Fermi National Accelerator Laboratory

FERMILAB-Pub- 90/i 1 I-A

May 1990

THE DYNAMIC INSTABILITY OF ADIABATIC BLASTWAVES

DONGSU RYU

NASA/Fermilab Astrophysics Center, Fermilab

and

ETHAN T. VISHNIAC_

Department of Astronomy, University o£ Texas

ABSTRACT 7-j 

Adiabatic blastwaves, which have a total energy injected from the center,

E o< _q, and propagate through a preshock medium with a density, PE ocr -W,

are described by a family o£ similarity solutions. Previous work has shown that adiabatic blastwaves with increasing or constant postshock entropy behind the shock

front are susceptible to an oscillatory instability, caused by the difference between

the nature of the forces on the two sides of the dense shell behind the shock front.

This instability sets in if the dense postshock layer is sufllciently thin. In this paper, we consider the stability of adiabatic blastwaves with a decreasing postshock

entropy. Such blastwaves, if they are decelerating, always have a region behind the

shock front which is subject to convection. Some accelerating blastwaves also have

such region, depending on the values of q, _, and 7 where 3, is the adiabatic index.

However, since the shock interlace stabilizes dynamically induced perturbations,

blastwaves become convectively unstable only if the convective zone is localized

around the origin or a contact discontinuity far from the shock front. On the other

hand, the contact discontinuity of accelerating blastwaves is subject to a strong

Raylcigh-Taylor instability. The frequency spectra of" the nonradial, normal modes

of adiabatic blastwaves have been calculated. The results have been applied to the

shocks propagating through supernovae envelopes. We show ti_at the metal/He and

He/H interlaces are strongly unstable against the R_ylelgh-Taylor instability. Tt_s

instability will induce mixing in supernovae envelopes. In addition we discuss the

implications of this work for the evolution of planetary nebulae.

,qtdJ]ect headingJ: hydrodynamics - instabilities - interstellar: matter - nebulae:

planetary - shock waves - stars: supernovae

Submitted to Ast_phyaica! Jo=rna!

PresidentialYoung Investig_or and Sloan Fellow

(NA_A-CR-IG2498) THE _YNA_[C I_STA_ILITY OF

ADIABATIC BLA_TNAVLS (Fermi Ndtionnl

Accelerator LAD.) 39 p CSCL 20_

@

N?O-Z3OT?

Unclds

G3/3_ 0281035

_p_r_dbyUn_erI_e$Re_e_r_hA=$_a_i_n_n_under_nt_a_i_h_he_ni_eds_a_e$Depa_men_fEnergy
I. INTRODUCTION

An understanding of the propagation of shock waves in astrophysical systems

can be essential for understanding the history and present state of such systems.

A blast shock wave (or a blastwave, in short) isgenerated, by an object releasing

energy, whether impulsively, as in a supernova explosion, or continuously, as in a

stellaror galacticwind. The resultingbiastwave can be described by a Sedov-Taylor

type similaritysolution,provided the totalinjected energy E isa power law in time,

and the density of the ambient medium isa power law in radius,

= po (1.2)

(Sedov 1946, 1959; Taylor 1950; Ostriker and McKee 1988). The similarity formallsm demands that the flow quantities depend only on a single dimensionless

parameter involving length and time scales. Such a similarityformalism provides

an excellent approximation to realisticsystems, when it is applied in the region

far from the spatial and temporal boundaries and when it contains an adequate

description of physics.

The stabilityof blastwaves has been a longstanding problem to physicistsand astrophyslcistsalike.Erpenbeck (1962) studied the stabilityof a steady plane-parallel

shock dividing space into two homogeneous regions and showed that it is not exponentlally unstable against a rippling of the shock front. Subsequently, several

workers, including Bertschinger (1986), demonstrated that the shock interface is

stable against dynamically induced perturbations because of the stabilizingtangentialvelocity in the postshock flow,produced by the obliquity of the rippled shock.

However, the contact discontinuitybetween the ejected material and the swept-up

gas resembles the classicalillustrationof a Rayleigh-Taylor instabilityin accelerating blastwaves. In fact, Bernstcin and Book (1978) showed that this contact

discontinuityin iaer,tropicblastwaves isindeed unstable.

Ryu and Vishnlac (1987, 1988) analyzed the stabilityof adiabatic blastwaves

propagating through a uniform medium (_a = 0) and derived the eigenvalues and

cigenfunctions for the nonradial, normal modes. Such blastwaves are susceptible

to an oscillatoryinstability,due to the differencebetween the nature of forces on

the two sides of the dense shellbehind the shock front. They become subject to

growing oscillationsif the shellis sufficientlythin. The physical mechanism for

this instabilitywas firstsuggested by Vishniac (1983) for an isothermal shellwith

infinitesimalthickness and laterconfzrmed by Bertschlnger (1986). In a subsequent

paper, Vishniac and Ryu (1989) gave an analyticalderivation of the eigenmodes of
3

isothermal blastwaves propagating through a uniform medium under the approximation that the evolutionary time scale of the perturbations is much shorter than

the evolutionary time scale of the shock itself. From the above work, it can be seen

that adiabatic blastwaves with increasing or constant postshock entropy from the

shock front, or with

6-2q (1.3)

- + - 1)q

where _, is the adiabatic index (see §II), are overstable for _/ sui_clently close to

1. For a realistic shock this condition is usually equivalent to requiring that the

postshock gas cool ei_ciently. In the absence of radiative cooling such shocks will

not be subject to these growing oscillations.

However, as pointed out by Bandiera (1984), some adiabatic blastwaves have a

region behind the shock front which is subject to convection and may be convectively unstable. The stability of such blastwaves was studied by Chevalier (1976)

and Bandiera (1984). Chevalier used the local Rayleigh-Taylor criterion to determine the existence of an instability and the dispersion relation for an incompressible

fluid under constant gravitational acceleration bounded by rigid boundaries to estimate the growth rate of the perturbations. On the other hand, Bandiera used the

