分子雲面狀結構之塌縮與穩定分析

全文

(1)國立臺灣師範大學理學院地球科學系碩士論文 Department of Earth Sciences College of Science. National Taiwan Normal University Master’s Thesis. 分子雲面狀結構之塌縮與穩定分析 Collapse and Stability Analysis of Sheet-like Structures in the Molecular Cloud. 沈孟嫺 Shen, Meng-Hsien. 指導教授：李悅寧博士 Advisor: Yueh-Ning Lee, Ph.D. 中華民國 109 年 8 月 August 2020.

(2)

(3) 中文摘要恆星形成發生在分子雲內的高密度區域，近年來的觀測發現分子雲裡普遍存在絲狀結構 (filamentary structure)，並且與前恆星核 (prestellar core) 形成有高度相關。Monocerous R2 (Mon R2) 分子雲區域的觀測資料呈現有類似面狀的結構 (sheet-like structure)，並且有數條的絲狀體 (filament) 往中心匯聚。我們的研究是想用物理模型來了解這種面狀結構的動力機制及其演化過程。我們考慮在其自身重力下 (self-gravity)，這個片狀結構的徑向塌縮動態，並解其自相似解 (self-similar solution)。我們畫出表面密度與速度的徑向剖面圖和Mon R2的觀測值做比對。接著我們使用微擾分析 (perturbative analysis) 研究絲狀體在這種結構下的發展，並比較生成絲狀體的速率與大尺度時間下塌縮的速率，並得到絲狀體的特徵數量與每條絲狀體的線質量密度。關鍵字 :恆星形成，星際介質，絲狀體，集中絲狀體系統，前恆星核，流體力學，自相似解，微擾分析. i.

(4) ii.

(5) Abstract Star formation occurs inside over-dense regions in molecular clouds. In recent years, observations have found filamentary structures in molecular clouds, which are highly associated with the prestellar cores. Observation of the Monocerous R2 region (Mon R2) suggests a sheet-like structure and shows a structure with several filaments converging into the central hub. We study the dynamics of this sheet-like structure and its evolution with a physical model. We consider the radial collapse of a sheet under self-gravity and solve the self-similar solutions. We produce surface density and velocity profiles which are compatible with observations of Mon R2. Finally, we study the development of filaments in this structure using perturbative analysis. By comparing the timescales of the filament growth rate and the global collapse, it will allow to obtain the characteristic number of filaments and the line mass of each filament.. Keywords: Star formation, Interstellar medium (ISM), Filament, Hub-filament system, Prestellar core, Hydrodynamics, Self-similar solution, Perturbative analysis. iii.

(6) iv.

(7) Contents 中文摘要. i. Abstract. iii. contents. vi. List of Figures. vii. 1. Introduction. 1. 1.1. Star formation and the interstellar medium . . . . . . . . . . . . . . . . . . . .. 1. 1.2. Fluid dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 5. 1.3. Collapse criterion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 7. 1.3.1. Sphere (γ < 34 ) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 7. 1.3.2. Cylinder (γ < 1) . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 7. Self-similar solution: example of a spherical collapse . . . . . . . . . . . . . .. 8. 1.4.1. Physical equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 8. 1.4.2. Similarity solution . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 9. 1.4. 2. Collapse solution of a self-gravitating sheet. 11. 2.1. Thin sheet approximation . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 11. 2.1.1. Vertical equilibrium . . . . . . . . . . . . . . . . . . . . . . . . . . .. 11. 2.1.2. Radial collapse equations . . . . . . . . . . . . . . . . . . . . . . . . .. 12. 2.1.3. Collapse criterion γ <. 5 4. . . . . . . . . . . . . . . . . . . . . . . . . .. 13. 2.2. The equation set of self-similar collapse . . . . . . . . . . . . . . . . . . . . .. 14. 2.3. Treatment of self-gravity in a flat system . . . . . . . . . . . . . . . . . . . . .. 15. 2.3.1. Monopole approximation . . . . . . . . . . . . . . . . . . . . . . . . .. 15. 2.3.2. Exact gravitational field from the surface density . . . . . . . . . . . .. 16. Mathematical analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 18. 2.4.1. The sonic line. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 18. 2.4.2. Classification of sonic line . . . . . . . . . . . . . . . . . . . . . . . .. 18. 2.4. v.

(8) 2.5 2.6 3. 4. 2.4.3. Eigensolutions through a sonic point . . . . . . . . . . . . . . . . . .. 19. 2.4.4. Asymptotic solutions . . . . . . . . . . . . . . . . . . . . . . . . . . .. 19. Numerical solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 22. 2.5.1. Monopole approximation . . . . . . . . . . . . . . . . . . . . . . . . .. 22. Physical domain . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 23. Perturbative analysis. 27. 3.1. Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 27. 3.2. Treatment of self-gravity in a flat system . . . . . . . . . . . . . . . . . . . . .. 30. 3.2.1. Approximation from local surface density . . . . . . . . . . . . . . . .. 30. 3.2.2. Exact gravitational field from the surface density . . . . . . . . . . . .. 30. Discussion and conclusions. 33. Bibliography. 35. vi.

(9) List of Figures 1.1. The column density maps of two subfields in Aquila and Polaris regions . . . .. 2. 1.2. A schematic diagram of the process of forming the prestellar core. . . . . . . .. 2. 1.3. Three-color image of the Mon R2 . . . . . . . . . . . . . . . . . . . . . . . .. 3. 1.4. The three map of the Mon R2 . . . . . . . . . . . . . . . . . . . . . . . . . . .. 4. 1.5. The Position-velocity diagrams along the filament F3 in the Mon R2 . . . . . .. 4. 1.6. The (x, y) plane showing the structure of complete solutions. . . . . . . . . . .. 10. 2.1. Iteration procedure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 15. 2.2. The velocity profile v(x) in monopole approximation . . . . . . . . . . . . . .. 22. 2.3. The density profile σ(x) in monopole approximation . . . . . . . . . . . . . .. 23. 2.4. The Σ profile for several moments of t . . . . . . . . . . . . . . . . . . . . . .. 24. 2.5. Compare with the velocity profile . . . . . . . . . . . . . . . . . . . . . . . .. 24. 2.6. The velocity profile for several moments of t . . . . . . . . . . . . . . . . . . .. 25. 2.7. Compare with the Σ profile . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 25. 3.1. Perturbative iteration procedure . . . . . . . . . . . . . . . . . . . . . . . . . .. 30. vii.

(10) viii.

(11) Chapter 1. Introduction. 1.1. Star formation and the interstellar medium. Star formation occurs inside over-dense region in molecular clouds. The interstellar medium (ISM) is composed of gas and dust, and the primary element is hydrogen, followed by helium. The ISM is affected by turbulence, gravity and magnetic field, forming complex structures. In particular, the molecular clouds are formed in high-density regions of the ISM. Within the molecular cloud, complex structures also develop as consequence of interaction among the physical mechanisms mentioned above. Early observations found that stars are born in molecular clouds. Moreover, in recent years, observations by space telescopes have found filamentary structures, which are highly associated with the prestellar cores (André et al. 2010) . Figure 1.1 shows two molecular cloud regions in the Gould Belt. In Aquila, there are many high density regions, which are characterized as prestellar cores and protostars. However, the Polaris Flare region does not show any sign of star formation activities. Two regions both have filamentary structures. This suggests that filaments form prior to star formation. 1.

