STAR FORMATION IN TURBULENT MOLECULAR CLOUDS WITH COLLIDING FLOW

The evolution of the high-density model with ${\cal M}_t = 10$ and ${\cal M}_f =0$ is described here as a fiducial model without the colliding flow. Figure 1 shows the evolution of model HT10F0 by the distributions of column density and sink particles. In the early phase, turbulence produces many shock waves, as shown in Figure 1(a). The density, which increases by a factor of ∼10² due to shock compression, indicates that the density contrast of the isothermal shock is equal to the square of the Mach number. The first star forms at t = 5.65 × 10⁵ yr, which corresponds 0.22 t_cross and 0.62 t_ff, as shown in Figure 1(b). Sink particle formation earlier than the freefall time is the result of the density increase due to the shock waves. At this stage, models HT10F0 and HT10F0wog show similar density distributions, which indicates that turbulence mainly controls the structure formation. Figure 1(c) shows the cloud at the stage where the maximum mass of sink particles reaches 10 M_☉ (hereafter referred to as "the stage of M_{⋆, max} = 10 M_☉"). By this stage, the density contrast increases more than that shown in Figure 1(b), and it is slightly higher than that of model HT10F0wog because of self-gravity. Sink particles are created in clusters along the dense filaments.

Zoom In Zoom Out Reset image size

Figure 2 shows the dependence of the high-density models (upper panels) and low-density models (lower panels) without colliding flow (${\cal M}_f=0$) on the initial strengths of the turbulence.

Zoom In Zoom Out Reset image size

The stronger turbulent models produce a fine structure of gas, whereas the weak turbulent models show parsec-scale filaments irrespective of the initial densities. Model HT30F0 exhibits small structures disturbed by strong turbulence, whereas model HT3F0 exhibits a network with long filaments. Similarly, model LT3F0 shows smoother filaments than does LT30F0.

For the high-density models, a model with stronger turbulence produces more sink particles in a shorter time (see also Table 1). This is attributed to strong compression due to a large amount of turbulence, as shown in the previous work of Klessen et al. (2000). The long loci of sink particles in the more turbulent model, HT30F0, indicate the high speeds of the sink particles. These sink particles move at velocities roughly similar to those of the gas where they formed. Details of the velocity distribution for sink particles are shown in Section 4.9.

The low-density models produce much smaller numbers of sink particles compared to the high-density models (see also Table 1). This indicates that the high-density models contain larger masses than the low-density models. Contrary to the high-density models, a model with stronger turbulence produces fewer sink particles by the stage of M_{⋆, max} = 10 M_☉, as shown in the lower panels of Figure 2 and Table 1. However, at a given time, a model with stronger turbulence produces more sink particles, showing that the turbulence promotes star formation even for the low-density models.

Model HT10F10 is the same as the previous model, HT10F0, but colliding flow with a Mach number of 10 is added. Figure 3 shows the evolution of model HT10F10. The upper and lower panels show the column densities along the x- and z-directions, respectively. Figures 3(a) and (d) show the cloud in the early phase. The column density along the x-direction (Figure 3(a)) is similar to that of model HT10F0 (Figure 1(a)), and the gas compression due to the colliding flow is shown in the column density along the z-direction (Figure 3(d)). Figures 3(b) and (e) show the cloud when the first sink particle forms. The colliding flow accumulates filaments to form a sheet at x ≃ 0, and the dense gas of the sheet leads to star formation at an earlier time than for the previous model.

Zoom In Zoom Out Reset image size

Figures 3(c) and (f) show the cloud at the stage of M_{⋆, max} = 10 M_☉. The thickness of the sheet increases because the accumulated filaments penetrate the sheet. In the comparison of Figure 1(c) and Figure 3(c), model HT10F10 exhibits more prominent fine filaments. The filamentary structure is enhanced by the self-gravity, while the structure formation is primarily dominated by the turbulence. The colliding flow increases the density in the sheet, increasing the effects of the self-gravity. We confirmed that the dense filaments with widths of ∼0.1 pc have the Jeans length less than 0.1 pc, indicating that the self-gravity also contributes to the structure formation mainly on a scale of the filament width. The number of sink particles is roughly the same as that of model HT10F0 (Table 1), but the elapsed time is shorter than that in model HT10F0 due to the high accretion rate, which indicates that the colliding flow promotes star formation.

The column density along the z-direction shows that dense filaments lie parallel to the sheet and low-density striations lie perpendicular to the sheet. The low-density striations are oriented parallel to the accretion flow on the sheet. The column density and the velocity are in agreement with recent observations of Palmeirim et al. (2013). After finding striations perpendicular to the filament in the B211/L1495 region in the Taurus molecular cloud, they suggest that gas may be accreting along the striations onto the main filament.

Figure 4 shows the high-density models with different strengths of colliding flows in column densities along the x-direction (upper panels) and the z-direction (lower panels). The model with a strong flow, HT10F30, exhibits strong compression of gas into the plane at x ≃ 0, as shown in the column densities along the z-direction. This leads to fine filaments and active formation of sink particles compared to the previous model, HT10F10, as shown in the column densities along the x-direction. Model HT10F10 produces 142 sink particles, while model HT10F30 produces 915 sink particles.

Zoom In Zoom Out Reset image size