local Schwarzschild criterion and a dispersion relation based on a local stability analysis. Unfortunately, in this case, convection is driven dynamically, not thermally,

and such a dynamically driven instability is generally a global instability, not a local

instability. The global stability analysis should include the overall dynamics of the

blastwaves and the correct boundary conditions, in addition to the local structure of

the blastwaves. Therefore the local satisfaction of the convection criterion does not

necessarily guarantee the existence of a global convective instability, and estimates

of the growth rate of perturbations based on a local stability analysis may be misleading. Here, we call this instability a convective instabillty, and reserve the name

Raylelgh-Taylor instability to indicate the instability in the contact discontinuity

(see §II).

In this paper, we consider the stability of adiabatic blastwaves with decreasing

entropy behind the shock front, or

6- 2q (1.4)

> + - 1)q"

They are subject to convective or Raylelgh-Taylor instabiRty depending on the

values of q, w, and 7- This problem is of more than purely formal interest. Realistic

astrophysical shocks often occur in the presence of density gradients and may involve

the injection of energy over time long compared to the evolutionary time scale of

the blastwave. Supernova explosions, where the shock propagates outward through

the envelope of a massive star, constitute an example of the former. The ejection of

planetary nebulae in the later stages of stellar evolution may be an example of both

significant density gradients and continuous energy injection. In this case, earlier

stellar winds may lead to a slgnificant, but rapidly decreasing, density enhancement

around the star. Moreover, the planetary nebula itself is probably driven by a
4

continuous, high-intensity stellar wind. Under these circumstances, instabilities

may result in global inhomogeneltles in the blastwaves at late times. An equally

important effect is the tendency of these instabilities to drive mixing in otherwise

chemically distinct layers in the postshock gas.

In this paper we have calculated thefrequency spectra of the nonradial, normal

modes of adiabatic blastwaves using a method similar to that used in Ryu and

Vishniac (1987, 1988). In §II of this paper we have reviewed the Sedov-Taylor type

similarity solutions for adiabatic blastwaves. In §III the perturbation equations

are derived, and in §IV the boundary conditions at the shock front and at the

contact discontinuity or the origin are discussed. In §V the numerical scheme used

to solve the perturbation equations is described and the results are discussed. In §VI

the results are applied to the shock propagating through the supernova envelopes.

Finally, in §VII we discuss the application of our results to the evolution of planetary

nebulae.

II. BLASTWAVES

We treat the postshock fluid as a one-component ideal ituld with an adiabatic

index "y. The flow is governed by the usual conservation equations:

8p

_ +V.(_') = 0, (2.1)

@,7

P_- _- P(_" _)_" A- _p = O, (2.2)

ap

÷ _. % + _p_. _ = 0. (_.3)

The energy conservation equation can be replaced by the entropy conservation

equation

+_'.e in_- - o (2.4)

except at the shock discontinuity. The shock is assumed to be very strong, so that

the fluid density, velocity, and pressure just behind the shock front are expressed in

terms of the preshock density PE and the shock velocity V, relative to the preshock

medium

P, -- PE_ A"- I1' (2.s)

2v°

vo - -- 1' (2.6)

2 pE Vs2

p. = --, (2.7)

_f+l
by the shock jump conditions. The flow pattern is entirely determined by two

quantities, PE and E, and from them a unique dimensionless parameter of similarity

flow can be formed,

( PE "_1/5

= _\_-_/ (2.s)

_

k _ ) "_-_"

The radius and velocity of the shock front change as

raoc_s--, (2.9)

v, = _-'---_, (2.10)

and the fluiddensity,velocity,and pressure just behind the shock front change as

pj cc _- s-h,

_3a OC _ 5-_o

Psoct- s-v .

(2.11)

(2.12)

(2.13)

It is convenient to introduce the dimensionless fluid density, velocity, and pressure, and the dimensionless radial coordinate for the flow throughout the region

behind the shock front as defined by

: P, (2.14)

P.

'U

'US

Ps

T

F----.

TB

Then, the conservation equations (2.1),(2.2),and (2.4)become

(2.15)

(2.16)

(2.1z)

2 d._ 4

+ + _ = O,

7+Id_ 7+lF

(2.1s)

(_-+-E2 ___ ) d'F+Tq-l_d_ d_ _-iid$ 3-w-q_=o 2-Fq '

-': (____) id_ _+ (_+i_--_)_ Id$ (6-2q+wq "2--T-q-'r,,,) = O.

(2.19)

(2.20)
Blastwaves with w < 5 admit a similarity solution where the radius and velocity of

the shock front remain finite, even though the total mass and energy contained in

a blastwave is infinite if _ > 3.

For q = 0 the above equations can be integrated analytically. The integration is

given in Sedov (1959), and also in Bandiera (1984) in a compact form. Blastwaves

with q = 0 are filled, i.e. the postshock flow extends to the origin, provided

< (7-7)/(7 + 1). In this case, near the origin, the density becomes

ocr _-l , • (2.21)

and the pressure approaches a constant value with

_

_-_ ocr _-1 (2.22)

On the other hand, blastwaves with q = 0 are hollow, i.e. the postshock flow ends

at a contact discontinuity, provided w > (7 - 7)/(7 + 1). In this case, near the

contact discontinuity, the density and pressure become

--re) _-'+_ , (2.23)

o¢ (_ -- _e)_. (2.24)

Here, _¢ is the position of the contact discontinuity. The special case with ¢# --

(7 - 7)/(7 + 1) is sometimes referred as the Primakoff blastwave, and the postshock

fluid quantities are expressed by a simple power law in radius:

oc _, (2.25)

i_ oc _3. (2.26)

The blastwaves with q - 0 can be classified into 7 characteristic types depending

on the values of u; and 7, and the ranges of w for each type are given in Bandiera

(1984). The plots of the postshock flow for each characteristic type of the blastwaves