(12) Ph. André et al.: The Herschel Gould Belt Survey. Fig. 1. Column density maps of two subfields in Aquila (left) and Polaris (right) derived from our SPIRE/PACS data. The contrast of the filaments with respect to the non-filamentary background has been enhanced using a curvelet transform as described in Appendix A. Given the typical width ∼10 000 AU of the filaments, these column density maps are equivalent to maps of the mass per unit length along the filaments. The color scale on the right of each panel is given in approximate unitsin of Aquila the critical (left) line massand of Inutsuka & Miyama (1997) as discussed in Sect. 4. Figure 1.1:shown Column density maps of two subfields Polaris (right) derived from SPIRE/PACS The areas where the filaments have a mass per unit length larger than half the critical value and are thus likely gravitationally unstable have been data (Fig.1highlighted from Andr´ e et al. 2010) . The green stars and blue triangles show the candidate Class 0 protostars and in white. The maximum line mass observed in the Polaris region is only ∼0.45 × the critical value, suggesting that the Polaris filaments are stable and unable to form stars at the present time. The candidate Class 0 protostars and bound prestellar cores identified in Aquila by Bontemps bound prestellar cores, respectively, which were identified in Aquila by Bontemps et al. (2010). et al. (2010) and Könyves et al. (2010) are shown as green stars and blue triangles, respectively. Note the good correspondence between the spatial distribution of the bound cores/protostars and the regions where the filaments are unstable to gravitational collapse.. Prestellar cores form from a sequence of gravitational collapse and fragmentations, that will eventually lead to the formation of stars. The following process describes how a prestellar core forms, and Figure 1.2 is a schematic diagram. First, the ISM contracts under gravity to form a clump. Second, inside a clump, gravitational instability leads to fragmentation into sheets. Third, inside a two-dimensional sheet, gravitational instability leads to fragmentation into filaments. Finally, the fragmentation of these latter results in the formation of prestellar cores.. Fig. 2. Core mass functions (blue histograms with error bars) derived from our SPIRE/PACS observations of the Aquila (left) and Polaris (right) regions, which reveal of total of 541 candidate prestellar cores and 302 starless cores, respectively. A lognormal fit (red curve) and a power-law fit (black solid line) to the high-mass end of the Aquila CMF are superimposed in the left panel. The power-law fit has a slope of −1.5 ± 0.2 (compared to a Salpeter slope of −1.35 in this dN/dlogM format), while the lognormal fit peaks at ∼0.6 M$ and has a standard deviation of ∼0.43 in log10 M. The IMF of single stars (corrected for binaries – e.g., Kroupa 2001), the IMF of multiple systems (e.g., Chabrier 2005), and the typical mass spectrum of CO clumps (e.g., Kramer et al. 1998) are also shown for comparison. Note the remarkable similarity between the Aquila CMF and the stellar IMF, suggesting a ∼ one-to-one correspondence between core mass and star/system mass with M!sys = " Mcore and " ≈ 0.4 in Aquila.. one-to-one basis, with a fixed and relatively high local efficiency, i.e., "core ≡ M! /Mcore ∼ 20−40% in Aquila. This is consistent with theoretical models according to which the stellar IMF is in large part determined by pre-collapse cloud fragmentation, prior to the protostellar accretion phase (cf. Larson 1985; Padoan & Nordlund 2002; Hennebelle & Chabrier 2008). There are several caveats to this simple picture (cf. discussion in André et al. 2009), and detailed analysis of the data from the whole GBS will be required to fully characterize the CMF–IMF relationship and, e.g., investigate possible variations in the efficiency "core with. environment. It is nevertheless already clear that one of the keys to the problem of the origin of the IMF lies in a good understanding of the formation process of prestellar cores, even if additional processes, such as rotational subfragmentation of prestellar cores into binary/multiple systems (e.g., Bate et al. 2003), probably also play an important role. Our Herschel initial results also provide key insight into the core formation issue. They support an emerging picture (see also Myers 2009) according to which complex networks of long, thin filaments form first within molecular clouds, possibly as a result Page 3 of 7. Figure 1.2: A schematic diagram of the process of forming the prestellar core.. However, it is highly unlikely to form a high-mass star from fragmentation of one single filament, because there is not enough mass. Observations suggest that star formation occurs preferentially along the filaments, with high-mass stars forming in the highest density regions where several filaments converge, called ridges or hubs (Treviño-Morales et al. 2019). This structure is usually called the hub-filament system, where high-mass stars and star clusters form. 2.

(13) In the following, we will give an example of such structure (Treviño-Morales et al. 2019). The Monocerous R2 (Mon R2) is an excellent target to study hub-filament system, at a distance of 830 pc. Figure 1.3 is Mon R2, shows the cloud with several filaments converging into the central hub. In the next paragraph, we will briefly introduce the hub-filament system from the observations of the Mon R2.. Figure 1.3: Left: Three-color image of the Mon R2 cluster-forming hub-filaments system.Red: H2 column density map derived from Herschel SPIRE and PACS observations (Didelon et al. 2015), green: 1.65 µm band of 2MASS (Two micron all survey; Skrutskie et al. 2006), and blue: 560 nm band of DSS (Digitalized Sky Survey; Lasker et al. 1990). (Fig.1 from Treviño-Morales et al. 2019). In Figure 1.4, the left column shows the integrated intensity maps over the whole surveyed area for the (1 → 0) transition lines of. 13. CO and C 18O molecules. It is similar to the map of. column density, indicating that the mass is concentrated in the central hub and the filaments. The middle column shows the velocity gradient, in which we can see that the upper left is moving toward us and the lower right is moving away from us. This can be explained as a radially collapsing velocity field within a plane with an inclination angle with respect to the line of sight. The right column shows high velocity dispersion in the central hub. This is possibly a sign of gravitational collapse or outflow from massive stars. Figure 1.5 shows the relationship between velocity and position along the filament. The velocity increases along the filament toward the central hub, suggesting that the gas is falling into the central hub. Overall, Fig.1.4 and Fig.1.5 present the Mon R2 region with several filaments converging into the central hub, and the filaments lie within a 2D plane. 3.

(14) Figure 1.4: Left panels show the integrated intensity maps over the whole surveyed area for the (1 → 0) transition lines of the 13CO and C 18O molecules. Middle panels present the velocity centroid. Right panels show the linewidth. The yellow labels, and the dotted lines, indicate the main features identified in the region. (Fig.3 from Treviño-Morales et al. 2019). Figure 1.5: Position-velocity diagrams along the ‘skeleton’ of filament F3 obtained from the 13CO (left) and C 18O (right) data cubes. The vertical yellow dashed lines indicate the transition between the hub and the filaments, corresponding to radii 20000 , 25000 (Rhub), and 30000 .( Fig.B.3 from Treviño-Morales et al. 2019). In this study, we aim to describe the dynamics of this system with a physical model. Therefore, we first review the study related to self-similar collapse solutions have been studied for an isothermal spheres in the ISM (e.g. Hunter 1977; Larson 1969; Penston 1969; Shu 1977), polytropic cylinder/filament (e.g. Kawachi & Hanawa 1998), and magnetized molecular cloud cores which are flattened (e.g. Contopoulos et al. 1998; Krasnopolsky & Königl 2002). None of the existing studies has considered the self-similar collapse of a thin sheet in the ISM. Therefore, 4.

(15) we consider the radial collapse of a sheet with self-gravity and solve for the self-similar solution. Whitworth & Summers (1985) found a family of solution (by Hunter 1977; Larson 1969; Penston 1969; Shu 1977), and did mathematical and physical analysis, as well as numerical analysis. Therefore, we follow their procedure in our work, and compare the obtained surface density and velocity profiles to observations in Chapter 2. In addition, we would like to study the development of filaments in this structure with perturbative analysis. Therefore, we review the study related to perturbative analysis. Hanawa & Matsumoto (1999) studied the spherical harmonics perturbation. They considered the stability of the Larson-Penston (Larson 1969; Penston 1969) solution against nonspherical perturbations. Ledoux (1951) studied the self-gravitating isothermal sheet plane perturbation. He did the horizontal plane wave perturbation in uniform sheet, that is, δΣ(x, z, t) in Σ(z), where z is the normal vector of the sheet. In contrast, Simon (1965) studied the vertical polytropic perturbation in uniform sheet, that is, δΣ(z, t) in Σ(z). There exists also perturbative studies of a rotating protoplanetary disk (e.g McKee 1991; Shi-Xue et al. 2000; Wu et al. 1995). However, a disk in Keplerian rotation is very different from our configuration, because we can have vφ ≈ 0 and. they have vr ≈ 0. In our case, we consider azimuthal perturbation in a thin sheet with the radial profile, that is, δΣ(r, θ, t) in Σ(r, t). This is described in Chapter 3. This study will give us the. numbers of filaments and the line mass of each filament. In the future, we can study the mass of cores in filaments from the fragmentation of these latter.. 1.2. Fluid dynamics. We consider Navier-Stokes Equations, which describes the dynamics of fluids. The conservation of mass equation, also called continuity equation, writes ∂ρ ~ + ∇ · (ρ~u) = 0, ∂t. (1.2.1). where ρ is the density, t is time, ~u is the velocity. The conservation of momentum equation is ∂(ρ~u) ~ ~ − ρ∇φ, ~ + ∇ · (ρ~u~u) = −∇P ∂t. (1.2.2). where P is pressure, φ is the gravitational potential. The conservation of energy equation is ∂ρε ~ ~ · ~u − ρ(∇φ) ~ · ~u + ∇ ~ · KH ∇T ~ + Q − Λ, + ∇ · (ρε~u) = −(∇P) ∂t 5. (1.2.3).

