© ESO, 2012

1. Introduction

The process of star formation has been a subject of active research over the past few decades. It involves the gravitational collapse of a cold dense core inside a molecular cloud, which heats up in its centre as the pressure and density increase from the compression, a problem which entails many physical processes (hydrodynamics, radiative transfer, magnetic fields, etc.) over a very wide range of spatial scales (over 5 orders of magnitude) (Larson 1969; Stahler et al. 1980; Masunaga et al. 1998). The collapsing material is initially optically thin to the thermal emission from the cold gas and dust grains and the entire energy gained from compressional heating is transported away by the escaping radiation, which causes the cloud to collapse isothermally in the initial stages of the protostar formation. When the optical depth of the cloud reaches unity, the radiation is absorbed by the system, which starts heating up, taking the core collapse through its adiabatic phase. The strong compression forms an accretion shock at the border of the adiabatic core (also known as Larson’s first core).

This radiative shock is believed to play a major role in the star formation process since it can either convert the kinetic energy of the material infalling onto the core into thermal energy, which is deposited inside the core, or radiate this energy away, depending on the strength of the shock and the opacity of the gas. This largely determines the key properties of the first core such as its size, mass, and entropy level as well as the magnitude of the radiative flux at the core border (Winkler & Newman 1980; Commerçon et al. 2011). The presence of this radiative shock requires calculations to include the full radiation hydrodynamics (RHD) system of equations, which makes computations demanding. Three-dimensional (3D) RHD simulations have only just recently become possible with modern computers. The first 3D study of the collapse of an isolated spherical cloud using smoothed particle hydrodynamics and flux-limited diffusion (FLD) was carried out by Whitehouse & Bate (2006) (see Whitehouse & Bate 2004, for a description of the methodology). Krumholz et al. (2007a) later used a grid-based adaptive mesh refinement RHD code to simulate star formation in a turbulent molecular cloud, using FLD in the mixed frame formalism (Krumholz et al. 2007b). More recently, with the continuous rise in computer power, several authors have been able to simulate star formation in increasingly larger turbulent clouds (Bate 2009; Offner et al. 2009), but these largescale 3D studies can never incorporate the same detailed physics as one-dimensional (1D) studies, which are still highly essential to uncover and isolate all the crucial physical processes at work.

In this paper, we improve on the recent 1D calculations of Commerçon et al. (2011) by including frequency-dependent radiative transfer in our model as opposed to the more simple grey approximation, which integrates the equations of radiative transfer over the entire frequency range. The motivation behind this extension lies in the strong dependence of the circumstellar material opacities as a function of frequency. Indeed, considerable changes in opacities (several orders of magnitude!) from the radio to the X-ray part of the electromagnetic spectrum are ubiquitous to interstellar gas and dust (see for example Ossenkopf & Henning 1994; Li & Draine 2001; Draine 2003; Semenov et al. 2003; and Fig. 1). This implies that the cloud core and circumstellar envelope can be optically thick to some parts of the spectrum and optically thin to others at the same time. This has potentially important consequences for the evolution and properties of the protostellar core.

Fig. 1 Spectral opacities for homogeneous spherical dust grains taken from Semenov et al. (2003) for five different temperature ranges. The grey colour bands represent the five different frequency groups that were used in the multigroup simulation. The dotted vertical lines indicate the additional splitting of groups 2 and 3 in the case of the 12-group simulation. The dashed line represents the Planck function for a black body with a typical temperature of 1000 K normalised to arbitrary units. The numbers in the lower part of the figure identify the different groups used in the 12-group simulation (see Sect. 3).