The model with weak colliding flow, HT10F3, does not exhibit a clear sheet structure at x ≃ 0, as shown in Figure 4(d). This is because the amplitude of the colliding flow is considerably lower than the mean velocity of the initial turbulence (${\cal M}_f < {\cal M}_t$). The low compression due to the colliding flow leads to a long elapsed time for the maximum mass of the sink particles to reach 10 M_☉. This long time brings about formation of many sink particles compared to the case for model HT10F10.

Figure 5 compares the low-density models with different strengths of colliding flows. In the low-density models, formation of sink particles begins after the colliding flows pass by each other. At the stage of M_{⋆, max} = 10 M_☉, the sheet structures due to colliding flows disappear.

Zoom In Zoom Out Reset image size

Figure 6 shows the decay of turbulence by measuring the standard deviations of the velocities (velocity dispersion), $\sigma _v = \langle ({\boldsymbol {v}}-\langle {\boldsymbol {v}}\rangle _m)^2 \rangle _m^{1/2}$, where 〈 · 〉_m denotes the mass-weighted average over the computational domain, i.e., $\langle {\boldsymbol {v}} \rangle _m = \int _V \rho {\boldsymbol {v}} \, dV/M_\mathrm{gas,tot}$ and M_{gas, tot} = ∫_VρdV. The total kinetic energy of gas is therefore evaluated as $E_K = M_\mathrm{gas,tot} \sigma _v^2/2$ when $\langle {\boldsymbol {v}}\rangle _m=0$. Our numerical scheme conserves the linear momentum of gas within a truncation error before sink particle formation, and the mean velocity $\langle {\boldsymbol {v}}\rangle _m$ therefore has small values compared to σ_v. Even after sink particle formation, the mean velocity remains small, $\langle {\boldsymbol {v}}\rangle _m / \sigma _v \lesssim 10^{-3}{\rm to}10^{-2}$ at most.

Zoom In Zoom Out Reset image size

For models without colliding flows (Figure 6(a)), models with a common ${\cal M}_t$ exhibit similar decay of turbulence irrespective of their initial densities. After the stages of formation of the first sink particle, as denoted by circles, the high-density models tend to have a slightly higher σ_v than the low-density models because of self-gravity. The models without colliding flows produce the first sink particle roughly around ∼ t_ff irrespective of the initial density, whereas the time of formation of the first sink particle tends to be early when assuming strong turbulence.

When the time and velocity dispersion are normalized by the crossing time and the effective Mach number, respectively, the evolutions for the different models converge to that for model MT10F0wog (non-gravitational model), as shown in Figure 7(a). This indicates that the turbulence decays in the crossing timescale, t_cross, are in agreement with Mac Low et al. (1998); Mac Low (1999), and Ostriker et al. (2001). The high-density models and the low-density model LT3F0 undergo gravitational collapse to form sink particles before the initially assumed turbulence decays considerably.

Zoom In Zoom Out Reset image size

Figure 6(b) and Figure 7(b) show velocity dispersions as functions of time for the models with colliding flows. The colliding flows contribute to σ_v, and σ_v at t = 0 depends on ${\cal M}_f$ when comparing models with the same ${\cal M}_t$, as shown in Figure 6(b). All the models tend to converge to the locus of the non-gravitational model as time proceeds (Figure 6(b)). This tendency indicates that colliding flows decay faster than do isotropic turbulent flows. For example, models with ${\cal M}_f = 30$ show rapid decreases in σ_v, and the timescales of the decreases are smaller than t_cross (Figure 7(b)), where σ_v decreases by a factor of 50% in a time of 0.3 t_cross. The rapid decreases in σ_v are responsible for the shock waves formed by the colliding flows. The colliding flow brings forward the first sink particle formation, and the elapsed time until the first sink particle formation is shorter than the freefall time irrespective of the initial densities.

A velocity field can be separated into solenoidal (transverse) and compressive (longitudinal) components by applying a Helmholtz decomposition. The decomposition of the velocity field has also been used to discuss turbulence (e.g., Kitsionas et al. 2009; Federrath et al. 2010b). We performed the velocity decomposition and estimated E_trans and E_long according to Equations (B8) and (B9), where E_trans and E_long are the powers of the transverse and longitudinal velocity components, respectively (see Appendix B). Figure 8 shows E_long/E_tot as functions of time, where E_tot = E_long + E_trans.

Zoom In Zoom Out Reset image size

The models without colliding flow show E_long/E_tot = 0 at the initial stages (t = 0) because the initial velocity field is incompressible. As time proceeds, E_long/E_tot increases. This increase indicates that the compressibility of the velocity field increases. In the early stages, models with the same turbulent Mach number and a different initial density exhibit similar increases in E_long/E_tot. Models HT10F0 and LT10F0 show increases in E_long/E_tot at a rate similar to the non-gravitational model HT10F0wog. Models with strong turbulence tend to show a rapid increase in E_long/E_tot, and the increase rates are in proportion to the crossing time, t_cross. These increases indicate that evolution in the early stages is controlled mainly by the turbulence.

In the later stages, the self-gravitational models HT10F0 and LT10F0 exhibit higher ratios of power than does the non-gravitational model HT10F0wog. This indicates that self-gravity increases the compressibility of the velocity field. The high-density models tend to exhibit more rapid increases in E_long/E_tot than do the low-density models, and this result underscores the differences in freefall times between them.

The non-gravitational model HT10F0wog, followed for a long time of 2.5 × 10⁷ yr, exhibits considerable oscillation of E_long/E_tot in the range of 0.15–0.2. After the velocity becomes subsonic σ_v ≲ c_s, the oscillation dumps and E_long/E_tot ≃ 0.15 is shown.

