# Smoothed Particle Hydrodynamics in Thermal Phases of a One Dimensional Molecular Cloud

###### Abstract

We present an investigation on effect of the ion-neutral (or ambipolar) diffusion heating rate on thermal phases of a molecular cloud. We use the modeling of ambipolar diffusion with two-fluid smoothed particle hydrodynamics, as discussed by Nejad-Asghar & Molteni. We take into account the ambipolar drift heating rate on the net cooling function of the molecular clouds, and we investigate the thermal phases in a self-gravitating magnetized one dimensional slab. The results show that the isobaric thermal instability criterion is satisfied in the outer parts of the cloud, thus, these regions are thermally unstable while the inner part is stable. This feature may be responsible for the planet formation in the outer parts of a collapsing molecular cloud and/or may also be relevant for the formation of star forming dense cores in the clumps.

## I Introduction

In most plasmas, the forces acting on the ions are different from those acting on the electrons, so naively one would expect one species to be transported faster than the other, by diffusion or convection or some other process. If such differential transport has a divergence, then it will result in a change of the charge density, which will in turn create an electric field that will alter the transport of one or both species in such a way that they become equal. In astrophysics, ambipolar diffusion refers specifically to the decoupling of neutral particles from plasma. The neutral particles in this case are mostly hydrogen molecules in a cloud, and the plasma is composed of ions (mostly protons) and electrons, which are tied to the interstellar magnetic field.

Nejad-Asghar (2007) has recently made the assumption that the molecular cloud is initially an uniform ensemble which then fragments due to thermal instability. He finds that ambipolar drift heating is inversely proportional to density and its value, in outer parts of the cloud, can be significantly larger than the average heating rates of cosmic rays and turbulent motions. His results show that the isobaric thermal instability can occur in outer regions of the cloud. The rapid growth of thermal instability results a strong density imbalance between the cloud and the low-density surroundings; therefore it may produce the cloud fragmentation and formation of the condensations.

Many authors have developed computer codes that attempt to model ambipolar diffusion. Black & Scott (1982) used a two-dimensional, deformable-grid algorithm to follow the collapse of isothermal, non-rotating magnetized cloud. The three-dimensional work of MacLow et al. (1995) treats the two-fluid model in a version of the zeus magnetohydrodynamic code. An algorithm capable of using the smoothed particle hydrodynamics (SPH) to implement the ambipolar diffusion in a fully three-dimensional, self-gravitating system was developed by Hosking & Whitworth (2004). They described the SPH implementation of two-fluid technique that was tested by modeling the evolution of a dense core, which is initially thermally supercritical but magnetically subcritical. Nejad-Asghar & Molteni (2007) have recently optimized the two-fluid SPH implementation to test the pioneer works on the behavior of the ambipolar diffusion in an isothermal self-gravitating molecular layer.

In this paper, we include the SPH equivalent of the energy equation, which follows from the first law of thermodynamics. Here, we include the ambipolar drift heating rate in the net cooling function, and we investigate the thermal phases in the self-gravitating magnetized molecular layer. For this purpose, the two-fluid SPH technique and its algorithm are given in section 2. Section 3 devotes to the basic equations of the one dimensional molecular layer and its evolution by SPH. The summary and conclusion are presented in section 4.

## Ii Numerical method

In two-fluid SPH technique of Hosking & Whitworth (2004), the initial SPH particles represent by two sets of molecular particles: magnetized ion SPH particles and non-magnetized neutral SPH particles. In this method, for each SPH particle we must create two separate neighbor lists: one for neighbors of the same species and another for those of different species. Consequently, each particle has two different smoothing lengths. In the following sections we refer to neutral particles as and , and ion particles as and ; the subscripts and refer to both ions and neutral particles. We adopt the usual smoothing spline-based kernel (Monaghan & Lattanzio 1985) and apply the symmetrized form proposed by Hernquist & Katz (1989)

(1) |

where , are positions of the particles and , respectively. is the smoothing length of particle when considering neighbors of kind . We update the adaptive smoothing length according to the octal-tree based nearest neighbor search algorithm (NNS).

The exact fraction of the total fluid that is ionized depends upon many factors (e.g. the neutral density, the cosmic ray ionization rate, how efficiently ionized metals are depleted on to dust grains). Here, we use the expression employed by Fiedler & Mouschovias (1992), which states that for ,

(2) |

where in standard ionized equilibrium state, and are valid. In reality, the gas in this case is very weakly ionized, thus, we adopt the approximation in fluid equations. In this case, the molecular cloud is considered as global neutral which consists of a mixture of atomic and molecular hydrogen (with mass fraction ), helium (with mass fraction ), and traces of and other rare molecules, thus, the mean molecular weight is given by .

