# Formation of Primordial Protostars

## 1 Introduction

According to the standard scenarios of cosmology, the hydrogen becomes neutral at $`z10^3`$. The absence of Ly$`\alpha `$ trough in quasar spectra, however, severely constrains the neutral fraction of the uniformly distributed intergalactic medium at $`z<5`$. This fact indicates that the hydrogen has to be reionized by energy injection from some objects in an intervening period. It is suggested that ultraviolet radiation from the known quasar population is insufficient to photoionize the universe (Shapiro & Giroux 1987). Recently, the role of the first stars in reionizing the universe has been emphasized by many authors (e.g., Fukugita & Kawasaki 1994; Ostriker & Gnedin 1996; Haiman & Loeb 1997).

However, star formation in the galaxy formation epoch is one of the most poorly understood processes involved in galaxy formation. Consequently, theories on the stellar reionization are undermined by many uncertainties, for example, the initial mass function, star formation rate and efficiency of the first stars. In order to resolve these uncertainties, detailed knowledge of the formation processes of a first star from primordial gas seems to be required in advance.

Primordial star formation can be viewed as the successive fragmentation and contraction processes of collapsed cosmological objects. In this paper we discuss the latter process, namely the contraction of primordial protostellar clouds <sup>1</sup><sup>1</sup>1We use the term “protostellar cloud” to identify a gravitationally bound contracting gas cloud that, in part or as a whole, becomes a star as it evolves. On the other hand, the term “protostar” means an accreting core in hydrostatic equilibrium formed at the center of a protostellar cloud into stars. This problem has been investigated with a one-zone approximation by e.g., Carlberg (1981); Palla, Salpeter, & Stahler (1983; see also Matsuda, Sato, & Takeda 1969; Yoneyama 1972).

On the other hand, Stahler, Palla, & Salpeter (1986) studied the main accretion phase of a primordial protostar assuming stationary accretion of ambient material onto a central hydrostatic core, which is a primordial counterpart of the present-day stellar core (Larson1969). In that paper, they assumed Shu-like collapse (Shu 1977), i.e., the cloud initially in an unstable hydrostatic equilibrium is overrun by an expansion wave that propagates outward from the center, leaving behind a growing central point mass and a free-fall velocity field. Further, the mass accretion rate is constant with time and its rate is $`=0.975c_\mathrm{s}^3/G`$, where $`c_\mathrm{s}`$ is the sound velocity in the cloud. Higher temperatures, $`10^3`$ K, in primordial clouds results in a higher mass accretion rate $`(4.41\times 10^3M_{\mathrm{}}\mathrm{yr}^1)`$ than the present-day value ($`10^5M_{\mathrm{}}\mathrm{yr}^1`$, e.g. Stahler, Shu, & Taam 1980). However, it remains unclear how the primordial protostellar cloud reaches the main accretion phase.

Villere & Bodenheimer (1987) investigated this problem by using hydrodynamical methods. They studied the spherically symmetric collapse of a primordial cloud of $`5\times 10^5M_{\mathrm{}}`$ into a protostar and pointed out the importance of H<sup>-</sup> cooling for the determination of the minimum mass of primordial stars. However, they did not treat the radiative transfer of H<sub>2</sub> lines, and their efficiency of cooling in optically thick H<sub>2</sub> lines was different from that of Palla et al. (1983). Therefore, in order to follow the collapse of primordial protostellar clouds accurately, we need to employ the hydrodynamical method together with solving the radiative transfer of H<sub>2</sub> lines.

In this paper we study the evolution of primordial protostellar clouds by performing hydrodynamical calculations, with radiative transfer of H<sub>2</sub> lines and continuous radiation for spherically symmetric clouds up to the main accretion phase.

The outline of this paper is as follows. In §2, we describe the method of our calculation. In §3, results of calculations are presented. We summarize our work and discuss its implication in §4. In the Appendix, our treatment of radiative transfer is described.

## 2 Method of Calculation

### 2.1 Hydrodynamics

We assume the cloud is spherically symmetric and neglect rotation, magnetic field, and external radiation field for simplicity. The deviation from spherical symmetry and the effects of rotation may be important. However, these effects make the problem two or three-dimensional and intractable so far. In primordial star formation, magnetic fields are considered to play a less important role than in present-day star formation. Since star formation takes place in the deep interior of a protogalactic cloud, and, moreover, in such a region the electron fraction is tiny, the effects of cosmic background radiation are not significant.

The hydrodynamic code we use is a usual spherical symmetry, explicit Lagrangian finite-difference scheme with von Neumann-Richtmyer artificial viscosity (e.g., Mihalas & Mihalas 1984; Thoul & Weinberg 1995).

The Lagrangian equations of hydrodynamics describing the collapse of a spherically symmetric cloud are the continuity equation

