2.1. The Initial Conditions of the Simulation Particles

2.2. Initial Energies

In a particle-based simulation, the thermal, kinetic and gravitational energies are given by

where Pi is the pressure and Φ_i is the gravitational potential at the location of particle i, with velocity vi and mass m_i. It should be emphasized that all of the SPH particles of a simulation must be used in the summation of equation 3.

Let α be defined as the ratio of the thermal energy to the gravitational energy and let β be the ratio of the kinetic energy to the gravitational energy, so that

and

The value of the speed of sound c₀ has been fixed at 225,000 cm/s, so that the initial turbulent cloud has the α₀ ratio given by 0.16, for all the collision models. The multiplicative constant mentioned in 2.1.2 has been adjusted so that the initial turbulent cloud has a β₀ ratio given by 0.5.¹

^{Higuchi et al. (2014)} presented line emission observations of the Brick using the Atacama Large Millimeter/Submillimeter Array and considered values of β₀ = 0.1 and α₀ = 0.02 taking into account a cloud mass of 2 × 10⁵ M_⊙, radius of 2.8 pc, a temperature of 20 K and a one-dimensional velocity dispersion of 4 km/s.

The virial parameter is very useful when characterizing the physical state of a gas structure, which is defined observationally by

where G is Newton’s gravitational constant, M and R are the mass and radius of the gas structure, and σ_1D is the intrinsic one-dimensional velocity dispersion of the hydrogen molecule. Assuming isotropic motions, a 3D velocity dispersion can be simply related by σ_3D = √3σ_1D. It should be noted that a gas structure in virial equilibrium would have β_vir = 1.

The empirical relation between the virial parameter is β_vir = 2 a β, where β is the dimensionless ratio defined in equation 5 and a is a numerical factor that is empirically included to take modifications of non-homogeneous and non-spherical density distributions into account. According to this empirical relation, the virial parameter of the simulation of this paper is approximately 1.

Later, the virial theorem will be useful to show the level of virialization of the simulation outcome. In general terms, for a gas structure in virial equilibrium, the energy ratios defined above in equations 4 and 5 satisfy the relation

It is expected that if a gaseous system has α + β > 1/2, then it will expand; in the other case, if α + β < 1/2, then the system will collapse. It must be mentioned that ^{Miyama et al. (1984)}, ^{Hachisu and Heriguchi (1984)} and ^{Hachisu and Heriguchi (1985)} obtained a criterion of the type α × β < 0.2 to predict the output of a given simulation.

2.3. The Azimuthal Velocity

Let us consider a massive agent, such as a dwarf spheroidal galaxy, which is located at the origin of a coordinate system. Let us place the molecular cloud of interest here to be in the z-axis, at a distance Z_C. This massive agent induces an escape velocity at each radius R (with respect to the center of the massive agent), so that

where M (R) is the mass contained up to the radius R and G is Newton’s gravitational constant. In spherical coordinates R, θ and φ, the velocity vector V→ has components V_R, V_θ and V_φ, respectively. These spherical components are related to the Cartesian components of velocity V_x, V_y and V_z by three simultaneous equations, whose coefficients are given in terms of the sine and cosine of the polar and azimuthal angles θ and φ, as follows:

Let us assume that Z_C and the radius of the cloud (which is defined in § 2 R₀ = 3 pc with respect to the center of the cloud), satisfy the relation R₀ / Z_C << 1, then as sin(θ) = R₀ / R ≈ R₀ / Z_C << 1 then θ ≈ 0, for which cos(θ) ≈ 1 and sin(θ) ≈ 0. In addition, in the particular case that the cloud follows a circular orbit around the massive center at radius R, then the velocity vector would only have a non-zero polar velocity component, V_θ, while the radial and azimuthal components V_R and V_φ are zero. Let us denote this velocity as V_θ = V_cir, where the minus sign indicates that the assumed rotation of the cloud is counter-clock-wise in its orbit around the massive agent. Under these simplifications, the relations 9 between Cartesian and spherical components of velocity are reduced to

so that the magnitude of the velocity of a particle located at any radius R is therefore always given by Vx2+Vy2=Vcir. These Cartesian components of the velocity will be added to the particle velocity defined in equation 1 to generate four new models in which this approximation will be implemented. It should be noted that the azimuthal angles have the same projection in both coordinate systems: the first is based on the cloud center and the second is located at the gravitational agent center.

We will call this velocity the “azimuthal velocity” because it is given only in terms of the azimuthal angle φ. Then, the approximation that replaces the tidal force by an azimuthal velocity, as described in equation 10, does not depend explicitly on the distance of the massive center to the cloud as long as ratio between the cloud radius to this distance is quite small.

Following with the model of a cloud of the CMZ, the mass of this massive agent has been fixed at M_B = 3.6 × 10⁶ M_⊙, which corresponds to the black hole located at the center of the Milky Way. In this case, the escape velocity at 500 pc is 5.57 km/s. To compare this velocity with those displayed at Figure 2, in terms of the speed of sound c₀ defined in § 2.2, we have a Mach number of V_cir = 2.47. It must be noted that the magnitude of this circular velocity V_cir is quite small as compared to that proposed by ^{Molinari et al. (2011)}, in which a model for the orbit of the gas stream near the massive center Sgr B2 is around 80 km/s. For this reason, we have also included a second set of models in which the massive agent is considered to be molecular gas concentrated in the nuclear bulge of the Milky Way, see ^{Launhardt et al. (2002)}, so that the total mass is M_H = 8.4 × 10⁸ M_⊙, see also ^{Mills (2017)}, for which the escape velocity at 500 pc is 85 km/s, such that the normalized velocity in terms of the speed of sound 37.82.

Fig. 1 Schematic diagram of the pre-collision geometry for a head-on collision (left-hand panel) and an oblique collision (right-hand panel). The color figure can be viewed online. 

Fig. 2 Distribution function for the radial component of the velocity at the initial snapshot at t/t_ff = 0, for the model (top left-hand) U with a low level of turbulence; (top right-hand) Ur with a small azimuthal velocity; (bottom left-hand) Ub with a high level of turbulence and (bottom right-hand) Urb with a large azimuthal velocity. f is the fraction of particles whose magnitude of the velocity vr/c0 (normalized with the sound speed) is smaller than that value shown in the horizontal axis. The color figure can be viewed online. 