The neutral density in place of neutral particles is estimated via usual summation over neighboring neutral particles

(3) |

while in place of ions, , is given by interpolation technique from the values of the nearest neighbors. The ion density is evaluated via equation (2) as follows

(4) |

for both places of ions and neutral particles. According to this new ion density, we update the mass of ion as

(5) |

in each time step so that the usual summation law

(6) |

for ions might being appropriate.

The SPH form of the drift velocity of ion particle is given by Hosking & Whitworth (2004) as

(7) |

where is the artificial viscosity between ion particles and

(8) |

where is an average density, and are the coefficients, is the mean sound speed, and is defined as its usual form

(9) |

with . Nejad-Asghar, Khesali & Soltani (2008) has recently considered the coefficients in the Monaghan’s standard artificial viscosity as time variable, and a restriction on them is proposed such that avoiding the undesired effects in the subsonic regions. Here, we use the Monaghan’s standard artificial viscosity, since the cloud contraction is quasi-hydrostatic and there is not supersonic motions and shock formation during this contraction.

It is better to keep in mind the second golden rule of SPH which is to rewrite formulae with the density inside operators (Monaghan 1992). In this case, the drift velocity is interpolated as

(10) |

where two extra density terms are introduced, one outside and one inside the summation sign. This comes as a result of the approximation to the volume integral needed to perform function interpolation (Nejad-Asghar & Molteni 2007).

There is no analytical expression that allows to calculate the value of drift velocity of the neutral particles. Here, we use the interpolation technique that starts at the nearest neighbor, then add a sequence of decreasing corrections, as information from other neighbors is incorporated (e.g., Press et al. 1992). These drift velocities at neutral places are used to estimate the drag acceleration

(11) |

instead the method of Hosking & Whitworth (2004) who used the expression of Monaghan & Kocharayan (1995). In the usual symmetric form, the self-gravitating SPH acceleration equation for neutral particle is

(12) |

where is the gravitational acceleration of particle .

The ion momentum equation assuming instantaneous velocity update so that we have

(13) |

where the first term on the right-hand side gives the neutral velocity field at the ion particle , calculated using a standard SPH approximation. Finally, the magnetic induction equation in SPH is

(14) |

where the usual notation is used.

In ambipolar diffusion process, the ion particles are physically diffused through the neutral fluid, thus, the ions will be bared in the boundary regions of the cloud (i.e. without any neutral particles in their neighbors). We check the position of ions before making a tree and nearest neighbor search, so that we send out the bared ion particles in the boundary regions of the simulation at next time-step. Both the cloud and boundary regions contain ion and neutral particles. We set up boundary particles ( up and down in ) using the linear extrapolation approach (from the values of the inner particles) to attribute the appropriate drift velocity, drag acceleration, pressure acceleration, and the magnetic induction rate to the boundary particles. The algorithm of the two-fluid SPH for simulation of thermal phases in a one dimensional molecular cloud is shown in Fig. 1.

## Iii layer evolution

### Iii-a basic equations

We write the continuity equation of neutral part as its common form

(15) |

while, the relation (2) is used to determine the ion density wherever it is required in the two-fluid equations. The momentum equation of neutral part then becomes

(16) |

where the gravitational acceleration obeys the poisson’s equation

(17) |

and the pressure is given by the ideal gas equation of state

(18) |

where is the mass of hydrogen atom and the temperature is approximated the same as for both neutral and ion fluids ().

The thermal energy per unit mass is

(19) |

where the mean internal (rotation and vibration) energy of an molecule is included. The energy equation follows from the first law of thermodynamics, that is

(20) |

where is the net cooling function

(21) |

where and are the heating rates due to cosmic rays and ambipolar diffusion, respectively, and and are the parameters for the gas cooling function that here we use the polynomial fitting functions, outlined by Nejad-Asghar (2007), as follows

(22) |

(23) |

where .

The magnetic fields are directly evolved by charged fluid component, as follows:

(24) |

where the last term outlines the ambipolar diffusion effect with drift velocity,

(25) |

which is obtained by assumption that the pressure and gravitational force on the charged fluid component are negligible compared to the Lorentz force because of the low ionization fraction. The notation in drift velocity represents the collision drag coefficient.

### Iii-B results

The chosen physical scales for length and time are , and , respectively, so that velocity unit is approximately . The Newtonian constant of gravitation is set for which the calculated mass unit is . Consequently, the derived physical scale for density, energy per unit mass, and drag coefficient are , , and , respectively. In this manner, the numerical values of and are and , respectively. The magnetic field is scaled in units such that the constant is unity. Since the magnetic flux density has dimensions

