# GROWTH OF PERTURBATIONS IN GRAVITATIONAL COLLAPSE AND ACCRETION

## 1 Introduction

Gravitational instability is responsible for a wide range of structures in the observable universe. Much effort has been devoted to structure formation in cosmological models, and it is well-established that the Hubble expansion slows the growth of perturbations (e.g., Peebles 1980). In this paper we are interested in the behavior of nonspherical perturbations in the background of spherical collapse or accretion.

Hunter (1962) studied the collapse of a homogeneous, pressureless, dust cloud, and showed that in the linear regime, perturbations of arbitrary shape and scale grow asymptotically as $`\delta \rho /\rho (t_0t)^1R^{3/2}`$, where $`t=t_0`$ denotes the time of complete collapse, and $`R`$ is the cloud radius. It was thought that this instability might be responsible for the fragmentation of collapsing protostellar clouds. Lin, Mestel and Shu (1965) studied the collapse of a homogeneous ellipsoidal cloud of dust, and showed that its ellipticity increases until a sheet or pancake forms. Homogeneity plays a key role in both examples. However, the presence of a slight initial central concentration significantly alters the evolution of a cloud. Since the dynamical time, $`(G\rho )^{1/2}`$, is shortest in the central region, a cuspy density and velocity profile develops. As a result, one may question the applicability of the Hunter and Lin-Mestel-Shu instabilities to realistic astrophysical situations. Goodman & Binney (1983) have already commented on problems with the Lin-Mestel-Shu instability for inhomogeneous clouds. In this paper, we show that perturbation growth is slowed by the central mass concentration in an inhomogeneous collapse.

A related problem concerns the stability of spherical accretion flow. Bondi (1952) found a class of solutions describing steady-state accretion onto a compact object from a homogeneous medium. There is a unique transonic “critical” flow, for which the mass flux is maximum; the other subcritical solutions describe subsonic flows. Bondi speculated that nature would prefer the critical flow (see Shu 1992 for a discussion) and underlined the importance of a linear stability analysis. After several attempts by a number of authors, the correct analysis was achieved by Garlick (1979) and Moncrief (1980), who showed that both the critical and subcritical flows are globally stable. However, this does not mean that perturbations can not grow spatially. Indeed, we show in this paper that nonspherical perturbations carried by fluid elements are amplified in the supersonic region of the critical flow (see also Goldreich, Lai & Sahrling 1996; Kovalenko & Eremin 1998).

The present study originated from our attempts to understand the origin of asymmetric supernova (Goldreich, Lai & Sahrling 1996; hereafter GLS). A large body of evidence suggests that Type II supernovae are globally asymmetric, and that neutron stars receive kick velocities of order a few hundred to a thousand km s<sup>-1</sup> at birth (see, e.g., GLS, Cordes & Chernoff 1998, and references therein). The origin of the kick is unknown. A class of mechanisms relies on local hydrodynamical instabilities in the collapsed stellar core (e.g., Burrows et al. 1995; Janka & Müller 1994, 1996; Herant et al. 1994) that lead to asymmetric matter ejection and/or asymmetric neutrino emission; but numerical simulations indicate these instabilities are not adequate to account for kick velocities $`>100`$ km s<sup>-1</sup> (Burrows & Hayes 1996; Janka 1998). Global asymmetric perturbations of presupernova cores may be required to produce the observed kicks (GLS; Burrows & Hayes 1996). GLS suggested that overstable g-modes driven by shell nuclear burning might provide seed perturbations which could be amplified during core collapse.<sup>1</sup><sup>1</sup>1See Lai & Qian (1998) and Arras & Lai (1999) for discussion/review on alternative mechanisms which rely on asymmetric neutrino transport induced by strong magnetic fields.

Hints regarding perturbation growth during collapse are obtained by considering the collapse of finite-mass fluid shells onto central point masses. In the context of core-collapse supernova, one might imagine that the spherical shell mimics the outer supersonic region, while the central mass represents the homologous core. A large number of growing modes can be identified (GLS; Appendix B). In particular, the dominant dipole mode grows according to $`\mathrm{\Delta }R(\theta ,t)/R(t)(t_0t)^s`$, with $`s<0`$, where $`R`$ is the mean radius of the shell and $`\mathrm{\Delta }R`$ is the perturbation. For a core-shell mass of order unity, the power-law index $`s1`$, corresponding to $`\mathrm{\Delta }R/RR^{3/2}`$.

The thin-shell model is too simplistic to be applicable to real presupernova collapse. One approach is to determine the stability of the self-similar collapse solution (Goldreich & Weber 1980; Yahil 1983). An analysis by Goldreich & Weber (1980) shows that the inner homologous core is stable against nonradial perturbations. This is not surprising given the significant role played by pressure in the subsonic collapse. Pressure is less important in the supersonically collapsing region, making it more susceptible to large scale instability. A stability analysis of Yahil’s self-similar solution, which extends the Goldreich-Weber solution to include a supersonically collapsing outer core, does not reveal any unstable global mode before the proto-neutron star forms (Lai 2000). However, it is perhaps more illuminating to treat the initial-value problem, and determine how initial perturbations evolve as the collapse proceeds. We carry out such an analysis in this paper.

The rest of the paper is organized as follows. The basic perturbation equations are summarized in §2, and in §3 we derive the asymptotic scaling relations for the perturbations in the regime where the collapse/accretion flow is supersonic. In §4 we present a numerical study of the evolution of perturbations during collapse; this numerical study not only confirms the analytic asymptotic relations, but also explores the regime where pressure is important. We discuss the implications of our results in §5.

## 2 Basic Equations

We consider barotropic fluid obeying the equation of state $`p=K\rho ^\gamma `$, where $`\gamma `$ is the adiabatic index and $`K`$ is a constant. The unperturbed flow is spherically symmetric, with velocity in the radial direction. The Eulerian perturbations of density, $`\rho `$, and velocity, $`𝐯`$, can be decomposed into different angular modes, each of which has the form

