Triggered star formation by shocks

Galactic star formation is classified two main processes. One is “spontaneous star formation,” in which the contraction of molecular clouds and subsequent star formation proceeds without any significant external disturbances. The other is “triggered star formation,” in which external factors promote the compression of molecular clouds, which induces star formation. The external compression is driven by shock waves generated by stellar winds, supernova explosions, cloud-cloud collision. For example, Elmegreen & Lada (1977) suggested that the compressed dense layer between shock and an ionization front can be gravitationally unstable. Cloud-cloud collisions have been proposed as an important mechanism in the formation of high-mass stars (e.g., Stolte et al., 2008; Furukawa et al., 2009; Wu et al., 2017). There is also observational evidence for star formation triggered by supernovae, ionization fronts, cloud-cloud and other shocks in the (e.g., Yokogawa et al., 2003; Ortega et al., 2004; Hester & Desch, 2005; Furukawa et al., 2009; Kinoshita et al., 2021).

, triggered star formation can occur when supersonic shocks sweep over clouds. The shock-cloud interaction is a fundamental and important physical phenomenon for star formation. , shock-cloud interactions are highly nonlinear hydrodynamic processes thus numerical simulation is tool to understand details. Over the past decades, the shock-cloud interaction has been explored by many groups (e.g., Klein et al., 1994; Xu & Stone, 1995; Boss, 1995; Nakamura et al., 2006; Pittard et al., 2009; Banda-Barragán et al., 2018). Most models assumed the two-dimensional (2D) axial-symmetric geometry. The first three-dimensional (3D) simulations were by Stone & Norman (1992). These simulations demonstrated that cloud destruction occurs faster in 3D the rapid growth of hydrodynamic instabilities in 3D. Later, the effects of various physical factors were explored by including turbulence (e.g.,Pittard et al., 2009), magnetic felds (e.g., Fragile et al., 2005), and radiative cooling (e.g., Mellema et al., 2002). incorporated self-gravity. These studies revealed the details of triggered star formation by shocks in the context of the formation of our System triggered by a supernova shock (e.g., Boss, 1995; Foster & Boss, 1996; Vanhala & Cameron, 1998; Vanhala & Boss, 2002 Boss et al., 2008; Leão et al., 2009; Boss & Keiser (2013); Li et al., 2014; Falle et al., 2017). These previous studies demonstrated that isothermal shocks can trigger gravitational collapse of stable clouds. Boss & Keiser (2013) that faster shocks destroy and disperse the cloud material before its collapse. , Falle et al. (2017) that slower shocks cannot induce collapse. These previous results indicate that only intermediate-speed shocks can trigger gravitational collapse.

Additionally, clouds and cores contain turbulent motions. Recently, Banda-Barragán et al. (2018) performed numerical experiments to investigate how cloud turbulence influences the shock-cloud evolution in the absence of self-gravity. They suggested that cloud turbulence cloud destruction and influences several ISM properties such as cloud porosity. Supersonic turbulence also enhances the acceleration of clouds to shocks.

In this , we 3D hydrodynamic numerical simulations to study shock-cloud . We an isothermal shock interacting with a Bonnor-Ebert sphere. In these simulations, self-gravity and sink particles included. This initial setup similar to those of Li et al. (2014) and Falle et al. (2017). Li et al. (2014) considered an application to the formation process of our Sun a shock-cloud interaction. Thus, their initial conditions very specific. Falle et al. (2017) investigated the early stages of shock-cloud and derived the condition for triggered collapse of stable Bonner-Ebert spheres. Here, we more general cases in a wider parameter range (from stable to unstable clouds). The inclusion of sink particles us to follow much longer evolution of the shock-cloud evolution. We also the effects of cloud turbulence on the shock-cloud interaction by including transonic cloud turbulence, which is often observed in the cloud cores in star-forming regions.

