A Drag force/frictional heating

A critical analysis of shock models for chondrule formation

Abstract

In recent years many models of chondrule formation have been proposed. One of those models is the processing of dust in shock waves in protoplanetary disks. In this model, the dust and the chondrule precursors are overrun by shock waves, which heat them up by frictional heating and thermal exchange with the gas.

In this paper we reanalyze the nebular shock model of chondrule formation and focus on the downstream boundary condition. We show that for large-scale plane-parallel chondrule-melting shocks the postshock equilibrium temperature is too high to avoid volatile loss. Even if we include radiative cooling in lateral directions out of the disk plane into our model (thereby breaking strict plane-parallel geometry) we find that for a realistic vertical extent of the solar nebula disk the temperature decline is not fast enough. On the other hand, if we assume that the shock is entirely optically thin so that particles can radiate freely, the cooling rates are too high to produce the observed chondrules textures. Global nebular shocks are therefore problematic as the primary sources of chondrules.

keywords:
Disks, Meteorites, Radiative transfer, Solar Nebula, Thermal histories
1\cortext

[cor1]Corresponding author

1 Introduction

The origin of chondrules is one of the biggest mysteries in meteoritics. These 0.11 mm sized silicate once molten droplets, abundantly found in chondritic meteorites, must have cooled and solidified within a matter hours (e.g. Hewins et al., 2005). From short-lived radionuclide chronology data (e.g. Villeneuve et al., 2009) it is known that this must have taken place during the first few million years after the start of the solar system, during the phase when the sun was still likely surrounded by a gaseous disk (the “solar nebula”). What makes chondrule formation mysterious is that this few-hour cooling time is orders of magnitude shorter than the typical few-million-year time scale of evolution of the solar nebula. Chondrules can thus not be a natural product of the gradual cooling-down of the nebula. Instead, chondrules must have formed during “flash heating events” of some kind – but which kind is not yet known. There exists a multitude of theories as to what these flash heating events could have been. Boss (1996) and Ciesla (2005) give nice overviews of these theories and their pros and cons. So far none of these theories has been universally accepted.

One of the most popular theories is that nebular shock waves can melt small dust aggregates in the nebula, causing them to become melt droplets and allowing them to cool again and solidify (Hood and Horanyi, 1991). The origin of such shocks could, for instance, be gravitational instabilities in the disk (e.g. Boley and Durisen, 2008) or the effect of a gas giant planet (e.g. Kley and Nelson, 2012). Detailed 1-D models of the structure of such radiative shocks, and the formation of chondrules in them, were presented by Iida et al. (2001), Desch and Connolly (2002), Ciesla and Hood (2002) or Morris and Desch (2010). These models show that such a shock, in an optically thick solar nebula, would lead to a temperature spike in the gas and the dust that lasts for only a few seconds to minutes with cooling rates of up to K/hour, followed by a more gradual cooling lasting several hours, with cooling rates of the order of K/hour. These appear to be the right conditions for chondrule formation, which is one of the reasons why this model is one of the favored models of chondrule formation nowadays.

In this paper we revisit this shock-induced chondrule formation model. Our aim is to investigate the role of the downstream boundary condition and the role of sideways radiative cooling.

If the shock is a large scale shock, e.g. due to a global gravitational instability, then on a small scale (the scale of several radiative mean free paths ) the shock can be modeled as a infinite 1-D radiation hydrodynamic flow problem. The “infinite 1-D” in this context means that the 3-D geometry of the problem only becomes important on scales , so that on scales a 1-D geometry can be safely assumed where the pre-shock boundary is set at and the post-shock boundary is set at . With the pressure scale height of a typical minimum mass solar nebula is several hundred times larger than the optical mean free path. The shock is assumed to be stationary in the -coordinate system, i.e. the coordinates move along with the shock and the shock is always at . In such an infinite 1-D system the full shock structure can be reconstructed when the physical variables at the upsteam boundary are all given. No downsteam boundaries at need to be given. In fact, the physical variables at follow uniquely by demanding that the mass-, momentum- and energy-flux at equals those of , but with subsonic gas velocity at . This gives a global Rankine-Hugoniot condition including all the radiative and dust physics. Note that right at the shock at the jump in the gas variables is given by a local Rankine-Hugoniot condition in which only the gas fluxes on both sides are set equal. We will work out these shock models in section 2 and show that after the temperature spike they lack a slower few-hour cooling phase. Instead, they stay at a constant temperature, i.e. the temperature in accordance with the global Rankine-Hugoniot condition. We will investigate in section 3.1 whether the chondrule peak temperature can be made high enough for melting while keeping the post-spike temperature low enough to retain volatile elements. According to Fedkin et al. (2012) and Yu et al. (1996) the high-temperature phase should be of the order tens of minutes rather than of hours. From textural constraints, Hewins et al. (2005) and Desch et al. (2012) conclude that the cooling rates have to be of the order of 10 – 10 K/hour.

