Abstract
By using the particle-based code Gadget2, we follow the evolution of a gas giant molecular cloud, in which a set of gas particles representing the wind are created by a Monte Carlo scheme and suddenly move outwards from the cloud’s center. The particles representing the gas cloud initially have a velocity according to a turbulent spectrum built in a Fourier space of 643 grid elements. The level of turbulence and the temperature of the cloud are both adjusted so that a gravitational collapse of the cloud is initially induced. All the winds are activated in a very early stage of evolution of the cloud. We consider only two kinds of winds, namely, one with spherical symmetry and the second one of a bipolar collimated jet. In order to assess the dynamical change in the cloud due to interactions with the winds, we show isovelocity and isodensity plots for all our simulations. We also report on the accretion centers detected at the last simulation time available for each model.
1. Introduction
Stars are born in large gas structures made of molecular hydrogen. According to [1] a cloud is a gas structure with mass and size around and 2–15 pc, respectively.
The physical process by which the molecular gas is transformed from a gas structure into stars is mainly gravitational collapse, whose main effect on the cloud is a reduction in size with the corresponding increase in density at various points of the cloud. These small overdensities within the larger gas structure are defined as cores, which have a typical mass and size around and pc, respectively, [1].
The basic idea of prompt fragmentation is that during gravitational collapse a molecular cloud may spontaneously break into several cores. At some point during the transformation process from clouds to cores and then to stars, the densest gas structures settle down in more or less dynamically stable gas objects called protostars. If the cores do not undergo further subfragmentation upon collapsing to higher densities, they will become protostars.
To achieve a better understanding of the huge scale changes due to this transformation process, we may recall that the number density of a typical star is around molecules per . Meanwhile, the number density of a typical typical cloud structure ranges around molecules per and that of a typical protostar ranges around molecules per .
Many other phenomena occurring in the interstellar medium can have a great influence on the evolution of the cloud, among others: (i) the highly ionized gas ejected by supernovae explosions; (ii) the bipolar collimated winds ejected by massive protostars; (iii) rapidly expanding regions which hit the slower and less dense gas structures.
For instance, in [2] a set of SPH simulations were conducted to study star formation triggered by an expanding region within a spherical gas structure of uniform density that was not globally collapsing. The expanding shock plays the role of a snowplow, which sweeps out the gas of the surrounding medium. The density of the swept out gas increases as a consequence of this agglomeration process, and a gravitational collapse may then be locally initiated. Small gas overdensities can be formed in this way, which may achieve the protostar stage. Because these two processes are complementary in forming protostars, this star formation scenario is known as the collect and collapse model. The one-dimensional model of freely expanding wind into an uniform ambient medium was analitically and numerically solved some time ago; see [3] and references there in.
In this paper, we investigate the change in the dynamical configuration of a typical turbulent cloud when a wind of particles is outwardly ejected from the central region of the cloud. The most important differences between the present paper and [2] are the turbulent nature of our cloud and its induced global tendency to gravitational collapse. Turbulence makes a big change in the spatial distribution as the cloud becomes filamentary and flocculent. The tendency to collapse favors the condensation of gas towards the central region of the cloud. Thus, a gas wind flowing through this collapsing and turbulent gas cloud is a very difficult hydrodynamical phenomenon that can only be solved by numerical simulations. It is shown that the presence of winds (i) keeps the collapsing nature of the turbulent cloud, but they cause a delay in the runway collapse; (ii) changes the number and location of accretion centers detected in the densest central region of the cloud.
Furthermore, in order to study the effects of protostellar outflows on the turbulence of a star forming region, in [4, 5], magneto-hydrodynamics simulations (MHD) were conducted with a mesh-based code which implements the adaptive mesh refinement technique (AMR). As we use here the fully parallelized Gadget2 code which implements the SPH technique, this will make it possible to carry out an interesting comparison between the two computational techniques.
2. The Physical System
In this section we briefly describe the physics of the cloud and the winds which will be considered in the following sections.
2.1. The Cloud
We consider here a typical spherical cloud with a radius pc and mass . Initially, it has a radially uniform density distribution with an average density given by gr , which is equivalent to a number density molecules for molecular hydrogen with molecular mass gr/mol. The size and mass of this cloud are chosen here to be typical in statistical sense, in accordance with [1].
The free fall time is defined as the time needed for an external particle to reach the center of the cloud when gravity is the only force pulling the particle. In this idealized gravitational collapse, we havewhere is Newton’s universal gravitational constant. For our cloud, we have years.
For a spherical cloud, the approximate total gravitational potential energy is . The average total thermal energy (kinetic plus potential interaction terms of the molecules) is , where is the Boltzmann constant, is the equilibrium temperature, is the total number of molecules in the gas, and is the speed of sound; see Section 3.5 for a more precise definition.
The kinetic energy can be estimated by , where is the average translational velocity of the cloud. In order to have both energies of the same order of magnitude, , the gas elements of the cloud must attain average velocities within the rangeor km/sec, for a speed of sound given byso that the corresponding temperature associated with the cloud is K.
It is possible to define the crossing time by means ofwhich sets a time scale for a sound wave to travel across the cloud. To make the crossing time comparable in magnitude to the free fall time of (1), the front wave must have velocities around or km/sec, which are velocities a little bit slower than the ones estimated above; see (2). In this paper, we will treat propagation velocities of gas particles ranging around Mach.
2.2. The Wind
The dynamical characteristics of the wind strongly depend on its type of source. All stars eject winds of the first kind, which are driven by stellar radiation. For instance, in cool stars, like the ones observed in the AGB (asymtotic giant branch) of the Galaxy, the winds cause a mass loss in the range yr whereas the terminal wind velocities are around km/sec. In OB stars, the mass loss ranges over /yr and the terminal wind velocities can go up to thousands of km/sec.
Supernovae dump around joules of thermal and kinetic energy into the interstellar medium. But there are several types of supernovas, so that the mass losses and terminal velocities are very different; see [6]. For example, for a supernova whose progenitor was a He star, the mass loss and terminal velocities are within the ranges /yr and km/sec, respectively. When the progenitor was a RSG star, their values range over /yr and km/sec, respectively.
It seems that all protostars eject highly collimated jets of gas during their formation process by gravitational collapse. The origin of these jets is still unclear, but it is very likely that the accretion disk and magnetic field around protostars play a crucial role in determining the velocities and degree of collimation of the jets. For the molecular winds associated with protostarts of Class 0 and Class 1, the characteristic velocities are around 20 km/sec. However, for optical jets of highly ionized gas, the typical jet velocities are a few hundred km/sec; see [7] and the references therein.
3. The Computational Methods and the Models
In this section we briefly describe the way we set up the physical system outlined in Section 2 in computational terms as well as the models we consider in this paper.
3.1. The Initial Configuration of Particles
We set million SPH particles for representing the gas cloud. By means of a rectangular mesh we make the partition of the simulation volume in small elements each with a volume ; at the center of each volume we place a particle (the th, say), with a mass determined by its location according to the density profile being considered, that is, with . Next, we displace each particle from its location by a distance on the order of in a random spatial direction.
As stated earlier, in this paper we only consider a uniform density cloud, for which , for all the simulations (see Section 2.1). Therefore, all the cloud particles have the same mass.
3.2. The Initial Turbulent Velocity of Particles
To generate the turbulent velocity spectrum for the cloud particles we follow a procedure based on the papers [8, 9]. We set a second mesh with the size of each element given by , , and . In Fourier space, the partition is , , and . Each Fourier mode has the components , , and , where the indices , , and take values in , , and , respectively. The wave number magnitude is , and so and . The Fourier wave can be equally described by a wave length , and then we see that and .
The components of the particle velocity arewhere the spectral index was fixed at and thus we have . The vector whose components are denoted by takes values obeying a Rayleigh distribution. The wave phase vector, , given by takes random values on the interval . The components of the vector are calculated by means of , where is a random number in . is a fixed parameter with value .
3.3. The Set-Up of the Particle Wind
We generate the initial wind particles by a traditional Monte Carlo scheme, in which particles are initially located randomly in the central volume of the cloud. Let , , and be random uniform variables taking real values within the interval . According to the fundamental probability conservation law for a system with spherical symmetry, we must have . By means of integration we obtain that the spherical coordinates of the particles are related to the uniform random variables by the following equations:where is the total mass contained in the central cloud region of radius and its mass density. It must be noticed that the first relation should be integrated numerically to obtain once has taken an allowed uniform random value.
In order to assign a velocity vector to each wind particle, let us consider the equation of mass conservation for a set of particles moving radially outwards; that is,In order to determine the wind velocities according to (7), where the radius ranges in the interval , the mass loss and the mass density will be considered as free parameters in the simulations; see Section 3.8. As the velocity magnitude diverges for particles around , we must then set up a cut velocity value, . Therefore, the velocity vector for the wind particle is given by , where is the speed of the particle obtained by (7) and is the unit vector pointing radially outwards.
Of course, there are other possibilities, which will be considered elsewhere: one is to fix the radial density and/or the velocity profile in order to obtain the mass loss as a result. For example, for modeling an expanding region, the authors of [2] proposed a complicated velocity function, which gives a constant expansion velocity at the last stages of time evolution, so that the average velocity of the shocked shell considered by [2] is or km/sec.
3.4. Resolution and Thermodynamical Considerations
Following [10, 11], in order to avoid artificial fragmentation, the code must fulfil certain resolution criteria, imposed on the Jeans wavelength , which is given bywhere is the instantaneous speed of sound and is the local density. To obtain a more useful form for a particle based code, the Jeans wavelength is transformed into the Jeans mass given byIn this paper, the values of the density and speed of sound are updated according to the following equation of state:as proposed by [12], where , and for the critical density we assume the value .
For the turbulent cloud under consideration, we have g, where we take .
In this paper, the mass of an SPH particle is g, so that and therefore the Jeans resolution requirement is very easily satisfied.
3.5. Initial Energies
Following [13], the dynamical properties of the initial distribution of gas are usually characterized by , the ratio of the thermal energy to the gravitational energy, and , the ratio of the rotational energy to the gravitational energy.
In a particle based code, we approximate the thermal energy of the cloud by calculating the sum over all the particles described in Section 3.1; that is,where and are the pressure and density associated with particle by means of the equation of state given in (10). In a similar way, the approximate potential energy iswhere is the gravitational potential of particle . For all the simulations considered in this paper, the values of the speed of sound (see (3)) and the level of turbulence are chosen so that the energy ratios have the numerical values
3.6. The Evolution Code
We carry out the time evolution of the initial distribution of particles with the fully parallel code, which is described in detail by [14]. is based on the method for computing the gravitational forces and on the standard SPH method for solving the Euler equations of hydrodynamics. incorporates the following standard features: (i) each particle has its own smoothing length ; (ii) the particles are also allowed to have individual gravitational softening lengths , whose values are adjusted so that for every time step is of order unity. fixes the value of for each time-step using the minimum value of the smoothing length of all particles; that is, if for , then .
The code has an implementation of a Monaghan-Balsara form for the artificial viscosity; see [15, 16]. The strength of the viscosity is regulated by setting the parameters and ; see in [14]. Here we fix the Courant factor to .
3.7. The Detection of Accretion Centers
Let us now describe the code implemented to detect accretion centers in the models. Any gas particle with density higher than is a candidate to be an accretion center. We localize all the candidate particles for a given snapshot at time .
We then test the separation between candidate particles: if there is one candidate with no other candidates closer than , then this particle is identified as an accretion center at time . We define as the neighbor radius for an accretion center, given by , where is the minimum value of the smoothing length of all particles; that is, if with being the smoothing length of each particle , for in this way determines a set of particles which are within the sphere having this radius and whose center is the accretion center itself. All those particles will contribute to the mass and momentum of the accretion center.
We carry out the detection of accretion centers for all the models, using a fixed value of given by , as was done by [5].
The code described here is a first step towards the full sink-particle technique implementation already in use worldwide. Our code is still incomplete as we have not yet included all the particle tests suggested in [17, 18].
3.8. The Models
The initial configuration of particles has many parameters and their values determine the models to consider as follows. To sum up the meaning of these parameters in order of importance, we make a list in Table 1, without any reference to some specific physical system. In the first line, determines the time when the wind appears within the turbulent cloud. In the second and third lines, and determine the velocity by (7). Next, the radius determines the initial spatial extension of the wind particles. () is the total number (the total mass) of wind particles within . is the peak wind velocity allowed at the initial simulation time to avoid divergences in the velocity given by (7). Finally, we need to define two angles in order to characterize the collimated jets: the polar angular aperture and the azimuthal angular aperture , centered around the -axis with , and the equator with , respectively.
In this paper, we fix and km/sec for all the simulations while and for the collimated jet simulations. The choice of a wind density value must not be so different to the cloud uniform value , to avoid discontinuity problems between the gas cloud and the gas wind when Gadget2 calculates the tree. Besides, it should be noticed that values of have been chosen to represent a core (with typical values of size and mass) within the cloud.
In Table 2 we present the values used to define the models calculated below as well as some results described in the coming sections. In this Table and throughout the paper, we use the following labels: the label “T” to indicate the turbulent gas cloud when the winds appear within the cloud, the label “TW” to refer to the spherically symmetric case, and the label “TWC” for the bipolar collimated jet.
The first six columns of Table 2 were already defined in the preceding paragraphs. In the seventh column of Table 2, we show the total mass swept out by the winds outside by the end of the simulation. In column 8 of Table 2, we show the velocities attained by those particles located far outside the initial cloud. In columns 9 and 10 of Table 2, we show the maximum evolution time and peak density reached in the simulations.
4. Results
To present the results of our simulations, we consider a slice of particles around the equatorial plane of the cloud; with these particles (around 10000) we make plots containing two complementary panels: one to show colored regions of isodensity and the other one to show the velocity field of the particles. Later, we present plots of the velocity distributions and the radial velocity profile, for which we use all the simulation particles. Finally, in order to make a comparison of the changes produced in the dynamical state of the central region of the cloud, we present and plots with the accretion centers detected by locating the densest particles of each simulation.
4.1. The Evolution of the Turbulent Cloud
Because of the random nature of the initial velocity distribution, a huge number of particle collisions occur across the cloud. In the first two panels of Figure 1 one can already notice one of the main characteristics of turbulence: the appearance of a filamentary and flocculent structure across the cloud, a spatial structure similar to “Swiss cheese.”
The particle collisions dissipate very quickly much of the kinetic energy of the cloud, so that the cloud goes to a physical state similar to a free fall collapse; then the gravitational force pulls the cloud mass to a global collapse towards its central region, as can be seen in the last panels of Figure 1.
As one can see in the top panel of Figure 2, very few particles can attain velocities much higher than those average velocities provided by turbulence alone. But almost all the simulation particles have velocities . In the middle panel of Figure 2 we show the radial profile of the velocity, which clearly indicates that almost all particles have the same average velocity across the cloud, which is around Mach. At the time , the fully flocculent cloud reaches a the peak in its density curve, as we show in the bottom panel of Figure 2.
(a)
(b)
(c)
We emphasize now a very important fact occurring at the outer regions of the turbulent cloud. Near the cloud boundary, the interpolating kernel of the Gadget2 code introduces small errors in the calculation of the particle density that makes the boundary density lower than it should be. So to speak, the mass of each SPH particle weighs less for particles located near the boundary. There are methods for correcting this density deficiency problem; see, for example, [19, 20]. Unfortunately, such corrections have not been incorporated so far into the public Gadget2 code. Furthermore, the turbulent cloud is not in hydrodynamic equilibrium and there is no external pressure acting upon the cloud, and then the outermost particles have nonequilibrated thermal pressure. All these features act together in such a way that the outer cloud particles expand outwards.
So, we have to keep in mind this expansion effect for the problem at hand, as we shall quantify the mass of the cloud swept out by the winds; see, for instance, Figure 3 and the next section.
(a)
(b)
(c)
4.2. The Effects of Winds in the Evolution of the Turbulent Cloud
Let us again consider the top and middle panels of Figure 2, which show that the wind particles near the cloud’s center can attain velocities much higher than . It must be kept in mind that these curves are smoothed; thus they appear to steep even for much greater than .
Because the wind particles are located at the center of the cloud by a Monte Carlo method, the particle density there significantly increases at time , when the winds are activated, as can be seen in the bottom panel of Figure 2. It should also be noticed by looking at this panel that the global collapse of the cloud does not change even when the winds particles are present, as the peak density curve of each run goes to higher values by the end of the simulation time.
We emphasize that for all the wind models the gravitational collapse of the cloud is slightly slower as their peak the curves begin to slope upwards a little later the middle up to the final stages of their evolution time. The reason to claim this delay is that the peak density values reached by the wind models are always lower than the turbulent cloud compared at a similar evolution time. For instance, let us mention that the turbulent cloud alone needs to achieve a central peak density of while the models TW1 and TWC3 have peak densities several orders of magnitude smaller at a evolution time slightly greather than . The model TWC1 needs to reach a peak density of , as can be seen in the last two columns of Table 2.
The isodensity and velocity plots for all the simulations are shown in Figures 4, 5, 6, 7, and 8. It is recommended to make zooming in the panels of these figures in order to notice in all the model the occurrence of a common transformation process: from a very disordered velocity distribution to a very clear tendency to condensate toward the central cloud region. There are still some differences in the final morphology and in the accretion centers detected in the central region of the cloud.
Let us first consider the models with spherically symmetric winds, TW1 and TW2. The wind particles for model TW1 push the cloud particles outwards. Thus, by the time , we see a void created in the central region of the cloud. However, gravity and viscosity act together in such a way that the particles quickly fill the void, as can be seen in the last panel of Figure 4. Instead, the wind particles for model TW2 diffuse through the cloud without creating a void. As there is no much more mass swept out than that of the turbulent cloud alone, it is possible that all the wind particles reach the outer cloud boundary, but they do not push gas particles outwards.
Let us now consider the bipolar collimated jet models TWC1 and TWC2 illustrated in Figures 6 and 7. Essentially, we see either collision (for model TWC1) or diffusion (for model TWC2) of the wind particles taking place in the directions in which the collimated jets are initially emitted.
For model TWC1, the lack of wind particles in the perpendicular directions of the emitted jet produces a slightly elongated void, as can be seen in the second panel of Figure 6. As expected, the cloud mass located in the perpendicular directions to the emitted jets is strongly attracted by the gravitational force of cloud mass located in the other side. When the central cloud region is filled again with mass, it adopts a slightly elongated condensation along the direction of the emitted jets. With the purpose of emphasizing this effect, we include another collimated model in which we increase the initial mass, the model TWC3, which is illustrated in Figure 8.
We emphasize that the fraction of gas swept out by the winds is quite significant: from nothing to dozens of solar masses moved even far beyond the cloud radius , as can be seen in the seventh column of Table 2. We also report in the eighth column the expansion velocities reached during the time we follow the simulations. The most evolved simulation we have is for model TWC1. This fact alone explains why it has the highest values reported in Table 2 for the the masses and .
As expected, all the wind models show a strong tendency to form accretion centers in the central region of the cloud. It must be noticed that the great differences appear when we compare the number and location of the accretion centers obtained for the turbulent cloud and for the wind models, as can be seen in Figures 9 and 10. In these plots, it is possible to see the spatial distribution of the accretion centers: every point represents a SPH particle, and then we show all the particles captured by the formed accretion centers.
(a)
(b)
(c)
(d)
(e)
(f)
5. Discussion
In this paper we carry out a fully 3D set of numerical hydrodynamical simulations within the framework of the technique, aimed at following the interaction of winds with a turbulent cloud. We comment now about three technical limitations we ougth to keep in mind when looking at the results reported here.
Firstly, the Monte Carlo method can generate couples of very close particles. These couples should be removed because the Gadget2 code gets into troubles when the construction of the tree is iniciatted. We romove close particles by running a code that searches for those particles which are closer than an allowed minimun distance. Thus more particles than really needed are removed by this procedure. The net effect is that there is no control on the total number of wind particles to be evolved in a simulation. This explains why the number of particles entering in the wind models are so different, as can be seen in Table 2.
Secondly, the Gadget2 code used here was written in the so-called density-entropy formulation, for which the adaptive smoothing length of each particle is determined by creating a kernel volume containing a constant mass, in which a number of neighbour particles is kept roughly constant for each particle . For this reason, we can not have a wind density so different from the cloud density as the smoothing lengths may become too large for those less clustered wind particles; see [21].
Thirdly, it must be taken into consideration that (i) we have used a low density threshold for a particle to become an accretion center and (ii) there is no restriction on the minimum number of particles that must be captured for a particle to really become an accretion center: so we have seen that there are a large number of accretion centers having just a very few captured particles. For these two reasons we must expect overproduction of accretion centers in some of our simulations. However, as we have applied the same code to all the models, a comparison is possible. Finally, we recall that accretion centers are made of gas particles which can collide or merge. In view of this, the future of the accretion centers is still unclear and requires taking the calculations further in time.
6. Concluding Remarks
A star formation scenario based only on the collapse of turbulent gas structures provides a very highly efficient transformation of the gas into protostars. However, this point is in contradiction with observations.
Because of this, another scenario must be considered, or at least a theoretical complement to the turbulent model is needed. As we have shown in this paper, winds must be considered an additional ingredient to complement the turbulent model with the hope that these models TW and TWC can alleviate the problem mentioned above as they make a delay in the runaway collapse of the cloud.
We have obtained a physical system showing simultaneously in-fall and out-flow motions of particles. We emphasize that these simultaneous motions are observed in clusters like NGC 1333 and NGC 2264; see [22, 23], respectively.
Another interesting objet is IRAS 16547-4247. It has a very dense single central region surrounded by a cloud of about . A couple of wind jets have been detected recently in this star forming region and they are belived to be generated by a massive protostar; see [24]. The model studied here can be modified in order to take into account the new findings.
This paper must be considered a first step towards a more complete study of winds within a collapsing turbulent cloud. A forthcoming publication will treat a different mass and density for the wind particles and the different directions of the emitted jets to model system like the IRAS source 16547-4247.
Conflict of Interests
The authors declare that there is no conflict of interests regardng the publication of this paper.
Acknowledgments
The authors would like to thank ACARUS-UNISON for the use of their computing facilities in the making of this paper.