Of all the star formation studies cited above, Masunaga et al.’s (1998) is the only one not to use the grey approximation. They implemented frequency-dependent radiative transfer in their 1D simulations by solving the transfer equation directly. They looked at exactly the same problem as the one discussed in this paper, and already addressed the astrophysical questions 15 years ago. However, their method cannot be realistically incorporated in future 3D simulations because of the overly complex nature of the full seven-dimensional radiative transfer equation. Yorke & Sonnhalter (2002) used a simpler multi-frequency FLD in 2D in the context of high mass (>30  M_⊙) star formation through gravitational collapse, but once again, solving the radiative transfer at every frequency would be too computationally expensive in 3D. On the other hand, the multigroup approach, for which the frequency domain is divided into a finite number of bins or groups and the opacities are averaged within each group, enables us to limit the number of frequency bins. This proves to be a good compromise between reproducing the important aspects of frequency-dependent opacities and saving computational cost. Moreover, using limited numbers of frequency groups allows us to use a more accurate radiative transfer model than the FLD, such as the M₁ model (see below).

We first describe the multigroup radiative transfer model, the opacities and the numerical method used in the simulations. The results are then presented and compared to grey simulations, highlighting the differences between the two schemes.

Fig. 2 Planck (left) and Rosseland (right) mean opacities as a function of temperature for the grey (thick grey) and multigroup (colours) simulations. The sharp discontinuities seen in the curves are due to transitions from one temperature range to the other, corresponding to the destruction of ice, silicates, etc.

2. The multigroup RHD collapse simulations

2.1. The equations of multigroup radiation hydrodynamics

The numerical method employed to solve the equations of RHD is the one described in Vaytet et al. (2011). We use the M₁ moment model for radiative transfer (Levermore 1984; Dubroca & Feugeas 1999; González et al. 2007) in its multigroup formalism, where the equations of radiative transfer are integrated into a finite number of frequency groups and the opacities are averaged over the same frequency ranges. This allows us to account for frequency dependence of the absorption and emission coefficients. The system of multigroup RHD equations in the comoving frame (all radiative quantities, including the frequencies, are expressed in the frame of reference moving with the fluid) is where ρ, u, p and e are the gas density, velocity, pressure and total energy, respectively, and Φ is the gravitational potential. We also define (6)where X_g = E_g, F_g, P_g which represent the radiative energy, flux and pressure inside each group g which holds frequencies between ν_g − 1/2 and ν_g + 1/2. Q is the third moment of the radiation specific intensity (representing a heat flux), N_g is the total number of groups and Θ_g(T) is the energy of the photons with a Planck distribution at temperature T inside a given group. The absorption coefficients σ_Pg, σ_Eg and σ_Fg are the means of σ_ν inside a given group weighted by the Planck function, the radiative energy and the radiative flux, respectively. The opacity κ_ν of the gas is defined as σ_ν = ρκ_ν. As in Mihalas & Mihalas (1984), the notation P:∇u represents the tensorial product P_ij∂ⁱu^j.

2.2. The opacities

At low temperatures (below 1500 K), the opacities of the interstellar gas are dominated by the one percent dust grains and the atomic and molecular gas opacities can be neglected. We used the monochromatic opacities from Semenov et al. (2003)¹, who provided dust opacities for various types of grains in five different temperature ranges (we assumed that the dust opacities are independent of gas density and that the dust is in thermal equilibrium with the gas; see Galli et al. 2002 for example). The spectral opacities for homogeneous spherical dust grains and normal iron content in the silicates (Fe/(Fe + Mg) = 0.3) are shown in Fig. 1; the different line colours correspond to the different temperature ranges explicated in the legend. The grey-shaded areas represent the five frequency groups used in the multigroup simulation: radio + microwave, far-IR, IR, visible and UV (see Table 1).

Planck (κ_P) and Rosseland (κ_R) mean opacities were calculated over the entire frequency range for the grey simulation and inside each frequency group for the multigroup simulation. We show in Fig. 2 the mean opacities that were tabulated for temperatures between 5 K and 1500 K using ~10 K-wide linear transition between the five different temperature regimes, which correspond to the destruction of ice, silicates, etc. (see Semenov et al. 2003). In the RHD Eqs. (2) − (5), common practice is to set σ_Eg = σ_Pg and σ_Fg = σ_Rg, and this is also what we have used in the present work (see Larsen & Lane 1994; González et al. 2007; Offner et al. 2009, for example). However, we wish to point out that the inaccuracies that arise from this approximation are reduced as the number of frequency groups used increases, since in the limit of infinite frequency resolution, all these quantities simply reduce to σ_ν. This approximation is thus less crude in a multigroup model than in a grey model.