$$\frac{m}{r}=4\pi r^2\rho ,$$
(1)   
the momentum equation

$$\frac{Dv}{Dt}=4\pi r^2\frac{p}{m}\frac{Gm}{r^2},$$
(2)   
the energy equation

$$\frac{Dϵ}{Dt}=p\frac{D}{Dt}(\frac{1}{\rho })\frac{\mathrm{\Lambda }}{\rho },$$
(3)   
and the equation of state

$$p=(\gamma _{ad}1)\rho ϵ.$$
(4)   
In the above equations, the quantities $`m,\rho ,v,p,\epsilon `$, and $`\mathrm{\Lambda }`$ are the mass within radius $`r`$, the density, velocity, pressure, thermal energy per unit mass, and net cooling rate per unit volume, respectively.

The adiabatic exponent $`\gamma _{\mathrm{ad}}`$ is given by

$$\gamma _{\mathrm{ad}}=1+\underset{i=1}{\overset{8}{}}y_i/\underset{i=1}{\overset{8}{}}\frac{y_i}{\gamma _i1},$$
(5)   
where the summation is over all chemical species, which include $`i=\mathrm{H}^0`$,$`\mathrm{H}^+`$,$`\mathrm{H}^{}`$,$`\mathrm{H}_2`$,$`\mathrm{He}^0`$,$`\mathrm{He}^+`$,$`\mathrm{He}^{++}`$, and $`e^{}`$. The adiabatic exponent for H<sub>2</sub> is given by

$$\frac{1}{\gamma _{\mathrm{H}_2}1}=\frac{1}{2}[5+2x^2\frac{e^x}{(e^x1)^2}],x\frac{6100\mathrm{K}}{T}$$
(6)   
for H<sub>2</sub>, where the second term in equation (6) counts the vibrational degrees of freedom and by $`1/(\gamma _i1)`$=1.5 for the other species. The concentration of the $`i`$ th species $`y_i`$ is defined by

$$y_i=\frac{n(i)}{n_\mathrm{H}},$$
(7)   
where $`n(i)`$ is the number density of the $`i`$ th species, and $`n_\mathrm{H}=n(\mathrm{H}^0)+n(\mathrm{H}^+)+n(\mathrm{H}^{})+2n(\mathrm{H}_2)`$ is the number density of hydrogen nuclei. We do not include the radiation force term in the momentum equation (2), since it is always negligible compared with the gas pressure in our runs.

The net cooling rate $`\mathrm{\Lambda }`$ consists of two parts, that is, the radiative cooling $`\mathrm{\Lambda }_{\mathrm{rad}}`$ and the chemical cooling $`\mathrm{\Lambda }_{\mathrm{chem}}`$. The former can be written by using the luminosity $`L(m)`$ as

$$\mathrm{\Lambda }_{\mathrm{rad}}(m)=\rho \frac{L(m)}{m}.$$
(8)   
The luminosity $`L(m)`$ is obtained by solving radiative transfer as described in §2.3 and the Appendix. We consider H<sub>2</sub> lines and continuum components that consist of H<sub>2</sub> collision-induced absorption (CIA), H<sup>-</sup> bound-free absorption, Lyman continuum absorption, etc., as sources of opacity. The net chemical cooling rate $`\mathrm{\Lambda }_{\mathrm{chem}}`$ is given by

$$\mathrm{\Lambda }_{\mathrm{chem}}(m)=\rho \frac{ϵ_{\mathrm{chem}}}{t},$$
(9)   
where $`ϵ_{\mathrm{chem}}`$ is the chemical binding energy per unit mass.

We keep the external pressure constant with time as the boundary condition. The timestep is the smaller of that determined by the Courant condition (Courant number=0.2) and $`10^3`$ of the local cooling time.

### 2.2 Chemical Reactions

We consider the following eight species: $`\mathrm{H}^0`$,$`\mathrm{H}^+`$,$`\mathrm{H}^{}`$,$`\mathrm{H}_2`$,$`\mathrm{He}^0`$,$`\mathrm{He}^+`$,$`\mathrm{He}^{++}`$, and $`e^{}`$. The abundance of helium atoms is assumed to be 28 % by mass.

At low density ($`n_\mathrm{H}<10^{15}\mathrm{cm}^3`$), reaction equations between $`\mathrm{H}^0,\mathrm{H}^+,\mathrm{H}^{},\mathrm{H}_2`$, and $`e^{}`$ are solved using an implicit difference scheme. Considered reactions are the same as Palla et al. (1983), while some reaction coefficients are altered (Table 1). The H<sup>-</sup> fraction is assumed to be the equilibrium value. Helium atoms are considered to be all neutral, because in our case this low density range corresponds to low temperature, namely $`T<2100\mathrm{K}`$.

