text:
# Very high column density and small reddening towards GRB 020124 at z = 3.201footnote 11footnote 1Based on observations with the Nordic Optical Telescope, which is operated on the island of La Palma jointly by Denmark, Finland, Iceland, Norway, and Sweden, at the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofísica de Canarias, and on observations collected by the Gamma-Ray Burst Collaboration at ESO (GRACE) at the European Southern Observatory, Paranal, Chile (ESO Large Programme 165.H–0464)

## 1 Introduction

Spectroscopy of the optical afterglows of cosmological GRBs allows detailed studies of the chemical and kinematical properties of gas along the lines of sight. Independently of the brightness of the host galaxy, the resulting insight into the GRB environment and its dust and gas content can provide clues to the nature of GRB progenitors as well as the properties of high-redshift galaxies. In particular, constraints on the column density of neutral hydrogen, N(H I), can be obtained through Ly$`\alpha `$ absorption if the burst is sufficiently distant to redshift the Ly$`\alpha `$ line into the near-UV/optical domain. Unlike X-ray studies where one infers N(H I) only indirectly based on an assumed metallicity, the derived N(H I) will be independent of metallicity.

In this regard, spectroscopy of GRB optical afterglows bears some resemblance to that of Damped Ly$`\alpha `$ Absorbers (DLAs, Wolfe et al. (1986)), which are gas-rich absorption systems intervening the lines of sight to background sources, such as QSOs. DLAs are chemically enriched and are therefore most likely caused by gas in, or close to, galaxies (e.g., Lu et al. (1996) and references therein; Pettini et al. (1997a)). In some cases galaxy counterparts of DLAs have been detected directly (Møller & Warren, 1993, 1998; Djorgovski et al., 1996; Fynbo, Møller & Warren, 1999; Møller et al., 2002). The important and interesting difference between QSO DLAs and GRB absorption systems is that the lines of sight through QSO DLAs are absorption cross-section selected whereas GRB absorbers (at the GRB redshift) are ‘GRB progenitor site’ selected. Therefore, GRB absorbers are likely to probe sightlines through different parts of DLA galaxies than QSO absorbers do, and hence contribute to a more complete picture of the properties of high-redshift galaxies.

The dust content of DLAs has been constrained mainly by two methods: i) by comparing the colors of background QSOs having intervening DLAs with the colors of a control sample of QSOs without intervening DLAs (Pei et al., 1991), and ii) by studying the metal abundance ratios, most notably Zn vs. Si and iron group elements, and inferring the dust depletion by comparison with interstellar clouds in the Galaxy (Lu et al., 1996; Kulkarni, Fall & Truran, 1997; Pettini et al., 1997b; Ledoux, Bergeron & Petitjean, 2002). The main conclusion of these studies is that DLAs in general are metal poor and that the fraction of the metals bound in dust grains is roughly the same as, or somewhat lower than, that found in the Galaxy.

This paper presents optical/near-IR photometry and optical spectroscopy of the afterglow of GRB 020124 (catalog ). GRB 020124 (catalog ) (burst trigger on Jan 24.44531 2002 UT) was the second HETE-II burst (Ricker et al., 2002) to be localized to arcsec precision and is the only burst for which a host galaxy has not been detected when searched for with HST (Berger et al., 2002; Bloom, Kulkarni & Djorgovski, 2002). We determine the redshift of the burst to be 3.20, which is the fourth highest redshift known for a GRB to date, and the third highest (after GRB 000131 (catalog ) at $`z=4.50`$ (Andersen et al., 2000) and GRB 030323 (catalog ) at $`z=3.37`$ (Vreeswijk et al., 2003a)) based on afterglow absorption lines rather than host galaxy emission lines. We also show that it has a very high column density, qualifying it as a DLA (the highest value of N(H I) of known GRB afterglows and higher than almost all QSO DLAs). We discuss the implications for the dust-to-gas ratio of the GRB surroundings, including the effects of metallicity and destruction of dust by the GRB itself.

We assume a cosmology where $`H_0=70`$ km s-1 Mpc-1, $`\mathrm{\Omega }_m=0.3`$, and $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$. For these parameters, a redshift of 3.20 corresponds to a luminosity distance of 27.47 Gpc and a distance modulus of 47.19. One arcsecond corresponds to 7.55 proper kpc, and the lookback time is 11.5 Gyr.

## 2 Observations

The observing log and photometric results are reported in Table 1 and summarized below.

### 2.1 Optical imaging

The position of the burst was observed at RIKEN with the 0.25 m f/3.4 hyperboloid astrograph equipped with unfiltered CCD camera AP7p. The field of view was $`50\mathrm{}\times 50\mathrm{}`$ which covered the entire $`12\mathrm{}`$ radius HETE-II error circle (Ricker et al., 2002). The observation started at Jan 24.530 UT; 126 frames of 20-s exposure were acquired by Jan 24.574 UT. PSF photometry was applied to each of the 126 frames. The resultant photometric measurements were combined to yield a 3.0$`\sigma `$ detection at the position of the afterglow reported by Berger et al. (2002).

R-band images of the optical afterglow were obtained with StanCam at the 2.56 m Nordic Optical Telescope (NOT) between Jan 26.035 and Jan 26.076 2002 UT. I-band images were obtained at VLT-Melipal (UT3) between Jan 26.311 and 26.318 UT. R-band images were obtained with FORS1 on VLT-Melipal between Jan 27.278 UT and 27.292 UT. The latter images have a pronounced background pattern due to instrument reflections of the very bright background (the moon was $`75`$ % illuminated). The variation in the pattern across the field is approximately 3–4 %. To flatten each individual image we used SExtractor to create a smoothed background image which was subsequently subtracted. The five background-subtracted images were then combined.

All images were preprocessed using standard tools and calibrated using nearby comparison stars from the field photometry of Henden (2002). The derived magnitudes are listed in Table 1.

### 2.2 Near-infrared imaging

Near-IR observations were obtained with ISAAC on VLT-Antu (UT1). The data consist of a total of 30 min in the Ks band and 30 min in the Js band obtained on Jan 26.2 UT and 30 min in the Ks band obtained on Jan 27.2 UT. The reduction was carried out with the “eclipse” software package (Devillard, 1997) and the IRAF “Experimental Deep Infrared Mosaicing Software” xdimsum. Eclipse was used to remove effects of electrical ghosts from science and calibration frames and to construct flat fields and bad-pixel maps from a series of twilight-sky flats. Sky subtraction and combination of the dithered science frames were carried out with xdimsum. To account for the rapidly varying background, the sky value in each pixel was calculated as the running median of the pixel value in the exposures taken immediately before and after a given exposure. To minimize the effects on the running median of bright objects in neighboring frames, the subtraction was carried out in two iterative steps: First, a cosmic-ray cleaned, background-subtracted combined image was constructed from which a mask of all the objects was created. This mask was then registered to the individual frames and taken into account in the following running median background subtraction. Residual background signatures were subsequently removed by fitting high-order polynomials to the rows and columns of the individual (masked) frames. Finally, the individual background-subtracted images were registered and median combined. Excerpts of the resulting images are shown in Fig. 1.

Photometric calibration to J and Ks magnitudes was performed using standard star observations of 9106/S301–D and 9149/S860–D (Persson et al., 1998) on the nights of observation. The images were astrometrically calibrated using comparison stars from the Guide Star Catalogue II (GSC-II)222The Guide Star Catalogue II is a joint project of the Space Telescope Science Institute and the Osservatorio Astronomico di Torino. Space Telescope Science Institute is operated by the Association of Universities for Research in Astronomy, for the National Aeronautics and Space Administration under contract NAS5-26555. The participation of the Osservatorio Astronomico di Torino is supported by the Italian Council for Research in Astronomy. Additional support is provided by European Southern Observatory, Space Telescope European Coordinating Facility, the International GEMINI project and the European Space Agency Astrophysics Division. using the WCStools software package by D. Mink333http://tdc-www.harvard.edu/TDC.html. The absolute astrometry has an rms of $`0\stackrel{}{\mathrm{.}}43`$ in the Ks band. The resulting coordinates of the infrared afterglow of GRB 020124 (catalog ) are R.A. (J2000.0) = 09h32m50.82s, Dec. (J2000.0) $`=11^{}31\mathrm{}11\stackrel{}{\mathrm{.}}0`$, consistent with the position of the radio and optical afterglow (Berger et al., 2002).

### 2.3 Spectroscopy

Spectroscopic observations were obtained on Jan 26.3 with FORS1 on VLT-Melipal using the 300V+10 grism and the order-separation filter GG375 in long-slit mode with a 1.0 arcsec wide slit. This configuration provides a wavelength coverage 3650 – 7500 Å at a resolution of 13 Å. The total integration time was 1200 sec, evenly divided between two exposures. Standard data reduction procedures were performed using MIDAS. The two separate exposures allowed us to reliably reject cosmic-ray events before co-addition. The spectrophotometric standard LTT9491 (Hamuy et al., 1994) was observed using the same configuration.

## 3 Results

### 3.1 Spectral energy distribution: spectral index and redshift estimate

The extracted, flux calibrated spectrum is shown in Fig. 2. Due to the faintness of the object spectrum and the resulting dominance of systematic errors from sky subtraction close to bright airglow lines, the signal-to-noise fluctuates rapidly in the red part of the spectrum. What may look like absorption or emission features redwards of 5500 Å should therefore not be trusted. The strong absorption line around 5100 Å is clearly significant, and several other absorption features bluewards of 5500 Å are marginally significant at the 3–4 $`\sigma `$ level.

One or more of the absorption lines seen in the blue part of Fig. 2 could be due to Ly$`\alpha `$, which would indicate a redshift in excess of 2.3. Our first objective was therefore to construct the spectral energy distribution (SED) in the optical to determine its slope, and to look for possible large-scale signatures (spectral drops) which might indicate the onset of the Lyman Forest and/or the Lyman Valley as expected in an object with redshift significantly in excess of 2.3. For this purpose we carefully defined a series of intervals along the spectrum where systematic errors from night-sky emission-line subtraction did not pose any problems. The total number of counts was determined in each interval of the afterglow spectrum as well as of the spectrum of the spectrophotometric standard. Uncertainties were calculated directly via propagation of photon statistics.

The resulting fluxes and corresponding uncertainties (after correction for airmass) in each bin are plotted in Fig. 3. In Fig. 3a the minimum-$`\chi ^2`$ pure power-law fit ($`f_\nu \nu ^\beta `$; $`\beta =2.36\pm 0.23`$; reported uncertainties are 1$`\sigma `$ errors throughout this paper) is overplotted, but is seen to be a very poor fit to the data points. In particular the data point at 5100 Å falls many standard deviations below any power-law fit. Such a large drop in that wide a bin can effectively only be caused by a damped Ly$`\alpha `$ line at a redshift of about 3.2. At this redshift one expects to see not only the Ly$`\alpha `$ line but also a significant drop due to the Lyman Forest, and the signature of the onset of the red slope of the Lyman Valley (Møller & Jakobsen, 1990). In Fig. 3b we therefore plot the same data points, but here we overlay a model power law including the predicted absorption due to the intervening intergalactic medium at $`z=3.2`$. The power law again represents the minimum-$`\chi ^2`$ fit and has an index of $`\beta =1.32\pm 0.25`$. Comparing the $`\chi ^2`$ of the fits in Fig. 3a and Fig. 3b, we find that even if we ignore the presence of the Ly$`\alpha `$ line (i.e. consider only the large-scale SED), the $`z=3.2`$ model SED represents an acceptable fit while the pure power law is rejected at the 98.6 % confidence level ($`\mathrm{\Delta }\chi ^2=16.4`$ for 20 degrees of freedom (dof)).

### 3.2 Redshift and H I column density of the Ly$`\alpha `$ line