In the context of astronomical objects, we these simulations to represent the interaction of ISM shocks with dense cores and Bok globules in star-forming regions. Although the observed dense cores and Bok globules not in perfect equilibrium, some that the density structures of the dense cores and Bok globules in reasonable agreement with those of Bonner-Ebert spheres (e.g., Bacmann et al., 2000; Kandori et al., 2005; Alves et al., 2001). In some star-forming regions such as rho Oph and Orion A, the majority of the cores is likely to be pressure-confined (Maruta et al., 2010; Kirk et al., 2017). Therefore, our initial setup of the simulations is expected to represent reasonable conditions are in the ISM.

This paper is organized as follows. In Section 2, we describe the numerical methods, initial conditions, and simulation models which we in this . In Section 3, we present the numerical results. , we discuss interpretations of simulation results in Section 4. We derived a simple analytic condition for triggered gravitational collapse. Finally, we summarize our results and conclusion in Section 5.

In this , the simulations conducted using Enzo¹¹1http://enzo-project.org, adaptive mesh refinement (AMR) code (Bryan et al., 2014). We used Version 2.6 of the Enzo code in a 3D Cartesian coordinate system $(x,y,z)$. We numerically the following hydrodynamic equations for mass, momentum and energy conservation:

and

where $\rho$ is the density, $\bm{v}$ is the velocity, $p$ is the pressure, $\epsilon$ is the internal energy, $h=\epsilon+p/\rho$ is the enthalpy, $\phi$ is the gravitational potential. We the ideal gas law:

$\phi$ determined by solving the following Poisson’s equation

where $G$ is the gravitational constant $\rho_{\rm particle}$ is the density of sink particles assigned onto the finest grids by using a second-order cloud-in-cell interpolation technique (Hockney & Eastwood, 1988). See Section 2.2.5 for the details of the sink particles.

We a mean molecular weight $\mu$= 2.3, and $\gamma$ set to 1.00001 for an approximate isothermal assumption. In this , we the purely hydrodynamic problem ignoring radiative cooling, heating, magnetic fields, and thermal conduction.

The hydrodynamic equations solved a Runge Kutta second-order based (van Leer, 1977). The Riemann problem solved the , a two-wave, three-state solver with no resolution of contact waves, while the reconstruction method for the MUSCL solver . Bryan et al. (2014) for more details.

The maximum density is a good indicator of cloud stability. Figure 3 shows the maximum density normalized to the initial central cloud density as functions of time after a shock at cloud. Hereafter, in all figures, $t=0$ indicates the time when the shock front first the surface of the cloud. The vertical dashed line one free fall time $t_{\rm ff}=(3\pi/32G\langle\rho_{\rm cl}\rangle)^{1/2}$, where $\langle\rho_{\rm cl}\rangle$ is the mean density of the initial cloud. The solid lines the cases for which gravitational collapse is triggered by the shocks. The cases the clouds not collapse after shock passage are dashed lines. Below, we these two cases “triggered-collapse case” and “no-collapse case” respectively.

One important feature is that only intermediate shocks of $M_{\rm sh}=1.414.00$ can trigger cloud contraction. For no-collapse cases, the maximum density at the beginning but to lower values after rebounding. For example, for the model with a weak shock of $M_{\rm sh}=1.2$, the maximum density by one cloud free-fall time. However, it gradually with time. For the strong shock of $M_{\rm sh}=5.64$, the maximum density at $t=0.1$ Myr. However, the maximum density to . Even for $M_{\rm sh}=1.41$, the cloud once before collapsing at t = 0.43 Myr (point (1) in Figure 3). In triggered-collapse cases, the rate of density-increase becomes at the beginning. After some time, the rate of density-increase and the maximum density more than $\rho_{\rm th}$ and sink particles are formed.

As in Section 3.1, the density evolution depends on the Mach number $M_{\rm sh}$. Here, we present density distribution results for four cases: (1) a weak shock with $M_{\rm sh}$=1.41, in which the cloud slowly collapses, (2) an intermediate shock with $M_{\rm sh}$=3.15, in which the cloud rapidly collapses after shock passage, (3) a strong shock with $M_{\rm sh}$=5.0, in which the cloud does not collapse, and (4) a weakest shock with $M_{\rm sh}$=1.2, in which the cloud does not collapse.