with 7 = 4/3 are shown in Figure 2 of Chevalier (1976).

For q > 0 equations (2.18), (2.19), and (2.20) can not be integrated analytically,

but the behavior of the postshock flow can be induced from the equations. All the

blastwaves with q > 0 are hollow. Near the contact discontinuity,

-- rc d_ 6 - 2q + wq

cc - 7_ for w < 3

2+q

<0 for w > 3.

The density converges to zero or a constant if

¢0 < 6-2q

- 27+(7-1)q'

(2.27)

(2.28)
and divergesotherwise. Also near the contact discontinuity,

i__ _ 3- _ - q, (2.29)

_a

and the pressure approaches a constant value if _ < 3 and diverges otherwise. On

the other hand, near the shock front,

d-_cc d_ 6(3 ....q) 42V-_(2 + q) [2(2+ 7) + (7 1)q]_, (2.3O)

d$ 2(27-1)(3-q)-47!_-_11(2+q)-[27-(7-1)q]_. (2.31)

(The above relationsalso hold for the cases with q = 0.) The blastwaves with q > 0

can be classifiedinto 5 characteristictypes depending on the values of q, _a,and 7,

and the ranges of _afor each type are following:

figure Ca)

0 - 2q (2.32)

ca< 27+(7-1)q'

figure (b)

6(3 - q) ? - I 4(2 + q) 0-2q <_ <

27+ (7- 1)q 2(2 +7) + (7-- 1)q 7-1- 12(2 +7) +(7- 1)q'

(2.33)

figure(c)

figure (d)

and fiffare (e)

0(3- q) 7- 1 4(2+ q)