Since the infinite 1-D shock solution is a geometric approximation we will implement effects of sideways cooling (i.e., from top and bottom of the disk) in section 3.2 to improve the realism of the model. This will re-introduce the slow cooling phase after the temperature spike, but we find that this slow cooling phase is of the order of weeks/months/years rather than hours, because 2-D/3-D radiative diffusion will take place on -scales of the same order as the vertical scale height of the disk.

Finally we discuss in section 3.3 two versions of the shock-induced chondrule formation where the cooling time can be rapid. One is a locally induced one, for instance due to a supersonic planetesimal bow shock (Hood, 1998; Morris et al., 2012). The other is if the disk is fully optically thin, allowing the post-shock material to cool straight to infinity.

2 The Model

Our one-dimensional shock model is built on the work by Desch and Connolly (2002). We used their approach, generalized it to arbitrarely many gas species, particle populations and chemical reactions and modified it where we think modifications or corrections are needed.

In this model we assume that all the parameters (densities, temperatures, velocities, etc.) only change along one direction, the -direction. Therefore, the model consists of infinitely extended, plane-parallel layers of constant temperatures, velocities and densities.

2.1 Radiative Transfer

Even though the model is one-dimensional, we must allow the photons to move into all three directions. But fortunately, because of the plane-parallel assumption the radiative transfer equations here can be vastly simplified.

In this case the optical depth is a monotonic increasing function of the distance from the post-shock boundary and therefore a meassure of the position inside the computational domain:

(1)

where is the location of the post-shock boundary. This means the optical depth increases from at the post-shock boundary to a maximum value of at the pre-shock boundary. is the absorption coefficient, which is the sum of the contribution of the gas and the particles:

(2)

The absoption coefficient of the gas is the product of the gas mass density and the temperature-dependent Planck-mean opacity , which we took from Semenov et al. (2003) (see figure 1). The absorption coefficient of the particles is the product of their number densities , their geometrical cross-section and their absoprtion efficiency , which we adopted from Desch and Connolly (2002) as

(3)

where is the particle radius.

Figure 1: Planck-mean opacity from Semenov et al. (2003). At K (vertical line) the fine-grained dust associated with the gas gets evaporated. Therefore, the opacity drops by one order of magnitude at this temperature. This is responsible for the “opacity knee” in the pre-heating phase, where the gas becomes optically thin.

To calculate the thermal histories of the particles one has to calculate the mean intensity at every position Our radiative transfer is grey, i.e. wavelength independent. The mean intensity is defined as the average intensity per solid angle coming from all directions

(4)

In the plane-parallel assumption the intensity depends only on the position and the angle between the incoming ray and the -axis. is defined as . Therefore, the mean intensity can be simplified to

(5)

by using the exponential integrals

(6)

(see Mihalas and Weibel-Mihalas, 1999). and are the incoming radiations from both boundaries. The source function is defined as

(7)

which are the wave-length integrated Plank functions of the gas and the particles at their given temperatures weighted by their respective emission efficiencies, which are the same as their absoprtion efficiencies according Kirchoff’s law. is the Stefan-Boltzmann constant.

The radiative flux is defined as the net flux of radiative energy in -direction

(8)

Later, we need the derivative of the radiative flux with respect to the -direction. In the plane-parallel assumption this is simply

(9)

where we now included different particle populations. Therefore, we have to take the sum over all populations.

2.2 Hydrodynamics

To calculate the evolution of the gas and particle parameters we used the one-dimensional stationary Euler equations, which are

(10)
(11)
(12)

with the total mass density , the velocity , the pressure and kinetic and thermal energy density . These equations are simply conservation laws of mass, momentum and energy respectively. S, S and S are the source terms of the designated quantities.

Particle dynamics

In this model we have different particle populations with radius , temperature , velocity and number density . The particles are not allowed to evaporate completely. This means the continuity equation (10) states

(13)

The particles are accelerated or decelerated by the drag force acting on population . The drag force is given by

(14)

(Gombosi et al., 1986), with the gas’ mass density and the drag coefficient of population (see A for futher details). If the gas velocity is higher than the particle velocity, the particles get accelerated and vice versa. With that, the force equation of the particles is

(15)

with the particles’ mass and the comoving derivative , which is defined as

(16)

The last equality is for the one-dimensional, stationary case. Using this, the force equation (15) yields