Models with colliding flows exhibit E_long/E_tot > 0 at the initial stages because the colliding flows are compressible and contribute to E_long. The high-density models tend to have a high ratio of E_long/E_tot than do the low-density models. Moreover, in comparison with the non-gravitational models, the self-gravitational modes have a high ratio. These comparisons also indicate that self-gravity enhances the compressibility of the velocity field.

The non-gravitational model HT10F10wog exhibits an oscillation of E_long/E_tot with decaying amplitude in the range of 0.2–0.4 when σ_v ≳ c_s (the early part of the epoch is shown in Figure 8), and it decreases to ≃ 0.15 when σ_v ≲ c_s. The amplitude of the oscillation is larger than that for the non-colliding model HT10F0wog. The oscillation is responsible for the oscillation on E_long, while E_trans decreases monotonically.

It is known that the probability distribution function (PDF) of density in isothermal turbulent gas without self-gravity is well approximated by a lognormal distribution (e.g., Vazquez-Semadeni 1994). The lognormal distribution is expressed as (Ostriker et al. 2001; Padoan & Nordlund 2002)

$$\begin{equation} p(s) ds = \frac{1}{ \left(2\pi \sigma _s^2\right)^{1/2} } \exp \left[ -\frac{1}{2} \left(\frac{s-\bar{s}}{\sigma _s} \right)^2 \right] ds, \end{equation} \tag{ 3 }$$

$$\begin{equation} s = \ln \left(\frac{\rho }{\rho _0}\right), \end{equation} \tag{ 4 }$$

$$\begin{equation} \bar{s} = -\frac{\sigma _s^2}{2}, \end{equation} \tag{ 5 }$$

where $\bar{s}$ and σ_s denote the average and the standard deviation of the logarithmic density, s, respectively. The standard deviation, σ_s, is found to be a function of the Mach number for the gas, ${\cal M}$, as (Padoan et al. 1997)

$$\begin{equation} \sigma _s^2 = \ln (1+ \gamma ^2{\cal M}^2), \end{equation} \tag{ 6 }$$

where γ ≃ 0.5 and ${\cal M} = \langle \left(v- \langle v\rangle \right)^2 \rangle ^{1/2}/c_s$, and 〈 · 〉 denotes a volume-weighted average.

At the initial stage, the PDFs are expressed by a delta function because the initial density is constant. As time proceeds, the profiles of the PDFs are broadened, and they are well expressed by the lognormal distributions for all models without colliding flow in the early phase, typically before sink particles form. Figure 9 shows the volume-weighted PDFs for models without colliding flow at the stages of M_{⋆, max} = 10 M_☉. For models with ${\cal M}_t = 3$, the PDFs have power-law tails in the high densities and so deviate from the lognormal functions. This indicates that the cloud structure in high density is affected by self-gravity. A deviation of PDFs from the lognormal distribution was also reported by Klessen (2000); Federrath et al. (2008), and Federrath & Klessen (2013) for self-gravitational isothermal turbulence. During the evolution, the profiles of the low-density tails oscillate around the lognormal distributions. The PDFs of the snapshots shown in Figure 9 therefore deviate from the lognormal distribution at low densities.

Zoom In Zoom Out Reset image size

In Figure 9, two fitting curves are plotted for each model based on the lognormal functions. The standard deviation of the logarithmic density is derived from the density distribution (green curve), and the standard deviation is derived by using Equation (6) (blue curve). In the case of non-gravitational driven turbulence, Konstandin et al. (2012) reported γ = 1/3 for solenoidal forcing and γ = 1 for compressive forcing (see Price et al. 2011). We adopt γ = 1/2 here because it provides a better fit to the data than do the other values. These two lognormal fitting curves roughly coincide and indicate that Equation (6) is satisfied well in the stages examined here, even in the case of self-gravitating gas.

Figure 10 shows the volume-weighted PDFs for models with colliding flow. In the cases with strong colliding flow, especially in the high-density models, the PDFs deviate considerably from the lognormal distributions around the peaks. This deviation is attributed to the sheet structures of the clouds caused by the colliding flow, as shown in Figures 4(e) and (f). Moreover, the two fitting curves, which show significantly different distributions, indicate that Equation (6) is no longer satisfied in the cases with strong colliding flow.

Zoom In Zoom Out Reset image size

PDFs for the column density are shown because they are often examined in observational studies, (e.g., Lombardi et al. 2011). Ostriker et al. (2001) showed that the PDF for the column density can also be fitted by lognormal functions.

Figure 11(a) shows the PDFs for the column density for model LT10F0. The density PDF for Model LT10F0 has a little excess from the lognormal function in the high densities (the bottom panel in the middle of Figure 9). The profiles of the column density PDFs are bumpy compared to the density PDFs because the sampling numbers for the column density are less than those for the density. For the first three stages shown in Figure 11(a), each PDF can be fitted well by a lognormal function (dotted curves). To examine the long-time evolution of the PDF for this model, we follow the cloud evolution until the stage of M_{⋆, max} = 50 M_☉ (t = 3.4 Myr). The PDF at this stage has a wing at high column density, which can be fitted by the power law of PDF∝Σ⁻². The estimated power indexes depend on the models; the power index is approximated as −3 for models HT3F0 and HT10F0, and seems to change from −2 to −3 as the column density increases for model LT3F0. Such a power-law wing in PDFs has often been observed in several star-forming regions (i.e., Kainulainen et al. 2009), as will be discussed in Section 5.1.1.