2.4. The Collision Models

It is important to emphasize that the cloud entirely contains the two subsets of particles that are going to collide. Let us call these subsets the precollision sub-clouds. They are located initially along the x-axis, so that the centers are: for the left-hand clump (-2.55, 0, 0) pc and for the right-hand clump (2.55, 0, 0) pc.

The radius of the pre-collision sub-clouds are chosen to be equal for two models, and different for another two models. The former models are head-on collisions. An impact parameter b has also been considered for the latter models, so that they are oblique collisions, in which b takes the value 1.5 pc along the y-axis. ^{Bekki and Couch (2004)} have demonstrated observationally that the most likely impact parameter b in cloud collisions in the Large Magellanic Cloud and the Small Magellanic Cloud is 0.5 D < b < D where D is the diameter of the cloud. For the radius R₀ of the cloud considered in this paper, b has been chosen such that b = 0.25 D.

We show all these models in Table 1. The label is shown in Column one. The impact parameter value is shown in Column two. In Columns three and four, the relationship of the radius and translational velocities are shown, for the left-hand and right-hand sub-clouds, respectively. It is important to emphasize that the relative pre-collision velocity of the sub-clouds is 29 km/s and is formed for nonidentical velocities for the left-hand and right-hand sub-clouds. In Column five, the value of the ratio of the kinetic energy to the gravitational energy for the initial configuration of particles is shown, see equation 5. Finally, Column six gives the value of the azimuthal velocity added to the cloud particles, see § 2.3. It must be clarified that these models are a sample from a larger set of models that was considered in a first manuscript, so that the numbers of the labels do not show any ordering.

2.5. Characterization of the Initial Turbulence

To show the nature of the turbulence that is implemented in § 2.1.2, we consider the distribution functions of the initial velocity.

In Figure 2, we show the distribution functions of the radial component of the velocity at the initial snapshot, so that in the vertical axis the fraction f of the simulation particles whose magnitude of the velocity is smaller than that value shown in the horizontal axis. The radial component has been calculated with respect to the origin of the coordinates located in the center of the cloud, which is located at the center of the simulation box.

According to the left-hand column panels of Figure 2 that is, for models U and Ub, half of the simulation particles have a negative radial velocity component, while the other half have a positive radial component. This symmetry is expected from the random process described in § 2.1 to generate the direction of the velocity vectors.

It must be emphasized that for the right-hand column panels of Figure 2 that is, for models Ur and U_rb, which include an azimuthal velocity, the symmetry of the curves with respect to the positive and negative radial components has been lost. These panels indicate that the azimuthal velocity favors that 80 percent of the particles have negative radial component of the velocity.

It must also be emphasized that both types of models U, Ur and Urb have the same initial turbulent velocity spectrum with the same level of energy, as defined in § 2.1 and § 2.2, and have the same translational velocity. The only difference between them is whether or not they include the azimuthal velocity, as described in § 2.3. It is observed in Figure 2 that the distribution function of the models U shows a magnitude of the velocity 12 percent smaller than that of models Ur.

Later, we will compare these curves at the initial snapshot with curves obtained for an snapshot of the final evolution stage.

A brief description of the dynamics of the isolated cloud is given in § 5.1.

2.6. The Evolution Code

The simulations of this paper are evolved using the particle-based Gadget2 code, which implements the SPH method to solve the Euler equations of hydrodynamics; see ^{Springel (2005)}. Gadget2 has a Monaghan-Balsara form for the artificial viscosity, see ^{Balsara (1995)}, so that the strength of the viscosity is regulated by setting the parameter α_ν = 0.75 and β_ν = ½. α_v, see Equations 11 and 14 in ^{Springel (2005)}. The Courant factor has been fixed at 0.1.

The SPH sums are evaluated using the spherically symmetric M4 kernel and so gravity is spline-softened with this same kernel. The smoothing length h establishes the compact support, so that only a finite number of neighbors to each particle contribute to the SPH sums. The smoothing length changes with time for each particle, so that the mass contained in the kernel volume is a constant for the estimated density. Particles also have gravity softening lengths ∈, which change step by step with the smoothing length h, so that the ratio ∈ / h is of order unity. In Gadget2, ∈ is set equal to the minimum smoothing length h_min, which is calculated over all particles at the end of each time step.

2.7. Resolution

^{Truelove et al. (1997)} demonstrated that the resolution requirement of a hydrodynamic simulation can be expressed in terms of the Jeans wavelength λ_J which is given by

where G is Newton’s gravitation constant, c is the instantaneous sound speed and ρ is the local density, so that a mesh-based simulation must always have its grid length scale l such that l < λ_J / 4.

^{Bate and Burkert (1997)} demonstrated that the resolution requirement for a particle-based code, the Jeans wavelength λ_J is better written in terms of the spherical Jeans mass M_J, which is defined by

so that an SPH code will produce correct results as long as the minimum resolvable mass m_r is always less than the Jeans mass M_J . The mass m_r is given by m_r ≈ M_J / (2N_neigh), where N_neigh is the number of particles included in the SPH kernel (i.e., the number of neighbors). Therefore, our simulations will comply with this resolution requirement if the particle mass m_p is such that m_p /m_r < 1.

As we mentioned in § 2.1, we have N = 13366240 SPH particles in each simulation and therefore m_p = 7.4 × 10⁻³ M_⊙. Now, if we consider that the highest peak density in our collision models is ρ_max = 5.0 ×10⁻¹² g/cm³, then the minimum Jeans mass would be given by (M_J )_min ≈ 0.594048M_⊙, so that we obtain m_r = 7.4 × 10⁻³ M_⊙. Thus, for that peak density the ratio m_p /m_r 1, and the Jeans resolution requirement is satisfied. In this sense, we are sure to avoid the growth of numerical instabilities or the occurrence of artificial fragmentation in all our simulations, up to densities smaller or equal than ρ_max. In the next sections, we will present our results in terms of a normalized density, so that log (ρ_max/ρ₀) is given by 7.9, where ρ₀ is the average density of the initial cloud, as we mentioned in § 2.

2.8. Equation of State

Most simulations in the field of collapse used an ideal equation of state or a barotropic equation of state (BEOS), as was proposed by ^{Boss et al. (2000)}:

where γ ≡ 5/3 and ρ_crit is a critical density, a parameter which we explain now. This BEOS takes into account the increase in temperature of the gas as it begins to heat once gravity has produced a substantial contraction of the cloud. In this paper, we also use this BEOS scheme for simplicity with a critical density ρ_crit = 5.0 × 10⁻¹⁴ g/cm³, which is 100 times smaller than the peak density considered in § 2.7 for the resolution requirement estimate; that is, ρ_max. However, it should be emphasized that it is only an approximation. Consequently, to describe correctly the transition from the ideal to the adiabatic regime, one needs to solve the radiative transfer problem coupled to gravity in a self-consistent way.

3.1. Evolution of the Density Peak

In Figure 3 we show the time evolution of the global density peak, irrespective of where the particle with the highest density is located in the simulation volume. As can be seen in this figure, all models collapse at different times (as expected).

Fig. 3 Evolution of the density peak for all the models (top left-hand) U with low level of turbulence; (top right-hand) Ur with a small azimuthal velocity (bottom left-hand) Ub with high level of turbulence and (bottom right-hand) Urb with a large azimuthal velocity. The color figure can be viewed online. 

Let us consider the top left-hand panel of Figure 3, the models with a low level of turbulence. The fastest collapse is that of Model U 13; followed by the collapse of Models U 9 and U 5, respectively. The slowest collapse is that of Model U 11.

This ordering seems to be a consequence of the gas dragging, which is caused by the asymmetry in radius and velocity. Consequently, as the gas flows, it is more difficult to condense by the action of the gravity. Model U 13 does not show any sign of gas dragging, because the gas remains of the head-on collision is still around of the pre-collision center. This is the reason why this model collapses first.

The density peak curve of Model U 9 very closely follows that of Model U 13; This happens because the right-hand sub-cloud acts as a primary member in the binary system formed, which is only slightly perturbed by the left-hand sub-cloud, that acts as a secondary member of the binary. The collapse takes place first in the primary, which is more massive.

The density peak curve of Model U 5 is the third to reach the collapse. This happens because the right-hand sub-cloud entirely swallows the left-hand sub-cloud during the collision, so it produces a mass perturbation in the central region, which must first settle down for the collapse to continue. The slowest collapse is that of Model U 11. This happens because the dragging of the colliding sub-clouds is maximum, so the mass does not stack easily.

Let us now consider the bottom left-hand panel of Figure 3, the Models Ub with a high level of turbulence. In this case, the behavior of all of the curves is the same as that explained earlier for the top panel, but for the panel of Model Ub there is a slight shift to the right to longer evolution times, above all for the Models U 5b, U 9b and U 13b. The collapse time of Model U 11b is almost similar to that observed in Model U 11.

We observe a very significant change in the time scale for the right column panels of Figure 3; in the top right-hand panel, Models Ur that include a low azimuthal velocity, the collapse is accelerated, so that the collapse time is shorter than the previous models U and Ub by 40 percent, approximately. In the bottom right-hand panel, Models Urb that include a high azimuthal velocity, the collapse is extremely fast and the time scale has been reduced to a range from 0.1 to 0.15 t/t_ff. In addition, the curves for Models Urb do not show any sign of the early random collisions between the SPH particles, due to the turbulence spectrum induced on each particle velocity.

3.2. Column Density Plots

The main outcome of the collision models is shown by means of column density plots of a thin slice of gas, parallel to the x-y plane. To make a proper comparison between the different models, we have selected for each model a snapshot whose peak density is such that log (ρ_max/ρ₀) ≈ 5.

In Model U 5, once the head-on collison between the left-hand and right-hand sub-clouds has taken place, the asymmetry in the original sub-clouds in both radius and translational velocity, makes the right-hand sub-cloud (the biggest and the fastest cloud) swallow and drag the left-hand sub-cloud (the smallest and slowest cloud). The resulting stirred gas oscillates from the left-hand side to the righthand side along the x-axis. The spatial symmetry in the original configuration of the right-hand subcloud is translated to the symmetry in the arms developed around the central cloudlet, as can be seen in the top left-hand panel of Figure 4. In § 3.1, this phenomenon was simply referred to as gas dragging, which was a useful way to explain the time of collapse by means of the density peak curves.

When the asymmetry in the radius is removed, the final result is a central cloudlet, that is elongated along the x-axis, with a strong bipolar outflow along the y-axis, as can be seen in the top right-hand panel of Figure 4, which is the outcome of Model U 13. Because of this result, we note that the asymmetry in the velocities of the collision model is not as important as the asymmetry in the radii with respect to the outcome of the simulations.

Fig. 4 Column density plots of the collision models U, for a thin slice of gas parallel to the x-y plane. The unit of length is one parsec. The models are shown in panels as follows: (top left-hand) U 5 (at time t/t_ff = 1.65 and peak density log (ρ _max/ρ₀) = 5.0); (top right-hand) U 13 (at time t/t_ff = 1.34 and peak density log (ρ_max/ρ₀) = 4.3); (bottom left-hand) U 9 (at time t/t _ff = 1.42 and peak density log (ρ_max/ρ₀) = 5.0); (bottom right-hand) U 11 (at time t/t_ff = 2.4 and peak density log (ρmax/ρ0) = 5.0). The color figure can be viewed online. 

Different results are obtained when an impact parameter is taken into account. In the case of Model U 9, as illustrated in the bottom left-hand panel of Figure 4, the asymmetry in the radii and velocities together with an impact parameter produce a weak binary system, in which the the left-hand subcloud (the smallest and slowest cloud) passes by and is attracted by the right-hand sub-cloud (the biggest and fastest cloud), so that part of the mass of the former is pulled out. Nevertheless, an arm is still visible around the remains of the pre-collision lefthand sub-cloud and a long arm also develops around the central cloudlet, which is analogous in origin to the arm formed in Model U 5.

In the case of Model U 11, when the symmetry in the radii of the pre-collision sub-clouds is restored but still in the presence of the impact parameter as in Model U 9, a binary system is formed as the main outcome, in which several gas bridges are seen to be strongly connecting the remains of the two colliding clumps, as can be seen in the bottom right-hand panel of Figure 4.

In Figure 5 we show the column density plots of Models Ub, with a high level of turbulence. As expected, there is a lot of similarity with the previous Models U, with a low level of turbulence, because the initial structure of the velocity spectrum is the same for both Models U and Ub. The only difference is the magnitude of the velocity. The larger magnitude of the velocity for Models Urb makes the arms and tails larger and better defined than those in Models U.