2.3. Initial and boundary conditions

The initial setup for the dense core collapse is identical to the one of Commerçon et al. (2011) and Masunaga et al. (1998). A uniform density sphere of mass M₀ = 1  M_⊙, temperature T₀ = 10 K (c_s0 = 0.187  km   s^-1) and radius R₀ = 10⁴ AU collapses under its own gravity. The ratio of thermal to gravitational energy in the cold gas cloud is (7)and its free-fall time is t_ff ~ 0.177 Myr. The radiation temperature is in equilibrium with the gas temperature (the energy of a black body with T = 10 K is divided among the frequency groups according to the Planck distribution) and the radiative flux is set to zero everywhere. The boundary conditions are reflexive at the centre of the grid (r = 0) and have imposed values equal to the initial conditions at the outer edge of the sphere. We used an ideal gas equation of state with the ratio of specific heats γ = 5/3 and the mean molecular weight μ = 2.375 for molecular hydrogen and a helium concentration of 0.27.

We also performed simulations of the collapse of 10 and 0.1 M_⊙ clouds (see Sect. 3.2), where the setup is almost identical to the 1 M_⊙ case. The thermal to gravitational energy ratio α was kept the same and the parent cloud radii were now R₀ = 10⁵ and 10³ AU, respectively. The free-fall times were then 10 times longer for the 10 M_⊙ case and 10 times shorter for the 0.1 M_⊙ case.

2.4. Numerical method

The numerical method is the spherically symmetric version of the one described in Vaytet et al. (2011); a 1D Lagrangean second-order Godunov RHD code. The hydrodynamics equations are solved using a classical MUSCL-Hancock scheme while the radiative transfer equations are solved using a HLLC solver with an asymptotic preserving correction (Berthon et al. 2007). The RHD equations are integrated with a scheme implicit in time, using a standard Raphson-Newton iterative method and the lapack library for the Jacobian matrix inversion. The radial (spherically symmetric) grid comprises 2000 cells and is logarithmically regular. The regridding scheme of Dorfi & Drury (1987) was also used (decoupled from the hydrodynamics) to maintain a good resolution at the accretion shock.

3. Results

3.1. The 1 M_⊙ case

A grey and a multigroup simulation of the collapse of a dense core were performed. In the grey run, the radiative quantities were integrated over the entire frequency range (0 to 3 × 10¹⁵ Hz) while for the multigroup run, five frequency groups were used (see Table 1). The grey run should be more or less equivalent to the simulation in Commerçon et al. (2011), because a very similar method was used (the opacities and the computational grid differ slightly and the solver for the radiative transfer is different). The simulations were run until the central density reached ρ_c = 10^-8  g   cm^-3. Figure 3 shows the radial profiles of the density, temperature, velocity, entropy, optical depth, luminosity, radiative flux and mass for the grey (grey thick solid line) and multigroup simulations (colours).

We first note that the results for the grey simulation are again very consistent with the results of Masunaga et al. (1998), as was noted in Commerçon et al. (2011). An adiabatic core has formed at the centre of the grid with a shock clearly visible at its border, at a radius of ~7 AU (see Figs. 3a and 3c for instance). From the temperature plot (b), we see that this is a supercritical radiative shock; the pre- and post-shock temperatures are identical and the radiative precursor is visible in between 7 and 100 AU. The properties of the first core are listed in Table 2; these are the first core radius (R_fc), the mass of the first core (M_fc), the mass accretion rate at the accretion shock (8)the accretion luminosity (9)the total luminosity (L_tot), the gas temperature at the centre of the core (T_c), the gas temperature at the accretion shock T_fc and the entropy at the centre of the core (10)These values are very close to those listed in Commerçon et al. (2011), which is once again expected since we used almost the same computational method.