Zoom In Zoom Out Reset image size

Figure 11(b) shows the column density PDFs for the colliding model HT10F10. The line of sight is the z-direction, which is perpendicular to the colliding flow. Along this line of sight, the sheet cloud appears in the edge-on view. The PDFs highly deviate from lognormal functions and exhibit two peaks, which are responsible for the sheet cloud with a high column density and the colliding flow with a low column density. We confirmed that, when the line of sight is taken to be the face-on direction, the PDFs can be fitted by a lognormal function. In summary, the features of the PDFs for the column density depend on the line of sight, and they can highly deviate from a lognormal function when the cloud is disturbed by a large-scale flow.

Figure 12 shows the mass distribution for sink particles in the histograms for all models. The high-density models exhibit clear features in the histograms, whereas the low-density models do not because of the small number of sink particles produced. For the high-density models, the mass distribution has a peak around ∼1 M_☉ and has tails in both the high masses and low masses. For all the high-density models except model HT10F30, the tails can be roughly fitted by the relationship $N_\mathrm{star}(M_\mathrm{star}) \propto M_\mathrm{star}^{-1.35}$, which is a classical initial mass function (IMF) of Salpeter (1955). Similar tails at high masses were reported for a cloud with a different mass of 500 M_☉ by Bate (2009) and Bate (2012). Model HT10F30 shows considerably steeper tails than does the fitting line at high mass. This is because this model produces 3.6–6.5 times more sink particles than do the other models (models HT10F3 and HT10F10) and the masses of the peaks and the maximum masses are roughly common among the models.

Zoom In Zoom Out Reset image size

The histograms presented here have a turnover near ∼1 M_☉, whereas the IMFs derived from observations in several star-forming regions tend to have a turnover near ∼0.07–0.5 M_☉ (e.g., Kroupa et al. 2013). A possible interpretation for the discrepancy is that the our simulations do not spatially resolve the formation of very low mass stars whose masses are ≲ 0.1 M_☉, e.g., brown dwarfs. As shown in Figure 13, some sink particles have ∼0.1 M_☉ when they are produced. If a finer resolution is adopted, the initial masses of the sink particles could be lower. This could shift the turnover to lower mass.

Zoom In Zoom Out Reset image size

Figure 13 shows the evolution of sink particle mass for the four representative models. As also shown in Table 1, the number of sink particles is larger in the high-density models than in the low-density models. The high-density models and the models with colliding flows take a shorter time for the sink particles to reach 10 M_☉ after the initial sink particle formation than do the other models. However, the typical accretion rates (slopes of the lines) do not exhibit a significant difference between models. Among all sink particles, those reaching 10 M_☉, as indicated by the thick lines, show slightly higher accretion rates (steeper slopes) for the high-density and/or colliding models.

Figure 14 shows the accretion rates for the sink particles explicitly. The accretion rates are spread over a wide range of $10^{-7}\,M_\odot \,\mathrm{yr}^{-1} \lesssim \dot{M} \lesssim 10^{-4}\,M_\odot \, \mathrm{yr}^{-1}$, which is a common range for the models. The accretion rate for each particle oscillates within this range. This range includes the typical accretion rate for gas of 10 K, i.e., $c_s^3/G = 2 \times 10^{-6}\,M_\odot \, \mathrm{yr}^{-1}$, while the majority of sink particles exceed this value. The most massive particle tends to have the highest accretion rate in each model (see thick lines in Figures 13 and 14). There is also a tendency that a model producing more stars has higher maximum mass of stars when observed at a given epoch. For example, among the four models shown in Figure 13, model HT10F10 has the highest maximum stellar mass, and this model produces the most number of stars at 0.1 Myr after the first sink particle formation. This indicates that the maximum mass of stars correlates with the number of stars formed there. A similar implication is proposed by the observational research of Dobashi et al. (2001), in which the maximum stellar luminosity of protostars correlates with the mass of the parent clouds. Moreover, some other observations have shown that the number of protostars also correlates with the mass of the parent cloud (Dobashi et al. 1996; Yonekura et al. 1997; Kawamura et al. 1998). When we simply assume that the stellar luminosity correlates with the stellar mass, these observations also show that the maximum stellar mass correlates with the number of stars in agreement with our simulations.

Zoom In Zoom Out Reset image size

Figure 15 shows the star formation efficiency as a function of the freefall time for each model. The star formation efficiency is defined as the ratio of the mass of all sink particles to the total mass of all particles and gas. Star formation in models with colliding flows (models HT10F10 and LT10F10) begins in earlier stages than for models without colliding flow (models HT10F0 and LT10F0). In comparing models with the same colliding flow, the high-density models tend to begin star formation earlier than do the low-density models. After star formation begins, the star formation efficiencies increase with a similar trend in all models. This similar trend indicates that the high-density models have higher accretion rates in total than do the low-density models because the clouds of the high-density models are more massive than those of the low-density models. These high accretion rates are mainly due to the large number of sink particles formed in the high-density models.

Zoom In Zoom Out Reset image size

Sink particles are produced with initial velocities so that the linear momenta of the material gas and the formed sink particles are conserved. Sink particles change their velocities via gas accretion onto the sink particles and gravitational interaction with the gas and other particles. Figure 16 shows the velocity distributions for sink particles in bar charts for the four representative models HT30F0, HT10F30, HT10F10, and LT10F0 at the stage of M_{⋆, max} = 10 M_☉. The high-density models have many sink particles, and each model exhibits a profile of velocity distribution with a peak and tails at both low and high velocities. In contrast, the low-density models do not have sufficient numbers of sink particles to show such a clear profile.