Following the determination of the underlying power-law index and approximate redshift, we proceeded to fit the Ly$`\alpha `$ line itself. The part of the spectrum useful for line-fitting is shown again in Fig. 4, here normalized by division with the fitted power-law. Overplotted are model absorption spectra with $`z=3.198`$. The model spectra include the Lyman series, Si II, Si III, and O I lines. On the blue side of the Ly$`\alpha `$ line we see the expected line-blanketing by the Lyman Forest. At this resolution, lines in the Lyman Forest are therefore not useful for a precise redshift determination, and we shall rely only on the Si II $`\lambda 1260`$ and O I $`\lambda 1302`$ lines at 5294 Å and 5468 Å for this purpose. The first line is an unblended Si II $`\lambda 1260`$ line providing a best-fit redshift of $`3.200\pm 0.004`$. The O I $`\lambda 1302`$ line is blended with a much weaker Si II $`\lambda 1304`$ line and provides a redshift of $`3.196\pm 0.004`$. The combined redshift of these two lines, including wavelength calibration errors, is $`z_{\mathrm{metal}}=3.198\pm 0.004`$ which we shall adopt as the systemic redshift of the absorber.

Because of the line-blanketing on the blue side of the Ly$`\alpha `$ line, only the red part of the line profile could be used to constrain the H I column density of the absorber. Within the uncertainty of the redshift we found that acceptable fits could be obtained in the range $`\mathrm{log}(N_{HI})=21.7\pm 0.2`$. The estimated $`1\sigma `$ range is marked in Fig. 4 as dotted lines. This column density is among the very highest values observed for QSO DLAs (Storrie-Lombardi & Wolfe, 2000). The absorption system is unlikely to be “intervening” (i.e. unrelated to the GRB host) because of the very small probability of such a high-column-density intervening absorber. More importantly, a significantly higher redshift of the GRB is ruled out by the lack of even stronger Lyman Forest absorption redward of 5200 Å (Fig. 3).

### 3.3 Characteristics of the optical decay

In addition to the data presented here, we include the data reported by Berger et al. (2002) to constrain the R-band lightcurve. In Fig. Very high column density and small reddening towards GRB 020124 at z = 3.20111Based on observations with the Nordic Optical Telescope, which is operated on the island of La Palma jointly by Denmark, Finland, Iceland, Norway, and Sweden, at the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofísica de Canarias, and on observations collected by the Gamma-Ray Burst Collaboration at ESO (GRACE) at the European Southern Observatory, Paranal, Chile (ESO Large Programme 165.H–0464) we plot all the R-band data points vs. the logarithm of time since the GRB. The power-law decay index derived from a weighted fit to all the data is $`\alpha =1.64\pm 0.03`$ ($`f_\nu (tt_{\mathrm{GRB}})^\alpha `$). This fit is formally rejected with a reduced $`\chi ^2`$/dof = 26/14. Excluding the late-time HST data we derive an acceptable fit with a decay index of $`\alpha =1.49\pm 0.04`$ and $`\chi ^2`$/dof = 6/12. Fitting only to the data prior to 1 day after the burst we derive a decay index of $`\alpha =1.45\pm 0.06`$ with $`\chi ^2`$/dof = 5/10. These results confirm the conclusion of Berger et al. (2002) that the lightcurve is becoming progressively steeper (a ‘break’) (see Fig. Very high column density and small reddening towards GRB 020124 at z = 3.20111Based on observations with the Nordic Optical Telescope, which is operated on the island of La Palma jointly by Denmark, Finland, Iceland, Norway, and Sweden, at the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofísica de Canarias, and on observations collected by the Gamma-Ray Burst Collaboration at ESO (GRACE) at the European Southern Observatory, Paranal, Chile (ESO Large Programme 165.H–0464)).

### 3.4 Optical/near-infrared SED

Assuming a power-law decay index $`\alpha =1.49\pm 0.04`$ we can interpolate the RIJKs-band data taken around Jan. 26 UT to a common epoch at Jan 26.20 2002 UT. After dereddening the contemporaneous RIJKs-band points for Galactic extinction (Schlegel, Finkbeiner & Davis, 1998), the colors are $`RKs=2.93\pm 0.26`$ and $`JKs=1.33\pm 0.20`$. These colors are consistent with the ones expected of a GRB afterglow (see the shaded region of the color–color diagram in Fig. 2 of Gorosabel et al. (2002a)).

We have used the intrinsic spectral index obtained from the VLT spectrum ($`\beta =1.32\pm 0.25`$, see § 3.1) to obtain a mock B-band point by a power-law extrapolation of the R-band measurement. In order to infer information on the extinction ($`A_V`$, E(B$``$V)), the BRIJKs-band SED was fitted with a functional form $`f_\nu \nu ^\beta \times 10^{0.4A_\nu }`$, where $`A_\nu `$ is the extinction in magnitudes at rest-frame frequency $`\nu `$. $`A_\nu `$ was parametrized in terms of $`A_V`$ using the three extinction laws reported by Pei (1992), i.e., for the Small Magellanic Cloud (SMC), Large Magellanic Cloud (LMC) and for the Milky Way (MW). For comparison purposes we also considered the unextinguished case (pure power law spectrum given by $`f_\nu \nu ^\beta `$).

As shown in Fig. 6 and Table 2 (third column) the fits of the three extinction laws are consistent with the data. The derived extinctions are consistent with the no-extinction case (in which case $`\beta =0.91\pm 0.14`$). Based on these fits we can set an upper limit of $`\mathrm{E}(\mathrm{B}\mathrm{V})<0.18`$. This upper limit however corresponds to an unrealistic value of $`\beta `$. Constraining the spectral index to be $`\beta >0.5`$ and the extinction law to be that of the SMC (which provides the best fit to this and most other afterglow SEDs) we find that $`\mathrm{E}(\mathrm{B}\mathrm{V})<0.065`$. The $`\beta `$ values derived for the three extinguished cases agree with the pure power-law case.

## 4 Discussion

### 4.1 Jet–wind fireball model

Within the framework of the afterglow synchrotron model (Sari, Piran & Narayan, 1998), the decay index, $`\alpha `$, and the spectral index, $`\beta `$, are related through the slope of the electron energy distribution, $`p`$. Given the low extinction inferred from the SED analysis we find that a jet ($`\nu <\nu _c`$, the synchrotron cooling frequency) expanding into a wind medium ($`nr^2`$) provides the best fireball model for the data. For $`\alpha _1=1.49\pm 0.04`$, this model requires $`p=2.32\pm 0.05`$, and $`\beta =0.66\pm 0.03`$, consistent with the observational results with modest extinction. Expansion into a homogeneous medium ($`\nu <\nu _c`$) is not strictly ruled out but requires an unusually high value of $`p`$: the predicted values are $`p=2.99\pm 0.05`$ and $`\beta =1.00\pm 0.03`$ (for $`\alpha _1=1.49\pm 0.04`$), consistent with the observational values in the absence of extinction.

### 4.2 How robust is the redshift determination?

The spectrum has low resolution and signal-to-noise ratio, but there are five independent pieces of evidence for the redshift. First, assuming a pure power-law for the spectrum provides a spectral index of $`2.36\pm 0.23`$ which is incompatible with $`\beta =0.91\pm 0.14`$ determined from the photometry. Second, the strong absorption line at 5100 Å cannot be explained by anything else than Ly$`\alpha `$ at a redshift around 3.2. Third, assuming $`z=3.2`$ and fitting a model with free spectral index but now including the predicted Lyman Forest and Valley absorption, one obtains a minimum-$`\chi ^2`$ fit for a spectral index of $`1.32\pm 0.25`$ which is fully compatible with the photometric determination. Fourth, the $`z=3.2`$ fit leads to an improvement in $`\chi ^2`$ over a single power-law alone (for the same degrees of freedom) which rejects the single power-law at the 98.6 % level. Fifth, there are two identified metal lines, and none of the lines predicted from the explicit fit (e.g. Ly$`\beta `$) are inconsistent with the observed spectrum (Fig. 4).

### 4.3 Dim or bright?

The fact that GRB 020124 (catalog ) was a fairly dim burst (Berger et al., 2002), but fitted by a normal fireball model without excessive extinction reinforces the conclusion of Hjorth et al. (2002) (GRB 980613 (catalog )), Berger et al. (2002) (GRB 020124 (catalog )), and Fox et al. (2003) (GRB 021211 (catalog )) that there is a population of dim, unextinguished GRB afterglows that can account for at least some of the large fraction of (‘dark’) bursts for which no optical afterglow is detected (Fynbo et al., 2001a). In the case of GRB 020124 (catalog ) this is mainly due to the fairly high redshift (unlike GRB 980613 (catalog ) and GRB 021211 (catalog ) which are at redshifts 1.096 (Djorgovski et al., 1999) and 1.006 (Vreeswijk et al., 2003b), respectively): If an afterglow is redshifted from $`z_1`$ to $`z_2`$ it is dimmed by $`\mathrm{\Delta }\mu +2.5(\beta \alpha 1)\mathrm{log}\left((1+z_2)/(1+z_1)\right)`$. For a median GRB redshift of $`z_1=1`$ and $`z_2=3.198`$ the relative dimming amounts to 1.82 mag (for $`\alpha =1.49`$ and $`\beta =0.91`$). If the afterglow had been 1.82 mag brighter it would have been a typical bright afterglow: $`R16.7`$ after two hours and $`R20.7`$ after one day (see e.g. Fig. 3 of Gorosabel et al. (2002b) or Fig. 2 of Fox et al. (2003)).

### 4.4 N(H I) vs. E(B$``$V)

The Galactic relation between the column density of neutral hydrogen and the reddening is N(H I)/E(B$``$V) = $`4.93\pm 0.28\times 10^{21}`$ cm2 mag-1 (Diplas & Savage, 1994). For the LMC and SMC the corresponding ratios are N(H I)/E(B$``$V) = $`2\pm 0.5\times 10^{22}`$ cm2 mag-1 and N(H I)/E(B$``$V) = $`4.4\pm 0.7\times 10^{22}`$ cm2 mag-1, respectively (Koornneef, 1982; Bouchet et al., 1985). In Table 3 we list this ratio for known GRBs for which the H I column density and E(B$``$V) have been measured or constrained from optical spectroscopy and optical/near-infrared photometry (the values given are for an SMC extinction law). In addition to GRB 020124 (catalog ) these are GRB 000301C (catalog ) (Jensen et al., 2001) and GRB 000926 (catalog ) (Fynbo et al., 2001b, c). We exclude GRB 021004 (catalog ) which has several intervening absorbers at different redshifts. It is evident that the values for N(H I)/E(B$``$V) are not larger than expected for galaxies like the LMC or SMC. Taken at face value these results indicate that the GRB surroundings are low in dust, because of low metallicity (e.g., for a universal dust-to-metals ratio) and/or because of a low dust-to-metals ratio (a large fraction of the metals being in the gas phase). This conclusion is fully consistent with what is found in QSO DLAs (Pettini et al., 1997b).

It has been suggested that the intense UV and X-ray flux from a GRB may photoionize the surrounding gas and destroy dust grains (Waxman & Draine, 2000; Fruchter, Krolik & Rhoads, 2001; Perna & Lazzati, 2002), resulting in changes in the derived N(H I), $`A_V`$, E(B$``$V), and their ratios. Galama & Wijers (2001) studied the inferred absorption of soft X-rays in BeppoSAX GRB afterglows and compared it to the optical extinction as inferred from a fireball model fit to the available data. They claimed evidence for a high N(H I)/$`A_V`$ ratio and concluded that dust destruction may be the cause. The most convincing case for dust destruction along these lines was made by Galama et al. (2003) who found a high value of N(H I)/$`A_V`$ for GRB 010222 (catalog ), based on Chandra X-ray Observatory data. At the same time they interpreted a strong (i.e., non-grey) far-ultraviolet (restframe wavelength $`1100`$ Å) component in the extinction of GRB 010222 (catalog ) as evidence for dust destruction.