The top panel in Figure 12 shows the time evolution of the mixing rate defined in Equation (12). The larger the shock Mach number is, the faster the mixing rate . For $M_{\rm sh}$ = 4.46$-$5.64, at $t=0.1$ Myr when the maximum density rebound , the mixing rate $0.2-0.3$, i.e. $20\%-30\%$ cloud gas mixed with the ambient gas. , for lower $M_{\rm sh}$ cases, the mixing rate not become higher than $0.2-0.3$ even when the cloud . A similar trend is indicated in the bottom panel in Figure 12, which shows the time evolution of living rate defined in Equation (14). For cases, rebounding, living rates $\sim$ 0.8, , for lower $M_{\rm sh}$ cases, living rates than $\sim$ 0.8 when maximum densities . During the cloud contraction, more gas removed in the larger $M_{\rm sh}$ cases than in the lower $M_{\rm sh}$ .

For triggered-collapse cases, the evolutions of mass and mass accretion rates of sink particles are shown in Figure 13. rates gradually by a few orders of magnitude and the mass of sink particles asymptotically. The top panel in Figure 13 indicates that the higher the Mach number , the lower the asymptotic sink particle mass . For $M_{\rm sh}$=1.41, 1.99, 3.15, and 4.00, periods when mass accretion rates above $10^{-6}M_{\odot}\rm yr^{-1}$, 0.22, 0.18, 0.13, and 0.07 Myr, respectively. The higher the Mach number , the shorter the time with high accretion rate . This trend of accretion timescale would affect the asymptotic sink particle mass.

Here, we results of turbulent cloud models. We followed models of turbulent ($M_{\rm tub}=1.00$) clouds with a critical Bonnor-Ebert radius by changing the shock Mach number of $1.41-5.64$. Figure 15 shows the time evolution of $\rho_{\rm max}/\rho_{\rm c}$, $M_{\rm mix}/M_{\rm cl}$, $M_{\rm live}/M_{\rm cl}$, $M_{\rm sink}/M_{\rm cl}$, and accretion rates of sink particles. Figure 15 (a) indicates that $M_{\rm sh}=1.41-3.15$ shocks cloud collapse, relatively stronger $M_{\rm sh}=4.46$ and $5.64$ shocks not induce cloud collapse. Figure 16 and 17 show the slices of the mass density distribution for $M_{\rm sh}$=3.15 and 4.46 as Figure 4. $M_{\rm sh}$=3.15, a strongly compressed layer formed and then the cloud and . , $M_{\rm sh}$=4.46, a strongly compressed layer formed initially the cloud later destroyed by the shock and mixed with the ambient gas. Figure 15 (b) and (c) indicate that the larger the Mach number , the shorter the mixing timescale with the ambient gas , higher living rates. Similar to the non-turbulent models, intermediate shocks cloud collapse, strong shocks destroyed clouds before collapse and mixed clouds with ambient gas faster. Figure 15 (d) and (e) show the evolution of mass and mass accretion rates of sink particles. Cross points in Figure 14 shows results of turbulent cloud models. The larger the Mach number, the shorter the span of effective accretion rate, and the lower the asymptotic mass of the sink particle.

To investigate the effect of cloud turbulence on shock-cloud evolution, here, we will compare turbulent models with corresponding non-turbulent models for the same cloud radius and shock Mach number. Figure 18 shows time of some physical quantities in both turbulent and non-turbulent models for $\xi=6.45$ and $M_{\rm sh}=3.15$. Figure 18 (a) maximum density evolution $\rho_{\rm max}/\rho_{\rm c}$. Up to $t\sim 0.2$ Myr density evolution in two cases similar. , the turbulent cloud a slower increase in density. As shown in Appendix D, for the turbulent model, the sink particle was formed by 0.49 Myr after the shock touched the cloud, for the non-turbulent clouds model, it was formed by 0.31 Myr. That is, the turbulent cloud a slower increase in density and slower sink particle formation than the non-turbulent cloud. Figure 18 (b) and (c) show the evolution of mixing rates and living rates. For the turbulent cloud, the living rate faster and mixing rate faster than that of non-turbulent . That is, the turbulent cloud faster with the surrounding ambient gas and destroyed. Figure 18 (d) and (e) show the evolution of sink particle mass and accretion rates. The asymptotic sink particle mass in the turbulent cloud model lower than non-turbulent counterparts. The accretion rate $10^{-5}M_{\odot}$ by $t=0.06\rm Myr$ the turbulent cloud, this occurs at $t=0.18$ Myr for the non-turbulent . For the turbulent cloud, accretion time shorter and the asymptotic mass of the sink particle lower.