Fig. 5 Column density plots of the collision models Ub, for a thin slice of gas parallel to the x-y plane. The unit of length is one parsec. The models are shown in panels as follows: (top left-hand) U 5b (at time t/t_ff = 2.15 and peak density log (ρ_max/ρ₀) = 4.34); (top right-hand) U13b (at time t/t_ff = 1.73 and peak density log (ρ_max/ρ₀) = 8.63); (bottom left-hand) U9b (at time t/t_ff = 2.08 and peak density log (ρ_max/ρ ₀) = 5.20); (bottom right-hand) U11b (at time t/t_ff = 2.42 and peak density log (ρ_max/ρ₀) = 5.14). The color figure can be viewed online. 

Let us now consider the iso-density plots for Models Ur, which are shown in the Figure 6. In this case, we show two columns of density plots. In the righthand column of Figure 6, we show the density plots for the snapshots with almost the same density peak shown in Figure 4 to allow comparison with Models U and Ub. We only observe an homogeneous and spherical collapse as the final result of simulations Ur.

In the left-hand column of Figure 6, we choose snapshots of the first stage of evolution, to show the early development of a central lump of gas, which is a direct consequence of the azimuthal velocity added to the particle velocity, see § 2.3. This central lump of gas makes the collapse of the cloud faster, as can be seen in Figure 3, because it acts as a centrally located mass attractor.

Fig. 6 Column density plots of the collision models Ur, with an azimuthal velocity included, see § 2.3. The models are shown in panels (from top to bottom) as follows: (first line) Model U 5r, (left) at time t/t_ff = 0.24, peak density log (ρ_max/ρ ₀) = 0.79 and (right) at time t/t_ff = 1.05, peak density log (ρ_max/ρ ₀) = 5.2; (second top line) Model U9r, (left) at time t/t_ff = 0.25, peak density log (ρ_max/ρ₀) = 0.8 and (right) at time t/t_ff = 1.05, peak density log (ρ_max/ρ ₀) = 4.8; (third line) Model U 11r, (left) at time t/t_ff = 0.23, peak density log (ρ_max/ρ ₀) = 0.78 and (right) at time t/t_ff = 1.05, peak density log (ρ_max /ρ₀) = 4.9; (fourth line) Model U13r, (left) at time t/t_ff = 0.22, peak density log (ρ_max/ρ₀) = 0.74 and (right) at time t/t_ff = 0.92, peak density log (ρ_max/ρ ₀) = 4.7. The color figure can be viewed online. 

Finally, let us consider the iso-density plots for Models Urb, which are shown in Figure 7. We observe the formation of a massive lump of gas at the cloud center, at the beginning of the evolution, similarly to those shown in the left-hand column of Figure 6. However, for each Model Urb, the massive lump of gas is quite bigger than for models Ur. This is a direct consequence of the higher azimuthal velocity added to the particle velocity, see § 2.3. This massive central lump of gas makes the collapse quite faster than that observed for Models Ur, as we notice happens in the bottom right-hand panel of Figure 3.

Fig. 7 Column density plots of the collision models Urb, for a thin slice of gas parallel to the x-y plane. The unit of length is one parsec. The models are shown in panels as follows: (top left-hand) U 5rb (at time t/t_ff = 0.12 and peak density log (ρ_max/ρ₀) = 5.0); (top right-hand) U13rb (at time t/t_ff = 0.12 and peak density log (ρ_max/ρ₀) = 4.62); (bottom left-hand) U9rb (at time t/t_ff = 0.09 and peak density log (ρ_max/ρ ₀) = 4.63); (bottom right-hand) U11rb (at time t/t_ff = 0.15 and peak density log (ρ_max/ρ ₀) = 4.62). The color figure can be viewed online. 

3.3. 3D Rendered Plots

In Figure 8, we show the spatial structure of the models U using 3D plots, for the same time and density chosen for the snapshot shown in Figure 4. Until now, the figures displayed in § 3.2, have been cuts parallel to the equatorial plane of the initial sphere, so that around of 10,000 particles are included in the slice shown. For the 3D plots, all of the particles with a density greater than log (ρ_max /ρ ₀) ≈ 0.7 and located within the region [-2.5,2.5] in the three Cartesian coordinates x, y, z, entered in the 3D plots.

In this case, the number of particles that are used to make the 3D plots ranges from 181266 to 5616884. The log of the density is rendered in the 3D plots by assigning a color and a vertical bar located in the bottom right-hand corner of each panel. It must be noted that an arbitrary rotation is done on the Cartesian coordinates to show some of the details of the spatial structure.

In the top left-hand panel of Figure 8, we see the remains of Model U 5, view from the rear (along the positive x-axis), in which one can see an elongated bulb. In the region where the unequal sub-cloud collision takes place, on the negative side of the x-axis, one can see a thick disk of gas surrounding the elongated bulb. This structure looks like a mushroom pointing toward the negative x-axis.

In the top right-hand panel of Figure 8, we see the remains of Model U 13; recall that this collision is head-on along the x-axis between two equal sized sub-clouds. For this reason, one can see an elongated solid tube of gas along the x-axis, surrounded by an almost spherical gas region, which is formed by the particles bounced from the collision, most of which escape away along the y-axis in both directions, positive and negative.

In the bottom left-hand panel of Figure 8, we see the remains of Model U 9, in which two unequal sized sub-clouds have an oblique collision. One can see only the remains of each separate sub-clouds, after their close encounter, so one can notice that the bottom sub-cloud is almost destroyed by the tidal force caused by the top sub-cloud.

In the bottom right-hand panel of Figure 8, we see the remains of Model U 11, in which two equalsized sub-clouds have an oblique collision. For this reason, the symmetry is evident between the top and bottom gas tubes, which are formed as the tracks of the original colliding sub-clouds. There is a complex bridge of gas connecting these top and bottom tubes. A disk of gas is surrounding the bridge.

Fig. 8 3D-plots of the collision Models U for the same snapshots of Figure 4, rotated arbitrarily to reveal more details of the configurations, in a window of spatial dimension from -2.5 to 2.5, in each of the axis-xyz, in which only those particles with log (ρ_max/ρ₀) > 0.7 have entered to make the plot. The panels and the number of particles entered are shown as follows: (top left-hand) U 5 with 243455 particles; (top right-hand) U 13 with 287037 particles; (bottom lefthand) U 9 with 181266 particles and (bottom right-hand) U 11 with 5616884 particles. The color figure can be viewed online. 