Fig. 3 Radial profiles of various properties during the collapse of a 1  M⊙ dense clump at a core central density ρc = 10-8  g   cm-3 for the grey (grey thick solid line) and multigroup (colours) models. For the multigroup run, the gas and total radiative (summed over all groups) quantities are plotted in red, while the other colours represent the individual groups; see legend in e). From top left to bottom right: a) density; b) gas (solid) and radiation (dashed) temperature; c) velocity; d) entropy; e) optical depth; f) luminosity; g) Planck average opacity; and h) reduced radiative flux.

Fig. 4 Same as Fig. 3 but for the 12-group simulation.

Fig. 5 Comparison of the radial profiles during the collapse of a 0.1  M⊙, 1 M⊙, and 10 M⊙ dense cloud for a core central density ρc = 10-8  g   cm-3. From top left to bottom right: a) density; b) gas temperature; c) velocity; d) entropy; e) radiative flux; and f) enclosed mass. The legend in the f) panel explains the curve colour-coding: the top row is the mass of the initial cloud and the bottom row is the number of frequency groups used in each simulation.

There are significant variations in radiative flux (h) (and hence luminosity) inside the first core, which occur where there are large jumps in dust opacity (g) when changing from one temperature range to the next (destruction of ice, pyroxene and other dust grain constituents, see Semenov et al. 2003), because the variation in temperature is important in this part of the system. The jump in luminosity at the position of the shock is due to the accretion luminosity, which is added to the internal luminosity of the first core to constitute the total luminosity. The power radiated per unit area is (11)and therefore the radiative luminosity just inside of the shock is for in the case of one frequency group (we recall that in this case). The accretion luminosity is L_acc = 2.16 × 10^-2  L_⊙ (see Table 2) and the sum of the two approximately corresponds to the total radiative luminosity observed just outside of the shock L_tot = 4.2 × 10^-2  L_⊙ ≃ L_rad + L_acc. The total luminosity comprises an additional 3 × 10^-3  L_⊙ because the optical depth per cell falls abruptly at the shock because of the drop in gas density, and more radiation thus escapes from the core, causing a mini-burst in luminosity.

Secondly, and much more importantly, we note that there are virtually no differences between the grey and multigroup simulations. All profiles are identical to their grey counterparts. We see that the IR frequency group (3) dominates inside the first core (radiative temperature and flux) while the far-IR group (2) takes over in the outer envelope. One also notes that the grey Planck mean opacity (grey thick solid line in g) is driven towards the values in these groups by the weighting of the Planck function for the temperatures of the core (200 − 1200 K) and the outer envelope (10 − 200 K). As listed in Table 2, the first core radius is slightly larger and the central temperature is marginally higher in the multigroup simulation, but these differences are of the order of one percent, hence insignificant for the global properties of the cores.

There is no jump in gas temperature at the accretion shock in either the grey or the multigroup simulations, illustrating that the radiative shock is supercritical in both cases. This and the upstream Mach number of the shock (M ~ 2) implies that the vast majority of the kinetic energy from the infalling material is radiated away at the shock (Commerçon et al. 2011).

In light of these small differences between the grey and five-group simulation, we ran a third simulation, this time using 12 frequency groups. Since groups 2 and 3 hold the majority of the radiative energy (and show strong opacity variations), they were split again while groups 1, 4 and 5 were kept unchanged. Group 2 was split into three logarithmically equal parts and group 3 into six parts, as illustrated by the dotted vertical lines in Fig. 1. The results for the 12-group simulation are compared to the grey run in Fig. 4 and the first core properties are listed in Table 2.

There are again no major differences between the 12-group and the two previous simulations. The first core radius is slightly larger in the 12-group case (difference of 10% compared to the grey run) and the central temperature is also marginally higher, but the other first core properties are very similar to the five-group simulation (the time needed to reach ρ_c was also 0.194 Myr; identical to the five-group run). Panels (b), (f) and (h) illustrate well the sharing of the radiative energy and flux among the different frequency groups. Groups 5−10 hold the majority of the radiative energy inside the core, where for the most part, temperatures exceed 600 K. Figure 4g shows that groups 5, 6 and 10 have a significantly higher opacity than the grey opacity, and that the opacities groups 7 and 8 are lower. This is also illustrated by Fig. 1 where groups 5 and 6 correspond to peaks in spectral opacity (light blue curve) while 7 and 8 are in the low-opacity region immediately following. The spectral opacities increase again in groups 9 and 10. This means that inside the core, there is a combination of high and low opacities, and that the grey opacity lies somewhere in the middle. This explains why the central temperatures differ between the grey and the 12-group case, and in this particular instance, the higher opacities in groups 5, 6 and 10 “win” over the lower opacities in 7 and 8, preventing a little more radiation from escaping and causing the core to be somewhat hotter. The heating of the core works against the gravitational contraction, allowing the core to have a larger radius.