As above, the turbulence cloud contraction, the destruction of the cloud, and the mass of the formed stars. We can that turbulence in the dense cloud has the effect of suppressing star formation by shocks.

As shown in Figure 14, shocks that are strong or weak cannot induce cloud collapse. This is because strong shocks tend to destroy the cloud hydrodynamic instabilities and weak shocks do not compress the cloud to become unstable. This implies the existence of parameter windows for cloud collapse. In this subsection, we the conditions for gravitational collapse triggered by shocks.

Figure 19 shows the for the $\xi=3.22$ models as Figure 5. The vertical axis shows the value of the velocity normalized by estimated propagating shock velocity at the cloud surface $v_{\rm sur}=v_{\rm sh}\chi_{\rm s}^{-1/2}$ (see Equation (B1)). with the color variable $C>0.99$ in green, while those of $0.1<C<0.99$ were shown in red. The panels in the upper row are for the $M_{\rm sh}=1.20$ cases corresponding to $t=0.27t_{\rm cc}$ and the rebounding phase. For $M_{\rm sh}=1.20$, gas accounting for more than 50$\%$ of the mass compressed almost without being stripped. However, as mentioned in Section 3.2.4, the clouds not sufficiently compressed to develop gravitational instability. The panels in the middle row are for the $M_{\rm sh}=3.15$ cases corresponding to $t=0.27t_{\rm cc}$ and before sink particle creation. In panel (c), the distribution is shaped like a horizontal V with the root at $v/v_{\rm sur}\sim 1$. The gas distributed from the root of V toward the origin of the figure corresponds to the gas compressed by the shock progressing to the center of the cloud. Gas distributed from the base of V to the upper left of the figure corresponds to gas stripped and dissipated by the shock. For $M_{\rm sh}=3.15$, compressed gas dense, gravitational collapse. The panels in the lower row are for the $M_{\rm sh}=5.00$ cases corresponding to $t=0.27t_{\rm cc}$ and the rebounding phase. In panel (e), the distribution is shaped like a V as in panel (c), but the V is wider and more gas is dissipated. As shown in Section 3.2.4, eventually, the entire cloud destroyed.

these results, two conditions must be for the cloud to contract. The first condition the degree of compression by shocks. For low , most of the gas does not become dense and gravitational collapse cannot be induced. It is assumed that for cloud collapse, compression to make Bonnor-Ebert spheres unstable is required. This condition would determine the lower limit of the Mach number. The second condition the destruction of clouds. Even if clouds are strongly , shocks with Mach numbers that the destruction of clouds. to collapse, the timescale for the development of gravitational instability must be shorter than that of destruction. This condition would determine the upper limit of the Mach number. Hereafter, based on the above two conditions, we will estimate the Mach number parameter window for the cloud collapse.

First, we consider the lower limit of $M_{\rm sh}$. Considering an isothermal shock with a Mach number $M_{\rm sh}$, the ambient gas pressure behind the shock is $M_{\rm sh}^{2}$.

Considering conditions for a Bonnor-Ebert sphere to become unstable, we expect

where $P_{\rm crit}$ is the critical pressure of Bonnor-Ebert sphere (see Equation (A9)). This relation implies

Substituting Equation (A2),(A4) and (A6), we a dimensionless expression:

Since the side of Equation (18) is a function of $\xi$, it can also be expressed a function of $\rho_{\rm c}/\rho$. the third degree polynomial regression of the side of Equation (18) in the range of $2.0\leq\rho_{\rm c}/\rho\leq 14.1$, we the approximate equation of lower limit of $M_{\rm sh}$