(16) where KH is heat conduction coefficient, T is the temperature, Q is radiative heating, Λ is a cooling function, and ε is specific energy equaling to u2 P + , γ−1 2 where γ is adiabatic index. Because the ISM is almost isothermal at temperature T equal to 10 Kelvin, we can replace the energy equation with an Equation of state P = nkT =. ρkT = ρc2 , µm p. (1.2.4). where k is Boltzmann constant, µ is mean molar weight, m p is proton mass, c = 200 m/s is the thermal sound speed. Finally, the Poisson equation describes the relation between the gravitational potential and the density: ~ · ∇φ ~ = 4πGρ. ∇. (1.2.5). Since we consider sheet-like structures in this study, we express Navier-Stokes equations in the cylindrical coordinate (r, θ, z). The continuity equation ∂ρ 1 ∂ 1 ∂ ∂ (rρur ) + (ρuθ ) + (ρuz ) = 0. + ∂t r ∂r r ∂θ ∂z. (1.2.6). The momentum equations ∂ur ∂ur ρuθ ∂ur ∂ur ρu2θ −∂P ∂φ ρ + ρur + + ρuz − = −ρ , ∂t ∂r r ∂θ ∂z r ∂r ∂r ∂uθ ∂uθ ρuθ ∂uθ ∂uθ ρur uθ −1 ∂P ρ ∂φ ρ + ρur + + ρuz + = − , ∂t ∂r r ∂θ ∂z r r ∂θ r ∂θ ∂uz ∂uz ρuθ ∂uz ∂uz −∂P ∂φ ρ + ρur + + ρuz = −ρ . ∂t ∂r r ∂θ ∂z ∂z ∂z. (1.2.7) (1.2.8) (1.2.9). The Poisson equation 2 2 ~ · ∇φ ~ = 1 ∂ (r ∂φ ) + 1 ∂ φ + ∂ φ = 4πGρ. ∇ r ∂r ∂r r2 ∂θ2 ∂z2. 6. (1.2.10).

(17) 1.3 1.3.1. Collapse criterion Sphere (γ < 43 ). If we consider a sphere of polytropic (P ∝ ργ ), self-gravitating gas, the mass of sphere is 4 M = πr3 ρ ∝ r3 ρ 3. ⇒ ρ ∝ r−3 .. (1.3.1). The internal energy is Eint = PV ∝ r3 ργ .. (1.3.2). The gravitational energy is Egrav =. GM 2 M 2 ∝ ∝ r5 ρ2 . r r. (1.3.3). In order to compare the relationship between internal energy and gravitational energy, we divide internal energy by gravitational energy, and substitute the relationship of (1.3.1) into (1.3.2) and (1.3.3) r3 ργ r−3γ+3 Eint ∝ 5 2 ∝ −1 = r−3γ+4 = rn . Egrav r ρ r. (1.3.4). If n > 0, the gravitational potential energy has a chance to dominate over pressure support when r decreases. On the contrary, if n < 0, the internal energy will eventually dominate when r shrinks to very small value. The collapse will be stopped by the thermal pressure in this case. There exists a critical value for the polytropic index, γ, such that collapse is possible : 4 −3γ + 4 > 0 ⇒ γ < γcrit = . 3 An isothermal gas has γ = 1. Therefore, the collapse solution exists and we can solve the self-similar equation. 1.3.2. Cylinder (γ < 1). If we consider a sphere of polytropic (P ∝ ργ ), self-gravitating gas, the mass of cylinder is M = ρπr2 L ∝ r2 ρ,. ⇒ ρ ∝ r−2 ,. (1.3.5). where L is the height of cylinder. The internal energy is Eint = PV = ργ c2 πr2 L ∝ r2 ργ .. (1.3.6). The gravitational energy is Egrav =. GM 2 ∝ M 2 ∝ r 4 ρ2 . L 7. (1.3.7).

(18) In order to compare the relationship between internal energy and gravitational energy, we divide internal energy by gravitational energy, and substitute the relationship of (1.3.5) into (1.3.6) and (1.3.7) r2 ργ r−2γ+2 Eint ∝ 4 2 ∝ = r2−2γ = rn . Egrav r ρ r0. (1.3.8). If n > 0, the gravitational potential energy has a chance to dominate over pressure support when r decreases. On the contrary, if n < 0, the internal energy will eventually dominate when r shrinks to very small value. The collapse will be stopped by the thermal pressure in this case. There exists a critical value for the polytropic index, γ, such that collapse is possible : 2 − 2γ > 0 ⇒ γ < 1. In the cylinder, the solutions of γ < 1 can be collapsed.. 1.4. Self-similar solution: example of a spherical collapse. In the following, we briefly describe the self-similar solution of a collapsing sphere (Shu 1977). The self-similar formalism allows us to significantly simplify the complex fluid dynamics equations. 1.4.1. Physical equation. The independent physical variables are radius r and time t. The dependent physical variables are mass interior to M(r, t); density ρ(r, t); outward radial flow velocity u(r, t); and pressure P(r, t). The conservation of mass equation is ∂M ∂M +u = 0, ∂t ∂r. ∂M = 4πr2 ρ. ∂r. (1.4.1). The relations (Eq.1.4.1) are equivalent to the continuity equation, ∂ρ 1 ∂ 2 + r ρu = 0, ∂t r2 ∂r. (1.4.2). ∂u ∂u 1 ∂P GM +u =− − 2 , ∂t ∂r ρ ∂r r. (1.4.3). and the momentum equation is. where G is the universal gravitational constant. The gas is isothermal, so we can express P = c2 ρ, with constant, uniform isothermal sound speed c. 8. (1.4.4).

(19) 1.4.2. Similarity solution. Following Shu (1977), we define the independent dimensionless similarity variable x as x=. r . ct. (1.4.5). We also define the dependent dimensionless similarity variables α, m and v as α(x) ρ(r, t) = , 4πGt2. c3 t M(r, t) = m(x), G. u(r, t) = cv(x).. (1.4.6). The substitution of equations (1.4.6) into equations (1.4.1) now yields the ordinary differential equations m−x. dm dm +v = 0, dx dx. dm = x2 α. dx. (1.4.7). The term dm/dx can be eliminated from the above relations to give m = x2 α(x − v).. (1.4.8). This formula plus some straightforward manipulation now allows us to express equations (1.4.2) and(1.4.3) as the coupled set of ordinary differential equations (ODE) # h i dv " 2 2 (x − v) − 1 = α(x − v) − (x − v), dx x # h i 1 dα " 2 2 = α − (x − v) (x − v). (x − v) − 1 α dx x. (1.4.9) (1.4.10). The coupled set of ODEs (1.4.9 and 1.4.10) are solved by numerical integration in the study of Whitworth & Summers (1985). Figure 1.6 shows the velocity profile for a family of solutions : ρ0 , where y is equivalent to v. This figure is divided into three parts with three asymptotic solutions of Whitworth & Summers (1985). 1. The asymptotic form as x → −0 (t → −∞) is v(x) = 2x/3 + (2 − 3α0 ) x3 /135 + O x5 , α(x) = α0 + α0 (2 − 3α0 ) x2 + O x4 ,. (1.4.11). with a flat central density α0 . When x → −0, the value of velocity v goes to zero and the value of density α is flat at the center.. 2. The asymptotic form as x → +0 (t → +∞) is v(x) = − (2m0 /x)1/2 − (3 ln [ε0 x] /2) (x/2m0 )1/2 + O(x ln[kx]), 1/2 h i α(x) = m0 /2x3 − 3 ln e2 ε0 x /4 (2m0 x)−1/2 + O(ln[kx]),. (1.4.12). with central point mass m0 . When x → +0, the value of velocity v is proportional to −x−1/2 and the value of density α is proportional to x−3/2 .. 9.

(20) 3. The asymptotic form as x → ∓∞ (t → ∓0) is v(x) = v∞ − (A − 2) /x + v∞ /x2 + [4v∞ + (A − 2) (A − 6)] /6x3 + O x−4 , α(x) = A/x2 − A (A − 2) /2x4 + O x−6 ,. (1.4.13). with uniform radial velocity v∞ at large radius and A = αx2 . When x → ∓∞, the value of velocity v is constant and the value of density α is proportional to x−2 .. They started by integrating from −0 to minus infinity with the asymptotic solutions at x →. −0. However, the integration cannot be done directly through the sonic point, because it is a mathematical singular point. They found several slopes after the sonic point, and integrated to minus infinity. Finally, they used the transformations x = −x∞ → +x∞ , v = v∞ + (A − 2)/x∞ → 198 5MNRAS.214. v∞ − (A − 2)/x∞ to map minus infinity x directly to x at plus infinity, and integrated from plus 4. infinity to +0 .. A. Whitworth and D. Summers. H 1 T r~T —£ ‘saiojduiKs* Suisiruo. © Royal Astronomical Society • Provided by the NASA Astrophysics Data System. Section 6.1 for details.. Figure 1.6: The (x, y) plane showing the structure of complete solutions and their relationship to the sonic line.(Whitworth & Summers (1985) Fig.1). 10.