3.4. Distribution Function of the Radial Component of Velocity

In Figure 9 we show in the vertical axis the fraction of particles with a velocity smaller than that shown in the horizontal axis. This distribution function of the velocity is taken at the same time and density as the snapshots shown in Figures 4, 5, 6 and 7. We consider only the radial component of the velocity, which is calculated with respect to the center of the cloud.

Fig. 9 Distribution function of the the radial component of the velocity, calculated with respect to the center of the cloud. The snapshots are taken at the same time and density as the snapshots shown in Figure 4, 5, 6 and Figure 7. (top left-hand) Models U with a low level of turbulence; (top right)-hand Ur with a low azimuthal velocity (bottom left-hand) Ub with high level of turbulence and (bottom right-hand) Urb with a high azimuthal velocity. The color figure can be viewed online. 

By comparing with the panels of Figure 2 in § 2.5, one can see that the fraction of the simulation particles with a negative component of the radial velocity has increased from 0.5 at time t/t_ff = 0 for all the models, to 0.9 for Models U , at time of the Figure 4; to 0.8 for Models Ub and to 0.9 for Models Ur. This means that a high fraction of the simulation particles have likely reached already (Models U 5, U 9 and U 13) or are flowing towards (Model U 11) an accretion center.

4.1 Radial Profile of the Density and Mass

In Figure 11 we show the radial profile of the density (in the left-hand column) and the mass (in the right-hand column), calculated with respect to the center of each cloudlet as defined in Figure 10.

It must be clarified that the density ρ_bin(r) and mass M(r)_bin, shown in the vertical axis of Figure 11, are determined by taking into account only those particles located within the radii r and δr, where δr is given by 4/500 pc per bin, and r goes from r = 0 (the cloudlet center) to r_max = 4 pc (even further than the edge of the cloudlet, because we want to study the environment of the cloudlets as well).

Let us consider the top line of Figure 11 for Models U. For Model U 11, there are two cloudlets. Consequently, we label the cloudlet located to the left of the vertical axis and above the horizontal axis, as Cloudlet “a” (in the upper left-hand region). We label the cloudlet located to the right of the y-axis and below the x-axis as Cloudlet “b”, as can be seen in the bottom right-hand panel of Figure 10.

The curves of ρ _bin for Models U5, U9, U11a and U11b are not steep, as opposed to the curve for Model U13. This means that the cloudlets of these models must have a mass increasing with radius, as can be seen in the right-hand panel of the top line of Figure 11.

Fig. 11 (Left-hand panel) Radial profile of the density and (right-hand panel) radial profile of the mass. In the vertical axis we show the mass and density averaged over those particles located within the radial spherical shell defined by r, so that these functions are not cumulative. In the horizontal axis, r starts at 0, the center of each cloudlet, as illustrated in Figure 10 for Models U and analogously calculated for all the other models. The panels on the top line are for Models U; the panels on the middle line are for Models Ur and the panels on the bottom line are for Models Urb. The snapshots are taken at the same time and peak density as in Figure 4, Figure 5, Figure 6 and Figure 7, respectively. The color figure can be viewed online. 

Model U13 is the only one that shows a density curve ρ _bin(r) decreasing significantly with radius r. For this behavior, the mass M(r)_bin(r) is almost kept constant for a wide range of radii. This would be the standard behavior of a dense cloudlet formed by gravitational attraction in a simulation.

In the second and third lines of Figure 11, from top to bottom, we show the curves for Models Ub and Ur, respectively. The behavior observed here is quite similar to the one described earlier for Model U. In the bottom line, we show the curves for Models Urb, which include a high azimuthal velocity and as we have seen in § 3.2, the simulation outcome changed significantly. In this case, the density curves are kept constant and the mass curves are slightly increasing functions of the radius, above all in the range of radius 0-3 pc.

^{Rathborne et al. (2015)} found curves for the radial profile of the mass and density of the Brick, such that the mass curve is always an increasing function of their effective radius, while the density curve is always a decreasing function, both curves extending up to a effective radius of 2 pc. For instance, the mass contained within its central one pc is approximately 6 10³ M_⊙, while the density curve follows a power-law over radii r^−1.2.

The increasing mass curves shown in the righthand column of Figure 11 are of the same order of magnitude as those reported by ^{Rathborne et al. (2015)}, taking into account that the curves of Figure 11 are not cumulative (as we mentioned earlier). Therefore, the exterior layers (outside of the densest central region) of the cloud, contain a substantial amount of mass.

4.2. Radial Profile of the Radial and Tangential Components of the Velocity

In Figure 12, we show the radial profile of the radial (left-hand column) and tangential (right-hand column) components of the velocity. We apply here the same radial partition described in § 4.1; that is, from the cloudlet center up to 4 pc.

Let us clarify the meaning of the tangential component of the velocity. In spherical coordinates (r, θ, φ), a gas particle has a magnitude of the velocity vector v_p with the components v_r,v_θ and v_φ. Then, we split the components of the velocity into radial v_r and tangential v_t = (v_θ + v_φ)/2, so that we can follow both components separately. We do this separation because the radial component can be associated with a collapse trend while the tangential component can be considered as a manifestation of turbulence, see ^{Guerrero and Vázquez (2020)}.

The left-hand column panels Figure 12 indicate that many particles move to the cloudlet center, mostly from the innermost region of the cloud. The right-hand column panels Figure 12 indicate that there is a non-zero, almost constant, tangential component of the velocity. These observations indicate that a lot of particles are falling towards the cloudlet center in trajectories that are slightly curved (i.e., not from a purely radial direction, such as in a freefall).

Fig. 12 (Left-hand column) Radial profile of the radial component of the particle velocity and (right-hand column) the radial profile of the tangential component of the particle velocity, both calculated with respect the center of each cloudlet illustrated in Figure 10 and normalized with the speed of sound. From top to bottom, the first line of panels is for Models U; the second line for Models Ub and the third and four lines for Models Ur and Urb, respectively. The snapshots considered for these calculations are taken at the same time and peak density shown in Figure 4, Figure 5, Figure 6 and Figure 7, respectively. The color figure can be viewed online. 

