© 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 so-called 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 radiation-hydrodynamics (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).

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 free-fall 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 flux-limited 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 semi-analytical 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 first-Larson 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 semi-analytical model that allows a simple interpretation of the numerical results. Eventually, Sect. 5 summarizes our main results.

2. Radiative shock – A semi-analytic 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₂ rises, producing a radiative flux of , 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₂. A shock with T₋ < T₂ is called a subcritical shock. Because the material entering the shock is preheated, the postshock temperature T ₊  exceeds its asymptotic equilibrium value T₂, 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₂. The shock velocity at which the postshock and preshock temperatures are equal defines the critical shock. For higher shock velocities, T₋ cannot exceed T₂, the excess energy forces the radiative precursor further into the upstream region with a temperature close to T₂. Pre- and post-shock 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₂. 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 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 (4)Since matter and radiation are in equilibrium and the material is opaque, we have E_r = 3P_r = a_RT⁴. Defining the compression ratio r = ρ₂ /ρ₁ = u₁ /u₂, we can rewrite Eqs. (2) and (3) in the non-dimensional form (5)and (6)where Π = P₂ /P₁, and ℳ₁ is the hydrodynamic Mach number, ℳ₁ = u₁ /c_s1 = u₁ /(γP₁ /ρ₁)^1 /2. The coupled Eqs. (5) and (6) are solved numerically to get the variations of Π and r as function of ℳ₁ and α₁.

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⁴, 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 , whereas it rises as in a non-radiating shock. If P_r /P < 10^-5, the non-radiating fluid approximation becomes valid for an upstream Mach number ℳ₁ < 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 low-mass 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 where ΔF_r = F_r2 − F_r−. We consider the case of a non-zero net flux |F_r2 − F_r−| across the shock. The jump relations, derived from the conservation Eqs. (1)–(3), read (10)and (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₂ and then h₋ = h₂. Using the same adimensional parameters as for the opaque case, we have (12)and (13)Since the shock is isothermal, r = Π and (14)Eventually, the radiative flux discontinuity, normalized to the upstream kinetic energy, is given by (15)where X thus represents the amount of incident kinetic energy radiated away at the shock front.

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 non-zero 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: (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 sub-critical 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 (17)with . If we evaluate the flux just at the discontinuity, we find the maximum preheating temperature T₋ from (18)We can now determine the critical temperature T_cr at which T₋ equals T₂(19)so that (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 (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₁ = 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₂ since all the upstream kinetic energy is radiated away: (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) (23)with ℛ = k_B /μm_H. The preshock temperature T₋ is estimated as (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 with the barotropic equation of state (EOS) as a closure relation (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 (28)Using the barotropic EOS yields (29)where r₁ = ρ₁ /ρ_ad and r₂ = ρ₂ /ρ_ad. Eventually, we get r by solving (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 (31)where r and Π are given by (30).

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₁ = 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₁ increases (note that for r₁ = 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 (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₁ = 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 new-born protostars. The entropy is indeed a key quantity for pre-main sequence evolutionary models as it entirely determines the mass-radius relationship of an adiabatic object, like the first and second Larson cores.

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 frequency-dependent 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²ρu. In the theory, the accretion rate is generally defined as (e.g. Stahler et al. 1980) (33)where ξ is a dimensionless coefficient and c_s0 the isothermal sound speed. The value of ξ is estimated assuming that a non-magnetic, non-rotating cloud, whose initial state corresponds approximately to a balance between thermal support and self-gravity, has comparable free-fall time, , and sound-crossing 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 expansion-wave solution, whereas the dynamic Larson-Penston solution yields ξ ~ 50 (Larson 1969; Penston 1969).

Another quantity directly comparable to the theory is the accretion luminosity, usually expressed as (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 flux-limited diffusion approximation (FLD), and the M1 model. The RHD equations are integrated in their non-conservative 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 Raphson-Newton 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₂, 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: (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 flux-limited diffusion (FLD) approximation and the M1 model.

3.4.3. The flux-limited diffusion approximation

The diffusion approximation is the most widely used moment model of radiation transport. The diffusion limit is valid when the photon mean-free 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 (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 (37)where λ = λ(R) is the flux limiter.

In this study, we retain the flux limiter derived by Minerbo (1978)(38)with R = |∇E_r| /(κ_RE_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 diffusion-type equation on the radiative energy is thus obtained (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 ad-hoc 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 (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) (41)where χ is the Eddington factor, I the identity matrix, and 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 (42)It is then trivial to recover the two asymptotic regimes of radiative transfer. When f → 0, χ = 1 /3, and P_r = 1 /3E_rI, 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 free-streaming 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 (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 (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 (x-axis) and , with T6 = T /106, (y-axis) using logarithm scales.

• In the M1 model, we consider the equations governing the evolution of an inviscid, radiating fluid (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². 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₀ = 1 M_⊙, temperature T₀ = 10 K (c_s0 = 0.187 km s^-1), and radius R₀ = 10⁴ AU. This initial setup corresponds to a ratio α of thermal to gravitational energies of α = 5R₀k_BT₀ /(2GM₀μm_H) ~ 0.97 and to a free-fall time t_ff ~ 1.77 × 10⁵ 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.

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).

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 Larson-Penston solution than to Shu’s expansion wave solution.

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 (46)For the barotropic EOS, we have (47)At high density, . The limiting value for the entropy is then  ~2.03 × 10⁹ 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.

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 low-mass protostar determines its mass-radius relationship, hence its pre-main sequence subsequent evolution. According to the first principle of thermodynamics, the global entropy loss of the first core can be estimated as (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 (49)At ρ_c = 1 × 10^-10 g cm^-3, we get M_fc ~ 2.3 × 10^-2 M_⊙. The entropy loss is thus (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. ρ₁ ~ 7 × 10^-14 g cm^-3, u₁ ~ 2 × 10⁵ cm s^-1, and T₁ ~ 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 semi-analytical 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 super-critical 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 Larson-Penston and Shu self-similar 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 free-fall velocity (51)The preshock density is given assuming a profile  ∝ r^-2 in the free-falling envelope (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 ℳ₁ ~ 2 and Fig. 2 shows that X ~ 1. The shock temperature thus reads (cf. Eq. (22)) (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 ρ₁ = 5.8 × 10^-14 g cm^-3, u₁ ~ 2.3 × 10⁵ cm s^-1 and T₁ ~ 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 = ρ₂ /ρ₁ ~ 50, which gives ρ₂ ~ 2.9 × 10^-12 g cm^-3, and u₂ = 4.6 × 10³ 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: ρ₁ ~ 7 × 10^-14 g cm^-3, u₁ ~ 2 × 10⁵ cm s^-1 and T₁ = 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.15 × 10^-12 g cm^-3, u₂ ~ 10⁴ cm s^-1 ≪ u₁, 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 frequency-dependent 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 flux-limited 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.

This study confirms the need to solve the complete RHD equations, i.e. to correctly take the coupling between radiation and hydrodynamics into account when addressing strong shock conditions, as occurs during the formation of the first Larson core in the context of star formation. Such a complete, consistent treatment of radiation and hydrodynamics is necessary to correctly handle the cooling properties of the accreting gas and thus to obtain the correct mechanical and thermal properties of the first core. This, in turn, sets up the initial conditions for the second collapse and the formation of the second core. Indeed, since the

first core has a short lifetime (a few hundred years), it is important to pursue these calculations during the second collapse, including H₂ dissociation. Such work is in progress.

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 Max-Planck-Institut 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/2007-2013 Grant Agreement no. 247060).

References

Online material

Appendix A: Effect of the numerical resolution on the energy budget through the shock. Case of a 0.01 M_⊙ dense core

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 mean-free 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 well-suited to the hydrodynamical shocks, thanks to the artificial viscosity scheme, it may encounter difficulties in the case of radiation-hydrodynamical 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²Pu). The total energy equation reads as (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.

All Tables

All Figures