(21) Chapter 2. Collapse solution of a self-gravitating sheet 2.1. Thin sheet approximation. The goal of this study is to derive the collapse solution which describes a system similar to Mon R2. Since this region is flat, we search for solutions of radial collapse of a flat system, using the thin sheet approximation (e.g. Schmitz 1990, 1994). 2.1.1. Vertical equilibrium. If we consider vertical hydrostatic equilibrium, the density of a uniform sheet can be found from solving the Poisson equation in the z-direction ∂P + ρg = 0, ∂z. g = −2πGΣ(z),. z ρ(z) = ρ0 sech2 ( ), H. (2.1.1) (2.1.2). where the scale height c H= p . 2πGρ0. (2.1.3). We can integrate in the vertical direction, z and obtain the surface density of z Zz Σ(z) = −z. z ρ(z0 )dz0 = 2ρ0 H tanh( ). H. (2.1.4). When z is large, r Σ(z) ≈ Σ = 2ρ0 H =. 2cρ0 . πG. When the sheet is not uniform, the above equilibrium becomes invalid. However, if . ∂Σ rˆ + 1 ∂Σ θˆ H(r, θ), ∂r r ∂θ 11. (2.1.5). (2.1.6).

(22) we can assume local vertical hydrostatic equilibrium, the density of a thin sheet can then be approximate as: 2. ρ(r, θ, z) = ρ0 (r, θ) sech. ! z , H(r, θ). (2.1.7). where the local scale height H(r, θ) =. c c2 = p . πGΣ(r, θ) 2πGρ0 (r, θ). (2.1.8). The above assumption is valid only when we are close to mid-plane, that is, z . H. On the contrary, when z H, we can approximate the relation between ρ and the surface density Σ as ρ(r, θ, z) = Σ(r, θ)δ(z),. (2.1.9). where δ(z) is the Dirac delta function.. 2.1.2. Radial collapse equations. The thin sheet description allows to eliminate dependancies in z. We therefore have scale height H(r, t); density ρ(r, t); surface density Σ(r, t); total mass M(r, t); and velocity u(r, t) as functions of radius r and time t. The mass equation becomes Z∞ Zr. 0. 0. ρ2πr dr dz =. M(r, t) = M(0, t) + z=−∞ r0 =0. ∂M ∂M +u = 0, ∂t ∂r. Zr. 2πΣr0 dr0 ,. r0 =0. ∂M = 2πrΣ. ∂r. (2.1.10). Integrating the density ρ in z-direction in the equation (1.2.6). The integration of ∂(ρuz )/∂z becomes zero because ρ is small for z at infinity. The relations (Eq.2.1.10) are equivalent to the continuity equation ∂Σ 1 ∂ 1 ∂ (Σuθ ) = 0. + (rΣur ) + ∂t r ∂r r ∂θ. (2.1.11). In this section, we do not consider variation in theta. Therefore, we neglect dependencies in θ for the moment, and rewrite ur to u. The continuity equation becomes ∂Σ 1 ∂ + (rΣu) = 0. ∂t r ∂r. (2.1.12). We consider a non-rotating sheet, that is, uθ =0. Therefore, we only consider the momentum equations in the r direction (i.e. Eq. 1.2.7). The momentum equation becomes ∂u ∂u 1 ∂P ∂φ +u =− − . ∂t ∂r ρ ∂r ∂r 12. (2.1.13).

(23) We express pressure P with the relation (Eq.1.2.4)(i.e. P = ρc2 ), and the integration along z yields 1 ∂P c2 ∂P 2c2 ∂Σ = = . ρ ∂r ρ ∂r Σ ∂r. (2.1.14). We express gravitational acceleration −∂φ/∂r with g(r, t). The momentum equation becomes ∂u ∂u −2c2 ∂Σ +u = + g. ∂t ∂r Σ ∂r. (2.1.15). The Poisson still needs to be solved in the 3D space, with θ-dependent term dropped, and we express the relation (Eq.2.1.9) 1 ∂ ∂φ ∂2 φ (r ) + 2 = 4πGρ = 4πGΣδ(z). r ∂r ∂r ∂z 2.1.3. Collapse criterion γ <. (2.1.16). 5 4. If we consider a sphere of polytropic (P ∝ ργ ), self-gravitating gas, the mass of a thin sheet is 1. M = ρπr2 H ∝ ρ 2 r2. ⇒ ρ ∝ r−4 .. (2.1.17). where ρ is the surface density and H is function of r c H(r) = p 2πGρ0 (r). 1. ⇒ H ∝ ρ− 2 .. The internal energy is 1. Eint = PV ∝ ργ r2 H ∝ r2 ργ− 2 .. (2.1.18). The gravitational potential energy is Egrav =. GM 2 M 2 ∝ ∝ r3 ρ. r r. (2.1.19). In order to compare the relationship between internal energy and gravitational energy, we divide internal energy by gravitational energy, and substitute the relationship of (Eq.2.1.17) into (Eq.2.1.18) and (Eq.2.1.19) 1. Eint r2 ργ− 2 r−4γ+4 ∝ 3 ∝ −1 = r−4γ+5 = rn . Egrav rρ r. (2.1.20). If n > 0, the gravitational potential energy has a chance to dominate over pressure support. On the contrary, if n < 0, the internal energy will eventually dominate when r shrinks to very small value. The collapse will be stopped by the thermal pressure in this case. There exists a critical value for the polytropic index, γ, such that collapse is possible : 5 −4γ + 5 > 0 ⇒ γ < . 4 An isothermal gas has γ = 1. Therefore, the collapse solution exists and we can solve the self-similar equation. 13.

(24) 2.2. The equation set of self-similar collapse. In the following, we derive the self-similar formalism of the radial collapse of a thin sheet. Navier-Stokes equations (2.1.12) and (2.1.15), and Poisson equation (2.1.16) ∂Σ 1 ∂ + (rΣu) = 0, ∂t r ∂r ∂u ∂u −2c2 ∂Σ Momentum equation: +u = + g, ∂t ∂r Σ ∂r ∂2 φ 1 ∂ ∂φ (r ) + 2 = 4πGρ = 4πGΣδ(z). Poisson equation: r ∂r ∂r ∂z Continuity equation:. (2.2.1) (2.2.2) (2.2.3). Following Shu (1977), we define the independent dimensionless similarity variable x, ζ as x=. r , ct. ζ=. z . ct. (2.2.4). We also define the dependent dimensionless similarity variables α, σ, v, m, γ and ψ as ρ(r, t) =. α(x) , 4πGt2. u(r, t) = cv(x),. Σ(r, t) =. c σ(x), 2πGt. c g(r, t) = γ(x), t. M(r, t) =. c3 t m(x) G. φ(r, t) = c2 ψ(x).. (2.2.5). The substitution of equations (2.2.5) into equations (2.1.10) now yields the ordinary differential equations m(x) − x. dm dm + v(x) = 0, dx dx. dm = xσ(x). dx. (2.2.6). The term dm/dx can be eliminated from the above relations to give m(x) = x(x − v(x))σ(x).. (2.2.7). The substitution of equations (2.2.5) into equations (2.2.1 , 2.2.2 and 2.2.3) yields the ordinary equations dv (x − v(x)) dσ (x − v(x)) = + , dx σ dx x dσ σ(x)(x − v(x)) dv σ(x) = + γ(x), dx 2 dx 2 1 ∂ ∂ψ(x) ∂2 ψ(x) (x )+ = α(x) = 2σ(x)δ(ζ(x)). x ∂x ∂x ∂ζ(x)2. (2.2.8) (2.2.9) (2.2.10). Rearrange these two equations(2.2.8) and (2.2.9), we obtain the coupled set of ordinary differential equation (ODE) 2 = (x − v)( + γ), dx x h i 1 dσ (x − v)2 2 − (x − v)2 = + γ, σ dx x. h. 2 − (x − v)2. i dv. which should be solved together with the Poisson equation (2.2.10). 14. (2.2.11) (2.2.12).

(25) 2.3. Treatment of self-gravity in a flat system. In this section, we would like to solve the coupled set of ODEs (2.2.11 and 2.2.12) by numerically integrating from small radius to infinity. However, we cannot solve the Poisson equation directly without knowing the entire σ profile in a cylindrical geometry. Therefore, we first apply monopole approximation to the gravity field. This allows us to perform a first integration of an approximated solution. We can then calculate the exact solution by method of iteration. Once σ(x) is known, γ(x) can be calculated accordingly. We can then re-integrate Equations (2.2.11 and 2.2.12) to obtain a new profile of σ(x), which is closer to the real value. This process should be repeated until convergence is reached. The profile of γ(x) can be calculated from σ with the Hankel transform.. Figure 2.1: Iteration procedure in monopole approximation. 2.3.1. Monopole approximation. In this subsection, we solve the coupled set of ODEs using monopole approximation. The continuity equation is the same as equation(2.2.1). We approximate the gravitational acceleration by that of a point mass equaling to the mass enclosed within r: g(r, t) = −. GM(r, t) ∂φ(r, t) ≈− . ∂r r2. (∵ φ ≈. GM(r, t) ) r. (2.3.1). The momentum equation (2.2.2) becomes ∂u ∂u −2c2 ∂Σ GM +u = − 2 . ∂t ∂r Σ ∂r r. (2.3.2). The substitution of equations (2.2.5) and (2.2.7) into equations (2.1.12) and (2.3.2) yields the ordinary equations dv (x − v(x)) dσ (x − v(x)) = + , dx σ dx x dσ σ(x)(x − v(x)) dv σ(x)2 = − (x − v(x)). dx 2 dx 2x Rearrange these two equations, we obtain the coupled set of ODEs h i dv 1 2 − (x − v)2 = (x − v) [2 − σ(x − v)] , dx x h i 1 dσ 1 2 − (x − v)2 = (x − v) [(x − v) − σ] . σ dx x 15. (2.3.3) (2.3.4). (2.3.5) (2.3.6).

