2 Numerical techniques and gas dynamics

The present work adopts the traditional isotropically-adaptive SPH, where the smoothing length is proportional to the radius of a spherical region inside which there is a fixed number of nearest particles. However, more refined approaches of adaptive interpolation, regarding anisotropy of the mass distribution, may be seen in the literature (e.g., [Fulbright et al. 1995]; [Owen & Fisher 1994]; [Owen et al. 1998]; [Shapiro et al. 1996]).

We adopted the B-spline ([Monaghan & Gingold 1983]) to model the smoothing kernel, defined in 3-dimensions as:

$$W_h(r)={3\over 4\pi h^3} \cases {{{4\over 3}-{2}\frac{r^2}{h^... ...}}\bigr)^3}, & $h\le r< 2h$ ,\cr\cr 0, & $r\ge 2h$ . \cr\cr }$$ (1)

From the equation above, the smoothing length h can be calculated as half the radius of a spherical region, enclosing a fixed number of the nearest neighbors of a given particle, placed in position ${\vec r}$.

We calculate the smoothed particle gradient of any fluid parameter A, for a particle i, by using the traditional form:

$$\langle\nabla{A}\rangle_i=\sum_j\frac{m_j}{\rho_j}{\vec\nabla}_iW_{ij}A_j,$$ (2)

where any particle density, $\rho_i$, is estimated by

$$\rho_i=\sum_j{m_j}W_{ij},$$ (3)

and the index j points to the particles in a list of the nearest neighbors of particle i, separated by distances smaller than a given radius. For the present kernel, this radius is naturally 2h_i.

The two-particle smoothing kernel, W_ij, appearing in equations above is the symmetrized gather-scatter form, proposed by Hernquist & Katz (1989):

$$W_{ij}\equiv {W(\vert{\vec r}_i-{\vec r}_j\vert,h_i)+W(\vert{\vec r}_i-{\vec r}_j \vert,h_j)\over 2}$$ (4)

with $W(\vert{\vec r}\vert,h)$ given by Eq. (1).

2.1 Momentum conservation

The usual form of writing the self-gravitating SPH momentum equation requires a symmetrized two-particle combination of pressure and artificial viscosity (see, e.g., Hernquist & Katz 1989; Monaghan 1992; Benz 1990):

$${\frac{{\rm d}{\vec v}_i}{{\rm d}t}}= -\sum_jm_j {\vec\nabl... ...\frac{p_j}{\rho _j^2}}+\Pi_{ij} \biggr) -\vec\nabla_i\Phi_i,$$ (5)

where $\vec\nabla_i\Phi_i$ is the acceleration of gravity on particle i, and $\Pi_{ij}$ is the Monaghan-Gingold artificial viscosity (Monaghan & Gingold 1983):

