Formation of chondrules in radiative shock waves
Key Words.:meteorites, meteors, meteoroids – interplanetary medium – protoplanetary disks – shock waves
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 gas-particle 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 particle-rich environment.
Methods:We calculated the shock waves using one-dimensional, time-independent equations of radiation hydrodynamics (RHD) involving realistic gas- and dust opacities and gas-particle friction. The evolution of spherical chondrules was followed by solving the heat conduction equation on an adaptive grid.
Results:The results for the shock-heating 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 post-shock region is not fully consistent with experimental constraints on chondrule formation since the models do not include additional molecular cooling processes.
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 gas-particle suspension models with significant heating via radiation from the other chondrule-sized 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 shock-heated 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 X-ray 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 post-shock region.
In this work we investigate the reaction of dust particles to the structure of radiative shock waves. Based on full one-dimensional time-independent 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 gas-particle 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 one-dimensional 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 time-independent planar shock waves in a protoplanetary nebula in this exploratory study. We can neglect self-gravity 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 Eddington-approximation 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, plane-parallel 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
and the equation of motion with gas pressure and radiative (Eddington) flux
where is the so-called Rosseland-mean of the opacity, is the light speed. The third equation is the gas energy equation written as
containing the adiabatic index , the Planck-mean 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 0-moment we obtain the equation of the radiation energy density
including the Eddington radiation pressure . The 1-moment leads to the equation for the radiation flux density ,
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
The source function is given by
for LTE, which we assumed throughout. is the temperature and the Stefan-Boltzmann constant. The corresponding radiation temperature can be obtained from the radiation energy density by
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 well-known Eddington approximation
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 up-stream 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 up-stream region. The same conditions hold for the far down-stream 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
and can transform the system of stationary RHD equations into
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 up-stream region; is given for an adiabatic sound velocity by
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
has to be reached in the far down-stream region behind the shock wave.
A close view on the aforementioned stationary RHD-equations 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 time-dependent 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 gas-dust 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 solid-to-gas ration of about . However, the chondrules may be formed in a dust-enriched environment (e.g. Ciesla & Hood, 2002) where this solid-to-gas 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 particle-gas suspensions requires a more detailed treatment of the opacities (Desch et al., 2005). Increasing the dust-to-gas 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 frequency-integrated moments of the radiation field, the frequency-integrated gas and dust opacities that enter RHD equations, in particular the frequency-dependent opacity , have to be integrated to obtain either the Planck-mean for the optically thin case or the Rosseland-mean in the optical thick case. This use of opacity tables for a dust-gas 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 Planck-mean 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 Gas-particle 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 up-stream 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
The equation of motion for a chondrule can be written as
The simple formulae for the frictional force is related to the relative or drift-velocity by
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)
Both terms are calculated as follows
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
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 one-dimensional 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 volume-integrated 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
where specifies the heat conduction coefficient. This coefficient can be written as
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. 1-7) 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 up-stream values the equations were integrated towards the shock position and matched by the usual Rankine-Hugoniot conditions to the down-stream 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. 18-24. Numerical experiments show that the type of gas-particle 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 mm-sized 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.
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 up-stream 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 ()
which also defines the flow time through these regions by
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.
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 Planck-mean opacities with dust can meet the basic requirements of shock-based 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 dust-enriched 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 up-stream 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 non-equilibrium 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 down-stream gas temperature can be set to the far up-stream 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 post-shock region. However, these authors emphasize that the post-shock 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 down-stream 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 dust-enriched environments caused by fragmentation and collisional debris (e.g. Ciesla & Hood, 2002) may significantly increase solid-to-gas ratios, and values of or even times the solar ratio can be expected. Clearly, shock waves propagating in such a densely dust-populated 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 shock-heated 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).
- Amelin Y., Krot A. N., Hutcheon I. D., Ulyanov A. A., 2002, Science 297, 1678-1683
- Bell K. R., Cassen P. M., Klahr H. H., Henning Th., 1997, ApJ486, 372-387
- 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, 257-263
- Ciesla F. J., Hood L. L., 2002, Icarus 158, 281-293
- Ciesla F. J., Hood L. L., Weidenschilling S. J., 2004, Meteoritics & Planetary Science 39, Nr. 11, 1809-1821
- Connolly H. C., Love S. G., 1998, Science Vol. 280, No. 5360, p. 62-67
- 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,
- 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, 265-276
- 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, 877-882
- 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, SLGA2-Ein 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 High-Temperature Hydrodynamic Phenomena, Academic Press, New York