$`\delta \rho (𝐫,t)`$ $`=`$ $`\delta \rho (r,t)Y_{lm}(\theta ,\varphi ),`$ (1)   
$`\delta 𝐯(𝐫,t)`$ $`=`$ $`\delta v_r(r,t)Y_{lm}(\theta ,\varphi )\widehat{r}+\delta v_{}(r,t)\widehat{}_{}Y_{lm}(\theta ,\varphi )+\widehat{}_{}\times \left[\delta v_{\mathrm{rot}}(r,t)Y_{lm}(\theta ,\varphi )\widehat{r}\right],`$ (2)   
where

$$\widehat{}_{}\widehat{\theta }\frac{}{\theta }+\frac{\widehat{\varphi }}{\mathrm{sin}\theta }\frac{}{\varphi },$$
(3)   
and $`\widehat{r},\widehat{\theta },\widehat{\varphi }`$ are unit vectors in spherical coordinates. The perturbations of pressure, $`p`$, and gravitational potential, $`\psi `$, have the same angular dependence as $`\delta \rho (𝐫,t)`$. The perturbed mass continuity equation reads

$$\frac{d\delta \rho }{dt}+(𝐯)\delta \rho +\frac{1}{r^2}\left(r^2\rho \delta v_r\right)^{}l(l+1)\frac{\rho \delta v_{}}{r}=0,$$
(4)   
where $`\rho ,𝐯=v\widehat{r}`$ denote the unperturbed spherical flow variables, prime stands for $`/r`$, and $`d/dt=/t+𝐯`$ is the total time derivative. The perturbed radial Euler equation can be written as

$$\frac{d\delta v_r}{dt}+v^{}\delta v_r=\left(\frac{\delta p}{\rho }\right)^{}\left(\delta \psi \right)^{},$$
(5)   
and the tangential Euler equation reduces to

$`{\displaystyle \frac{d\delta v_{}}{dt}}+{\displaystyle \frac{v}{r}}\delta v_{}={\displaystyle \frac{1}{r}}\left({\displaystyle \frac{\delta p}{\rho }}\right){\displaystyle \frac{1}{r}}\delta \psi ,`$ (6)   
$`{\displaystyle \frac{d}{dt}}(r\delta v_{\mathrm{rot}})=0.`$ (7)   
The perturbed Poisson equation is

$$\frac{1}{r^2}\left[r^2(\delta \psi )^{}\right]^{}\frac{l(l+1)}{r^2}\delta \psi =4\pi G\delta \rho .$$
(8)   
The vorticity perturbation reads

$$\times \delta 𝐯=\frac{\delta v_T}{r}\widehat{r}\times \widehat{}_{}Y_{lm}+\frac{l(l+1)}{r}\delta v_{\mathrm{rot}}Y_{lm}\widehat{r}+\frac{1}{r}(r\delta v_{\mathrm{rot}})^{}\widehat{}_{}Y_{lm},$$
(9)   
where

$$\delta v_T\delta v_r\delta u^{},$$
(10)   
with

$$\delta u(r,t)r\delta v_{}(r,t).$$
(11)   
Thus $`\delta v_{\mathrm{rot}}`$ and $`\delta v_T`$ are related to the vorticity of the perturbed flow. Equation (6) is transformed to

$$\frac{d\delta u}{dt}=\frac{\delta p}{\rho }\delta \psi ,$$
(12)   
with the aid of equation (11). Combining equations (5) and (12), we find

$$\frac{d}{dt}\delta v_T=v^{}\delta v_T$$
(13)   
Equations (7) and (13) express the conservation of circulation in a barotropic fluid.<sup>2</sup><sup>2</sup>2For homogeneous collapse or expansion, as in an expanding universe, we have $`v=(\dot{R}/R)r`$, where $`R`$ is the scale factor. Equation (13) then becomes the familiar $`d(R\delta v_T)/dt=0`$. The former also reflects the conservation of angular momentum; following a fluid element, $`\delta v_{\mathrm{rot}}1/r`$. Since, as we shall prove shortly, $`vr^{1/2}`$ as $`r0`$, equation (13) implies the asymptotic relation $`\delta v_Tr^{1/2}`$. We shall focus on irrotational flows from here on. We neglect $`\delta v_{\mathrm{rot}}`$ because it is decoupled from the density perturbation, and $`\delta v_T`$ because it decays inward. The continuity equation (4) for irrotational flow simplifies to

$$\frac{d\delta \rho }{dt}+(𝐯)\delta \rho +\frac{1}{r^2}\left(r^2\rho \delta u^{}\right)^{}l(l+1)\frac{\rho \delta u}{r^2}=0.$$
(14)   
Note that the Poisson equation (8) has the following integral solution:

$$\delta \psi (r,t)=\frac{4\pi G}{2l+1}\left[\frac{1}{r^{l+1}}Q_l(r,t)+r^lS_l(r,t)\right],$$
(15)   
where

$$Q_l(r,t)=_0^rx^{l+2}\delta \rho (x,t)𝑑x,S_l(r,t)=_r^{\mathrm{}}x^{1l}\delta \rho (x,t)𝑑x.$$
(16)   
Equations (12), (14), and (15) determine the perturbed flow. This particular form of the perturbation equations is convenient for implementation as a Lagrangian numerical code (see §4).

## 3 Asymptotic Analysis

Consider a cloud in hydrostatic equilibrium with an initial density profile that decreases outward. As its pressure is depleted, the cloud starts to collapse. Since the dynamical time, $`(G\rho )^{1/2}`$, decreases outward, cuspy density and velocity profiles will be established after the core has collapsed. By contrast, a uniform density dust cloud collapses homologously, and remains uniform as the collapse proceeds. We now study the behavior of the flow and its perturbations in the asymptotic regime where the gas pressure is negligible compared to gravity.