$$\Pi_{ij}\;=\; \cases { {\alpha\bar c_{ij}\mu_{ij}+\beta\mu_{i... ...${\vec v}_{ij}\cdot{\vec r}_{ij} < 0$ ,\cr 0,&otherwise,\cr }$$ (6)

with

$$\mu_{ij}=-{\frac{{\vec v}_{ij}\cdot{\vec r}_{ij}}{\bar h_{ij}}}{\frac 1{r_{ij}^2/\bar h_{ij}^2+\eta ^2}},$$ (7)

and the following convention is assumed: $a_{ij}\equiv a_i-a_j$ and $\bar a_{ij}\equiv 0.5(a_i+a_j)$. The coefficient $\eta$ is usually 0.1, and the coefficients $\alpha\sim1$ and $\beta\sim2$. However, we heuristically found the values $\alpha =3$, $\beta =5$, to be a good choice to stabilize the hierarchical integration scheme. For lower values, the solutions diverged on the critical phase of the shock in an adiabatic experiment of 3-D spherical collapse, shown in Appendix C.

The coefficient $\bar c_{ij}$ in Eq. (6) is usually defined as the arithmetic mean of the adiabatic sound speed, $c_{_{\rm S}}$, for particles i and j respectively. However, we found that the artificial viscosity works efficiently if we use the exact definition of $c_{\rm S}$ in adiabatic shocks but it allows particle interpenetration if we adopt the present molecular mixture and cooling, mainly in the post-shocking layers where it has cooled down to lower temperatures as $\sim{10}$ K. The main reason for this effect is that a considerable reduction in the specific thermal energy occurs at lower temperatures, which reduces considerably the local sound speed. We minimized this effect replacing $c_{\rm S}^2$ by the equivalent specific thermal energy of a (non-physical) 100% dissociated hydrogen cloud:

$$u_{\textup{H\, {\mdseries\textsc{i}}}}={\frac 32}\biggl( X+{\frac Y4}\biggr) R_{\rm o}T+{\frac{X{D}_{\rm o}}{2m_{{\rm H}}}},$$ (8)

so that $\bar c_{ij}=0.5(u_{{\textup{H\, {\mdseries\textsc{i}}}}_i}+u_{{\textup{H\, {\mdseries\textsc{i}}}}_j})$, where $D_{\rm o}$, X and Y are discussed below (Sect. 2.2).

Monaghan (1989) discussed an alternative solution to the problem of particle penetration in SPH. However, we have not adopted such device since the integration scheme has shown numerical instability on supersonic shocks in our calculations.

  
2.2 Chemical equilibrium and thermal energy

The interstellar clouds in our experiments are considered a mixture of atomic and molecular hydrogen, helium and traces of CO. Only H, H₂ and He contribute effectively to the thermal energy; the role of CO is limited to the cooling of the gas, at temperatures lower than 10⁴ K and at high densities. We assume that the chemical abundance (in mass) is X = 0.75 and Y = 0.25. The chemical equilibrium between H and H₂ is the same of ML91, which is only valid in dark regions. The relative abundance, y, of atomic H obeys the equilibrium equation:

$${\frac y{1-y}}={\frac{2.11}{\rho Xy}}{\rm e}^{-{D}_{\rm o}/kT},$$ (9)

where

$$y={\frac{\rho ({\rm H})}{\rho X}},$$ (10)

and $D_{\rm o}$ = 4.477 eV is the dissociation energy of H₂ molecule. The specific thermal energy for the gas mixture is given by

$$u(T)={\frac 32}{\frac{{R}_{_{\rm o}}}\mu }T +{\frac{Xy{D}_{_{\rm o}}+X(1-y)E({\rm H}_2)}{2m_{_{{\rm H}}}}},$$ (11)

and the equation of state is given by

$${\frac p\rho}={\frac{{R}_{\rm o}}\mu }T.$$ (12)

with $\mu$ being the mean molecular weight:

$${\frac 1\mu }=\left( {\frac{1+y}2}\right) X+{\frac Y4}.$$ (13)

The function $E({\rm H}_2)$ is the mean internal (non-translational) energy of the H₂, given by

$$E({\rm H}_2)=f_{_{{\rm rot}}}+{\frac{k\theta_{_{{\rm vib}}}}{\exp(\theta _{_{{\rm vib}}}/T)-1}},$$ (14)

where $f_{_{{\rm rot}}}$ is the rotational internal energy, that we extracted from Eq. (3.9) of ML91, and $\theta _{_{{\rm vib}}}$ is 6100 K.

Temperatures are calculated by iterations from Eq. (11) predicting more accurate values by solving T from the linear term. In general, convergence is reached within a good tolerance ( $\sim 10^{-4}$) after two or three iterations.

  
2.3 Energy conservation. Cooling and heating

Figure 1: The total cooling rate in ergs s-1 cm-3 for the mixture CO, H2, HI and HII given in Eq. (17) (see text) for six different number densities of hydrogen species: 10-2, 100, 102, 104, 106, 108 cm-3

The particle form of the first law of thermodynamics is expressed as

$${\frac{{\rm d}u_{_i}}{{\rm d}t}}={\frac 12}\sum_jm_j{\vec v}_{{ij}}\cdot {\vec \nabla }_iW_{{ij}}\Pi _{ij}-{\cal L}_i,$$ (15)

where we define ${\cal L}_{_i}$ as the rate of specific energy loss:

$${\cal L}_{_i}={\frac{\Lambda _{_i}-\Gamma_{_i}}{\rho _{_i}}}.$$ (16)

$\Lambda _i$ and $\Gamma _i$ being the cooling and heating functions, respectively.

The cooling model takes into account two main processes that dominate in different ranges of temperature. Below 10⁴ K the cooling is mostly due to molecules; we adopt the model of Hollenbach & McKee (1979, hereafter HM79). Above 10⁴ K, free-free emission by electrons is the main cooling mechanism, which is well described by the DM72 cooling function. In the simulations presented in this work, the temperature reached values just above 10⁴ K.

The total expression for the specific cooling rate in erg cm^-3 s^-1 is:

$$\Lambda =8.5\ 10^{-5}n_{_{{\rm H}}}^2L({\rm CO})+n({\rm H}_2)^2L({\rm H}_2)$$

$$+n({\rm H})n({\rm H}_2)L({\rm H})+n({\rm H})^2L(\textup{H\, {\mdseries\textsc{ii}}}),$$ (17)

where $n_{\rm H}$ is the number density of hydrogen species, $n({\rm H})$ is the number density of atomic hydrogen, and $n({\rm H}_2)$ is the number density of molecular hydrogen:

$$n_{_{{\rm H}}}=n({\rm H})+n({\rm H}_2).$$ (18)

The functions $L({\rm CO})$, $L({\rm H})$, $L({\rm H}_2)$were taken from the work of HM79, which is well summarized by ML91. The cooling function L(H II) is the ionized hydrogen contribution that is given by the interpolation formula of DM72:

$$L(\textup{H\, {\mdseries\textsc{ii}}})=10^{-21}\biggl[10^{-0.1-1.88(5.23-\log_{10}T)^4}$$

$$+10^{-1.7-\Theta (T)}\biggr]({\rm erg~cm^3~s^{-1}}),$$ (19)

with

$$\Theta(T) = \cases{0.2\,(6.2-\log_{10}T)^4,&if $\log_{10}T\,\leq\,6.2$ ,\cr0,&otherwise.\cr}$$ (20)

This contribution is effective only in the initial stages of the collision experiments, where the gas temperature reaches values comparable to that one of an H II region.

A plot of the cooling rate $\Lambda$ as a function of the temperature, for different number densities of hydrogen species, is shown in Fig. 1. The peak around 10³ K is due to the cooling by H₂, expressed by the functions $L({\rm H})$ and $L({\rm H}_2)$ reported by HM79. These functions rise strongly between 10² and 10³ K, as shown in Figs. 7 and 8 of HM79. The steep decrease of the H₂ contribution at about 10³ K is due to the dissociation of H₂. According to Eq. (9), the relative abundance of molecular H, is affected by a strong transition in chemical equilibrium near 10³ K. This decrease does not appear in the figures of HM79 because they are interested in the cooling per H species and not in the total effective cooling per cm³. The CO contribution is evident in Fig. 1 for high number densities and T<10⁴ K, since its efficiency increases with the square of the density of H species. One can see in Fig. 1 that the cooling law is represented by a straight line, for T<10⁴ K and for $n_{\rm H}$ $\sim{10^8}$, except for a small peak due to H₂ cooling. The H₂-cooling profile is evident at low densities, at temperatures lower than 10⁴ K.

The heating mechanism is assumed to be a combination of the effect of cosmic rays and of H₂ production by grains, according to the expressions given by Spitzer (1978). Heating by H₂ formation on dust grains is effective below a critical temperature limit, of the order of 10² K, that depends on the surface texture of the grain. There are possibly many other heating processes in the interstellar medium. For instance, one possible mechanism we have not included here is the heating due to gas-grain collisions (e.g., [Spitzer 1978]; [Tielens & Hollenbach 1985]). We believe the mechanisms that we have taken into account are the dominant ones.

The heating by cosmic rays is given by

$$\Gamma_{_{{\rm CR}}} \approx 3.8\ 10^{-29}n_{_{{\rm H}}}{\rm (erg~cm^{-3}~s^{-1})},$$ (21)

and for the heating by H₂ molecules evaporated from dust grains,

$$\Gamma_{_{{\rm H}}d} \approx 2.2\ 10^{-28}n_{_{{\rm H}}}n({\rm H}){\rm (erg~cm^{-3}~s^{-1})}.$$ (22)

The total heating per unit volume is given by

$$\Gamma =\Gamma _{_{{\rm CR}}}+\Gamma _{_{{\rm H}}d}.$$ (23)

2.4 The computer code

2.4.1 Basic implementation

The present code belongs to a new class of hydrodynamic codes, called TreeSPH, which was originally developed by HK89. The main characteristic of TreeSPH is the computational reuse of the tree-descent on both gravity-calculations and neighbors-searching. Essentially, we followed the directions given by classical works in the literature (e.g., [Barnes 1990]; [Hernquist 1990]; HK89; [Makino 1990]) to construct the tree-subroutines.

One particular application of the SPH-interpolation technique is the program we developed to visualize the simulation results in gray scale or in pseudo-color images (cf. Figs. 6 and 7). To be visualized, each particle quantity (e.g., density, temperature etc.) is smeared out in a 2-D pixel-grid, and conveniently converted to pixel quantities (8 or 24-bits). As the out-put file of a given simulation stores relevant SPH information (e.g., the smoothing length h), the particle smoothing-length itself is used to cover some pixels inside a radius $r_i=2h_i/\sigma$, where $\sigma$ (in length unit per pixel) is the grid scale-factor. It may also be used to control the image resolution. The interpolation is performed with particles that obey some criterion (e.g., with densities or temperatures in a given range, z-coordinate tolerance to define a slide etc.).

A more complete version of the code, including a more complex cooling for the mixture H, H₂, OH, H₂O, CO, with the H, H₂ and H₂O chemistry shall be presented in a future paper.

2.4.2 Gravity forces and neighbors searching

Both gravity acceleration and potential field are calculated by means of the octal-tree scheme, called BH method (e.g., [Barnes & Hut 1986], 1989; [Hernquist 1987], 1988, 1990; [Makino 1990]), with the quadrupole approximation properly derived for the softened potential (see Appendix A).

Essentially, the same procedure to perform gravity calculations is applied to perform the nearest neighbors searching. The difference in both implementations, besides the specific calculations, is the stop-condition for tree-descents. On gravity computations, the stop condition is the essential of the BH method, which is described in several papers concerning this matter (e.g., Barnes & Hut 1986, 1989; Hernquist 1987, 1988).

In the nearest-neighbors searching (NNS), a tree-descent is performed in order to assemble for each particle a list of the nearest other particles that are separated with distances below or equal to a given radius R. The NNS tree-descent through a given cell stops if this one is unable to contain any neighbor inside, which may be translated to the following condition:

$$\max\{\vert x_i-x_{\rm c}\vert,\vert y_i-y_{\rm c}\vert,\vert z_i-z_{\rm c}\vert\}>R+s,$$ (24)

where (x_i,y_i,z_i) is the Cartesian coordinates of the fixed particle, with generic index i, $(x_{\rm c},y_{\rm c},z_{\rm c})$ is the Cartesian coordinates of cell's mass-center, and s is the cell's size.

If the above condition is not satisfied, the focused cell may contain neighbors, and the tree-descent is resumed until the newly pointed cell encloses a single particle inside. This new candidate, with index j, is one of the nearest neighbors if the following condition is true:

$$\vert{\vec r}_j-{\vec r}_{\rm c}\vert\le R,$$ (25)

and then, this particle index is stored in the neighbors-list for the fixed particle. After a complete search, each particle have their corresponding nearest-neighbors lists.

The NNS algorithm itself is applied to determine the smoothing lengths. The NNS routine returns, for particle with index i, a list of the nearest neighbors within an input radius R_i. If the number of neighbors found within R_iexceeds an expected number $N_{{\rm f}}$, the list is sorted according to the neighbors distances, and then the exceeding particles are discarded from the neighbors-list ([Heller 1993]). Thus, the corrected radius is the maximum distance in the remaining neighbors- list. The smoothing length, in the case of the adopted kernel model is then $h_i=0.5\,R_i$.

If the number of neighbors found by the NNS procedure is less than the expected $N_{{\rm f}}$, the new initial radius for neighbors search is twice the old value so that a new cycle of search-and-sort is required until a convergence criterion is satisfied. It has been adopted adopted a tolerance of $\pm 2$ neighbors in the present work.

For SPH computations, it is required a variant of the above algorithm to assemble a list of the nearest particles that give non-zero contribution to the smoothing kernel. This effective-neighboring condition may enlarge the list for some configurations, reaching increases of about $50\%~N_{\rm f}$. We call this algorithm the effective-neighbors-search (ENS), which consolidates the construction of the final neighbors list. Along with an ENS tree-descent for a given particle i, any resolved particle j is assigned as an effective neighbor if the relation $d\leq 2\max(h_i,h_j)$ is satisfied.

The adopted gravity softening is the classical Aarseth's approach, where each particle is is gravitationally equivalent to a Plummer's sphere with the same mass of the original particle and with a half-mass radius $\epsilon$. This approach does not reproduce Keplerian orbits even for a pair of particles separated far beyond the effective-radius ( $d>2\epsilon$). This problem is overcome by the softening scheme of HK89, where the B-spline kernel itself is applied to perform the softening.

At the initialization of a new SPH (or pure-gravitational) experiment, the program calculates iteratively the maximum softening length $\epsilon_{\rm max}$ as a function of the gravity linear-scale:

$$\epsilon_{\rm max}= \biggl({-{G}M_{\rm tot}\over E_{\rm G}}\biggr){N_{\rm tot}}^{-1/3},$$ (26)

where $N_{\rm tot}$is the total number of particles in the system, $M_{\rm tot}$ is the system's total mass, $E_{\rm G}$ is the total gravitational energy, calculated by means of the gravitational part of the treecode itself. $\epsilon$ is then estimated as the mean interparticle spacing for an equivalent homogeneous self-gravitating system with the above characteristics.

The energy $E_{\rm G}$ depends on the chosen value for $\theta$ and the initial value for $\epsilon$, which can be zero or, if a pertinent system's linear scale lis known, its initial value can be $\epsilon=l\,N^{-1/3}$ in order to accelerate the $\epsilon$-determination. Then, a new value for $\epsilon$ is returned from Eq. (26). This corrected value for $\epsilon$ can be used as a new initial value to repeat the above procedure, allowing then a more accurate value for $\epsilon$. In general, three iterations are sufficient to reach a good value for $\epsilon$.

In order to perform calculations with spatially varying softening length, we adopted the scatter form (HK89) since this allows quadrupole approximation without involving spatial derivatives of $\epsilon$. However, we shall impose an upper limit $\epsilon_{\rm max}$to the softening length, either being the one calculated from Eq. (26), at the initialization of the experiment, or an input value, so that the definitive softening length is written as

$$\epsilon=\min(\epsilon_{\rm max},\, s)\geq\epsilon_{\rm min},$$ (27)

with s being the cell's size. The lower limit $\epsilon_{\rm min}$is used to avoid clumps to collapse indefinitely in dissipative simulations, yielding numerically unsolvable scale-contrasts. Of course, the adaptive $\epsilon$ in Eq. (27) does not allow energy conservation since the corresponding Hamiltonian is no longer constant. However, this energy error is negligible in comparison to the high amount of energy loss due to the efficient cooling adopted in the present work. Particularly, it has been noticed that the propagated errors, adopting Eq. (27), are in general less significant than common conservation-errors of SPH simulations of adiabatic collapses. However, on simulations of pure gravitational collapse of low-velocity-dispersion spherical systems, the results were catastrophic at the critical passage of the collapse (splash). An alternative approach to perform tree-gravity calculations, with variable resolution, is presented by Nelson & Papaloizou (1994), where the tolerance parameter, $\theta$, is halved if a given cell has distance smaller than $\epsilon$.

Figure 2: The hierarchical scheme for the multiple time-step leap-frog. In the figure it is assumed the pre-fixed integration, so that upper level particles (coarser time-steps) are integrated before a descent

2.5 The integration scheme

2.5.1 Time-stepping

The multiple time-scale integration-scheme requires bookkeeping to control what particles are integrated with a given time-step, and how a time-step is selected regarding synchronization of the finite-differences solutions. One possible time scheme is the one illustrated in Fig. 2, where the finite-differences equations are solved along with a binary hierarchy of time-steps.

For each time-step, a list (time-bin) of dynamically similar particles is assigned. Each time-bin is indexed by a time-depth, $l=\{0,1,...\}$, which, in turn, labels a time-step $\Delta{t}_l=2^{-l}\Delta{t}_0$, where $\Delta{t}_0$ is the root time-step that is calculated in the program initializations. For this reason, a time-step must be an ease-to-handle floating-point, so that divisions by two will not compromise the time integrity. Thus, the chosen initial time-steps must be numbers like 0.0078125, 0.00390625 etc. Numbers such as $1/\pi$, $\exp{-23}$ etc. will round-off on divisions by two, inducing synchronization breaking, and the program will hang due to infinite loops or it will at least incur in wrong results.

An upper-limit time-step is calculated from stability/accuracy considerations, which may be associated to the local physical properties. Regarding synchronization, the definitive time-step is determined as the maximum power of two that is smaller than or equal to the initial value.

For pure gravitational purposes the maximum time-step is adopted as being half of the cell's size divided by the particle's velocity (e.g., Jernigan & Porter 1989), assuming that the octal-tree itself is the standard of rest. For SPH simulations, the maximum individual time-step is calculated by the Courant stability condition:

$$\tau^{\rm c}_{i}=\frac{h_i}{\max(c_{s_i},\sum_j\mu_{ij},...)},$$ (28)

 so that we determine the particle time-step (time-bin) as being the maximum time-step satisfying:

$$\tau_k=2^{-k}\tau_{\rm o}\leq\tau^{\rm c}_{i}.$$ (29)

The time-stepping algorithm descends recursively to the deepest time-bin. Before a descent, a time accumulator (indexed by the current time-bin) is set equal to zero. After the descent is exhausted, the current time-step is accumulated; the particles of the current time-bin are simultaneously integrated, until the returning condition of the accumulated time for the current level equals the time-step of the parent time-bin. The time-stepping order may be changed from pre-fixed integration the post-fixed one by placing the leapfrog subroutine-call at the end of the time-bin descent subroutine-call. However, we have noticed that in some tests the post-fixed integration increases the numerical cooling in compressive shock simulations.

2.5.2 Modified leapfrog

The SPH equations of motion are integrated by means of a time-adaptive leapfrog, so that corrections due to the changing of time-steps are explicitly included in the finite-differences equations. A particle may jump from a time-bin to another so that along its history the particle may have visited several time-bins. The finite-differences scheme is derived in an approach similar of that given by Hernquist & Katz (1989), so that positions and specific thermal energies are time-centered and velocities are synchronized with the time-level. Thus, if the SPH particles are dynamically described by their position vector, ${\vec x}_{n+1/2}$, velocity, ${\vec v}_n$, and specific thermal energy, u_n+1/2, we have the following finite difference equations:

$${\vec x}_{n+1/2}={\vec x}_{n-1/2}+\bar\tau_n{\vec v}_n+\bar\tau_n\delta\tau_n {\vec a}_{n-1/2}$$ (30)

$${\vec v}_{n+1}={\vec v}_n+\tau_{n+1/2}{\vec a}_{n+1/2}$$ (31)

and

$$u_{n+1/2}=u_{n-1/2}+0.5\bar\tau_n(\dot u_{n+1/2}+\dot u_{n-1/2})$$ (32)

where $\bar\tau_n=0.5(\tau_{n+1/2}+\tau_{n-1/2})$, $\delta\tau=0.25(\tau-\tau_{n-1/2})$, and the time-step for the current time-bin is $\tau_n+1/2$ and the old individual time-step was $\tau_{n-1/2}$.

Time-centered velocities, appearing in the acceleration and in the specific thermal-energy-rate, are predicted by

$${\vec v}_{n+1/2}=0.5({\vec v}_{n+1}+{\vec v}_n).$$ (33)

A more detailed derivation of the above leap-frog equations is given in Appendix B.