Zoom In Zoom Out Reset image size

The velocity distributions for gas are plotted for comparison in Figure 16. For the high-density models HT30F0 and HT10F30, the velocity distributions for sink particles are similar to those for gas; the peak positions and widths of the profiles coincide between the sink particles and the gas. For model HT10F30, the velocity distribution for gas has a wing at high velocities because of the accretion flow onto the sheet cloud. The accretion flow does not considerably affect the velocities of the sink particles because the accretion flow has low density and the sink particles exist around the high-density sheet clouds. Model HT10F10 also has an accretion flow onto the sheet, but the flow is denser and slower than that for model HT10F30 (see Figure 4 for comparison of the densities of accretion flows between these models). The profile of the velocity distribution for gas therefore appears to shift to the high velocities with respect to the velocity distribution for the sink particles.

For all low-density models, the velocity distributions for sink particles roughly agree with those for gas, even though the number of sink particles is small. Note that sink particle formation occurs after the sheet due to the colliding flow disappearing for the low-density models. Even for the colliding models, the velocities of sink particles follow the velocity distribution for gas.

The fitting curves for both sink particles and gas are also plotted in Figure 16 under the assumption that the velocities of sink particles and gas follow a Gaussian distribution. When each velocity component, (v_x, v_y, v_z), follows a Gaussian distribution with a velocity dispersion, σ, the distribution for the absolute value of the velocity $|{\boldsymbol {v}}|$ is given by

$$\begin{equation} p(|{\boldsymbol {v}}|; \sigma) = \frac{2 |{\boldsymbol {v}}|^2}{(2 \pi)^{1/2} \sigma ^3} \exp \left(-\frac{ |{\boldsymbol {v}}|^2 }{2 \sigma ^2}\right). \end{equation} \tag{ 7 }$$

The velocity dispersion for the sink particles is evaluated by

$$\begin{equation} \sigma _\mathrm{star} = \left(\frac{1}{ 3N_\star } \sum _{i=1}^{N_\star } |{\boldsymbol {v}}_i|^2 \right)^{1/2}, \end{equation} \tag{ 8 }$$

and that for gas is evaluated by using a mass-weighted integral,

$$\begin{equation} \sigma _\mathrm{gas} = \left(\frac{ \int \rho ({\boldsymbol {r}}) |{\boldsymbol {v}}({\boldsymbol {r}})|^2 dV}{ 3 \int \rho ({\boldsymbol {r}}) dV} \right)^{1/2}. \end{equation} \tag{ 9 }$$

The numerical factor of three in Equations (8) and (9) comes from the number of spatial dimensions.

The fitting curves indicate that the velocities of sink particles and gas are well described by a Gaussian distribution for model HT30F0. The velocity dispersion for sink particles is roughly the same as that for gas (σ_star ≃ σ_gas). Such a tendency is confirmed in all the high-density models without colliding flow. In the high-density models with colliding flow, gas has a significantly larger velocity dispersion than that for sink particles (σ_gas ≳ σ_star) because of the high-velocity wing due to the colliding flow. Due to the large velocity dispersion, the fitting curve for gas shifts to high velocities with respect to sink particles. In these models we also confirmed that, for sink particles, the velocity dispersion for v_x tends to be smaller than those for v_y and v_z. The dispersion difference indicates that the gravity of the sheet decelerates the motion of the sink particles along the x-direction, which is the direction normal to the sheet.

Sink particles and gas change their velocity distributions as time proceeds. For model HT30F0, a considerable decrease in the velocity dispersion for gas is seen in Figure 6, even after formation of the first sink particle. In this model, the velocity dispersion for gas decreases: σ_gas/c_s = 10.76 and 7.41 at the stages of M_{⋆, max} = 1 M_☉ and 10 M_☉, respectively. The velocity dispersion for sink particles also decreases: σ_star/c_s = 9.09 and 7.09 at the stages of M_{⋆, max} = 1 M_☉ and 10 M_☉, respectively. This indicates that stars with a velocity dispersion similar to that for gas are produced continuously. On the other hand, for model HT10F0, the velocity dispersion for gas decreases slightly, as seen in Figure 6, and both σ_star and σ_gas are roughly constant.

Reconstruction of channel maps by using the data of the numerical simulations is useful for comparison of observations and numerical simulations, as demonstrated by Dobashi et al. (2014). The channel maps here are reconstructed for four models: HT10F0, HT10F10, LT10F0, and LT10F10. In reconstructing a channel map, gas is assumed to be optical thin and both the non-thermal and thermal velocity components are considered; each numerical cell has a velocity along the line of sight, v_los, with a Gaussian profile, which is broadened by the sound speed, and the central velocity of the Gaussian is equal to the bulk velocity. The x-direction is adopted as the line of sight, which is parallel to the direction of the colliding flow.

Figure 17 shows composite maps for the two velocity channels. In each composite map, two representative velocity channels of −1.50 km s⁻¹ ⩽ v_lols ⩽ −0.51 km s⁻¹ and 0.51 km s⁻¹ ⩽ v_lols ⩽ 1.50 km s⁻¹ are shown in cyan and red colors, respectively. Components existing in both velocity channels are expressed in white because cyan and red are complementary colors.

Zoom In Zoom Out Reset image size