In a comparison of the properties of DLAs and GRB absorbers, Savaglio, Fall & Fiore (2003) claimed that GRB 990123 (catalog ) ($`z=1.60`$), GRB 000926 (catalog ) ($`z=2.04`$), and GRB 010222 (catalog ) ($`z=1.48`$) (of which only GRB 000926 (catalog ) has a measured N(H I) (Fynbo et al., 2001b)) have higher metal column densities and contain more dust than DLAs at similar redshifts. They further suggested that the relatively low reddening (restframe wavelength $`>1200`$ Å) observed towards these GRBs may be due to a grey extinction law, produced by dust destruction (Perna, Lazzati & Fiore, 2003).

Castro et al. (2003) argue that the chromium-to-zinc ratio towards GRB 000926 (catalog ) indicates that the host is depleted in dust relative to local values to a similar degree as QSO DLAs at the same redshift. This conclusion is consistent with what we have suggested above for GRB 020124 (catalog ), namely that the GRB surroundings have a low dust-to-gas ratio, similar to QSO DLAs (Pettini et al., 1997b). An additional argument in favor of this interpretation is that a significant amount of grey extinction would make the afterglow intrinsically even brighter than indicated above (§4.3). For example, assuming a Galactic value of N(H I)/$`A_V`$, log N(H I)= 21.7, destruction of 2/3 of the dust producing the extinction in the rest-frame V band (Perna, Lazzati & Fiore, 2003) and a flat extinction curve would dim the afterglow by 1 mag. Any non-grey component would significantly increase this value.

We conclude that the high value of the ratio between column density and optical extinction first found by Galama & Wijers (2001) from X-ray spectroscopy remains when N(H I) is estimated from optical spectroscopy in GRB afterglows with a DLA host and the extinction is estimated from the observed reddening of the optical/NIR afterglow. The dust destruction interpretation originally proposed by Galama & Wijers (2001) to account for the X-ray result could conceivably also be applied in this case. The observed small reddening would be due to a modification of the extinction law due to preferential destruction of small grains by the prompt UV and X-ray radiation from the GRB, leading to a grey extinction law. However, we have found that the straightforward alternative explanation, namely that GRBs occur in less chemically enriched environments, similar to those of the SMC and QSO DLAs, is fully consistent with the optical observations reported here. Future joint optical and X-ray spectroscopy of $`1.6z2.5`$ GRBs combined with accurate multi-wavelength observations of the afterglow, in particular in the optical/NIR regime, can be used to distinguish between these two possibilities and at the same time fix the metallicity of the burst environment.

We thank Jochen Greiner, Sylvio Klose, and Kristian Pedersen for useful comments. JPUF acknowledges support from the Carlsberg Foundation. JMCC acknowledges the receipt of a FPI doctoral fellowship from Spain’s Ministerio de Ciencia y Tecnología. The authors acknowledge benefits from collaboration within the EU FP5 Research Training Network “Gamma-Ray Bursts: An Enigma and a Tool”. This work was also supported by the Danish Natural Science Research Council (SNF).


format: text

text:
# Protostellar Disk Dynamos and Hydromagnetic Outflows in Primordial Star Formation

## 1 Introduction

Primordial star formation has attracted much attention recently, due to advances in numerical models of cosmic structure formation that can resolve down to solar system scales (Abel, Bryan, & Norman 2002, hereafter ABN), constraints on the reionization epoch from the Wilkinson Microwave Anisotropy Probe observations of the microwave background (Kogut et al. 2003), and the discovery of very metal poor Galactic halo stars (Aoki et al. 2002; Christlieb et al. 2002), whose elemental abundances may be the fingerprints of the earliest supernovae.

The typical mass of the first stars is an important, but unknown, quantity. It may be set by feedback processes local to the forming star. To investigate this possibility it is necessary to model the star formation process by including as much physics as possible, whilst making allowance for a realistic geometry of the accreting gas. Magnetic fields play crucial roles in present-day star formation, but have usually been assumed to be unimportant in the primordial case. It is the goal of this paper, which builds on the detailed, yet non-magnetic, star formation model of Tan & McKee (2003a), to critically examine this assumption.

At typical locations in the standard cosmological model, baryonic over-densities collapse to form the first stars at redshifts $`20`$ (ABN, Bromm, Coppi, & Larson 2002). These first “complex” baryonic systems, consisting of a protostar and accretion disk, live for many dynamical times so that amplification of seed magnetic field by dynamo action may be possible. If the field strength becomes dynamically significant, then an outflow can be generated.

The ubiquity and energy of molecular outflows from present-day protostars (Zuckerman, Kuiper, & Kuiper 1976; Kwan & Scoville 1976; Edwards & Snell 1982; Bally & Lada 1983) was a surprise and forced a radical change in theoretical models of star formation. Outflows are now thought to play a key role redistributing angular momentum in the accretion disk (Blandford & Payne 1982), setting the star formation efficiency (Levreault 1984; Myers et al. 1988; Matzner & McKee 2000), and depositing energy into the star-forming interstellar medium (Norman & Silk 1980). It is often speculated that magnetic fields and outflows may be unimportant for primordial star formation because the strength of seed fields is presumably very weak. Thus they have generally been ignored in studies of the collapsing gas cloud (e.g. Haiman, Thoul, & Loeb 1996; Tegmark et al. 1997; Nakamura & Umemura 1999; ABN; Bromm et al. 2002) and in models of the formation of the protostellar core (Stahler et al. 1986; Omukai & Nishi 1998; Ripamonti et al. 2002; Tan & McKee 2003a). However, as in the local Universe, the presence or absence of dynamically-important magnetic fields and outflows from primordial protostellar disks has important implications for the basic star formation process.

The amplification of seed magnetic field by dynamo action in Population III protostars and its ejection in outflows may also bear on the long-standing question of the origin of galactic and intergalactic field strengths (e.g. Kulsrud et al. 1997, and references therein), particularly if this mode of star formation was widely and evenly spread throughout the Universe. Such a model would be analagous to those which propose that fields are ejected from active galactic nuclei (AGN) (Daly & Loeb 1990; Colgate & Li 1999; Kronberg et al. 2001), though the AGN models would produce more ordered flux on large scales. If field amplification in Population III star formation is very efficient, then much of the intergalactic magnetic field may have been in place before or at the same time as the era of reionization, which has implications for models of field generation in cosmic ionization fronts (Subramanian, Narasimha, & Chitre 1994; Gnedin, Ferrara, & Zweibel 2000).

Starting from conditions representative of the end state of the numerical simulation of ABN, Tan & McKee (2003a, hereafter TM) predicted the mass infall rate, accretion disk structure, and protostellar evolution associated with primordial star formation. Using these results, in this paper we consider a wider range of possible accretion disk models, calculate their structure over the evolution of the protostar and assess their gravitational stability with respect to fragmentation (§2). We discuss the ability of these disks to generate large scale magnetic fields via a dynamo (§3), which is likely to be a key requirement (Lubow et al. 1994; Blackman 2003; Blackman & Tan 2003) for magnetically mediated outflows from thin disks (e.g. Blandford & Payne 1982; Lovelace et al. 1987; Pelletier & Pudritz 1992; Königl & Pudritz 2000). We show that strong, dynamically important fields can be produced even from the weak seed fields expected of primordial protostellar disks. We address the implications of feedback from hydromagnetic winds on the star formation process in §4, and conclude in §5.

## 2 Primordial Protostellar Accretion Disks

The dark matter halos and gas clouds relevant to primordial star formation form from hierarchical mergers of smaller structures, including the accretion of matter along filaments (ABN). Thus in general we expect that the initial conditions for star formation are gas cores with non-zero angular momentum. The collapse should proceed to a disk, which then channels mass to the protostar via viscous torques.

We use TM’s model of the rotating, freely falling accretion envelope that develops inside the sonic point, $`r_{\mathrm{sp}}`$, of the collapsing gas. The rotation is parameterized by

$$f_{\mathrm{Kep}}\frac{v_{\mathrm{rot}}(r_{\mathrm{sp}})}{v_{\mathrm{Kep}}(r_{\mathrm{sp}})}=\frac{v_{\mathrm{rot}}(r_{\mathrm{sp}})}{(GM_{\mathrm{sp}}/r_{\mathrm{sp}})^{1/2}},$$
(1)
which is the ratio of the mass-weighted rotational velocities, $`v_{\mathrm{rot}}`$, relative to Keplerian, $`v_{\mathrm{Kep}}`$, of infalling matter at a radial distance from the protostar corresponding to the mean location of the sonic point. Inside this region, which contains a total mass $`M_{\mathrm{sp}}`$, we assume conservation of angular momentum (Ulrich 1976) so that matter falls freely to the midplane and forms a disk of size

$$r_d=f_{\mathrm{Kep}}^2r_{\mathrm{sp}}3.44\left(\frac{f_{\mathrm{Kep}}}{0.5}\right)^2\left(\frac{1+f_d}{ϵ_d}\right)^{9/7}\left(\frac{m_{}}{M_{}}\right)^{9/7}K^{10/7}\mathrm{AU}.$$
(2)
For the numerical evaluation we have used the mass-radius relation expected at the sonic point of the collapsing cloud (TM), which involves the entropy parameter, $`K^{}`$, of the polytropic equation of state of the cloud, relative to the fiducial case with $`T=300\mathrm{K}`$ at a hydrogen density of $`10^4\mathrm{cm}^3`$ — larger values of $`K^{}`$ correspond to denser gas cores. We have also assumed an efficiency, $`ϵ_d`$, of collapse from the core to the disk, which is unity when feedback processes are unimportant, and related $`M_{\mathrm{sp}}`$ to the stellar mass, $`m_{}`$, at the time when $`M_{\mathrm{sp}}`$ has collapsed to the disk: $`m_{}(1+f_d)=M_{\mathrm{sp}}`$, where $`f_d`$ is the disk mass relative to the star’s and with the assumption that there is little diversion of the accretion flow from $`r_{\mathrm{sp}}`$ to $`r_d`$. From the results of the simulations of ABN we have taken a fiducial value of $`f_{\mathrm{Kep}}=0.5`$, but there is likely to be a significant dispersion depending on the formation history of the halo and cloud. The fiducial value of $`f_d`$ adopted by TM is $`1/3`$, because the disks are expected to be relatively massive since their angular momentum transport must initially be driven by gravitational instabilities.

Material falls onto the disk (and directly to the star) at all radii $`r<r_d`$ (Figure 2 of TM; Cassen & Moosman 1981). The characteristic disk scale is much greater than the stellar radius ($`r_{}100R_{}=0.5\mathrm{AU}`$, TM) so most matter accretes to the star via the disk. The evolution of $`r_{}`$ as the stellar mass increases is shown in Figure 1a, together with the linear growth of $`r_d`$. Although infall onto the disk provides some effective viscosity, this is not important in the very inner regions and so we assume a standard “$`\alpha _{\mathrm{ss}}`$” model for the viscosity (Shakura & Sunyaev 1973). We discuss the choices for the value of $`\alpha _{\mathrm{ss}}`$ below.

The overall rate of collapse and thus accretion rate of the star is given approximately by $`\dot{m}_{}m_{}/t_{\mathrm{ff}}`$, where $`t_{\mathrm{ff}}`$ is the local free-fall timescale evaluated at the density of the region that contains mass $`m_{}`$ at the moment it undergoes dynamical collapse. The accretion rate is thus set by the structure of the pre-stellar gas core, which in turn is set by the cooling properties of trace amounts of molecular hydrogen in the primordial gas. A more careful evaluation of the accretion rate to the star (TM) yields