Before the number density exceeds $`10^{15}\mathrm{cm}^3`$, the H<sub>2</sub> fraction easily reaches the equilibrium value. At $`n_\mathrm{H}>10^{15}\mathrm{cm}^3`$ we switch the computation to the equilibrium chemistry. Although the $`e^{}`$ fraction has not reached equilibrium yet, this does not cause any significant error, since $`e^{}`$ plays no role in H<sub>2</sub> formation/dissociation at that time.

The equilibrium chemistry ($`n_\mathrm{H}>10^{15}\mathrm{cm}^3`$) is computed by solving coupled Saha equations between $`\mathrm{H}^0,\mathrm{H}^+,\mathrm{H}_2,\mathrm{He}^0,\mathrm{He}^+,\mathrm{He}^{++}`$, and $`e^{}`$.

### 2.3 Radiative Transfer

We consider H<sub>2</sub> rovibrational lines and continuum components as sources of opacity and treat their transfer separately.

H<sub>2</sub> line transfer is dealt with as a frequency dependent transfer problem as is described in detail in the Appendix. We consider the first three vibrational states with rotational levels up to $`J=20`$ in each vibrational state following Palla et al. (1983). The population of each level is determined by assuming local thermodynamic equilibrium (LTE), since we consider the density range well above the critical number density $`n_{\mathrm{cr}}10^4\mathrm{cm}^3`$, where radiative and collisional deexicitation rates become equal for the molecular hydrogen. We consider both thermal and kinematic Doppler shifts as a line broadening mechanism.

Because of the lack of dust grains, primordial gas has only very weak continuum opacity. Continuum transfer is treated as a gray problem using mean opacity. When the cloud becomes substantially optically thick, it becomes difficult to keep numerical accuracy in solving the gray transfer equation. In such a case, we separate the cloud into two regions, namely, the opaque interior and transparent envelope. We solve the transfer equation only in the outer region, while the diffusion approximation is employed in the interior.

As the mean opacity, we take the continuum-Planck mean opacity for metal free gas at 1000K$`<T<`$7000K from Lenzuni, Chernoff & Salpeter (1991), which includes all the important continuum processes, specifically, bound-free absorption by $`\mathrm{H}^0`$ and $`\mathrm{H}^{}`$, free-free absorption by $`\mathrm{H}^0`$, $`\mathrm{H}^{}`$, $`\mathrm{H}_2`$, $`\mathrm{H}_2^{}`$, $`\mathrm{H}_2^+`$, $`\mathrm{H}_3`$, $`\mathrm{He}^0`$, and $`\mathrm{He}^{}`$, photodissociation of $`\mathrm{H}_2`$, and $`\mathrm{H}_2^+`$, Rayleigh scattering by $`\mathrm{H}^0`$, $`\mathrm{H}_2`$, and $`\mathrm{He}^0`$, Thomson scattering by $`e^{}`$, and collision-induced absorption by $`\mathrm{H}_2`$ due to collisions with $`\mathrm{H}_2`$, $`\mathrm{He}^0`$, and $`\mathrm{H}^0`$. They omit all lines, for example, Ly$`\alpha `$ or $`\mathrm{H}_2`$ rovibrational (quadrupole) transitions in their opacity. We treat the $`\mathrm{H}_2`$ lines as was described above. The Ly$`\alpha `$ lines are not so important in cooling in our case, because at such high temperatures where they are excited, the density is so high and their optical depth is so huge that these photons are always absorbed before escaping the cloud. At $`T<1000\mathrm{K}`$, we set the continuum opacity to be zero. At $`T>7000\mathrm{K}`$, it is supplemented by using the Rosseland mean opacity for metal free gas from Rogers & Iglesias (1992). Our adoption of the Planck mean opacity at $`T<7000\mathrm{K}`$ and the Rosseland mean opacity at higher temperature may deserve some comments. At low temperature, and then at low density, protostellar clouds are transparent to continuum emission, so that Planck mean opacity is appropriate. While the continuum optical depth exceeds unity $`(T>2000\mathrm{K})`$, both mean opacities have close values within a factor of 2 in our relevant range (Lenzuni et al. 1991). Actually, they match smoothly to each other at the interface.

We solve H<sub>2</sub> line transfer and continuum transfer independently, i.e., In solving H<sub>2</sub> line transfer, we neglect continuum opacity. However, errors caused by this are not significant. The reason is as follows: once the clouds become opaque to continuum opacity in a small central region where all the evolution occurs essentially, H<sub>2</sub> line cooling is negligible compared with continuum and chemical cooling there. There still exists the outer region where H<sub>2</sub> line radiation is the dominant cooling agent. However, in such a region the evolutionary timescale is so long in comparison with the central one that this region remains practically unchanged up to the formation of the hydrostatic core at the center (i.e., the end of our calculations). After the continuum optical depth exceeds unity, we set the cooling rate due to H<sub>2</sub> lines equal to zero.