Let us recall the behavior of a test particle and let its velocity magnitude be given by v_g. This v_g is determined by v_g = √2 GM (r) / r, where M(r) is the mass contained up to radius r and G is Newton’s gravitational constant. This v_g can be considered as the velocity when a test particle arrives at distance r from the central mass M, having started from rest at infinity, where its gravitational potential is zero. As is well known, for a spatially bounded mass of radius R_g, the velocity v_g must increase with the radius r, such that 0 < r < R_g. Once the radial coordinate r is outside the bounded mass, that is, for r > R _g the velocity v_g simply decreases.

It must be noted that a significant fraction of the total velocity magnitude v_p comes from its radial component v_r, though a minor fraction of v_p comes from the tangential components grouped in v_t. Then, we can think about a curve of v_p if we see a curve of v_r, because we only need to transform from v_r to v_p by changing the velocity sign, from negative to positive values, so that the behavior of both curves, v_r and v_p, would be very similar.

Let us consider the panels of the top line of Figure 12, which are for Models U. In terms of the “imagined curves” of v_p (r), the cloudlets U 11a and U 11b follow the behavior expected for the test particle velocity, indicating that the cloudlet radius R (the analog of the bonded mass) is around 1 pc for both cloudlets of Model U 11. The curves for Models U 9 and U 13 indicate that the bounded mass has a very small radius R, which is slightly smaller than 0.5 pc. For Model U 5, the resolution of the radial partition is not fine enough to indicate a radius R of the cloudlet found. This bounded mass can be identified with the size of the region from the center of the cloud, which is a particle reservoir, out of which the particles flow towards the cloud center. The core of the collapsing cloud, which is formed by the particles with higher density of the simulation, is located in these cloud centers.

The panels of the second line of Figure 12, from top to bottom, show the curves for the Models Ub, which are very similar to those already described for Models U.

The panels of the third line of Figure 12, which are for Models Ur, with a low azimuthal velocity, indicate a behavior very similar to that observed for Models U 9 and U 13, that is, a region of particle reservoir is about 1 pc in radius from the center of the cloud. These behaviors can be better seen in the panels on the bottom line of Figure 12, which are for Models Urb. Models Urb have a high azimuthal velocity, whose effect is more clearly seen because the in-fall velocity is quite higher than in the previous Models U, Ub and Ur. Instead of a bounded mass, the size of the region of strong in-fall gas is determined by the size of the centrally located lump of gas induced by the azimuthal velocity, see the lefthand column of Figure 6, for an illustration.

4.3. Integral Properties of the Cloudlets

In Figure 13 we show the values of the dimensionless ratios α and β, respectively, calculated only for a cloud region that includes the cloudlets and their surroundings. Two parameters are used: the first parameter is log (ρ_min), which is a lower bound for density; and the second parameter r_max is a maximum radius, which is taken with respect to the cloudlet’s center. To calculate the ratios α_f and β_f for the cloudlets, we consider only those particles that have a density greater than log (ρ_min) and are located at a radius smaller than r_max. To make a comparison between all the models and to calculate the properties of all the models shown in Figure 13, we have used the following values log (ρ_min) = 0.0 and r_max = 1.5 pc.

One can see in the top left-hand panel of Figure 13 that Models U 5, U 9 and U 13 have a similar value for the ratio α_f, which is around 0.1. For the cloudlets “a” and “b” of Model U 11 (i.e., U 11a and U 11b) α_f is around 0.05. The values of the ratio β_f ranges from 0.1 for Model U 13, around 0.13 for Model U 5, and a little higher than 0.2 for Model U 9. Cloudlets U 11a and U11b have the highest values of β_f ≈ 0.55. Cloudlets U 11a and U 11b are the only ones over-virialized, because their sum α_f + β_f >1/2; while Models U 5, U 9 and U 13 are sub-virialized, because their sum α_f +β_f < 1/2. The same behavior is observed for Models Ub, with a high turbulence, as can be seen in the bottom left-hand panel of Figure 13.

In both the top right-hand panel and the bottom left-hand panel of Figure 13, we see that curves for the Models Ur and Ub, respectively, show a behavior that is very similar to that already observed for the curves of Models U. For all Models Ur and Ub the ratio α_f is around the value 0.09, with a clear tendency to lower values. The ratio β_f for these models ranges from 0.1 to 0.5 Models Ur, Ur5, Ur9 and Ur13 are sub-virialized, because their sum α_f + β_f < 1/2. Meanwhile, the two cloudlets of Model Ur11 are over-virialized, because their sum α_f + β_f > 1/2. The same behavior is observed for Models Ub.

On the opposite side, all Models Urb, as shown in the bottom right-hand panel of Figure 13, are overvirialized, because the high initial azimuthal velocity is manifested in the excess of kinetic energy, such that the β_f ratio for these models is in the range 7-11. In spite of this excess of kinetic energy, all Models Urb have collapsed.

We show the mass associated with the cloudlets of the models on the vertical axis of Figure 14. We have considered the mass of only those particles that entered into the calculation of the physical properties shown in § 4. In the horizontal axis, we show Models from 1 to 4, ordering the models in the following way: U 5, U 9, U 11 and U 13, respectively.

The mass of the cloudlets for Models U 5, U 9 and U 13 are shown in the top left-hand panel, so that the mass ranges from log (M_f /M_⊙) =3.6 to 3.8. Two cloudlets of Model U 11 have the largest masses, around log (M_f /M_⊙) = 4.6.

In contrast to what we have observed when compared to other physical properties of Models U with Ub, the behavior is different; in the case of the mass of the cloudlets, without any trend.

The mass of the cloudlets for Model Ur is shown in the top right-hand panel of Figure 14. We see that Models Ur5, Ur9 and Ur11 all have similar cloudlet masses, which are around log (M_f /M_⊙) = 4.5. Meanwhile, the cloudlet mass for Model Ur13 is a little smaller than the mass of the other models, ≈ log (M_f /M_⊙) = 4.4.

Given that Models Urb have the highest in-fall radial velocity of the all the models, as can be seen in 12, then their mass accretion rate must be the highest too; for this, the mass enclosed is systematically higher than in the other models, although there is not much difference in the mass values observed in the panels of Figure 14 when compared to the big difference in the infall radial velocity.