$$\dot{m}_{}=0.026(1+f_d)^{10/7}ϵ_d^{10/7}K^{15/7}\left(\frac{m_{}}{M_{}}\right)^{3/7}M_{}\mathrm{yr}^1.$$
(3)
These high accretion rates make it likely that the disk will build itself up to a mass significant compared to the stellar mass (we have set $`f_d=1/3`$) and will become susceptible to gravitational instabilities. Two dimensional simulations of clumpy, self-gravitating disks show self-regulation with $`\alpha _{\mathrm{ss}}(\mathrm{\Omega }t_{\mathrm{th}})^1`$ up to a maximum value $`\alpha _{\mathrm{ss}}0.3`$ (Gammie 2001), where $`\mathrm{\Omega }`$ is the orbital angular velocity, $`t_{\mathrm{th}}\mathrm{\Sigma }kT_{\mathrm{c},\mathrm{d}}/(\sigma T_{\mathrm{eff},\mathrm{d}}^4)`$ is the thermal timescale, $`\mathrm{\Sigma }`$ is the surface density, $`T_{\mathrm{c},\mathrm{d}}`$ is the disk’s central (midplane) temperature, and $`T_{\mathrm{eff},\mathrm{d}}`$ the effective photospheric temperature at the disk’s surface. Fragmentation occurs when $`\mathrm{\Omega }t_{\mathrm{th}}3`$: this condition has the best chance of being satisfied in the outermost parts of the disk that are still optically thick (see below). In addition to local gravitational viscosity, global ($`m=1`$ mode) instabilities (Adams et al. 1989; Shu et al. 1990) may be efficient at driving inflow if the disk becomes too massive.

In addition to gravitational instabilities, the magneto-rotational instability (MRI) may become important, yielding viscous stresses that correspond to somewhat smaller values of $`\alpha _{\mathrm{ss}}`$. This process could be particularly important in the inner accretion disk, where gravitational instability is suppressed. Values often quoted for the MRI are $`\alpha _{\mathrm{ss}}0.01`$ (Balbus & Hawley 1998), but larger values have also been measured and neither simulations nor physical understanding have yet determined why it cannot differ by an order of magnitude (or whether the constant $`\alpha _{\mathrm{ss}}`$ formalism is really appropriate for MRI-driven turbulence). We evaluate the ability of the MRI to operate in §3. If the disk is stable with respect to the MRI, then $`\alpha _{\mathrm{ss}}`$ may be set by molecular viscosity and thus much smaller. In this case the accretion rate may effectively go to zero. However, the disk is still being fed at a rate given by equation (3), so that the surface mass density inevitably builds up until self-gravity (and the associated larger values of $`\alpha _{\mathrm{ss}}`$) becomes important.

Given these considerations, we calculate the radial structure of the inner disk at any given point in the protostellar evolution by assuming it is fed smoothly at a rate given by equation (3). We also use the standard theory of steady, thin, viscous accretion disks, with a spatially constant viscosity parameter, $`\alpha _{\mathrm{ss}}`$, and ignore energy injection from the star. The viscosity is assumed to be a function of the total pressure, though it is gas pressure that dominates over radiation pressure in the earlier stages ($`m_{}20M_{}`$), which we think are relevant for field growth (§3). We evaluate cases with $`\alpha _{\mathrm{ss}}=0.01,0.3`$.

The inner scale of the disk is set by $`r_{}`$, which is predicted from the protostellar evolution model of TM. It is somewhat sensitive to the choice of $`\alpha _{\mathrm{ss}}`$ for $`m_{}20M_{}`$ because this sets the temperature of gas in the inner disk and this thermal energy is advected into the star (Fig. 1a). We adopt primordial gas abundances $`X=0.76`$ and $`Y=0.24`$ so that in the atomic phase $`\mu =1.22m_\mathrm{H}`$ and $`n_{\mathrm{He}}=0.079n_\mathrm{H}`$, and use the opacities of zero metallicity gas of Rogers & Iglesias (1996) for $`T6000\mathrm{K}`$ and of Lenzuni, Chernoff, & Salpeter (1991) for $`T6000\mathrm{K}`$ (however, here we are interested primarily in the higher temperature regime). The method assumes that the disk’s central temperature, $`T_{\mathrm{c},\mathrm{d}}`$ is much greater than the surface temperature, $`T_{\mathrm{eff},\mathrm{d}}`$, i.e. that the disk is very optically thick: this condition is not well satisfied during the very earliest stages, or for the outer disk, since the opacity is set by $`\mathrm{H}^{}`$ and declines rapidly for $`T7000\mathrm{K}`$.

We calculate the ionization state of H and He from the Saha equation and include the energy needed for ionization of H and He and the thermal energy in the disk’s energy equation:

$$F=\frac{3Gm_{}\dot{m}_{}}{8\pi r^3}\left[1\left(\frac{r_{}}{r}\right)^{1/2}\right]+\frac{\dot{m}_{}}{4\pi r}\frac{d}{dr}\left(\frac{5}{3}\overline{ϵ}_{\mathrm{th}}+\overline{ϵ}_I\right),$$
(4)
where $`F`$ is the flux emitted from one side of the disk, and $`\overline{ϵ}_{\mathrm{th}}`$ and $`\overline{ϵ}_I`$ are the thermal energy and dissociation/ionization energy per unit mass, averaged over the thickness of the disk. The first term on the right hand side of eq. (4) is the viscous dissipation per unit disk face area. Note that usually $`d\overline{ϵ}_{\mathrm{th}}/dr<0`$ and $`d\overline{ϵ}_I/dr<0`$. The effect of the thermal energy term is small for thin disks such as those we consider, so we drop it from our numerical calculations. However, the dissociation/ionization term, which is often neglected, can be an order of magnitude larger than thermal term and can have a major impact on the disk structure. We therefore retain this term. We follow Frank, King, & Raine (1995), using an approximate solution based on radiative diffusion so that $`F4\sigma T_{\mathrm{c},\mathrm{d}}^4/3\tau `$, where $`T_{\mathrm{c},\mathrm{d}}`$ is the midplane temperature and $`\tau `$ is the optical depth from the midplane to the surface.

Dividing the disk into many discrete radial zones, we solve the full set of disk structure equations (e.g. Frank et al. 1995), starting from the outer regions where the gas is almost completely neutral. The high accretion rates of primordial protostars lead to the ionization energy being quite important, particularly in the earlier stages of evolution before the protostar starts contracting towards the main sequence. For complete ionization of atomic gas, the energy absorption rate is

$$\left(\frac{dE}{dt}\right)_{\mathrm{I},\mathrm{max}}=\dot{m}_{}ϵ_{\mathrm{I},\mathrm{max}}=2.4\times 10^3\left(\frac{\dot{m}_{}}{10^2M_{}\mathrm{yr}^1}\right)L_{},$$
(5)
calculated assuming a dissociation energy of 13.6 eV for H ionization, and 24.6 eV and 54.4 eV for the first and second stages of He ionization.111The results of ABN indicate that the first solar mass or so of material to be accreted would be almost fully molecular, which for a dissociation energy of 4.48 eV for $`\mathrm{H}_2`$ increases the coefficient of equation (5) by 11%. We do not include this minor correction factor in our numerical calculations.

Figure 2 shows three examples of disk structure for $`\alpha _{\mathrm{ss}}=0.01`$, with $`m_{}=1,10,100M_{}`$. At these masses the stellar sizes are taken to be $`r_{}=100,300,4R_{}`$ and the accretion rates (eq. 3) are $`\dot{m}_{}=(17,6.4,2.4)\times 10^3M_{}\mathrm{yr}^1`$, respectively. Figure 3 shows the equivalent models with $`\alpha _{\mathrm{ss}}=0.3`$.

The self-gravity in the disks can be gauged by considering the Toomre stability parameter for a Keplerian disk, $`Qc_s\mathrm{\Omega }/(\pi G\mathrm{\Sigma })`$, which is $`<1`$ in unstable regions. The outer parts of the disk are most susceptible to gravitational instability. Following a line of reasoning similar to that presented by Goodman (2003) concerning self-gravity in quasar accretion disks, consider a disk with $`\alpha _{\mathrm{ss}}=0.01`$ provided by the MRI: the model will be physically self-consistent inside the radius at which $`Q=1`$. At larger radii, we expect that gravitational viscosity will become important, allowing larger values of $`\alpha _{\mathrm{ss}}`$ so that $`Q`$ remains approximately equal to unity. The results of Gammie (2001) show that once $`\alpha _{\mathrm{ss}}0.3`$ then fragmentation occurs. This corresponds to the location in our $`\alpha _{\mathrm{ss}}=0.3`$ models where $`Q=1`$. From Figure 3 we can see that $`Q>1`$ for all the regions of the disks with $`T_{\mathrm{c},\mathrm{d}}7000\mathrm{K}`$, and so fragmentation should not occur in these regions. At larger radii we expect the disks to become optically thin because of the temperature dependence of $`\mathrm{H}^{}`$ opacity (approximately $`T^{14}`$), and have thermal timescales long compared to the orbital period. This appears to occur before $`Q<1`$, so that fragmentation would then not occur anywhere in the disk.

The infall of material onto the disk provides an effective viscosity and source of turbulence that is likely to be quite important in the outer disk, and that this will also help prevent fragmentation. Even if $`Q`$ became $`<1`$, fragmentation is not inevitable: global modes of gravitational instability, such as spiral waves, could develop instead. These modes, if present, are relatively efficient at transporting angular momentum and could then connect the outer infall region with the inner, viscous accretion disk, where first local gravitational viscosity, and then perhaps viscosity from the MRI, bring material all the way to the immediate vicinity of the protostar. The final accretion to the stellar surface may be mediated in part by angular momentum loss in protostellar outflows, the properties of which are discussed in §4.

We conclude that fragmentation of typical primordial protostellar accretion disks does not occur, if the results of Gammie’s (2001) simulations are applicable. However, further numerical studies are desirable to confirm this prediction. We note that in the most advanced simulations to date (ABN) of the global primordial star formation process, which probe conditions in the initial collapse down to scales of order astronomical units (just before formation of a hydrostatic protostellar core and an optically thick accretion disk), there is also no indication of fragmentation. Primordial star formation appears to favor single (or perhaps binary) stars rather that star clusters. The reservoir of gas available to the forming star is several hundreds to thousands of solar masses, but the typical mass of the stars is likely to be set by feedback from the protostar acting on its own accretion. Following TM, here we shall consider the uninterrupted growth of the star at rate given by eq. (3), accreting from a stable accretion disk, and use this as the framework for assessing magnetic field generation and subsequent feedback effects.

## 3 Dynamos and Protostellar Disks

Stahler et al. (1986) and Omukai & Palla (2001, 2003) considered the structure of primordial protostars in the limit of spherical accretion. They found that the protostar is radiatively stable (apart from a very thin outer layer, $`10^7`$ of the stellar mass; Omukai 2003, private comm.) until $`m_{}100M_{}`$, when nuclear H burning creates a convective core. Thus, except at very high masses, little dynamo amplification in the protostar is expected. We therefore study helical disk dynamos, focusing on the production of large scale ordered fields because these are the most susceptible to escape into the corona (Blackman 2003; Blackman & Tan 2003). There they can drive outflows that produce a feedback on the star formation process. A non-helical dynamo may also be present, but this produces disordered random fields that do not drive a strong outflow. For the helical dynamo to operate, turbulence in a $`stratified`$ rotator is the crucial requirement.

### 3.1 Basic Theory

We discuss likely sources of turbulence, the ability of MRI-driven turbulence to produce large scale fields, and then the helical dynamo.

#### 3.1.1 Sources of Turbulence