$\rho_{\rm c}/\rho=14.1$ corresponds to $\xi=\xi_{\rm crit}$. When $\rho_{\rm c}/\rho>14.1$, a Bonnor-Ebert sphere is unstable.

Next, we consider upper limit of $M_{\rm sh}$. Iwasaki & Tsuribe (2008) that the timescale of the gravitational instability of isothermal layers bounded by shock wave with the Mach number $M_{\rm sh}$ is an order of $t_{\rm ff}/\sqrt{M_{\rm sh}}$ where $t_{\rm ff}$ is the free-fall time scale of the preshock region. In our model, the sum of the timescale on which the shock propagates through the cloud and the timescale on which gravitational instability behind the shock is estimated

where $t_{\rm ff}$ is the free fall time scale of spherical cloud.

It is known that the timescale of the cloud destruction is the order of the cloud crushing time $t_{\rm des}$=$\alpha t_{\rm cc}$ (see Appendix B). Previous numerical hydrodynamical simulations that $\alpha\sim 2.5-4$ (e.g., Klein et al., 1994; Poludnenko et al., 2002). by shocks, gravitational instability must within the timescale of cloud destruction $t_{\rm des}$. Hence, we expect

This implies

This equation can be rewritten dimensionlessly as

the the minus third degree polynomial regression of the side of Equation (23) in the range of $2\leq\rho_{\rm c}/\rho\leq 100$, we the approximate equation of upper limit of $M_{\rm sh}$

Figure 20 shows the pressure ratio $(P_{\rm 0}/P_{\rm crit})$ vs $(M_{\rm sh})$ parameter space for an initially stable Bonnor-Ebert sphere. The solid line demarcates the Mach number estimate based on Equation (18) above which the cloud will become unstable by the shock. The light gray shaded area shows the upper limit estimate for Mach number based on Equation (23) applying $2.25<\alpha<2.50$ above which the cloud will be destroyed by hydrodynamic instabilities before the collapse. In the dark gray region, clouds are expected to be induced to collapse by shocks. Based on simulation results, red points indicate initial cloud parameter pairs corresponding to eventual collapse, while black points indicate cloud parameter pairs corresponding to predicted non-collapse. that red points are in dark or light gray regions and black points are outside of the light gray region. Hence, the simulation results are in agreement with our estimates of conditions. However, the upper limit estimate for Mach numbers has small range in terms of $\alpha$. This upper limit is derived based on estimate, and hydro instability is a complex process with uncertainty. It may be difficult to set an upper limit with a single line. the range of $\alpha$ used is not large. The criteria we have derived are useful for estimates of stable cloud collapse.

Figure 21 shows the density ratio $(\rho_{\rm c}/\rho_{\rm 0})$ vs $(M_{\rm sh})$ parameter space similar to Figure 20 with all parameter pairs shown in Figure 14. For initially unstable clouds (i.e. in the right-side region of the critical density line), when applying smaller $\alpha$ ($\sim$ 2.0), results of initially unstable clouds consist with the upper limit estimate. The reason for this gap in $\alpha$ range between stable and unstable clouds is that for unstable clouds, the propagating shock speed estimated Equation (B1) and actual one an exact match. For unstable clouds, the density ratio $\rho_{\rm c}/\rho_{\rm 0}$ is larger and a density gradient along the radial direction affects the shock velocity propagating in the cloud. Therefore, the timescale for shock propagation and cloud contraction expressed by the side of Equation (21) is affected, and the $\alpha$ range changes. However, our upper limit estimate is useful for an order discussion. Although the upper limit has some uncertainty for larger unstable clouds, the general tendency of results is consistent with our estimate.

In our calculation, for clouds initially at $\xi>6.45$, at the time when propagating shock to compress clouds, little gravitational collapse is in progress. That is, the density profile of the clouds is almost identical between the initial conditions and when the shocks arrive. We note that $\xi>6.45$ unstable clouds in which collapse has progressed further by the time the shock arrives, shock-cloud evolution may change and ranges of $M_{\rm upp}$ from our results.