(26) This equation set can be then integrated numerically given proper boundary conditions. 2.3.2. Exact gravitational field from the surface density. The following process is the use of the Bessel function to solve Poisson equation and the Hankel transform of these functions. The function Jk (sx)e−s|z| cos(kθ) satisfies the Poisson equation(1.2.10) for all values of s. We define −s|z| ˆ ψ s (x, θ, z) = sψ(s)J cos(kθ), k (sx)e. (2.3.7). such that Z∞ ψ(x, θ) =. ψ s (x, θ, z = 0)ds. (2.3.8). 0. is a linear combination with all possible values of s. Correspondingly, σ s (x, θ) = sσ(s)J ˆ k (sx) cos(kθ),. (2.3.9). such that Z∞ σ(x, θ) =. σ s (x, θ)ds.. (2.3.10). 0. We can verify that the above relations give us the Hankel transform pairs: Z∞. Hk. Z∞. sds σ(s) ˆ Jk (sx) ⇐⇒ σ(s) ˆ =. σ(x) = 0. (2.3.11). xdx ψ(x)Jk (sx).. (2.3.12). 0. Z∞. Hk. Z∞. ˆ ˆ = sds ψ(s)J k (sx) ⇐⇒ ψ(s). ψ(x) =. xdx σ(x) Jk (sx),. 0. 0. The Bessel solution for k = 0 inserts ψ s (Eq.2.3.7) and σ s (Eq.2.3.9) into the Poisson equation, we can find the relation between ψˆ and σ. ˆ When |ζ| > 0, ! 2 ∂J ∂ J s 0 (sx) 0 (sx) 2 2 ˆ ψ(s) +s + s J0 (sx) = 0. x ∂(sx) ∂(sx)2. (2.3.13). When ζ = 0, ! 2 s ∂J0 (sx) 2 ∂ J0 (sx) ˆ = 2σ(s)J +s − 2sJ0 (sx)δ(ζ) ψ(s) ˆ 0 (sx)δ(ζ). x ∂(sx) ∂(sx)2. (2.3.14). The first two items are small with respect to the third item, so we obtain ˆ −2sψ(s)J ˆ 0 (sx)δ(ζ) = 2σ(s)J 0 (sx)δ(ζ), 16. (2.3.15).

(27) such that σ(s) ˆ ˆ =− ψ(s) = s. Z∞ xdx ψ(x)J0 (sx).. (2.3.16). 0. Finally, the gravitational acceleration at the monopole approximation can be calculated: ∂ψ = γ(x) = − ∂x. Z∞. s2 ˆ [J1 (sx) − J−1 (sx)] ds ψ(s) 2. 0 Z∞. ˆ s2 ds ψ(s) J1 (sx). =. (∵ J−1 = −J1 ). 0. Z∞ =−. sds σ(s) ˆ J1 (sx).. (2.3.17). 0. We consider an axisymmetric thin sheet case (k = 0). We summarize the conversion from σ into γ(x) through the calculation of Hankel transform. −H1−. H0. σ(x) −→ σ(s) ˆ −→ γ(x). 17.

(28) 2.4. Mathematical analysis. In this section, we need to understand the mathematical analysis of monopole approximation before we would like to solve the coupled set of ODE (2.2.11 and 2.2.12) by numerically integrating. We refer to the method of mathematical analysis of Whitworth & Summers (1985). 2.4.1. The sonic line. We replace equations (2.3.5) and (2.3.6) with the following formula dv B(x, v, σ) = ; dx A(x, v). dσ C(x, v, σ) = ; dx A(x, v). (2.4.1). A = 2 − (x − v)2 ;. (2.4.2). 1 B = (x − v)(2 − σ(x − v)); x. (2.4.3). σ (x − v)((x − v) − σ). x. (2.4.4). C=. From equation (2.4.2)-(2.4.4) we see that there are two singular surface in (x, v, σ) space given √ by A = 0, or v = x ± 2. Physically acceptable solution can only cross these surfaces through the sonic point where A = B = C = 0. If we exclude parts of the sonic line corresponding to negative mass, that is, when x(x − v) < 0, this means that the coordinates of the sonic points must satisfy ±x s > 0,. vs = xs ∓. √. 2,. √ σ s = ± 2.. (2.4.5). The ±/∓ convention in equation (2.4.5) is defined such that the upper sign will always refer to. sonic points with x s > 0, and therefore t > 0; the lower sign will always refer to sonic with x s < 0, and therefore t < 0. 2.4.2. Classification of sonic line. Since the sonic points are singular, we can find the topological structure of paths through a sonic point by investigating the eigenvalues of the matrix  ∂A ∂A ∂A   √  √   ∓ 2 ± 2 0    ∂x ∂v ∂σ              1  1 1  ∂B ∂B ∂B   −  −  M =   =  x  x x s s s  ∂x ∂v ∂σ             ∂C ∂C ∂C   1 1 1     − −  xs xs xs ∂x ∂v ∂σ 18. (2.4.6).