We expect that the initial driver of turbulence is gravitational instability (Gammie 2001). This may then be followed by MRI-driven turbulence (Balbus & Hawley 1998).

For sufficiently ionized and magnetized disks, the combination of radially decreasing angular velocity in a rotating plasma and seed magnetic field can in principle drive the MRI (Balbus & Hawley 1991, 1998). The vertical stratification then provides the turbulence with the helical property from which a large scale helical dynamo feeds. However, the critical wave number below which the MRI grows, depends on the initial field strength: there is a minimum initial field for the instability to drive turbulence in the disk. Since primordial star-forming environments have only very weak seed fields (discussed later), it is important to evaluate this critical field strength. It is determined by requiring that the MRI growth time is shorter than the diffusion time at the scale of the initial critical vertical growth wavenumber associated with the initial seed field (Balbus & Hawley 1991;1998). This wavenumber is $`k_c^2=3\mathrm{\Omega }/v_A`$. Since the growth rate is of order $`\mathrm{\Omega }`$, the condition that a given dissipation process is slower than the growth rate requires

$$\chi k_c^2/\mathrm{\Omega }=3\chi \mathrm{\Omega }/v_A^2<1,$$
(6)
where $`\chi `$ is a particular microphysical diffusion coefficient to be considered. Note that if $`\chi `$ is an ambipolar diffusion coefficient, then no constraint on the field emerges since then $`\chi v_A^2`$. However, if $`\chi `$ is the Spitzer magnetic diffusivity (Spitzer 1967), then the condition becomes

$$\left(\frac{B}{\mathrm{G}}\right)^2>4.0\times 10^4T_3^{3/2}\left(\frac{m_{}}{M_{}}\right)^{1/2}\left(\frac{r}{1000R_{}}\right)^{3/2}\left(\frac{\mathrm{ln}\mathrm{\Lambda }}{10}\right)^{1/2}\left(\frac{\rho }{10^9\mathrm{g}\mathrm{cm}^3}\right),$$
(7)
where $`\mathrm{\Lambda }`$ is the usual Coulomb logarithm and $`T_3=T/10^3\mathrm{K}`$.

From this estimate we can see that while the MRI process may be able to operate given the expected initial conditions of present-day star formation, in the early universe the fields are likely to be too weak: we estimate in §3.2 that the initial fields that thread the disk have a maximum value of order $`10^{16}\mathrm{G}`$. Yet, even if the MRI does not work initially, gravitational instability in the disk can in principle provide local turbulence that is amenable to an $`\alpha _{\mathrm{ss}}`$ treatment (e.g. Gammie 2001; Wada et al. 2002). As we argued in §2, lacking other sources of viscosity, the accretion disk should establish a state in which the accretion is driven by gravitational instability. The angular momentum transport could also be global, via modes that are not describable with the $`\alpha _{\mathrm{ss}}`$ formalism. The issue of which mechanism operates and at what levels is not fully resolved (Balbus & Papaloizou 1999), requiring further numerical study. Here we assume that local turbulence does arise.

In this context and for the interpretation of simulations, it should be noted that the meaning of $`\alpha _{\mathrm{ss}}`$ is sometimes misunderstood: a turbulent disk incurs large fluctuations so that simulation data must be averaged temporally or spatially for appropriate comparisons to the $`\alpha _{\mathrm{ss}}`$ formalism to be made. The formalism for the helical dynamo that we discuss below, only requires that the stratified helical turbulence yields a viscosity that can be modeled by an $`\alpha _{\mathrm{ss}}`$ parameter in the mean statistical sense.

In a turbulent disk, the small scale field grows faster than the helical large scale field by a factor of order $`\alpha _{\mathrm{ss}}^{1/2}`$ (Blackman 1998, 2003). Thus, well before the large scale field growth saturates, the small scale field will grow to the critical value required for the MRI to operate. For this reason, even if the initial turbulence is driven by gravitational instability, we expect the later phase of field amplification to be mediated by the MRI and its associated $`\alpha _{\mathrm{ss}}`$.

#### 3.1.2 Comment on MRI and Large Scale Fields

Field amplification by the MRI in a non-stratified disk is sometimes said to produce large-scale fields (Balbus & Hawley 1998), but the meaning is subtle and easily misunderstood. In the simplest MRI simulations that use a non-stratified periodic box, the saturated magnetic energy spectrum emerges as Kolmogorov with a $`k^{5/3}`$ shape. No global flux can grow in a periodic box. The toroidal field has the largest coherent scale, but this is only of order the disk height. Most importantly, although the energy of this field is sustained, the sign of the field fluctuates rapidly, on a timescale of order the turbulent correlation time (about equal to a Keplerian rotation time for MRI unstable disks). The reason is that scales of order the box thickness are themselves unstable to the MRI, and thus are the largest turbulent scale. The largest scale fields produced by the non-stratified MRI have the same spatial and temporal properties of the largest scale of the turbulence (Hawley et al. 1995). Such turbulent scale fields can be shredded before buoyantly rising to the coronae (Blackman 2003; Blackman & Tan 2003), and are therefore not ideally suited for driving outflows.

However when stratification is present, the MRI can produce helical turbulence. “Helical” describes the property that rising turbulent eddies in the northern (southern) hemisphere expand and twist clockwise (counterclockwise) to conserve angular momentum, while falling turbulent eddies twist counterclockwise (clockwise). Both rising and falling plasma maintain the same sign of pseudoscalars such as $`𝐯\times 𝐯`$ and $`𝐛\times 𝐛`$, (where $`𝐯`$ and $`𝐛`$ are the turbulent velocity and turbulent magnetic field in Alfvén units, respectively), with each having opposite signs in opposite hemispheres. Such finite pseudoscalars enable amplification of global scale fields by helical dynamo action when the rising plasma is threaded by a magnetic field (e.g. Moffatt 1978; Parker 1979; Brandenburg 2001; Blackman & Field 2002). In a real non-periodic system the helical dynamo can maintain a global large scale flux over many rotation periods, as seen in non-periodic box simulations of disks (Brandenburg et al. 1995; Brandenburg & Donner 1997). Note that Stone et al. (1996) claimed that no helicity was observed in their stratified periodic box disk simulations, when averaged over hemispheres (Stone 2003, private comm.). Further studies are required to reconcile these results.

#### 3.1.3 The Helical Dynamo

Recently, Blackman & Field (2002) developed a dynamical nonlinear helical dynamo theory that accounts for the nonlinear backreaction of the magnetic field on large scale field growth and agrees with recent helical dynamo simulations (Brandenburg 2001). The driver of exponential large scale field growth, $`\alpha _d`$, turns out to be proportional to the residual helicity $`𝐯\times 𝐯𝐛\times 𝐛`$ (Pouquet, Frisch, & Leorat 1976), rather than just the kinetic helicity $`𝐯\times 𝐯`$ used in textbooks for similar dynamos. Large scale field growth occurs in two phases: an early, fast kinematic phase, and a late, slow resistive phase. Generation of the large scale field is associated with large scale magnetic helicity, whilst helicity conservation means that the small scale helicity grows with the opposite sign. As this happens, $`𝐛\times 𝐛`$ grows and offsets the $`𝐯\times 𝐯`$, shutting down the dynamo. These principles seem to survive in dynamos with shear (Blackman & Brandenburg 2002). Although the above results are for simplified dynamos, we shall use the saturation principles in what follows, by incorporating them into an $`\alpha \mathrm{\Omega }`$ accretion disk. Other approaches to mean field accretion disk dynamos have been developed by, for example, Pudritz (1981), Rüdiger, Elstner, & Stepinski (1995), Stepinski (1995), Reyes-Ruiz & Stepinski (1997; 1999) and Campbell & Caunt (1999).

To estimate the large scale field, $`\overline{𝐁}`$, produced by a helical dynamo, we need the growth rate. We use a mean field theory for large scale field growth in a turbulent disk. This approach presumes that the mean quantities vary on time scales long compared to that of the fluctuating turbulent quantities. Therefore, the derived mean field growth rate, $`\gamma `$, should be slower than the eddy turnover rate, i.e. $`\gamma <\mathrm{\Omega }`$, where the latter is also the growth rate of the random field (e.g. Balbus & Hawley 1998). We shall find that this condition is indeed satisfied.

We assume $`\overline{𝐁}=\overline{𝐁}_0\mathrm{Exp}[\gamma t+ik_zz]`$, where $`k_z`$ is the vertical wavenumber. The basic equations for the radial and toroidal field in an “$`\alpha _d\mathrm{\Omega }`$” dynamo for a Kelperian disk become (e.g. Ruzmaikin et al. 1988)

$$\gamma \overline{B}_r=ik_z(\alpha _d\overline{B}_\varphi )\beta k_z^2\overline{B}_r$$
(8)
and

$$\gamma \overline{B}_\varphi =(3/2)\mathrm{\Omega }\overline{B}_r\beta k_z^2\overline{B}_\varphi ,$$
(9)
where we ignore spatial dependencies of the pseudoscalar $`\alpha _d`$ and turbulent transport scalar $`\beta `$, and assume averages are taken within a given hemisphere. We assume $`\alpha _d=q\alpha _{d0}`$ where the time-dependent multiplicative factor $`q1`$ represents the backreaction from $`𝐛\times 𝐛`$, with $`q=1`$ in the kinematic regime.

There are two possible dominant contributions to $`\alpha _{d0}`$: a Coriolis force effect and a magnetic buoyancy effect. First consider the Coriolis force acting on the largest turbulent eddies, of typical scale $`l`$. A rising and expanding (falling and contracting) eddy rotates in the opposite (same) sense as that of the underlying disk rotation during its turbulent correlation time, $`t_{\mathrm{ed}}=l/v`$. Therefore, while conserving angular momentum about the disk axis, such eddies gain a specific angular momentum about their own symmetry axis such that $`𝐯\times 𝐯`$ is non-vanishing and of opposite sign in each hemisphere. Using the equation for $`\times 𝐯`$ evolution and the mass continuity equation, and assuming that the density varies on scale $`h`$, this leads to the estimate (Ruzmaikin et al. 1988)

$$\alpha _{d0}=t_{\mathrm{ed}}𝐯\times 𝐯t_{\mathrm{ed}}^2𝐯(\times (\mathrm{\Omega }\times 𝐯))t_{\mathrm{ed}}^2𝐯𝛀𝐯\frac{\mathrm{\Omega }l^2}{h}\times \mathrm{Min}[R_b,1],$$
(10)
where the last relation follows from mass continuity, and $`R_b\frac{v}{l\mathrm{\Omega }}`$ is the Rossby number.

For thin disks in the Shakura-Sunyaev formalism,

$$l\alpha _{\mathrm{ss}}^{1/2}h/R_b^{1/2}$$
(11)
and

$$v\alpha _{\mathrm{ss}}^{1/2}R_b^{1/2}c_s,$$
(12)
which follow since $`\nu vlv^2/(R_b\mathrm{\Omega })=\alpha _{\mathrm{ss}}c_sh`$ and $`\mathrm{\Omega }hc_s`$. Using (11) in (10), we have

$$\alpha _d=q\alpha _{d0}q\alpha _{\mathrm{ss}}\mathrm{\Omega }h.$$
(13)
It also follows that the maximum fractional kinetic helicity satisfies