The most prominent features are shown in model HT10F10 (Figure 17(b)), where red and cyan filaments are tangled up with each other on a 0.1 pc scale. Moreover, the red and cyan filaments show an anti-correlated distribution because of the collision of filaments, as shown by Dobashi et al. (2014). Similar anti-correlation between the different velocity components was also reported by Shimoikura et al. (2013). A considerable number of sink particles exist near the region of collisional interaction, and this indicates that the collision of filaments induces star formation. Stars are also formed along the dense filaments, e.g., the red filament running from (y, z) = (− 0.7, −1.3)pc to (− 0.2, −0.6)pc.

In contrast, models HT10F0, LT10F0, and LT10F10 exhibit coarser patterns of red and cyan colors on a parsec scale (Figures 17(a), (c), and (d)). These patterns indicate that velocity shear on the cloud scale exists. This is attributed to turbulence, which has a larger kinetic energy in a longer wavelength mode. Therefore, the cloud-scale gas motion should dominate over short wavelength modes when considering turbulence with the power spectrum examined here. We confirmed that similar anti-correlations on the cloud scale are also shown irrespective of the choice of line of sight for the models without colliding flows.

The low-density colliding model LT10F10 has a sheet cloud caused by the colliding flow before sink particles form. The sheet cloud disappears at the stage shown in Figure 17(d). When the sheet cloud exists, the channel map shows the anti-correlation pattern on a sub-parsec scale, as also seen for model LT10F10. However, the contrast for the anti-correlation pattern is much weaker than that for the high-density model HT10F10 because the filaments in the low-density model have a low density contrast.

Figure 18 shows a position–velocity (PV) diagram for examining the slit region in Figure 17(b). Every sink particle is associated with a dense gas component in the PV diagram. Figure 18(a) shows the PV diagram for the region of −2.5 pc ⩽ z ⩽ −1.25 pc, where the interaction between clouds with different velocities induces the star formation shown in Figure 17(b). The PV diagram clearly shows that dense gas components pass through each other with supersonic velocities in the x-direction (line of sight).

Zoom In Zoom Out Reset image size

In Figure 18(b), two dense components with different velocities exist in −0.5 pc ≲ z ≲ −0.4 pc without an associated sink particle. This figure shows that these components have not yet collided. The components are shown in white in the composite map of Figure 17(b). Indeed, we confirmed that these components collide 0.2 Myr later, and star formation is induced.

The regions shown in Figures 18(c) and (d) have only two sink particles, and active star formation is not induced there. In this region, anti-correlation between the dense cyan and the thin red components is observed at z ≃ 1.7 pc and a sink particle is associated (Figure 17(b)). The PV diagram suggests that the sink particle is associated with the cyan component, and it interacts with the ambient gas with low density (Figure 18(d)).

5.1.1. Probability Distribution Functions

5.1.2. Velocity Distributions

As shown in the analysis of the velocity power in Section 4.6, non-gravitational non-colliding model HT10F0wog keeps E_long/E_tot = 0.15–0.2. This ratio agrees with that obtained by Federrath et al. (2010b), who found that non-gravitational isothermal driven turbulence with ${\cal M}_t \simeq 5.5$ has a ratio of E_long/E_tot ≃ 0.2 when the driving force is 1/3 or less of the longitudinal mode, i.e., F_long/F_tot ⩽ 1/3 (see their Figure 8). The ratio of E_long/E_tot ≃ 1/3 is expected from geometrical considerations. (e.g., Nordlund & Padoan 2003). Federrath et al. (2010b) showed that the velocity Fourier spectra shows E_long(k)/E_tot(k) ≃ 1/3 in the inertia range of k for the solenoidal driving force. For the decaying turbulence presented here, the ratio of E_long(k)/E_tot(k) in the inertia range tends to decrease from 0.4 to 0.25 as time proceeds, and the ratio agrees with the case of the turbulence driven by solenoidal force.

In the case with colliding flow, model HT10F10wog shows a higher value of E_long/E_tot in the early stage. This result suggests that the colliding flow provides more compressibility. The compressible flow shows a rapid decay, and E_long/E_tot = 0.2–0.25 in a long period of (2–5) t_cross, which are slightly higher values than those in the non-colliding model. Compared between models HT10F0wog and HT10F10wog, E_trans shows almost the same values during the decay of turbulence, while E_long exhibits slightly higher values in model HT10F10wog than in model HT10F0wog after the rapid decay of E_long in the early stages of model HT10F10wog. The slightly higher value of E_long/E_tot in the colliding model is due to the slightly higher value of E_long. Federrath et al. (2010b) showed that non-gravitational turbulence shows E_long/E_tot = 0.6 when the driving force is purely compressive. This ratio is higher than the colliding model presented here. This is because compressibility is continuously supplied by the driving force in the case of Federrath et al. (2010b), while it is provided by the colliding flow only in the initial stage in our model.

In our simulations, the gravitational collapse leading to star formation occurs locally compared with the cloud scale, and the density structure on the cloud scale is mainly a result of the interaction of turbulent flows. However, self-gravity enhances the compressibility of the velocity field in the whole cloud. For example, models HT3F0 and LT3F0 exhibit E_long/E_tot ≃ 0.6 at the last stages. Such a high compressibility of the velocity field requires an extreme driving force of F_long/F_tot ≳ 0.9 in the case of non-gravitational turbulence (Federrath et al. 2010b).

