Formation of chondrules in radiative shock waves I. First results, spherical dust particles, stationary shocks
Key Words.:
meteorites, meteors, meteoroids – interplanetary medium – protoplanetary disks – shock wavesAbstract
Context:The formation of chondrules in the protoplanetary nebulae causes many questions concerning the formation process, the source of energy for melting the rims, and the composition of the origin material.
Aims:The aim of this work is to explore the heating of the chondrule in a single precursor as is typical for radiation hydrodynamical shock waves. We take into account the gasparticle friction for the duration of the shock transition and calculate the heat conduction into the chondrules. These processes are located in the protoplanetary nebulae at a region around , which is considered to be the most likely place of chondrule formation. The present models are a first step towards computing radiative shock waves occurring in a particlerich environment.
Methods:We calculated the shock waves using onedimensional, timeindependent equations of radiation hydrodynamics (RHD) involving realistic gas and dust opacities and gasparticle friction. The evolution of spherical chondrules was followed by solving the heat conduction equation on an adaptive grid.
Results:The results for the shockheating event are consistent with the cosmochemical constraints of chondrule properties. The calculations yield a relative narrow range for density or temperature to meet the requested heating rates of as extracted from cosmochemical constraints. Molecular gas, opacities with dust, and a protoplanetary nebula with accretion are necessary requirements for a fast heating process. The thermal structure in the far postshock region is not fully consistent with experimental constraints on chondrule formation since the models do not include additional molecular cooling processes.
Conclusions:
1 Introduction
Chondrules (Greek Chondros, grain) or chondren are silicate globules with sizes of to and consist mainly of olivine and pyroxene. They make up to of the meteorite volume, which is why a whole group of meteorites is called chondrites. The mineral structure of the chondrules together with their chemical composition tell a lot about their formation process. Based on their observed consistence and structure the heating has to occur very quickly at a rate of , which strongly favours melting of the origin material. A subsequent rapid cooling within a few hours is also required to produce the observed mineralogical features (see, e.g. Hood & Kring (1996) for a detail discussion of these basic properties). The analysis of various chondrules proves that multiple heating events have been present to generate their final structure. Based on these investigations the origin of the chondrules lies within regions of the early protoplanetary nebulae and their age is considered to be between 1 and 4 million years after the formation of the solar system (e.g. Amelin et al., 2002; Kita et al., 2005).
There exists a large number of models for chondrule formation. All of them contain more or less plausible assumptions. Owing to the aforementioned structural characteristics of the chondrules, the duration, temperature, and pressure at the heating event have to occur within a spatially localised zone of the protoplanetary nebulae as summarised, e.g. by Boss (1996). Shock waves appear to meet most constraints and this scenario is currently the most plausible model, as discussed, e.g. by Connolly & Love (1998), Desch et al. (2005). Since the models based on heating in shock waves make fewer assumptions on the physical environment, we will focus on these events, which have been investigated by a large amount of detailed modelling, e.g. Hood & Horanyi (1991, 1993), Ruzmaikina & Ip (1994), Hood (1998), Desch & Connolly (2002), Miura et al. (2002), and Miura & Nakamoto (2005, 2006). Ciesla & Hood (2002) and Ciesla et al. (2004) developed gasparticle suspension models with significant heating via radiation from the other chondrulesized particles. In the models of Desch & Connolly (2002) the heating process is caused either by frictional heating between gas and dust particles when traversing an adiabatic shock wave or by thermal conduction from the shockheated gas to the cooler dust particles. Desch and Connolly also took significant heating by radiation from other chondrules into account. Nevertheless, without detailed modelling of the physical processes within the protoplanetary disc the origin of these shock waves remains unanswered. In the literature several mechanisms have been suggested to explain the shock waves, i.e. accretion shocks on the surface of the protoplanetary disc (Ruzmaikina & Ip, 1994), infall of clumps of gas onto the nebulae (Tanaka et al., 1998), bow shocks produced by planetesimals on eccentric orbits (Weidenschilling et al., 1998), and Xray flares from the young Sun (Nakamoto et al., 2005). Clearly, all these scenarios give only a simplified picture of the numerous interactions between radiation field, and gas, and dust particles. In particular the treatment of adiabatic shock waves is mostly made with a simple equation of state and based on the equations of hydrodynamics without calculating a consistent radiative precursor region.
The recent work of Morris & Desch (2010) is based on a hydrodynamical shock model including a complete molecular line cooling due to , a treatment of the radiation field, and improved opacities of the solids. These authors found the effects of molecular cooling to be minimal because the combination of high column densities of water, hydrogen recombination/dissociation and radiation from the upstream regions reduce the rapid cooling times of chondrules in the postshock region.
In this work we investigate the reaction of dust particles to the structure of radiative shock waves. Based on full onedimensional timeindependent equations of radiation hydrodynamics (RHD) (e.g. Mihalas & Mihalas, 1984), we follow the evolution of a spherical dust particle that is swept through these shock transitions, which are parameterised by the Mach number. In this first study we include the effects of gasparticle drift but the chondrules are still assumed to preserve their spherical shape and are not directly coupled through their radiative properties to the equations of RHD. Their temperature structure is obtained by solving a heat conduction equation allowing a thermal expansion of the dust particles. Thanks to the full system of RHD, the overall structure of the radiative shock waves is calculated correctly together with the radiative precursor that heats the dust particles before they enter the shock jump. All models presented are located in a protoplanetary nebula at a typical distance of , the most likely origin of chondrules (Hood & Kring, 1996).
The next section deals with the basic equations describing the interaction of dust particles, gas, and radiation in onedimensional plane geometry together with the numerical method. In sect. 3 we present the results of several model calculations followed by a discussion of the possible formation process of chondrules in radiating shock waves.
2 Basic equations
2.1 Equations of radiation hydrodynamics (RHD)
The equations of radiation hydrodynamic (RHD) treat the interplay of matter and radiation by including the momentum and energy exchange between the two components. A detailed derivation of these equations can be found in Landau & Lifschitz (1981) and Mihalas & Mihalas (1984). To investigate the physical behaviour of dust particles we restricted the problem to timeindependent planar shock waves in a protoplanetary nebula in this exploratory study. We can neglect selfgravity of the gas, and the resulting system of ODEs can be solved without the need of artificial viscosity to broaden shock waves. We furthermore assumed the Eddingtonapproximation as closure condition for the radiation field as well as an ideal molecular gas with a constant adiabatic index of for the equation of state.
The following three equations are the stationary, planeparallel Euler equations for the gas interacting with the radiation field where denotes the gas density and the gas velocity. We begin with the equation of continuity
(1) 
and the equation of motion with gas pressure and radiative (Eddington) flux
(2) 
where is the socalled Rosselandmean of the opacity, is the light speed. The third equation is the gas energy equation written as
(3) 
containing the adiabatic index , the Planckmean opacity , the Eddington radiation energy density and the source function .
The next two equations are needed to describe the radiation field by writing down the Eddington moments of the radiation transport equation. Starting with the 0moment we obtain the equation of the radiation energy density
(4) 
including the Eddington radiation pressure . The 1moment leads to the equation for the radiation flux density ,
(5) 
This system of ODEs has to be solved for radiative shock waves together with additional closure conditions, and boundary conditions. Hence, we need an equation of state and we adopted in this first study an ideal gas
(6) 
The source function is given by
(7) 
for LTE, which we assumed throughout. is the temperature and the StefanBoltzmann constant. The corresponding radiation temperature can be obtained from the radiation energy density by
(8) 
Since the above equations contain three moments, , , and of the grey specific radiation intensity, we need an additional closure condition for the radiation field. Within the limit of isotropic intensity distributions we can use the wellknown Eddington approximation
(9) 
relating the second to the zeroth moment. For the radiating shock waves in a dusty environment such as the solar nebula with its high optical depths this approximation is usually valid.
2.2 Boundary conditions
According to the previous section our five basic variables are given by the gas density , the gas velocity , the gas pressure , the radiation energy density , and the radiation flux . Hence, we have to specify five far upstream values at . Denoting this gas density by and the gas velocity by , the gas pressure also fixes the gas temperature through the equation of state (6). From the equation of motion (2) we see that only leads to consistent boundary conditions, and the equation of radiation energy (4) requires and together with the source functions we have , i.e the radiation temperature is equal to the gas temperature for the far upstream region. The same conditions hold for the far downstream region for stationary problems (e.g. Zel’dovich & Raizer, 1967). As already discussed in the literature, these conditions appear to be not fulfilled in all shock computations used to study the formation of chondrules. Because the radiative flux is zero at both boundaries, a maximum (or minimum) has to occur in between and any physical possible solution has a maximum located at the shock position . Ignoring for this discussion the difference between the Rosseland and Planck mean, we introduce the optical depth by
(10) 
and can transform the system of stationary RHD equations into
(11)  
(12)  
(13)  
(14)  
(15) 
We see that the optical depth defines through a length scale and hence yields a typical time scale available for heating of particles in the precursor region.
The properties of the shock waves are parameterised by the Mach number in the upstream region; is given for an adiabatic sound velocity by
(16) 
A detailed analytical treatment of RHD shock waves can be found in Zel’dovich & Raizer (1967). We used the inverse compression ratio , and a value of
(17) 
has to be reached in the far downstream region behind the shock wave.
A close view on the aforementioned stationary RHDequations reveals that the model assumes a single shock wave with a precursor for the melting of the chondrules. However, the complexity observed in meteorites favours formation regions where many precursors may exist, e.g. due to the interaction with bow shocks of moving planetesimals and to accretion processes. A more definite answer to the evolution of chondrules therefore requires a more elaborate treatment of the various interaction processes between the gas, dust particles and radiation fields within timedependent flows. In particular, in another step it will be necessary to consider also the momentum and energy exchange between gas and particles because this study only takes into account the heating of the dust particles embedded in the gas without additional energy and momentum exchange. Hence, the cooling rates in the post shock region can disagree with data coming from furnace experiments (e.g. Hewins et al., 2005).
2.3 Gas and dust opacities
As mentioned before, we neglected the backreaction of the dust particles on the gas, e.g. the evolution of the dust particle’s size distribution function during the shock transition, which also influences the transparency of the gasdust mixture. The opacities adopted in our computations were computed for evolved planetary discs with a chemical composition like our solar system, i.e. also a solidtogas ration of about . However, the chondrules may be formed in a dustenriched environment (e.g. Ciesla & Hood, 2002) where this solidtogas ratios can be 100 or even 1000 times the solar ratio. In these zones the particle densities are much higher than the adopted values and the passage of a shock wave through this particlegas suspensions requires a more detailed treatment of the opacities (Desch et al., 2005). Increasing the dusttogas ratio increases the total opacity, and therefore all length scales defined through will become smaller, in particular the time a chondrule needs to pass the radiative precursors (cf. Fig. 2) is shortened.
Since we use frequencyintegrated moments of the radiation field, the frequencyintegrated gas and dust opacities that enter RHD equations, in particular the frequencydependent opacity , have to be integrated to obtain either the Planckmean for the optically thin case or the Rosselandmean in the optical thick case. This use of opacity tables for a dustgas mixture can only be a first approximation to a more realistic situation where the dust particles are processed by the interaction with larger bodies.
In particular, we used the Rosseland and Planckmean dust opacities for protoplanetary discs from Semenov et al. (2003), which are based on optical constants for spherical dust with porous aggregate particles and normal silicates.
2.4 Gasparticle friction
When the gas passes the shock wave, the velocity changes in the precursor region, the adiabatic shock, and the subsequent postshock cooling region. A particle changes its velocity through collisions with the moving gas and will therefore need some time to relaxate to zero relative velocity. The speed of a chondrule will be higher compared to the gas since it maintains its initial upstream velocities in the absence of friction. Consequently, the exposition of dust grains in the hot shock zone is shorter compared to the gas. To account for this effect, we solved the following system of equations for a moving chondrule at a position
(18) 
The equation of motion for a chondrule can be written as
(19) 
The simple formulae for the frictional force is related to the relative or driftvelocity by
(20) 
where is the mass of the particle and its radius. is the drag coefficient given for the numerical calculations according to (Sandin & Höfner, 2003)
(21) 
Both terms are calculated as follows
(22)  
(23) 
with
(24) 
defines the type of collision, i.e. diffusive collisions lead to and specular collisions to , denotes the temperature of the particles, is the temperature of the gas and is the thermal speed of the Maxwellian velocity distribution. Figure 1 shows the velocity difference between gas and particles and also the force on the particles plotted for a shock wave with Mach number , a preshock temperature , a preshock density and an adiabatic index . The radius of the particle is set to and the density corresponds to assuming forsterit MgSiO.
The stopping length of a typical chondrule with radius and density can be estimated by equating the mass of a dust particle with the gas mass contained in a cylinder of length and a cross section corresponding to the dust particle. This simple approximation leads roughly to
(25) 
and for the above values (also used in Fig. 1) we derived km. Hence, we expect relative velocities between gas and chondrules of about in the vicinity of the shock waves in agreement with Fig. 1.
2.5 Heat conduction in particles
As a first step we studied the heating of a spherical dust particle by heat conduction where the temperature changes because of the passage through the shock at the particle surface. In the simplest approximation we assumed a spherical shape of the particle as an initial condition in the upstream region. The particle is also assumed to be in thermal equilibrium with the surrounding gas, to be located many optical depths away from the radiating shock wave, and to move with the gas, i.e. no drift velocity .
For calculating the particle temperature we solved the onedimensional conduction equation in spherical geometry. Since the gas moves at high velocities and the diameter of the chondrule is small (at maximum a few millimeters), we can assume that the temperature is homogeneous over the spherical surface. For further applications we rewrite the heat conduction equation in a volumeintegrated form that also allows changes of the radial extension or the geometrical shape for heating (and melting). For the particle temperature the spherical heat equation reads
(26) 
where specifies the heat conduction coefficient. This coefficient can be written as
(27) 
is the heat conduction efficiency, is the chondrule density and its specific heat.
For the numerical solution of this PDE we used a standard implicit difference scheme (e.g. Richtmyer, 1967), which leads to a tridiagonal system of equations for the temperature distribution at the different particle radii. In this formulation the total radius can change in time due to a thermal expansion of the particle. The initial condition is specified by a constant temperature within the particle for , and later on we assumed an outer boundary condition of as well as at the centre. The particle position is determined by integrating Eq. (18).
2.6 Numerical calculations
For the numerical solution of the RHD system (Eqs. 17) together with the opacity tables we adopted a standard package SLGA2 (Raith, Schoenauer & Glotz, 1984) for ODEs, written in FORTRAN. We sought solutions where but . Together with the remaining upstream values the equations were integrated towards the shock position and matched by the usual RankineHugoniot conditions to the downstream values. A continuous transition of the gas temperature is not possible but the radiation energy density and the radiation temperature cross the shock wave without discontinuity.
On top of the RHD shock solutions we followed the path of a dust particle by integrating Eqs. 1824. Numerical experiments show that the type of gasparticle collision is negligible for our calculations and we adopted throughout. During the trajectory the particle temperature is not yet known, but should be in between radiation temperature and gas temperature , and for these computations we considered the gas temperature as the appropriate outer boundary condition for the chondrules, i.e. .
The calculation of the heat conduction in the dust particles was made by numerically solving the tridiagonal equation system, which we obtained from discretisation of the heat conduction equation in spherical symmetry (u.v. 2.5). We also included the thermal expansion and estimated the amount of radial increase by taking a linear expansion coefficient of together with a temperature change of . We obtained an increase of or about 1.5% for a mmsized chondrule. Clearly, this effect is negligible but we adopted a finite volume discrete version of Eq. (26) for subsequent applications, studying also the temporal evolution from fluffy to spherical particles.
3 Results
For all RHD shock wave models listed in Table 1 we considered a spherical dust particle with a radius of . According to Hood & Kring (1996), this size is close to the largest chondrules and we assumed a molecular gas with an adiabatic index of throughout. According to Hewins et al. (2005), the peak temperatures during the forming process of the chondrule are in a range between . These gas temperatures destroy all molecules and thereby also change the adiabatic index to a value of . However, to study the effects within the radiative precursor as well as the radiative cooling effects behind the shock waves, we tried to keep the equation of state as simple as possible for the moment. The chondrules were assumed to consist of forsterit (MgSiO) with the following particle characteristics. The particle density was set to , the heat conductivity and the specific heat were taken from the literature (e.g. Krishnaiah, Singh & Jadha, 2004; Yingwei, 2004). The linear heat expansion coefficient is .
The protoplanetary disc with accretion has been modelled by Bell et al. (1997), where an angular momentum transport coefficient leads to a mass accretion rate of . Outside of these discs are gravitationally unstable. Taking a typical solar distance of we found a gas density of in the mid plane at a temperature . According to Hood & Kring (1996), the temperature in the chondrule precursor region must not exceed so that in chondrules the existence of chemical compounds such as FeS can be explained. Furthermore, the melting time has to be in the range of about , which requires heating rates of and the subsequent temperature decline has to occur over a period of an hour or more (Hewins et al., 2005).
In Figure 2 a schematic curve illustrates the basic temperature structure as required from the cosmochemical evidence found within chondrules. Table 2 lists several models with a maximum temperature between and (see figure 3 and 4). The initial density and the temperature in the upstream region () correspond to different Mach numbers of the shock transitions. The adopted numbers vary around values suggested by simple models of the protoplanetary nebulae.
Since the spatial structures within a radiative shock can be approximated by exponential functions, see, e.g. Zel’dovich & Raizer (1967), we assumed for our analysis of the precursor region () and the relaxation zone ()
(28) 
which also defines the flow time through these regions by
(29) 
denotes the distance from the shock wave located at . The quantiles are also given in Table 1 allowing a more detailed characterisation of the shock structures with their typical length and time scales shaping the dust particles.
The temperature difference between the surface and the center of a chondrule is relative small, as seen from the values of . Figure 5 shows the maximum temperature in a spherical body in dependence of its radius with by transition of RHD shock fronts.
Model  1  2  3  4  5 