## 3 Results

### 3.1 Initial Conditions

We consider primordial protostellar clouds which are fragments of larger ones. Fragmentation of a collapsing cosmological object into primordial protostellar clouds can be viewed as follows. It is known that a spherical cloud in pressure-free collapse is unstable with respect to non-spherical perturbations (Lin, Mestel, & Shu 1965; Hutchins 1976; Susa, Uehara, & Nishi 1996). As a consequence, the collapsing cosmological object first becomes a pancaking disk, which is gravitationally unstable and fragments more easily into filamentary clouds rather than into spherical clouds (Miyama, Narita, & Hayashi 1987a, 1987b). Such a filamentary cloud is also gravitationally unstable and subject to fragmentation but does not fragment immediately, as long as the filamentary cloud collapses so fast that the density perturbation has not sufficient time to grow (Inutsuka & Miyama 1992, 1997). Because the virial temperature of a filamentary cloud is constant with time, the free-fall time $`t_{\mathrm{ff}}`$ decreases faster than the cooling time $`t_{\mathrm{cool}}`$ as contraction proceeds. Eventually, $`t_{\mathrm{ff}}t_{\mathrm{cool}}`$, i.e., the collapse becomes slow, then the filamentary cloud fragments.

At the time of fragmentation, fragments cannot be highly gravitationally unstable. Hence, we take hydrostatic equilibrium clouds with adiabatic stratification as the initial conditions for simplicity. Namely, the density distribution is represented by $`\rho =\rho _c(\theta _n)^n`$, where $`\rho _c`$ is the central density, $`n=1/(\gamma _{ad}1)`$ is the polytropic index, and $`\theta _n`$ is the Emden function of index $`n`$. The pressure distribution is given by $`p=K\rho ^{\gamma _{ad}}`$, where $`K`$ is a free parameter fixed by giving the number density at half mass radius $`n_h`$ and the total mass of the cloud $`M`$. These clouds are cut off at the radius $`r_s`$ where the density falls off $`10^3`$ from the central value. Then parameters needed to specify an initial state are the mass of the cloud $`M`$, the number density at the half mass radius $`n_h`$, and the concentrations of chemical species. The initial parameters for several models are given in Table 2. Here, $`N`$ is the number of radial grids. As shown in Table 2, the initial width of the radial zones is increased by constant ratio for each zone as distance from the center is increased. This ratio is determined so as to resolve the smallest structure found in test runs with a sufficient number of radial grids. The mass and initial parameters of run A are close to those of a primordial cylindrical cloud at the fragmentation epoch (Uehara 1998). Hereafter, we call run A the fiducial run. The parameters of run B are the same as those of run A, except the number of radial grids is reduced to 100 in order to check the accuracy of our calcualtion. The match of the evolutionary paths between run A and run B is fairly good. Therefore, radial grids of 100 cells are enough to follow the collapse as long as the resolution around the center is sufficient. Runs C and D correspond to smaller and denser fragments than the fiducial run, which result from the later fragmentation of cylindrical clouds (Uehara et al. 1996). The evolution in these runs becomes almost the same as that of the corresponding inner mass regions of the fiducial run. If the H<sub>2</sub> fraction is too small at the begining, heat generated by H<sub>2</sub> formation exceeds the radiative loss by H<sub>2</sub> lines and the cloud expands until a sufficient amount of H<sub>2</sub> is produced. However, this expansion is clearly a result of inappropriate initial conditions. We put as much H<sub>2</sub> fraction as necessary in these runs to avoid such expansion. In the above runs, we set the H<sub>2</sub> concentration $`f_{\mathrm{H}_2}=\mathrm{const}.`$ with radius initially, while in reality the $`f_{\mathrm{H}_2}`$ distribution has a gradient, that is, higher in the interior and lower in the envelope. To see its effect, we set $`f_{\mathrm{H}_2}=0`$ initially in run E. In this case the cloud does not collapse dynamically until a sufficient amount of H<sub>2</sub> is formed. It is found that, however, the evolution thereafter is the same as in the fiducial case.

As was mentioned above, the structure of collapsing primordial clouds becomes almost the same after initial transient phases that merely reflect inappropriate initial conditions. Therefore, we describe the evolution in the fiducial run in detail hereafter.

### 3.2 The collapse of a primordial protostellar cloud

Initially the cloud is in hydrostatic equilibrium, and the cooling time is somewhat longer than the free-fall time. For the fiducial run the cooling time $`t_{\mathrm{cool}}\rho ϵ/\mathrm{\Lambda }=5.1\times 10^5`$ yr and the free-fall time $`t_{\mathrm{ff}}(3\pi /32G\rho )^{1/2}=3.2\times 10^5`$ yr at the center of the cloud initially. After some quasi-static contraction, the cooling time at the center becomes shorter than the free-fall time. At this time, the cloud begins to collapse dynamically in a free-fall time scale. In reality, the evolution should only be physically relevant thereafter, since a fragment is in the state where $`t_{\mathrm{cool}}t_{\mathrm{ff}}`$ at the fragmentation epoch (Low & Lynden-Bell 1976; Rees 1976; Silk 1977; Uehara et al. 1996).