(17)

The energy budget of the particles is given by the balance between frictional heating by the gas drag, thermal contact with the gas and radiative heating on one hand and radiative cooling on the other hand. The effects of frictional heating and thermal contact with the gas are combined into a single heating rate per unit surface area of the particle given by

(18)

(Gombosi et al., 1986), with the heat transfer coefficient and the recovery temperature . In the limit of there can be still heat exchange if (see A for further details). The radiative heating rate per unit surface area of the paticles is given by , which is the balance between radiation received by the mean intensity and the energy radiated away according the Stefan-Boltzmann law. Combining the effects of frictional and thermal heating and radiative heating the net heating rate per unit surface area is .

This net heating rate can raise the particles’ temperatures and/or evaporate them. We assume that the fraction of the net heating rate that goes into evaporation is , which is a function of the particle temperature. increases monotonically from 0 to 1 within a temperature interval centered on the evaporation temperature of 2000 K used by Desch and Connolly (2002). We have arbitrarely chosen (see figure 2). We had to do this because of numerical reasons (strong if-statements can make it impossible to estimate the Jacobian for the implicit integration scheme), but the positive side effect is that we can account now for the heterogeneity of the particle material with different evaporation temperatures.

Figure 2: The fraction of evaporation yields for a smearing of the evaporation temperature of 2000 K within our choice of the interval .

Therefore, the change in the particles’ temperature is given by

(19)

with the specific heat capacity at constant pressure of the particle material. Using again the definition of the comoving derivative (16), the particles’ material mass density and using the change in the particles’ temperature can be written as

(20)

The other part of the net heating rate, that goes into evaporation, is given by

(21)

with the latent heat of evaporation of the particle material. The negative sign is needed because the heat received here goes into shrinking the particles’ mass. This yields for the change in the particle radius

(22)

If the net heating rate is negative then is set to zero to avoid artificial condensation in equation (22). We neglect condensation and nucleation here because we assume it takes place on much longer timescales than considered here. We want to point out that this implementation of evaporation is not the correct physical treatment. To do it in an correct way one has to integrate the Hertz-Knudsen equation and possibly a nucleation model, which makes the whole system even more complex. We also want to point out that in our conclusive simulations, the particle temperature is always safely below the evaporation temperature.

With the continuity equation of the particles (13) and the change in the particle velocity (17), the change in the particles’ number densities can be calculated

(23)

With equations (17), (20), (22) and (23) all the particle parameters can be calculated throughout the computational domain knowing the gas parameters.

Gas Dynamics

The gas consists of different gas species, which are assumed to be well-coupled and share the same temperature and velocity . Each species has a number density (do not confuse this with of the particles) and a mean molecular weight .

The model also includes different chemical reactions with their respective net reaction rates of reaction , which are the number of reactions per unit time per unit volume. In chemical equilibrium the net reactions rates are equal to zero (which does not mean that no chemical reaction takes place).

Every reaction costs or sets free energy. The definition is such, that a positive reaction rate sets free the energy . That means that the total energy per unit time and unit volume set free by reaction is . With chemistry the continuity equations of the different gas species are

(24)

where is the creation rate of gas species due to reaction (i.e., the change in number density due to reaction k). The component of the creation rates accounts for changes due to non-chemical processes, e.g. evaporation since the mass of the evaporated material has to be added to the gas.

The total mass loss per unit time and unit volume of all particle populations is

(25)

This mass has to be added to the different gas species via the component of their creation rates

(26)

where is the fraction of the evaporated mass added to gas species . For mass conservation has to hold.

To close the system of equations one has to solve for the gas temperature , the gas velocity and the number densities . To do so we use Euler’s momentum equation (11) and sum up the total momentum of the gas and the particles

(27)

The first term in the sum over is the momentum carried by the gas, the second term is the pressure according the ideal gas law. The sum over refers to the momentum carried by the particles. Since the particles are pressureless there is no pressure term. This formula can be manipulated to get to

(28)

(see B for detailed calculations).

The same kind of calculations can be done for Euler’s energy equation (12)

(29)

The first term in the sum over is the kinetic energy carried by the gas, whereas the second term is the sum of the internal energy of the gas and the pressure term according the ideal gas law. is the number of degrees of freedom of gas species . The first term in the sum over is the kinetic energy carried by the particles, the second term is the thermal energy of the particles. In addition to that the radiative flux has to be included. The right hand side is not equal to zero, because energy can be created or consumed by chemical reactions. This equation can be further manipulated to

(30)

(see B for detailed calculations).