It must be emphasized that the mass scale shown in Figure 14 is in agreement with the mass observed for an open cluster of stars, which is around 10⁴ M_⊙, while the mass scale of a globular cluster of stars is around 10⁵ M_⊙, see ^{Kumai et al. (1993)}.

The results of the collapse of a gas core (e.g., in the so called “standard isothermal simulation”) are identified as protostars (^{Boss (1995)}, ^{Boss et al. (2000)}, ^{Burkert and Alves (2009)} and ^{Arreaga (2007)}). In the same sense, the gas structures obtained from the simulations of the present paper and whose properties are shown in § 4 can be called a proto-cluster of protostars. As is well-known for simulations of the collapse of a gas core, the mass of the proto-stars depends on the mass of the parent cloud, see for instance ^{Arreaga (2016)}. The same is expected to be true for proto-clusters.

Finally, it must be recalled that the results displayed in Figure 13 and Figure 14 are taken when the collapse is still ongoing, so that the peak density is around log (ρ_max/ρ₀) ≈ 5. As we have seen in § 3.1, the final state of the collapse reaches a peak density around log (ρ_max/ρ₀) ≈ 8.

Fig. 13 The ratio of the thermal energy to the gravitational energy α versus the ratio of the kinetic energy to the gravitational energy β, of the cloudlets defined in § 4.3, at the same time as the snapshots shown in Figure 4, Figure 5, Figure 6 and Figure 7, respectively. (top left-hand) U with a low level of turbulence; (top right-hand) Ur with a low azimuthal velocity; (bottom left-hand) Ub with a high level of turbulence and (bottom right-hand) Urb with a high azimuthal velocity. The color figure can be viewed online. 

Fig. 14 The log of the mass of the cloudlets defined in § 4.3, at the same time as the snapshots shown in Figure 4, Figure 5, Figure 6 and Figure 7, respectively. (Top left-hand) U with a low level of turbulence; (top right-hand) Ur with a small azimuthal velocity; (bottom left-hand) Ub with a high level of turbulence and (bottom right-hand) Urb with a large azimuthal velocity. The color figure can be viewed online. 

5.1. The Collapse of the Isolated Turbulent Cloud

We mentioned in § 1 and § 2.2 that the initial conditions of the isolated cloud are chosen to favor its gravitational collapse. The curve of the peak density for the isolated cloud develops a small peak at a time smaller than t/t_ff = 0.1. This increase of density happens because of the multitude of gas lumps formed by the collisions between gas particles that occur simultaneously throughout the cloud. This density peak does not appear for low levels of initial kinetic energies, which is measured by the β ratio defined in equation 5; for high values of the β ratio, this early peak is quite noticeable.

In contrast, the time required by the isolated cloud to reach the highest density values, for instance log (ρ_max/ρ₀) ≈ 6, does not depend significantly on the level of the initial energy, at least for a wide range of the initial β ratio. This is due to the fact that almost all of the kinetic energy available is equally dissipated by the random collision of particles, as described in the previous paragraph. The time required for the cloud to reach its highest peak density is around t/t_ff = 2.5, which is of the same order of time that can be seen in Figure 3 for the collapse of the collision models.

In the time interval between 0.1 < t/t_ff < 2.0, a relaxation of the small gas lumps occurs throughout the cloud, so that the peak density curve decreases quickly. From there, it increases very slowly, up to times t/t_ff >2.0, at which the final collapse takes place very quickly.

From the point of view of the column density plots, the occurrence of collisions between gas particles, as a consequence of the turbulent velocity field implemented initially, is seen as a random formation of many over-dense lumps of gas, which are homogeneously distributed across the entire cloud volume, see the left-hand panel of Figure 15. Later, when the initial kinetic energy of the cloud is dissipated, the cloud reaches a physical state similar to a free-fall collapse, which is seen as a clear tendency to a global collapse towards its central region. However, at the final evolution stage that could be followed in this paper, the mass accretion with spherical symmetry is lost, so that a central dense filamentary structure forms that is highly anisotropic and with a high possibility of fragmenting, see the right-hand panel of Figure 15. The behavior described in this section is paradigmatic of turbulence.

Fig. 15 Column density plots of the isolated turbulent cloud, for a thin slice of gas parallel to the x-y plane. The unit of length is one parsec. The plots are shown in panels as follows: (left-hand) at time t/t_ff = 0.02 and peak density log (ρ_max/ρ₀) = 0.53; (right-hand) at time t/t_ff = 2.5 and peak density log (ρ_max/ρ₀) = 8.0. The color figure can be viewed online. 

5.2. Does the Turbulence make a Difference in the Collision Simulations?

The occurrence of the collision between the subclouds induced by the translation velocity v_L : v_R prevents the gas particles from forming small lumps of gas throughout the cloud by means of early random collisions; as explained in § 5.1.

There is over-dense gas in the contact region between the colliding sub-clouds. This over-density accelerates the collapse of the remaining gas of the cloud, so that the turbulence does not have time enough to get relaxed by dissipation of the kinetic energy. For this reason, the turbulence does not play a fundamental role in the outcome of the simulations, such as U and Ub. In fact, if one turns off the turbulence and keeps only the collision process of the sub-clouds in these models, then the results are basically the same.

5.3. Does the Level of Turbulence make a Difference in the Collision Simulations?

We recall that the level of turbulence in the simulations can be modified by introducing an arbitrary multiplicative constant in equation 1. As we mentioned in § 2.2.1, there is interest in considering models of turbulent clouds with extreme initial kinetic energy in addition to those clouds with low-level turbulence, which are more favored statistically. Consequently, we have studied the effect of the level of turbulence on the simulations (i.e., Models Ub), as can be seen in Table 1.

The average Mach velocity ℳ_p of the gas particles for the low-level of turbulence Models U is around ℳ_p ≈ 2.9. For the radial component of the velocity (calculated with respect to the origin of coordinates of the simulation box), the average Mach number is ℳ_r ≈ 0.13. For the tangential component of the velocity, the average Mach number is ℳ_t ≈ 0.01. In the meantime, the translational velocities (around 15 km/s) given to the particles that are to collide are of order ℳ_c ≈ 6.6. Then, for Models U, ℳ_c >> ℳ_p,r,t.