The evolutionary sequences of number density, temperature, velocity, and H<sub>2</sub> concentration distributions are illustrated for the fiducial run in Figure 1 and for run C, in Figure 3, as a function of the radial distance. As we can see from Figures 1 a, the collapse proceeds almost self-similarly and is very reminiscent of Larson-Penston’s similarity solution (Penston 1969; Larson 1969). Since the free-fall time depends inversely on the density, the collapse proceeds most rapidly in the center where the density is the highest, leaving the outer part practically unchanged. The dash in the figure indicates the slope of density gradient for the isothermal case, namely $`\rho r^2`$ (Bodenheimer & Sweigart 1968). The slope in our case is slightly steeper, and its value is about $`2.2`$. For the Larson-Penston similarity solution in the case of a polytropic equation of state, i.e., $`p=K\rho ^\gamma `$, the density gradient may be written as (Larson 1969)

$$\frac{\mathrm{ln}\rho }{\mathrm{ln}r}=\frac{2}{2\gamma }.$$
(10)   
Therefore $`\gamma 1.1`$ in our case (see §3.3).

The evolution of the central temperature and the ratio of specific heats $`\mathrm{\Gamma }=\frac{\mathrm{ln}p/t}{\mathrm{ln}\rho /t}`$ at the center are illustrated in Figure 3 a and 3 b. According to Figure 3 b, the central value of $`\mathrm{\Gamma }`$ is about 1.1 in a large range of density, which coincides with the value expected from the density gradient. An initial low value of $`\mathrm{\Gamma }`$ at the center does not mean that the central part of the cloud cools efficiently at that time. The reason is the following. In quasi-static contraction, although $`\mathrm{\Gamma }`$ becomes nearly equal to the critical value 4/3 on average over a whole cloud, this needs not be $`\mathrm{\Gamma }4/3`$ at the center. Since the inner region cools faster than the outer part, a specific entropy gradient (i.e., $`\frac{s}{r}>0`$) is established in the course of contraction. In the case that the cloud contracts from the initial adiabatic stratification (i.e., $`\frac{s}{r}=0`$ ) to such a state, $`\mathrm{\Gamma }`$ tends to become lower in the inner region, while larger than 4/3 in the outer region in quasi-static contraction (see Omukai et al. 1998). Our initial low value of $`\mathrm{\Gamma }`$ at the center merely reflects this fact. After the cloud starts dynamical collapse, efficient H<sub>2</sub>-line cooling makes the ratio of specific heats $`\mathrm{\Gamma }`$ continue to decrease toward unity (corresponding to the isothermal case) on average, until the cloud becomes optically thick to several lines when the central number density $`n_\mathrm{c}10^{11}\mathrm{cm}^3`$.

Cooling due to molecular hydrogen formed by efficient three-body processes, which become the dominant mode of H<sub>2</sub> production for number density $`n>10^8\mathrm{cm}^3`$, makes the temperature in the central region even lower than in the outer part, as we can see in Figure 1 (b). These processes completely transform hydrogen in the central part of about $`1M_{\mathrm{}}`$ into molecular form before the temperature becomes so high that efficient H<sub>2</sub> dissociation begins. However, the outer region, i.e., $`m>\mathrm{several}M_{\mathrm{}}`$ in mass coordinate, remains totally atomic to the end (Figs. 1 and 2 d).

By the time the central number density reaches about $`10^{11}\mathrm{cm}^3`$ and the central temperature $`T_c900\mathrm{K}`$, more than 60 % of hydrogen atoms have been converted into molecular form in a small region around the center, and some of the H<sub>2</sub> lines begin to become optically thick. Eventually the cooling rate per unit mass around the center becomes lower than that in the outer region. Although this raises the ratio of specific heats $`\mathrm{\Gamma }`$ steadily, it does not exceed the critical value 4/3, as we can see in Figure 3 (b). The reason is that there are always sufficient lines that are efficient in cooling, namely, those whose optical depth is nearly unity.

When the central number density reaches about $`3\times 10^{13}\mathrm{cm}^3`$ and the central temperature about $`1600\mathrm{K}`$, the hydrogen molecules begin to dissociate gradually. However, this is only temporary. At the same time, H<sub>2</sub> CIA-continuum cooling comes into play. Continuum cooling is so strong as to stop and even to reverse dissociation, until it becomes optically thick and radiative cooling is no longer efficient. After a small central part of the cloud becomes opaque ($`3\times 10^{16}\mathrm{cm}^3,2000\mathrm{K}`$), the full-scale dissociation begins subsequently.