Overall, however, the variations in dust opacity as a function of frequency are not sharp enough to induce a strong impact on the results, compared to the grey simulation, which weighs the Planck opacity using the Planck function (and its derivative with respect to temperature in the case of the Rosseland opacity) over the entire frequency range and adequately reproduces the variation of the opacity as a function of gas temperature. Significant differences between grey and multigroup simulations would appear if strong variations in opacities (several orders of magnitude) occured on scales much smaller that the typical width of the Planck function.

3.2. The 0.1 M_⊙ and 10 M_⊙ cases

To check the universality of the results above, we also performed simulations of the collapse of a 0.1 M_⊙ and a 10 M_⊙ cloud. As mentioned in Sect. 2.3, the initial setup is identical to the 1 M_⊙ case, except for the system’s initial radius. The same opacities and frequency groups were also used. The results for the simulations using 1, 5 and 12 frequency groups are shown in Fig. 5, along with the previous results from the 1 M_⊙ case.

The first aspect to notice is that the results are strikingly similar to the 1 M_⊙ case, for the 0.1 M_⊙ and 10 M_⊙ clouds. The properties of the first core for the new cloud masses are listed in Table 2, and this illustrates the point further. The shock radii and masses of the first core are very similar (differences of ~20%). Some trends do appear in the results; for all core properties, except for the entropy, we note some dependence upon the mass of the parent cloud, although weak, this effect is the most obvious for the core’s central temperature (see Fig. 5b). Figure 5e also shows the more massive parent clouds producing slightly more luminous cores, as was noted by Masunaga et al. (1998).

The second trend in the results listed in Table 2 is that the 12-group simulations always appear to yield hotter, larger and more massive cores compared to the grey runs, but once again those differences remain very small. Overall, the properties of the first core, and especially the entropy inside the core, appear to be independent of the size and mass of the initial cloud (cf. Masunaga et al. 1998). There are here again no major differences between the grey and multigroup simulations for all the runs.

4. Conclusions

We have performed multigroup RHD simulations of the gravitational collapse of a 1 M_⊙ cold dense cloud core up to the formation of the first Larson core, reaching a central density of ρ_c = 10^-8  g   cm^-3. Five groups were first used to sample the opacities in the frequency domain and the results were compared to a grey simulation. Only small differences were found between the two runs, with no major structural or evolutionary changes. The main properties of the resulting first core formed in the centre of the grid such as its mass, size and entropy were not changed significantly. We then performed a third simulation using 12 frequency groups, which again did not affect the results in any major fashion. The grey simulation was found to slightly underestimate the first core’s radius, mass and central temperature, however.

We also performed simulations of the collapse of a 0.1 and 10 M_⊙ parent cloud to confirm the robustness of the results stated above. The properties of the first core were found to be quasi-independent of the initial conditions in the cloud; all cores have a radius of ~7 AU, a mass of ~ 2 × 10^-2  M_⊙ and an entropy of ~2 × 10⁹  erg   K^-1   g^-1 at their centre. As for the dependency on the radiative transfer model used, there were once again no significant differences between the grey and multigroup simulations in the 0.1 and 10 M_⊙ cases.