Equations (28) and (30) are a set of two coupled differential equations, which can be solved for and . Together with the continuity equation (24) of the gas and the particle differential equations (17), (20), (22) and (23) this closes the system of equations.

Setting and using the gas species , , this reduces in principle to the model of Desch and Connolly (2002), but correcting for some sign errors and implementing the smooth transition of the evaporation temperature.

2.3 Numerical Method

The system of equations (17), (20), (22), (23), (24), (28) and (30) is extremely stiff. Therefore it would require a very small step size and a very large number of grid points to numerically integrate it with an explicit integration scheme. We used here the implicit integrator DVODE (Brown et al., 1989) to integrate the equations simultaneously through the computational domain.

We performed our calculations in the comoving frame of the shock front, which is set to be at . At the shock front the gas parameters are changed according the Rankine-Hugoniot jump conditions, while the particle parameters remain unchanged.

After every complete integration the radiative transfer calculations have to be done again with the new particle and gas parameters. Then the integration is repeated with these new radiative parameters. This is repeatedly done until convergence is reached.

Another crucial point is the calculation of the post-shock boundary temperature , which is needed for the radiative transfer calculation. Desch and Connolly (2002) calculated by using incorrect jump conditions adopted from Hood and Horanyi (1991). This was corrected by Morris and Desch (2010) by calculating their own jump conditions. They found post-shock temperatures on the order of and concluded that chondrule formation is not possible in strictly one-dimensional models. Since disks are not one-dimensional objects they will eventually cool by radiation. Therefore Morris and Desch (2010) loosened the one-dimensional assumption by setting the post-shock temperature to the pre-shock temperature . This approach raises a problem, however, since it forces the gas and dust to radiatively cool through the downstream boundary. Since this downstream boundary is not a physical boundary, but just a computational boundary, this does not appear to be justified. The cooling is then dependent on the distance between shock front and post-shock boundary.

In our model we do not set the post-shock temperature a-priori since in a one-dimensional stationary model all downstream parameters are completely set by the upstream conditions. Therefore we perform right before the first iteration an additional integration without radiation. The temperature reached there at the post-shock boundary is then used as a first approximation of the post-shock temperature. After every further iteration, we check the radiative flux through the post-shock boundary. If the flux is positive i.e., the final temperature is higher than the post-shock temperature, we increase slightly and vice versa. If convergence is reached then the radiative flux at the post-shock boundary is equal to zero. We want to point out that is only true at the boundaries, which have to be at large enough .

The post-shock temperatures we found by this approach are even higher than those calculated by Morris and Desch (2010). We think this is due to some approximations performed in deriving their jump conditions. Later in section 3 we introduce a method to perform simulations with a vertical energy loss.

3 Results

For reasons of comparison we used the same input parameters as Desch and Connolly (2002), which we want to repeat here.

We have one particle population with initial radius of 300 m and an material mass density of 3.3 g/cm. Their initial number density can then be calculated by assuming a gas density of g/cm, a dust-to-gas ratio of 0.005 and assuming that 75 % of the dust mass are in the chondrule precursors. The particles have a heat capacity of erg/g/K. Between 1400 K and 1820 K melting takes place. This results in an effective heat capacity of erg/g/K within this temperature interval (see Desch and Connolly, 2002, and references therein for detailed descriptions). The latent heat of evaporations is erg/g.

The intial temperature of both gas and dust is 300 K. The gas consists of four species: atomic hydrogen, molecular hydrogen, helium and silicon monoxide (SiO). The only chemical reaction we consider is hydrogen dissociation and recombination, which consumes 4.48 eV for every H molecule that breaks apart. We used the reaction rates given in Desch and Connolly (2002) adopted from Cherchneff et al. (1992). The reaction can also go in reverse direction; then setting free energy. If the gas exceeds 1350 K, then the dust associated with it evaporates. The mass of the evaporated dust gets added to the silicon monoxide.

3.1 The standard shock

Figure 3: Left: Particle temperature in a shock with shock speed km/s. Right: Same as left but with gas temperature. See text for further description.

A standard shock with a shock speed of km/s is shown in figure 3. Most of the time the gas and the particles are well-coupled and share the same temperature. Already 3-4 hours before the shock front, the particles receive radiation from the hot gas and particles behind the shock. The temperature increases rapidly until roughly 1350 K. From that point on the temperature only slightly increases towards the shock front. This change in slope is related to the opacity (see figure 1). As the gas reaches temperatures of 1350 K the fine grained dust associated with the gas gets evaporated and the opacity drops by one order of magnitude. From that point on the opacity is purely caused by molecular lines. Therefore the pre-heating region can be sub-divided into an optically thick and an optically thin region. It is important that the opacity is not set to zero because then the gas could not actively cool anymore. The only chance to lose energy would then be thermal contact with the particles. This would artificially slow down the cooling process.