In the case of the star formation in present-day molecular clouds, after the central region becomes opaque to dust continuum, a transient core in hydrostatic equilibrium is formed with central density and temperature $`n_\mathrm{c}=10^{14}\mathrm{cm}^3`$, and $`T_\mathrm{c}=170\mathrm{K}`$ immediately after the formation (Larson 1969). In the primordial case, no transient core is formed when the central part of the cloud becomes optically thick to continuous radiation, since H<sub>2</sub> dissociation follows successively. The density and temperature at the center at the begining of the full-scale dissociation are nearly the same as in the present-day case. Then the evolution also becomes nearly the same.

When most of the hydrogen molecules are dissociated, the ratio of specific heats $`\mathrm{\Gamma }`$ rises above the critical value 4/3 (Fig. 3 b). The partial ionization of the hydrogen atoms reduces $`\mathrm{\Gamma }`$ temporarily but not below 4/3, since the density is already sufficiently high for compressional work to dominate cooling owing to ionization. After the central part of the cloud contracts almost adiabatically to some extent, the hydrostatic core, whose mass is about $`5\times 10^3\mathrm{M}_{\mathrm{}}`$, forms at the center when the central number density reaches $`10^{22}\mathrm{cm}^3`$ and the central temperature $`3\times 10^4\mathrm{K}`$. These physical dimensions are almost the same as those of the stellar core found in the calculation of present-day star formation. The formation of the core can be seen more clearly in the velocity distributions (Figs. 1 c and 2 c). When the central hydrostatic core forms, a shock front develops at the surface. The core grows in mass by accretion of envelope material and eventually will become an ordinary star.

### 3.3 The Mass Accretion Rate

Because of our explicit method of hydrodynamics, the Courant condition makes time steps extremely short and it becomes numerically very costly to pursue the evolution after the core formation. Moreover our artificial viscosity method smears the accretion shock front, thus producing a small temperature gradient, which results in a low luminosity and low radiative cooling when we use the diffusion approximation (Winkler & Newman 1980). We then stop our calculation at some arbitrary time after the core formation. Therefore the evolution thereafter cannot be known exactly.

However, we found that the Larson-Penston-type similarity solution for the equation of state $`p=K\rho ^\gamma `$ (Yahil 1983; Suto & Silk 1988), where $`K=4.2\times 10^{11}`$ (in cgs) and $`\gamma =1.09`$ reproduces well our results of the evolution in the fiducial run before the core formation (see Figure 4). Here we assume that the evolution of envelope after the central core formation is also described by the same similarity solution. Larson-Penston type similarity solutions can be extended naturally even after the central density becomes infinite (i.e., core formation) (Hunter 1977; Yahil 1983). The solutions after the central density becomes infinite have the finite point mass at the center (i.e., protostar), which grows in mass with time as (Yahil 1983; Suto & Silk 1988)

$`M_{}`$ $`=`$ $`K^{3/2}G^{(13\gamma )/2}t^{43\gamma }m_0`$ (11)   
$`=`$ $`0.11M_{\mathrm{}}({\displaystyle \frac{t}{1\mathrm{y}\mathrm{r}}})^{0.73},`$   
where $`M_{}`$ is the mass of the central protostar, $`t`$ is the time elapsed since the core formation, and $`m_0`$ is a non-dimensional constant whose value is 20 for $`\gamma =1.09`$ Larson-Penston type solution. We used $`K=4.2\times 10^{11}`$ and $`\gamma =1.09`$ in the second expression of equation (11). Thus the mass accretion rate to the central protostar is given by

$`\dot{M}_{}`$ $`=`$ $`(43\gamma )K^{3/2}G^{(13\gamma )/2}t^{3(\gamma 1)}m_0`$ (12)   
$`=`$ $`8.3\times 10^2M_{\mathrm{}}\mathrm{yr}^1({\displaystyle \frac{t}{1\mathrm{y}\mathrm{r}}})^{0.27}.`$   
Using the relation in equation (11), we can write the mass accretion rate as a function of an instantaneous mass of the protostar as