2(2 -k 7) -{- (7 -- 1)q 7 + 1 2(2 -{-7) + (7-- 1)q

2(27-1)(3-q) 7-1 47(2+q)

< 27:(7-1)q-7+127-(7-1)q'

(2.34)

2(27 - 1)(3 - q) 7-1 47(2 + q)

27- (7 - 1)q 7+127- (7- l)q

< w < 3 -- q, (2.35)

3-- q < w < 5. (2.36)

Figure 1 shows the plots of the postshock flow for each characteristic type of the

blastwaves with 7 = 5/3 and q = I. Each plot ismeant to illustrativeof a range of

values of ca which yield similar results.

It iswell known that a region issubject to a convective motion ifthe entropy

increases in the direction of gravity (see,for example, Landau and Lifshitz1959).

Denoting by 8 the entropy per unit mass, which isdefined as
8

it can be seen from equation (2.20) that the entropy increases or stays constant

from the shock front toward the origin or the contact discontinuity provided

6 - 2q

< 23" + (3"- 1)q' (2.38)

and decreases otherwise. On the other hand, even though we are neglecting gravity,

the presence of a pressure gradient implies the flow is accelerated and gives rise to

an effective gravity

l d_

gaff - _ d_'" (2.39)

The postshock region around the contact discontinuity is decelerated (or accelerated

in the direction of increasing entropy) provided

< 3-q (2.40)

(see equation (2.29)), and the postshock region behind the shock front is decelerated

provided

2(23"- 1)(3- q) 3"- 1 4_(2+ q)

< 23" - (3" - 1)q 3" + 1 23" - (3" - 1)q (2.41)

(see equation (2.31)). Hence, all or a portion of postshock region is subject to a

convective motion if

6 - 2q

<_

23" + (3" - 1)q

2(23"- 1)(3- q) 3"- i 4_(2+ q) ]

< re=I(3 - q), _- (7- _q - _ _i 24- _4-_)q "

(2.42)

In such region, the local strength of the convective motion may be estimated by the

local growth rate given by Bandiera (1984),

o" _g eff--_rd_ . (2.43)

Prom equations (2.18), (2.19), and (2.20), it can be seen that 0.2 approaches a

constant near the origin in filled blastwaves, or if

7-3"

q = 0 and t_ < _ (2.44)

- 3"+1'

and 0.2 oc 1/(r - re) near the contact discontinuity in the hollow blastwaves, or if

q > 0 or q = 0 and t_ > 7 - 3". (2.45)

3'-t-1

Hence, if convection occurs in the postshock flow, it is expected to be highly localized

in hollow blastwaves but not in filled blastwaves. However, the question whether the

blastwaves locally subject to a convective motion become dynamically unstable is a

different problem, as discussed in §I. Since the shock interface stabilizes dynamically
9

induced perturbations, we expect that only those blastwaves with a highly localized

convective zone far from the shock interface,or with

7--/<_<3 and q=O

7+1

6 - 2q

or < _ < q- 3 and q > O,

2"),+ (_ - 1)q

are convectively unstable.

(2.46)

The accelerating contact discontinuitywith

t#>3-q (2.47)

and the decelerating shock front with

2(27 - 1)(3 -- q) q,-1 4'7(2 -{-q) (2.48)

2_ - (_ - 1)¢ -y+ 12_ - (_ - 1)q

have configurations subject to the Rayleigh-Taylor instability. However, again

since the shock front tends to suppress the instability,we expect that only those

blastwaves with the acceleratingcontact discontinuityare subject to the RayleighTaylor instability,i.e.only those blastwaves where the instabilityislocalized near

the contact discontinuityand far from the shock front.

- In the following sections, we will provide a proof of these arguments and find

the relevant frequency spectra by performing a linear stabilityanalysis.

III. PERTURBATION EQUATIONS

In this paper we willfollow the notation used in Ryu and Vishniac (1987, 1988)

except that w and q have been added to the llstof parameters.

Let 6_, 6v_ and d/_ be the normalized, perturbed variables defined in the Eulerian

coordinates as

6_,(,',0,¢,,t) -- #(,',0,¢,,t) - #o(,',',), (3.1)

Ps

_(,',O,¢,t)- ,,oCr,t)_

6_,(,,0,_,_)-- , (3.2)

,s

6f_(r,O,_,t) -- p(r,O,4_,t) - po(r,t), (3.3)

Ps

where po, vo, and po are the unperturbed fiuid quantities considered in the previous

section. Then, by linearlzingthe hydrodynamic conservation equations (2.1),(2.2),

and (2.3)about the perturbed variables,we obtain the perturbation equations

a6,_ 3' + 1 .... @6,_ Oq_ _ _

+ 15 - _ _op.,_- + _6_ + 2_6p+ _$. 6_+ °_6_, = 0, (3.4) 2 2+q @t 2 8_
10

-=o

where V is the normalized spatial derivativeand _, and 6_'T are the radial and

tangential components of the normalized perturbed °velocity.

Since the unperturbed flow is self-similar,we have seen in the previous section

that the normalized, unperturbed/_uid quantities _, #, and @ depend only on a

single,normalized, dimensionless radius 9. Hence, we assume that the norma/ized,

perturbed variables can be written in a self-similarform consistingof a term which

depends only on the normalized radius,the usual expansion in spherical harmonics,

and a power law in time:

_(r,0,4,t)= _ _z,_(_)y_,_(0, 4)t' (3.z)

l,rrt

l,rn,

_f,(,.,o,_,O- F_.@_(O_,,,(o,_)e,

where the operator V T is defined as

, (3.s)

(3.9)

" 1O^ 1 0.

(3.i0)

Here, s is the dimensionless perturbation growth rate. It depends on I and may

be either real or a complex number° The Erst-order perturbation equations, which

are the partial differentialequations, are reduced to a set of ordinary differential

equations:

( ) [ --] 7 +19 d6_ _d6$, d$ 2_- 7+i 7+15

2 -_--{'P-_-+ d-_+ 9 2 w+ 2 2+qa 6_

d- _-_ -}- 2 6@, -- _(_ + 1)_6_ T "- 0,

(_____2,)_,_,+_-i_@__ _-_@___ _°_

(3.11)

+ [d_d_ 7+13-_-q+7+lS--w 2 2+q 2 2+qs ] _6_,=0, (3.i2)
11

(_ 7+I_ )_dg%T+ [:=_7+13-w-q " 2 dF .2 2+q

"7- 1 161B= 0 ,

2

+'7

7+15-wsl

+ --7_ j A6_T

[ id_ i@ l

+ L- _ + _J _6_,

2 _d_ 2 2+qS 5_

7+2 1_)_ld_ _ 7+2 15-W2+q.s]P=6_p = 0.

(3.13)

(3.14)

Here, we have dropped the subscript Im from 6plm, 6v,lm, 6VTin.,, and 6pl m for

clarity. Although we do not distinguish notationally between, say, 6_(r, 0, _b,t) and

6_(_), it is clear that the latter form is applied in the above equations. In deriving

the equations (3.11) to (3.14), we have used the relations

a6_ 2 + q _ 06_ s 6. (3.15) --_'-= 5--wt a_ +_ p'

-_- = 5- ,,,e a_ + _'_"' (3.1e)

(96_ 2 + q _ a6_ s . (3.17)

IV. BOUNDARY CONDITIONS

a) Outer Boundary Condition8

The outer boundary conditions at the shock frontforthe perturbation equations

(3.11) to (3.14) are derived from the requirement of the mass, momentum, and

energy flux conservations across the shock front:

b",] = O, (4.1)

[p+ =O, (4.2)

[_r] = O, (4.3)

_-'i =0.

Here, ul. is the fluid velocity relative to and perpendicular to the shock front and

u T is the fluid velocity parallel to the shock front. The square brackets denote the

difference in the enclosed quantities across the shock front.
12

Let _Ta be the perturbation in the position of the shock front, ra, and _/a be the

normalized perturbation def_ued as

_.(o,g,,0 - ,7o(0,q,,0 (4.5)

,.°

FoUo_ingequations(3.7)to (3.10),weassumethat _, ca_be written as:

_7.(0,_,0= _. _a,,,_,,,(o,_)_°. (4.6)

l,w'.,

Then, the equations (4.1) to (4.4) give the normalized boundary conditions at the

shock front for the perturbation equations (3.11) to (3.14)

a_= -_,_ - _,_., (4.7)

d_+ 2-%_q'+ 1 _,, (4.s)

6_ -- --_° --w_, + 2 2--_9, + 1 _°, (4.9)

6"_T = --'_o, (4.10)

where the subscript Im has been dropped for clarity. We have assumed that the

ambient medium is unperturbed in deriving the above conditions. We will set _o - 1

for convenience.

b) Inner Boundarp Conditions

The derivation of the inner boundary conditions is subtler than that of the outer

boundary conditions. The physical arguments for the inner boundary conditions

are different for filled and hollow blastwaves. The inner boundary conditions at the

origin for filled blastwaves, i.e., for

q = 0 and _ < 7 - _f (4.11)

--'y+l'

were considered in Ryu and Vishniac (1987). From the requirement that the fluid

should not undergo divergent perturbations in the Lagrangian sense at the origin,

or the fractional density change of the displaced volume element should go to zero

at the origin, the boundary condition

_p = 0 (4.12)

is obtained. The same boundary condition is also obtained from the requirement

that the energy perturbation should vanish at the origin.
13

The inner boundary conditions at the contact discontinuity for hollom blastwaves, 1.e, for

q>O or q=O and w> _,7-7 (4.13)

7+1

were considered in Ryu and Vishniac (1988). They are derived from the requirement

that there are no mass, momentum and energy fluxes across the contact discontinuity, or

u, = 0 (4.14)

at the contact discontinuity and

[p]=o (4.15)

across the contact discontinuity. Here, u, is the fluid velocity relative to and

perpendicular to the contact discontinuity. Let 77c be the perturbation in the

position of the contact discontinuity, re, and _c be the normalized perturbation

defined as

_c(0,_,t) -- _c(0,_,t). (4.16)

Ta

Then, the above two equations give

_,+_c:: 7+I2 5-w(1) k2--C_q _+ _°= o,

(4.17)

@+ _ = o. (4.18)

In deriving the above conditions, it has been assumed that there is no perturbation

inside the contact discontinuity, where a vacuum forms if q = 0 and w > (7-7)/(7+

1) or the sound speed is so large to smooth out immediately all the perturbations

if q > 0. In the case that d_/d_ = 0 near the contact discontinuity, or for

6

q=O w< --

7+1

6 --2q

q>O w<

27 + (7 --1)q'

the second condition becomes

(4.19)

6_=0 (4.20)

which is completely sui_cient as the inner boundary condition. In this case the

other one may be used to compute @c. However, in the case that d_/d_" # 0 near

the contact discontinuity, both conditions must be invoked to obtain a solution.

V. NUMERICAL SOLUTIONS

The _st-order perturbation equations (3.11)to (3.14)with the outer boundary

conditions (4.7) to (4.11) and the inner boundary conditions (4.12) or (4.17) and
14

(4.18) constitute elgenvalue equations for the eigenvalue s, which we have solved

using the following iteration method. First, the derivatives of the normalized,

perturbed variables with respect to the dimensionless radius have been calculated

from the equations (3.11) to (3.14)

dF

-o

where _Q is a vector with

(5.2)

and ._.is a 4 × 4 matrix whose elements are functions of p, v, p, d_ld_, d_ld_, d_ld_,

and _. Second, for given values of q, _, 7, and l, the eigenvalue s has been estimated.

Using the outer boundary conditions at the shock front, the above equation has

been integrated numerically up to the origin or to the contact discontinuity. The

integration has been carried out using the Runge-Kutta-Verner _th-order and sixthorder method. In order to improve the accuracy around the origin or the contact

discontinuity where some quantities diverge, i) ln_, ln_, or ln_ instead of _ has

been used as the independent varlable and il) the unperturbed fluid quantities near

the origin or the contact discontinuity have been replaced by the leading power-law

solutions discussed in §II when possible. Finally, the integrated perturbed quantities

at the origin or at the contact discontinuity have checked against the inner boundary

conditions. If they have not satisfed the inner boundary conditions, a has been

mocllfed to a better value s + As by the formula

As = 6_

(5.3)

or

As =

(5.4)

Here, _6_/8s and 06_;,/c_s have been calculated from

The iteration has been continued unt_ ]As[ < 10 -4, and the resulting s has been

taken as the eigenvalue for given q, w, 7, and l. Typically, s converges in 3 or 4

iterations.

In the stability analysis, the blastwaves with q = 0 can be classified into 4

characteristic types, as discussed _n §II (refer, for instance, equations (2.38), (2.40),
, 15

(2.44), and (2.45)). The ranges of w for which each type applies are given below:

type (a)

3

< - and q = 0, (5.6)

`7

type (b)

type (c)

type (d)

3 7-'7

- < _ < _ and q = O, (5.7) `7 -`7+i

7 -`7 < _ < a and q = o, (5.8)

`7-F1

w > 3 and q = O. (5.9)

The plots of the resulting frequency spectra for each characteristic type of blastwaves

with `7 = 4/3 are shown in Figure 2. Each plot is meant to illustrative of a range

in _ which yields similar results.

The blastwaves with q > 0 can be classified into 3 characteristic types, and the

ranges of _ for which each type applies are given below:

type (a)

6- 2q (5.i0)

(__ 2`7+ (`7- l)q'

type (b)

type (C)

6 - 2q

2`7+ (`7-i)q< _ _< q- 3, (5.11)

_o > q- 3. (5.12)

The plots of the resulting frequency spectra for each characteristic type of blastwaves

with `7 = 5/3 and q = 1 are shown in Figure 3. Figure 3b actually corresponds to the

shock wave structure given in Figure lc, but the frequency spectra corresponding

to perturbations of Figures lb and ld are only slightly di_erent. Consequently, we

have omitted the corresponding plots and allowed the remaining figure to stand for

all three cases.

The frequency spectra

3a are the same frequency

and have been reproduced

growing oscillatory modes

front is sufflciently thin, or

of the 8lied decelerating blastwaves in Figures 2a and

spectra as those in Ryu and Vishniac (1987 and 1988)

here for completeness. Such blastwaves are unstable to

only when the thickness of the shellbehind the shock

for `7sufficientlyclose to one.

Blastwaves with _a in the range given in (5.7) are apparently subject to a

convective motion but have a convective region extending from the origin up to the
16

shock front. The corresponding frequency spectra, such as Figure 2b, indicate that

such blastwaves are actually stable against convection because the shock interface

stabilizes the dynamically induced perturbations. On the other hand, blastwaves

with _ in the ranges given in (5.8) and (5.11) have a convective region highly

localized around the contact discontinuity. The corresponding frequency spectra,

such as Figures 2c and 3b, indicate that such blastwaves are unstable against the

convective instability. In the blastwaves subject to convection, the sound wave

modes at long wavelengths bifurcate progressively into growing modes and decaying

modes. The first mode to split, which we call the primar_ mode, is also the mode

which becomes unstable first. The modes which split at higher l, which we call the

higher harmonic modea, become unstable at progressively shorter wavelengths. The

primary mode is always the most unstable mode. The bifurcating wavelength, _b

at which a given mode splits (_ becomes real), and the critical wavenumber, kc, at

which a given mode becomes unstable (8 becomes positive), depend on the pressure

and density distributions in the postshock shell. The detailed dependence of k b

and kc on the internal structure will not be given here, but will be investigated

in the subsequent paper (Ryu 1990). In order for the convective instability to be

effective, the pressure gradient should be large (see Figure 3b). However, the region

initially subject to convection may become less convective as the instability grows,

since effect of convection will be to reduce the gradient of the entropy. Clearly,

the detailed nonlinear evolution of the convective instability in shocks should be

investigated by high-resolution, three-dimensional numerical simulations.

The accelerating contact discontinuity in blastwaves with _ > 3 - q is unstable

against the Rayleigh-Taylor mode over all wavelengths. Figures 2d and 3c show the

frequency spectra of a growing mode, which is the only non-oscillatory mode we can

identify. Of course, it is expected that there are an infinite number of oscillatory

modes corresponding to sound waves propagating throughout the postshock region.

However, blastwaves with _ > 3 - q are generally stable against such modes over all

wavelengths. We have actually identified several such modes but the real part of s

of such modes is always negative. In the subsequent paper (Ryu 1990), we will show

that, for perturbations with wavelengths smaller than the thickness of the shell, the

growth rate of the ltaylelgh-Taylor instability is given by

r ~ (5.13)

where g is the effective gravity exerted on the shell.

Figure 4 shows the postshock flow of the Primakoff blastwave with -y = 5/3,

q = 0, and _ = 2 and its frequency spectra. This is particularly interesting case,

since this corresponds to a supernova explosion in a preexisting supersonic stellar

wind. Bernsteln and Book (1980) considered the the stability of the Primakoff

blastwave and showed that it is stable against perturbations on all wavelengths.

However, Gai_et (1984a, b) argued that Bernstein and Book had neglected the entropy perturbations and concluded the Primakoff blastwave can be proved to be

stable only against radial modes. Our analysis shows that the Primakoff blastwave
17

is stable also against nonradial, normal modes, even though it is formaUy subject

to a convective motion.

Figure 5 shows the distributions of the perturbed fluid quantities 6_, 6_,, 6_,

and 6_ T of blastwaves which are unstable against the convective instability or the

Rayleigh-Taylor instability. The perturbation wavelength corresponds to

A 2_r

ra l

= 0.628,

(5.14)

which is larger than the thickness of the shells. In both cases, the perturbations are

highly localized around the contact discontinuity.

VI. MIXING IN SUPERNOVA ENVELOPES

The advent of Supernova 1987A in the Large Magellanic Cloud provided the

direct observational evidence that hydrodynamic processes may play an important

role in the early stage of a supernova explosion. A large scale mixing of the

ejecta were suggested from the following observations: i) Both x-ray and 7-ray

were discovered earlier than expected (see, e.g. , Sunyaev et al. 1987; Dota_ui et

al. 1987; Matz et aI. 1988). ii) After a rapid initial decline, the bolometric light

curve increased very slowly and then reached a plateau-like maximum (see, e.g. ,

Catchpole et al. 1987). iii) Broad lines with FWHM v __ 2000 - S000_m s -1 of

Fe II, Ni II, Ar II, and Co II were observed in late 1987 and early 1988 (see, e4.

, Erickson et aI. 1988; Rank et al. 1988; Witteborn et aI. 1989). Several possible

mechanisms for mixing were suggested (see Arnett, Fryxeil, and M_dler 1989 for

more discussion and references), and among them some hydrodynamic instability

has been considered the most probable solution.

Ml_Jler e_ al. (1989) considered the stability of blastwave produced by a point

explosion in a n = S polytrope and applied the stability analysis of an incompre_8-

ible fluid confined between two rigid boundaries discussed in Chandrasekhar (1961).

They concluded that the initial perturbations in their model do not grow significantly in the shock propagation time scale. However, Ebisuzaki, Shigey_, and

Nomoto (1989), and Benz and Thlelemann (1990) considered the stability of blastwave propagating through a realistic progenitor star and applied a local stability

analysis based on the Itayleigh-Taylor criterion of an incompressible ituid (Ebisuzaki,

Shigeyama, and Nomoto) and the Schwarzschild criterion of a compressible fluid

(Benz and Thielemann). Both groups concluded that the two composition interfaces between the hydrogen-rich and helium zones and the helium and metal zones

are eztremeIy unstable. However, as mentioned in §II, an analysis of the dynamic

instability of blastwaves should include the overall dynamics, the boundary conditions, and the local background structure. Therefore, estimates of the growth rate

of the perturbations based on local criteria may be misleading.
18

On the other hand, Nagasawa, Nakamura, and Miyama (1988) performed numerical simulations of supernova explosions in n = 3 and n = 1.5 polytropic stars

using a three-dimensional smoothed particle hydrodynamic (SPH) code and found

the formation of clumpy structures produced by a compressible Rayleigh-Taylor instability. In more accurate numerical simulations using a hydrodynamic code based

on the piecewise parabolic method (PPM), Arnett, Fryxell, and M_iller (1989) considered a supernova explosion in a realistic progenitor star with M = 15._® by

Arnett (1987). They found that the metal/He interface is strongly unstable against

the ltayleigh-Taylor instability and heavier metals penetrate into the He zone. On

the other hand, they found that the He/H interface stays relatively stable.

We have applied the results of §V to a blastwave propagating through the supernova envelope. The metal/He and He/H interfaces in Ebisuzaki, Shigeyama, and

Nomoto (1989) (see their Figures 1 and 2) form contact discontinuities with large

density jumps and appear subject to the Rayleigh-TayIor instability, as suggested

by Chevalier (1976). As a result, they are expected to be strongly unstable over

all wavelengths. However, different models for the supernova progenitor may result

in composition interfaces which are stable against the Rayleigh-Taylor instability.

The stability of the composition interfaces may be tested by

(6.1)

The growth time scale of the perturbations is comparable to the dyn_cal time

scale of the contact discontinuity for perturbation wavelengths comparable to the

thickness of the unstable zone, (at least in the linear regime). Consequently, the long

wavelength perturbations will grow about as fast asthe supernova explosion evolves.

According to the linear stability analysis, perturbations with smaller wavelengths

grow faster. However, perturbations with wavelengths smaller than the thickness

of the unstable zone are expected to be strongly localized around the composition

interfaces. This should prevent them from dominating the long wavelength perturbations. We conclude that the modes that most effectively promote mixing will be

the long wavelength modes, i.e. those with wavelengths comparable to the local

pressure scale height. On the other hand, the reglons between the two composition interfaces and the shock front are expected to be 0nly weakly, or not, unstable

against the convective motions. This is because the small pressure gradient in both

regions makes the convective instability ineffective compared to the Rayleigh-Taylor

instability of the composition interfaces, and also because convective motions in the

region between the He/H interface and the shock front will become stabilized if they

extend to the shock front.

VII. PLANETARY NEBULAE

In the case of a stellar wind driving a shock into its surroundings, it is necessary

to take q _ 1. For planetary nebulae the standard model is that the bright shell
19

of ionized material represents a shock being driven into gas which originates from

an earlier, and gentler, steUar wind emitted by the same star. When the shock

emerges from this neutral material the shock wave will expand rapidly into the low

density interstellar medium and the nebula will eventually become less conspicuous.

At earlier times the shock will encounter a density gradient that depends on the

history of mass loss from the central star. In general, we have _ > 2, the value we

would have for a steady mass loss at a constant velocity.

In an earlier paper Breitschwerdt and Kahn (1990) have pointed out that under

certain reasonable assumptions concerning the evolution of a planetary nebula and

its central star the nebula will be subject to the Rayleigh-Taylor instability. They

pointed out that this may explain the small scale features seen in most planetary

nebulae. This result is consistent with the work we have presented here, but we

wish to make a more general point as well. Taking q >_ 1 and _a > 2 we see

immediately from equation (2.48) that we are always in the regime where the

contact discontinuity is accderating and the dense shell is subject to RayleighTaylor instabilities. In particular, we see that this result does not depend on

the value of "r in the dense shah (or equivalently, on the importance of radiative

cooling). In other words, regardless of the exact modal one uses for the evolution