At the position of the shock front the parameters of the gas only are changed (i.e., without radiation or dust) according the Rankine-Hugoniot jump conditions wheareas the particle parameters remain unchanged. Therefore, the gas reaches temperatures as high as 2700 K. The particles suddenly find themselves surrounded by gas that is much hotter than before the shock and that has velocities much smaller than those of the particles. The particles are now heated up by thermal contact with the hot gas and by frictional heating. But since the radiative cooling into the pre-shock region is relatively effective, both the gas and the particles quickly reach an equilibrium post-shock temperature before the particles are able to adapt to the high gas temperatures. In that way the maximum temperature the particles reach is K smaller than the maximal gas temperature.

The constant post-shock equilibrium temperature is a consequence of the one-dimensionality of the model. As soon as the gas and the particles are a few optical depths behind the shock front there is no way for them to lose energy into the pre-shock region by radiation. Since the model consists of infinitely extended plane-parallel layers, they can not cool in lateral directions.

Therefore our model has in principle two different jump conditions: one is local at the position of the shock front, where only the gas parameters are changed according the Rankine-Hugoniot jump conditions. The other global jump condition forces the gas and particles to be at rather high post-shock temperatures after the relaxation of the temperature spike at the shock front.

From now on the maximum temperature the particles experience in the spike is referred to as peak temperature , whereas the equilibrium post-shock temperature is .

The results for different shock speeds look similar. The lower the shock speed, the lower and and the smaller the pre-heating region. Only if the shock speed is too low for the gas to reach 1350 K, then there’s also a lack of the “opacity knee” noticable.

The chemical reaction of hydrogen dissociation works as energy sink in this case. Instead of raising the temperature, energy is consumed in breaking molecular bonds. Simulations without chemical reactions show an overall increase in temperature downstream by K. We want to point out that only a small fraction of the H is dissociated. This demonstrates that it is very important to include chemistry in such simulations: the temperature deviation can make the difference between losing particles via evaporation or not. Whether or not other chemical reactions are equaly important has to be investigated.

Figure 4: Cooling rates of the particles at particle temperatures T K in the standard shock

Figure 4 shows the cooling rates of the particles in the standard shock at temperatures T K which is the important temperature regime for crystallization. As seen here the peak of the cooling rates is at the upper limit of what is allowed by experimental constraints. The standard shock itself is too weak to completely melt the particles. In addition to that the cooling rates drop to zero at already  1600 K. This is discussed later in this section.

Figure 5: Velocity of the particles.

Another interesting feature is shown in figure 5: the velocity of the particles in the simulation. As soon as the pre-heating sets in the particles’ velocity decreases because the gas velocity decreases. From initially km/s to roughly km/s shortly before the shock front. After the gas velocity is changed according the Rankine-Hugoniot conditions at the shock front the particles are rapidly decellerated by the drag force within minutes to km/s.

The Mach number of the gas is initially of the order of 5. But due to the increase in temperature and decrease in velocity it is only roughly 2 at the position of the shock, where the Mach number is applied in the Rankine-Hugoniot conditions.

Figure 6: The peak and post-shock temperatures of the standard shock for different shock speeds. The upper solid line denotes a temperature of 1820 K, which we assume the particles have to reach to be completely molten. The lower solid line is at 1400 K. This is the temperature the particles have to reach at least after the shock to retain their volatiles. This is not possible in the standard case.
Figure 7: Same as figure 6 but with higher gas mass densities of 10 g/cm (left) and 10 g/cm (right).

The problem of chondrule formation within this scenario is summarized in figure 6. If we assume that the particles have to reach a temperature of at least 1820 K to be completely molten and at the same time to cool rapidly back down below at least 1400 K to retain volatile materials which are observed in them, this is not possible in a one-dimensional model of an optically thick disk.

If we increase the gas mass density (see figure 7) the problem improves slightly and the two lines of peak and post-shock temperature are further apart. But only at unrealistic densities of g/cm we could have both requirements fulfilled at the same time. According the disk model by Bell et al. (1997) these high midplane densities can only be found at radii of AU or closer to the sun, where the temperatures are already on the order of 10 K to 10 K depending on the accretion rate.

But we want to point out that post-shock temperatures of 1400 K might already be too high to retain the volatiles. In addition to that the particles are already at high temperatures in the pre-heating phase for a prolonged time ( 3 h for the standard shock, cf. figure 3), which is already too long to be consistent with chondrule formation (Fedkin et al., 2012).