2.0  1.0  3.0  2.0  2.0  
300  500  100  
3.97  3.97  3.98  2.92  7.23  
Model  6  7  8  9  10 
5.83  5.7  6.0  5.83  
300  400  200  
2100  
4.48  4.48  4.48  3.80  5.61  
26.8  26.2  28.8  20.7  33.7 
4 Discussion
For model 1 with initial temperature and peak temperature an initial density was chosen so that the requirement of a heating rate is fulfilled. Model 2 shows that at initial density of the necessary heating rate is not achieved. A higher initial density , like in model 3, leads to significantly higher heating rates. Models 4 and 5 show variations of the initial temperature , based on model 1. A temperature decrease leads to higher heating rates (model 5). The necessary heating rates cannot be reached any more by model 4. Models 1 to 5 are almost subcritical shock waves, models 6 to 10 with peak temperatures constitute supercritical shock waves. With density and temperature variations models 6 to 10 react in the same way as models 1 to 5, but the changes have a greater effect on the duration of the preheating wave and thereby on the heating rate. The increase of temperature in supercritical shock waves at the beginning of the preheating waves fulfills the necessary heating rates. The particles stay for a long time in the preheating wave until the shock occurs.
Investigating the range of possible initial conditions as constrained by models of protoplanetary discs, we found that the gas can be heated up to and that the Rosseland and Planckmean opacities with dust can meet the basic requirements of shockbased models. These RHD shock waves provide the requested heating rates. Calculations using a monoatomic gas with , opacities without dust (Alexander opacities of Grevesse & Noels (1993)) and a disc models without accretion (e.g. the minimum mass nebulae of Black & Metthews (1985)) reveal that the particles are propagating within the preheating wave for days to years, which does not match the required heating rates. Molecular gas and the dustenriched opacities result in better agreements with the observational facts of chondrule research. The models of a protoplanetary disc including accretion (Bell et al., 1997) have much higher gas densities in the mid plane around . This gas densities are essential to compute the requested heating rates. Summarising we emphasize that the increase of the upstream density together with a decrease of the temperature makes molecules and therefore as well as higher opacities possible. All these effects reduce the duration of the dust particle in the preheating wave and increase the heating rates.
The calculated temperature differences between surface and core of the chondrules (Fig. 5) are correct as long as the chondrules remain solid particles. But up to now we have not incorporated phase transitions in this chondrule model. Clearly, strong shock waves with high Mach numbers as well as shock waves with smaller preheating waves will produce larger temperature differences but are unlikely to ensure the stringent cosmochemical constraints. As depicted in model 10, the preheating wave vanishes completely and the sudden increase of the temperature at the outer boundary yields the largest temperature differences within chondrules. For these conditions it is easier to explain partial or shell melting of the particles.
Immediately after the shock passage we observe a rapid cooling on a time scale of minutes owing to radiative losses determined by the opacities. This effect is followed by a longer time scale of several days where nonequilibrium processes further reduce the post shock temperature of the gas and the chondrules. The rapid cooling rates of our models of about K/hr seem to disagree with furnace experiments, which indicate values of around K/hr. The reason for this discrepancy is due to different time scales involved during cooling processes in the post shock region. The calculations of, e.g. Desch et al. (2005) include the energy and momentum exchange between the chondrules and the gas. Therefore the temperature structure in the post shock region is controlled over longer time scales by an additional cooling term , which includes the radiative cooling by molecules. Consequently, the far downstream gas temperature can be set to the far upstream temperature . This leads to overall cooling rates of K/hr, which cannot be obtained within our physical description, which neglects this additional cooling term. The boundary condition demands for the gas temperature as specified, e.g. in Desch et al. (2005), and the post shock region cools down to the initial temperature by radiative looses into the surrounding protosolar nebula. However, this additional cooling occurring on a longer spatial and temporal scale reduces the overall cooling rates to less than K/hr. From furnace experiments (Hewins et al., 2005) it has become clear that the exact shape of cooling path has no influence on the final cooling rate.
The recent work of Morris & Desch (2010) reveals how additional physical mechanisms like molecular line cooling due to water influence the cooling history in the postshock region. However, these authors emphasize that the postshock region is only slightly changed because several effects partly cancel each other out, i.e. radiative heating from the shock front as well as recombination and dissociation of molecules. Summarising these recent results of Morris & Desch (2010) we note that the improved nebular shock models can meet all cosmochemical constraints of condrule formation. The shock speeds necessary to melt chondrules have to be increased, e.g. from to .
Several models have been developed during the past years to calculate the behaviour of dust particles within shock waves (e.g. Desch & Connolly, 2002; Miura & Nakamoto, 2005). However, none of the proposed models have included the full set of RHD equations with the most recent gas and dust opacities, which are essential for the correct structure of the precursor region. This zone is shaped by the radiation transmitted from the shocked downstream gas and the dust particles are advected through this region towards the shock wave. As inferred from the cosmochemical evidence, this region provides a basic heating process and will therefore restrict the possible Mach numbers. As seen from the models (cf. Table 1), already moderate Mach numbers around (for ) meet the constraints of chondrule formation with parameters reasonably within protoplanetary nebula models.
Opacities changes due to dustenriched environments caused by fragmentation and collisional debris (e.g. Ciesla & Hood, 2002) may significantly increase solidtogas ratios, and values of or even times the solar ratio can be expected. Clearly, shock waves propagating in such a densely dustpopulated environment requires a more detailed treatment of the opacities (Desch et al., 2005), as adopted in these simulations. In particular, all time and length scales controlled by the opacity are then changed.
As mentioned already, the shape and/or the mass of the particles can be modified as they propagate through the shock region. First, the partial or total melting could lead to a mass loss by evaporation or, depending of the surface tension, also to a rearrangement of the surface, i.e. the particles can become more spherical. Secondly, drift velocities acting on the surface can result in an erosion and thereby decrease the particles mass. Thirdly, the radiative interactions with the dust particles together with the additional heat needed to melt the particle should be included into these models of radiative shockheated chondrules. The dissociation of molecules also demands a more realistic treatment of the equation of state. These structural effects of the dust particles will be discussed and included in a forthcoming paper in more detail.
Acknowledgements.
The publication is supported by the Austrian Science Fund (FWF).References
 Amelin Y., Krot A. N., Hutcheon I. D., Ulyanov A. A., 2002, Science 297, 16781683
 Bell K. R., Cassen P. M., Klahr H. H., Henning Th., 1997, ApJ486, 372387
 Black D. C., Metthews M. S. (Eds.), 1985, Protostars and Planets II, University of Arizona Press, Arizona
 Boss A. P., 1996, in: Chondrules and the Protoplanetary Disk, Cambridge University Press, 257263
 Ciesla F. J., Hood L. L., 2002, Icarus 158, 281293
 Ciesla F. J., Hood L. L., Weidenschilling S. J., 2004, Meteoritics & Planetary Science 39, Nr. 11, 18091821
 Connolly H. C., Love S. G., 1998, Science Vol. 280, No. 5360, p. 6267
 Desch S. J., Connolly Jr. H. C., 2002, Meteorit. Planet. Sci. 37, 183
 Desch S. J., Ciesla F.J, Hood L.L., Nakamoto, T. 2005, in: Chondrites and the Protoplanetary Disk, Eds. A.N. Krot, E.R.D. Scott, B. Reipurth, ASP, San Francisco, Vol. 341, p. 849

