Issue 
A&A
Volume 530, June 2011



Article Number  A13  
Number of page(s)  15  
Section  Interstellar and circumstellar matter  
DOI  https://doi.org/10.1051/00046361/201016213  
Published online  27 April 2011 
Physical and radiative properties of the firstcore accretion shock^{⋆}
^{1}
Max Planck Institute for Astronomy,
Königstuhl 17,
69117
Heidelberg,
Germany
email: benoit@mpiahd.mpg.de
^{2}
Laboratoire AIM, CEA/DSM  CNRS  Université Paris Diderot,
IRFU/SAp,
91191
GifsurYvette,
France
^{3}
École Normale Supérieure de Lyon, CRAL, UMR 5574 CNRS, Université
de Lyon, 46 allée
d’Italie, 69364
Lyon Cedex 07,
France
^{4}
School of Physics, University of Exeter,
Exeter
EX4 4QL,
UK
Received:
26
November
2010
Accepted:
14
February
2011
Context. Radiative shocks play a dominant role in star formation. The accretion shocks on first and second Larson cores involve radiative processes and are thus characteristic of radiative shocks.
Aims. In this study, we explore the formation of the first Larson core and characterize the radiative and dynamical properties of the accretion shock, using both analytical and numerical approaches.
Methods. We developed both numerical radiationhydrodynamics calculations and a semianalytical model that characterize radiative shocks in various physical conditions, for radiating or barotropic fluids. Then, we performed 1D spherical collapse calculations of the first Larson core, using a grey approximation for the opacity of the material. We considered three different models for radiative transfer: the barotropic approximation, the flux limited diffusion approximation, and the more complete M1 model. We investigate the characteristic properties of the collapse and of the first core formation. Comparison between the numerical results and our semianalytical model for radiative shocks shows that the latter reproduces the core properties obtained with the numerical calculations quite well.
Results. The accretion shock on the first Larson core is found to be supercritical; i.e., the post and preshock temperatures are equal, implying that all the accretion shock energy on the core is radiated away. The shock properties are described well by the semianalytical model. The fluxlimited diffusion approximation is found to agree quite well with the results based on the M1 model of radiative transfer, and is thus appropriate for studying the star formation process and allows a tractable and relatively correct treatment of radiative transfer in multidimensional radiationhydrodynamics calculations. In contrast, the barotropic approximation does not correctly describe the thermal properties of the gas during the collapse.
Key words: stars: formation / methods: analytical / methods: numerical / hydrodynamics / radiative transfer
Appendix A is available in electronic form at http://www.aanda.org
© ESO, 2011
1. Introduction
Star formation involves a wide variety of complex physical processes. Among them, radiative transfer appears to be one of the most important ones, since it governs the behaviour of the earliest phases of the collapse and the formation of the socalled first core (Larson 1969). It has been well established that the molecular gas experiences different thermal regimes from isothermal to adiabatic during the earliest phase of the collapse via the coupling between gas, dust, and radiation (e.g. Larson 1969; Tscharnuter & Winkler 1979; Masunaga et al. 1998).
Since the pioneering works of Winkler & Newman and Tscharnuter (Winkler & Newman 1980; Tscharnuter 1987), relatively few studies have been devoted to the accretion shock on the first Larson core. This shock is a radiative shock; i.e., radiation can escape from the infalling material. Once this infalling gas has been shocked and has reached an optical depth of 1, it becomes more or less adiabatic, at a given entropy level characteristic of the first core’s initial energy content. This entropy content is kept roughly constant during the nearly adiabatic subsequent stages of the collapse and the formation of the second core. The accretion shock on the first core is thus of prime importance, since it determines the entropy level of the subsequent stages of star formation up to the formation of the protostar itself.
A radiative shock is a shock that is so strong that it emits radiation that in turn affects the hydrodynamic behaviour of the flow. The equations of radiationhydrodynamics (RHD) must thus be solved consistently to properly explore this kind of a process. Radiative shocks are common in astrophysics, e.g. in star formation regions or around supernovae remnants. They have been studied well from the theoretical, experimental, and numerical points of view (e.g. Bouquet et al. 2000; Drake 2007; González et al. 2009). The wide variety of radiative shocks makes it difficult to classify them (Drake 2005; Michaut et al. 2009).
In this work, we study radiative shocks both numerically and analytically. We focus on the properties of the accretion shock during the formation of the first prestellar core. In the first part of the paper, we recall the main properties of radiative shocks. Then, considering the jump relations for a radiating fluid, we focus on some particular kinds of radiative shocks, with upstream and downstream material characterized by different optical depth properties. We also consider the case of a shock taking place in a barotropic material, since the barotropic approximation is the simplest way to characterize the thermal evolution of the gas during the collapse (e.g. Commerçon et al. 2008).
Fig. 1 Temperature and density profiles in a subcritical shock (left) and a supercritical shock (right). Adapted from Zel’Dovich & Raizer (1967). 
In the second part, we study the impact of various models of radiative transfer on the protostellar collapse. Radiation hydrodynamics plays a crucial role here, for instance by evacuating the compressional energy, leading to a nearly isothermal freefall collapse phase. Radiation transfer also has a dramatic impact on the accretion shock, with the infalling gas kinetic energy either converted into internal energy in the static adiabatic core or radiated away. This ratio between the energy accreted onto the star and the one radiated away is a key quantity, as it ultimately determines the entropy content of the forming star. Using a 1D spherical code, we compare results of calculations done with a barotropic EOS, a fluxlimited diffusion approximation, and a M1 model for radiative transfer, assuming grey opacity for the infalling material. The 1D calculations retain the important virtue of allowing detailed and complex physical processes to be considered in dense core collapse calculations, which is not the case in a multidimensional approach. In consequence, a 1D approach enables us to characterize the impact of these various physical processes on the results, allowing a better characterization of the validity of simplified multidimensional studies (e.g. Commerçon et al. 2010, 2011).
This paper is organized as follows. In the first section, we recall the basic properties of radiative shocks and develop semianalytical models that can be applied in some particular cases for a radiating fluid and a barotropic fluid. In Sect. 3, we detail the numerical method and the physics input used in our calculations. In Sect. 4, we derive the firstLarson core properties from the numerical calculations, using various approximations for the radiative transfer, and compare the results with the ones obtained by Masunaga et al. (1998). We also present an original semianalytical model that allows a simple interpretation of the numerical results. Eventually, Sect. 5 summarizes our main results.
2. Radiative shock – A semianalytic model
2.1. A qualitative picture of radiative shocks
Depending on the shock’s strength, radiative shocks belong to two groups: the subcritical shocks and the supercritical shocks. As the strength of a shock increases, the postshock temperature T_{2} rises, producing a radiative flux of $\mathit{\sigma}{\mathit{T}}_{\mathrm{2}}^{\mathrm{4}}$, which increases very rapidly. This flux penetrates the upstream material and preheats it to a temperature T_{−} immediately ahead of the shock front (radiative precursor), which is proportional to the incident flux, then T_{−} increases rapidly with the shock strength and eventually becomes equal to T_{2}. A shock with T_{−} < T_{2} is called a subcritical shock. Because the material entering the shock is preheated, the postshock temperature T_{ + } exceeds its asymptotic equilibrium value T_{2}, and decays downstream as the material cools by emitting photons that propagate across the shock (see Fig. 1 for various shock characteristics).
For stronger shocks, the preheating becomes so strong that the preshock temperature equals the postshock equilibrium temperature T_{2}. The shock velocity at which the postshock and preshock temperatures are equal defines the critical shock. For higher shock velocities, T_{−} cannot exceed T_{2}, the excess energy forces the radiative precursor further into the upstream region with a temperature close to T_{2}. Pre and postshock temperatures are equal, so that the supercritical shock is isothermal. Radiation and matter are still out of equilibrium in some part of the precursor, but balance out when temperature approaches T_{2}. The upstream kinetic energy is radiated at the shock.
Most of the early work on radiative shocks has been described in Zel’Dovich & Raizer (1967) and Mihalas & Mihalas (1984), where readers can find the basic equations of the front structure, radiative precursor extension, etc.
2.2. Jump relations for a radiating material
Consider the jump relations (Rankine Hugoniot) for a radiating flow (see Mihalas & Mihalas 1984; Zel’Dovich & Raizer 1967; Drake 2007). Conservation of mass, momentum, and energy yields $\begin{array}{ccc}& & {\mathit{\rho}}_{\mathrm{1}}{\mathit{u}}_{\mathrm{1}}\mathrm{=}{\mathit{\rho}}_{\mathrm{2}}{\mathit{u}}_{\mathrm{2}}\mathrm{\equiv}\mathit{m\u0307}\mathit{,}\\ & & {\mathit{\rho}}_{\mathrm{1}}{\mathit{u}}_{\mathrm{1}}^{\mathrm{2}}\mathrm{+}{\mathit{P}}_{\mathrm{1}}\mathrm{+}{\mathit{P}}_{\mathrm{r}\mathrm{1}}\mathrm{=}{\mathit{\rho}}_{\mathrm{2}}{\mathit{u}}_{\mathrm{2}}^{\mathrm{2}}\mathrm{+}{\mathit{P}}_{\mathrm{2}}\mathrm{+}{\mathit{P}}_{\mathrm{r}\mathrm{2}}\mathit{,}\\ & & \mathit{m\u0307}\mathrm{(}{\mathit{h}}_{\mathrm{1}}\mathrm{+}{\mathit{\rho}}_{\mathrm{1}}{{\mathit{u}}_{\mathrm{1}}^{\mathrm{2}}}^{\mathrm{)}}\mathrm{+}{\mathit{F}}_{\mathrm{r}\mathrm{1}}\mathrm{+}{\mathit{u}}_{\mathrm{1}}\mathrm{(}{\mathit{E}}_{\mathrm{r}\mathrm{1}}\mathrm{+}{\mathit{P}}_{\mathrm{r}\mathrm{1}}\mathrm{)}\mathrm{=}\\ & & \mathit{m\u0307}\mathrm{(}{\mathit{h}}_{\mathrm{2}}\mathrm{+}{\mathit{\rho}}_{\mathrm{2}}{{\mathit{u}}_{\mathrm{2}}^{\mathrm{2}}}^{\mathrm{)}}\mathrm{+}{\mathit{F}}_{\mathrm{r}\mathrm{2}}\mathrm{+}{\mathit{u}}_{\mathrm{2}}\mathrm{(}{\mathit{E}}_{\mathrm{r}\mathrm{2}}\mathrm{+}{\mathit{P}}_{\mathrm{r}\mathrm{2}}\mathrm{)}\mathit{,}\end{array}$where subscripts “1” and “2” respectively denote the upstream and downstream states, all radiation quantities are estimated in the comoving frame, h corresponds to the gas specific enthalpy (h = ϵ + P /ρ, with e = ρϵ = P /(γ − 1)), and ṁ is the mass flux through the shock. In comparison with the hydrodynamical case, the pressure is now the total pressure; i.e., the gas plus radiation contributions, is P + P_{r}, and the specific enthalpy is the total specific enthalpy, ϵ + (P + P_{r} + E_{r}) /ρ. The radiation energy E_{r} and radiation pressure P_{r} are important only at high temperatures or low densities, whereas radiation flux F_{r} plays a fundamental role in all radiative shocks.
Contrarily to the hydrodynamical case, the system of Eqs. (1)–(3) can not be solved explicitly. We need to make some assumptions on both the upstream and downstream materials. In what follows, we distinguish two cases: i) the upstream and downstream materials are opaque; ii) the upstream material is optically thin and the downstream material remains opaque.
2.2.1. Radiative shock in an optically thick medium
This is the most studied case (e.g. Mihalas & Mihalas 1984; Drake 2007). It occurs, for instance, in shocks in stellar interiors, where matter is both dense and hot. At a sufficiently large distance from the front, matter and radiation are in equilibrium and both the downstream and upstream materials are opaque (F_{r1} = F_{r2} = 0). Any radiation crossing the front from the hot downstream material into the cool upstream material is reabsorbed in the radiative precursor, into which it propagates by diffusion. Outside this diffusion layer, the flux vanishes, and we have $\begin{array}{ccc}\mathit{m\u0307}\mathrm{(}{\mathit{h}}_{\mathrm{1}}\mathrm{+}{\mathit{\rho}}_{\mathrm{1}}{{\mathit{u}}_{\mathrm{1}}^{\mathrm{2}}}^{\mathrm{)}}\mathrm{+}{\mathit{u}}_{\mathrm{1}}\mathrm{(}{\mathit{E}}_{\mathrm{r}\mathrm{1}}\mathrm{+}{\mathit{P}}_{\mathrm{r}\mathrm{1}}\mathrm{)}& \mathrm{=}& \mathit{m\u0307}\mathrm{(}{\mathit{h}}_{\mathrm{2}}\mathrm{+}{\mathit{\rho}}_{\mathrm{2}}{{\mathit{u}}_{\mathrm{2}}^{\mathrm{2}}}^{\mathrm{)}}\\ & & \mathrm{+}{\mathit{u}}_{\mathrm{2}}\mathrm{(}{\mathit{E}}_{\mathrm{r}\mathrm{2}}\mathrm{+}{\mathit{P}}_{\mathrm{r}\mathrm{2}}\mathrm{)}\mathit{.}\end{array}$(4)Since matter and radiation are in equilibrium and the material is opaque, we have E_{r} = 3P_{r} = a_{R}T^{4}. Defining the compression ratio r = ρ_{2} /ρ_{1} = u_{1} /u_{2}, we can rewrite Eqs. (2) and (3) in the nondimensional form $\mathit{\gamma}{\mathrm{\mathcal{M}}}_{\mathrm{1}}^{\mathrm{2}}\left(\frac{\mathit{r}\mathrm{}\mathrm{1}}{\mathit{r}}\right)\mathrm{=}\mathrm{(}\mathrm{\Pi}\mathrm{}{\mathrm{1}}^{\mathrm{)}}\mathrm{+}{\mathit{\alpha}}_{\mathrm{1}}\left({\left(\frac{\mathrm{\Pi}}{\mathit{r}}\right)}^{\mathrm{4}}\mathrm{}\mathrm{1}\right)\mathit{,}$(5)and $\frac{\mathit{\gamma}}{\mathrm{2}}{\mathrm{\mathcal{M}}}_{\mathrm{1}}^{\mathrm{2}}\left(\frac{{\mathit{r}}^{\mathrm{2}}\mathrm{}\mathrm{1}}{{\mathit{r}}^{\mathrm{2}}}\right)\mathrm{=}\frac{\mathit{\gamma}}{\mathit{\gamma}\mathrm{}\mathrm{1}}\left(\frac{\mathrm{\Pi}}{\mathit{r}}\mathrm{}\mathrm{1}\right)\mathrm{+}\mathrm{4}{\mathit{\alpha}}_{\mathrm{1}}\left(\frac{{\mathrm{\Pi}}^{\mathrm{4}}}{{\mathit{r}}^{\mathrm{5}}}\mathrm{}\mathrm{1}\right)\mathit{,}$(6)where Π = P_{2} /P_{1}, ${\mathit{\alpha}}_{\mathrm{1}}\mathrm{=}\mathrm{(}\mathrm{1}\mathit{/}\mathrm{3}\mathrm{)}{\mathit{a}}_{\mathrm{R}}{\mathit{T}}_{\mathrm{1}}^{\mathrm{4}}\mathit{/}{\mathit{P}}_{\mathrm{1}}$ and ℳ_{1} is the hydrodynamic Mach number, ℳ_{1} = u_{1} /c_{s1} = u_{1} /(γP_{1} /ρ_{1})^{1 /2}. The coupled Eqs. (5) and (6) are solved numerically to get the variations of Π and r as function of ℳ_{1} and α_{1}.
The compression ratio rises from r = 4 for a pure hydrodynamic shock with γ = 5 /3 to r = 7, which corresponds to the limiting compression ratio for a gas with γ = 4 /3, i.e pure radiation, a gas made of photons. The stronger the shock, the greater the radiative effects, even for a very low initial ratio of radiative pressure to gas pressure. When the upstream radiative pressure increases, the shock becomes immediately radiative since radiative pressure increases as T^{4}, whereas the gas pressure only increases linearly with temperature. As shown in Mihalas & Mihalas (1984), for a strong radiating shock, since the compression ratio is fixed, the temperature ratio only increases as ${\mathrm{\mathcal{M}}}_{\mathrm{1}}^{\mathrm{1}\mathit{/}\mathrm{2}}$, whereas it rises as ${\mathrm{\mathcal{M}}}_{\mathrm{1}}^{\mathrm{2}}$ in a nonradiating shock. If P_{r} /P < 10^{5}, the nonradiating fluid approximation becomes valid for an upstream Mach number ℳ_{1} < 10.
2.2.2. Radiative shock with an optically thin upstream material
If we now suppose that the shock is propagating in an optically thin material with an opaque downstream region, as in the case of star formation (e.g., Calvet & Gullbring 1998), the radiative energy and pressure in the lowmass star formation context can be neglected compared to the gas internal energy for the first core accretion shock.
We then consider the material and radiative quantities at the discontinuity, outside the spike in gas temperature (i.e., between regions with subscripts “2” and “–”). We have $\begin{array}{ccc}& & {\mathit{\rho}}_{\mathrm{}}{\mathit{u}}_{\mathrm{}}\mathrm{=}{\mathit{\rho}}_{\mathrm{2}}{\mathit{u}}_{\mathrm{2}}\mathrm{\equiv}\mathit{m\u0307}\mathit{,}\\ & & {\mathit{\rho}}_{\mathrm{}}{\mathit{u}}_{\mathrm{}}^{\mathrm{2}}\mathrm{+}{\mathit{P}}_{\mathrm{}}\mathrm{=}{\mathit{\rho}}_{\mathrm{2}}{\mathit{u}}_{\mathrm{2}}^{\mathrm{2}}\mathrm{+}{\mathit{P}}_{\mathrm{2}}\mathit{,}\\ & & \mathit{m\u0307}\mathrm{(}{\mathit{h}}_{\mathrm{}}\mathrm{+}{\mathit{\rho}}_{\mathrm{}}{{\mathit{u}}_{\mathrm{}}^{\mathrm{2}}}^{\mathrm{)}}\mathrm{=}\mathit{m\u0307}\mathrm{(}{\mathit{h}}_{\mathrm{2}}\mathrm{+}{\mathit{\rho}}_{\mathrm{2}}{{\mathit{u}}_{\mathrm{2}}^{\mathrm{2}}}^{\mathrm{)}}\mathrm{+}\mathrm{\Delta}{\mathit{F}}_{\mathrm{r}}\mathit{,}\end{array}$where ΔF_{r} = F_{r2} − F_{r−}. We consider the case of a nonzero net flux F_{r2} − F_{r−} across the shock. The jump relations, derived from the conservation Eqs. (1)–(3), read $\frac{{\mathit{\rho}}_{\mathrm{2}}}{{\mathit{\rho}}_{\mathrm{}}}\mathrm{=}\frac{\left[\mathrm{(}\mathit{\gamma}\mathrm{+}\mathrm{1}\mathrm{)}{\mathit{P}}_{\mathrm{2}}\mathrm{+}\mathrm{(}\mathit{\gamma}\mathrm{}\mathrm{1}\mathrm{)}{\mathit{P}}_{\mathrm{}}\right]{\mathit{u}}_{\mathrm{2}}\mathrm{+}\mathrm{2}\mathrm{(}\mathit{\gamma}\mathrm{}\mathrm{1}\mathrm{)}\mathrm{\Delta}{\mathit{F}}_{\mathrm{r}}}{\left[\mathrm{(}\mathit{\gamma}\mathrm{+}\mathrm{1}\mathrm{)}{\mathit{P}}_{\mathrm{}}\mathrm{+}\mathrm{(}\mathit{\gamma}\mathrm{}\mathrm{1}\mathrm{)}{\mathit{P}}_{\mathrm{2}}\right]{\mathit{u}}_{\mathrm{2}}}\mathit{,}$(10)and $\frac{{\mathit{T}}_{\mathrm{2}}}{{\mathit{T}}_{\mathrm{}}}\mathrm{=}\left(\frac{{\mathit{P}}_{\mathrm{2}}}{{\mathit{P}}_{\mathrm{}}}\right)\frac{\left[\mathrm{(}\mathit{\gamma}\mathrm{+}\mathrm{1}\mathrm{)}{\mathit{P}}_{\mathrm{}}\mathrm{+}\mathrm{(}\mathit{\gamma}\mathrm{}\mathrm{1}\mathrm{)}{\mathit{P}}_{\mathrm{2}}\right]{\mathit{u}}_{\mathrm{}}\mathrm{}\mathrm{2}\mathrm{(}\mathit{\gamma}\mathrm{}\mathrm{1}\mathrm{)}\mathrm{\Delta}{\mathit{F}}_{\mathrm{r}}}{\left[\mathrm{(}\mathit{\gamma}\mathrm{+}\mathrm{1}\mathrm{)}{\mathit{P}}_{\mathrm{2}}\mathrm{+}\mathrm{(}\mathit{\gamma}\mathrm{}\mathrm{1}\mathrm{)}{\mathit{P}}_{\mathrm{}}\right]{\mathit{u}}_{\mathrm{}}}\mathrm{\xb7}$(11)These two relations show that radiative energy transport (i.e. the radiative flux) across a shock can significantly alter the density, temperature and velocity profiles of the flow. Both upstream and downstream materials are affected over distances that depend on the material opacity. As mentioned before, the structure of a radiative shock is as follows: the upstream material is preheated by a radiation precursor, while the downstream material is cooled by radiative losses.
As for the previous opaque case for which we use an iterative solver, the downstream quantities cannot be derived analytically from the conservation relations for given upstream conditions. The radiative flux has to be known, and the result depends on the upstream flow.
The simplest case to study is a supercritical shock, where T_{−} = T_{2} and then h_{−} = h_{2}. Using the same adimensional parameters as for the opaque case, we have $\mathit{\gamma}{\mathrm{\mathcal{M}}}_{\mathrm{}}^{\mathrm{2}}\left(\frac{\mathit{r}\mathrm{}\mathrm{1}}{\mathit{r}}\right)\mathrm{=}\mathrm{\Pi}\mathrm{}\mathrm{1}\mathit{,}$(12)and $\frac{\mathit{\gamma}}{\mathrm{2}}{\mathrm{\mathcal{M}}}_{\mathrm{}}^{\mathrm{2}}\left(\frac{{\mathit{r}}^{\mathrm{2}}\mathrm{}\mathrm{1}}{{\mathit{r}}^{\mathrm{2}}}\right)\mathrm{=}\frac{\mathit{\gamma}}{\mathit{\gamma}\mathrm{}\mathrm{1}}\left(\frac{\mathrm{\Pi}}{\mathit{r}}\mathrm{}\mathrm{1}\right)\mathrm{+}\frac{\mathrm{\Delta}{\mathit{F}}_{\mathrm{r}}}{{\mathit{P}}_{\mathrm{}}{\mathit{u}}_{\mathrm{}}}\mathrm{\xb7}$(13)Since the shock is isothermal, r = Π and $\mathrm{\Pi}\mathrm{=}\mathit{\gamma}{\mathrm{\mathcal{M}}}_{\mathrm{}}^{\mathrm{2}}\mathit{.}$(14)Eventually, the radiative flux discontinuity, normalized to the upstream kinetic energy, is given by $\mathit{X}\mathrm{=}\frac{\mathrm{\Delta}{\mathit{F}}_{\mathrm{r}}}{\mathrm{0.5}{\mathit{\rho}}_{\mathrm{1}}{\mathit{u}}_{\mathrm{}}^{\mathrm{3}}}\mathrm{=}\frac{{\mathit{\gamma}}^{\mathrm{2}}{\mathrm{\mathcal{M}}}_{\mathrm{}}^{\mathrm{4}}\mathrm{}\mathrm{1}}{{\mathit{\gamma}}^{\mathrm{2}}{\mathrm{\mathcal{M}}}_{\mathrm{}}^{\mathrm{4}}}\mathit{,}$(15)where X thus represents the amount of incident kinetic energy radiated away at the shock front.
Fig. 2 Compression ratio r and fraction of kinetic energy radiated away at the shock, X, as a function of the upstream Mach number for a supercritical shock, in the case of an optically thin upstream material. 
Figure 2 shows the evolution of the adimensional compression ratio, r, and flux discontinuity, X, as a function of the upstream Mach number for a supercritical shock. As seen in the figure, the compression ratio does not saturate, in contrast to the case of an opaque material. As shown by Eqs. (10) and (11), a nonzero flux across the front, with F_{r2} > F_{r−}, increases the density jump, whereas it decreases the temperature jump. To have an isothermal shock at low Mach number, ℳ < 2, the downstream velocity is determinant since in that case the upstream kinetic energy is not entirely radiated away. In protostellar collapse calculations, the characteristic Mach number at the first core accretion shock is ℳ ~ 3, implying that almost all the kinetic energy (>99%) is radiated away at this stage.
Including radiative energy and radiative pressure in Eqs. (8) and (9) does not affect the dependence of the compression ratio upon the upstream Mach number. On the other hand, the dependence of the amount of kinetic energy radiated away becomes stronger: $\mathit{X}\mathrm{=}\frac{\mathrm{\Delta}{\mathit{F}}_{\mathrm{r}}}{\mathrm{0.5}{\mathit{\rho}}_{\mathrm{}}{\mathit{u}}_{\mathrm{}}^{\mathrm{3}}}\mathrm{=}\frac{{\mathit{\gamma}}^{\mathrm{2}}{\mathrm{\mathcal{M}}}_{\mathrm{}}^{\mathrm{4}}\mathrm{}\mathrm{1}}{{\mathit{\gamma}}^{\mathrm{2}}{\mathrm{\mathcal{M}}}_{\mathrm{}}^{\mathrm{4}}}\mathrm{(}\mathrm{1}\mathrm{+}\mathrm{4}{\mathit{\alpha}}_{\mathrm{1}}\mathrm{)}\mathit{.}$(16)The same analysis cannot be carried out for a subcritical shock with an optically thin upstream medium, since, in that case, no constraint allows us to close the system of Eqs. (8) and (9). One needs a prescription on the radiation flux or luminosity in the upstream material to close the system.
2.3. Super or subcritical shock?
In this section, we characterize in which regime the radiative shock takes place, depending on the upstream and downstream material properties.
2.3.1. Opaque material
In an opaque material, the simplest case is the one of a supercritical shock, in which the preshock gas ahead of the discontinuity is heated up to the postshock temperature (see Sect. 2.1.2). The radiative energy absorbed in the upstream region is only used to raise the gas temperature (Zel’Dovich & Raizer 1967, p. 536). Assuming for simplicity that the gas is neither compressed nor slowed down, that the radiative pressure is negligible compared to the gas pressure (valid for early low mass star formation stages), and that the upstream internal energy is negligible far from the shock (region with subscript “1”), the equation of energy reads as $\mathrm{}{\mathit{F}}_{\mathrm{r}}\mathrm{=}{\mathit{u}}_{\mathrm{}}{\mathit{\rho}}_{\mathrm{}}\mathit{\u03f5}\mathrm{\left(}{\mathit{T}}_{\mathrm{}}\mathrm{\right)}\mathit{,}$(17)with $\mathit{\u03f5}\mathrm{\left(}\mathit{T}\mathrm{\right)}\mathrm{\equiv}\mathit{\u03f5}\mathrm{=}\frac{\mathrm{1}}{\mathit{\gamma}\mathrm{}\mathrm{1}}\frac{{\mathit{k}}_{\mathrm{B}}\mathit{T}}{\mathit{\mu}{\mathit{m}}_{\mathrm{H}}}$. If we evaluate the flux just at the discontinuity, we find the maximum preheating temperature T_{−} from $\mathrm{\left}{\mathit{F}}_{\mathrm{r}}\mathrm{\right}\mathrm{~}\mathit{\sigma}{\mathit{T}}_{\mathrm{2}}^{\mathrm{4}}\mathrm{=}{\mathit{u}}_{\mathrm{}}{\mathit{\rho}}_{\mathrm{}}\mathit{\u03f5}\mathrm{\left(}{\mathit{T}}_{\mathrm{}}\mathrm{\right)}\mathit{.}$(18)We can now determine the critical temperature T_{cr} at which T_{−} equals T_{2}$\mathit{\sigma}{\mathit{T}}_{\mathrm{cr}}^{\mathrm{4}}\mathrm{=}{\mathit{u}}_{\mathrm{}}{\mathit{\rho}}_{\mathrm{}}\mathit{\u03f5}\mathrm{\left(}{\mathit{T}}_{\mathrm{cr}}\mathrm{\right)}\mathit{,}$(19)so that ${\mathit{T}}_{\mathrm{cr}}\mathrm{=}{\left(\frac{{\mathit{u}}_{\mathrm{}}{\mathit{\rho}}_{\mathrm{}}{\mathit{k}}_{\mathrm{B}}}{\mathrm{(}\mathit{\gamma}\mathrm{}\mathrm{1}\mathrm{)}\mathit{\mu}{\mathit{m}}_{\mathrm{H}}\mathit{\sigma}}\right)}^{\mathrm{1}\mathit{/}\mathrm{3}}\mathit{,}$(20)which defines the supercritical shock condition.
2.3.2. Optically thin upstream material
For an optically thin upstream material, the aforederived criterion for a supercritical shock is simply derived by assuming that all the upstream kinetic energy is radiated away at the shock, which yields ${\mathit{T}}_{\mathrm{cr}}\mathrm{=}{\left(\frac{\mathrm{0.5}{\mathit{\rho}}_{\mathrm{}}{\mathit{u}}_{\mathrm{}}^{\mathrm{3}}}{\mathit{\sigma}}\right)}^{\mathrm{1}\mathit{/}\mathrm{4}}\mathit{.}$(21)
2.4. Estimate of the preshock temperature
In this section, we calculate the preshock temperature as a function of the upstream quantities, for the various natures of the shock.
2.4.1. Supercritical shock with an optically thin or thick upstream material
This is the simplest case since the results are independent of the optical depth of the upstream material. Once the upstream density is known, the velocity and the temperature are set. When the upstream gas is not compressed (u_{1} = u_{−}) and when there is a strong shock (ℳ > 2 so that X ~ 1), it is easy to get the shock temperature T_{s} = T_{−} = T_{2} since all the upstream kinetic energy is radiated away: ${\mathit{T}}_{\mathrm{s}}\mathrm{=}{\left(\frac{\mathrm{0.5}{\mathit{\rho}}_{\mathrm{1}}{\mathit{u}}_{\mathrm{1}}^{\mathrm{3}}}{\mathit{\sigma}}\right)}^{\mathrm{1}\mathit{/}\mathrm{4}}\mathit{.}$(22)
2.4.2. Subcritical shock with an optically thick upstream material
For a subcritical shock with an optically thick upstream material, the shock temperature is fixed by the upstream velocity. The equilibrium postshock temperature is (Mihalas & Mihalas 1984) ${\mathit{T}}_{\mathrm{2}}\mathrm{=}\frac{\mathrm{2}\mathrm{(}\mathit{\gamma}\mathrm{}\mathrm{1}\mathrm{)}{\mathit{u}}_{\mathrm{}}^{\mathrm{2}}}{\mathrm{\mathcal{R}}\mathrm{(}\mathit{\gamma}\mathrm{+}\mathrm{1}{\mathrm{)}}^{\mathrm{2}}}\mathit{,}$(23)with ℛ = k_{B} /μm_{H}. The preshock temperature T_{−} is estimated as $\frac{\mathit{m\u0307}\mathrm{\mathcal{R}}{\mathit{T}}_{\mathrm{}}}{\mathit{\gamma}\mathrm{}\mathrm{1}}\mathrm{=}\mathrm{2}\mathit{\sigma}{\mathit{T}}_{\mathrm{2}}^{\mathrm{4}}\sqrt{\mathrm{3}}\mathit{,}$(24)which indicates that all the radiative energy going across this location is absorbed and heats up the upstream gas.
2.5. Jump relations for a barotropic gas
A widely used approximation to handle radiative transfer during the first stages of star formation is the barotropic approximation (e.g. Commerçon et al. 2008). For a barotropic material, the temperature and pressure depend only on the density. The total energy is thus not conserved. In this case, the jump relations are simply $\begin{array}{ccc}& & {\mathit{\rho}}_{\mathrm{1}}{\mathit{u}}_{\mathrm{1}}\mathrm{=}{\mathit{\rho}}_{\mathrm{2}}{\mathit{u}}_{\mathrm{2}}\mathrm{\equiv}\mathit{m\u0307}\mathit{,}\\ & & {\mathit{\rho}}_{\mathrm{1}}{\mathit{u}}_{\mathrm{1}}^{\mathrm{2}}\mathrm{+}{\mathit{P}}_{\mathrm{1}}\mathrm{=}{\mathit{\rho}}_{\mathrm{2}}{\mathit{u}}_{\mathrm{2}}^{\mathrm{2}}\mathrm{+}{\mathit{P}}_{\mathrm{2}}\mathit{,}\end{array}$with the barotropic equation of state (EOS) as a closure relation $\mathit{P}\mathrm{\propto}\mathit{\rho}\left[\mathrm{1}\mathrm{+}{\left(\frac{\mathit{\rho}}{{\mathit{\rho}}_{\mathrm{ad}}}\right)}^{\mathit{\gamma}\mathrm{}\mathrm{1}}\right]\mathit{,}$(27)where ρ_{ad} is the critical density at which the gas becomes adiabatic (see Sect. 3.4.1). Using the same adimensional parameters as previously, we get $\mathit{\gamma}{\mathrm{\mathcal{M}}}_{\mathrm{1}}^{\mathrm{2}}\left(\frac{\mathit{r}\mathrm{}\mathrm{1}}{\mathit{r}}\right)\mathrm{=}\mathrm{\Pi}\mathrm{}\mathrm{1.}$(28)Using the barotropic EOS yields $\mathrm{\Pi}\mathrm{=}\frac{{\mathit{P}}_{\mathrm{2}}}{{\mathit{P}}_{\mathrm{1}}}\mathrm{=}\frac{{\mathit{\rho}}_{\mathrm{2}}}{{\mathit{\rho}}_{\mathrm{1}}}\frac{\mathrm{1}\mathrm{+}{\mathit{r}}_{\mathrm{2}}^{\mathrm{(}\mathit{\gamma}\mathrm{}\mathrm{1}\mathrm{)}}}{\mathrm{1}\mathrm{+}{\mathit{r}}_{\mathrm{1}}^{\mathrm{(}\mathit{\gamma}\mathrm{}\mathrm{1}\mathrm{)}}}\mathrm{=}\mathit{r}\frac{\mathrm{1}\mathrm{+}{\mathit{r}}_{\mathrm{1}}^{\mathrm{(}\mathit{\gamma}\mathrm{}\mathrm{1}\mathrm{)}}{\mathit{r}}^{\mathrm{(}\mathit{\gamma}\mathrm{}\mathrm{1}\mathrm{)}}}{\mathrm{1}\mathrm{+}{\mathit{r}}_{\mathrm{1}}^{\mathrm{(}\mathit{\gamma}\mathrm{}\mathrm{1}\mathrm{)}}}\mathit{,}$(29)where r_{1} = ρ_{1} /ρ_{ad} and r_{2} = ρ_{2} /ρ_{ad}. Eventually, we get r by solving $\mathrm{1}\mathrm{+}{\mathit{r}}_{\mathrm{1}}^{\mathrm{(}\mathit{\gamma}\mathrm{}\mathrm{1}\mathrm{)}}\mathit{\gamma}{\mathrm{\mathcal{M}}}_{\mathrm{1}}^{\mathrm{2}}\frac{\mathit{r}\mathrm{}\mathrm{1}}{\mathit{r}}\mathrm{=}\mathrm{(}\mathit{r}\mathrm{}\mathrm{1}\mathrm{)}\mathrm{+}{\mathit{r}}_{\mathrm{1}}^{\mathit{\gamma}\mathrm{}\mathrm{1}}\mathrm{(}{\mathit{r}}^{\mathit{\gamma}}\mathrm{}\mathrm{1}\mathrm{)}\mathit{.}$(30)One can also derive an equivalent luminosity, assuming a shock in a perfect gas with the same jump properties. From Eq. (9), we have $\mathit{X}\mathrm{=}\frac{\mathrm{\Delta}{\mathit{F}}_{\mathrm{r}}}{\mathrm{0.5}{\mathit{\rho}}_{\mathrm{1}}{\mathit{u}}_{\mathrm{1}}^{\mathrm{3}}}\mathrm{=}\frac{\mathrm{2}}{\mathit{\gamma}\mathrm{}\mathrm{1}}\frac{\mathrm{1}}{{\mathrm{\mathcal{M}}}_{\mathrm{1}}^{\mathrm{2}}}\left(\mathrm{1}\mathrm{}\frac{\mathrm{\Pi}}{\mathit{r}}\right)\mathrm{+}\left(\mathrm{1}\mathrm{}\frac{\mathrm{1}}{{\mathit{r}}^{\mathrm{2}}}\right)\mathit{,}$(31)where r and Π are given by (30).
Fig. 3 Evolution of the compression ratio r and of the amount of kinetic energy radiated away, X, as a function of the upstream Mach number for a barotropic material, with density ratios r_{1} = 0.1 (black), 1 (red), and 10 (blue). The dotted lines represent the evolution for a supercritical shock with an optically thin upstream material (cf. Fig. 2). 
Figure 3 illustrates the evolution of the adimensional quantities r and X as a function of the upstream Mach number for a barotropic material, for density ratios r_{1} = 0.1, 1, and 10, which represent transition values between the optically thin (isothermal) and thick regimes (adiabatic). For comparison, we again plot the evolution of r and X for a supercritical shock with an optically thin upstream material. The compression ratio is not bound by a limiting value as the upstream Mach number increases, as for the case of a supercritical shock with an optically thin upstream material, but in contrast to the case of an opaque material. The compression ratio decreases as r_{1} increases (note that for r_{1} = 10, the material should be totally optically thick in the context of star formation). The shock is never found to be supercritical (X ~ 1), even though almost all the incident kinetic energy is assumed to be radiated away at a high upstream Mach number. This clearly shows that a barotropic EOS approximation cannot treat radiative shocks properly. Compared to the case of a supercritical shock with an optically thin upstream that is relevant for the accretion shock on the first Larson core in the transition regions between the isothermal and adiabatic regimes, the compression ratio and the amount of energy radiated away are underestimated.
The specific entropy jump at the shock is expressed as $\mathrm{\Delta}\mathit{s}\mathrm{=}{\mathit{s}}_{\mathrm{2}}\mathrm{}{\mathit{s}}_{\mathrm{1}}\mathrm{=}{\mathit{C}}_{\mathrm{v}}\left[\mathrm{ln}\left(\frac{{\mathit{P}}_{\mathrm{2}}}{{\mathit{\rho}}_{\mathrm{2}}^{\mathit{\gamma}}}\right)\mathrm{}\mathrm{ln}\left(\frac{{\mathit{P}}_{\mathrm{1}}}{{\mathit{\rho}}_{\mathrm{1}}^{\mathit{\gamma}}}\right)\right]\mathit{,}$(32)where s is the specific entropy and C_{v} = k_{B} /(μm_{H}(γ − 1)). The entropy jump for a barotropic fluid is thus proportional to Π /r^{γ}. Figure 4 shows the evolution of log(Π /r^{γ}) for a barotropic fluid (with r_{1} = 0.1, 1, and 10) and for a shock with an optically thin upstream material. The entropy jump at the shock is overestimated with a barotropic law. This would lead to an incorrect initial entropy level and profile for newborn protostars. The entropy is indeed a key quantity for premain sequence evolutionary models as it entirely determines the massradius relationship of an adiabatic object, like the first and second Larson cores.
Fig. 4 Evolution of the entropy jump, Δs ∝ Π /r^{γ}, as a function of the the upstream Mach number for a barotropic material, with density ratio r_{1} = 0.1 (black), 1 (red), and 10 (blue). The dotted line represents the evolution for a supercritical shock with an optically thin upstream material. 
In this section, we have shown that the downstream quantities of radiative shocks strongly depend on the nature of the shock and on the model used to describe radiation transport. In the following sections, we study the case of a particular radiative shock in detail: the accretion shock on the first Larson core, in the context of star formation.
3. The 1D spherical numerical calculations – Method
3.1. Introduction and previous work
In this section, we introduce our numerical method and basics concepts, for a good understanding of the following part of this work. The two companion papers Masunaga et al. (1998), and Masunaga & Inutsuka (2000) present a very extensive study of 1D protostellar collapse. The first paper focusses on the formation and the properties of the first core, using a frequencydependent RHD model, while the second paper addresses the formation of the second core. These authors use a 1D spherical code in the Lagrangean frame. Moment equations of radiation in the comoving frame are solved following Stone et al. (1992). Masunaga et al. (1998) use a variable Eddington tensor factor (VETF) that retains a frequency dependence, whereas grey opacities are used to calculate the coupling between matter and radiation and the work of the radiative pressure.
In Masunaga et al. (1998), the cores initially have uniform density and temperature, and a radius adjusted so that the cloud should be initially slightly more massive than the Jeans mass. Initial masses range from 0.1 M_{⊙} to 3 M_{⊙}. They find that, regardless of the initial cloud mass, the first core radius R_{fc} is almost constant with R_{fc} ~ 5 AU, and the first core mass is M_{fc} ~ 0.05 M_{⊙}. In their study, R_{fc} is defined as the point where the gas pressure is balanced by the ram pressure of the infalling envelope.
To compare our numerical results with analytical theories, it is useful to define some measurable quantities. A useful one is the mass accretion rate Ṁ. For a 1D spherical model, the mass accretion rate is simply Ṁ = 4πr^{2}ρu. In the theory, the accretion rate is generally defined as (e.g. Stahler et al. 1980) $\mathit{M\u0307}\mathrm{=}\mathit{\xi}\frac{{\mathit{c}}_{\mathit{s}\mathrm{0}}^{\mathrm{3}}}{\mathit{G}}\mathit{,}$(33)where ξ is a dimensionless coefficient and c_{s0} the isothermal sound speed. The value of ξ is estimated assuming that a nonmagnetic, nonrotating cloud, whose initial state corresponds approximately to a balance between thermal support and selfgravity, has comparable freefall time, ${\mathit{t}}_{\mathrm{ff}}\mathrm{~}{\mathit{R}}_{\mathrm{c}}^{\mathrm{3}\mathit{/}\mathrm{2}}\mathit{/}\mathrm{(}\mathit{G}{\mathit{M}}_{\mathrm{c}}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}$, and soundcrossing time, t_{sc} ~ R_{c} /c_{s0}. The mass accretion rate is thus given by Ṁ ~ M_{c} /t_{ff}. Shu (1977) obtained ξ = 0.975 for the expansionwave solution, whereas the dynamic LarsonPenston solution yields ξ ~ 50 (Larson 1969; Penston 1969).
Another quantity directly comparable to the theory is the accretion luminosity, usually expressed as ${\mathit{L}}_{\mathrm{fc}}^{\mathrm{acc}}\mathrm{=}\frac{\mathit{G}{\mathit{M}}_{\mathrm{fc}}\mathit{M\u0307}}{{\mathit{R}}_{\mathrm{fc}}}\mathit{,}$(34)where M_{fc} corresponds to the mass of the accreting core. This accretion luminosity thus corresponds to the case where all the infalling kinetic energy is radiated away, i.e. to a supercritical radiative shock. The accretion luminosity can be directly measured in numerical calculations from the mass accretion rate, and compared with the intrinsic effective luminosity emerging from the first core.
3.2. Numerical method
We use a 1D full Lagrangean version of the code developed by Chièze and Audit (see Audit et al. 2002), which integrates the equations of grey radiation hydrodynamics under three different assumptions (in order of increasing complexity order): a barotropic EOS, the fluxlimited diffusion approximation (FLD), and the M1 model. The RHD equations are integrated in their nonconservative forms using finite volumes and an artificial viscosity scheme in tensorial form. The RHD equations are integrated with an implicit scheme in time, using a standard RaphsonNewton iterative method.
3.3. General assumptions for the radiation field
The first approximation in our calculations is to consider a grey radiation field, i.e. only one group of photons integrated over all frequencies. We also neglect scattering and assume local thermodynamical equilibrium (LTE) everywhere. The second assumption concerns the coupling between radiation and hydrodynamics, for which we write the RHD equations in the comoving frame.
3.4. Models for the radiative transfer
According to previous studies using an accurate model for the radiation field (e.g. Masunaga et al. 1998), it is well established that the molecular gas follows two thermal regimes during the first collapse. At low density, the gas is able to radiate freely and to couple with the dust. This is the isothermal regime, where the Jeans mass decreases with density. Then, when the gas becomes denser, typically ρ > ρ_{ad} ~ 1 × 10^{13} g cm^{3}, the radiation is trapped within the gas that begins to heat up adiabatically. The Jeans mass then increases with density until the gas becomes hot enough to dissociate H_{2}, leading to the onset of the second collapse.
3.4.1. The barotropic equation of state approximation
As mentioned earlier, the easiest way to describe the thermal evolution of the gas without solving the radiative transfer equation is to adopt a barotropic EOS that reproduces both the isothermal and adiabatic limiting regimes as a function of the density. We use the following barotropic EOS: $\frac{\mathit{P}}{\mathit{\rho}}\mathrm{=}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathrm{=}{\mathit{c}}_{\mathit{s}\mathrm{0}}^{\mathrm{2}}\left[\mathrm{1}\mathrm{+}{\left(\frac{\mathit{\rho}}{{\mathit{\rho}}_{\mathrm{ad}}}\right)}^{\mathrm{2}\mathit{/}\mathrm{3}}\right]\mathit{,}$(35)where ρ_{ad} is the critical density at which the gas becomes adiabatic. This critical density is obtained by more accurate calculations and depends on the opacities, the composition, and the geometry of the molecular gas. At low densities, ρ ≪ ρ_{ad}, c_{s} ~ c_{s0} = 0.19 km s^{1} (for T ~ 10 K), and the molecular gas is able to radiate freely by thermally coupling to the dust and remains isothermal at 10 K. At high densities, ρ > ρ_{ad}, and we assume that the cooling due to radiative transfer is trapped by the dust opacity. Therefore, P ∝ ρ^{5 /3}, which corresponds to an adiabatic, monoatomic gas with adiabatic exponent γ = 5 /3. Molecular hydrogen behaves like a monoatomic gas until the temperature reaches several hundred Kelvin since the rotational degrees of freedom are not excited at lower temperatures (Whitworth & Clarke 1997; Masunaga & Inutsuka 2000).
3.4.2. Moment models
In RHD calculations, the radiative transfer equation should be formally integrated over six dimensions at each time step. This process is too computationally demanding for multidimensional numerical analysis. Using angular moments of the transfer equation is thus very useful by allowing a large reduction of the computational cost. However, each evolution equation of a moment of the transfer equation involves the next higher order moment of the intensity. Consequently, as for the kinetic theory of gases, the system must be closed by using an ad hoc relation that gives the highest moment as a function of the lower order moments. For radiation transport, the closure theory is usually limited to the two first moments of the transfer equation. A closure relation for the system is needed, and this relation is prime. Many possible choices for the closure relation exist. In the following, we present two models based on more or less accurate closure relations, which we use in this work: the fluxlimited diffusion (FLD) approximation and the M1 model.
3.4.3. The fluxlimited diffusion approximation
The diffusion approximation is the most widely used moment model of radiation transport. The diffusion limit is valid when the photon meanfree path is short compared with other length scales in the system. On the contrary, the approximation is no longer accurate in the transport regime. In the diffusion limit, photons diffuse through the material in a random walk. Readers can find an accurate derivation of the diffusion limit in Mihalas & Mihalas (1984), Sect. 80.
In the diffusion limit, the radiative energy E_{r} and radiative flux F_{r} are related simply by ${F}\mathrm{r}\mathrm{=}\mathrm{}\frac{\mathit{c}}{\mathrm{3}{\mathit{\kappa}}_{\mathrm{R}}}\mathrm{\nabla}{\mathit{E}}_{\mathrm{r}}\mathit{,}$(36)where κ_{R} is the Rosseland mean opacity. The radiative flux is expressed directly as a function of the radiative energy and is proportional and colinear to the radiative energy gradient. Equation (36) has no upper limit, but for optically thin flows, the effective propagation speed of the radiation must be limited to c (F_{r} ≤ cE_{r}). We thus have to limit the propagation speed of the radiation by means of a flux limiter. Equation (36) is then expressed as ${F}\mathrm{r}\mathrm{=}\mathrm{}\frac{\mathit{c\lambda}}{{\mathit{\kappa}}_{\mathrm{R}}}\mathrm{\nabla}{\mathit{E}}_{\mathrm{r}}\mathit{,}$(37)where λ = λ(R) is the flux limiter.
In this study, we retain the flux limiter derived by Minerbo (1978)$\mathit{\lambda}\mathrm{=}\{\begin{array}{c}\\ \mathrm{2}\mathit{/}\mathrm{(}\mathrm{3}\mathrm{+}\sqrt{\mathrm{9}\mathrm{+}\mathrm{12}{\mathit{R}}^{\mathrm{2}}}\mathrm{)}& \mathrm{if}& \mathrm{0}\mathrm{\le}\mathit{R}\mathrm{\le}\mathrm{3}\mathit{/}\mathrm{2}\mathit{,}\\ \mathrm{(}\mathrm{1}\mathrm{+}\mathit{R}\mathrm{+}\sqrt{\mathrm{1}\mathrm{+}\mathrm{2}\mathit{R}}{\mathrm{)}}^{1}& \mathrm{if}& \mathrm{3}\mathit{/}\mathrm{2}\mathit{<}\mathit{R}\mathrm{\le}\mathrm{\infty}\mathit{,}\end{array}$(38)with R = ∇E_{r} /(κ_{R}E_{r}). The flux limiter has the property that λ → 1 /3 in optically thick regions and λ → 1 /R in optically thin regions.
In the FLD approximation, a unique diffusiontype equation on the radiative energy is thus obtained $\frac{\mathit{\partial}{\mathit{E}}_{\mathrm{r}}}{\mathit{\partial t}}\mathrm{}\mathrm{\nabla}\mathrm{\xb7}\left(\frac{\mathit{c\lambda}}{{\mathit{\kappa}}_{\mathrm{R}}}\mathrm{\nabla}{\mathit{E}}_{\mathrm{r}}\right)\mathrm{=}{\mathit{\kappa}}_{\mathrm{P}}\mathrm{(}\mathrm{4}\mathit{\pi B}\mathrm{}\mathit{c}{\mathit{E}}_{\mathrm{r}}\mathrm{)}\mathit{.}$(39)
3.4.4. M1 model
In the M1 model, the radiation transport is described by the first two moments of the radiative transfer equation (radiative energy and flux). We use a closure relation introduced by Dubroca & Feugeas (1999). The M1 method is used in the RHD code HERACLES (González et al. 2007). Based on a minimum entropy principle, which is a method able to account for the large anisotropy of the radiation, as well as for the correct diffusion limit. The main advantage of the M1 system is that the underlying photon distribution function is not isotropic, but has a prefered direction of propagation. The M1 model also allows explicitly getting rid of the adhoc limitation of the flux in the transport regime.
We now consider the evolution equations of the zeroth and first moments of the specific intensity in the laboratory frame $\{\begin{array}{c}\\ {\mathit{\partial}}_{\mathit{t}}{\mathit{E}}_{\mathrm{r}}\mathrm{+}\mathrm{\nabla}\mathrm{\xb7}{F}\mathrm{r}& \mathrm{=}& {\mathit{\kappa}}_{\mathrm{P}}\mathrm{(}\mathrm{4}\mathit{\pi B}\mathrm{}\mathit{c}{\mathit{E}}_{\mathrm{r}}\mathrm{)}\\ {\mathit{\partial}}_{\mathit{t}}{F}\mathrm{r}\mathrm{+}{\mathit{c}}^{\mathrm{2}}\mathrm{\nabla}\mathrm{\xb7}{P}_{\mathrm{r}}& \mathrm{=}& \mathrm{}{\mathit{\kappa}}_{\mathrm{R}}\mathit{c}{F}\mathrm{r}\mathit{,}\end{array}$(40)where κ_{P} is thePlanck mean opacity. We use the Rosseland mean in the first moment equation in order to yield the diffusion limit. As a closure relation, the radiative pressure is often expressed as P_{r} = DE_{r}, where D is the Eddington tensor. Assuming that the direction of the radiative flux is an axis of symmetry of the local specific intensity, the Eddington tensor is given by (Levermore 1984) $D\mathrm{=}\frac{\mathrm{1}\mathrm{}\mathit{\chi}}{\mathrm{2}}I\mathrm{+}\frac{\mathrm{3}\mathit{\chi}\mathrm{}\mathrm{1}}{\mathrm{2}}{n}\mathrm{\otimes}{n}\mathit{,}$(41)where χ is the Eddington factor, I the identity matrix, and ${n}\mathrm{=}\frac{{f}}{\mathit{f}}$ a unit vector aligned with the radiative flux. In the M1 approximation, the Eddington factor χ is then obtained by minimizing the radiative entropy (Dubroca & Feugeas 1999). The Eddington factor is then found to be $\mathit{\chi}\mathrm{=}\frac{\mathrm{3}\mathrm{+}\mathrm{4}{\mathit{f}}^{\mathrm{2}}}{\mathrm{5}\mathrm{+}\mathrm{2}\sqrt{\mathrm{4}\mathrm{}\mathrm{3}{\mathit{f}}^{\mathrm{2}}}}\mathrm{\xb7}$(42)It is then trivial to recover the two asymptotic regimes of radiative transfer. When f → 0, χ = 1 /3, and P_{r} = 1 /3E_{r}I, which corresponds to the diffusion limit with an isotropic radiative pressure. On the other hand, if f = 1, χ = 1, and P_{r} = n ⊗ nE_{r}, which corresponds to the freestreaming limit. Between these two limits, the closure relation ensures that energy remains positive and that the flux is limited (F_{r} < cE_{r}).
3.4.5. Systems of RHD equations

In the barotropic EOS approximation, the thermal behaviour of the gas is determined by the choice of the EOS. The energy equation becomes superfluous, and the Euler equations reduce to $\{\begin{array}{c}\\ {\mathit{\partial}}_{\mathit{t}}\mathit{\rho}\mathrm{+}\mathrm{\nabla}\left[\mathit{\rho}{u}\right]\mathrm{=}\mathrm{0}\\ {\mathit{\partial}}_{\mathit{t}}\mathit{\rho}{u}\mathrm{+}\mathrm{\nabla}\left[\mathit{\rho}{u}\mathrm{\otimes}{u}\mathrm{+}\mathit{P}I\right]\mathrm{=}\mathrm{}\mathit{\rho}\mathrm{\nabla \Phi}\mathit{,}\end{array}$(43)where Φ the gravitational potential.

In the FLD approximation, the system of RHD equations corresponds to the Euler equations plus the equation on the radiative energy $\{\begin{array}{c}\\ {\mathit{\partial}}_{\mathit{t}}\mathit{\rho}\mathrm{+}\mathrm{\nabla}\left[\mathit{\rho}{u}\right]\mathrm{=}\mathrm{0}\\ {\mathit{\partial}}_{\mathit{t}}\mathit{\rho}{u}\mathrm{+}\mathrm{\nabla}\left[\mathit{\rho}{u}\mathrm{\otimes}{u}\mathrm{+}\mathit{P}I\right]\mathrm{=}\mathrm{}\mathit{\rho}\mathrm{\nabla \Phi}\mathrm{+}\mathit{\lambda}\mathrm{\nabla}{\mathit{E}}_{\mathrm{r}}\\ {\mathit{\partial}}_{\mathit{t}}\mathit{E}\mathrm{+}\mathrm{\nabla}\left[{u}\mathrm{(}\mathit{E}\mathrm{+}{\mathit{P}}^{\mathrm{)}}\right]\mathrm{=}\mathrm{}\mathit{\rho}{u}\mathrm{\xb7}\mathrm{\nabla \Phi}\mathrm{+}\mathit{\lambda}\mathrm{\nabla}{\mathit{E}}_{\mathrm{r}}\mathrm{\xb7}{u}\\ \mathrm{}{\mathit{\kappa}}_{\mathrm{P}}\mathrm{(}\mathrm{4}\mathit{\pi B}\mathrm{}\mathit{c}{\mathit{E}}_{\mathrm{r}}\mathrm{)}\\ {\mathit{\partial}}_{\mathit{t}}{\mathit{E}}_{\mathrm{r}}\mathrm{+}\mathrm{\nabla}\left[{u}{\mathit{E}}_{\mathrm{r}}\right]\mathrm{=}\mathrm{}{P}_{\mathrm{r}}\mathrm{:\nabla}{u}\mathrm{+}\mathrm{\nabla}\mathrm{\xb7}\mathrm{(}\frac{\mathit{c\lambda}}{{\mathit{\kappa}}_{\mathrm{R}}}\mathrm{\nabla}{\mathit{E}}_{\mathrm{r}}\mathrm{)}\\ \mathrm{+}{\mathit{\kappa}}_{\mathrm{P}}\mathrm{(}\mathrm{4}\mathit{\pi B}\mathrm{}\mathit{c}{\mathit{E}}_{\mathrm{r}}\mathrm{)}\mathit{.}\end{array}$(44)This system is closed by the perfect gas relation and the simplest FLD approximation for the radiative pressure, P_{r} = λE_{r}.
Fig. 5 Rosseland opacity made of a mix of Semenov et al. (2003) model at low temperature and Ferguson et al. (2005) model at high temperature. The Rosseland mean opacity is plotted as a function of temperature (xaxis) and $\mathit{R}\mathrm{=}\mathit{\rho}\mathit{/}\mathrm{\left(}{\mathit{T}}_{\mathrm{6}}^{\mathrm{3}}\mathrm{\right)}$, with T_{6} = T /10^{6}, (yaxis) using logarithm scales.

In the M1 model, we consider the equations governing the evolution of an inviscid, radiating fluid $\{\begin{array}{c}\\ {\mathit{\partial}}_{\mathit{t}}\mathit{\rho}\mathrm{+}\mathrm{\nabla}\left[\mathit{\rho}{u}\right]\mathrm{=}\mathrm{0}\\ {\mathit{\partial}}_{\mathit{t}}\mathit{\rho}{u}\mathrm{+}\mathrm{\nabla}\left[\mathit{\rho}{u}\mathrm{\otimes}{u}\mathrm{+}\mathit{P}I\right]\mathrm{=}\mathrm{}\mathit{\rho}\mathrm{\nabla \Phi}\mathrm{+}{\mathit{\kappa}}_{\mathrm{R}}{F}\mathrm{r}\mathit{/}\mathit{c}\\ {\mathit{\partial}}_{\mathit{t}}\mathit{E}\mathrm{+}\mathrm{\nabla}\left[{u}\mathrm{(}\mathit{E}\mathrm{+}{\mathit{P}}^{\mathrm{)}}\right]\mathrm{=}\mathrm{}\mathit{\rho}{u}\mathrm{\xb7}\mathrm{\nabla \Phi}\mathrm{+}{\mathit{\kappa}}_{\mathrm{R}}{F}\mathrm{r}\mathit{/}\mathit{c}\mathrm{\xb7}{u}\\ \mathrm{}{\mathit{\kappa}}_{\mathrm{P}}\mathrm{(}\mathrm{4}\mathit{\pi B}\mathrm{}\mathit{c}{\mathit{E}}_{\mathrm{r}}\mathrm{)}\\ {\mathit{\partial}}_{\mathit{t}}{\mathit{E}}_{\mathrm{r}}\mathrm{+}\mathrm{\nabla}\left[{u}{\mathit{E}}_{\mathrm{r}}\right]\mathrm{+}\mathrm{\nabla}\mathrm{\xb7}{F}\mathrm{r}\mathrm{=}\mathrm{}{P}_{\mathrm{r}}\mathrm{:\nabla}{u}\mathrm{+}{\mathit{\kappa}}_{\mathrm{P}}\mathrm{(}\mathrm{4}\mathit{\pi B}\mathrm{}\mathit{c}{\mathit{E}}_{\mathrm{r}}\mathrm{)}\\ {\mathit{\partial}}_{\mathit{t}}{F}\mathrm{r}\mathrm{+}\mathrm{\nabla}\left[{u}{F}\mathrm{r}\right]\mathrm{+}{\mathit{c}}^{\mathrm{2}}\mathrm{\nabla}\mathrm{\xb7}{P}_{\mathrm{r}}\mathrm{=}\mathrm{}\left({F}\mathrm{r}\mathrm{\xb7}\mathrm{\nabla}\right){u}\mathrm{}{\mathit{\kappa}}_{\mathrm{R}}\mathit{c}{F}\mathrm{r}\mathit{,}\end{array}$(45)where ρ is the material density, u the velocity, P the isotropic thermal pressure, Φ the gravitational potential, and E the fluid total energy E = ρϵ + 1 /2ρu^{2}. This system is closed by the perfect gas relation and the M1 relation (42).
3.5. Initial and boundary conditions
We use the same model as Masunaga et al. (1998), i.e. a uniform density sphere of mass M_{0} = 1 M_{⊙}, temperature T_{0} = 10 K (c_{s0} = 0.187 km s^{1}), and radius R_{0} = 10^{4} AU. This initial setup corresponds to a ratio α of thermal to gravitational energies of α = 5R_{0}k_{B}T_{0} /(2GM_{0}μm_{H}) ~ 0.97 and to a freefall time t_{ff} ~ 1.77 × 10^{5} yr. Boundary conditions are very simple: for hydrodynamics, we impose a constant thermal pressure equal to the initial pressure (other quantities are free); and we impose a vanishing gradient on radiative temperature for the radiation field. Calculations were performed using a Lagrangean grid containing 4500 cells.
Summary of first core properties for ρ_{c} = 1 × 10^{10} g cm^{3}.
3.6. The opacities
For the M1 model and the FLD approximation, we use the set of opacities given by Semenov et al. (2003) for low temperature (<1000 K) and Ferguson et al. (2005) for high temperature (>1000 K), which we compute as a function of the gas temperature and density. In Fig. 5, the table for the Rosseland opacity is plotted as a function of temperature and density. In Semenov et al. (2003), the dependence of the evaporation temperatures of ice, silicates, and iron on gas density are taken into account. In this work, we use spherical composite aggregate particles for the grain structure and topology and a normal iron content in the silicates, Fe/(Fe + Mg) = 0.3.
4. The 1D numerical calculations of the first collapse
As mentioned before, all calculations presented here were performed with 4500 cells, distributed according to a logarithmic scale in mass in our initial setup. These mass and spatial resolutions are high enough to resolve quantitatively the accretion shock energy budget (see Appendix A). Calculations with the barotropic EOS were performed using a critical density of ρ_{ad} = 3.7 × 10^{13} g cm^{3}. We derived the value of ρ_{ad} from the intersection of extrapolated lines of the isotherm and the adiabat in the log(ρ)log(T) plane, obtained with calculations using the M1 model (see Fig. 6).
Fig. 6 Sketch of the estimate of ρ_{ad} used in the barotropic EOS, according to the densitytemperature distribution obtained with the M1 model. 
4.1. Results during the first collapse. Formation of the first core
Table 1 gives results of barotropic EOS (Baro), M1, and FLD calculations, when the central density reaches ρ_{c} = 1 × 10^{10} g cm^{3}. For all calculations, the dynamical times are very close to ~0.189 Myr ~1.07 t_{ff}. We defined the first core radius as the radius at which the infall velocity is the highest (i.e., the position of the accretion shock). The accretion luminosity was estimated at the first core border, according to equation (34). The mass accretion rate Ṁ and the accretion parameter ξ were evaluated at R_{fc}. Table 1 also displays the central temperature T_{c}, the central specific entropy s_{c}, and the first core temperature at the border T_{fc}.
As seen from Table 1, the first core radius and mass agree with Masunaga et al. (1998) results. The results from the M1 or FLD calculations are very similar. The central temperature with the FLD is slightly higher (by ~2%) than the one in the M1 calculations, which confirms that the diffusion approximation tends to slightly overestimate the cooling. The central entropy obtained with the M1 and FLD models is lower than the one obtained with the barotropic EOS, and this difference increases with time. Indeed, taking radiative transfer into account allows the energy to be radiated away during the collapse, a process which is not properly handled with a barotropic EOS. This leads to a lower entropy level of the first core with the FLD or the M1 models. The temperature at the border of the first core is higher in the M1 and FLD cases, since the photons escaping the accretion shock heat up the infalling material. On the other hand, the mass accretion rate, the mass, the radius, and the accretion parameter ξ of the first core are in good agreement between the three models. The values of ξ obtained in our calculations are closer to the LarsonPenston solution than to Shu’s expansion wave solution.
Fig. 7 Evolution of the central temperature (left) and the central entropy (right) as a function of the central density during the first collapse for the M1 (solid line), FLD (dashed red line) and barotropic (dasheddotted line). 
Figure 7 shows the evolution of the central temperature and the central entropy as a function of the central density. From Fig. 7a, we notice the perfect bimodal thermal behaviour of the gas from isothermal to adiabatic. The critical density at which the gas begins to heat in both the M1 and FLD calculations is the same. The difference with the barotropic case stems from our choice of the barotropic EOS and critical density. For FLD and M1, we find a slope of γ_{eff} − 1 ~ 0.61 < 2 /3 at high density. This means that the first core is not completely adiabatic, but it experiences some heat loss, yielding a significant cooling.
In Fig. 7b, the central entropy is plotted as a function of the central density. At high density, we recover the adiabatic regime, and the M1 and FLD calculations settle at the same entropy level in the centre, which is lower than the one reached by the barotropic model. In the barotropic case, the entropy is determined by the EOS. The value of the “adiabat” at which the gas settles is $\mathit{s}\mathrm{\propto}\mathrm{ln}\left(\frac{\mathit{P}}{{\mathit{\rho}}^{\mathit{\gamma}}}\right)\mathrm{\xb7}$(46)For the barotropic EOS, we have $\mathit{s}\mathrm{\propto}\mathrm{ln}\left[{\mathit{\rho}}^{\mathrm{}\mathrm{2}\mathit{/}\mathrm{3}}{\mathit{c}}_{\mathit{s}\mathrm{0}}^{\mathrm{2}}\left(\mathrm{1}\mathrm{+}{\left(\frac{\mathit{\rho}}{{\mathit{\rho}}_{\mathrm{ad}}}\right)}^{\mathrm{2}\mathit{/}\mathrm{3}}\right)\right]\mathrm{\propto}\mathrm{ln}\left(\frac{{\mathit{c}}_{\mathit{s}\mathrm{0}}^{\mathrm{2}}}{{\mathit{\rho}}_{\mathrm{ad}}^{\mathrm{2}\mathit{/}\mathrm{3}}}\mathrm{+}\frac{{\mathit{c}}_{\mathit{s}\mathrm{0}}^{\mathrm{2}}}{{\mathit{\rho}}^{\mathrm{2}\mathit{/}\mathrm{3}}}\right)\mathrm{\xb7}$(47)At high density, $\mathit{s}\mathrm{\to}\mathrm{ln}\mathrm{(}{\mathit{c}}_{\mathit{s}\mathrm{0}}^{\mathrm{2}}\mathit{/}{\mathit{\rho}}_{\mathrm{ad}}^{\mathrm{2}\mathit{/}\mathrm{3}}\mathrm{)}$. The limiting value for the entropy is then ~2.03 × 10^{9} erg K^{1} g^{1} for the barotropic EOS, which is higher than the value obtained with FLD and M1. Moreover, as the slope of the thermal profile (Fig. 7a) is not exactly equal to γ = 5 /3 in the M1 and FLD calculations, as mentioned above, the gas tends to cool and to decrease its entropy level at a rate ~10^{3}–10^{4} erg K^{1} g^{1} s^{1} at the centre. This is one of the immediate and most important consequences of correctly taking radiative transfer into account in the collapse: the accretion shock becomes a real radiative shock and radiation is transported outwards in the core.
Fig. 8 Evolution of the central temperature, density, entropy, and optical depth with time for the M1 (solid line), FLD (dashed red line), and barotropic (dasheddotted line) models. 
Fig. 9 Radial profiles of various properties during the collapse of a 1 M_{⊙} dense clump for a core central density ρ_{c} = 1 × 10^{10} g cm^{3}, for the M1 (solid line), FLD (dashed red line) and barotropic (dasheddotted line) models. From top left to bottom right: a) density; b) gas temperature; c) velocity; d) entropy; e) optical depth; f) luminosity; g) radiative flux; and h) integrated mass. 
Figure 8 shows the central temperature, density, entropy, and optical depth evolution during the collapse. Variations in all variables are quite similar for temperature lower than T ~ 100 K. At 100 K, there is a first discontinuity in the opacity, owing to the destruction of icy grains. This discontinuity moves to higher temperatures as density increases, so that density affects the opacity to some extent. However, these differences are really small. From Fig. 8c, we clearly see that the entropy level remains constant with the barotropic model, whereas the first core keeps cooling to lower entropy levels with the M1 and FLD calculations. This is an important result, since, as mentioned previously, the entropy level of a lowmass protostar determines its massradius relationship, hence its premain sequence subsequent evolution. According to the first principle of thermodynamics, the global entropy loss of the first core can be estimated as ${\left(\frac{\mathit{\partial L}}{\mathit{\partial m}}\right)}_{\mathit{t}}\mathrm{=}\mathrm{}\mathit{T}{\left(\frac{\mathit{\partial s}}{\mathit{\partial t}}\right)}_{\mathit{m}}\mathrm{=}\mathrm{}{\left(\frac{\mathit{\partial \u03f5}}{\mathit{\partial t}}\right)}_{\mathit{m}}\mathrm{}\mathit{P}{\left(\frac{\mathit{\partial v}}{\mathit{\partial t}}\right)}_{\mathit{m}}\mathit{,}$(48)where ϵ is the specific internal energy and v = 1 /ρ is the specific volume. If the first core was perfectly adiabatic, the compressional work would be converted entirely into internal energy, and the radiative loss would be zero. As shown in Fig. 9f, the luminosity increases significantly within the first core in the FLD and M1 calculations. This increase corresponds to the radiative loss, which amounts to ~4 × 10^{3} L_{⊙}. The corresponding entropy loss of the first core is then $\frac{\mathrm{\Delta}{\mathit{s}}_{\mathrm{fc}}}{\mathrm{\Delta}\mathit{t}}\mathrm{~}\mathrm{}\frac{\mathrm{1}}{\mathit{T}}\frac{\mathrm{\Delta}\mathit{L}}{{\mathit{M}}_{\mathrm{fc}}}\mathrm{\xb7}$(49)At ρ_{c} = 1 × 10^{10} g cm^{3}, we get M_{fc} ~ 2.3 × 10^{2} M_{⊙}. The entropy loss is thus $\frac{\mathrm{\Delta}{\mathit{s}}_{\mathrm{fc}}}{\mathrm{\Delta}\mathit{t}}\mathrm{~}\mathrm{}\frac{\mathrm{3.3}\mathrm{\times}{\mathrm{10}}^{1}}{\mathit{T}}\mathrm{erg}{\mathrm{K}}^{1}{\mathrm{g}}^{1}\mathit{.}$(50)For a typical first core temperature of a few 100 K, the entropy loss rate is ~10^{3} erg K^{1} g^{1} yr^{1}, which is consistent with the value found at the centre. Because the first core lasted about a few hundred years, the total entropy loss thus represents a few percent of the core’s initial entropy content, and this entropy loss increases with time, as shown in Fig. 1 of Masunaga et al. (1998). Such an entropy loss cannot be handled with a barotropic EOS.
In Fig. 8d, we see that the first core quickly becomes optically thick once it begins to heat up. The optical depth at the centre is so high that observations are not able to catch the central evolution of the cores at this stage of the evolution.
4.2. Mechanical and thermal profile of the first prestellar core
Figure 9 displays the radial profiles of the density, gas temperature, velocity, specific entropy, optical depth, luminosity, radiative flux, and integrated mass, once the central density reaches ρ_{c} = 1 × 10^{10} g cm^{3}, for calculations with M1, FLD, and the barotropic EOS. As mentioned before, the first core radius is smaller with the barotropic model, whereas both M1 and FLD results are similar. For these models, we find only small differences around τ ~ 1; i.e., at the transition between the optically thick and thin regimes. Since the M1 model naturally recovers the diffusion and transport regimes and since the FLD model is defined to recover these limits as well, it is natural that we get similar results when either the transport or the diffusion regimes are well established. Although the accretion shock is located within the transition region around τ ~ 10, the aforementioned small differences do not affect the first core properties. The entropy jump at the accretion shock is much higher with M1 and FLD than with the barotropic approximation. We also see from the temperature and entropy profiles that the barotropic EOS cannot correctly handle the transition from the isothermal to the adiabatic regime, as pointed out earlier. The radiative precursor in front of the shock is not reproduced with the barotropic EOS. In that case, the gas becomes rapidly isothermal after the shock, whereas it is heated by photons escaping from the shock in the FLD and M1 models. As a consequence, the core temperature is higher in the barotropic case. The differences in the behaviour of the radiative flux at small radius come from the differences in temperature between the M1 and FLD models.
The emergent luminosity, on the other hand, is the same for both the M1 and FLD models. The luminosity jump, ~0.013 L_{⊙} is consistent with the accretion luminosity estimated from Eq. (34), 0.014 L_{⊙}. This means that all the infalling kinetic energy is radiated away at the shock boundary; i.e., that the shock is supercritical at the formation of the first core.
If we apply the criterion derived in Eqs. (20) and (21) to the upstream quantities estimated in Fig. 9, i.e. ρ_{1} ~ 7 × 10^{14} g cm^{3}, u_{1} ~ 2 × 10^{5} cm s^{1}, and T_{1} ~ 57 K, the corresponding critical temperatures are T_{cr} ~ 23 K with (20) and T_{cr} ~ 57 K with (21). This confirms that the accretion shock on the first core is supercritical. However, since accretion takes place in the transition region between optically thin and thick regimes and since most of the upstream region is optically thin (see luminosity and optical depth profiles in Fig. 9), the most relevant model is the one we developed for a supercritical shock with a upstream optically thin material.
4.3. Comparison with an analytical model
In Sect. 2, we developed a semianalytical solver that can be applied to the accretion radiative shock on the first core and that covers three cases: the case of either sub or supercritical shocks in an optically thick material and the case of a supercritical shock in an optically thin material. In the following, we develop a simple model for the protostellar collapse, where the upstream quantities are estimated under some basic assumptions, and we compare this model with our numerical results. We need to estimate the density, the velocity, and the temperature in the preshock region in the context of the first core formation. We can easily get the density (and the velocity) from the LarsonPenston and Shu selfsimilar solutions. The tricky part is to estimate the temperature before the shock. However, as mentioned in Sect. 2.4, the preshock temperature is determined by the upstream velocity. Assuming that the accretion shock is supercritical and that the upstream material is optically thin, it is then easy to get this temperature and to recover all fluid variables. Omukai (2007) proposes an alternative model that fully describes the first core characteristics but does not use jump relations for a radiating fluid. In our model, we consider the characteristics at the first core border and the jump conditions for a radiating fluid.
We approximate the upstream velocity by the freefall velocity ${\mathit{u}}_{\mathrm{1}\mathit{,}\mathrm{ff}}\mathrm{=}\sqrt{\frac{\mathrm{2}\mathit{G}{\mathit{M}}_{\mathrm{fc}}}{{\mathit{R}}_{\mathrm{fc}}}}\mathrm{\xb7}$(51)The preshock density is given assuming a profile ∝ r^{2} in the freefalling envelope ${\mathit{\rho}}_{\mathrm{1}}\mathrm{=}\frac{{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}}{\mathrm{2}\mathit{\pi G}}{\mathit{R}}_{\mathrm{fc}}^{2}\mathit{,}$(52)where c_{s} = (P /ρ)^{1 /2} is the isothermal sound speed (at 10 K). The temperature is estimated by assuming a supercritical shock, i.e. all the upstream kinetic energy being radiated away. Moreover, our calculations using the FLD or M1 models yield ℳ_{1} ~ 2 and Fig. 2 shows that X ~ 1. The shock temperature thus reads (cf. Eq. (22)) $\mathit{\sigma}{\mathit{T}}_{\mathrm{s}}^{\mathrm{4}}\mathrm{=}\frac{{\mathit{\rho}}_{\mathrm{1}}{\mathit{u}}_{\mathrm{1}\mathit{,}\mathrm{ff}}^{\mathrm{3}}}{\mathrm{2}}\mathrm{\xb7}$(53)We take the first core properties obtained in our numerical calculations, similar to those in Masunaga et al. (1998), i.e. M_{fc} ~ 2.3 × 10^{2} M_{⊙} and R_{fc} ~ 8 AU. If we apply Eqs. (51)–(53) to these quantities, we get ρ_{1} = 5.8 × 10^{14} g cm^{3}, u_{1} ~ 2.3 × 10^{5} cm s^{1} and T_{1} ~ 50 K.
We now apply these preschock quantities to the jump properties derived in Sect. 2.2.2 for a supercritical shock with an optically thin upstream material. According to Eq. (14), we get r = ρ_{2} /ρ_{1} ~ 50, which gives ρ_{2} ~ 2.9 × 10^{12} g cm^{3}, and u_{2} = 4.6 × 10^{3} cm s^{1}. Using Eq. (15), we get X ~ 0.99: all the infalling kinetic energy is radiated away at the shock.
We can compare these analytical results with the numerical ones given in Fig. 9 and in Table 1: ρ_{1} ~ 7 × 10^{14} g cm^{3}, u_{1} ~ 2 × 10^{5} cm s^{1} and T_{1} = 57 K. Thus, the analytical estimates for the preshock quantities agree quite closely with the numerical values. We also read from the figure and table ρ_{2} ~ 2.15 × 10^{12} g cm^{3}, u_{2} ~ 10^{4} cm s^{1} ≪ u_{1}, again very similar to our analytical estimates. These comparisons validate our simple analytical model fot infering the properties of the first Larson core.
5. Summary and perspectives
In this paper, we have conducted RHD calculations aimed at describing the physical and radiative properties of radiative shocks. We applied these calculations to the formation of the first Larson core, during the early phases of protostellar collapse. We also derived simple analytical models and compared the results with the ones obtained in the simulations. The main results of our study can be summarized as follows.