### 3.1 Unperturbed Spherical Flow

Consider how the velocity, $`v_m`$, and density, $`\rho _m`$, of a fluid shell with enclosed mass $`m`$ change as the shell collapses from its initial radius $`r_{m0}`$ to a smaller radius $`r_m`$. The pressure is negligible in the supersonic region, so for $`r_mr_{m0}`$, we have

$$v_m\left(\frac{2Gm}{r_m}\right)^{1/2}r_m^{1/2}.$$
(17)   
Then from the continuity equation, we obtain 

$$\rho _mr_m^{3/2}.$$
(18)   
Note that these relations do not require the mass of the collapsed object to be fixed. Indeed, with $`m_c(t)`$ the mass of the collapsed core at time $`t`$, we have

$$v(r,t)\left[\frac{2Gm_c(t)}{r}\right]^{1/2}r^{1/2},\rho (r,t)\frac{\dot{m}_c(t)}{4\pi r^2v(r,t)}r^{3/2},$$
(19)   
where $`\dot{m}_c`$ is the mass accretion rate onto the core.<sup>3</sup><sup>3</sup>3If the central region of a pressureless cloud is non-singular to begin with, then at the moment when the center reaches infinite density, the central density and velocity profiles are $`\rho r^{12/7}`$ and $`vr^{1/7}`$ (Penston 1969). However, after the core has formed, the profiles are given by equation (19). When the accretion time, $`m_c/\dot{m}_c`$, is much longer than the dynamical time of the flow, equation (19) describes the inner region of a steady-state Bondi flow. In a dynamical collapse, when the accreted mass becomes much larger than the original core mass, dimensional analysis implies $`m_c(t)=K^{3/2}G^{(3\gamma 1)/2}t^{43\gamma }\overline{m}_c`$, where $`\overline{m}_c`$ is a dimensionless number, and $`t`$ is measured from the moment when the center collapses. For $`\gamma =1`$, this reduces to the familiar $`m_c=\overline{m}_c(c_s^3/G)t`$ (see Shu 1977). Equation (19) describes the central region of Shu’s expansion-wave solution (Shu 1977), and the post-collapse extension of the Larson-Penston solution (Larson 1969; Penston 1969; Hunter 1977) in the context of star formation, as well as Yahil’s post-collapse solution in the context of core-collapse supernova (Yahil 1983).

Note that the above asymptotic scaling solution assumes supersonic flow for $`r0`$, i.e., $`vc_s\rho ^{(\gamma 1)/2}r^{3(\gamma 1)/4}`$. This requires $`\gamma <5/3`$. The special case of $`\gamma =5/3`$ is considered in Appendix A.

### 3.2 Perturbations

We investigate asymptotic power-law solutions to equations (12), (14) and (15). Let $`\delta \rho r^a`$ and $`\delta ur^b`$. The Poisson equation has a general solution of the form $`\delta \psi Q_c/r^{l+1}+r^2\delta \rho `$, where the first term arises from a central multipole moment, $`Q_c`$, and the second term is due to the density perturbation outside the the central core. In most astrophysical situations, $`Q_c`$ is zero or close to zero; possible exceptions are discussed separately in §3.3. For example, in accretion onto a star, the supersonic flow may be stopped by a standing shock near the stellar surface. Inside the shock, any inhomogeneity carried in by the gas will be smeared out on a local dynamical timescale. In accretion onto a Schwarzschild black hole, the event horizon defines the inner boundary of the supersonic flow, and the no-hair theorem ensures that mass multipole moments are not retained by the black hole. Thus we have $`\delta \psi r^{a+2}`$.

In the asymptotic regime, $`d/dtv(/r)`$, so equations (12) and (14) reduce to:

$`{\displaystyle \frac{bv\delta u}{r}}+\gamma K\rho ^{\gamma 2}\delta \rho +\delta \psi =0,`$ (20)   
$`\left(a+{\displaystyle \frac{3}{2}}\right){\displaystyle \frac{v\delta \rho }{r}}+\left[b\left(b{\displaystyle \frac{1}{2}}\right)l(l+1)\right]{\displaystyle \frac{\rho \delta u}{r^2}}=0.`$ (21)   
Equation (21) implies $`b=a+2`$ for $`a3/2`$. Equation (20) has the scaling form $`𝒪(br^{b3/2})+𝒪(Kr^{a3(\gamma 2)/2})+𝒪(r^{a+2})=0`$, from which we see immediately that $`b=0`$ and $`a=2`$ (for $`\gamma <5/3`$). To obtain the scaling behavior for $`\delta v_r=\delta u^{}`$, we need a higher order correction for $`\delta u`$. Let $`\delta u=\delta u_0+\delta u_1r^{b_1}`$, where $`\delta u_0`$ and $`\delta u_1`$ are constants independent of $`r`$. For $`K=0`$ (the pressureless case) or $`\gamma <2/3`$ so that $`\delta p/\rho <\delta \psi `$ asymptotically, equation (12) gives $`𝒪(b_1v\delta u_1r^{b_11})𝒪(r^{a+2})`$. We then have $`b_1=a+7/2=3/2`$. For $`K0`$ and $`\gamma >2/3`$, equation (12) reduces to $`𝒪(b_1v\delta u_1r^{b_11})+𝒪(r^{a3(\gamma 2)/2})=0`$, which gives $`b_1=(53\gamma )/2`$. To summarize, the asymptotic scaling relations for the perturbations are:

$`\delta u`$ $``$ $`\delta u_0+\delta u_1r^{b_1},`$ (22)   
$`\delta v_r`$ $``$ $`r^{b_11},\delta v_{}r^1,`$ (23)   
$`\delta \rho `$ $``$ $`2l(l+1){\displaystyle \frac{\rho \delta u}{rv}}r^2,`$ (24)   
where