Hansen et al. (2011) examined the decay of anisotropic turbulence and showed that the decay rate of the turbulence depends on the crossing time of the isotropic component only, and the decay rate is much lower for anisotropic turbulence than for isotropic turbulence. For their highly anisotropic turbulence, $\langle |\nabla \times {\boldsymbol {v}}|^2\rangle$ is much larger than $\langle |\nabla \cdot {\boldsymbol {v}}|^2\rangle$, and so it seems that the solenoidal flow dominates over the compressive flow. In our case, models with colliding flows can be considered as anisotropic turbulent models, in which the colliding flow corresponds to a Fourier component of k_x = 1 in the compressive mode. In contrast to the anisotropic models of Hansen et al. (2011), the velocity dispersion due to the colliding flow decays faster than that for the isotropic turbulence (see Figure 6(b)). The anisotropic turbulence therefore decays faster than the isotropic turbulence because anisotropy is responsible for the compressible mode.

The magnetic field is ignored in this paper, although it is known that the magnetic field plays an important role in star formation.

Many simulations have been performed for turbulent clouds while taking into account magnetic fields with and without self-gravity (e.g., Mac Low 1999; Ostriker et al. 2001). When considering a relatively strong magnetic field (plasma beta ≃ 0.2), a turbulent cloud tends to collapse along the magnetic field lines, and then a sheet structure oriented perpendicular to the magnetic field tends to form (e.g., Nakamura & Li 2008). Even when strong turbulence of ${\cal M}\sim 10$ is imposed on the initial stage, rapid decay of the turbulence leads to formation of a sheet structure perpendicular to the magnetic field.

In the magnetically critical case, the magnetic field strength is given by B_cr = 2πG^1/2Σ (Nakano & Nakamura 1978; Tomisaka et al. 1988), where Σ denotes the column density for a typical scale and is estimated here as Σ = ρ₀λ_J. The magnetic energy of the critical field strength is therefore given by $B_\mathrm{cr}^2/(8\pi) = (\pi ^2/2) \rho _0 c_s^2 = 4.9\, \rho _0 c_s^2$. This indicates that a flow faster than ∼3c_s dominates over the magnetic field, and a sheet forms due to compression by the flow rather than the magnetic field.

Chen & Ostriker (2014) performed MHD simulations of turbulent clouds with a colliding flow, taking into account ambipolar diffusion. Their model settings are similar to ours except for the magnetic field; they assumed an initial number density of n = 10³ cm⁻³ and a colliding flow of Mach 10. Because of the high Mach number of the colliding flow, a sheet forms due to the compression of the flow, and the gas and the magnetic field are accumulated into the sheet. The sheet is perturbed little compared to our cases because they assume relatively weak turbulence. The evolution of the sheet cloud is highly dependent on the ionization rate and the direction of the magnetic field. The models with a low ionization rate exhibit several features similar to the hydrodynamical models, e.g., the post-shock densities, collapse times, masses of cores. Thus, several features of our hydrodynamical models may hold when our models are extended to strongly diffusive MHD.

When a high-density portion undergoes gravitational collapse, a sink particle is created there. Sink particle creation is implemented mainly according to Federrath et al. (2010a); the conditions for particle creation are that (1) the density is higher than the threshold density, ρ_sink, at the position of the new particle, (2) the gravitational potential takes the local minimum there, (3) the divergence of the velocity is negative ($\nabla \cdot {\boldsymbol {v}} < 0$), (4) all the eigenvalues of the symmetric parts of the velocity gradient tensor $\nabla {\boldsymbol {v}}$ are negative, (5) the total energy of the gas within the sink radius r_sink is negative (E_th + E_kin + E_grav < 0), and this indicates that the gas is gravitationally bound there. In addition, sink particle creation is forbidden at a position within a distance of 2r_sink from other sink particles in order to avoid an overlap of sink regions, where r_sink denotes a sink radius, which is set at r_sink = 4Δx in the finest grid level (cf. Krumholz et al. 2004). These conditions are checked at every hydrodynamical timestep, and a sink particle is created when the conditions are satisfied. The merger of sink particles can also be implemented, but this function is switched off in the simulations presented in this paper.

The method for accretion is the same as that of Machida et al. (2010). A sink particle accretes gas that is located within the distance of the sink radius, r_sink, from the sink particle. The accretion onto sink particles is performed at each hydrodynamical timestep. The mass accreted onto the ith sink particle is given by

$$\begin{equation} \dot{M} \Delta t = \int _{|{\boldsymbol {r}} - {\boldsymbol {r}}_i| < r_\mathrm{sink}}\Delta \rho ({\boldsymbol {r}}) dV, \end{equation} \tag{ A1 }$$

$$\begin{equation} \Delta \rho ({\boldsymbol {r}}) = \max \left[ \rho ({\boldsymbol {r}}) - \rho _\mathrm{sink}, 0\right], \end{equation} \tag{ A2 }$$

where $\dot{M}$ and Δt denote the mass accretion rate and the hydrodynamical timestep, respectively. The integration of Equation (A1) is a volume integration within a sphere with a radius of r_sink and a center of ${\boldsymbol {r}}_i$, which coincides with the position of the ith sink particle. The sink particle obtains the mass given by Equation (A1), and gas reduces the density by $\Delta \rho ({\boldsymbol {r}})$ according to Equation (A2). In the accretion process, the gas velocity remains unchanged and the sink particle therefore obtains a linear momentum of