(26) |

while has dimensions

(27) |

specifying therefore scales the magnetic field equal to . With aforementioned units, the thermal energy per unit mass (19) is represented by , the heating rates due to cosmic rays and ambipolar diffusion are and

(28) |

respectively, and the parameters for the gas cooling function are

(29) |

(30) |

The initial conditions for this simulation are a parallel magnetic field directed perpendicular to the -axis so that the initial ratio of magnetic to gas pressure is everywhere a constant (), and a density profile given by

(31) |

where is the initial central density and is a length-scale parameter according to the initial sound speed (see Fig. 1 of Nejad-Asghar & Molteni 2007). The magnetic field is assumed to be frozen in the fluid of charged particles and the central density is assumed to be . We choose a molecular cloud with the mass fraction of molecular hydrogen and helium are and , respectively, and have initial uniform temperature of . We assume that the cloud layer is spread from to (according to Fig. 1 of Nejad-Asghar & Molteni 2007). The initial values of the cooling and heating functions are shown in Figure 2a. As presented in this figure, the isobaric thermal instability criterion,

(32) |

is satisfied in the outer parts of the cloud, thus, these regions are thermally unstable while the inner part is stable as outlined by Nejad-Asghar (2007).

The present SPH code has the main features of the TreeSPH class so that the nearest neighbors searching are calculated by means of this procedure. We integrate the SPH equations using the one order simple integration scheme. The selection of time-step, , is of great importance. There are several time-scales that can be defined locally in the system. For each particle , we calculate the smallest of these time-scales using its smallest smoothing length, , i.e.

(33) |

where is the Alfvén speed of ion fluid and is the Courant number which in this paper is adopted equal to (for numerical stability). The evolution were carried out to time so that the stability of the outer parts of the cloud may be revealed. The values of the cooling and heating rates versus position (and neutral density) are shown in Fig. 2b-d, at times , and , respectively.

## Iv Summary and conclusions

In the preceding work of Nejad-Asghar & Molteni (2007), two-fluid SPH implementation was used to check the pioneer works on behavior of the ambipolar diffusion in an isothermal self-gravitating layer. In this paper, we include the SPH equivalent of the energy equation that can be obtained from the first law of thermodynamics. Here, we incorporate the ambipolar drift heating in the net cooling function, and we investigate the thermal phases in the self-gravitating magnetized molecular layer.

The values of the cooling and heating functions are shown in Figure 2. According to this figure, the isobaric thermal instability criterion is satisfied in the outer parts of the cloud, thus, these regions are thermally unstable while the inner part is stable. The SPH equations were integrated, and the evolution were carried out to time so that the stability of the outer parts of the cloud is revealed.

It is obvious that the instability of the cloud at the outer parts, causes the formation two relative cool regions in that area. The rapid growth of thermal instability results a strong density imbalance between the cloud and the low-density surroundings. This feature may be responsible for the planet formation in the outer parts of a collapsing molecular cloud and/or may also be conscientious for the formation of star forming dense cores in the clumps.

## Acknowledgment

The MNA would like to thank the Research Council of the Damghan University of Basic Sciences.

## References

- [1] D.C. Black and H. Scott, Astrophys. J., vol. 263, p. 696, 1982.
- [2] R.A. Fiedler and T.C. Mouschovias, Astrophys. J., vol. 391, p. 199, 1992.
- [3] L. Hernquist and N. Katz, Astrophys. J. Sup., vol. 70, p. 419, 1989.
- [4] J.G. Hosking and A.P. Whitworth, Mon. Not. Royal Aston. Soc., vol. 347, p. 994, 2004.
- [5] M.M. MacLow, M.L. Norman, A. Konigl and M. Wardle, Astrophys. J., vol. 442, p. 726, 1995.
- [6] J.J. Monaghan, Ann. Rev. Astron. Astrophys., vol. 30, p. 543, 1992.
- [7] J.J. Monaghan and A. Kocharyan, Comp. Phys. Commun., vol. 87, p. 225, 1995.
- [8] J.J. Monaghan and J.C. Lattanzio, Astron. Astrophys., vol. 149, p. 135, 1985.
- [9] M. Nejad-Asghar, Mon. Not. Royal Aston. Soc., vol. 379, p. 222, 2007.
- [10] M. Nejad-Asghar and D. Molteni, Joint Eroup. National Astron. Meeting, Yerevan, Armenia, 20-25 August, 2007, asto-ph0802.2437.
- [11] M. Nejad-Asghar, A.R. Khesali and J. Soltani, Astophys. Space Sci., vol. 313, p. 425, 2008.