1.
The properties of the first collapse and the characteristics of the first core in our calculations are in good agreement with the ones obtained by Masunaga et al. (1998). We find that the first core has a typical radius of ~8 AU and a mass ~2 × 10^{2} M_{⊙}. This sets up the initial conditions for the second collapse and the formation of the second Larson core.

2.
We showed that, at the first core stage, the accretion shock is a supercritical radiative shock, at which all the infalling kinetic energy is radiated away, and that a barotropic EOS cannot reproduce the correct jump conditions at the shock. The FLD and M1 calculations show that there is a substantial entropy loss during the formation of the first core, owing to the radiative loss. A barotropic EOS cannot handle this cooling mechanism correctly. In consequence, the first core’s entropy content obtained with a barotropic approximation is overestimated, compared with the calculations that solve the radiative transfer. Such a cooling effect can have a strong impact on the core fragmentation process and, eventually, on the initial properties of the future protostar (Commerçon et al., in prep.).

3.
We confirmed that, when radiative cooling is properly taken into account, the transition from an isothermal to an adiabatic phase during the first collapse and the formation of the first Larson core does not necessarily correspond to an optical depth of unity, as initially shown by Masunaga & Inutsuka (1999).

4.
We developed a simple analytical model for supercritical shocks within an optically thin medium, which reproduces the jump conditions obtained with the numerical calculations well. We show that the compression ratio in such a shock can become very high (r ~ 50). We plan to keep exploring this issue with a frequencydependent radiative transfer model. Indeed, strong shocks on (massive) protostars are known to be optically thick for hard photons, while optically thin for UV radiation (e.g. Stahler et al. 1980).