As shown in Section 3.4, comparing results with the same initial radii, the higher Mach number of the shock is, the lower asymptotic sink particle mass and the shorter accretion time is. This would be because the higher the Mach number is, the faster the shock the cloud around sink particles. As shown in Figure 6, the cloud around the sink particle is stripped and mixed with ambient gas. Cloud material is accelerated and displacement of cloud material and the sink particle . , the accretion timescale would be determined by how the incoming shock can accelerate and strip accreting clouds from sink particles.

Here, we will the moving of stripped clouds and sink particles. We use the $\langle z\rangle$ defined in Equation (15) to analyze the moving of stripped clouds.

Figure 22 shows the $\langle z\rangle$ and position of sink particles in the z-direction for $\xi=3.22$ and $M_{\rm sh}=3.15$ $4.00$. Displacement of the cloud and sink particle is initially small because of the low mass of sink particles. , the sink particle drifts behind and displacement becomes larger because accretion progresses and the sink particle mass .

Figure 23 shows the relative and velocity $\langle z\rangle$ and those of sink particles for $\xi=3.22$ and $M_{\rm sh}=3.15$ $4.00$. The higher the Mach number of the shock is, the faster the relative and velocity . This trend is also for other radii (see Appendix F). Therefore, relative displacement and velocity of the particle and cloud determined by the shock speed and the timescale of accretion and asymptotic mass.

Figure 24 shows the as Figure 19 for both turbulent and non-turbulent ($M_{\rm sh}=3.15$ and $\xi=6.45$). Figure 24 (a) and (d) are at the same time. Figure 24 (b) and (e) also at the same time, and Figure 24 (e) corresponds to the before sink particle. Comparing Figure 24 (a) and (d), cloud mass less on the high density side compared . As in Section 3.5, the turbulent cloud has a slower increase in density and a slower sink particle formation than the non-turbulent cloud. This slowness would be due to the effective pressure from turbulence in the cloud. Internal turbulence increases effective cloud pressure to the order of $\sim\rho\sigma^{2}$. This increased pressure enhances cloud diffusivity and suppresses cloud contraction by the external ram pressure due to shock. Therefore, turbulence prevents rapid gravitational collapse and high density region is formed slowly.

Figure 24, throughout the evolution the turbulent cloud, the distribution of the gas shown in red ($C<0.99$) from that of the non-turbulent cloud. For example, in Figure 24 (d), red and green gas shows near V-shaped distribution, while in Figure 24 (a), red gas shows a fan-like distribution. In other words, more gas is dispersed and mixed with the ambient gas. This trend is consistent with that of Figure 18 (b) and (c). Clouds with a turbulent velocity are more prone to Kelvin-Helmholtz instabilities than non-turbulent counterparts. Therefore, that turbulent vortical motions in the clouds diffuse the cloud, and the timescale of the mixture becomes shorter.

Figure 23 shows the relative and velocity $\langle z\rangle$ and those of sink particles for both turbulent and non-turbulent ($M_{\rm sh}=3.15$ and $\xi=6.45$). For the turbulent case, throughout the entire evolution, the relative is greater than that of the non-turbulent case. Most of the time, the relative velocity is also greater for the turbulent case. One reason for this trend is that the cloud with a turbulent velocity is destroyed by the shocks and the gas around the sink particles is stripped faster downstream of the shock. Another reason is that, for the turbulent case, the sink particle formed later, and by the time mass accretion , more gas been accelerated. , the asymptotic sink particle mass in the turbulent clouds models is lower than the non-turbulent counterparts.

Comparing No.5 and No.8 models results in Appendix D, the upper limit of the Mach number for cloud collapse . In turbulent cloud cases, the upper limit of the Mach number is $3.15<M_{\rm upp}<4.46$ for the not turbulent cases, it is $5.64<M_{\rm upp}<6.00$. , parameter window for the cloud collapse narrower for turbulent cloud cases.

, turbulence in the cloud makes triggered star formation by shocks more difficult. Pressure due to turbulence cloud contraction and facilitates mixing with the ambient gases. The asymptotic stellar mass . that our simulation models address magnetic fields. In a realistic ISM, magnetic fields exist and can alter the physical process of the shocked cloud. In this , we purely hydrodynamic cases and investigate the effects of turbulence.