For the high-level of turbulence Models Ub, we have an average ℳ_p ≈ 25. Despite this significant increase of the magnitude of the velocity, the average radial and tangential components do not change appreciably with respect to those of Models U; that is, ℳ_r ≈ −0.12 and ℳ_t ≈ 0.01. In contrast, for Models Ub we have ℳ_c << ℳ_p.

In spite of the opposite features in the relation of Mach numbers for Models U and Ub with respect to the translational velocity, the outcome of the simulations U and Ub do not show any significant difference with respect to the final configuration of Models U 5, U 9 and U 13. The only differences that can be observed are: (i) that the double bridge of gas formed in Model U 11 becomes only one bridge in Model U 11b; and (ii) that the density peak curves, shown in Figure 3, for the Ub models are displaced to the right-hand side at large free-fall times, so that the collapse takes a little longer than for Models U.

5.4. How Useful is the Approximation of an Azimuthal Velocity to Mimic Tidal Forces in the Collision Simulations?

There is an obvious problem with the approximation of a velocity instead of a tidal force, as described in § 2.3, which is that the azimuthal velocity entered only once in the simulations, as an initial condition of the gas particles. Obviously, this is a severe limitation of the model, similar to that of the turbulence, which is not replenished continually during a simulation. Therefore, the effect of the tidal interaction must be activated during all the simulation time.

However, we observe in Figure 6 that the immediate effect of the azimuthal velocity on the simulated cloud is a strong tendency for the gas to be accumulated quickly at the cloud’s center. In this case, if the azimuthal velocity terms given in equation 10 were implemented at every time step of the simulation, then to model more appropriately the tidal force over all the simulation time one would expect this tendency to accelerate the central collapse of the cloud.

It should be emphasized that the previous statement is based on the results of a very naive model, in which the only information about the massive center exerting a gravitational force on the cloud is by means of the circular velocity, which is given by √2 GM (R)/R, as described in § 2.3.

According to § 2.3, the approximation of the azimuthal velocity is valid as long as the ratio between the cloud radius to the distance to the gravitational center is quite small. A way to check the applicability of this approximation is obviously to make the calculation without the approximation. However, this is not an immediate calculation. The main difficulty is the difference in length and mass scales when considering a small cloud (with very few parsecs of radius) near a massive object (probably with a scale of kpc in radius, separated from the cloud by several hundreds of pc, or even a kpc, and whose mass can quite greater than that of the cloud), both of which must have evolved together in the same simulation code.

For instance, ^{Gnedin (2003)} resorts to a resimulation technique, so that a low-resolution simulation of the massive object (e.g., a central dwarf galaxy) is first carried out to obtain an approximate gravitational potential. This is then used in a second high-resolution simulation of the cloud, in which this potential is taken into account as an external timevarying field on the gas particles. However, applying this technique to the problem presented in this work would require a future paper.

5.5. A Brief Review of the Literature on this Subject

Many papers have simulated isolated clouds and followed their collisions. However, simulations of clouds under the influence of an external gravitational potential are limited in number.

Let us now mention briefly some results of more accurate calculation methods of the tidal effects on clouds, which is a subject that has a long history. For instance, ^{Sigalotti and Klapp (1992)} used a time-varying gravitational potential to calculate the equal-sized cloud-cloud tidal interaction of clouds that are in an elliptic orbit, and reported configuration transformations on the clouds in their course to collapse.

More recently, ^{Longmore et al. (2013)} proposed that the collapse of the Brick is a progenitor of a star cluster, whose collapse was triggered as a consequence of the tidal compression exerted by Sgr B2 during the most recent peri-center passage of the Brick.

^{Kruijssen et al. (2015)} determined a realistic orbit of a dense gas streams in the CMZ. In a subsequent paper, ^{Kruijssen et al. (2019)} calculated the tidal interaction of the galactic center on the orbit followed by the dense gas streams of the CMZ. They found that the tidal interaction acting upon the clouds makes a compression on the vertical direction, which causes the clouds to become pancake-like structures.

Later, ^{Dale et al. (2019)} simulated the evolution of turbulent clouds in orbit at the CMZ. The authors assumed a similar magnitude of the kinetic energy to the gravitational energy and found that the clouds collapse rapidly. This paper is a mature way of simulating tidal force in cloud dynamic evolutions by introducing the tidal force potential.

5.6. Applicability of these Simulations to Represent the Evolved Clouds of the CMZ

In view of § 5.2 and § 5.3, the collision process (and its parameters) is clearly the dominant physical mechanism in shaping the appearance of the cloud in the simulation outcome. It is possible that this collision of sub-clouds, with the cloud’s self-gravity, is the dominant process of the cloud evolution, even over the tidal interaction with the massive center, and above all, for the small scale of the circular velocity that is induced, as compared to the magnitude of the other velocities involved, which are the turbulent and the translational velocities, see § 2.3.

For Models U 5, U 9, and U 13, the geometry of the resulting configuration can be well characterized by defining a center and a radius of a cloudlet. For Model U 11, this spherical structure does not make sense, as can be seen in Figures 4 and 8. The outcome of Model U 11 is a complex, structured molecular gas cloud that exhibits an interconnected network of components. This is the only model that can be compared or approximated to the cloud configuration called the Brick.

In fact, for the Brick, a shell-like structure with radius of 1.3 pc has been revealed from observations in the integrated intensity map of SO. For instance, see Figures 1 and 3 of ^{Longmore et al. (2012)}; Figures 1 and 3 of ^{Kauffmann et al. (2013)}; Figures 1 and 2 of ^{Higuchi et al. (2014)}. ^{Kruijssen et al. (2019)} presented three panels in their Figure 6, to compare the results of ALMA observations of the Brick to a synthetic observation obtained from a numerical simulation. A complex gas structure can be seen in these panels, in which the gas condenses in a persistent diagonal direction with many twisted and bending filaments connected in a messy way. Model U 11 of this paper clearly shows a similar diagonal direction of the dense gas.

It may seem that there is a huge problem with Models Ur, given that all of the different structures obtained as a result of the collision process in Model U are destroyed because of the azimuthal velocity of Models Ur. However, the configuration obtained in Models Ur can be well recognized as the final outcome of the formation process of a YMC, which is observed to be a strong central condensation of gas, with an enclosed mass of stars of about 10⁴ M⊙ with a size of one pc; see for instance ^{Rathborne et al. (2015)}. The Arches cloud is an example of this kind of observed configuration; see ^{Portegies al. (2010)}.