5.
We showed that, at least in 1D spherical calculations, the fluxlimited diffusion approximation is appropriate to studying the earliest stages of star formation, as it gives very similar results to the more complete M1 model for radiative transfer. However, our 1D spherical geometry code cannot account for multidimensional effects like the anisotropy of the radiation field.
first core has a short lifetime (a few hundred years), it is important to pursue these calculations during the second collapse, including H_{2} dissociation. Such work is in progress.
Online material
Appendix A: Effect of the numerical resolution on the energy budget through the shock. Case of a 0.01 M_{⊙} dense core
Fig. A.1 Radial profiles of various first core properties during the collapse of a 0.01 M_{⊙} clump for a core central density ρ_{c} = 1.8 × 10^{11} g cm^{3}, for calculations done with 4500 cells (dotted red line) and 18 000 cells (solid black line). From top left to bottom right: a) density; b) gas temperature; c) velocity; d) entropy; e) optical depth; f) luminosity. 
Fig. A.2 Normalized energy balance as a function of the mass for the calculations with 18 000 cells (left) and 4500 cells (right), when the central density reaches ρ_{c} = 1.8 × 10^{11} g cm^{3}. Scales are logarithmic and normalized to the rate of change of the total energy d(E_{k} + E_{p} + U_{i}) /dt. 
In contrast to the hydrodynamical case, the structure of a radiative shock can extend over large distances, depending on the optical properties of the material. For an optically thin material, the photon meanfree path is long, so the shock structure is very extended compared to the viscous scale. In this work (see Sect. 3), we present numerical calculations of dense core collapse, using a fixed resolution in mass; i.e., the mesh is not refined in the large gradient zones. Although this Lagrangean description is wellsuited to the hydrodynamical shocks, thanks to the artificial viscosity scheme, it may encounter difficulties in the case of radiationhydrodynamical flows, in particular in the optically thin region (upstream region, outside the first core).
In this appendix, we present the results of the collapse of a 0.01 M_{⊙} dense core, using the same initial ratio of thermal energy over gravitational energy as in Sect. 3 (α ~ 0.97). To investigate the effect of the numerical resolution, we performed calculations with 4500 cells and 18 000 cells, using the FLD model.
Figure A.1 shows the density, gas temperature, velocity, entropy, optical depth, and luminosity radial profiles for the two calculations at a central density ρ_{c} ~ 1.6 × 10^{11} g cm^{3}. Although there are some significative differences in the radiative precursor region (i.e. the transition region between optically thin and thick regions, where 2 < τ < 0.5) and in the estimate of the first core radius (~10%), the entropy, density, velocity, and luminosity jumps are about the same. In both calculations,
the shock is supercritical, and the amount of energy radiated away is about the same (L = 1.5 × 10^{2} with 4500 cells, and L = 1.44 × 10^{2} with 18 000 cells). This means that the overall properties of the first accretion shock, including its global energy budget remain correctly calculated even at low resolution. However, using 18 000 cells, we see that the spike in the gas temperature is resolved and that the radiative precursor length is much shorter. On the other hand, the central entropy within the first core is the same in both cases, indicating that the cooling of the first core by radiation is not affected by the lack of resolution in the radiative shock.
Figure A.2 shows the normalized energy balance at ρ_{c} = 1.8 × 10^{11} g cm^{3} for the calculations run with 18 000 cells (left) and 4500 cells (right). The figures display the rate of change of potential energy dE_{p} /dt, kinetic energy dE_{k} /dt, internal energy dU_{i} /dt, total energy d(E_{k} + E_{p} + U_{i}) /dt, and the work done by thermal pressure and radiative flux (L_{rad} + 4πr^{2}Pu). The total energy equation reads as $\frac{\mathrm{d}}{\mathrm{d}\mathit{t}}\mathrm{(}{\mathit{U}}_{\mathrm{i}}\mathrm{+}{\mathit{E}}_{\mathrm{k}}\mathrm{+}{\mathit{E}}_{\mathrm{p}}\mathrm{+}{\mathit{E}}_{\mathrm{r}}\mathrm{)}\mathrm{+}\frac{\mathit{\partial}}{\mathit{\partial r}}\mathrm{\left[}\mathrm{4}\mathit{\pi}{\mathit{r}}^{\mathrm{2}}\mathrm{\right[}\mathit{u}\mathrm{(}\mathit{p}\mathrm{+}{\mathit{P}}_{\mathrm{r}}\mathrm{)}\mathrm{+}{\mathit{F}}_{\mathrm{r}}{\mathrm{\left]}}^{\mathrm{\right]}}\mathrm{=}\mathrm{0.}$(A.1)First, we see from Fig. A.2 that the radiative pressure exerts a negligible work compared to the thermal pressure one. Comparing the energy balance of the two calculations, we see that it is globally the same, which confirms that the calculations done with 4500 cells receive the correct features of the first core accretion shock, even though the radiative structure of this shock is not resolved.
Acknowledgments
We thank the anonymous referee for comments, which helped to improve the clarity of the manuscript. Calculations were performed on the GODUNOV cluster at SAp/CEA. The research of B.C. is granted by the postdoctoral fellowships from the MaxPlanckInstitut für Astronomie. The research leading to these results received funding from the European Research Council under the European Community’s Seventh Framework Programme (FP7/20072013 Grant Agreement no. 247060).
References
 Audit, E., Charrier, P., Chièze, J. P., & Dubroca, B. 2002, unpublished [arXiv:astroph/0206281] [Google Scholar]
 Bouquet, S., Teyssier, R., & Chièze, J. P. 2000, ApJS, 127, 245 [NASA ADS] [CrossRef] [Google Scholar]
 Calvet, N., & Gullbring, E. 1998, ApJ, 509, 802 [NASA ADS] [CrossRef] [Google Scholar]
 Commerçon, B., Hennebelle, P., Audit, E., Chabrier, G., & Teyssier, R. 2008, A&A, 482, 371 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Commerçon, B., Hennebelle, P., Audit, E., Chabrier, G., & Teyssier, R. 2010, A&A, 510, L3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Commerçon, B., Teyssier, R., Audit, E., Hennebelle, P., & Chabrier, G. 2011, A&A, 529, A35 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Drake, R. P. 2005, Ap&SS, 298, 49 [NASA ADS] [CrossRef] [Google Scholar]
 Drake, R. P. 2007, Phys. Plasmas, 14, 043301 [NASA ADS] [CrossRef] [Google Scholar]
 Dubroca, B., & Feugeas, J. N. 1999, Comptes Rendus de l’Académie des Sciences, 329, 915 [Google Scholar]
 Ferguson, J. W., Alexander, D. R., Allard, F., et al. 2005, ApJ, 623, 585 [NASA ADS] [CrossRef] [Google Scholar]
 González, M., Audit, E., & Huynh, P. 2007, A&A, 464, 429 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 González, M., Audit, E., & Stehlé, C. 2009, A&A, 497, 27 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Larson, R. B. 1969, MNRAS, 145, 271 [NASA ADS] [CrossRef] [Google Scholar]
 Levermore, C. D. 1984, J. Quant. Spect. Rad. Transf., 31, 149 [NASA ADS] [CrossRef] [Google Scholar]
 Masunaga, H., & Inutsuka, S.i. 2000, ApJ, 531, 350 [NASA ADS] [CrossRef] [Google Scholar]
 Masunaga, H., Miyama, S. M., & Inutsuka, S.I. 1998, ApJ, 495, 346 [Google Scholar]
 Michaut, C., Falize, E., Cavet, C., et al. 2009, Ap&SS, 322, 77 [NASA ADS] [CrossRef] [Google Scholar]
 Mihalas, D., & Mihalas, B. W. 1984, Foundations of radiation hydrodynamics, ed. D. Mihalas, & B. W. Mihalas [Google Scholar]
 Minerbo, G. N. 1978, J. Quant. Spect. Rad. Transf., 20, 541 [Google Scholar]
 Omukai, K. 2007, PASJ, 59, 589 [NASA ADS] [Google Scholar]
 Penston, M. V. 1969, MNRAS, 144, 425 [NASA ADS] [CrossRef] [Google Scholar]
 Semenov, D., Henning, T., Helling, C., Ilgner, M., & Sedlmayr, E. 2003, A&A, 410, 611 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Shu, F. H. 1977, ApJ, 214, 488 [NASA ADS] [CrossRef] [Google Scholar]
 Stahler, S. W., Shu, F. H., & Taam, R. E. 1980, ApJ, 241, 637 [NASA ADS] [CrossRef] [Google Scholar]
 Stone, J. M., Mihalas, D., & Norman, M. L. 1992, ApJS, 80, 819 [NASA ADS] [CrossRef] [Google Scholar]
 Tscharnuter, W. M. 1987, A&A, 188, 55 [NASA ADS] [Google Scholar]
 Tscharnuter, W. M., & Winkler, K. 1979, Comp. Phys. Commun., 18, 171 [Google Scholar]
 Whitworth, A. P., & Clarke, C. J. 1997, MNRAS, 291, 578 [NASA ADS] [Google Scholar]
 Winkler, K., & Newman, M. J. 1980, ApJ, 236, 201 [NASA ADS] [CrossRef] [Google Scholar]
 Zel’Dovich, Y. B., & Raizer, Y. P. 1967, Physics of shock waves and hightemperature hydrodynamic phenomena, ed. W. D. Hayes & R. F. Probstein (New York: Academic Press) [Google Scholar]
