Abstract
This work presents three-dimensional hydrodynamical simulations with the fully parallel GAGDET2 code, to model a rotating source that emits wind in order to study the subsequent dynamics of the wind in three independent scenarios. In the first scenario we consider several models of the wind source, which is characterized by a rotation velocity and an escape velocity , so that the models have a radially outward wind velocity magnitude given by 1, 2, 4, 6, and 8 times . In the second scenario, we study the interaction of winds emitted from a binary system in two kinds of models: one in which the source remains during the wind emission and a second one in which all the source itself becomes wind. In the third scenario we consider the interaction of a rotating source that emits wind within a collapsing and rotating core. In this scenario we consider only wind models of the second kind built over a new initial radial mesh, such that the angular velocity of the wind is 1, 100, and 1000 times the angular velocity of the core .
1. Introduction
The star formation process is initiated when random density fluctuations trigger gravitational collapse in large gas structures (called clumps) made of molecular hydrogen. According to [1], the first dynamically important overdensities formed by this process (called cores) have a typical mass of 0.5–5 and size of 0.03–0.2 pc, respectively. These overdensities can undergo further gravitational collapse to reach higher densities, until they become protostars. For instance, the number densities of a typical gas clump is around molecules per cm−3 while that of a typical core is around molecules per cm−3 and that of a typical protostar varies around molecules per cm−3.
Over the last two decades, a large number of papers have focused on numerical studies of the gravitational collapse of isolated cores; see the review given in [2] and references therein. The outcome of many of these numerical experiments was the formation of a set of low-mass protostars, sometimes grouped in binary systems; see the review in [3]. Furthermore, all these results have been validated by using different codes, such that a very interesting comparison of results obtained in the recent history of collapse calculations of cores with codes based on different numerical techniques is given by [4].
Recent technical improvements on the observational side have made it possible to observe a few protostars in their Class 0 evolution stage. It is thus known that all protostars eject highly collimated jets of gas during their formation process by gravitational collapse. In fact, protostars of Class 0 and Class 1 emit winds with characteristic velocities around 20 km/s; see [5] and the references therein.
The winds emitted by an already formed star can also have an important feedback effect on the surrounding medium up to large scales. In this case, the ejected winds are driven by stellar radiation. In cool stars, the winds cause a mass loss in the range of , whereas the terminal wind velocities are around km/s. In OB stars, the mass loss ranges over /yr and the terminal wind velocities can be in thousands of km/s. In fact, observational evidence of large scale gas outflows has been provided by studying the systems M17 and RCW49i, in which a stellar bow shock around O-stars is seen outside the central cluster; see [6].
On the theoretical side, in [7] a set of hydrodynamical simulations was conducted considering three types of stars ejecting winds into a turbulent and heterogeneous clump. Although these simulations did not include gravity or radiative losses, they found a complex interaction in which the winds carved channels in the clump; see also [8, 9].
Recently, many papers have been published with regard to the study of winds and their effects on the different stages of the star formation process. For instance, [10, 11] studied the effects of protostellar outflows on the turbulence of a star forming region. These simulations were done with the AMR technique, and like other simulations also based on grids the information of the wind is entirely codified in only one or two cells, while the remaining grid elements host the information of the turbulent cloud. This situation is similar for simulations done with a particle-based technique: either the wind source is skipped or the source information is contained in very few particles.
In [12], a toy model of a compact and rotating star is presented in which an attractive force acts upon the SPH particles in order to produce a static gas structure; see also [13]. In this paper we present a hydrodynamical model of a rotating source of wind without any compensating force included. Despite the fact that the dynamical characteristics of the emitted wind strongly depend on its type of source, we want to keep the source unchosen (as much as possible) in building three scenarios where both the source and the emitted winds are relevant. For instance, in the first scenario we consider several wind models, in which the wind velocity magnitude is given by 1, 2, 4, 6, and 8 times the rotation velocity of the source . We find that the rotation of the source induces the formation of a mass shell through rear-end collisions among those wind particles not gravitationally escaped from the mass of the source.
In [14] a hydrodynamical model of a binary merger is presented in order to explain the three rings observed around the system SN 1987A. Furthermore, in [15], hydrodynamical models of rotating stars are presented and used to study close binary interactions by means of compressible ellipsoidal polytropes. Motivated by these examples, in the second scenario we study the interaction of winds by constructing binary systems in which one wind source collides against another, so that they are initially in front of the other without translation motion between them.
By applying a variant of the sink particle technique described in [16], which was implemented to the public GADGET2 code for dynamically identifying accretion centers, we observe that the first overdensities are formed not in the collision interface of the emitted winds but at the rear of the colliding sources.
These types of colliding events are not uncommon, as it has been observed that a huge number of massive stars and protostars are forming binary systems with orbital periods short enough for them to be able to exchange mass or even merge during their lifetime; see [17].
Recently, there have been several very interesting papers dedicated to the interaction of winds in binary systems; see, for instance, [18–20], where the binary components are orbiting around each other; in [21] a binary system was simulated in which each wind source had a constant mass loss rate; in [22, 23], binary systems formed by two very dissimilar components were considered; in [24] two velocity wind models were studied.
In the third scenario, we consider a rotating source that emits wind into a rotating and collapsing core. In our previous papers on collapse reported in [25–27], a regular Cartesian grid was used to make the partition of the simulation domain in small volume elements, each of which has an SPH particle located not necessarily at its center. Considering that the scale of the source is one million times smaller than the scale of the core, this Cartesian grid is inadequate when wind particles are going to be emitted from the core center, as very few grid elements will be dedicated to hosting wind particles. We also mention the technical problems caused by the very different spatial scales involved; when a source is hosted into the core, the time-step becomes prohibitively small with regard to the achievement of the evolution of the system.
For this reason, in this paper we have changed the geometry of the initial grid from Cartesian to spherical, so that a set of concentric shells are created and populated with SPH particles. This grid change allows us to mimic the mass and dynamics of the wind source, as the first concentric shell is dedicated to hosting all the wind particles.
In order to demonstrate the correctness of our spherical shell mesh to make collapse calculation, here we reconsider the collapse of a variant of the isothermal test core. Later, we go on to study the effects of the wind on this collapsing core by making use of the spherical shell mesh.
In [28], a set of SPH simulations were conducted in order to study a star formation scenario known as the collect-and-collapse model, in which an expanding region sweeps out the gas of a surrounding gas structure of uniform density. In our case, the winds are considered a small perturbation propagating radially outward into a uniform density and rotating core.
In addition, [29] reported a set of particle-based 3D hydrodynamical simulations in order to study the effects of a wind source located far from a core, modeled as a collapsing Bonnor-Ebert sphere. They implemented the HEALPix algorithm in order to achieve an isotropic wind from their source. We generate this isotropic wind in a natural way from our initial spherical mesh.
In [30], a filamentary cloud is shown to fragment into small overdensities because of the presence of turbulent perturbations. The results of this study showed that wind propagation mainly affects the early evolution of the core, such that the presence of winds does not change the collapsing nature of the rotating core but does slightly change the nature of the filament detected in the densest central region of the core and the resulting products as well.
2. Initial Systems and Models
Let us briefly describe here the physics of the systems considered in this paper and also the way in which we set up these systems in computational terms.
2.1. Setup of the First and Second Scenarios
In the first scenario we consider a hydrodynamical model of a wind source, which will be involved in a binary collision in the second scenario.
2.1.1. Initial Configuration of Source Particles
For the first and second scenarios, we consider that the wind source has a spherical and uniform distribution with radius and mass . The average source density is then given by g . We consider that the source is rigidly rotating with a linear velocity of km/s, so that its angular velocity is given by rad/s.
The escape velocity from an analogous but nonrotating mass distribution (with mass and radius ) can be estimated to be km/sec., where is Newton’s gravitational constant.
By means of a regular Cartesian mesh we make the partition of the simulation volume in small elements each with volume . At the center of each volume we place a particle (the th, say), with a mass determined by the density by means of with . Next, we displace each particle from its location by a small distance in a random spatial direction. For the first and second scenarios considered in this paper, we set thousand SPH particles to represent the source.
2.1.2. Winds
To represent the winds for the first and second scenarios, we put 650 thousand SPH particles in a thin layer around the source surface, built as explained in Section 2.1.1. The initial configuration of particles is shown in Figure 1: (a) is a colored iso-density plot and (b) is a velocity distribution plot. The former illustrates the very thin layer of high density, where the wind particles are located; the latter shows the magnitude of the velocity vector by locating an arrow at each particle site.
(a)
(b)
The wind particles initially have two velocity components; one is the rotation velocity of the source, , and the other is the radial velocity outwards from the source center, . In Figure 1(b), the wind particles seem to occupy a larger area than the source particles, as the arrow length is proportional to the magnitude of the velocity vector. The presence of the rotation component of the velocity must be noted in this panel.
As an extension to the thin layer model introduced here, for the second scenario, we will consider the case in which all the source particles become wind particles; that is, all the particles, irrespective of their initial location, will have a rotation velocity and simultaneously a radial outward velocity.
2.1.3. Thermodynamics
As a first approximation, we consider the gas of the source to be described by an ideal equation of state,where is the sound speed, which we take to be in the range ofso that the corresponding temperature associated with the source is around K.
2.2. Setup of the Third Scenario
In the third scenario, a rotating wind source is located within a collapsing and rotating core. When the source is present, the particle time-step becomes prohibitively small with regard to the achievement of the evolution of this system. To overcome this technical problem, we test here a radial mesh to be populated with particles, such that those particles located in the innermost radial shell are considered wind particles while all the rest of the particles represent the core.
2.2.1. Initial Radial Mesh for Core Particles
For the third scenario, we consider a spherical core with radius cm and mass . This core is rigidly rotating around the -axis with an angular velocity given by rad/s, so that the initial velocity of the th particle is given by .
The free-fall time is defined as the time needed for an external particle to reach the center of the core when gravity is the only force acting on this test particle. In this idealized gravitational collapse, we havewhere g cm−3 is the average density. For the uniform core in consideration here, s.
We divide the total volume of the sphere of radius into a given number of bins, , such that is the volume of a spherical shell. Each shell can be characterized by a radial interval , such that its initial and final radii are and , respectively. The radius is determined by the condition that is kept constant. Thus, we haveIn this way, the first shell is determined by the radial interval , while the second shell will be delimited by , and so on. Let us define the radius of a given shell as the average radius .
Next, by means of a Monte Carlo scheme, we populate each concentric shell with a given number of particles, , so that the particles are located randomly on all the available surface of each spherical shell. The spherical coordinates of the particles of a given shell are related to uniform random variables and (taking real values within the interval ) by the following equations: Thus, we have a total of particles distributed in the spherical volume of the core, such that the total mass in each shell is constant and therefore the global density of the core is also constant. To achieve a constant density distribution in a local sense, we further apply a radial perturbation to all the particles of a given shell such that any particle can be displaced radially outward or inward randomly, but a perturbed particle cannot cross to a different shell.
2.2.2. Initial Mass Perturbation of Core Particles
We implement a spectrum of mass perturbation on the initial distribution of core particles, such that at the end of the simulation it might result in the formation of a binary system.
If the unperturbed mass of the simulation particle is , then the mass of the particle is given by where the perturbation amplitude is set to and the mode is fixed to ; is the azimuthal spherical coordinate.
It is known that this kind of mass perturbation is very important in the formation of a filament as a result of the collapse of the isolated core. This kind of mass perturbation was successfully applied by us in our previous papers regarding collapse (see [25, 26]) and noticeably also by other authors; see, for example, [31].
It is also known that there are other different methods for the implementation of density perturbations, for example, those in which the particle mass remains unchanged; see for instance [32].
2.2.3. Winds
In order to implement the wind particles for the third scenario, we take advantage of the radial mesh described in Section 2.2.1, so that only the particles located in the first shell become wind particles. In order to increase the number of wind particles, more particles are created randomly in the neighborhood of each wind particle. We are able to give to all the wind particles a mass and angular velocity different from those of the core particles.
In all the wind models considered in this paper, the mass of the wind particle is times smaller than other involved particles, depending on the scenario: either source particle or core particle; see Section 3. Furthermore, for the third scenario, the first radial shell is populated with 10000 wind particles while the rest of the concentric shells contain 2 million core particles.
2.2.4. The Equation of State
In the third scenario, to take into account the heating of the core gas caused by core contraction and energy dissipation from artificial viscosity, we use the barotropic equation of state proposed in [33] to represent the thermodynamics:where and the critical density g cm−3; , the sound speed, is given by cm s−1, so that the corresponding temperature associated with the core is K.
2.2.5. Resolution
In all the simulations of this third scenario, we use a total of two million SPH particles which, according to the convergence study described in [25], is high enough to fulfill the resolution requirements described below.
In order to ensure the correctness of a collapse simulation, preventing the occurrence of the artificial fragmentation discovered by [34], an SPH run must fulfill some resolution criteria, expressed in terms of the Jeans mass, given bywhere is the Jeans wavelength, is the instantaneous sound speed, is Newton’s universal gravitational constant, and is the local peak density. It was shown in [35] that there is also a mass limit resolution criterion which needs to be met besides that of [34], for an SPH run to produce correct results involving self-gravity. It was then shown that both criteria are met as long as the minimum resolvable mass is always less than the Jeans mass .
For the core under consideration, and using a typical peak density of order g cm−3 and taking the number of neighbor particles included in the SPH kernel as , then we take to be given by . In this paper, the mass of an SPH core particle is , so that and therefore the resolution requirements are satisfied.
Moreover, in order to empirically demonstrate the convergence of the simulations of the third scenario, we present in the Appendix a new run made with much more particles.
2.3. The Evolution Code
We carry out the time evolution of the initial distribution of particles with the fully parallel Gadget2 code, which is described in detail by [31]. Gadget2 is based on the tree-PM method for computing the gravitational forces and on the standard SPH method for solving Euler hydrodynamic equations. Gadget2 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 , the values of which are adjusted such that for every time-step is of order unity. Gadget2 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 Gadget2 code has an implementation of a Monaghan-Balsara form for the artificial viscosity; see [36, 37]. The strength of the viscosity is regulated by setting the parameter and ; see equation (14) in [31]. Here we set the Courant factor to 0.1.
The public Gadget2 code used in this paper presents some potential technical problems when wind particles are simultaneously evolved with core particles; these problems, mainly that the particle time-step becomes prohibitively small with regard to the achievement of the overall evolution of the system, are caused by the different mass scales involved, as the discrete version of the Navier-Stokes hydrodynamical equations are written in the so-called density-entropy formulation; see [32].
Let us now briefly describe the modifications implemented into the Gadget2 code for detecting accretion centers, which is going to be used only in the second scenario. Any particle with a density higher than is an accretion center candidate. We localize all candidate particles for a given time . We then test the separation between candidates; if there is one candidate with no other candidate 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 . Thus, determines a set of particles that are within the sphere of radius and the center of which is the accretion center itself. All these particles will give up their mass and momentum to the accretion center. We then change the Gadget2 type for all these particles, and therefore they will no longer be advanced in time.
It is important to clarify that here our parameter plays the role of the inner accretion radius defined by [38] and this is why we decided to skip all the tests done on the accreted particles. The application of the limited technique described here makes sense for us, as our main aim is only to detect the places where the densest regions are formed. We have tested these modifications applied to the public Gadget2 code in [16].
2.4. Models
Because the escape velocity from a nonrotating source similar in mass and size to the one described in Section 2.1.1 is around 1000 km/s, we consider here 5 models of wind particles with a radial velocity ranging from to 2000 km/sec. We label the models A or T; either the source either is or is not left behind when the wind particles are emitted radially outward, respectively. Furthermore, we use a number to indicate the radial velocity of the model considered; see Table 1.
In order to investigate the dynamics of the wind particles when the sources are forming a binary system, we run two models, labeled A3 + A3 and T3 + T3, which according to the source enters directly in the calculation and does not do so, respectively. The binary systems are formed by locating two wind sources separated by and with no precollision velocity given, so that they are left in free fall.
In the third scenario, as we do not observe any significant change in the evolution of the core for the wind models with increasing radial velocity , then we go on to consider only wind models in which the wind particles have higher angular velocities, as shown in Table 2.
3. Results
To present the results of our simulations we consider a slice of approximately 10000 particles around the equatorial plane of the sphere: either of the source or the core. With these particles we make plots of two complementary types: one to show colored regions of iso-density and the other one to show the velocity field of the particles.
3.1. Wind Models in the First Scenario
Let us now consider the simulation outcomes obtained for the first scenario. We first notice that, for all the wind models labeled A1–A5 (of the first kind where the source is kept behind), the wind particles are separated very quickly from the source particles. However, for two low models, we see that some wind particles return to the central region, as the gravitational attraction from the source is strong enough to pull them back; see Figures 2(a) and 2(b).
(a)
(b)
(c)
(d)
The peak density increases in all these wind models because the source tends to collapse further; see Figure 3(a). For this, we make a second density plot which includes only the wind particles, so that the peak density drops as the wind particles move radially outwards, as can be seen in Figure 3(b).
(a)
(b)
We now consider the models labeled T1–T5 (in which no source is left behind). We notice that the inertia carried by the wind particles due to their rotation velocity plays a crucial role in determining the subsequent dynamics of the wind particles. This is because the particles initially located in the inner region have (linear) rotation velocity with a smaller magnitude than those particles located farther out. As the radial component of the velocity must be added with the rotation velocity to give the total wind particle velocity, then the radial motion of the exterior particles is more deflected by the rotation velocity than for inner particles.
A very important dynamical consequence of this different deflection in the radial direction of particles is that the inner particles catch up those particles initially located farther out, and, for the occurrence of this rear-end collision, a thin and dense layer of gas is formed in the shocked region of wind particles, as can be seen in Figure 4.
(a)
(b)
(c)
(d)
In order to compare the radial mass distribution in T models, in which the instantaneous radius at time is so different because of the wide range of radial velocities considered, we apply a variant of the procedure described in Section 2.2.1 for constructing the radial mesh. Thus, we first determine the maximum radius achieved at the last snapshot available for each T model; we then divide the total volume of that sphere of radius in a given number of bins, , such that is the volume of a spherical shell characterized by a variable radial interval or by a number of bins. Next, we count the number of particles located in the shell corresponding to the first bin, given by , in the second bin , and so on, such that each with is determined by the condition that is kept constant for each model. We emphasize that the radial interval of each bin is different for each T model.
In Figure 5 we show the radial mass distribution for all the T models in terms of these bin numbers. We see that the mass is dispersed around a well-defined bin number. We note that for models T1 and T2 the mass is located in the first bin while, for models T3–T5, the mass is located in subsequent bins, indicating that the higher the radial velocity of the wind models the farther out the mass shell is formed.
The time evolution of the peak density for all the T models can be seen in Figure 6. We see that for models with greater than the escape velocity, the density drops, as expected for a simple gas expansion. For models T1 and T2, we see that the density increases due to a particle condensation in the formed mass shell.
3.2. Binary Models in the Second Scenario
In this section we present and compare the results observed when a binary system is formed with two wind sources of the same type, either A3 or T3, respectively, which are the models considered for the second scenario; see Section 2.4. In these models, labeled A3 + A3 or T3 + T3, the colliding sources are of the same mass and size; thus we observe that a small bridge of particles is formed directly connecting the two A3 sources, as can be seen in Figures 7(a) and 7(b). Shortly afterwards, we see that the gravitational pull between the sources is manifested by their shape deformation while the mass of the sources feeds the longitudinal overdensity formed at the collision front: while a fraction of mass squirts out at the ends of this collision front, another fraction of mass bounces back. Therefore, there are two more slight lines of dense particles formed at the rear side of each source, as can be seen in Figures 7(c) and 7(d).
(a)
(b)
(c)
The central longitudinal overdensities seem to be fragmenting in their center while the other two overdensity lines almost disappear, as can be seen in Figures 7(e) and 7(f). It should be noted that the rotation velocity of the sources plays no role whatsoever during the interaction of winds from a binary system.
By applying a variant of the sink particle technique explained in Section 2.3, we detect the places where the densest particles are located; see [16]. Thus, we observe that the first more noticeable particle overdensities are seen to be formed in the rear of the extended front region, as can be seen in Figure 8.
When no sources are taken into account during the emission of wind particles, a contact region is quickly formed where the wind particles coming from each T3 source collide. Many of these shocked wind particles are bounced out by the opposite ends of a perpendicular line along the colliding front. The formed longitudinal overdensity is quickly fed by the sources as can be seen in Figures 9(a), 9(b), 9(c), and 9(d). Just at the end of the longitudinal overdensity, the bounced particles catch up the particles which are emitted first without any collision. We notice the formation of two well-defined overdensities; see Figures 9(e) and 9(f).
(a)
(b)
(c)
3.3. Collapse of the Isolated Core in the Third Scenario
The temporal evolution of an isolated core such as the one described in Section 2.2.1 has been reproduced by many research groups around the world using different codes based on either grids or particles; see for instance [39] and the references therein. Thus, it has been proven that an isolated rotating core contracts itself to an almost flat configuration approximately within a free-fall time of dynamical evolution, so that the central densest region of the core begins to lengthen while forming a well-defined filament that is surrounded by a halo.
In our simulation started from a new radial mesh formed by a set of concentric shells, we also see the formation of a well-defined pair of overdensities connected by a filament, as illustrated in Figures 10(a) and 10(b). It must be noted that due to the mass perturbation described in Section 2.2.2, the core has reached this resultant configuration. Next, we notice that the filament in the central region of the collapsed core fragments such that the two overdensities become separated and start orbiting around each other, as can be seen in Figures 10(c) and 10(d).
(a)
(b)
This familiar binary system resulting from the collapse of a rotating and uniform density core demonstrates the correctness of our approach by means of a radial mesh.
3.4. Collapse of the Core in the Presence of Wind
Here we follow the collapse of the core when all the particles of the first shell become wind particles. In this case, three free parameters have to be considered: the wind mass , the wind angular velocity , and the wind radial outward velocity . Only for simplicity, we fix the mass of the wind particle to be times smaller than the mass of a core particle; that is, , for every particle; see Section 2.2.2.
For a radial velocity of km/s, we first consider the simplest model, labeled Rp in Table 2, such that . That is, the wind particles rotate as a rigid body just as the core particles do. The results of this model are shown in Figure 11.
(a)
(b)
(c)
(d)
The initial concentration of wind particles can be easily seen in Figure 11(a). By the time , the wind particles have created a small hole in the central region of the core. A spherical shell of shocked particles is still visible by the time , as can be seen in Figure 11(c). For later times, the evolution of this wind model proceeds in the same manner as does the isolated core model described in Section 3.3 and illustrated in Figure 10.
As we do not observe any significant change in the evolution of this system by increasing the radial velocity , we go on to only consider wind models with higher angular velocities.
The second wind model, labeled R2, now has an angular velocity given by . The results can be seen in Figure 12. We see that the central hole is bigger and is developed earlier than the previous Rp model. At later times, we again observe the same evolution as that seen for the isolated core, but this time the filament appears to be thinner and it shows signs of fragmentation.
(a)
(b)
(c)
(d)
In the third wind model, labeled R2Om2, in which we have , we observe the formation of the central hole very early in the core’s evolution, as can be seen in Figure 13. The shell formed by the shocked particles is already expanded until a radius of at time . Another difference with previous wind models is the appearance of channels formed by the wind particles in their way through the core. The initially planted overdensities (by means of the mass perturbation described in Section 2.2.2) still appear but are now more elongated. We do not observe the formation of a bridge of dense particles connecting these initial overdensities. Nevertheless, at later evolution times, we still observe the formation of a weak filament connecting the binary system as in previous wind models.
(a)
(b)
(c)
(d)
In the third scenario considered here we observe that all the wind models do not change the collapsing nature of the rotating core, as can be seen in Figure 14; see Sections 3.3 and 3.4. Despite the fact that wind model with the highest angular velocity destroys the central core region, the core particles quickly reform this region, indicating that the overdensity seeds planted as explained in Section 2.2.2 are also unaffected, so that the binary system is formed as expected in all the models.
Despite the fact that the nature of the filament detected in the densest central region of the core is almost unchanged, the breakage of these filaments is slightly affected by the winds, for it occurs at different times and with different outcomes, as can be seen in Figure 15, where we present a comparison of the last snapshot available for each model.
(a)
(b)
(c)
(d)
It should be noted that if we further increase the initial value of the angular velocity given to the wind particles, then they would diffuse so quickly into the core that they would go across the core in a time much shorter than the time needed for the planted overdensities to be developed, and, in this case, we would not observe any change in the final results of the simulation. In order to illustrate this process of diffusion of the wind particles through the core, in Figure 16, we present another version of the colored density plots shown in Figures 13(b), 13(c), and 13(d), in which one can distinguish between the wind particles and the core particles.
(a)
(b)
4. Concluding Remarks
Very complex physical phenomena are caused by the interaction of stellar winds with their surrounding ISM, which has been observed recently, for instance, in [40], who reported distributions in the Orion-Eridanus superbubble, suggesting that the emission of a large amount of energies from massive stars can heat their surrounding material, such that large cavities can be formed; in [41] the interaction of an expanding bubble into the N49 region was discussed, likely caused by wind emission from a first generation of massive stars, so that this interaction triggers the formation of a second generation of massive stars.
In view of this complexity, the purely hydrodynamical model of a wind source considered in this paper is quite simple. However, as we mentioned earlier, most of the simulations so far cannot resolve individual stars, so that their feedback effects are collectively modeled, not individually as we try to do in this paper, where the wind source considered is dynamically characterized by its mass, its angular velocity, and the radial velocity of its emitted winds. The values shown in Section 2.1.1 for the mass and size of the source are only for purpose of illustration. In this paper we only investigate the effects and importance of considering or not considering these dynamical characteristics explicitly in simulations for three independent scenarios.
In the first scenario we find that the rotation of the source induces the formation of a mass shell formed by those emitted wind SPH particles not gravitationally escaped from the mass of the source; see Section 3.1. We also observe that the source of A models collapses further, a behavior that would be good to avoid, as a stable source would be more useful for this kind of wind model. We only mention here that the desired stability of the source may be achieved by adding a repulsive force between the source particles. For the wind T models, we observe that the mass shell of low models also collapses further, as their local density increases.
For the high models of both kinds A and T, we only observe a simple expansion of the wind particles. In this first scenario the radial velocity of the emitted winds and the angular velocity, either of the source for A models or of the wind particles for T models, are more relevant than the mass of the source in determining the final results.
In the second scenario considered here, we observe that both kinds of binary models, A3 + A3 and T3 + T3, present a dense filament formed at the collision front, formed along the perpendicular direction of the precollision source motion; see Section 3.2.
We clearly see that the filament formed in model A3 + A3 is fed by the source particles. Many of these particles are bounced back, in such a way that the first noticeable overdensities are formed not in the collision front but at the rear side of each source. We think that this occurs because the geometry of each source is deformed, first by a local conservation of angular momentum and second by the collision itself; thus, the spherical geometry is transformed to a thick disk and then to a deformed disk. For later times, as can be seen in Figure 7(c), we observe that the particles in the collision front are expelled from the ends, where we also observe the formation of some overdensities in the collision front.
This effect of forming the first overdensities some distance from the collision front is not visible in model T3 + T3. Here we observe that the filament formed along the collision front expels all its particles such that a pair of huge overdensities are visible at the ends of the filament. We are not able to evolve the simulation A3 + A3 further to confirm this observation because the computational cost of evolving the A3 + A3 model is much higher than T3 + T3 model. In this second scenario the angular velocity of the source is not at all relevant in determining the final results.
In the third scenario, we focus on a wind source centrally hosted in a rotating core. For this particular configuration, there is more evidence of the feedback effects from numerical simulations, for instance, in [42], in which the disruptive effects of gas outflows on the host cloud were demonstrated.
Here we find that the angular velocity of the wind source appears to be the most important dynamical characteristic in determining the main results of the simulation. We find that in general the emitted winds do little damage to the gas core.
If we take into account that the wind source considered in the third scenario is a massive star already in the main sequence, then it simultaneously emits photons (as ionizing radiation) and material (as stellar winds). Undoubtedly, a more complete simulation must include both kinds of emissions. When radiation emission is also considered (see, e.g., [43]), then we expect two events to occur: (i) the radiation emitted by the source creates a radiation pressure acting upon the surrounding core gas, which will then be quickly heated and thus the fragmentation of the core will be significantly diminished or even suppressed; (ii) most of the radiation particles will leak out of the core very quickly, as they travel much faster than the stellar wind particles and because the core is formed by neutral particles.
It should be noted that [44] presented the first 3D collapse calculation using an SPH code with radiative transfer included. It was demonstrated that the barotropic equation of state, like the one that we use in this paper, is suitable for describing the thermodynamics of the core collapse until the formation of the first pressure-supported condensations; beyond this limit, there are differences found in density and temperatures, between the simulation using the barotropic equation of state and that which includes radiative transfer.
In this paper, as a first approximation, the wind emission is a unique event, as it is emitted for only one time. It would be interesting to produce a situation in which the winds are created and emitted in each time-step, so that a continuous emission rate could be simulated. In this case, one would expect the collapse delay to be longer and the disruption wind effects to be more significant.
Appendix
A Higher-Resolution Run of the R2Om2 Model
An investigation on the needed resolution in SPH with regard to the achievement of convergence in the collapse of cores is presented in [25]. Because the models of [25] were constructed on a Cartesian grid while in this paper a radial grid has been used to locate the particles in the third scenario, then a demonstration of convergence is needed in order to show the reliability of the calculations presented here.
Thus, we recalculate the R2Om2 model using two times more particles to represent the wind and four times more particles to represent the core, such that a total of 8 million SPH particles are used in this test run. The convergence is achieved, as can be seen in Figure 17.
(a)
(b)
(c)
(d)
Conflict of Interests
The author declares that there is no conflict of interests regarding the publication of this paper.
Acknowledgment
Guillermo Arreaga-García would like to thank ACARUS-UNISON for the use of their computing facilities in the preparation of this paper.