(29) with the partial derivatives evaluated at the sonic point. Because there is a continuous line of sonic points, the determinant of this matrix is zero, i.e. √ √ 2 2 − 2 ) = 0, λI - M = λ(λ2 ± 2λ ± xs xs. (2.4.7). one of the eigenvalues is zero. And then, we neglect this eigenvalue because it does not represent a solution of the physical equations. The other two eigenvalues are given by the quadratic equation λ2 + bλ + c = 0 with √. b = ± 2, and hence. √ √ − 2( 2 ∓ x s ) c= , x2s. √ 2(x s ∓ 2)2 4 d = b − 4c = + 2. 2 xs xs 2. Since d ≥ 0 for all x s , none of the sonic points is a spiral or centre. 2.4.3. Eigensolutions through a sonic point. Linearizing the similarity equation (2.4.1) through (2.4.4) in the vicinity of a sonic point by substituting dv/dx = v0s , dσ/dx = σ0s , x = x s + δ, v = v s + v0s δ, σ = σ s + σ0s δ , and letting δ → 0. Straightforward algebra then yields the gradients of the eigensolutions: q q √ √ √ 2 x s + x s ± 2 2x s + 4 x s + 2 2 + x2s ± 2 2x s + 4 type1, v0s = , σ0s = , 2x s 2x s. type2, 2.4.4. v0s =. xs −. q. √ x2s ± 2 2x s + 4 2x s. , σ0s =. q √ √ x s + 2 2 − x2s ± 2 2x s + 4 2x s. .. (2.4.8). (2.4.9). Asymptotic solutions. 1. x → −0. Substitute x close to −0 into equations(2.3.5) and (2.3.6), and simplify to dv dv dv dv 2v σv2 − x2 + 2xv − v2 = 2 − σx + 2σv − − , dx dx dx dx x x dσ dσ dσ dσ σv2 σ2 v 2 − x2 + 2xv − v2 = σx − 2σv − σ2 + + . dx dx dx dx x x 2. (2.4.10) (2.4.11). Suppose v and σ are equation(2.4.12) and substitute into equations (2.4.10)&(2.4.11)    v = vp xp,      (2.4.12) σ = σq x q ,        p > 0, q ≥ 0. 19.

(30) The following tables show how to solve these equations. 2. 0. x x1. p−1 2v1 2v2. 2. 0. x x1. dv dx. dσ dx. q−1 2σ1 2σ2. −x2. dv dx. p+1. −x2. dσ dx. q+1. +2xv. dv dx. p + p0. 2xv. dσ dx. p + q0. −v2. dv dx. p + p0 + p00 − 1. −v2. dσ dx. p + p0 + q − 1. −σx. +2σv. −. 0 2. q+1. p+q. −σ0. 2v1 σ0. p−1 −2v1 −2v2. σx. −2σv. −σ2. q+1. p+q. q + q0 −σ20 −2σ0 σ1. σ0. 2v x. 2. −2v1 σ0. −. −. p + p0 + q − 1 −v21 σ0. σv2 x. −. q + p + p0 − 1 v21 σ0. σv2 x. σ2 v x. q + q0 + p − 1 σ20 v1 2v1 σ0 σ1 + v2 σ20. We obtain the solution of each term through the above calculation process  σ0 1       v1 = 2 , v2 = − 24 , · · ·   σ2 σ0 5σ30    + ,···  σ1 = − , σ2 = 4 16 96 and we obtain solutions as   1 1   x − σ0 x 2 + · · · v =    2 24     (2.4.13)   !    1 1 2 5 3 2      σ = σ0 − 4 σ0 x + 16 σ0 − 96 σ0 x + · · · However, the monopole approximation is not valid when σ is almost constant at x → −0. There is the gamma actually goes to zero. Substitute x close to −0 into equations(2.2.11) and (2.2.12), and simplify to dv dv dv dv 2v 2 − x2 + 2xv − v2 = 2 + xγ − vγ − , dx dx dx dx x dσ dσ dσ dσ σv2 2 − x2 + 2xv − v2 = σx − 2σv + + σγ. dx dx dx dx x In the same way, we obtain solutions as !   1 1 σ0 3    v= x+ + x + ···    2 64 32        !    σ20 2 1    + x + ···  σ = σ0 +  16 8. (2.4.14) (2.4.15). (2.4.16). 2. x → +0 Substitute x close to +0 into equations (2.3.5) and (2.3.6), and simplify to dv −v = (2 + σv), −v2 dx x dσ σv −v2 = (v + σ). dx x 20. (2.4.17) (2.4.18).

(31) Suppose v and σ are equation(2.4.19) and substitute into equations (2.4.17) and (2.4.18)   −p    v = ax , (2.4.19)     σ = bx−q . and we obtain   2 −3p−1   = −2x−p−1 − abx−2p−q−1 ,  a px     aqx−2p−q−1 = ax−2p−q−1 + bx−p−2q−1 .. (2.4.20). If p = q, the equation (2.4.20) will be established and substituted into equation (2.2.7), and will get the solution as  1    p=q= ,    2  √      a = − 2m0. r ,. b=. m0 . 2. and we obtain solutions as  √ 1   2m0 x− 2 , v = −         r     m0 − 1    σ= x 2.  2. (2.4.21). (2.4.22). 3. x → ∞. Substitute x to ±∞ into equations (2.3.5) and (2.3.6), and simplify to dv = 2 − σx + 2vσ dx dσ = σx − 2σv − σ2 . (2 − x2 + 2xv) dx (2 − x2 + 2xv). (2.4.23) (2.4.24). Suppose v and σ are equation(2.4.25) and substitute into equations (2.4.23)&(2.4.24)    v = v∞ + A0 x−p + A1 x−p−1 + A2 x−p−2 + · · ·      (2.4.25)        σ = B0 x−q + B1 x−q−1 + B2 x−q−2 + · · · and we get relationship as    p = q = 1,       A20 + 10A0 + 16v2∞   A =n , A = 2v , A = ··· (2.4.26) 0 ∞ 1 ∞ 2    6   2  −A0 + 2A0    ···  B0 = 2 − A0 , B1 = 0, B2 = 2 We replace the relationship (2.4.26) back to equation (2.4.25), so we obtain solution as   n∞ 2v∞ 16v2∞ + n2∞ + 10n∞    v = v + + 2 + + ··· ∞  3   x x 6x   (2.4.27)       2 − n∞ n∞ (n∞ − 2)    σ= − + ··· x 2x3 21.

(32) 2.5. Numerical solution. In the section 2.3, we mentioned two steps to solve the coupled set of ODE (2.2.11 & 2.2.12) by numerically integrating from small radius to infinity. Here, we use the asymptotic solutions at x = −0 in the section 2.4 as boundary conditions to initialize the numerical integration. 2.5.1. Monopole approximation. We use the asymptotic solutions at x → −0 (Eq.2.4.13) as boundary conditions to solve the. coupled set of ODE (2.2.11 and 2.2.12). The central surface density, σ0 , is a parameter that. we can choose, but not any value can be successfully integrated. We start by integrating from zero to minus infinity. However, when we integrate to the sonic point, we cannot continue to integrate, because this is a mathematical singular point. There for we use the slopes found in (2.4.8 and 2.4.9) to restart the integration after the sonic point. After that, we map minus infinity x directly to x at plus infinity, and then we integrate from plus infinity to zero. This corresponds to the moment when t passes through zero. We show the velocity profile v(x) in Fig. 2.2, and the density profile σ(x) in Fig. 2.3. Figure 2.2 shows the results of v(x) profile for several σ0 values. When x → −0, the value of. velocity goes to zero. When x → ±∞, the value of velocity is almost constant. When x → +0, the value of velocity is proportional to −x−1/2 , and it becomes singular at x = 0.. x-. 101 100 10 1 10 2. 0=-2.00e+00 0=-3.98e+00 0=-1.00e+03 0=-1.00e+09. 0 10 2 10 1 100 101 102 103. -102. -100. -100 -10. 2. x. 10. 2. 100. 100. 102. Figure 2.2: The velocity profile v(x). Red: σ0 = −2; Yellow: σ0 = −3.98; Blue: σ0 = −103 ; Green: σ0 = −109 . The dotted line is the sonic line. 22.

(33) Figure 2.3 shows the results of σ profile for several σ0 values. When x → −0, the value. of surface density σ is flat at the center. When x → ±∞, the value of surface density σ is. proportional to 1/x . When x → +0, the value of surface density σ is proportional to x−1/2 , and it becomes singular at x = 0.. x-. 105. 0=-2.00e+00 0=-3.98e+00 0=-1.00e+03 0=-1.00e+09. 103 101 10. 1. 10. 3. 10. 5. -102. -100. -100 -10. 2. x. 10. 2. 100. 100. 102. Figure 2.3: The density profile σ(x). Red: σ0 = −2; Yellow: σ0 = −3.98; Blue: σ0 = −103 ; Green: σ0 = −109 . The dotted line is the sonic line.. 2.6. Physical domain. In the subsection 2.5.1, we have found the solutions in the similarity domain. In this section, we would like to know how σ looks like in the physical coordinate of r and t, so we transform back to function σ of r and t. The Fig. 2.4 show the Σ profile for several moments of t with σ0 = 103 , and how a sheet-like structure collapses toward the center. Before time zero (t < 0), the central density increases until a central mass singularity forms at time 0. After time zero (t > 0), the central mass accretes the infalling gas, so the density decreases with time. 23.

(34) ( 0 = 103, t < 0). t=-1.00e+07 yr t=-3.16e+06 yr t=-1.00e+06 yr t=-3.16e+05 yr t=-1.00e+05 yr. (H2/cm2). 1025 1023 1021. ( 0 = 103, t > 0). t=1.00e+05 yr t=3.16e+05 yr t=1.00e+06 yr t=3.16e+06 yr t=1.00e+07 yr. 1025. 1019 1017 10. r. 1027. (H2/cm2). r. 1027. 1023 1021 1019. 10. 10. 8. 10. 6. 10. 4. 10. r (pc). 2. 100. 102. 1017 10. 104. (a) t < 0. 10. 10. 8. 10. 6. 10. 4. 10. r (pc). 2. 100. 102. 104. (b) t > 0. Figure 2.4: The Σ profile for several moments of t. We compare the above result to the observation results of Figure 2.5. We can see that reproduce the right order of magnitude for the column density, around 1021 at 1 pc.. (a). (b). Figure 2.5: (a) Herschel H2 column density (in cm−2 , Didelon et al. (2015)). The black polygon shows the area surveyed with the IRAM-30m telescope, while the white box corresponds to the inner 0.7 pc × 0.7 pc around the central hub and zoomed in the right panel. (b) The 13CO concentric annular surface density. The yellow dotted lines mark the different slopes in the surface density profiles. (Fig.1 and Fig.5 from Treviño-Morales et al. 2019). The Fig. 2.6 shows the u profile for several moments of t in σ0 = 103 . Before time zero (t < 0), the central velocity is small and the outer parts are almost constant. After time zero (t > 0), the velocity increases toward the center and the velocity increases with time. These is because the central mass is increasing. 24.

(35) r u ( 0 = 103,t < 0). t=-1.00e+07 yr t=-3.16e+06 yr t=-1.00e+06 yr t=-3.16e+05 yr t=-1.00e+05 yr. 103 102 1001 10. r u ( 0 = 103, t > 0). 104. u (m/s). u (m/s). 104. t=1.00e+05 yr t=3.16e+05 yr t=1.00e+06 yr t=3.16e+06 yr t=1.00e+07 yr. 103. 102. 2. 10. 1. 100. r (pc). 101. 101 10. 102. 2. 10. 1. (a) t < 0. 100. r (pc). 101. 102. (b) t > 0. Figure 2.6: The velocity profile for several moments of time.. We compare the above result to the observation results of Figure 2.7. We can see Fig. 2.7 that the velocity increases toward the central hub. It is consistent with the result of positive time in Fig. 2.6(b). The infall velocity is a few kilometers per second.. Figure 2.7: Position-velocity diagrams (color panels) along the ‘skeleton’ of filament F1 obtained from the 13CO (left) and C18 O (right) data cubes. The vertical yellow (dashed) lines indicate the transition between the hub and the filaments, corresponding to radii 20000 , 25000 (Rhub), and 30000 . The middle panels show the variation of velocity against the offset along the filament in two different manners: The dotted black line corresponds to the velocity obtained at the central pixel that constitute the skeleton of the filament, while the blue line shows the velocity along the skeleton after averaging over the velocity range shown in the top panels. The green lines indicate the velocity range where most of the emission of the filament resides. (Fig.B.3 from Treviño-Morales et al. 2019). 25.

(36) 26.

(37) Chapter 3. Perturbative analysis In this chapter, we would like to further understand how the filament structure of Mon R2 is formed. Therefore, we study the development of filaments in the sheet-like structure with perturbative analysis.. 3.1. Formulation. In this section, we would like to study how filaments grow inside a radially collapsing sheet. Therefore, we re-introduce the dependences in θ into our equations. We express the relation (2.1.10) in the continuity equation of Navier-Stokes equations (1.2.6) ∂Σ 1 ∂ 1 ∂ (rΣur ) + (Σuθ ) = 0. + ∂t r ∂r r ∂θ. (3.1.1). We express the relation (2.1.14)in the momentum equations of Navier-Stokes equations (1.2.7) and(1.2.8), the momentum equations become ∂ur ∂ur uθ ∂ur u2θ 2c2 ∂Σ ∂φ + ur + − =− − , ∂t ∂r r ∂θ r Σ ∂r ∂r. (3.1.2). ∂uθ ∂uθ uθ ∂uθ ur uθ 2c2 ∂Σ 1 ∂φ + ur + + =− − . ∂t ∂r r ∂θ r Σr ∂θ r ∂θ. (3.1.3). 2 2 ~ · ∇φ ~ = 1 ∂ (r ∂φ ) + 1 ∂ φ + ∂ φ = 4πGρ. ∇ r ∂r ∂r r2 ∂θ2 ∂z2. (3.1.4). The Poisson equation. We define the independent dimensionless similarity variable x, θ, τ and ζ as  r   x= ,    ct      θ = θ,            log t, t > 0.    τ = log |t| sgn(t) =        − log(−t), t < 0.      z    ζ= . ct 27. (3.1.5).

(38) Here τ represents a time relative to the stationary solution in the self-similar coordinate. We also define the dependent dimensionless similarity variables σ, vr , vθ and φ as Σ(x, θ, τ) =. c σ(x), 2πGt. ur (x, θ, τ) = cvr (x),. uθ (x, θ, τ) = cvθ (x),. φ(x, θ, τ) = c2 ψ(x). (3.1.6). The substitution of equations (3.1.5) and (3.1.6) into equations (3.1.1)-(3.1.4) yields the ordinary equations ∂vr vθ ∂σ σ ∂vθ σvr ∂σ ∂σ − (x − vr ) +σ + + −σ+ = 0, (3.1.7) ∂τ ∂x ∂x x ∂θ x ∂θ x ∂vr 2 ∂σ ∂ψ ∂vr vθ ∂vr v2θ Momentum equations: ± − (x − vr ) + − =− − , (3.1.8) ∂τ ∂x x ∂θ x σ ∂x ∂x ∂vθ ∂vθ vθ ∂vθ vr vθ 2 ∂σ 1 ∂ψ ± − (x − vr ) + + =− − , (3.1.9) ∂τ ∂x x ∂θ x xσ ∂θ x ∂θ 1 ∂ψ ∂2 ψ 1 ∂2 ψ ∂2 ψ Poisson equation: + + + = α = 2σδ(ζ). (3.1.10) x ∂x ∂x2 x2 ∂θ2 ∂ζ 2 Continuity equation: ±. We consider small perturbation and express them as    σ(x, θ, τ) = σb (x, θ, τ) + δσ(x, θ, τ),        vr (x, θ, τ) = vb (x, θ, τ) + δvr (x, θ, τ),      vθ (x, θ, τ) = δvθ (x, θ, τ),         ψ(x, θ, τ) = ψb (x, θ, τ) + δψ(x, θ, τ),       α(x, θ, τ) = αb (x, θ, τ) + δα(x, θ, τ).                               . (3.1.11). δσ(x, θ, τ) = δσ(x, θ)eωτ = δσ(x)eωτ cos(kθ), δvr (x, θ, τ) = δvr (x, θ)eωτ = δvr (x)eωτ cos(kθ), sin(kθ) , δvθ (x, θ, τ) = δvθ (x, θ)eωτ = δvθ (x)eωτ k δψ(x, θ, τ) = δψ(x, θ)eωτ = δψ(x)eωτ cos(kθ), δα(x, θ, τ) = δα(x, θ)eωτ = δα(x)eωτ cos(kθ).. The base solutions are simply the solutions that we previously found, which are stationary (∂/∂τ = 0) and axisymmetric ( ∂/∂θ = 0). They satisfy the equations: σb. ∂vb ∂σb ∂σb σb vb −x + vb − σb + = 0, ∂x ∂x ∂x x. (3.1.12). ∂vb ∂vb 2 ∂σb ∂ψb + vb + + = 0, ∂x ∂x σb ∂x ∂x. (3.1.13). 1 ∂ψb ∂2 ψb ∂2 ψb + + = αb = σb δ(ζ). x ∂x ∂x2 ∂ζ 2. (3.1.14). −x. 28.

(39) We rewrite (x − vb ) to vr0 , and rewrite formulas related to ψ with the following relationship ∂ψ(x) ∂δψ(x) ωτ , δγ x (x, θ) = − e cos(kθ), ∂x ∂x k 1 ∂ψ(x) , δγθ (x, θ) = δψ(x)eωτ sin(kθ). γθ (x) = − x ∂θ x. γ x (x) = −. (3.1.15). Substituting equations (3.1.11) into equations (3.1.7)-(3.1.10), and we use base solutions (3.1.12)(3.1.14) to simplify the perturbation equations to ! ! ∂δσ vr0 ∂vb σb ∂σb σb ∂δvr vr0 = ±ω − + δσ + + δvr + δvθ + σb , ∂x x ∂x x ∂x x ∂x ! ! ∂δvr 2 ∂σb ∂vb 2 ∂δσ vr0 =− 2 δσ + ±ω + δvr + − δγ x , ∂x ∂x σb ∂x σb ∂x 2k2 vb ∂δvθ =− δσ + ±ω + vr0 δvθ − kδγθ , ∂x xσb x ∂2 δψ 1 ∂δψ ∂2 δψ k2 + − δψ + = δα = 2δσδ(ζ). x ∂x ∂x2 x2 ∂ζ 2. (3.1.16) (3.1.17) (3.1.18) (3.1.19). Rearrange these three equations (3.1.16)-(3.1.18), we obtain the coupled set of ODEs ! v2r0 2 ∂σb ∂vb 2 ∂δσ 2 − vr0 = − ±ωvr0 − − + vr0 δσ ∂x x σb ∂x ∂x ! ∂σb ∂vb vr0 σb (3.1.20) − ±ωσb + + vr0 + σb δvr x ∂x ∂x vr0 σb δvθ + σb δγ x , − x ! 2ω 2vr0 vr0 ∂σb 2 ∂vb 2 ∂δvr 2 − vr0 =− ± − − + δσ ∂x σb xσb σ2b ∂x σb ∂x ! 2 ∂σb ∂vb 2 + vr0 − ±ωvr0 + + δvr x σb ∂x ∂x 2 − δvθ + vr0 δγ x , x vr0. ∂δvθ 2k2 vb =− δσ + ±ω + δvθ − kδγθ . ∂x xσb x. which should be solved together with the Poisson equation (3.1.19).. 29. (3.1.21). (3.1.22).

(40) 3.2. Treatment of self-gravity in a flat system. Similarly, we cannot calculate δψ(x) without knowing the complete profile of δσ(x). We therefore look for an approximation in order to do the first integration, and then try to find the exact solution by iteration. Once δσ(x) is known, δγ(x) can be calculated accordingly. We can then re-integrate Equations (3.1.20, 3.1.21 and 3.1.22) to obtain a new profile of δσ(x), which is closer to the real value. This process should be repeated until convergence is reached. The profile of δγ(x) can be calculated from δσ(x) with the Hankel transform. The following process is the use of the Bessel function to solve Poisson equation and the Hankel transform of these functions.. Figure 3.1: The iteration procedure for the perturbative solutions. 3.2.1. Approximation from local surface density. In this subsection, we solve the coupled set of ODEs using approximation from local surface density. We approximate the gravitational acceleration by using the local surface density profile. x δψ(x) ≈ − δσ(x), k. 3.2.2. δγ x (x) = −. ∂ψ(x) 1 ∂δσ(x) ≈ (δσ(x) + x ), ∂x k ∂x. k δγθ (x) = δψ(x) ≈ −δσ(x). x (3.2.1). Exact gravitational field from the surface density. The function Jk (sx)e−s|z| cos(kθ) satisfies the Poisson equation(1.2.10), for all values of s. We define δψ(x, θ, z) = δψ(x, z) cos(kθ),. (3.2.2). −s|z| ˆ δψ s (x, z) = sδψ(s)J , k (sx)e. (3.2.3). such that Z∞ δψ(x) =. δψ s (x, z = 0)ds 0. 30. (3.2.4).

(41) is a linear combination with all possible values of s. Correspondingly, δσ(x, θ) = δσ(x) cos(kθ),. (3.2.5). δσ s (x) = sδσ(s)J ˆ k (sx),. (3.2.6). such that Z∞ δσ(x) =. δσ s (x)ds.. (3.2.7). 0. We can verify that the above relations give us the Hankel transform pairs: Z∞. Z∞. Hk. sds δσ(s) ˆ Jk (sx) ⇐⇒ δσ(s) ˆ =. δσ(x) = 0. (3.2.8). xdx δψ(x)Jk (sx).. (3.2.9). 0. Z∞. Hk. Z∞. ˆ ˆ = sds δψ(s)J k (sx) ⇐⇒ δψ(s). δψ(x) =. xdx δσ(x) Jk (sx),. 0. 0. The Bessel function of order k, Jk , is the solution for Poisson equation of an infinitely thin sheet ∂2 δψ 1 ∂δψ ∂2 δψ k2 + − δψ + = δα = 2δσδ(ζ). x ∂x ∂x2 x2 ∂ζ 2. (3.2.10). Inserting δψ s cos(kθ) (Eq.3.2.2) and δσ s cos(kθ) (Eq.3.2.5) into the Poisson equation, we can find the relation between δψˆ and δσ. ˆ When |ζ| > 0, ! 2 s ∂Jk (sx) k2 2 ∂ Jk (sx) 2 ˆ +s δψ(s) − 2 Jk (sx) + s Jk (sx) = 0. x ∂(sx) ∂(sx)2 x. (3.2.11). When ζ = 0, ! 2 k2 s ∂Jk (sx) 2 ∂ Jk (sx) ˆ = 2δσ(s)J ˆ +s − 2 Jk (sx) − 2sJk (sx)δ(ζ) δψ(s) k (sx)δ(ζ). x ∂(sx) ∂(sx)2 x. (3.2.12). The first three terms between the parentheses are smaller than the fourth item, so we obtain ˆ −2sJk (sx)δψ(s)δ(ζ) = 2δσ(s)J ˆ k (sx)δ(ζ),. (3.2.13). such that ˆ ˆ = − δσ(s) δψ(s) = s. Z∞ xdx δψ(x)Jk (sx), 0. 31. (3.2.14).

(42) where δ(ζ) is the Dirac delta function. Finally, the gravitational acceleration for the infinity thin sheet then can be calculated as δγ x (x, θ) = −. ∂δψ(x, θ) ∂(δψ(x) cos(kθ)) =− ∂x ∂x. Z∞. = 0 Z∞. =. (3.2.15). s2 ˆ [Jk+1 (sx) − Jk−1 (sx)] cos(kθ) ds δψ(s) 2 s [Jk−1 (sx) − Jk+1 (sx)] cos(kθ) ds δσ(s) ˆ 2. 0. = δγ x (x) cos(kθ), 1 ∂δψ(x, θ) 1 ∂(δψ(x) cos(kθ)) =− x ∂θ x ∂θ Z∞ k ˆ = s2 ds δψ(s) Jk (sx) sin(kθ) x. δγθ (x, θ) = −. (3.2.16). 0. Z∞ =. s2 ˆ [Jk−1 (sx) + Jk+1 (sx)] sin(kθ) ds δψ(s) 2. 0. Z∞ =−. s [Jk−1 (sx) + Jk+1 (sx)] sin(kθ) ds δσ(s) ˆ 2. 0. = δγθ (x) sin(kθ). The following process is to convert σ into γ(x) through the calculation of Hankel transform. Hk. 1 −1 −1 2 (Hk−1 −Hk+1 ). Hk. −1 +H −1 ) − 12 (Hk−1 k+1. σ(x) −→ σ(s) ˆ σ(x) −→ σ(s) ˆ. −→. −→. 32. δγ x (x) δγθ (x). (3.2.17) (3.2.18).

(43) Chapter 4. Discussion and conclusions Throughout this work, we have looked at self-similar collapse of sheet-like structures and the filament growth in the molecular cloud. This thesis focuses on the understanding of the sheetlike structure formation: how it collapses, under which physical mechanisms it evolves, and how filaments grow inside it. In Chapter 1, we introduced the star formation in molecular clouds and the relevant physical mechanisms. We presented the observations of Mon R2. Chapter 2 focuses on the radial collapse of a thin sheet under self-gravity in the ISM. We derived the self-similar formulation from Navier-Stokes equations. And then, we solved the self-similar solution for radial collapse with numerical integration. We studied the behavior of the solutions for varying values of σ0 , and our results successfully describe the collapse of a sheet-like structure toward the center. Our model describes the global morphology and dynamics. This will allow to calculate the total mass, velocity, thickness and infer the inclination. We compared our results with the observations of Mon R2. They are indeed compatible. Furthermore, we studied the development of filaments in this structure with azimuthal perturbation in Chapter 3. We derived the perturbation formulation on top of the radial collapse self-similar solutions, and formulated the method of iteration for solving the self-gravity. In the future, we will solve this perturbation solution to understand the development of filaments. From the growth rate of the perturbative analysis, we will be able to find the most unstable mode k. This will give us the most probable number of filaments, and in turn, give us the line mass of each filament. Furthermore, we can study the fragmentation within these filaments to obtain the characteristic mass of prestellar cores. Two scenarios of cluster formation are well-known: nature vs nurture (Longmore et al. 2014). Our model describes a system that is forming stars while it is still globally collapsing. Stars form inside filaments, which are like conveyor belts (nurture). On the contrary, we cannot see that the gas first collapses to form a star-forming clump, and then the cluster forms in situ (nature). From the observation results, we see that the timescale of filament growth should be shorter than that of the global collapse. In the scenario described by our model, the massive 33.

(44) stars tend to form in the central hub. On the other hand, low-mass stars form in the filaments, and then fall into the central hub. This suggests a possible mechanism for primordial mass segregation of the cluster. To sum up, this study described the radial collapse of a thin sheet under self-gravity in the ISM, and solved the self-similar solution with numerical integration. Our results were consistent with the observations of Mon R2. In the future, our model can be further studied for the development of filaments in the hub-filament system.. 34.

(45) Bibliography André, P., Men’shchikov, A., Bontemps, S., et al. 2010, A&A, 518, L102 Bontemps, S., André, P., Könyves, V., et al. 2010, A&A, 518, L85 Contopoulos, I., Ciolek, G. E., & Königl, A. 1998, ApJ, 504, 247 Didelon, P., Motte, F., Tremblin, P., et al. 2015, A&A, 584, A4 Hanawa, T. & Matsumoto, T. 1999, ApJ, 521, 703 Hunter, C. 1977, ApJ, 218, 834 Kawachi, T. & Hanawa, T. 1998, PASJ, 50, 577 Krasnopolsky, R. & Königl, A. 2002, ApJ, 580, 987 Larson, R. B. 1969, MNRAS, 145, 271 Lasker, B. M., Sturch, C. R., McLean, B. J., et al. 1990, AJ, 99, 2019 Ledoux, P. 1951, Annales d’Astrophysique, 14, 438 Longmore, S. N., Kruijssen, J. M. D., Bastian, N., et al. 2014, in Protostars and Planets VI, ed. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning, 291 McKee, M. R. 1991, A&A, 251, 689 Penston, M. V. 1969, MNRAS, 144, 425 Schmitz, F. 1990, A&A, 230, 509 Schmitz, F. 1994, A&A, 284, 1007 Shi-Xue, D., Lan-Tian, Y., Xiao-Ci, L., & Xue-Bing, W. 2000, Ap&SS, 272, 383 Shu, F. H. 1977, ApJ, 214, 488 Simon, R. 1965, Annales d’Astrophysique, 28, 40 Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, 131, 1163 35.

(46) ´ et al. 2019, A&A, 629, A81 Treviño-Morales, S. P., Fuente, A., Sánchez-Monge, A., Whitworth, A. & Summers, D. 1985, MNRAS, 214, 1 Wu, X.-B., Li, Q.-B., Zhao, Y.-H., & Yang, L.-T. 1995, ApJ, 442, 736. 36.

(47)