All Tables
All Figures
Fig. 1 Temperature and density profiles in a subcritical shock (left) and a supercritical shock (right). Adapted from Zel’Dovich & Raizer (1967). 

In the text 
Fig. 2 Compression ratio r and fraction of kinetic energy radiated away at the shock, X, as a function of the upstream Mach number for a supercritical shock, in the case of an optically thin upstream material. 

In the text 
Fig. 3 Evolution of the compression ratio r and of the amount of kinetic energy radiated away, X, as a function of the upstream Mach number for a barotropic material, with density ratios r_{1} = 0.1 (black), 1 (red), and 10 (blue). The dotted lines represent the evolution for a supercritical shock with an optically thin upstream material (cf. Fig. 2). 

In the text 
Fig. 4 Evolution of the entropy jump, Δs ∝ Π /r^{γ}, as a function of the the upstream Mach number for a barotropic material, with density ratio r_{1} = 0.1 (black), 1 (red), and 10 (blue). The dotted line represents the evolution for a supercritical shock with an optically thin upstream material. 

In the text 
Fig. 5 Rosseland opacity made of a mix of Semenov et al. (2003) model at low temperature and Ferguson et al. (2005) model at high temperature. The Rosseland mean opacity is plotted as a function of temperature (xaxis) and $\mathit{R}\mathrm{=}\mathit{\rho}\mathit{/}\mathrm{\left(}{\mathit{T}}_{\mathrm{6}}^{\mathrm{3}}\mathrm{\right)}$, with T_{6} = T /10^{6}, (yaxis) using logarithm scales. 