3.2 Vertical energy loss

It is clear that a one-dimensional, plane-parallel model does not match the situation in actual protoplanetary disks, since such disks are not one-dimensional objects. They have a vertical extent with decreasing densities at higher altitudes above the midplane. At AU (the region of the asteroid belt where most of the chondrules can be found today) the pressure scale height of a disk around a solar mass star with a midplane temperature of K is AU, with the mean molecular weight of the gas amu and the proton mass .

If the disk is heated up by a shock whose propagation direction lies in the plane of the disk it can cool in vertical direction by radiative diffusion. But to cool down to the pre-shock temperature the gas and particles have to travel at least a distance comparable to the disk height. This is because in a vertically optically thick disk the radiative diffusion is a photon random walk: before it finds its way up a distance , it has an equal chance of moving a distance downstream. Radiative diffusion cooling will thus not create a steeper temperature gradient in the downstream (in-plane) direction than in vertical direction where the radiation escapes. A downstream cooling length of at least (assuming AU) amounts, with a gas velocity of 6 km/s, to at least one month.

To estimate this we added an energy loss term to our equations

(31)

We transferred this energy loss into a temperature change and included it into the differential equation of our gas temperature

(32)

where the denote the terms in the standard equation. The parameter is a length scale, which determines the strength of the energy loss. The larger , the smaller is the energy loss. Therefore is also a rough estimate of the vertical extent of the disk assumed. With time the temperature should approach the pre-shock temperature. To get a definitive answer one has to perform fully three-dimensional radiative transfer calculations. But this is beyond the scope of this paper.

Figure 8: The standard shock (black) compared to simulations with a vertical energy loss with different loss parameters of 4 % (red), 2 % (green) and 1 % (blue) of the pressure scale height AU at 3 AU. The shock speed is 6.5 km/s in all cases.

We have done this for different length scale parameters essentially simulating disks in which all the dust has settled into a very thin midplane layer. The results are compared to the standard case, which corresponds to , in figure 8. As seen here, the higher the energy loss (the lower the length scale ), the smaller are the pre-shock regions and the peak temperatures. The length scale parameters chosen here are extremely small, only a few percent of the disk’s pressure scale height. It is very questionable if such disks exist. And even in those thin disks the particles are cooked at temperature above 1400 K for hours. The shock speed here is not even enough to completely melt the particles. With the vertical energy loss the particles still spend 2-3 hours in the high-temperature pre-shock phase. At higher shock speeds and therefore temperatures the situation is even worse.

3.3 Optically thin case

To investigate an optically thin case, where the particles can freely lose energy by radiation, we did another run where we set the mean intensity to at every position. This means the gas and the particles are always in a radiation field with an ambient temperature of . At the shock front the gas parameters are changed as usual and the particles can adapt to it. This could correspond for example to bow shocks created by planetesimals on eccentric orbits (Hood, 1998). Here we assume the shocked volume of the disk is small compared to the unshocked medium such that the particles mostly see the radiation field from the unshocked gas. Of course, this assumption breaks down close to the shock front, where fully three-dimensional calculations are needed. But it could also represent the case for which the shock loses all opacity due to dust evaporation.

Figure 9: Particle temperature of an optically thin shock with a shock speed of km/s. See text for details.

The result is shown in figure 9. As expected there is no pre-heating. Therefore the destination temperature for the jump conditions is lower and therefore also the target temperature. To reach the melting temperature in an optically thin case, we needed to increase the shock speed to km/s, instead of km/s in the optically thick cases.

Right after the shock the gas and particles cool down rapidly and approach asymptotically the ambient temperature. The particles are only for a few minutes at critical temperatures and are back below 500 K after h.

Figure 10: Cooling rates at in the optically thin simulation (solid line) and with the Planck-mean opacity decreased by a factor of 10 (dashed line). The discontinuity at is due to the change of the heat capacity regime introduced earlier in section 3.

Unfortunately, the cooling rates (figure 10, solid line) are at least two orders of magnitude too high in the crystallization regime at T K to produce the observed chondrule textures (see e.g. Hewins et al., 2005; Desch et al., 2012). Since the Planck mean opacity is always an upper limit on the opacity we arbitrarily decreased it by a factor of to investigate the effect of lower opacities. The gas has then a lower ability to directly cool by radiation. The major cooling channel is then by thermal contact with the particles, which slows down the overall cooling process. The result is shown in figure 10 (dashed line). The cooling rates are still too large to be consistent with chondrule formation. Decreasing the opacity even further does not have any effect, since it is already low enough that the cooling via the particles is the dominant process.