Grevesse N., Noels A., 1993, Opal Opacity Tables, active: 2007,
Hypertext: http://wwwphys.llnl.gov/Research/OPAL/existing.html  Hewins R. H., Connolly Jr. H. C., Lofgren G. E., Libourel G., 2005, in: Chondrites and the Protoplanetary Disk, ASP Conference Series, Vol. 341, p. 286
 Hood L. L., 1998, Meteorit. Planet. Sci. 33, 97
 Hood L. L., Horanyi M., 1991, Icarus 93, 259
 Hood L. L., Horanyi M., 1993, Icarus 106, 179
 Hood L. L., Kring D. A., 1996, in: Chondrules and the Protoplanetary Disk, Cambridge University Press, 265276
 Kite et al, 2005, in: Chondrites and the Protoplanetary Disk, Eds. A.N. Krot, E.R.D. Scott, B. Reipurth, ASP, San Francisco, Vol. 141, p. 558
 Krishnaiah S., Singh D.N., Jadha G.N., 2004, International Journal of Rock Mechanics and Mining Sciences, Vol. 41, 877882
 Landau L.D., Lifschitz E.M., 1981, Lehrbuch der Theoretischen Physik VI, Hydrodynamik, Akademischer Verlag, Berlin
 Mihalas D., Mihalas B.W., 1984, Foundations of Radiation Hydrodynamics, Oxford University Press
 Miura H., Nakamoto T., Susa H., 2002, Icarus 160, 258
 Miura H., Nakamoto T., 2005, Icarus 175, 289
 Miura H., Nakamoto T., 2006, ApJ651, 1272
 Morris M.A., Desch S.J., 2010 ApJ722, 1474
 Nakamoto T., Hayashi M. R., Kita N. T., Tachibana S., 2005, ASP Conference Series, Vol. 341, p. 883
 Raith K., Schoenauer W., Glotz G., 1984, SLGA2Ein selbststeuerndes Lösungsverfahren für Anfangswertprobleme bei gewöhnlichen Differentialgleichungen, Interner Bericht 13/79 des Rechenzentrums der Universität Karlsruhe
 Richtmyer, R.D., Morton, K.W., 1967, Difference Methods for Initial Value Problems, Interscience Publishers New York
 Ruzmaikina T. V., Ip W. H., 1994, Icarus 112, 430
 Sandin C., Höfner S., 2003, A&A 398, 253
 Semenov D., Henning Th., Helling Ch., Ilgner M., Sedlmayr E., 2003, A&A 410, 611
 Tanaka K. K., Tanaka H., Nakazawa K., Nakagawa Y., 1998, Icarus 134, 137
 Weidenschilling S. J., Marzari F., Hood L. L., 1998, Science 279, 681

Yingwei Fei, 2004, Thermal Expansion, active: 2007,
Hypertext: http://www.agu.org/reference/minphys/6 fei.pdf  Zel’dovich Ya. B., Raizer Yu. P., 1967, Physics of Shock Waves and HighTemperature Hydrodynamic Phenomena, Academic Press, New York