for planetary nebulae, the implication of this paper is that small scale instabilities

and mixing will necessarily be part of their evolution. Whether this results in large

scale features depends on whether or not there are other stabilizing influences at

work. For example, a sufficiently intense ionizing flux may stabilize large amplitude

ripples in the inner boundary by evaporating Iagglng dense regions, an effect which

is beyond the scope of this paper.

We would like to thank Roger A. Chevalier and Hyesung Kang for comments

on manuscripts. The calculations were carried out using fadlities of the Minnesota

Supercomputer Institute and Fermilab. DR thanks the Minnesota Supercomputer

Institute for the hospitality, where some of this work was done. The work of DR was

supported in part by the DOE and by NASA at Fermilab through grant NAGW1340.

2go_e - During the completion of this paper, we received a preprint from

J. Goodman. He considered the stability of the blastwaves with q = 0 and

(7 - _/)/('), + 1) < _a < 3. We found that his results agree with ours for that

particular case.
2O

REFERENCES

Arnett, W. D. 1987, Ap. 3, 319, 136.

Arnett, W. D., Fryxell, B. A., and Mftller, E. 1989, Ap. 3. (Letters), 341, L63.

Bandiera, R. 1984, Astr. Ap., 139, 368.

Benz, W., and T_elemann, F.-K. 1990, Ap. J. (Letters), 348, L17.

Bernstein, I. B., and Book, D. L. 1978, Ap..L, 22]$, 633.