$$\begin{equation} \int _{|{\boldsymbol {r}} - {\boldsymbol {r}}_i| < r_\mathrm{sink}}\Delta \rho ({\boldsymbol {r}}) {\boldsymbol {v}}({\boldsymbol {r}}) dV. \end{equation} \tag{ A3 }$$

The accretion of the linear momentum changes the velocity of the sink particle. The mass and the linear momentum are conserved between the gas and the sink particles in the accretion process.

The sink particles are Lagrangian particles, moving according to the accretion of linear momentum described above and the gravitational interaction. The gravitational force of the ith sink particle is

$$\begin{equation} {\boldsymbol {g}}_{\mathrm{sink}, i}({\boldsymbol {r}}) = \left\lbrace \begin{array}{@{}ll}\displaystyle -\frac{G m_i}{|{\boldsymbol {r}}-{\boldsymbol {r}}_i|^3} ({\boldsymbol {r}}-{\boldsymbol {r}}_i) \;\; & {({\rm for} |{\boldsymbol {r}}-{\boldsymbol {r}}_i| \ge r_\mathrm{soft})}\\ \displaystyle -\frac{G m_i}{r_\mathrm{soft}^3} ({\boldsymbol {r}}-{\boldsymbol {r}}_i) \;\; & {({\rm for} |{\boldsymbol {r}}-{\boldsymbol {r}}_i| < r_\mathrm{soft})} \end{array} \right., \end{equation} \tag{ A4 }$$

where m_i denotes the mass of the ith sink particle. Equation (A4) represents the gravitational force of a uniform sphere with a radius of r_soft, which corresponds to a softening radius and is set at r_soft = r_sink here. The gas is accelerated due to the self-gravity of the gas and the sum of the gravitational forces of all the sink particles:

$$\begin{equation} {\boldsymbol {g}}({\boldsymbol {r}}) = -\nabla \Psi ({\boldsymbol {r}}) + \sum _i {\boldsymbol {g}}_{\mathrm{sink}, i}({\boldsymbol {r}}). \end{equation} \tag{ A5 }$$

When computing the gravitational forces for a sink particle acting on the cells within the softening radius ($|{\boldsymbol {r}}-{\boldsymbol {r}}_i| < r_\mathrm{soft}$), each cell is subdivided into 8³ subcells, and the gravitational forces of the sink particle acting on the subcells are summed (Krumholz et al. 2004).

The gravitational force of gas acting on the sink particle is evaluated as a reaction force of ${\boldsymbol {g}}_{\mathrm{sink},i}$ (see Krumholz et al. 2004). This ensures conservation of the linear momentum through the gravitational interaction between the gas and the sink particle within a round-off error. An alternative method for evaluation of the gravitational force of gas acting on the sink particle is taking the average of the gravity of the gas over the spherical region with the softening radius:

$$\begin{equation} {\boldsymbol {g}}_{\mathrm{gas}, i} = - \frac{ \int _{|{\boldsymbol r} - {\boldsymbol r}_i| \, <\, r_{\rm soft}} \nabla \Psi ({\boldsymbol {r}}) dV }{ \int _{|{\boldsymbol r} - {\boldsymbol r}_i| \, <\, r_{\rm soft}} dV }, \end{equation} \tag{ A6 }$$

where each cell is subdivided into 8³ cells for integration. This alternative method is valid because Equation (A4) coincides with the gravity of a uniform sphere. Finally, the gravitational force acting on the ith sink particle is given by

$$\begin{equation} {\boldsymbol {g}}_i = {\boldsymbol {g}}_{\mathrm{gas},i} + \sum _{j\ne i} {\boldsymbol {g}}_{\mathrm{sink},j}({\boldsymbol {r}}_i). \end{equation} \tag{ A7 }$$

The second-order leapfrog scheme is adopted for the time integration of sink particles. Sink particles and hydrodynamics generally have a common timestep. The timestep Δt is restricted so that the travel distances of the particles in a step are less than the smallest cell width of the numerical grid,

$$\begin{equation} \Delta t = \min (\Delta t_\mathrm{CFL}, \Delta t_\mathrm{vs}) \end{equation} \tag{ A8 }$$

$$\begin{equation} \Delta t_\mathrm{vs} = C_\mathrm{vs} \min _{i}\left[ \frac{\Delta x }{ v_{i,x} }, \frac{\Delta y }{ v_{i,y} }, \frac{\Delta z }{ v_{i,z} } \right], \end{equation} \tag{ A9 }$$

where Δt_CFL denotes a timestep determined by the Courant–Friedrichs–Lewy (CFL) condition, and ${\boldsymbol {v}}_i = (v_{i,x}, v_{i,y}, v_{i,z})$ denotes the velocity of the ith sink particle. The constant C_vs is set at 0.5 here. When a sink particle is accelerated by strong gravity, the timestep is divided into a sub-timestep of

$$\begin{equation} \Delta t _\mathrm{sub} = C_\mathrm{gs} \min _{i, j \ne i} \left[ \left(\frac{ \min (|{\boldsymbol {r}}_i - {\boldsymbol {r}}_j|, \Delta x, \Delta y, \Delta z) }{ |{\boldsymbol {g}}_i | } \right)^{1/2} \right]. \end{equation} \tag{ A10 }$$

The constant C_gs is set at 0.1 here. To reduce the computational cost, the gravity of gas acting on sink particles is approximated to be unchanged during the sub-timesteps. Similar sub-cycling is also adopted in Krumholz et al. (2004) and Federrath et al. (2010a).