$$f_h=\mathrm{Min}[1,\alpha _{d0}/v\alpha _{\mathrm{ss}}^{1/2}/R_b^{1/2}].$$
(14)
It is important to note that while equation (10) is a standard estimate for $`\alpha _{d0}`$, Brandenburg & Donner (1997) found that for large scale field generation in local shearing box simulations, the contribution seemed to have the opposite sign of that required to produce the measured mean field. It was later suggested (Brandenburg 1998) that a magnetic buoyancy term of the form $`t_{\mathrm{ed}}b_xb_y(g_z/p)`$ where $`b_x,b_y`$ are the fluctuating components of the magnetic field, $`g_z`$ is the vertical gravity, and $`p`$ is the pressure, actually has the correct sign. Since this term is $`(1/\mathrm{\Omega })(v_A^2/c_s^2)\alpha _{\mathrm{ss}}(v_\varphi ^2/r)(h/r)\alpha _{\mathrm{ss}}c_s\mathrm{\Omega }l^2/h`$, which has the same order of magnitude and scalings as equation (10), it does not significantly affect our calculations based on the use of equation (13), unless the two cancel to high order. The cycle period analyzed by Brandenburg & Donner (1997) suggests some cancellation. However, they studied the late time evolution of the dynamo, after the kinematic regime saturates.

With our assumed form for $`\overline{𝐁}`$, we solve equations (8) and (9), use equation (13), and set $`k_z1/h`$ and $`\beta \nu =\alpha _{\mathrm{ss}}c_sh`$, to derive the growth rate

$$\gamma =(3\alpha _dk_z\mathrm{\Omega }/4)^{1/2}\beta k_z^2=\alpha _{\mathrm{ss}}^{1/2}c_s(\sqrt{3}q/2\alpha _{\mathrm{ss}}^{1/2})/h.$$
(15)
Note that with $`q=1`$, for growth this requires a dynamo number $`D|(3/2)\mathrm{\Omega }\alpha _0h^3/\beta ^2|>2=D_c`$, the critical dynamo number. Critical dynamo numbers for disks are sometimes derived to be $`>8`$ (e.g. Ruzmaikin et al. 1988) but this value depends on the assumed profile of $`\alpha (z)`$. We have crudely assumed a nearly constant profile with $`z>0`$, and thus our critical dynamo number is low.

At a given stage in the evolution when the stellar mass is $`m_{}`$, we estimate $`\overline{B}`$ as the smaller of (i) the saturation value at the end of the kinematic regime and (ii) the value reached given the growth time that the disk dynamo has been able to operate,

$$t_{\mathrm{grow}}=t(m_{})t_{\mathrm{init}},$$
(16)
where $`t_{\mathrm{init}}`$ is the time at which the protostellar disk becomes susceptible to turbulence. We define $`t=0`$ at the time when the protostellar core has just started to form, then from equation (3) we have

$$t(m_{})=27(1+f_d)^{10/7}ϵ_d^{10/7}K^{15/7}(m_{}/M_{})^{10/7}\mathrm{yr}.$$
(17)
Several factors potentially determine $`t_{\mathrm{init}}`$. Firstly, we require a critical level of ionization so that the field’s growth time due to the dynamo is substantially shorter than the diffusion time. Brandenburg et al. (1995) found a critical fractional ionization of disk material of $`10^{14}`$ for growth of the MRI, corresponding to a neutral-ion collision frequency $`10`$ times the orbital frequency. The MRI has a growth time $`\alpha _{\mathrm{ss}}`$ times that of the large scale helical field and so our critical ionization fraction is $`\alpha _{\mathrm{ss}}^110^{14}10^{13}`$. This condition is well-satisfied for the inner regions of primordial protostellar disks, once the star has a mass that is even just a small fraction of a solar mass. Thus $`t_{\mathrm{init}}`$ (and the corresponding $`m_{,\mathrm{init}}`$) is set by the time for disk formation in the collapsing flow, i.e. when the centrifugal radius of the flow becomes significantly larger than the protostellar radius. For the fiducial case of $`K^{}=1`$ and $`f_{\mathrm{Kep}}=0.5`$, this occurs very early on, when the star has mass $`m_{,\mathrm{init}}0.4M_{}`$.

One might think that the time available for field growth at a given radius is the accretion timescale at that radius, rather than the disk age. However, since the mean field formalism is statistical, amplification occurs as long as the system remains in quasi-statistical equilibrium (i.e. as long as helically turbulent plasma is steadily supplied to that location). A related point is that as long as $`\beta >(L^2/r^2)\nu `$, where $`Lh`$ is the vertical variation scale of the large scale field, a poloidal mean field is likely not frozen into the plasma on an accretion timescale and would diffuse away if the $`\alpha _d`$ helicity effect is ignored. This then motivates the need for the $`\alpha _d`$ effect and large scale dynamos in the first place (Lubow et al. 1994; Blackman 2003; Blackman & Tan 2003). For completeness however, we will also provide the estimated field when the growth time is limited to the local accretion time.

After the kinematic regime, subsequent field growth becomes resistively limited as the kinetic and current helicities nearly balance. The magnitude of the helical component of the mean field when this occurs is given by (Blackman & Field 2002; Blackman 2003)

$$\overline{B}_H(4\pi \rho f_h)^{1/2}v(l/h)^{1/2}(4\pi \rho )^{1/2}\alpha _{\mathrm{ss}}c_s,$$
(18)
where the latter similarity follows from equations (11), (12), and (14). We take this to be an estimate of the poloidal field, i.e. we assume $`\overline{B}_r\overline{B}_H`$. We assume that the strength of the toroidal field at the surface $`\overline{B}_\varphi `$, is comparable to its value within the disk. This toroidal field is larger than $`\overline{B}_r`$ because it is linearly amplified above the value of $`\overline{B}_r`$ by shear as it rises into the coronae during a buoyancy or diffusion timescale. These timescales are $`t_{\mathrm{buo}}h(4\pi \rho )^{1/2}/\overline{B}_\varphi `$ and $`t_{\mathrm{dif}}=h/\alpha _{\mathrm{ss}}c_s`$, respectively. Linear growth of $`\overline{B}_\varphi `$ above the level predicted by equation (18) in a time $`t_{\mathrm{buo}}`$ gives $`\overline{B}_\varphi =\overline{B}_r\mathrm{\Omega }t_{\mathrm{buo}}=\alpha _{\mathrm{ss}}^{1/2}c_s(4\pi \rho )^{1/2}`$, which from equations (18) and (12) implies $`\overline{B}_r/\overline{B}_\varphi =\alpha _{\mathrm{ss}}^{1/2}`$. Similarly, if the growth of $`\overline{B}_\varphi `$ is limited by turbulent diffusion, then $`\overline{B}_\varphi =\overline{B}_r\mathrm{\Omega }t_{\mathrm{dif}}`$, so that $`\overline{B}_r/\overline{B}_\varphi \alpha _{\mathrm{ss}}`$.

To summarize: for an initial seed field $`\overline{𝐁}_0`$ we expect the final field strengths to be

$$\overline{B}_r\mathrm{Min}\{\overline{B}_0\mathrm{Exp}[\gamma _0t_{\mathrm{grow}}],\alpha _{\mathrm{ss}}c_s(4\pi \rho )^{1/2}\},$$
(19)
where $`\gamma _0\alpha _{\mathrm{ss}}^{1/2}c_s/h=\alpha _{\mathrm{ss}}^{1/2}\mathrm{\Omega }`$ is the kinematic $`\gamma `$, obtained from equation (15) with $`q1`$, and

$$\overline{B}_\varphi =\alpha _{\mathrm{ss}}^1\overline{B}_r$$
(20)
if $`\overline{B}_\varphi /(4\pi \rho )^{1/2}<\alpha _{\mathrm{ss}}c_s`$, (i.e. when $`t_{\mathrm{buo}}>t_{\mathrm{dif}}`$) or

$$\overline{B}_\varphi =\alpha _{\mathrm{ss}}^{1/2}\overline{B}_r$$
(21)
if $`\overline{B}_\varphi /(4\pi \rho )^{1/2}>\alpha _{\mathrm{ss}}c_s`$, (i.e. when $`t_{\mathrm{buo}}<t_{\mathrm{dif}}`$). The latter case applies in the saturated state.

### 3.2 Expected Values for Star-Forming Disks

To determine the expected $`\overline{B}`$ in star-forming disks, we must estimate the seed field, assess flux freezing for the initial collapse of the cloud to form the disk, and estimate the growth timescale to determine if the field has saturated. The seed field for the primordial star-forming environment can be estimated from the Biermann battery mechanism for non-barotropic flows. On scales much larger than that of an indvidual star, numerical simulations indicate that by $`z=18`$, seed fields of order $`\overline{B}_010^{26}`$ Gauss can be generated (Kulsrud et al. 1997). We assume such a field is generated inside the pre-stellar core region in the simulation of ABN, which has a size $`100\mathrm{pc}`$ and a density $`10^{26}\mathrm{g}\mathrm{cm}^3`$. The infall speeds within this region are of order km/s, vastly greater than the drift speeds of the ions through the neutrals (the ionization fractions are much greater than the equilibrium values for the gas, which has a minimum temperature of $`200\mathrm{K}`$). Thus field is frozen into the contracting gas: going from parsec to AU scales the field strength increases by factors of $`10^{10}`$, so the seed field in the disk may then be $`10^{16}\mathrm{G}`$.

It may actually be that much larger seed fields could arise in the protostar from the Biermann battery term if the pressure and density gradients are misaligned inside the protostar as well. The protostellar Biermann battery (e.g. Kulsrud et al. 1997) value reached after a rotation time (at which point other growth terms become important) can be estimated to be $`\overline{B}_0(m_Hc/e)\frac{P\times \rho }{\rho ^2\mathrm{\Omega }}2\times 10^{11}(T/10^5\mathrm{K})^{1/2}(r/10^{13}\mathrm{cm})^1(h/r)`$ Gauss, where we have used $`\mathrm{\Omega }=c_s/h`$ and have conservatively taken $`1/r`$.

In the disk, flux freezing and field advection are not guaranteed since, as discussed above, turbulent diffusion of the mean field may be relatively efficient and the accretion timescale may be relatively long. Thus, in addition to a fiducial initial seed field strength of $`\overline{B}_0=10^{16}\mathrm{G}`$, we also consider a conservative value of $`\overline{B}_0=10^{26}\mathrm{G}`$. Since the disk-dynamo amplifies the initial field exponentially, the uncertainties in the strength of the seed field do not significantly affect the time to reach saturation.

Using $`\gamma _0\alpha _{\mathrm{ss}}^{1/2}c_s/h=\alpha _{\mathrm{ss}}^{1/2}\mathrm{\Omega }`$ from (15) with $`q=1`$, and approximating $`t_{\mathrm{grow}}=t(m_{})=0.7m_{}/\dot{m}_{}`$, the first term on the right of equation (19) becomes

$$\overline{B}_r8.9\times 10^8\mathrm{G}\left(\frac{\overline{B}_0}{10^{26}\mathrm{G}}\right)\mathrm{Exp}\left[\left(\frac{\alpha _{\mathrm{ss}}}{0.01}\right)^{1/2}\left(\frac{m_{}}{M_{}}\right)^{3/2}\left(\frac{r}{10^{13}\mathrm{cm}}\right)^{3/2}\left(\frac{\dot{m}_{}}{0.01M_{}\mathrm{yr}^1}\right)^1\right],$$
(22)
while the second term gives

$$\overline{B}_r49\mathrm{G}(\alpha _{\mathrm{ss}}/0.01)^{1/2}(m_{}/M_{})^{1/4}(\dot{m}_{}/0.01M_{}\mathrm{yr})^{1/2}(r/10^{13}\mathrm{cm})^{3/4}(h/10^{12}\mathrm{cm})^{1/2},$$
(23)
where mass continuity at $`r=2r_{}`$, $`\rho =\dot{m}_{}/(4\pi h^2\alpha _{\mathrm{ss}}c_s)`$, was used, and we choose parameters commensurate with the results of §2. Although primordial protostars are thought to have accretion rates somewhat higher than the $`10^2M_{}\mathrm{yr}^1`$ adopted in equation (22), it is clear from these equations that we expect the field to saturate when the protostar has a relatively small mass. Note that the seed fields for present-day Galactic protostellar disks would be many orders of magnitude higher than the primordial case: the initial Galactic interstellar material possesses a field of $`10^6`$ G. Thus for a wide range of parameters equation (23) is the relevant estimate for the dynamo-amplified large scale field strengths in present-day, as well as primordial, protostellar disks. For the primordial case, we evaluate the saturation point more accurately in the $`m_{}`$ versus $`\dot{m}_{}`$ plane, and show the results in Figure 4.