1980, Ap. J., 240, 223.

Bertschlnger, E. 1986, Ap..L, 304, 154.

BreRschwerdt, D., and Kahn, F. D. 1990, preprlnt.

Catchpole, R. M. et al. 1987, M. N. R. A. S., 229, 15p.

Chandrasekhar, S. 1961, Hydrodynamic and Hydromagnetic Stability (Oxford: Clarendon Press).

Chevalier, R. A. 1976, Ap. J., 207, 872.

Dotani, T. et al. 1987, .Nature, 330, 230.

Eblsuzaki, T., Shlgeyama, T., and Nomoto, K. 1989, Ap. J. (Letters), 344, L65.

Ericks0n, E. F., Haas, M. R., Colgan, S. W. J., Lord, S. D., Burton, M. G., Wolf,

J., Hollenbach, D. J., and Werner, M. 1988, Ap. J. (Letters), 330, L39.

Erpenbeck, J. J. 1962, Phys. Fluids, 5, 1181.

Gaffet, B. 1984a, Ap. J., 279, 419.

. 1984b, Asia. Ap., 135, 94.

Landau, L. D., and Lifshitz, E. M. 1959, Fluid Mechanics (New York: Pergamon).

Matz, S. M., Share, G. H., Leising, M. D., Clmpp, E. L., Vestrand, W. T., PurceLl,