$$\dot{M}_{}=3.7\times 10^2M_{\mathrm{}}\mathrm{yr}^1(\frac{M_{}}{M_{\mathrm{}}})^{0.37}.$$
(13)   
According to equation (12) or (13), the mass accretion rate is huge and is diminishing with time. From the first expression of equation (12), we can see that there are three parameters that determine the mass accretion rate, i.e., $`\gamma ,K`$ and $`m_0`$. The decrease of the accretion rate with time is the common nature of self-similar collapse with $`\gamma >1`$ as can be seen by the first expression of equation (12). The parameter $`K`$ is related to the temperature of protostellar clouds, and the parameter $`m_0`$ concerns the type of the similarity solution (i.e., the Larson-Penston-type or the Shu-type solution, etc.) as well as $`\gamma `$. Stahler et al. (1986) pointed out the higher mass accretion rate $`4.41\times 10^3M_{\mathrm{}}\mathrm{yr}^1`$ of primordial protostars than present-day value because of higher temperature (then the higher $`K`$ value) in primordial protostellar clouds, although they thought of Shu-like collapse. In our analysis, however, the collapse proceeds like the Larson-Penston-type solution, whose parameter $`m_0`$ is about an order of magnitude higher than that of the Shu-type solution. For example, in the isothermal ($`\gamma =1`$) Larson-Penston collapse the density is 4.4 times higher than the hydrostatic equilibrium value and the fluid velocity is 3.3 times than the sound speed at the time of central core formation. This results in the value of $`m_0`$ and the mass accretion rate 48 times higher than that in isothermal Shu collapse of the same temperature (Hunter 1977).

## 4 Conclusion

In this paper we have investigated the collapse of primordial clouds into protostars using a spherically symmetric Lagrangian hydrodynamic code combined with radiative transfer of H<sub>2</sub> lines/continuum and chemical reactions.

The collapse proceeds self-similarly like the Larson-Penston similarity solution until the central region reaches stellar density. During the collapse, the central part of about $`1M_{\mathrm{}}`$ becomes fully molecular before the dissociation, while the outer part ($`m>\mathrm{several}M_{\mathrm{}}`$ in mass coordinate) remains almost fully atomic. When the small central part of the cloud reaches stellar densities ($`10^2\mathrm{g}\mathrm{cm}^3`$), the hydrostatic core, whose physical dimensions are roughly the same as that in present-day star formation, is formed at the center of the cloud. The mass accretion rate is typically $`10^2M_{\mathrm{}}\mathrm{yr}^1`$, which is about 3 orders of magnitude higher than in the present-day case, and it declines with time.

Our results suggest a stellar core of $`5\times 10^3M_{\mathrm{}}`$ at its formation. On the other hand, Palla et al. (1983) argued the minimum Jeans mass in the primordial gas clouds to be $`0.05M_{\mathrm{}}`$ at H ionization. This mass can be identified with the mass of the stellar core at the time of formation (Carlberg 1981). This descrepancy seems to come not only from Palla et al. (1983)’s one-zone treatment, but also from the fact that Palla et al. (1983) did not take into account H<sub>2</sub> CIA continuum.

We are grateful to H. Masunaga, H. Susa, and H. Uehara for fruitful discussions, to H. Sato for his continuous encouragement, to T. Chiba and N. Sugiyama for reading our manuscript carefully, to S. Hayward for checking the English, and to the referee, F. Palla, for improving this manuscript. We also would like to thank participants in the workshop “Molecular Hydrogen in the Early Universe” (Firenze, 1997 December 4-6) for invaluable comments. We owe thanks to the YITP computer system for the numerical analyses. This work is supported in part by a Grant-in-Aid of Scientific Research from the Ministry of Education, Science, Sports and Culture No. 09740174(RN).     

## Appendix A Radiative Transfer in Spherical Symmetry

In this Appendix, we describe our schemes of calculations of radiative transfer. The specific intensity $`I_\nu (\mathrm{ergs}\mathrm{sec}^1\mathrm{cm}^2\mathrm{sr}^1\mathrm{Hz}^1)`$ along a ray is calculated by solving the radiative transfer equation (e.g., Rybicki & Lightman 1979),

$$\frac{dI_\nu }{ds}=\alpha _\nu (S_\nu I_\nu ),$$
(A1)   
where $`s`$ is the displacement along the ray, $`\alpha _\nu `$ is the absorption coefficient, and $`S_\nu `$ is the source function. If the matter is in LTE, $`S_\nu =B_\nu (T)`$.

In a dynamically collapsing cloud, supersonic fluid velocity develops, then Doppler shift owing to bulk motion becomes comparable with thermal Doppler width, and we have to take it into account. For H<sub>2</sub> lines, absorption coefficients can be written by using Einstein $`B`$-coefficients,

$$\alpha _\nu =\frac{h\nu _{co}}{4\pi }\varphi (\nu _{co})(n_1B_{12}n_2B_{21}).$$
(A2)   
Here, $`\nu _{co}=\nu (1v_t/c)`$ is the frequency in the comoving frame, $`v_t`$ is the bulk velocity tangent to the ray, $`n_1`$ and $`n_2`$ are the number densities of molecules in the lower and upper level of the transition, respectively, and $`\varphi (\nu )`$ is the line profile function. We only consider line broadening due to thermal Doppler effects. Then