In the fiducial case the field in the inner disk saturates by about $`1M_{}`$, and this does not depend very sensitively on $`\alpha _{\mathrm{ss}}`$. In lower angular momentum cores the disk emerges later, at which point the dynamo rapidly causes the field to saturate. Although we expect that mean field dynamo can grow field at a given location for longer than the local accretion time, we also evaluate the cases where $`t_{\mathrm{grow}}=t_{\mathrm{acc}}2\pi n_{\mathrm{orb}}/\mathrm{\Omega }r^2/(\alpha _{\mathrm{ss}}c_sh)`$. The resulting field strength must always be smaller than that from equation (22) for disks for which $`\alpha _{\mathrm{ss}}`$ models are self-consistent. The saturation point is delayed until masses $`m_{}1020M_{}`$ in the fiducial case.

## 4 Implications of Large Scale Fields for Primordial Star Formation

### 4.1 The Power of Magnetically Mediated Outflows

The primary effect of the magnetic field that we consider is the generation of a magnetically mediated jet or wind. Integrating the vertical Poynting flux through the disk in each vertical half gives an estimate of the magnetic luminosity available for driving each outflow:

$$L_{\mathrm{mag},1/2}\overline{B}_{\varphi ,s}\overline{B}_{z,s}\mathrm{\Omega }r^3,$$
(24)
where $`\overline{B}_{\varphi ,s}`$ and $`\overline{B}_{z,s}`$ are the toroidal and vertical fields just inside the disk surface (defined at one scaleheight above the midplane). The field strengths we estimated in the previous section represent fields inside the disk, so the use of equation (24) really requires a solution of the global boundary value problem (e.g. Rüdiger, Elstner & Stepinski 1995; Reyes-Ruiz & Stepinski 1997;1999). For present purposes, we make crude estimates for the surface fields. We take $`\overline{B}_{\varphi ,s}\overline{B}_{\varphi ,\mathrm{ave}}=\alpha _{\mathrm{ss}}^{1/2}\overline{B}_{r,\mathrm{ave}}`$ as estimated in the previous section, where $`\overline{B}_{\varphi ,\mathrm{ave}}(3/4)\overline{B}_\varphi `$ and $`\overline{B}_{r,\mathrm{ave}}(3/4)\overline{B}_r`$ are the respective average values over a scaleheight, using flux freezing in which the density falls a factor of $`1/e`$ between the midplane and surface. Note that for a constant $`\alpha _{\mathrm{ss}}`$ disk, $`v_A/c_s`$ is a constant (Blackman 1998). Then with $`c_sr^{1/2}`$ and $`\rho r^{3/2}`$, we have $`B_rr^{5/4}`$, so that $`L_{\mathrm{mag}}`$ is dominated by the contribution generated at the inner disk radius.

Now we use the above relations to estimate $`\overline{B}_{z,s}`$ by appealing to $`\overline{𝐁}=0`$. We integrate this over an annular wedge of the disk extending from the midplane to the surface, with azimuthal extent $`\mathrm{\Delta }\varphi `$, and with radial extent $`\delta r`$, where $`\delta rr`$ is the characteristic radial scale over which the vertical and radial fields change significantly. The toroidal flux does not contribute to this integral, and by the assumption that the vertical field vanishes at the midplane, no net flux penetrates the equatorial plane. Then integrating $`\overline{𝐁}=0`$ over the described volume gives

$$\mathrm{\Delta }\varphi \left[B_{r,\mathrm{ave},i}h_ir_ir_oB_{r,\mathrm{ave},i}\left(\frac{r_i}{r_o}\right)^{5/4}\left(\frac{h_i}{r_i}\right)r_o\right]=\mathrm{\Delta }\varphi \frac{4}{3}\left[B_{z,i}r_i^{5/4}(r_o^{3/4}r_i^{3/4})\right],$$
(25)
where the subscripts $`i`$ and $`o`$ refer to inner and outer radii, and we assume $`h/r=h_i/r_i`$, that is, the ratio is independent of $`r`$, which is a reasonable approximation (see Figs. 2 and 3). In the limit $`r_or_i`$ and using $`B_{r,\mathrm{ave}}=(3/4)B_r`$ from above, we obtain

$$B_{z,i}/B_{r,i}=(9/16)(h/r).$$
(26)
From equation (21), $`\overline{B}_{\varphi ,s}=\alpha _{\mathrm{ss}}^{1/2}B_{r,\mathrm{ave}}=(3/4)\alpha _{\mathrm{ss}}^{1/2}B_{r,i}`$ so that equation (24) for the outflow power becomes

$$L_{\mathrm{mag},1/2}(27/64)(h/r)\alpha _{\mathrm{ss}}^{1/2}\overline{B}_r^2\mathrm{\Omega }r^3.$$
(27)
The procedure for using equation (27) to assess the value of $`L_{\mathrm{mag}}`$ for star-forming disks, is to pick a particular star formation model (which means choosing $`\alpha _{\mathrm{ss}}`$, $`m_{}`$ and $`\dot{m}_{}`$) and inferring from its solution whether the field has saturated or not (from eqs. 22 and 23). The appropriate field strength is then used in equation (27) to estimate the wind luminosity.

The case of a saturated field (eq. 23) is the most interesting:

$`L_{\mathrm{mag},1/2}`$ $``$ $`(27/64)\alpha _{\mathrm{ss}}^{1/2}{\displaystyle \frac{Gm_{}\dot{m}_{}}{r}}=0.84\alpha _{\mathrm{ss}}^{1/2}L_{\mathrm{max},1/2},`$ (28)
$`=`$ $`92\left({\displaystyle \frac{\alpha _{\mathrm{ss}}}{0.01}}\right)^{1/2}{\displaystyle \frac{m_{}}{M_{}}}{\displaystyle \frac{\dot{m}_{}}{0.01M_{}\mathrm{yr}^1}}\left({\displaystyle \frac{r}{10^{13}\mathrm{cm}}}\right)^1L_{},`$ (29)
where $`L_{\mathrm{max},1/2}=Gm_{}\dot{m}_{}/2r`$ is the upper limit on the magnetic luminosity through each disk plane available from accretion. Note that $`L_{\mathrm{mag},1/2}<L_{\mathrm{max},1/2}`$ as expected, for all reasonable values of $`\alpha _{\mathrm{ss}}`$ (formally $`\alpha _{\mathrm{ss}}<1`$). The evolution of $`L_{\mathrm{mag}}`$ with $`m_{}`$ is shown in Figure 5a, together with the radiative luminosity from the fiducial protostellar evolution model of TM. The wind luminosity is a fixed fraction of the accretion luminosity from the inner disk.

Fields produced by a helical dynamo are ordered such that their smallest gradient scale exceeds that which would otherwise result in their turbulent destruction before rising to the corona. But there is an additional process in the generation of global large scale fields in the corona that we have not discussed above: the fields supplied by the disk dynamo must incur further evolution to even larger scales once in the corona. Their scale in the disk is a few times $`h`$, but the needed large scale in the corona is $`r`$. While this further opening or “dynamical magnetic relaxation” process requires further study with direct application to the disk-corona system, there is ample reason to believe that it happens sufficiently fast when the supplied fields from the dynamo are helical (as they would be from a helical disk dynamo, and further supplemented by foot point shear): It has long been known that helical fields relax to the largest scale available to them in magnetically dominated environments (Taylor 1986) and it has been further shown recently that when helical magnetic fields are injected into such a magnetically dominated region, this relaxation occurs on an Alvfén crossing time (Blackman & Field 2003). The speed is independent of the value of the magnetic diffusivity, as long as the latter is not identically zero. For this reason the dynamical relaxation in the corona would naturally allow further opening of disk-supplied fields to global scales, not unlike the “magnetic carpet” evolution in the solar corona (Schrijver & Zwaan 2000).

Given the magnetic luminosity, we estimate the mass-loading and momentum of the outflow. We assume that at large distances from the star, the energy content of the flow is dominated by the kinetic energy of the matter, so that $`L_{\mathrm{mag}}=\dot{m}_wv_w^2/2`$, where $`\dot{m}_w`$ is the mass outflow rate and $`v_w`$ is the terminal velocity. Theoretical models of magneto-centrifugally driven outflows suggest that $`v_w`$ is about equal to the escape velocity from the protostar, $`v_{\mathrm{esc},}`$ (e.g. Shu et al. 2000; Königl & Pudritz 2000). We write $`v_w=f_vv_{\mathrm{esc},}`$, where $`f_v`$ is a factor of order unity. Then the mass flux of the outflow is

$$\dot{m}_w=\frac{27}{32}\frac{\alpha _{\mathrm{ss}}^{1/2}}{f_v^2}\dot{m}_{}0.084\dot{m}_{}.$$
(30)
We note that $`\alpha _{\mathrm{ss}}`$ is uncertain, as are the numerical prefactors of equations (29) and (30). For a massive protostar that is approaching the Eddington limit, the mass-loading of the outflow may be significantly enhanced.

### 4.2 Feedback Effects from Outflows

In present-day star formation, protostellar outflows contribute the following feedback effects: (1) they prevent the accretion of gas in the polar regions of the star-forming core that would otherwise have joined the star and they divert a fraction of the accretion flow in the disk (e.g. Matzner & McKee 2000, hereafter MM); (2) they inject mechanical energy into giant molecular clouds, which helps to keep the clouds turbulent (Norman & Silk 1980); (3) through the presence of dense, outflowing gas near the star, they inhibit the reach of other feedback effects, such as ionizing radiation (Tan & McKee 2003b). We discuss the significance of these processes in the context of primordial star formation.

#### 4.2.1 Star Formation Efficiency

First, consider the expulsion of gas from the polar regions of the core. By balancing magnetic tension and magnetic pressure gradients, Matzner & McKee (1999) (see also Shu et al. 1995; Ostriker 1997) showed that far from the star the angular distribution of the momentum in a radial hydromagnetic wind can be approximated by

$$\frac{dp_w}{d\mathrm{\Omega }}=\frac{p_w}{4\pi }\frac{1}{\mathrm{ln}(2/\theta _0)(1+\theta _0^2\mu ^2)},$$
(31)
where $`\mu =\mathrm{cos}\theta `$ and $`\theta _0`$ is a small angle, which is estimated to be $`0.01`$ for winds from low-mass stars (Matzner & McKee 1999), but may be somewhat larger for winds from massive stars. We adopt this value in the primordial case (our results are not particularly sensitive to this choice). Equation (31) is also a reasonable description for the angular distribution of the mass flux in the wind, assuming the terminal velocity (far from the star, but before deceleration by the ambient gas) is independent of $`\theta `$.