$$b_1=\{\begin{array}{cc}3/2\hfill & \text{for }K=0\text{ or }\gamma <2/3\text{,}\hfill \\ (53\gamma )/2\hfill & \text{for }K0\text{ and }\gamma 2/3\text{.}\hfill \end{array}$$
(25)   
Note that the above results apply for $`\gamma <5/3`$. The special case of $`\gamma =5/3`$ Bondi accretion is discussed in Appendix A.

### 3.3 Physical Interpretation

Equations (22)-(24) describe the Eulerian perturbations. To derive scaling relations for the Lagrangian displacement, $`𝝃(𝐫,t)=\xi _r(r,t)Y_{lm}\widehat{r}+\xi _{}(r,t)\widehat{}_{}Y_{lm}`$, and for the Lagrangian density perturbation, $`\mathrm{\Delta }\rho `$, we note that $`\mathrm{\Delta }𝐯=\delta 𝐯+(𝝃)𝐯=d𝝃/dt=(𝝃/t)+(𝐯)𝝃`$, which yields

$`\delta v_r`$ $`=`$ $`{\displaystyle \frac{\xi _r}{t}}+v\xi _r^{}v^{}\xi _r,`$ (26)   
$`\delta v_{}`$ $`=`$ $`{\displaystyle \frac{\xi _{}}{t}}+v\xi _{}^{}{\displaystyle \frac{v}{r}}\xi _{}.`$ (27)   
In the asymptotic regime, we find using equation (23) that

$$\xi _rr^{b_1+1/2},\xi _{}r^{1/2}.$$
(28)   
The Lagrangian density and velocity perturbations are

$$\frac{\mathrm{\Delta }\rho }{\rho }\frac{\mathrm{\Delta }v_{}}{v}r^{1/2},\frac{\mathrm{\Delta }v_r}{v}r^{b_11/2}.$$
(29)   
The scaling relations (28) and (29) can also be derived directly from Lagrangian perturbation theory.

Equations (22)-(24) and (29) are the main results of this paper.<sup>4</sup><sup>4</sup>4This result was discussed without derivation in GLS. Similar results were also obtained by Kovalenko & Eremin (1998) in the context of Bondi accretion. They can be understood from the following simple consideration: The time that a fluid element spends near radius $`r`$, of order $`t_r=dt/d\mathrm{ln}rr^{3/2}`$, decreases rapidly as $`r`$ decreases. Unless the specific torque, $`(\delta \psi +\delta p/\rho )`$, increases inward sufficiently rapidly to compensate for the decreased time, the specific angular momentum, $`\delta u`$, of the fluid element will be independent of $`r`$ at small radii. Thus $`\mathrm{\Delta }v_{}r^1`$. This leads to a tangential displacement $`\xi _{}t_r\delta v_{}r^{1/2}`$, which induces a density perturbation $`\mathrm{\Delta }\rho /\rho \xi _{}/rr^{1/2}`$.

Equation (29) indicates that the fractional density and tangential velocity perturbations grow in supersonic collapse/accretion. The radial velocity perturbation, however, can be affected by pressure even in the regime where the pressure has negligible effect on the unperturbed flow. We see that $`\mathrm{\Delta }v_r`$ grows with decreasing $`r`$ when $`\gamma >1`$, and $`\mathrm{\Delta }v_r/v`$ grows only when $`\gamma >4/3`$ (but recall that the scaling relations apply only for $`\gamma <5/3`$).

### 3.4 Special Cases: Possibility of Faster Growth

Now consider a hypothetical situation in which the inner boundary of the flow is a “sticky sphere”: Once a fluid element enters the sphere, it gets stuck on the spot where it enters. In this case, a finite multipole moment, $`Q_c`$, will accumulate at the center (see eq. ). In the region where $`Q_c`$ dominates the gravitational perturbation, we have $`\delta \psi Q_c/r^{l+1}r^2\delta \rho `$. It is easy to show from equations (20)-(21) that the perturbations have the following scaling behavior (for $`\gamma <5/3`$):

$`\delta u`$ $``$ $`{\displaystyle \frac{2}{2l1}}{\displaystyle \frac{r}{v}}\delta \psi r^{1/2l},`$ (30)   
$`\delta \rho `$ $``$ $`{\displaystyle \frac{3}{2}}{\displaystyle \frac{\rho \delta u}{rv}}r^{l3/2},`$ (31)   
$`\delta \psi `$ $``$ $`\left({\displaystyle \frac{4\pi G}{2l+1}}\right){\displaystyle \frac{Q_c}{r^{l+1}}}r^{l1}.`$ (32)   
These scalings depend on $`l`$. Even for $`l=1`$, the growth is faster than the case discussed in §3.2 and §3.3. These scalings can be understood as follows: The central multipole moment exerts a torque $`r^{l1}`$ and a radial force $`r^{l2}`$ on a fluid element. The angular momentum grows as $`t_r\delta \psi r^{l+1/2}`$, and the radial velocity perturbation grows as $`\delta v_rr^{l1/2}`$. Thus we have $`\delta v_r/v\delta v_{}/vr^l`$. The Lagrangian displacement is $`\xi _r\xi _{}t_r\delta vr^{l+1}`$, and the resulting density perturbation is $`\mathrm{\Delta }\rho /\rho \xi /rr^l`$.

As we have argued in §3.2, it is unlikely the scaling relations derived in this section will apply in general situations since it is hard to imagine that the central core can retain its multipole moments for a time much longer than the dynamical time of the flow at the inner boundary. There are two possible exceptions:

(i) For $`l=1`$ modes: The core has an extra degree of freedom, i.e., it can have linear motion in response to flow perturbations. This gives rise to a net central dipole moment and the possibility of faster growth of dipolar perturbations. Let $`m_c`$ be the core mass and $`Z_c`$ be the position (along the $`z`$-axis) of its center-of-mass relative to the the origin of coordinates. The displaced core induces a potential perturbation $`\delta \psi (Gm_cZ_c/r^2)\mathrm{cos}\theta `$. From equations (30)-(32), the flow perturbations are given by

$$\frac{\delta \rho }{\rho }\frac{3Z_c}{2r}\mathrm{cos}\theta ,\frac{\delta v_r}{v}\frac{Z_c}{2r}\mathrm{cos}\theta ,\frac{\delta v_{}}{v}\frac{Z_c}{r}\mathrm{sin}\theta ,$$
(33)   
where we have used $`v(2Gm_c/r)^{1/2}`$. However, the asymptotic perturbations given by equation (33) simply represent a spherical flow centered at $`Z_c`$. Thus it is not surprising that the sum of the gravitational force exerted on the core plus the rate at which momentum flows across a surface surrounding it sum to zero.

(ii) For $`l=2`$ modes: Suppose the central core consists of a rotating star, or a star with a massive circumstellar disk. This could result from an early phase of accretion/collapse with significant angular momentum. Subsequent spherical accretion, with no net angular momentum, will be affected by the central quadrupole.

## 4 Numerical Calculations of Linear Perturbation Growth

The asymptotic scaling relations derived in §3 apply only in the supersonic regime. To determine the behavior of perturbations under general conditions, we numerically follow the collapse of a self-gravitating cloud and evolve the nonspherical perturbation carried by each fluid element. Because of the large disparity in the timescales involved in the central region and the outer region, it is essential for the code to have a wide dynamical range. Previous multidimensional simulations of Bondi accretion (e.g., Ruffert 1994) did not achieve high enough resolution in the central region to reveal the growth of perturbation. We restrict our calculations to the linear regime. Thus perturbations associated with different $`Y_{lm}`$ evolve independently, and our calculations involve only one spatial dimension.

We have constructed a one-dimensional Lagrangian finite-difference code. The unperturbed flow variables $`(r,v,\rho )`$ are followed with a standard scheme (Bowers & Wilson 1991), and the Eulerian perturbations $`(\delta \rho ,\delta u,\delta \psi )`$ are evolved using equations (12), (14) and (15). The flow is covered by a uniform mass grid. The quantities $`r,v,\delta u,\delta \psi `$ are zone-edge-centered, while $`\rho ,\delta \rho `$ are zone-centered. A staggered leapfrog integration scheme is adopted to ensure second-order accuracy in time.

### 4.1 Pressureless Collapse

To calibrate our code and check the asymptotic scaling relations of §3, we study the collapse of a centrally concentrated dust cloud. Figure 1 shows an example of such a collapse calculation. The cloud, of total mass $`m=1`$ and radius $`r=1`$, is initially at rest with density profile $`\rho r^1`$; i.e. the radius of a mass shell with enclosed mass $`m`$ is $`r_m=m^{1/2}`$. We initialize an $`l=2`$ perturbation with $`\delta \rho /\rho =1`$ in arbitrary units and $`\delta u=0`$. The Poisson equation is solved to give the potential perturbation, $`\delta \psi `$. We impose an inner boundary at $`r_c=10^3`$: Once a mass shell enters this boundary, it is removed from the simulation domain, and its perturbations are immediately smeared out; i.e., the central quadrupole moment, $`Q_c`$, is maintained at zero. The profiles $`\rho r^{3/2}`$ and $`vr^{1/2}`$ for the unperturbed flow are established near the center as the collapse proceeds. Figure 1 depicts the evolution of the perturbation carried by three different mass shells. We see that as the mass shells collapse to small radii, the analytic asymptotic scalings $`\delta u`$ constant and $`\delta \rho /\rho r^{1/2}`$ are achieved. Calculations with other initial conditions confirm that these scalings are generic features of perturbation growth in the absence of central multipole moments.

We have also studied the case where the central multipole moment, $`Q_c`$, is nonzero, and confirmed the steeper scalings derived in §3.4.

### 4.2 Collapse with Finite Pressure

Next we study the collapse of clouds having finite pressure. For definiteness, we choose the initial cloud to be a $`\gamma =4/3`$ spherical polytrope in hydrostatic equilibrium. The collapse is initiated by reducing $`\gamma `$ to $`1.3`$, and by reducing $`K`$ by $`10\%`$. Either one of these reductions alone is adequate to induce the collapse. This mimics Type II supernova collapse, where a white dwarf core of a massive star collapses to a neutron star. However, to focus on the growth of perturbations during the collapse, we do not include a shock in our calculation. We define an inner boundary at $`r_c=0.005`$; the original cloud radius is $`r=1`$ and its mass $`m=1`$. When the central density becomes greater than $`10^6`$, corresponding to a few times nuclear density if the initial cloud has mass and radius typical of a Chandrasekhar mass white dwarf, we cut out the flow inside $`r_c`$ from the computational domain. This enables us to follow the collapse and accretion of the rest of the cloud.

Figure 2 shows the unperturbed density and velocity of several different mass shells as functions of their Lagrangian radii. The inner region collapses homologously; the Mach numbers of individual shells remain below unity outside $`r_c`$. The outer region of the flow goes through a transonic point, and eventually attains the free-fall asymptotics, with $`vr^{1/2}`$ and $`\rho r^{3/2}`$. Inspecting the flow profiles at different time slices (not shown), we confirm that our numerical results agree with the self-similar solution (Yahil 1983) in the regime where it applies. 

Figure 3 and Figure 4 show two examples of the evolution of $`l=2`$ perturbations during the collapse depicted in Fig. 2. In Fig. 3, the initial perturbation is chosen to be $`\delta \rho /\rho =1`$ and $`\delta u=0`$, while in Fig. 4, the initial perturbation corresponds to the eigenfunctions of the first g-mode of a $`\gamma =4/3`$ polytrope, with adiabatic index $`\gamma _1=5/3`$. This value of $`\gamma _1`$ is used only for setting up the initial perturbations; after the collapse starts, the adiabatic index is set to $`\gamma =1.3`$. We see that the perturbations carried by the inner mass shells ($`m=0.2,0.6`$), which never become supersonic, vary in an oscillatory manner with no increase in amplitude. This is consistent with the result of Goldreich & Weber (1980) that the homologous inner core of a collapsing $`\gamma =4/3`$ polytrope is stable against non-radial perturbations. However, the outer region of the cloud ($`m>0.85`$) attains a high Mach number and eventually approaches free-fall. We see from Figs. 3 and 4 that the density perturbation grows in this outer region, and that the asymptotic scaling relations derived in §3.2 are recovered.

## 5 Discussion

Nonspherical perturbations are amplified during the supersonic collapse or accretion of a centrally concentrated gas cloud. We have derived asymptotic scaling relations for their growth. These general results have implications for several different astrophysical problems which we now discuss.

In spherical accretion onto a black hole, we expect the radiative efficiency to increase as a result of nonradial perturbations in the flow. In the asymptotic regime, the tangential velocity perturbation scales as $`\delta v_{}r^1`$, and the corresponding Mach number scales as $`\delta v_{}/c_sr^{(3\gamma 7)/4}`$. For $`\gamma <5/3`$, the Mach number grows faster than $`r^{1/2}`$. Similarly, the Mach number associated with the radial velocity perturbation (see eq.  for $`\gamma >2/3`$) scales as $`\delta v_r/c_sr^{3(\gamma 1)/4}`$. One might expect the formation of shocks which lead to thermalization of the flow and high radiative efficiency (see Chang & Ostriker 1985 for previous discussion on the formation of shocks in spherical flows). This point has also been noted recently by Kovalenko & Eremin (1998), who derived similar scaling relations for Bondi accretion.

In the context of core-collapse supernovae, our results indicate that perturbations in the homologous inner core do not grow, but those in the outer core, involving $`15\%`$ of the core mass, and in the envelope are amplified. Since $`\delta \rho /\rho `$ scales as $`r^{1/2}`$, we expect that the amplification factor is at most $`10`$ for $`r`$ decreasing from $`1500`$ km to $`15`$ km. It is possible that dipole perturbations obey the $`r^1`$ scaling (see §3.4), in which case the amplification factor could be larger. Interestingly, if overstable g-modes driven by shell nuclear burning are responsible for the seed of presupernova perturbations (see GLS), it is exactly at the outer core where the perturbation amplitude is expected to be the largest. The asymmetric density perturbation may lead to asymmetric shock propagation and breakout, which then give rise to asymmetry in the explosion and a kick to the neutron star (e.g., Burrows & Hayes 1996).

Finally, we note that our analysis has neglected rotation (i.e., the flow does not have a net angular momentum). In the perturbative regime, rotation (around the $`z`$-axis) is represented by the last term of equation (2) with $`m=0`$, i.e., $`\delta 𝐯_{\mathrm{rot}}=\delta v_{\mathrm{rot}}(r,t)(Y_{l0}/\theta )\widehat{\varphi }`$. As shown in §2, this rotational perturbation is decoupled from the density perturbation. Therefore, provided the rotational velocity is small in comparison to the radial velocity, i.e., $`|\delta v_T||v|`$, we expect our scaling relations for the growth of perturbations to be valid.

A major portion of this research was done between 1995 and 1997, when D.L. was a postdoc in theoretical astrophysics at Caltech; support from a Richard C. Tolman fellowship is gratefully acknowledged. This work is supported in part by NASA grant NAG 5-8356, by a research fellowship from the Alfred P. Sloan foundation, and by NSF grant 94-14232.    

## Appendix A Perturbation of Bondi Accretion for $`\gamma =5/3`$

The Bondi solution for $`\gamma =5/3`$ is special because the flow remains subsonic for $`r>0`$. We can think of the sonic point as being located at $`r=0`$ (in Newtonian theory). For $`rGM/c_{\mathrm{}}^2`$, where $`M`$ is the central mass and $`c_{\mathrm{}}`$ is the sound speed at infinity, the flow velocity and density are given by

$`|v|`$ $``$ $`c_s\left({\displaystyle \frac{GM}{2r}}\right)^{1/2},`$ (A1)   
$`\rho `$ $``$ $`\rho _{\mathrm{}}\left({\displaystyle \frac{GM}{2c_{\mathrm{}}^2}}\right)^{3/2}r^{3/2}.`$ (A2)   
Although these have the same scalings as equations (17) and (18), the asymptotic behavior of the perturbations is quite different. Indeed, for $`\gamma =5/3`$, eqs. (20) and (21) yield

$$\delta \rho r^{l1/2},\delta ur^{l+3/2},\delta \psi r^l.$$
(A3)   
Thus all perturbations decrease as $`r`$ decreases in keeping with the subsonic nature of the unperturbed flow. If a central multipole moment is present (see §3.4), the following asymptotics become dominant:

$$\delta \rho r^{l3/2},\delta ur^{l+1/2},\delta \psi r^{l1}.$$
(A4)   
## Appendix B Perturbations in Collapsing Spherical Shells

Consider a spherical fluid shell of mass $`M_s`$ falling from infinity onto a central point mass $`M_c`$. The shell radius, $`R_0(t)`$, evolves in time according to

$$\frac{d^2R_0}{dt^2}=\frac{GM}{R_0^2},\mathrm{with}M=M_c+\frac{M_s}{2},$$
(B1)   
which gives

$$R_0(t)=\left(\frac{9GM}{2}\right)^{1/3}(t)^{2/3},$$
(B2)   
where we have set $`t=0`$ at the point of complete collapse. The surface density and radial velocity are given by

$$\mathrm{\Sigma }_0(t)=\frac{M_s}{4\pi R_0^2}(t)^{4/3},V_0(t)=\frac{dR_0}{dt}(t)^{1/3}.$$
(B3)   
### B.1 Perturbation Equations

The dynamical variables for a perturbed shell are its surface density $`\mathrm{\Sigma }(\theta ,t)`$, radius $`R(\theta ,t)`$, radial velocity $`V_r(\theta ,t)=\dot{R}`$, where dot indicates $`/t`$, and tangential velocity $`𝐕_{}(\theta ,t)`$. We assume spherical coordinates and axisymmetry so that there is no $`\varphi `$-dependence. The core mass, $`M_c`$, is free to move. We use $`Z_c(t)`$ to denote its displacement from the coordinate origin. Note that $`\mathrm{\Sigma },V_r,𝐕_{}`$ are rigorously defined from the three-dimensional fluid variables via

$$\mathrm{\Sigma }\frac{1}{R^2}_R^{R+}\rho r^2𝑑r,V_r\frac{1}{\mathrm{\Sigma }R^2}_R^{R+}\rho r^2v_r𝑑r,𝐕_{}\frac{1}{\mathrm{\Sigma }R^2}_R^{R+}\rho r^2𝐯_{}𝑑r,$$
(B4)   
where the integration runs through the thickness of the shell. Using these definitions and the standard hydrodynamical equations, we derive the continuity and Euler equations for the shell:

$`{\displaystyle \frac{\mathrm{\Sigma }}{t}}={\displaystyle \frac{2\dot{R}}{R}}\mathrm{\Sigma }\mathrm{\Sigma }_{}𝐕_{},`$ (B5)   
$`{\displaystyle \frac{V_r}{t}}=\ddot{R}={\displaystyle \frac{1}{2}}\left[\left({\displaystyle \frac{\mathrm{\Phi }}{r}}\right)_{R+}+\left({\displaystyle \frac{\mathrm{\Phi }}{r}}\right)_R\right],`$ (B6)   
$`{\displaystyle \frac{𝐕_{}}{t}}={\displaystyle \frac{\dot{R}𝐕_{}}{R}}{\displaystyle \frac{1}{2}}\left[\left(_{}\mathrm{\Phi }\right)_{R+}+\left(_{}\mathrm{\Phi }\right)_R\right],`$ (B7)   
where we have assumed that the nonspherical perturbation is small. Note that here $`_{}(1/R)\widehat{}_{}=(1/R)\left[\widehat{\theta }(/\theta )+(\widehat{\varphi }/\mathrm{sin}\theta )(/\varphi )\right]=(1/R)\widehat{\theta }(/\theta )`$. The gravitational potential $`\mathrm{\Phi }=\mathrm{\Phi }_c+\mathrm{\Phi }_s`$ includes contributions from the core, $`\mathrm{\Phi }_c`$, as given by

$$\mathrm{\Phi }_c(𝐫,t)=\frac{GM_c}{|𝐫Z_c\widehat{z}|}=GM_c\underset{l}{}\frac{Z_c^l}{r^{l+1}}\left(\frac{4\pi }{2l+1}\right)^{1/2}Y_{l0}(\theta ,\varphi ),$$
(B8)   
for $`R>Z_c`$. The potential produced by the shell satisfies

$$^2\mathrm{\Phi }_s(𝐫,t)=4\pi G\mathrm{\Sigma }(\theta ,t)\delta [rR(\theta ,t)].$$
(B9)   
Finally we need the equation of motion for the core mass:

$$\ddot{Z}_c=\left(\frac{\mathrm{\Phi }_s}{r}\right)_{r=Z_c,\theta =0}.$$
(B10)   
Consider linear perturbation modes associated with spherical harmonics $`Y_{l0}`$:

$`R(\theta ,t)=R_0(t)\left[1+a_l(t)Y_{l0}\right],`$ (B11)   
$`𝐕_{}(\theta ,t)=\dot{R}_0(t)b_l(t)\widehat{}_{}Y_{l0},`$ (B12)   
$`\mathrm{\Sigma }(\theta ,t)=\mathrm{\Sigma }_0(t)\left[1+c_l(t)Y_{l0}\right].`$ (B13)   
Solving equation (B9), the shell potential to linear order in $`a_l,c_l`$ is

$$\mathrm{\Phi }_s(𝐫,t)=\frac{GM_s}{r}\frac{1}{2l+1}\left(\frac{GM_s}{R_0}\right)\left(\frac{R_0}{r}\right)^{l+1}\left[c_l+(l+2)a_l\right]Y_{l0},$$
(B14)   
for $`r>R(\theta ,t)`$ (outside the shell), and

$$\mathrm{\Phi }_s(𝐫,t)=\frac{GM_s}{R_0}\frac{1}{2l+1}\left(\frac{GM_s}{R_0}\right)\left(\frac{r}{R_0}\right)^l\left[c_l(l1)a_l\right]Y_{l0},$$
(B15)   
for $`r<R(\theta ,t)`$ (inside the shell). Using the unperturbed solution (eqs. \[B2\] and \[B3\]), the perturbation equations (B5)-(B7) for the shell reduce to

$`2t\dot{a}_l+t\dot{c}_l{\displaystyle \frac{2}{3}}l(l+1)b_l=0,`$ (B16)   
$`9t^2\ddot{a}_l+12t\dot{a}_l2a_l={\displaystyle \frac{4M_c}{M}}\delta _{l1}z_c+{\displaystyle \frac{4M_c}{M}}a_l{\displaystyle \frac{2M_s}{(2l+1)M}}\left[l(l1)a_l+{\displaystyle \frac{1}{2}}c_l\right],`$ (B17)   
$`3t\dot{b}_l+b_l={\displaystyle \frac{M_c}{M}}\delta _{l1}z_c+{\displaystyle \frac{M_s}{(2l+1)M}}\left(c_l+{\displaystyle \frac{3}{2}}a_l\right),`$ (B18)   
where $`\delta _{l1}`$ is Kronecker delta, and

$$z_c\left(\frac{4\pi }{3}\right)^{1/2}\left(\frac{Z_c}{R_0}\right).$$
(B19)   
The equation of motion for the core mass, equation (B10), becomes

$$9t^2\ddot{z}_c+12t\dot{z}_c2z_c=\frac{2M_s}{3M}c_1.$$
(B20)   
Setting

$$a_l,b_l,c_l,z_c(t)^s,$$
(B21)   
equations (B16)-(B20) reduce to a set of algebraic equations from which the eigenvalue, $`s`$, and the corresponding eigenmode can be determined. We discuss these eigenmodes below.

### B.2 $`l=1`$ Modes

There are six roots for $`s`$. Two of these, $`s=1/3,2/3`$, are trivial modes which do not involve any surface density perturbation ($`c_1=0`$), and for which the core experiences no acceleration ($`\ddot{Z}_c=0`$). These correspond to the collapse of a uniform shell onto a displaced core; for $`\mathrm{\Delta }R=RR_0(t)`$, $`\mathrm{\Delta }R/R_0(t)^{1/3}`$, and for $`\mathrm{\Delta }R=`$ constant, $`\mathrm{\Delta }R/R_0(t)^{2/3}`$.

Two of the remaining four roots correspond to stable modes ($`s>0`$). The two unstable modes are shown in Fig. 5. Both lead to the growth of the separation between the center of mass of the shell and the position of the core: (i) “Bending” Mode: the surface density perturbation grows because one side of the shell collapses faster than the other, and the geometric center of the shell moves in opposition to the motion of the central mass; when $`M_s0`$, the mode has $`s=1`$ and $`(a_1,b_1,c_1)=(1,0,2)`$; (ii) “Jeans” Mode: the surface density perturbation grows due to the internal tangential flow in the shell, while the shell’s geometric center suffers little displacement with respect to the position of the central mass; when $`M_s0`$, the mode has $`s=1/3`$ and $`(a_1,b_1,c_1)=(0,1,4)`$. The “bending” mode grows more rapidly, with $`s`$ ranging from $`1`$ for $`M_s0`$ to $`s=(\sqrt{17}+1)/6=0.854`$ for $`M_c0`$.

### B.3 $`l=2`$ Modes

For $`l=2`$, the core mass experiences no acceleration. There are two types of unstable modes as shown in Fig. 6: (i) “Bending” Mode: the north pole and south pole of the shell collapse faster and have higher density than the equator, leading to a quadrupolar density perturbation; (ii) “Jeans” Mode: the density perturbation is mainly due to tangential fluid motion within the shell. The “bending” mode is the more rapidly growing mode, with $`s`$ in the range between $`1`$ and $`0.85`$;

### B.4 Large-$`l`$ Limit

Our results for general $`l`$ will not be presented here. But it is of interest to consider the $`l1`$ limit.

(i) “Bending” Mode: The dispersion relation of bending waves on a pressureless, nonrotating surface is $`\omega ^2=2\pi G\mathrm{\Sigma }_0|k|`$, with $`|k|l/R_0`$. This gives

$$\omega =\pm \frac{s_l}{t},\mathrm{with}s_l\frac{1}{3}\left(\frac{lM_s}{M}\right)^{1/2}.$$
(B22)   
The wave evolves as $`\mathrm{exp}(i\omega 𝑑t)=(t)^{\pm is_l}`$. To obtain the amplitude evolution, we use the conservation of wave action (energy per unit mass divided by frequency) $`\omega |\mathrm{\Delta }R|^2`$, which gives $`|\mathrm{\Delta }R|(t)^{1/2}`$, and $`|\mathrm{\Delta }R|/R_0(t)^{1/6}`$. Thus

$$\frac{\mathrm{\Delta }R}{R_0}(t)^s,s=\frac{1}{6}\pm is_l.$$
(B23)   
This agrees with our numerical results.

(ii) “Jeans” Mode: The dispersion relation of density wave is $`\omega ^2=2\pi G\mathrm{\Sigma }_0|k|`$, which gives $`\omega =\pm is_l/t`$. Similar to (i), we find

$$\frac{\mathrm{\Delta }R}{R_0}(t)^{1/6}\mathrm{exp}(i\omega 𝑑t)(t)^s,s=\frac{1}{6}\pm s_l.$$
(B24)   
This also agrees with our numerical results.

### B.5 Effect of Internal Pressure

We can include the effect of fluid pressure by adding a term $`(1/\mathrm{\Sigma })_{}(\mathrm{\Sigma }C_s^2)`$ to the tangential Euler equation (B7); the radial equation (B6) is not affected by pressure. For definiteness, we parametrize the shell-averaged sound speed, $`C_s`$, by

$$C_s=\beta \left(\frac{GM}{R_0}\right)^{1/2},$$
(B25)   
where $`\beta `$ is a constant. This amounts to adding a term $`\beta ^2c_l`$ to the right-hand-side of equation (B18).

Figures 5 and 6 show the effect of pressure on the eigenvalues of the $`l=1`$ and $`l=2`$ modes, respectively. We see that pressure always tends to stablize the modes. However, for $`M_s/M_c`$ less than a few, the “bending” mode is only slightly affected.

In the large-$`l`$ limit, the “bending” mode is unaffected by the pressure, thus equation (B23) still applies. For the “Jeans” mode, the dispersion relation is $`\omega ^2=k^2C_s^22\pi G\mathrm{\Sigma }_0|k|`$, which gives

$$\omega ^2=\left(2\beta ^2l^2\frac{lM_s}{M}\right)\frac{1}{9t^2}.$$
(B26)   
For $`lM_s/(2\beta ^2M)`$, we have $`\omega \pm (\sqrt{2}/3)\beta l/t`$. Using similar precedure as in §B.4, we obtain

$$\frac{\mathrm{\Delta }R}{R_0}(t)^s,\mathrm{with}s\frac{1}{6}\pm i\frac{\sqrt{2}}{3}\beta l.$$
(B27)