Nevertheless, we did note that the more frequency groups used, the larger the (small) differences were, and that the multigroup simulations always appeared to yield hotter (differences of 1%), larger (differences of ~ 10%) and more massive (differences of 3%) cores compared to the grey runs. The multigroup simulations also revealed some interesting facts on the distribution of the radiation energy and flux among the different frequency groups. Indeed, the vast majority of the radiative energy and flux inside the first core is contained within the IR frequency group, while the far-IR radiation dominates in the outer envelope. This could help identify first cores in IR observations of giant molecular clouds. We wish to point out here that Masunaga et al. (1998) were already using frequency-dependent radiative transfer by solving the transfer equation directly, but their method cannot be realistically used for future 3D simulations, unlike our new multigroup scheme.

The grey solution appears to work well for the first stage of collapse since the medium is optically thick for virtually all frequency groups behind the accretion shock (see Fig. 4e) and the gas and radiation temperatures are very closely coupled. However, differences are expected to arise if the gas and radiation are decoupled, and if radiation of very different energies is present at the same time in the system, as is the case for the second phase of the collapse. Indeed, once the temperatures are high enough in the core, the dissociation of molecular hydrogen enables further rapid collapse, giving rise to the situation pictured in Stahler et al. (1980) where a very hot core mainly produces UV and X-ray radiation, which travels through an opacity gap (created from the sublimation of the circumstellar dust at temperatures above 1500 K) and hits the dust front, where it is absorbed (the dust is optically thick to UV radiation) and re-emitted into the IR. This situation can evidently only be realistically modelled with a frequency-dependent treatment of the radiative transfer.

The second phase of collapse requires using a more sophisticated equation of state compared to the ideal gas equation and we will make the use of the Saumon-Chabrier-Van-Horn equation of state (Saumon et al. 1995) in an upcoming paper to model the complicated physics involved. The frequency-dependent radiative transfer will also allow us to study the radiative flux at the accetion shock in great detail and confirm or invalidate the predictions on the inward/outward flux budget of Stahler et al. (1980), which have important consequences on the rate at which star formation occurs.

Acknowledgments

The research leading to these results has received funding from the European Research Council under the European Community’s Seventh Framework Programme (FP7/2007-2013 Grant Agreement No. 247060). B.C. greatfully acknowledges support from his postdoctoral fellowship at the Max-Planck-Institut für Astronomie and the support of the CNES postdoctoral fellowship program. The authors would like to thank the referee for useful comments, which have helped to improve this paper.

References

All Tables

All Figures

Fig. 1 Spectral opacities for homogeneous spherical dust grains taken from Semenov et al. (2003) for five different temperature ranges. The grey colour bands represent the five different frequency groups that were used in the multigroup simulation. The dotted vertical lines indicate the additional splitting of groups 2 and 3 in the case of the 12-group simulation. The dashed line represents the Planck function for a black body with a typical temperature of 1000 K normalised to arbitrary units. The numbers in the lower part of the figure identify the different groups used in the 12-group simulation (see Sect. 3).

In the text

	Fig. 2 Planck (left) and Rosseland (right) mean opacities as a function of temperature for the grey (thick grey) and multigroup (colours) simulations. The sharp discontinuities seen in the curves are due to transitions from one temperature range to the other, corresponding to the destruction of ice, silicates, etc.
In the text

Fig. 3 Radial profiles of various properties during the collapse of a 1  M⊙ dense clump at a core central density ρc = 10-8  g   cm-3 for the grey (grey thick solid line) and multigroup (colours) models. For the multigroup run, the gas and total radiative (summed over all groups) quantities are plotted in red, while the other colours represent the individual groups; see legend in e). From top left to bottom right: a) density; b) gas (solid) and radiation (dashed) temperature; c) velocity; d) entropy; e) optical depth; f) luminosity; g) Planck average opacity; and h) reduced radiative flux.

In the text

	Fig. 4 Same as Fig. 3 but for the 12-group simulation.
In the text

Fig. 5 Comparison of the radial profiles during the collapse of a 0.1  M⊙, 1 M⊙, and 10 M⊙ dense cloud for a core central density ρc = 10-8  g   cm-3. From top left to bottom right: a) density; b) gas temperature; c) velocity; d) entropy; e) radiative flux; and f) enclosed mass. The legend in the f) panel explains the curve colour-coding: the top row is the mass of the initial cloud and the bottom row is the number of frequency groups used in each simulation.

In the text