To evaluate the efficiency of star formation from a core (i.e. the fraction of mass that forms a star) $`ϵ_{\mathrm{core}}`$, we follow MM and find the angle, $`\theta _{\mathrm{esc}}\mathrm{cos}^1\mu _{\mathrm{esc}}`$, where the wind sweeps up core material to the surface escape speed of the core, $`v_{\mathrm{esc},\mathrm{c}}`$. The analysis assumes thin, radiative shocks, purely radial motion and monopole gravity. We consider cores with angular mass distributions of the form $`dM/d\mathrm{\Omega }=(1/4\pi )Q(\mu )M`$ with $`Q(\mu )=(1\mu ^2)^n/_0^1(1\mu ^2)^n𝑑\mu `$, and evaluate models with $`n`$, ranging from $`0`$, isotropic, to $`4`$, which describes a flattened distribution that mimics the effect of some rotational support. In the notation of MM, the escape condition is given by $`(1+\theta _0^2\mu _{\mathrm{esc}}^2)Q(\mu _{\mathrm{esc}})=m_{}/(XM)`$, where $`X=0.132c_g[\mathrm{ln}(2/\theta _0)/\mathrm{ln200}]v_{\mathrm{esc},\mathrm{c},5}[(p_w/m_{})/40\mathrm{k}\mathrm{m}\mathrm{s}^1]^1`$, and $`c_g`$ is a factor of order unity that accounts for the effects of gravity on shock propagation. Using the results of MM, we estimate $`c_g4.6`$, for steady winds that decouple from the swept-up shell at the core edge in a core with $`k_\rho =20/9`$ (which is the fiducial value of TM). These cores have $`v_{\mathrm{esc},\mathrm{c}}=3.22(M/1000M_{})^{1/7}K^{5/7}\mathrm{km}/\mathrm{s}`$. We ignore the influence of material beyond the core “edge”, which is partly balanced by our neglect of the wind’s influence on the shell beyond this point. The efficiency is given by

$$ϵ_{\mathrm{core}}=\frac{1}{1+(\dot{m}_w/\dot{m}_{})(1\varphi _w)}_0^{\mu _{\mathrm{esc}}}Qd\mu ,$$
(32)
where $`\varphi _w_0^{\mu _{\mathrm{esc}}}Pd\mu `$ and $`P(\mu )`$ is the dimensionless force distribution in angle. Note that unlike $`ϵ_d`$, which was used in equation (3), this definition of efficiency includes the effect of diversion of material from the disk into the outflow. However, since $`f_w\dot{m}_w/\dot{m}_{}1`$, $`ϵ_{\mathrm{core}}ϵ_d`$.

In summary, we are assuming that the properties of the protostellar outflow, and how they evolve with $`m_{}`$, do not depend on $`ϵ_{\mathrm{core}}`$. In reality, as $`ϵ_{\mathrm{core}}`$ and $`\dot{m}_{}`$ become smaller, then the protostellar evolution changes, e.g. the stellar growth time becomes longer relative to the Kelvin-Helmholz time, so that the star would tend to be smaller. However, since $`p_w`$ is integrated over the stellar evolution, it should not be so sensitive to these changes.

We show the evolution of $`p_w/m_{}`$ with $`m_{}`$ for our fiducial model in Figure 5b. Initially $`r_{}`$ is approximately constant and the specific momentum of the wind grows only because of the increasing mass and escape velocity from the star. Then above $`m_{}20M_{}`$, the star starts shrinking towards the main sequence and there is a jump in the momentum input. This marks the start of the transition to lower efficiency star formation from the core (Fig. 5c). This transition may have important implications for the final mass of the star, and will be examined in more detail in a future paper. Note that this efficiency is not particularly sensitive to the parameter $`\theta _0`$ when the critical angle for ejection of gas from the core is $`\theta _0`$. That is, mass limits to star formation resulting from this process do not depend on $`\theta _0`$ so long as $`\theta _01`$.

#### 4.2.2 Feedback on the Larger-Scale Cloud and Break-Out from the Dark Matter Halo

We have seen that the protostellar outflow is eventually able to break free of the core. With the distribution of momentum given by equation (31), this is particularly true near the rotation axis, since about half of the initial outflow mass is concentrated inside $`\theta _0`$. A $`100M_{}`$ star releases a wind with total momentum of about $`10^4M_{}\mathrm{km}/\mathrm{s}`$, with about a quarter of this inside $`\theta _0`$ of each jet, for which the injected mass is about $`2M_{}`$. Each central jet region can then sweep up $`500M_{}`$ to about $`5\mathrm{km}/\mathrm{s}`$, which is the approximate sound speed and escape speed of the gas in the vicinity of the dark matter mini-halo. The fiducial core of TM has a total virialized baryonic mass of $`10^3M_{}`$ (based on the simulation of Abel et al. 2002), with a few $`\times 10^4M_{}`$ of gas still infalling within a couple of hundred parsecs. For reasonable values of $`\theta _00.1`$, and even with allowance for turbulent entrainment of matter at the sides of the jets and lateral spreading, the outflows should easily be able to eject some gas from the dark matter halo. Thus we expect the extent of mechanical feedback to reach beyond the immediate environment of the star-forming core.

The physical implications of this feedback are difficult to assess, since radiative feedback from the massive protostar is also acting on the cloud during this time. If outflows do escape from the halo, then they will start to pollute the intergalactic medium with their magnetic fields. This is only one possible mechanism for producing such fields. Others include field generation in shocks caused by supernovae or ionization fronts (Subramanian et al. 1994; Gnedin et al. 2000), and these mechanisms should have a longer reach than the protostellar outflows. Outflows from AGN are also likely to be important, although these may be more sparsely distributed than Population III stars.

#### 4.2.3 Confinement of Radiative Feedback by Protostellar Outflows

The presence of dense, outflowing gas near the protostar can suppress the influence of radiative feedback. For example, the propagation of ionizing photons should be impeded. The characteristic density at the base of the flow region near the star is

$$n_{\mathrm{H},\mathrm{base}}=4.0\times 10^{11}f_v^1\left(\frac{f_r}{2}\right)^2\left(\frac{f_w}{0.1}\right)\left(\frac{\dot{m}_{}}{10^3M_{}\mathrm{yr}^1}\right)\left(\frac{r_{}}{10\mathrm{R}_{}}\right)^{3/2}\left(\frac{m_{}}{100M_{}}\right)^{1/2}\mathrm{cm}^3,$$
(33)
where $`f_w\dot{m}_w/\dot{m}_{}`$ and we have assumed the wind is launched from the inner part of the disk from an annulus of inner radius $`f_rr_{}=2r_{}`$, and outer radius $`2f_rr_{}`$.

Assuming an isotropic flux, the ionizing luminosity necessary to balance recombinations and fully ionize the gas, is given by

$$\frac{S_{\mathrm{crit}}}{4\pi }=_{f_rr_{}}^{2f_rr_{}}\alpha ^{(2)}n_{\mathrm{H},\mathrm{base}}n_er^2𝑑r,$$
(34)
where $`\alpha ^{(2)}=2.62\times 10^{13}\mathrm{s}^1`$ is the radiative recombination rate to the excited states of hydrogen at $`10^4\mathrm{K}`$ and $`n_e`$ is the electron number density. We assume singly ionized He, so that

$$S_{\mathrm{crit}}=2.9\times 10^{47}f_v^2\left(\frac{f_r}{2}\right)^1\left(\frac{f_w}{0.1}\right)^2\left(\frac{\dot{m}_{}}{10^3M_{}\mathrm{yr}^1}\right)^2\left(\frac{m_{}}{100M_{}}\right)^1\mathrm{ph}\mathrm{s}^1.$$
(35)
Models of Population III stellar structure (Tumlinson & Shull 2000; Baraffe, Heger, & Woosley 2001; Bromm, Kudritzski, & Loeb 2001; Schaerer 2002) yield zero age main sequence (ZAMS) ionizing luminosities that can be approximated via $`S=9.42\times 10^{49}(m_{}/30M_{})^{2.5}\mathrm{ph}\mathrm{s}^1`$ for ($`12<m_{}/M_{}<40`$), $`S=9.17\times 10^{49}(m_{}/100M_{})^{1.7}\mathrm{ph}\mathrm{s}^1`$ for ($`40<m_{}/M_{}<120`$) and $`S=3.75\times 10^{50}(m_{}/300M_{})^{1.2}\mathrm{ph}\mathrm{s}^1`$ for ($`120<m_{}/M_{}<1000`$). The critical ZAMS mass that can ionize the outflow from a protostar accreting at the fiducial rate (Eq. 3) is then given by

$$m_{,\mathrm{crit}}=34f_v^{0.46}\left(\frac{f_r}{2}\right)^{0.23}\left(\frac{f_w}{0.1}\right)^{0.46}ϵ_d^{0.66}K^{0.98}M_{},$$
(36)
so long as $`12<m_{,\mathrm{crit}}/M_{}<40`$. Similar relations can be derived for when $`m_{,\mathrm{crit}}`$ is in a higher mass range.

Primordial protostars only reach the ZAMS at masses around $`100M_{}`$ (Omukai & Palla 2001; 2003; TM). Before this the ionizing luminosity is initially very small because the protostar is large and relatively cool. In the fiducial case the ionizing luminosity, including the contribution from accretion luminosity, then increases rapidly to become equal to the ZAMS value at about $`40M_{}`$ (Figure 7 of TM). Therefore, we expect the outflow to be able to confine ionizing radiation for much of the pre-main sequence phase of protostellar evolution up to about $`40M_{}`$, but not after that. These results have implications for feedback models that are based on the ionization of the protostellar disk (e.g. Hollenbach et al. 1994). In present-day massive star formation this process of confinement by the outflow can create extremely compact H II regions around massive protostars (Tan & McKee 2003b). A more complete analysis of the feedback effects in the primordial case, including those from FUV photons, is beyond the scope of this paper and is deferred to a future study.

## 5 Conclusions

We have presented calculations of the radial structure of primordial protostellar disks, using realistic models for the rate of accretion and the size of the protostar (Tan & McKee 2003a). We have followed the evolution of disk properties as the protostar grows in mass and changes in size. These disks differ from those around present-day protostars, because of the lack of cooling due to dust and because of the high rates at which they are fed matter from the collapsing core. The high accretion rates and large protostellar sizes require inclusion of the ionization energy in the structure equations, the effect of which is to reduce the disk temperatures and luminosities, particularly for disks around protostars with masses $`20M_{}`$. We argue that fragmentation in the disks is unlikely to occur, if local gravitational instabilities can provide an effective viscosity of the magnitude seen in the simulations of Gammie (2001).

We have shown that in situ dynamos can produce dynamically significant large scale magnetic fields in primordial star-forming disks. For present-day disks, such fields would likely be important for driving outflows in all star-forming systems. For the primordial case, the disks of the lowest mass protostars would not have enough time to generate strong fields, but the majority of the star-forming parameter regime at higher masses leads to strong fields even for very conservative assumptions about the initial seed magnetic field. Thus models of primordial star formation should include the dynamical influence of magnetic fields.

We emphasize that for these conclusions to be valid, gravitational instability, rather than the MRI, must provide the initial source of turbulence in the disk. Even though the ionization fraction is high, the seed field is below the critical value required at the mode of fastest growth. This contrasts with the situation in present-day protostellar disks. A second requirement is that the turbulence in the disk must generate helicity. It is important for this stage of the evolution to be investigated with numerical simulations that can resolve the helicity.

If large scale magnetic fields can be amplified to their saturation strength, then their main effect is to drive a collimated wind. This reduces the star formation efficiency by diverting material from the accretion disk and by sweeping-up ambient gas from the polar regions, that otherwise would have collapsed. The star formation efficiency is significantly reduced once the protostar reaches about $`100M_{}`$. The outflows are probably strong enough to eject matter from the dark matter mini-halo and to start magnetizing the intergalactic medium. The presence of outflows also increases the density distribution of gas near the star in the region above the accretion disk. The full effects of this have not yet been worked out, but it appears likely that the propagation of ionizing photons is inhibited, at least while the protostar is in the earlier stages of contraction towards the main sequence. These results will be applied to a more detailed study of protostellar feedback in a forthcoming paper.

JCT is supported by a Spitzer-Cotsen fellowship from Princeton University and NASA grant NAG5-10811. EGB acknowledges support from DOE grant DE-FG02-00ER54600. We thank C. McKee, V. Pariev, and K. Omukai for discussions.