$$\varphi (\nu )=\frac{1}{\mathrm{\Delta }\nu _\mathrm{D}\sqrt{\pi }}e^{(\nu \nu _0)^2/(\mathrm{\Delta }\nu _\mathrm{D})^2}.$$
(A3)   
Here, $`\nu _0`$ is the (comoving) line-center frequency and the Doppler width $`\mathrm{\Delta }\nu _\mathrm{D}`$ is defined by

$$\mathrm{\Delta }\nu _\mathrm{D}=\frac{\nu _0}{c}\sqrt{\frac{2kT}{\mu m_\mathrm{H}}},$$
(A4)   
where $`T`$ is the temperature at each radius in the cloud, $`m_\mathrm{H}`$ is the mass of a hydrogen atom, and $`\mu =2`$ is the molecular weight for the molecular hydrogen. We include the first three vibrational states, with rotational levels up to $`j=20`$ in each vibrational state following Palla et al. (1983). We use $`A`$-coefficients of molecular hydrogen from Turner, Kirby-Docken, & Dalgarno (1977). B-coefficients can be known from A-coefficients by the Einstein relations

$$g_1B_{12}=g_2B_{21},$$
(A5)   
and

$$A_{21}=\frac{2h\nu _{0}^{}{}_{}{}^{3}}{c^2}B_{21},$$
(A6)   
where $`g_i`$ is the statistical weight of the level $`i`$. The critical number density $`n_{\mathrm{cr}}`$, where radiative and collisional de-excitation rates become equal for the molecular hydrogen, is $`10^4\mathrm{cm}^3`$. In our calculations, the number density is well above the critical value, so that almost all of excited molecules would be deexcited by collision with other molecules. Therefore, collisional LTE is established, and the level populations are easily determined.

Figure 5 shows the geometry of grids of radiative transfer calculation in the case of spherical symmetry (Hummer & Rybicki 1971; see also Masunaga, Miyama, & Inutsuka 1998). To know the luminosity at the $`i`$-th radial mesh, $`r_i`$, we solve the transfer equation (A1) at a frequency $`\nu `$ along from the 0-th to $`(i1)`$-th rays, and obtain the specific intensity at $`r_i`$ in the direction $`\theta _{ij}`$, $`I_\nu (r_i,\theta _{ij})(i+1ji1)`$. This procedure will be repeated for each frequency mesh. We consider the frequency range for each line within 6 times the Doppler width corresponding to the central temperature from the line center of static media. Each line is separated into 30 frequency meshes in our calculations, which is indicated to be sufficient by some experiments. Integrating over $`\theta `$ and $`\nu `$, we obtain the monochromatic energy flux $`F_\nu `$ and the luminosity $`L`$:

$$F_\nu (r)=2\pi _0^\pi I_\nu (r,\theta )\mathrm{cos}\theta 𝑑\theta ,$$
(A7)   
$$L(r)=4\pi r^2_0^{\mathrm{}}F_\nu (r)𝑑\nu .$$
(A8)   
The cooling rate by lines whose optical depth at line center is less than 0.1 is calculated by using an optically thin approximation.

We approximate continuum transfer as an grey transfer problem. For grey material $`\alpha _\nu =\alpha `$, integrating over frequency we have

$$\frac{dI}{ds}=\alpha (SI),$$
(A9)   
where

$$I_0^{\mathrm{}}I_\nu 𝑑\nu ,S_0^{\mathrm{}}S_\nu 𝑑\nu .$$
(A10)   
For LTE matter, $`S=B(T)=\frac{ac}{4\pi }T^4`$. We solve the frequency integrated transfer equation (A9) along the same rays as in Figure 5, and integrating over the polar angle $`\theta `$ we obtain the continuum luminosity.

As the cloud contracts, the increment of optical depth between adjacent radial grids $`\mathrm{\Delta }\tau `$ around the center exceeds unity. In such a case, we can not resolve the temperature gradient precisely, and the numerical accuracy may be severely damaged in solving the gray transfer equation. Then we treat the inner opaque region and the outer transparent region separately. In the inner region the diffusion equation is solved, i.e., the continuum energy flux is given by

$$F_{\mathrm{cont}}=\frac{16\sigma T^3}{3\alpha }\frac{T}{r},$$
(A11)   
and in the outer region the transfer equation (A9) is solved as described above. These solutions for continuum energy flux are matched at the interface, which is chosen so as to satisfy the following conditions (Narita, Nakano, & Hayashi 1970):

$$|\frac{B(r)I(r,\pi )}{I(r,\pi )}|<0.2,$$
(A12)   
$$\tau (r)>4,$$
(A13)   
and

$$\mathrm{\Delta }\tau >1,$$
(A14)   
where $`I(r,\pi )`$ is the inward intensity, and $`\tau `$ is the optical depth that is measured from the surface in the radial direction. Under these conditions, the error in $`F_{\mathrm{cont}}(r)`$ is less than 4 percent.