aa r X i v : . [ a s t r o - ph ] J un Critical Accretion Rate for Triggered Star Formation
Tomoyuki Hanawa and Akihito Soeda [email protected] ABSTRACT
We have reexamined the similarity solution for a self-gravitating isothermalgas sphere and examined implication to star formation in a turbulent cloud.When parameters are adequately chosen, the similarity solution expresses anaccreting isothermal gas sphere bounded by a spherical shock wave. The massand radius of the sphere increases in proportion to the time, while the centraldensity decreases in proportion to the inverse square of time. The similaritysolution is specified by the accretion rate and the infall velocity. The accretionrate has an upper limit for a given infall velocity. When the accretion rate isbelow the upper limit, there exist a pair of similarity solutions for a given set ofthe accretion rate and infall velocity. One of them is confirmed to be unstableagainst a spherical perturbation. This means that the gas sphere collapses toinitiate star formation only when the accretion rate is larger than the upper limit.We have also examined stability of the similarity solution against non-sphericalperturbation. Non-spherical perturbations are found to be damped.
Subject headings: accretion — hydrodynamics — shock waves — stars: formation
1. INTRODUCTION
Similarity solutions have contributed very much to our understanding of star formationprocess. The classical similarity solution by Larson (1969) and Penston (1969) elucidatedthe runaway nature of gravitational collapse. The density increases in proportion to theinverse square of the time during the runaway collapse phase. We learned from the similaritysolutions of Shu (1977) and Hunter (1977) that the accretion rate of a protostar is of theorder of c s /G where G and c s denote the gravitational constant and the isothermal soundspeed of gas. Center for Frontier Science, Chiba University, Inage-ku, Chiba, 263-8522, Japan II region.Tsai & Hsu (1995) found two classes of similarity solution; the first class describesaccretion onto protostar while the second one does failure of star formation. The centraldensity decreases in proportion to the inverse square of the time, ρ c ∝ t − , in the second classsolution. Although this solution has not gained much attention thus far, it provides an insighton dynamical compression of a molecular cloud core. If a dense clump of gas is compressedby an external force, the temporal increase in the density may trigger gravitational collapseand star formation. One can surmise existence of threshold of gravitational collapse. If thedynamical compression is either weak or short, the clump will bounce back to expansion.A shock wave will be formed when accreting gas is stoped by the expansion (Adams & Shu2007, see, e.g.). The similarity solution of Tsai & Hsu (1995) demonstrated that a sphericalcloud can expand even when it is steadily compressed by a shock wave. On the other hand,the shock compressed gas sphere will collapse owing to its self gravity if the shock is strongand lasts for a long enough period.In this paper we reexamine the similarity solution of Tsai & Hsu (1995) while keepingits negative implication in mind. We find the condition for existence of similarity solutiondescribing expansion of a gas sphere. Conversely it will tell us condition for a shock com-pressed gas sphere to collapse by its self gravity. We also study stability of the similaritysolution. The similarity solution denies collapse due to the self gravity only when it is stable.We review the similarity solution in § § § § § §
4. 3 –
2. Model and Methods of Computation2.1. Similarity Solution
We consider an isothermal gas of which distribution is spherically symmetric. Then thehydrodynamical equations are expressed as ∂ρ∂t + 1 r ∂∂r ( r ρv ) = 0 , (1) ∂v∂t + v ∂v∂r + 1 ρ ∂P∂r + GM r r = 0 , (2)and ∂M r ∂r = 4 πr ρ (3)where P = c ρ . (4)Here, the symbols, ρ , v , P , M r , G , and c s denote the density, velocity, pressure, massinside the radius r , gravitational constant, and the isothermal sound speed, respectively. Asoriginally shown by Larson (1969) and Penston (1969), the hydrodynamical equations havesimilarity solutions, ρ ( r, t ) = ̺ ( ξ )4 πGt , (5) v ( r, t ) = u ( ξ ) c s , (6) M r ( r, t ) = c tG µ ( ξ ) (7)where ξ = rc s t . (8)We restrict ourselves in the case of t > ∂̺∂ξ = − [ µ − ξ − u ) ] ̺ [( u − ξ ) − ξ (9) ∂u∂ξ = ( µ −
2) ( ξ − u )[( u − ξ ) − ξ (10)and µ = ξ ( ξ − u ) ̺ . (11) 4 –We assume that the density is finite at the center since we are interested in applicationto a starless core, i.e., case of no star formation. Then the density and velocity should beexpressed as ̺ = ̺ − ρ (cid:18) ̺ − (cid:19) ξ + ρ (cid:18) ̺ − (cid:19) (cid:18) ̺ − (cid:19) + O ( ξ ) (12) u = 2 ξ − ̺ (cid:18) ̺ − (cid:19) ξ + O ( ξ ) (13)near the center. Following Tsai & Hsu (1995) we assume that the flow has a shock wave at ξ = ξ sh . Then the Rankine-Hugonio relation gives us the condition, ̺ + ̺ − = u + − ξ sh u − − ξ sh , (14)and ( u + − ξ sh ) ( u + − ξ sh ) = 1 . (15)where ̺ + and ̺ − denote the the densities of the pre- and post-shocked gases, respectively,and u + and u − do the velocities of them, respectively.In the region far from the origin, a solution of Equations (9) and (10) approaches theasymptotic solution, ̺ = ˙ Mv inf ξ + O ( ξ − ) , (16) u = − v inf + O ( ξ − ) . (17)The symbol, v inf , denotes the infall velocity at the infinity while ˙ M denotes the accretionrate. The similarity solution can be specified either by ( ̺ c , ξ sh ) or by ( v inf , ˙ M ).Note that Larson-Penston solution has the same asymptotic form. Only when v inf and˙ M vanish, the solution Equations (9) and (10) have a different asymptotic form, “plussolutions”of Shu (1977).We integrate equations (9) and (10) by the 4th order Runge-Kutta method from ξ = 0for a given set of ̺ c and ξ sh to obtain a similarity solution. The infall velocity and accretionrate are obtained numerically as a function of ̺ c and ξ sh . We also obtain the mass enclosedin the shock front, M c = Z ξ sh ̺ ξ dξ . (18) 5 – We have performed a normal mode analysis to examine the stability of the similaritysolution. In the analysis, the density is assumed to be expressed as ρ ( r, t ) = ̺ ( ξ ) + t σ δ̺ ( ξ ) Y mℓ ( θ, ϕ )4 πGt , (19)where Y mℓ ( θ, ϕ ) denotes the spherical harmonic function. This particular form is chosen be-cause the similarity solution has no specific timescale. An eigenmode has a growth timescaleand an unstable perturbation grows exponentially, when the unperturbed state is stationaryand has a specific timescale. When the physical quantities vary according to a power lawin time in an unperturbed state, the eigenmode grows (or decay) in proportion to a powerof time, | t | σ . See Hanawa & Matsumoto (1999) for the justification of Equation (19). Theyanalyzed linear stability of the Larson-Penston solution against a non-spherical perturba-tion, using the coordinates, ( ξ, θ, ϕ, ln t ) instead of the ordinary spherical coordinates,( r, θ, ϕ, t ). It is shown that the similarity solution can be expressed as a steady state andan unstable mode grows in propotion to exp ( σ ln t ) = t σ in the coodinates.The power index, σ , is obtained as an eigenvalue of the perturbation equations as shownin Appendix. When the real part of σ is positive, the mode grows in time and the similaritysolution is unstable. We call the real part of σ , the growth index in the following accordingto the suggestion by F. H. Shu, the referee of this paper. When σ is complex, the modegrows (or decays) while oscillating. The imaginary part of σ is related with the frequency ofoscillation, although the physical oscillation period increases with time.
3. Results3.1. Similarity Solutions
First we obtained a series of similarity solutions for a given ̺ c by increasing ξ sh . Theinfall velocity, v inf , decreases monotonically with increase in ξ sh . It reaches v inf = 0 at some ξ sh and the similarity solution terminates. By compiling the similarity solutions, we obtained˙ M as a function of v inf and ξ sh . The curves denote ˙ M as a function of ξ sh for given v inf inFigure 1.The accretion rate, ˙ M , is maximum at a certain ξ sh for a given v inf . It increases withincrease in ξ sh at a small ξ sh while it decreases with increase in ξ sh at a large ξ sh . The formeris denoted by thin curves while the latter is by thick ones. This means that there exists twosimilarity solutions for a given set of v inf and ˙ M . 6 –Figure 2 shows two similarity solutions having v inf = 1 . M = 1 . ξ sh = 0 .
986 while the dashed curves do that of ξ sh = 0 . ξ ≥ . § M for a given v inf .The upper limit is highest ˙ M max = 1 .
312 at v inf = 1 .
64 . Similarity solutions do not existfor ˙
M > . § Figure 4 shows the growth index of the spherical perturbation, σ r , as a function of ξ sh for a series of similarity solutions having a given v inf . The solid curves denote the growthindex of modes having real index. The dashed lines denote the real part of complex indecies.One of the index is positive ( σ r >
0) and the similarity solution is unstable only when theshock radius ( ξ sh ) is smaller than a critical value. The condition of neutral stability coincideswith that of maximum accretion rate for a given v inf as expected. When ξ sh is a little largerthan the critical, the solution is stable and the most slowly damping mode has a real index.When ξ sh is large enough, the solution is stable and the most slowly damping mode has acomplex index.Our survey is limited to modes having low indecies, i.e., in the range | σ i | ≤ .
0. Itis however unlikely that we have missed unstable spherical perturbations. A spherical per-turbation induces only sound waves and they are confined within the gas sphere since it isbounded by the shock wave. A high frequency sound wave has a shorter wavelength and isunlikely to be Jeans unstable. ℓ = 2 Mode We have studied ℓ = 2 mode as a typical non-spherical perturbation. This is in partbecause the dipole ( ℓ = 1) mode is unlikely to be excited. If the dipole mode grows, theinner and outer parts of the gas sphere should move in the opposite direction each other tokeep the center of gravity. 7 –Figure 6 denotes the eigenfrequencies of the ℓ = 2 modes for similarity solutions of v inf = 1.0. The abscissa denotes ξ sh while the ordinate denotes the index, σ . All the modesare damping ( σ i < − σ i ≃ ± . − . ≤ σ r ≤ − . v inf = 4.0. Again, all the modes are damping.The damping index is larger than unity ( σ r < − v inf = 1.0.One might ask the reason why the similarity solution is stable against a bar ( ℓ = 2)mode. We think that the bar mode is damped by expansion of the gas sphere. Since theradius of the shock front increases with the time, the asphericity of the shock compressed gassphere decreases unless the displacement grows faster than the radius. When a gas sphere iscollapse, the bar mode can be excited as shown by Lin, Mestel, & Shu (1965) for the pressureless gas and Hanawa & Matsumoto (1999) for an isothermal gas.We have not yet studied the modes of ℓ ≥
3. However they are also unlikely to beunstable, since the self-gravity does not excite a short wave perturbation.
4. Discussion
We obtained the critical accretion rate above which there exits no similarity solution.The critical rate can be interpreted as the minimum accretion rate for a high density clumpto initiate self-gravitational collapse. The critical accretion rate can be rewritten as dMdt (cid:12)(cid:12)(cid:12)(cid:12) cr = 3 . c s Gv , (20)for v & c s in the dimensional form.Equation (20) gives us an estimate for a converging flow to initiate gravitational collapse.We shall consider a spherical region of which surface is surrounded by a converging flow. Theradius, inflow velocity, and density are assumed to be r , v and ρ , respectively. Then thegravitational collapse will be initiated when the mass accretion rate exceeds the critical,4 πr ρ v > . c s G v . (21) 8 –Equation (21) can be rewritten as 2 rλ J > . c s v , (22)where λ J = 2 π c s p πGρ . (23)Since λ J denotes the Jeans length, Equation (22) means that the effective Jeans lengthreduces in proportion to the inverse of the Mach number.The Jeans mass is proportional to the cube of the Jeans length for a given density. Thusthe effective Jeans mass should reduce to M J , eff = M J (cid:18) λ J , eff λ J (cid:19) (24)= M J (cid:18) v . c s (cid:19) − . (25)This implies that compression of sub Jeans mass clump may result in gravitational collapsein the region of flow convergence. Note that the effective Jeans mass is several order ofmagnitude smaller than the classical one when v & c s .The compression should continue for a certain timescale for a dynamically compressedclump to collapse by its self gravity. If we evaluate the minimum timescale to be the effectiveJeans length divided by the flow velocity, it is shorter than the free-fall timescale by a factorof the Mach number squared, τ comp ≃ λ J , eff v ≃ τ ff (cid:18) vc s (cid:19) − . (26)The timescale can be translated into the wavelength of perturbation. A compressedclump can collapse by the self gravity if the wavelength of velocity perturbation is longer thanthe effective Jeans length. If turbulence contains velocity perturbations of long wavelengths,gravitational collapse due to dynamical compression will take place somewhere in the cloud.In such case we can expect a number of clumps of which masses are much smaller than theclassical Jeans mass.We thank T. Matsumoto and F. Nakamura for discussion and valuable comments onthe original manuscript. We also thank F. H. Shu for his valuable comments as the referee.We have added comparison of Tsai & Hsu (1995) solution with Bonnor-Ebert solution inthe revised manuscript according to the comments given to the original manuscript. This 9 –study is financially supported in part by the Grant-in-Aid for Scientific Research on PriorityArea (19015003) of The Ministry of Education, Culture, Sports, Science, and Technology(MEXT). A. Linear Stability Analysis
Since the density is given by Equation (19) in our linear stability analysis, the gravita-tional potential should be expressed asΦ ( r, t ) = c s [ ψ ( ξ ) + t σ δψ ( ξ ) Y mℓ ( θ, ϕ )] . (A1)The velocity, v = ( v r , v θ , v ϕ ), is assumed to be expressed as v r ( r, t ) = c s [ u ( ξ ) + t σ δu r Y mℓ ( θ, ϕ )] , (A2) v θ ( r, t ) = − c s t σ δu θ ( ξ ) ℓ + 1 ∂∂θ Y mℓ ( θ, ϕ ) , (A3) v ϕ ( r, t ) = − c s t σ δu θ ( ξ )( ℓ + 1) sin θ ∂∂ϕ Y mℓ ( θ, ϕ ) , (A4)in the spherical coordiantes. We assume that the perturbation is vanishingly small in theregion very far from the center. In other words, we restrict ourselves to search for aninstability due to the internal structure of the shock compressed gas sphere. Thus the flowis assumed to have no vorticity,( ℓ + 1) δu r + ddξ ( ξ δu θ ) = 0 . (A5)The assumption of no vorticity is based on the Kelvin’s circulation theorem. It says thecirculation, Γ ≡ I C u · d u (A6)is conserved for an isothermal (barotropic) gas (see, e.g., Chapter 6 of Shu 1992). Since theLagrangian loop C expands, the circulation velocity should decrease with the time to conserveΓ, although it can grow with the time if measured at a given ξ . See Hanawa & Matsumoto(2000) for more details on the mathematical treatment of perturbations having vorticity.After some manipulation, the perturbation equations are written as( σ + 1) δ̺ + 1 ξ ddξ (cid:0) ξ δf (cid:1) + ℓ ̺ δu θ ξ = 0 , (A7) 10 –( σ + 2) δf + 1 ξ (cid:8) ξ (cid:2) w δf + (1 − w ) (cid:3)(cid:9) − (cid:18) ̺ w + 2 ξ (cid:19) δ̺ + ̺ δγ + ℓ̺ wδu θ ξ = 0 , (A8) ddξ δψ = δγ , (A9) ddξ δγ = δ̺ − δγξ − ℓ ( ℓ + 1) ξ δψ (A10)where w = u − ξ , (A11) δu θ = ℓ + 1( σ + 1) ξ (cid:18) w δu r + δ̺̺ + δψ (cid:19) , (A12) δf = ̺ δu r + w δ̺ , (A13) δγ = ∂∂ξ δψ . (A14)The perturbation equations do not contain the azimuthal wavenumber, m , and accordinglythus the growth rate does not depend on m . This is because the unperturbed state (similaritysolution) is spherically symmetric.We must consider the jump condition at the shock front as well as the boundary condi-tions to solve the perturbation equation. The density perturbation is discontinuous betweenthe pre- and post-shocked flows and contains change due to shift of the shock front. The lat-ter introduces the Dirac’s delta function proportional to the shift and the density differencebetween the pre- and post flows. After some manipulation we obtain the jump condition,( w + − w − ) [ δf ] − ( ̺ + − ̺ − ) (cid:20) wδu r + δ̺̺ (cid:21) = 0 , (A15)( σ + 1) [ δγ ] + [ δf ] = 0 (A16)where the bracket denotes the difference between ξ = ξ sh ± ε . The symbols with thesubscripts, + and − , denote the values at ξ = ξ sh ± ε , respectively. The perturbation in thepotential, δψ , is continuous even at ξ = ξ sh . These jump conditions enable us to connect theperturbation in the post-shocked flow ( ξ < ξ sh ) and that in the pre-shocked flow ( ξ > ξ sh ).When the perturbation is spherically symmetric ( ℓ = 0), Equations (A7) and (A10) arelinearly dependent and we obtain δγ = − δfσ + 1 (A17) 11 –for ℓ = 0. Thus we need two independent boundary conditions for ℓ = 0 and three for ℓ = 0.The spherical perturbation should vanish in the pre-shocked since the unperturbed flowis supersonic. Otherwise the perturbation diverges when the mode is unstable. Thus theeffective outer boundary is set at ξ = ξ sh and the jump condition is applied there. The otherboundary condition is set so that the velocity perturbation vanishes at the origin, ξ = 0and is proportional to ξ near the origin. We have integrated the perturbation equationnumerically with the Runge-Kutta method from the origin and searched the eigenvalue, σ ,by try and error.Non-spherical perturbations should also be regular at the origin and vanishingly small atinfinity. We calculated asymptotic solutions around the origin and those around the infinityto examine the condition for perturbations to be regular according to Hanawa & Matsumoto(1999). Some asymptotic solutions are divergent since the perturbation equations are sin-gular at the origin and infinity. When ℓ = 0, the perturbation should be expressed by theasymptotic solution, δ̺ = B̺ ξ ℓ , (A18) δu r = − Aℓξ ℓ − − C ( ℓ + 2) ξ ℓ + 1 , (A19) δψ = (cid:20) A (cid:18) σ + 1 − ℓ (cid:19) − B (cid:21) ξ ℓ , (A20)where C = 14 ℓ + 6 (cid:20) ℓ (cid:18) ̺ − (cid:19) A + (cid:18) σ − ℓ − (cid:19) B (cid:21) . (A21)Here the symbols, A and B , denote arbitrary constants. Equations (A18) and (A20) meanthat the density and potential perturbations can be expanded by a series of polynomialsstarting from the ℓ -th power in the Cartesian coordinates. Otherwise the perturbationdiverges at the origin. The velocity perturbation is also expanded by another series ofpolynomials starting from the ( ℓ − δ̺ = 2 D ( ℓ + 1) ̺ ( σ + ℓ + 2)( σ + ℓ + 3) ξ ℓ + 3 , (A22) δu r = D ( ℓ + 1)( σ + ℓ + 2) ξ ℓ + 2 , (A23) δψ = Dξ ℓ + 1 , (A24) 12 –for ξ ≫
1. Here, the symbol D denotes an arbitrary constant. Equation (A24) meansthat the potential perturbation in the region ξ ≫ ξ = 10 − ) and a verylarge ξ ( ≃ σ , by examiningwhether the numerically solutions from both ends satisfy the jump condition at the shockfront. REFERENCES
Adams, F. C., & Shu, F. H. 2007, ApJ, 671, 497Bonnor, W. B. 1956, MNRAS, 116, 351Ebert. R. 1955, ZAp, 37, 217Hanawa, T., & Matsumoto, T. 1999, ApJ, 521, 703Hanawa, T., & Matsumoto, T. 2000, ApJ, 540, 962Hunter, C. 1977, ApJ, 218, 834Krasnopolsky, R., & K¨onigl, A. 2002, ApJ, 580, 987Larson, R. B. 1969, MNRAS, 145, 271Larson, R. B. 2003, Reports on Progress in Physics, 66, 1651Lin, C. C., Mestel, L., & Shu, F. H. 1965, ApJ, 143, 1431Narita, S., Hayashi, C., Miyama, S. M. 1984, Prog. Theor. Phys., 72, 1118Penston, M. V. 1969, MNRAS, 144, 425Saigo, K., & Hanawa, T. 1998, ApJ, 493, 342Shu, F. H. 1977, ApJ, 214, 488 13 –Shu, F. H. 1992, The Physics of Astrophysics Volume II Gas Dynamics, (Mill Valley, Uni-versity Science Books)Shu, F. H., Lizano, S., Galli, D., Cant´o, J. & Laughlin, G. 2002, ApJ, 580, 969Tsai, J. C., Hsu, J. J. L. ApJ, 448, 774
This preprint was prepared with the AAS L A TEX macros v5.2.
14 – ξ sh d M / d t v inf = 1.02.03.04.05.06.07.0 Fig. 1.— The accretion rate, ˙ M , is shown as a function of xi sh for a given v inf . 15 – ξ = r /( c s t )-2-10123 l og ρ , v log ρ v d M /dt = 1.204 v ∞ = 1.4 Fig. 2.— The density and velocity distributions are shown for two similarity solutions havingthe same v inf and dM/dt . 16 – v inf d M / d t c r i v inf Fig. 3.— The critical accretion rate is shown as a function of v inf . 17 – ξ sh -2-1012 σ r v inf =7.06.05.0 1.02.03.04.0 Fig. 4.— Each curve denotes the real part of eigenfrequency, σ r , of a spherical perturbationas a function of the shock radius, ξ sh for a given v inf . 18 – ξ sh σ i v inf =7.0 6.0 5.0 4.0 3.0 2.01.0 Fig. 5.— The same as Fig. 4 but for the imaginary part, σ i . 19 – ξ sh -2.5-2.0-1.5-1.0-0.50.00.5 σ σ i σ r v inf =1.0 l =2 Fig. 6.— The eigenfrequency, σ , is shown as a function of ξ sh for similarity solutions having v inf = 1 .
0. The solid curves denote the real part of σ for non-spherical perturbations having ℓ = 2, while dashed curves do that of imaginary part. Mode a has complex eigenfrequencieswhile modes b and c have real eigenfrequencies. 20 – ξ sh -2.5-2.0-1.5-1.0-0.50.00.5 σ σ i σ r v inf =4.0 l =2 Fig. 7.— The same as Fig, 6 but for v infinf