W. R., Strlckman, M. S., and Reppin, C. 1988, .Nature, 331,416.
21

Mfiller, E., Hillebrandt, W., Orlo, M., HSftich, P., MSnchmeyer, R., and Fryxell,

B. A. 1989, Aatr. Ap., 220, 167.

Nagasawa, M., Nakamura, T., and Miyama, S. M. 1988, Pub. Astr. Soc. Japan, 40,

691.

Ost_ker, J. P., and McKee, C. F. 1988, Rev. Mo& Phys., 60, 1.

Rank, D. M. et aL 1988, Nature, 331,505.

Ryu, D. 1990, in preparation.

Ryu, D., and Vishniac, E. T. 1987, Ap. J., 313, 820.

• 1988, Ap. J., 331,350.

Sedov, L. I. 1946, Prikl. Mat. Mekh., 10, 241.

1959, Similarity and DimenJionaI Method8 in Mechanic8 (New

York: Academic).

Sunyaev, R. et al. 1987, Nature, 330, 227.

Taylor, G. I. 1950, Proc. Royal Soc. London, 210A, 175.

Vislmiae, E. T. 1983, Ap. ?., 274, 152.

Vishniae, E. T., and Ryu, D. 1989, Ap..L, 337, 917.

Wittebora, F. C., Bregman, J. D., Wooden, D. H., Pinto, P. A., Rank, D. M.,

Woosley, S. E., and Cohen, M. 1989, Ap. J. (Lettera), 338, LY.
22