In the bow shock scenario the optically thin approximation breaks down close to the shock front. Therefore, in reality there will be a pre-heating just before the shock front and less cooling after the shock. Whether this can produce the desired cooling rates has to be investigated by fully two-dimensional simulations. Further – but still one-dimensional – simulations by Morris et al. (2012) suggest that this could be the case.

4 Conclusions

Our 1-D radiative nebular shock model treats the variables at the downstream boundary as output of the model instead of boundary conditions that can be set. The iteration of the model then automatically finds the right downstream state of the matter, given the upstream boundary conditions. We find that this procedure prevents a post-shock slow (few minutes) cooling process. Instead, after the temperature spike and super-rapid cooling right after the main shock, the temperatures stay virtually constant. Only at distances from the shock comparable to the scale height of the disk will the 1-D approximation break down and will sideways (i.e. upward and downward) cooling set in. In our model we mimic this with a simple sideways cooling term.

For cases where the shock structure is local, for instance the bow shock of a planetesimal (see Morris et al., 2012), the shock scenario for chondrules might work because then the 1-D plane parallel assumption breaks down and sideways cooling can commence in a matter of hours. If the disk is optically thin then the cooling rates are too high to be consistent with the constraints on chondrule formation.

We conclude that while the nebular shock model for chondrules may work for local shocks (e.g. planetesimal bow shocks) the scenario has difficulties for global shocks in an optically thick nebula. This is because such global radiative shocks do not produce sufficient cooling after the temperature spike to be consistent with meteoritic constrains.

5 Acknowledgements

This work has been supported by the Deutsche Forschungsgemeinschaft Schwerpunktprogramm (DFG SPP) 1385 “The first ten million years of the solar system”.

We would like to thank Guy Libourel and an unknown referee for their fruitful comments on this paper.

Appendix A Drag force/frictional heating

If the particles are in a gas flow, which has a different velocity, they are accelerated or decellerated, respectively, by exchanging momentum with the gas. During this process they also get heated up by frictional heating of the gas molecules, just like a spacecraft is heated up while re-entering the Earth’s atmosphere. This was described by Gombosi et al. (1986, and references therein). If the particle radius is smaller than the mean free path of the gas molecules – which is valid up to a gas mass density of the order of g/cm for millimeter-sized particles (Desch and Connolly, 2002) – the drag coefficient used in equation (14) is given by

(33)

The parameter is the absolute value of the difference of the particle and the gas velocity measured in units of sound speeds

(34)

with the mean molecular weight . The error function is defined as

(35)

The recovery temperature used in equation (18) is given by

(36)

If , which means particles and gas share the same velocity, then , which means the particles adapt the gas temperature with time. The heat transfer coefficient is given by

(37)

Appendix B Hydrodynamic calculations

Equations (27) and (29) have to be extensively manipulated to get to the two coupled differential equations (28) and (30) of the gas temperature and velocity. Since this calculus is somewhat obscure we want to present it here in all details.

Equation (27), which represents the momentum conservation consists of the derivative of three terms. We want to do the derivative for every term separately beginning with

(38)

The first term here can be replaced using the continuity equation of the gas (24) while the second term is already in the needed form only replacing the sum with the total gas mass density

(39)

Similar transformations can be applied to the second term of equation (27)

(40)

As with equation (38) the continuity equation of the gas (24) can be used to manipulate the first term, while the second is already in its final form

(41)

The third term of equation (27) with the particle momentum can be transformed as

(42)

The first term is equal to zero because of the continuity equation of the particles (13) and the second term can be replaced by the drag force (17). The change in mass in the third term can be transformed to a change in particle radius

(43)

Adding all terms up this leads to equation (28).

Similar transformations can be applied to equation (29). We do them here for all terms in the sums separately. The first term in the sum over represents the kinetic energy of the gas due to its flow

(44)

While the second term of this result is already in its final form, the first term can be replaced by the continuity equation of the gas (24)

(45)

The second term in the sum over in equation (29) is the combined internal energy of the gas and the pressure. It can be transformed as follows

(46)

where again the continuity equation of the gas (24) was used.

The first term in the sum over in equation (29) – the kinetic energy of the particles – can be written as

(47)

The first term here is again equal to zero because of equation (13). The second term can be replaced by the drag force (17), while the third term can again be transformed into a derivative of the particle radius

(48)

The second term in the sum over in equation (29) is the internal energy of the particles

(49)

where again the continuity equation of the particles (13) was used. The remaining term can be further transformed to

(50)

The derivative of the radiative flux is already given in equation (9). Summing up every term, this yields equation (30). The formulas for and can be inserted from equations (20) and (22). Equations (28) and (30) are two coupled differential equations of the form