In the text 
Fig. 6 Sketch of the estimate of ρ_{ad} used in the barotropic EOS, according to the densitytemperature distribution obtained with the M1 model. 

In the text 
Fig. 7 Evolution of the central temperature (left) and the central entropy (right) as a function of the central density during the first collapse for the M1 (solid line), FLD (dashed red line) and barotropic (dasheddotted line). 

In the text 
Fig. 8 Evolution of the central temperature, density, entropy, and optical depth with time for the M1 (solid line), FLD (dashed red line), and barotropic (dasheddotted line) models. 

In the text 
Fig. 9 Radial profiles of various properties during the collapse of a 1 M_{⊙} dense clump for a core central density ρ_{c} = 1 × 10^{10} g cm^{3}, for the M1 (solid line), FLD (dashed red line) and barotropic (dasheddotted line) models. From top left to bottom right: a) density; b) gas temperature; c) velocity; d) entropy; e) optical depth; f) luminosity; g) radiative flux; and h) integrated mass. 

In the text 
Fig. A.1 Radial profiles of various first core properties during the collapse of a 0.01 M_{⊙} clump for a core central density ρ_{c} = 1.8 × 10^{11} g cm^{3}, for calculations done with 4500 cells (dotted red line) and 18 000 cells (solid black line). From top left to bottom right: a) density; b) gas temperature; c) velocity; d) entropy; e) optical depth; f) luminosity. 

In the text 
Fig. A.2 Normalized energy balance as a function of the mass for the calculations with 18 000 cells (left) and 4500 cells (right), when the central density reaches ρ_{c} = 1.8 × 10^{11} g cm^{3}. Scales are logarithmic and normalized to the rate of change of the total energy d(E_{k} + E_{p} + U_{i}) /dt. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.