FIGURE CAPTIONS

Fig. 1.-- The characteristic types of the postshock flow of the blastwaves with

7 --"5/3 and q - 1. The ranges of w for which each particular type applies are

given in the text.

Fig. 2. m The characteristic types of the frequency spectra of the blastwaves with

7 = 4/3 and q = 0. (a) Theblastwavewithw = 0 which is subject to an

oscillatory instability but stable against the dynwnlcally induced perturbations.

(b) The blastwave with _ = 2.3 which is subject to a convective motion but

stable against the dynamically induced perturbations. (c) The blastwave with

w = 2.7 which is subject to a convective motion and unstable against the

dynamically induced perturbations. (d) The blastwave with _ = 3.2 which is

subject to a Raylelgh'Taylor instability and unstable against the dynamically

induced perturbations. The solid lines show the real eigenvalues 8 and the

dashed Hnes show the real parts of the complex elgenvalues S. In (b) and (c),

at the points where two solid lines meet, a dashed line starts and extends up

to l = 1. Some dashed lines are omitted for the clarity.

Fig. 3.-- The characteristic types of the frequency spectra of the blastwaves with

_/= 5/3 and q - 1. (a) The blastwave with w = 0 which is subject to an oscillatory instability but stable against the dynamically induced perturbations. (b)

The blastwave with w = 1.4 which is subject to a convective motion and unstable against the dynamically induced perturbations. (c) The blastwave with

= 2.5 which is subject to a P_ylelgh-Taylor instability and unstable against

the dynamically induced perturbations. The solid lines show the real eigenvalues s and the dashed lines show the real parts of the complex elgenvalues

8.

Fig.

Fig.

4.-- The Primakoff blastwave with _f = 5/3, q = 0, and 0_ -- 2 (a) and

its frequency spectra (b). The Primakoff blastwave is subject to a convective

motion but stable against the dynamically induced perturbations. The solid

lines show the real eigenvalues a and the dashed lines show the real parts of

the complex elgenvalues a. In (b), at the points where two solid lines meet, a

dashed line starts and extends up to l -- 1. Some dashed lines are omitted for

the clarity.

5.-- The perturbations of the postshock flow. (a) The blastwave with V = 4/3,

q = 0, and w = 2.7, which is convectively unstable, for l -- 10 and a = 0.239.

(b) The blastwave with _f -- 4/3, q = 0, and _ = 3.2, which is Raylelgh-Taylor

unstable, for I = 10 and a - 0.713.
. 23

Author's Postal Addresses

Dongsu Ryu:

NASA/Fermilab Astrophysics Center

Fermi National Accelerator Laboratory

MS 209, PO Box 500

Batavia, IL 60510-0500

Ethan T. Vishniac:

Department of Astronomy

University of Texas

Austin, TX 78712
I I I I

I I I i I I I

0

.8 .85 .9

y = 5/3

q= 1.0

(,.) - O0

.95 1
0'-

1.2

1

.8

.6

.4

.2

0

.8

I I I

m

m

m

m

I

w

m

I

'7 = 5/3-

1.0--

1.05-

.85 .9 .95 1
3

2

1

i I I

' I ' ' ' ' I ' ' ' ' I ' ' ' '

"y = 5/3-

q= 1.0 -

_= 1.4-

m

0.8',,,IIw,ll=,,,IPIIz .85 .9 .95 1
3

2

1

i

n I

l , l l. I t...,J I I I lll I l l I l

0.8 .85 .9 .95 i
3

2

c_

1

m

B

0

.8 .85 .9 .95 I

R/R_
cl:S

l,-0

4

2

0

B

m

m

B

m

n

m

B

I I I I

T = 4/3

q = 0.0

= 0.0

i I I

D

i

-4 k

0

%

%

i i ! i n I

log(l)

I i I I I

I

m

u

m

m

m

m

u

m

_' 'nil lul oil um aul nil am | | um mum anl nm am |

I I I I I n

2 3
60

(i,,)

, , , , I ' ' ' ' I I I i I

- T = 4/3

4_lq = O. 0 q

- _=2.3

2

0

-2

-4

0 I 2 3
420

-2

-4

0

12

3
cO

v

r--,q 0

-2

-4

, , , , I , , I l I , /

[_ "y = 4/3 /

q = 0.0 /

/

i

n

m

m

B

B

m

m

m

0

! I

1

1 1 ' l l

2

I

i

a

m

a

m

m

3

_wJ
c_

V

m.,-m¢

(I)

4

2

0

-2

-4

m

m

m

m

I I I I

"7 = 5/3

q = 1.0

= 0.0

i I I I I i i i I

, I, ,, , I,, ,,

1 2

log(O

I

w

m

m

w

i

m

I

3

m

J

/
4

2

_,= _/3

_ q=l.0

_= 1.4

I

v

0

-2

-4

0 I 2 3

log(Z)
r_

V

(D

0

-2

-4

-__ = _/3 /

_ _=i.o /.

m

m

0

m

I

I

I I I I

log(/,)

I

2

I I I I

I

m

u

m

u

m

w

3
l'eI'''l'''l'7 - 5/3

I

.8

.6

.4

.2

0

0

q = 0.0

c_ = 2.0

I I

J I

.2 .4 .6 .8 1

F,o_rc,.
i,,,,-I

(D

I I i I J I I I 1 i I I I l

- _,= 5/3

4-q=o. °

- a_ = 2.0

2

0

-2

-4 " , I , l l i

0 1 2 3

: ..Z-L

r_f4 i

d
'-0

400

2OO

0

I

_,= 4/3

q = 0.0

_=2.7

l= 10

s = 0.239

.5 .6 .7 .8 .9 1
c.o

30

20

10

0

-10

m

m

i

.7

I I I i

T = 4/3

q = 0.0

= 3.2

l= 10

s = 0.713

_vT/v_

I ,1 I , I , , , , I

.8 .9

n

q

m

i

I

q -:_