with the respective and . Decoupled they can be written as

Footnotes

  1. journal: Icarus

References

  1. Bell, K. R., Cassen, P. M., Klahr, H. H., Henning, T., 1997. The Structure and Appearance of Protostellar Accretion Disks: Limits on Disk Flaring. The Astrophysical Journal 486, 372 – 387.
  2. Boley, A. C., Durisen, R. H., 2008. Gravitational Instabilities, Chondrule Formation, and the FU Orionis Phenomenon. The Astrophysical Journal 685, 1193 – 1209.
  3. Boss, A. P., 1996. A concise guide to chondrule formation models. In: Chondrules and the Protoplanetary Disk. Cambridge University Press, pp. 257 – 263.
  4. Brown, P. N., Byrne, G. D., Hindmarsh, A. C., 1989. VODE: A variable-coefficient ODE Solver. SIAM Journal on Scientific Computing 10, 1038–1051.
  5. Cherchneff, I., Barker, J. R., Tielens, A. G. G. M., 1992. Polycyclic aromatic hydrocarbon formation in carbon-rich stellar envelopes. The Astrophysical Journal 401, 269 – 281.
  6. Ciesla, F. J., 2005. Chondrule-forming processes – An Overview. In: Chondrites and the Protoplanetary Disk. Astronomical Society of the Pacific, p. 811.
  7. Ciesla, F. J., Hood, L. L., 2002. The Nebular Shock Wave Model for Chondrule Formation: Shock Processing in a Particle-Gas Suspension. Icarus 158, 281 – 293.
  8. Desch, S. J., Connolly, Jr., H. C., 2002. A model of the thermal processing of particles in solar nebular shocks: Application to the cooling rates of chondrules. Meteorites & Planetary Science 37, 183–207.
  9. Desch, S. J., Morris, M. A., Connolly, Jr., H. C., Boss, A. P., 2012. The importance of experiments: Constraints on chondrule formation models. Meteorites & Planetary Science 47, 1139 – 1156.
  10. Fedkin, A. V., Grossman, L., Ciesla, F. J., Simon, S. B., 2012. Mineralogical and isotopic constraints on chondrule formation from shock wave thermal histories. Geochimica et Cosmochimica Acta 87, 81 – 116.
  11. Gombosi, T. I., Nagy, A. F., Cravens, T. W., 1986. Dust and neutral gas modeling of the inner atmospheres of comets. Reviews of Geophysics.
  12. Hewins, R. H., Connolly, H. C., Lofgren, G. E., J., Libourel, G., 2005. Experimental Constraints on Chondrule Formation. In: Chondrites and the Protoplanetary Disk. Vol. 341. ASP Conference Series, p. 286.
  13. Hood, L. L., 1998. Thermal processing of chondrule precursors in planetesimal bow shocks. Meteoritics & Planetary Science 33, 97 – 107.
  14. Hood, L. L., Horanyi, M., 1991. Gas dynamic heating of chondrule precursor grains in the solar nebula. Icarus 93, 259 – 269.
  15. Iida, A., Nakamoto, T., Susa, H., Nakagawa, Y., 2001. A Shock Heating Model for Chondrule Formation in a Protoplanetary Disk. Icarus 153, 130 – 150.
  16. Kley, W., Nelson, R. P., 2012. Planet-Disk Interaction and Orbital Evolution. Annual Review of Astronomy and Astrophysics 50, 211–249.
  17. Mihalas, D., Weibel-Mihalas, B., 1999. Foundations of Radiation Hydrodynamics. Dover.
  18. Morris, M. A., Boley, A. C., Desch, S. J., Athanassiadou, T., 2012. Chondrule formation in bow shocks around eccentric planetary embryos. The Astrophysical Journal 752.
  19. Morris, M. A., Desch, S. J., 2010. Thermal Histories of Chondrules in Solar Nebula Shocks. The Astrophysical Journal 722, 1474 – 1494.
  20. Semenov, D., Henning, T., Helling, C., Ilgner, M., Sedlmayr, E., 2003. Rosseland and Planck mean opacities for protoplanetary discs. Astronomy & Astrophysics 410, 611–621.
  21. Villeneuve, J., Chaussidon, M., Libourel, G., 2009. Homogeneous Distribution of Al in the Solar System from the Mg Isotopic Composition of Chondrules. Science 325, 985 – 988.
  22. Yu, Y., Hewins, H. H., Zanda, B., 1996. Sodium und Sulfur in Chondrules: Heating Time and Cooling Curves. In: Chondrules and the Protoplanetary Disk. Cambridge University Press, pp. 213 – 219.