A&A 420, 957-974 (2004)
DOI: 10.1051/0004-6361:20035915

Observations of L1521F: A highly evolved starless core

A. Crapsi^1,2 - P. Caselli³ - C. M. Walmsley³ - M. Tafalla⁴ - C. W. Lee⁵ - T. L. Bourke⁶ - P. C. Myers²

1 - Università degli Studi di Firenze Dipartimento di Astronomia e Scienza dello Spazio, Largo E. Fermi 5, 50125 Firenze, Italy
2 - Harvard-Smithsonian Center for Astrophysics, 60 Garden Street, Cambridge, MA 02138, USA
3 - INAF-Osservatorio Astrofisico di Arcetri, Largo E. Fermi 5, 50125 Firenze, Italy
4 - Observatorio Astronómico Nacional (IGN), Alfonso XII, 3, 28014 Madrid, Spain
5 - Taeduk Radio Astronomy Observatory, Korea Astronomy Observatory, 61-1 Hwaam-dong, Yusung-gu, Daejon 305-348, Korea
6 - Harvard-Smithsonian Center for Astrophysics, Submillimeter Array Project, 645 N. A'ohoku Place, Hilo, HI 96720, USA

Received 19 December 2003 / Accepted 19 March 2004

Abstract
We observed the pre-stellar core L1521F in dust emission at 1.2 mm and in two transitions each of  ${\rm N}_2{\rm H}^+$, ${\rm N}_2{\rm D}^+$, ${\rm C}^{18}{\rm O}$ and  ${\rm C}^{17}{\rm O}$ in order to increase the sample of well studied centrally concentrated and chemically evolved starless cores, likely on the verge of star formation, and to determine the initial conditions for low-mass star formation in the Taurus Molecular Cloud. The dust observation allows us to infer the density structure of the core and together with measurements of CO isotopomers gives us the CO depletion. ${\rm N}_2{\rm H}^+$ and ${\rm N}_2{\rm D}^+$ lines are good tracers of the dust continuum and thus they give kinematic information on the core nucleus. We derived in this object a molecular hydrogen number density n(${\rm H}_2$ $)\sim 10^6$  $\hbox{{\rm cm}}^{-3}$ and a CO depletion factor, integrated along the line of sight, $f_{\rm D} \equiv 9.5\times10^{-5}/x_{\rm obs}({\rm CO})\sim 15$ in the central 20 $^{\prime \prime }$, similar to the pre-stellar core L1544. However, the  $N(\hbox{${\rm N}_2{\rm D}^+$ })/N(\hbox{${\rm N}_2{\rm H}^+$ })$ column density ratio is $\sim$0.1, a factor of about 2 lower than that found in L1544. The observed relation between the deuterium fractionation and the integrated CO depletion factor across the core can be reproduced by chemical models if  ${\rm N}_2{\rm H}^+$ is slightly (factor of $\sim$2 in fractional abundance) depleted in the central 3000 AU. The  ${\rm N}_2{\rm H}^+$ and  ${\rm N}_2{\rm D}^+$ linewidths in the core center are $\sim$0.3  $\hbox {{\rm km s}}^{-1}$, significantly larger than in other more quiescent Taurus starless cores but similar to those observed in the center of L1544. The kinematical behaviour of L1521F is more complex than seen in L1544, and a model of contraction due to ambipolar diffusion is only marginally consistent with the present data. Other velocity fields, perhaps produced by accretion of the surrounding material onto the core and/or unresolved substructure, are present. Both chemical and kinematical analyses suggest that L1521F is less evolved than L1544, but, in analogy with L1544, it is approaching the "critical'' state.

Key words: radio lines: ISM - radio continuum: ISM - submillimeter

1 Introduction

The pre-stellar core L1544 has recently been the subject of much study because while apparently in hydrostatic equilibrium, there are indications that it is close to the "critical'' state at which it will become gravitationally unstable and from which it will dynamically collapse (see discussions by Tafalla et al. 1998; Ciolek & Basu 2000; Caselli et al. 2002a,b). If correct, this is of fundamental importance because it defines the initial conditions for the formation of a protostar, which affects many theoretical studies of low mass star formation.

Clearly trying to extrapolate general trends from a single object is difficult and a larger number of L1544-like cores (preferably with the same external environment) should be studied. Unfortunately there are rather few other objects with similar properties due to the short timescale of this phase. According to the Ciolek & Basu (2000) model, for example, (contraction of a disk driven by ambipolar diffusion) L1544-like properties fit the model structure of a core at times  $3{-}8\times 10^4$ years prior to the collapse, after an evolution of  $2.6\times10^6$ years. Thus in that particular model, L1544 finds itself in the last few percent of its evolution prior to becoming a protostar. While this may be a somewhat too literal interpretation of the model results, it shows that "L1544-type cores'' should be relatively rare.

Further progress requires the definition of what is a L1544-like core. One answer is to use current estimates of dust emission and absorption selecting cores of dust extinction upwards of 50 mag. Another approach is to say that cores which show signs of infalling gas (as does L1544, see Williams et al. 1999; Tafalla et al. 1998) are "L1544-twins''. This latter indicator is complicated by the fact that at the high densities found in the nuclei of cores similar to L1544, many molecular species and in particular CO and CS freeze-out onto dust grain surfaces (see Kramer et al. 1999; Caselli et al. 1999; Bacmann et al. 2002; Bergin et al. 2002; Jørgensen et al. 2002; Tafalla et al. 2002); observing such tracers implies observing the low density surrounding envelope. However, recent studies indicate that species whose abundance is linked to that of molecular nitrogen such as N₂H⁺ and NH₃ (as well as their deuterated counterparts) do not condense out in the same fashion and hence can be used as tracers of the dense gas (Bergin & Langer 1997). The extent to which this is true is debatable but it is a useful hypothesis and substantiated by the general similarity of the spatial distributions seen for example in dust emission and in maps of  ${\rm N}_2{\rm H}^+$ (Tafalla et al. 2002, 2004).

Caselli et al. (2002a,b) have used  ${\rm N}_2{\rm H}^+$ and  ${\rm N}_2{\rm D}^+$ to derive the physical, chemical and kinematical properties of L1544. They found that L1544 has a central  ${\rm N}_2{\rm H}^+$ column density of  $1.5 \times 10^{13}$ cm^-2 and a column density ratio  $N(\hbox{${\rm N}_2{\rm D}^+$ })/N(\hbox{${\rm N}_2{\rm H}^+$ })$ of 0.24. The ${\rm N}_2{\rm H}^+$ linewidths towards the nucleus (the dust emission peak) are roughly 0.3 km s^-1 and decrease as one goes to positions away from the center. The line of sight velocity measured in  ${\rm N}_2{\rm H}^+$ (1-0) and  ${\rm N}_2{\rm D}^+$ (2-1) shows a gradient along the minor axis of the elliptical structure seen in 1.3 mm dust emission but no clear gradient along the major axis.

In this paper, we will study another core in the Taurus complex,  L1521F (at an assumed distance of 140 parsec) using the same approach as in our study of L1544. Repeating the L1544 study carried out by Caselli et al. (2002a,b) is important because it allows us to check to what extent L1544 is an exceptional case. In order to do this we need another source which has the same general characteristics as L1544. The source selection was made using some preliminary results we obtained in a survey carried out at the IRAM-30 m telescope.  L1521F stood out as being the only core in Taurus, besides L1544, with strong  ${\rm N}_2{\rm D}^+$ (2-1) emission compared to  ${\rm N}_2{\rm H}^+$ (1-0). This suggests enhanced deuterium fractionation implying an advanced evolutionary state (Caselli et al. 2002b). Previous observations of this object have been carried out by Mizuno et al. (1994), Onishi et al. (1996), Codella et al. (1997), and Lee et al. (1999a). Onishi et al. (1999) also studied L1521F (which they call MC 27), and found a high central density, suggesting that this is the most evolved starless condensation in Taurus. L1521F was also noted by Lee et al. (1999b) as a strong infall candidate, in their survey of CS and ${\rm N}_2{\rm H}^+$ lines in starless cores, although later mapping of the two tracers has shown extended "red'' asymmetry in the CS(2-1) profiles (Lee et al. 2001).

In Sect. 2 of this paper, we describe our observational procedure. In Sect. 3 we present the observational results deriving the physical characteristics of the source and analysing its chemical and kinematical properties. In Sect. 4 we discuss the observational results and the summary can be found in Sect. 5.

2 Observations

The observations were carried out between April 2002 and January 2003 at the IRAM-30 m in three different runs.

In April 2002, we observed the core in ${\rm N}_2{\rm H}^+$ (1-0), ${\rm N}_2{\rm H}^+$ (3-2), ${\rm N}_2{\rm D}^+$ (2-1) and  ${\rm N}_2{\rm D}^+$ (3-2). In general, we used the symmetric frequency switch mode and the facility autocorrelator; in Table 1, we summarize the main observational parameters. The frequencies of the ${\rm N}_2{\rm H}^+$ (1-0), ${\rm N}_2{\rm D}^+$ (2-1) and  ${\rm N}_2{\rm D}^+$ (3-2) have been updated following the recent determinations of Dore et al. (2004); the values in the table refer to the F₁   F = 2 3  $\rightarrow$ 1 2, 2 3  $\rightarrow$ 1 2, and 4 5  $\rightarrow$ 3 4 hyperfine components of the ${\rm N}_2{\rm H}^+$ (1-0), ${\rm N}_2{\rm D}^+$ (2-1) and ${\rm N}_2{\rm D}^+$ (3-2) transitions, respectively. For  ${\rm N}_2{\rm H}^+$ (3-2) we used the frequency of the 2 1  $\rightarrow$ 1 0 component as determined by Caselli et al. (2002a). In the case of the April 2002 ${\rm N}_2{\rm D}^+$ (3-2) observations, in order to improve the baseline quality, we used also the "Wobbler switching'' mode with a 240 $^{\prime \prime }$ throw. We reached an rms sensitivity in main-beam brightness units of about 100 mK in all lines except ${\rm N}_2{\rm H}^+$ (3-2) ($\sim$400 mK). The pointing was checked every 2 h by means of a 3 or 2 mm continuum scan on nearby quasars and was accurate to within $\sim$4 $^{\prime \prime }$.

In order to refine the maps, originally taken with a 20 $^{\prime \prime }$ spacing, we observed in Nov. 2002 with a 10 $^{\prime \prime }$ grid (but 5 $^{\prime \prime }$ spacing in the inner 20 s).

Table 1: Telescope settings and parameters.

Between January 2003 and March 2003, we obtained continuum data at 1.2 mm together with observations of  ${\rm C}^{18}{\rm O}$ and  ${\rm C}^{17}{\rm O}$ (1-0) and (2-1). These data were taken in service mode by the IRAM staff.

The continuum data were obtained using MAMBO II, the 117-channels bolometer available at the 30 m. We mapped the core within an area of  $150\hbox{$^{\prime \prime}$ }\times 150\hbox{$^{\prime \prime}$ }$ scanning in azimuth with a 5 $^{\prime \prime }$/s speed and an interval between the subscans of 8 $^{\prime \prime }$. The atmospheric attenuation was measured to be 0.14 based on tipping curves measured after the map. The data were calibrated using the sources HL Tau and L1551 for which we assumed fluxes of 0.9 Jy and 1.4 Jy respectively and the final sensitivity was 5 mJy per 10.5 $^{\prime \prime }$ beam. The calibration error inherent in this comparison is likely to be at least ten percent due to both atmospheric fluctuations and calibration errors.

The ${\rm C}^{18}{\rm O}$ data were taken using the on-the-fly technique. We simultaneously mapped the ${\rm C}^{18}{\rm O}$ (1-0) and  ${\rm C}^{18}{\rm O}$ (2-1) using both polarizations for each line. The area covered was  $150\hbox{$^{\prime \prime}$ }\times 150\hbox{$^{\prime \prime}$ }$ and was scanned in the Right Ascension direction; the distance between the subscans was 5 $^{\prime \prime }$ as was the angular separation between two successive dumps. We also obtained a 9 points map in ${\rm C}^{17}{\rm O}$ (1-0) and  ${\rm C}^{17}{\rm O}$ (2-1) centered at the dust peak and spaced by 20 $^{\prime \prime }$.

   
3 Results

   
3.1 Integrated intensity maps and continuum emission

	Figure 1: Dust emission from L1521F. Levels are 30, 55 and 80 mJy/beam. Reference position is RA: 04:28:39.8 DEC: 26:51:35 (J2000). The dashed ellipse best fits the core structure with a 2D Gaussian. The dotted ellipse is the result of a 2D-Gaussian fit to the whole map, including the more extended emission. The black rectangle shows the area mapped in line observations (see Fig. 2).
Open with DEXTER

Figure 2: L1521F integrated intensity maps in the observed molecular transitions. Contour levels are 45, 70, 95% of the relative peak in each map (whose values are: 5.9, 0.74, 1.1, 0.33, 2.5 and 2.3  $\hbox {{\rm K km s}}^{-1}$ in  ${\rm N}_2{\rm H}^+$ (1-0), ${\rm N}_2{\rm H}^+$ (3-2), ${\rm N}_2{\rm D}^+$ (2-1), ${\rm N}_2{\rm D}^+$ (3-2), ${\rm C}^{18}{\rm O}$ (1-0) and  ${\rm C}^{18}{\rm O}$ (2-1) respectively). The ${\rm N}_2{\rm H}^+$ (3-2) data were smoothed to 26 $^{\prime \prime }$ resolution to increase S/N and help the comparison with ${\rm N}_2{\rm H}^+$ (1-0); for the same reasons  ${\rm N}_2{\rm D}^+$ (3-2) was smoothed to 16 $^{\prime \prime }$ (the HPBW of the 30 m antenna at the frequency of the the ${\rm N}_2{\rm D}^+$ (2-1) line). We overlaid on each map the ellipse which best fits the dust emission structure (in dashed line) and the 80 mJy/beam contour (dotted line).

Open with DEXTER

We show the observed map of the 1.2 mm dust continuum in Fig.  1, and the map of integrated line intensity, obtained in the observed transitions of  ${\rm N}_2{\rm H}^+$, ${\rm N}_2{\rm D}^+$, and  ${\rm C}^{18}{\rm O}$, in Fig. 2. The reference position for these maps is (04:28:39.8, 26:51:35) in J2000 coordinates.

The bolometer map was reduced with the IRAM standard reduction program NIC. From Fig. 1, we see that the observed emission has a "cometary'' structure in the sense that the low level contours are well-fitted by an ellipse with the maximum offset from the center. We can fit the general elliptical structure of the core with a 2D-Gaussian centered at offset (-30 $^{\prime \prime }$, 20 $^{\prime \prime }$) with full width half-power dimensions of  $274\hbox{$^{\prime \prime}$ }\times 170\hbox{$^{\prime \prime}$ }$ and position angle 25$^\circ$. If the more extended emission is not included in the fit, the 2D-Gaussian is centered on the dust peak position at (-10 $^{\prime \prime }$, 0 $^{\prime \prime }$), has half-power dimensions of  $127\hbox{$^{\prime \prime}$ } \times 77 \hbox{$^{\prime \prime}$ }$, a position angle of 25$^\circ$ and an aspect ratio equal to 1.6. The peak intensity is 90 mJy/beam.

The spectral line data were reduced using the standard IRAM package CLASS. A summary of line parameters at the dust peak is given in Table  2 and the corresponding spectra are shown in Fig. 3. The spectra shown in Fig. 3 (as well as the values in Table 2) have been derived from data Gaussian-smoothed to a resolution of 26 $^{\prime \prime }$ in the cases of  ${\rm N}_2{\rm H}^+$ and  ${\rm N}_2{\rm D}^+$, but unsmoothed in the case of  ${\rm C}^{18}{\rm O}$. These are also the effective resolutions of the ${\rm N}_2{\rm H}^+$ maps shown in Fig. 2 (the two  ${\rm N}_2{\rm D}^+$ maps are smoothed to 16 $^{\prime \prime }$, the angular resolution at the (2-1) frequency).

Table 2: Line parameters derived from hyperfine structure fitting at position (-10, 0). The values refer to spatially averaged spectra (26 $^{\prime \prime }$ beam).

Figure 3: Spectra of ${\rm N}_2{\rm H}^+$ (1-0), ${\rm N}_2{\rm H}^+$ (3-2), ${\rm N}_2{\rm D}^+$ (2-1), ${\rm N}_2{\rm D}^+$ (3-2), ${\rm C}^{18}{\rm O}$ (1-0), and ${\rm C}^{18}{\rm O}$ (2-1) towards the dust emission peak in L1521F. A fit of the hyperfine pattern (or in the case of  ${\rm C}^{18}{\rm O}$ a Gauss fit) of the lines assuming identical excitation temperatures for all satellites is also shown (dotted) for comparison.

Open with DEXTER

One clear result from Fig. 3 is that in L1521F there is evidence for an "infall signature'' although less marked than in L1544, where it is attributed to extended infall onto the core (Williams et al. 1999) or to central infall in a depleted core nucleus (Caselli et al. 2002a). In fact, we see evidence for asymmetric profiles, with the blue peak brighter than the red peak, in ${\rm N}_2{\rm H}^+$ (1-0) (where we derive an optical depth of order 4 in the main component based on our fit to the hyperfine satellites), but the two peaks are not clearly separated as in the case of L1544.

Figure 4: Enlargement of the ${\rm N}_2{\rm H}^+$ (3-2) spectrum of Fig. 3 toward offset (-10, 10), which shows tentative evidence for a double-peaked feature. The top panel show the line profile (grey histogram) and the hfs fit (black curve) assuming two velocity components; the bottom panel show the hfs fit assuming one velocity component. The velocities of the two components concide with the LSR velocities of clumps N1 and N2 observed by Shinnaga et al. (2004) in their interferometric observations, see Fig. 11 for reference.

Open with DEXTER

Table 3: Angular dimensions derived from two dimensionsal Gauss fits to maps in the 1.2 mm continuum, ${\rm N}_2{\rm H}^+$, and ${\rm N}_2{\rm D}^+$ emission.

It is interesting to note that two peaks are present in the spectrum of ${\rm N}_2{\rm H}^+$ (3-2) line toward the offset (-10, 10), as shown in Fig. 4. The hyperfine structure (hfs) fit to the ${\rm N}_2{\rm H}^+$ (3-2) line, assuming the presence of two velocity components along the line of sight, gives $V_{\rm LSR ~ 1} = 6.35\pm0.1$  $\hbox {{\rm km s}}^{-1}$ and  $V_{\rm LSR ~ 2}= 6.55\pm 0.04$  $\hbox {{\rm km s}}^{-1}$. The line widths of the two components are $\Delta v_1 = 0.15\pm 0.02$  $\hbox {{\rm km s}}^{-1}$ and  $\Delta v_2 = 0.19 \pm0.06$  $\hbox {{\rm km s}}^{-1}$, marginally (factor $\la$ 1.5) broader than the thermal  ${\rm N}_2{\rm H}^+$ line width at 10 K. Although this result should be confirmed with higher sensitivity data, the velocities are consistent with those of the N1 and N2 ( ${\rm N}_2{\rm H}^+$) clumps observed by Shinnaga et al. (2004) with interferometric observations. The hfs fit to the ${\rm N}_2{\rm H}^+$ (1-0) line, assuming two velocity components and fixing the velocities to the values found with the (3-2) line, is consistent with the observed spectrum (if the two velocities are not fixed, the hfs fitting procedure applied to the ${\rm N}_2{\rm H}^+$ (1-0) spectrum does not converge, probably because the two components are not well separated). This suggests that the two clumps observed by Shinnaga et al. (2004) in ${\rm N}_2{\rm H}^+$ (1-0) are also present in our single dish data, although the ${\rm N}_2{\rm H}^+$ (1-0) extended emission partially hides the features arising from the central region.

Figure 2 shows that L1521F has in common with several other cores (see also Tafalla et al. 2002; Caselli et al. 2002a,b) the property that while the maps of the nitrogen bearing molecules have a similar appearance to those observed in dust emission, maps in rare CO isotopomers such as  ${\rm C}^{18}{\rm O}$ show no correlation with the dust. In fact, the distribution of  ${\rm C}^{18}{\rm O}$ (1-0) and (2-1) is essentially flat over the area we have mapped, suggesting perhaps that the layer sampled in  ${\rm C}^{18}{\rm O}$ is either in the background or foreground relative to that seen in dust emission.

We have been able to make some estimate of the optical depths of the observed  ${\rm C}^{18}{\rm O}$ lines by comparison with our 9-point ${\rm C}^{17}{\rm O}$ map (20 $^{\prime \prime }$ spacing centred on (-5, -5)). From the integrated intensity ratio of the ${\rm C}^{18}{\rm O}$ and  ${\rm C}^{17}{\rm O}$ (1-0) lines and assuming an intrinsic abundance ratio [ ${\rm C}^{18}{\rm O}$]/[ ${\rm C}^{17}{\rm O}$] of 3.65 (from [¹⁸O]/[¹⁷O] = 3.65, Penzias 1981; see Kramer et al. 1999 for a discussion of the validity of these techniques), we can put in general upper limits on the optical depth of ${\rm C}^{18}{\rm O}$ (1-0) of order unity. At the (-5, -5) and (-25, -5) offsets, there is however evidence that ${\rm C}^{18}{\rm O}$ (1-0) is optically thick and in fact we derive ${\rm C}^{18}{\rm O}$ (1-0) optical depths of 1.0 and 1.5 respectively at these two offsets (errors 50%). In these positions we used the ${\rm C}^{17}{\rm O}$ (1-0) data and, assuming $T_{\rm ex} = 10$ K (consistent with the observed ${\rm C}^{17}{\rm O}$ (2-1)/(1-0) line ratio) and [¹⁸O]/[¹⁷O] = 3.65, we found an average ${\rm C}^{18}{\rm O}$ column density of $\sim$ $2\times 10^{15}$ cm^-2 in the central 25 $^{\prime \prime }$.

It is also interesting that we find that ${\rm N}_2{\rm H}^+$ (1-0) and dust emission maps have similar sizes, whereas ${\rm N}_2{\rm D}^+$ maps are somewhat smaller. Fitting a 2D Gaussian to the ${\rm N}_2{\rm H}^+$ and ${\rm N}_2{\rm D}^+$ maps in the same fashion as for the dust emission, we find the parameters given in Table 3 for the angular sizes of  ${\rm N}_2{\rm H}^+$ and ${\rm N}_2{\rm D}^+$ in L1521F. We deduce from these data a linear size for L1521F seen in ${\rm N}_2{\rm H}^+$ (1-0) of roughly 14 000 AU and dimensions in ${\rm N}_2{\rm D}^+$ (2-1) of  $13~000 \times 6000$ AU. It is interesting to note the smaller sizes derived in the higher J transitions of both species, suggesting that they are sampling somewhat higher density gas. The 1.2 mm dust continuum on the other hand is more extended suggesting that either an increase in excitation or in abundances is causing the ${\rm N}_2{\rm D}^+$ and ${\rm N}_2{\rm H}^+$ (3-2) line emission to increase in relative strength toward the dust peak. Notable also are the large aspect ratios (1.6-1.9) observed in ${\rm N}_2{\rm D}^+$ suggesting an ellipsoidal or triaxial form for the high density core.

   
3.2 Density and mass distribution

We expect the observed continuum emission at mm wavelengths to be optically thin and can relate the 1.2 mm flux to the ${\rm H}_2$ column density under the approximation of constant dust emissivity and temperature:

$$N(\hbox{${\rm H}_2$ }) = \frac{S_{\rm 1.2~mm}}{\Omega_{\rm beam} ~ B_{\nu}(T)~ \kappa_{\rm 1.2~mm} \; m } ,$$

where $N(\hbox{${\rm H}_2$ })$ is the averaged ${\rm H}_2$ column density, $S_{\rm 1.2~mm}$ is the emitted 1.2 mm flux density from the core, $\Omega_{\rm beam}$ is the telescope beam solid angle, $\kappa_{\rm 1.2~mm}$ is the dust absorption coefficient per gram of gas at 1.2 mm, $B_{\nu}(T)$ is the Planck function at temperature T, and m the mean molecular mass. In our calculations we used a dust temperature $T_{\rm dust} = 10$ K, close to the gas temperature found by Codella et al. (1997) using ammonia data (9.1 K), a dust opacity $\kappa_{\rm 1.2~mm} = 0.005$ cm² g^-1 (André et al. 1996) and m = 2.33 amu (with uncertainties of typically 50%; see e.g. Bianchi et al. 2003 and Gonçalves et al. 2004). This expression can be integrated over the 1.2 mm map to derive the total gas mass.

In this fashion, we derive a core mass from the continuum measurements, within the large ellipse in Fig. 1, of  $5.5 \pm0.5$ ${M}_\odot$ corresponding to an integrated total flux of  $3.0\pm 0.3$ Jy. (The error is dominated by errors in calibration and baseline removal.) Within the 26 $^{\prime \prime }$ beam, used for our line measurements, the enclosed mass is 0.7 ${M}_\odot$. It is interesting to note that 5 ${M}_\odot$ and 0.8 ${M}_\odot$ are the average total mass and 1 Jeans mass in the Taurus Molcular Cloud complex (Goodwin et al. 2004).

For the purpose of comparison with model calculations, we have used the 1.2 mm continuum data to estimate the density distribution under the assumption of spherical symmetry. This inevitably involves a rough approximation since the L1521F core clearly is not spherically symmetric (the aspect ratio in the 1.2 mm continuum emission is 1.6, see Table 3) and would be better approximated with an ellipsoid.

Nevertheless, we have followed the technique adopted by Tafalla et al. (2002) and fit our data with a model of the form:

$$n (r) = \frac{n_{0}}{1 + \left( \frac{r}{r_0} \right)^\alpha} \cdot$$

We thus circularly averaged the 1.2 mm data around the peak and used a $\chi^2$ routine to fit the dust profile using the above model integrated along the line of sight and convolved with a 10.5 $^{\prime \prime }$ Gaussian. The result of the fitting procedure can be seen in Fig. 5. The parameters that best fit our data are  n₀= 10⁶ cm^-3, r₀ = 20 $^{\prime \prime }$ (corresponding to 0.014 parsec or 2800 AU) and $\alpha = 2.0$.

	Figure 5: Circularly averaged dust emission fitted with the expression n0/(1+(r/r0)$^{\alpha }$) convolved with a 10.5 $^{\prime \prime }$ Gaussian. The parameters that best fit the L1521F data are n0= 106 cm-3, r0 = 20 $^{\prime \prime }$and  $\alpha = 2$. Dots with error bars are binned data and the dashed curve shows the observational beam.
Open with DEXTER

   
3.3 CO freeze-out

The comparison between the dust continuum emission and the ${\rm C}^{18}{\rm O}$ integrated intensity map allows the determination of the amount of CO freeze-out onto dust grain surfaces, integrated along the line of sight. This is possible because the millimeter continuum data furnish $N(\hbox{${\rm H}_2$ })_{\rm 1.2~mm}$, the column density of molecular hydrogen across the core, assuming optically thin conditions. The same quantity is obtained from ${\rm C}^{18}{\rm O}$ lines ( $N(\hbox{${\rm H}_2$ })_{\rm CO}$), again under the assumption of optically thin conditions, and adopting the relation between CO and ${\rm H}_2$ valid in undepleted conditions ([CO]/[${\rm H}_2$ $] \simeq 9.5\times10^{-5}\equiv x({\rm CO})_{\rm can}$, the "canonical'' CO abundance value; Frerking et al. 1982). From the $N(\hbox{${\rm H}_2$ })_{\rm 1.2~mm}$/ $N(\hbox{${\rm H}_2$ })_{\rm CO}$ column density ratio, the integrated CO depletion factor ($f_{\rm D}$), or the ratio between the canonical and the observed CO abundance ( $x({\rm CO})_{\rm can}$/ $x({\rm CO})_{\rm obs}$), is easily derived.

In practice, this process requires the division of the 1.2 mm dust continuum emission ( ${\cal F}_{\rm 1.2~mm}$[mJy/beam $] = S_{\rm 1.2~mm}$[mJy]/ $\Omega_{\rm beam}$) map by the ${\rm C}^{18}{\rm O}$ integrated intensity ( $W_{\hbox{${\rm C}^{18}{\rm O}$ }}$) map. The map-division has been carried out using the IRAM image manipulation software GRAPHIC, after degrading the continuum map to the angular resolution of the ${\rm C}^{18}{\rm O}$ (1-0) observations (22 $^{\prime \prime }$) so that the integrated depletion factor can be expressed by:

 

$$f_{\rm D} x({\rm CO})_{\rm can} \times \frac{N(\hbox{${\rm H}_2$ })_{\rm 1.2~mm}}{N(\hbox{${\rm C}^{18}{\rm O}$ })} \frac{[\rm ^{18}O]}{[\rm ^{16}O]}$$ =

$$8.5\times 10^{-2} \frac{{\cal F}_{{\rm 1.2~mm}}({\rm mJy/22\hbox{... ...me}$ }\ beam})} {W_{\hbox{${\rm C}^{18}{\rm O}$ }}({\rm K ~ km ~ s^{-1}})}\cdot$$ = (1)

To derive Eq. (1) we have assumed a dust temperature $T_{\rm dust} = 10$ K, $\kappa_{\rm 1.2~mm} = 0.005$ cm² g^-1, and an abundance ratio  $[\rm ^{16}O]/[\rm ^{18}O] = 560$ (Wilson & Rood 1994).

There are several caveats to the above procedure. One is that the "canonical abundance'' appears to vary from cloud to cloud (Lacy et al. 1994; Alves et al. 1999; Kramer et al. 1999) and is roughly one third of the diffuse cloud carbon gas phase abundance (Sofia et al. 1997). Given that CO is expected to represent essentially all the gas phase carbon in molecular clouds, this suggests depletion of carbon in some form (not necessarily as solid CO) even on the outskirts of cores (we note that the direct study of CO solid state features in Taurus (Chiar et al. 1995) shows that solid CO towards 4 field stars in Taurus is less than 40% of the canonical gas phase CO abundance and is observed for extinctions $A_{\rm V}$ greater than 6 mag). Thus we conclude that in particular cases such as L1521F, it is quite possible that we are using a value of  $x({\rm CO})_{\rm can}$ which is a factor of order 2 too large or small thus influencing the values of $f_{\rm D}$ which we derive but not the trends over our map.

Another problem is that the values of $f_{\rm D}$ which we derive are integrated along the line of sight in a situation where the observed ${\rm C}^{18}{\rm O}$ emission forms in an outer shell whereas the dust emission mainly emanates from the dense core nucleus. As a consequence, we observe ${\rm C}^{18}{\rm O}$ mainly from foreground and background gas which (see Fig. 2) has an essentially flat distribution over the region mapped by us. One concludes that our values of $f_{\rm D}$ are strict lower limits to the CO depletion in the core nucleus from which dust emission is observed. It is also the case that in this situation, the form of our map of $f_{\rm D}$ is essentially that of the dust emission (as we indeed find, see Eq. (2)). The local distribution of the CO depletion factor (which we call  $F_{\rm D}$, see below) may differ greatly and be much more highly peaked than in our map. Nevertheless, the $f_{\rm D}$ values derived by us are direct observables and we have therefore used them for the purpose of correlating with parameters such as the observed deuterium fractionation. We have also used them for model comparisons (Sects. 4.1.1 and 4.1.2) bearing in mind the above discussion.

	Figure 6: CO depletion factor (full contours) against dust emission smoothed to a 22 $^{\prime \prime }$ beam (grey scale). $f_{\rm D}$ contours range between 6 and 15, in steps of 3. The maximum value of $f_{\rm D}$ is 18 and is located at offset (0, 9). Note the good correspondence between the distribution of  ${\rm N}_2{\rm D}^+$ (2-1) (dashed contours) and the CO depletion map.
Open with DEXTER

In Fig. 6 the $f_{\rm D}$ map is shown (thick contour) overlapped with the smoothed 1.2 mm map (grey scale) and the ${\rm N}_2{\rm D}^+$ (2-1) map (dashed contours). The CO depletion factor increases between 6 at the lowest contour of the 1.2 mm map ( ${\cal F}_{\rm 1.2mm}= 95$ mJy/22 $^{\prime \prime }$ beam) to 18 at the peak position (offset [0, 9]), which is 11 $^{\prime \prime }$ off the millimeter dust emission peak (offset [-8, +2]), where ${\cal F}_{\rm 1.2~mm} = 320$ mJy/22 $^{\prime \prime }$ beam). From Fig. 6 we immediately see that $f_{\rm D}$ correlates with the ${\rm N}_2{\rm D}^+$ emission (and deuterium fractionation, see Sect. 3.4) and the 1.2 mm dust flux. The good correlation between $f_{\rm D}$ and  ${\cal F}_{\rm 1.2~mm}$ is also evident in Fig. 7, where $f_{\rm D}$ is plotted as a function of  $N(\hbox{${\rm H}_2$ })_{\rm 1.2~mm}$ ($\equiv$ $4.25 \times 10^{20}\ {\cal F}_{\rm 1.2~mm}$ mJy/(22 $^{\prime \prime }$ beam), with the assumptions on  $T_{\rm dust}$ and  $\kappa_{\rm 1.2~mm}$ as described above). Thus, in L1521F, and with the caveats discussed above, $f_{\rm D}$is linearly dependent on the ${\rm H}_2$ column density ( $N(\hbox{${\rm H}_2$ })$); the best fit to the data in Fig. 7 (dashed curve) gives:

 

$$f_{\rm D} 1.5 ~ \left[ \frac{N(\hbox{${\rm H}_2$ })}{10^{22}\ {\rm cm}^{-2}} \right]^{0.9}\cdot$$ (2)

In retrospect, this is not surprising as $N(\hbox{${\rm C}^{18}{\rm O}$ })$ does not vary greatly over the map. We note that $f_{\rm D}$ is always $\ga$2 in the region traced by the present observations, suggesting that even at the core edges (where $N(\hbox{${\rm H}_2$ }) \sim 10^{22}$ cm^-2) CO molecules are partially ($\simeq$30%) frozen onto grain mantles. This result is consistent with the average value of CO depletion ($\simeq$25%) gauged from observations of solid CO features in the direction of background stars in the Taurus complex (Chiar et al. 1995). The thick curve in Fig. 7 is the result of a simple chemical model of a centrally concentrated cloud, where CO depletion is taken into account, and where the CO abundance has been integrated along the line of sight to obtain the observed $f_{\rm D}$ ( $f_{\rm D} = \int {F}_{\rm D}(r) ~ {\rm d}l / \int {\rm d}l$, with  $F_{\rm D}(r)$being the CO depletion factor within the core, thus function of the cloud radius r; see Sect. 4.1). We anticipate here that the $f_{\rm D}$ value at the cloud center is only a very small fraction ($\sim$1%) of $F_{\rm D}$ in the central few thousand AU of the core.

	Figure 7: Integrated CO depletion factor $f_{\rm D}$ against ${\rm H}_2$ column density as derived from the data in Fig. 6. The dashed curve is the linear least fit to the data, whereas the solid curve is the $f_{\rm D}$vs. $N({\rm H_2})$ relation found in the simple model described in Sect. 4.1.2.
Open with DEXTER

   
3.4 N₂H⁺ - N₂D⁺ column densities, deuterium fractionation, and volume densities

The ${\rm N}_2{\rm H}^+$ and ${\rm N}_2{\rm D}^+$ column densities have been determined using the "constant- $T_{\rm ex}$'' (CTEX) approximation, which reduces to simple analytic expressions (see Appendix A in Caselli et al. 2002b), and the Large Velocity Gradient (LVG) approximation. Both approaches give reasonable column density estimates as long as optical depths are small. When, as for example for ${\rm N}_2{\rm H}^+$ (1-0), the optical depths in the main hyperfine satellites are large, one is best (independent of method) to use the weakest of the satellites or alternatively the optical depth inferred from the intensity ratio of the weakest satellite to the strong main components. The errors in any case stem from the difficulties in estimating the optical depth of thick lines compounded with possible non-LTE effects for the hyperfine satellites (Caselli et al. 1995). Errors due to the estimate of the partition function (i.e. the fraction of the species in unobserved levels) appear to be less. We in general report column densities for ${\rm N}_2{\rm H}^+$ using the CTEX approach based on the integrated intensity of the weakest satellite of ${\rm N}_2{\rm H}^+$ (1-0) and for ${\rm N}_2{\rm D}^+$ assuming optically thin conditions. From comparison between the different approaches employed by us, we estimate the column density errors to be 30 percent.

We have also used LVG estimates to infer the density towards the peak and edges of our ${\rm N}_2{\rm H}^+$ map. Here, we assume a temperature of 10 K and have used rates from Green (1975) for collisions of He with ${\rm N}_2{\rm H}^+$. Based on the values in Table 2 (data smoothed to 26 $^{\prime \prime }$ spatial resolution) for the integrated intensities of  ${\rm N}_2{\rm H}^+$, we find $n(\hbox{${\rm H}_2$ }) = 4 \times 10^5$  $\hbox{{\rm cm}}^{-3}$ towards the dust peak of L1521F and  $2.5\times 10^5$  $\hbox{{\rm cm}}^{-3}$ 30 $^{\prime \prime }$ North (offset [-10, 30]). For ${\rm N}_2{\rm D}^+$ (2-1) and (3-2), the corresponding numbers are  $6\times10^5$  $\hbox{{\rm cm}}^{-3}$ and  $3.5\times 10^5$  $\hbox{{\rm cm}}^{-3}$, at the peak and (-10, 30) offset positions, respectively. These values are similar to the corresponding values for L1544 consistent with the idea that they have similar density distributions. The density estimates are somewhat smaller than estimates based on the dust emission, probably due to the different (factor of 2.4 lower) spatial resolution, and consistent perhaps with the idea that ${\rm N}_2{\rm H}^+$ is abundant in a shell outside but close to the dust peak. However, given that the LVG method assumes homogeneous conditions, the LVG-derived densities are averages along the line of sight, thus lower values are expected when compared to the continuum data analysis, which takes into account the core density structure.

	Figure 8: [ ${\rm N}_2{\rm D}^+$]/[ ${\rm N}_2{\rm H}^+$] in the central region as estimated from CTEX calculations (see Sect. 3.4). Levels are 0.066, 0.078, and 0.091. The dashed contour is the small ellipse of Fig. 2.
Open with DEXTER

The deuterium fractionation is directly estimated from the $N(\hbox{${\rm N}_2{\rm D}^+$ })/N(\hbox{${\rm N}_2{\rm H}^+$ })$ column density ratio ($\equiv$ $R_{\rm deut}$), and the $R_{\rm deut}$ map in L1521F, assuming CTEX conditions, is shown in Fig. 8. $R_{\rm deut}$ ranges between $\sim$0.02 at the core edge to 0.1 in a region about 20 $^{\prime \prime }$ in size and centered on the dust peak position. The peak value of  $R_{\rm deut}$ is about a factor of 2 smaller than that found in L1544 (Caselli et al. 2002b). We note that the column density of  ${\rm N}_2{\rm H}^+$ is similar in the two cores, and that the factor of 2 of difference in deuterium fractionation is due to the (factor of 2) larger  ${\rm N}_2{\rm D}^+$ column density in L1544. This suggests that, although the two cores are similar in structure, L1521F is probably slightly less evolved than L1544 (see Sect. 4).

   
3.5 Line width vs. impact parameter

Quiescent starless cores mapped in ${\rm NH}_{3}$(1,1) and ${\rm N}_2{\rm H}^+$ (1-0) lines typically show a "velocity coherent'' zone of nearly constant line width ($\Delta v$) within the half-maximum contour map, followed by a $\Delta v$ rise at larger distances from core center (e.g. Goodman et al. 1998). There are however exceptions to this general trend, as pointed out by Caselli et al. (2002c). In particular, L1544 shows a significant increase of  ${\rm N}_2{\rm H}^+$ and  ${\rm N}_2{\rm D}^+$ (but not H¹³CO⁺ and DCO⁺) line widths toward the center (factor of 1.5 within $\sim$50 $^{\prime \prime }$; Caselli et al. 2002a). This increase has been interpreted as evidence of infall in the central few thousands AU, where CO and related species (such as H¹³CO⁺ and DCO⁺) are heavily depleted. Indeed, the line-width increase is consistent with models of ambipolar diffusion (see Sect. 4.2).

Figure 9: ${\rm N}_2{\rm H}^+$ (1-0) intrinsic line-width vs. impact parameter. The big dots are averages of the data points within 15 $^{\prime \prime }$ bins. The solid curve is the $\Delta v - b$ relation derived from the velocity field predicted by the t3 model of Ciolek & Basu (2000) (see Caselli et al. 2002a), where a constant turbulent field across the cloud has been included, and the dashed curve is from the same model but accounting for the presence of a central molecular "hole'' of 2000 AU in size (see Sect. 4.2).

Open with DEXTER

The same trend has been observed in L1521F using ${\rm N}_2{\rm H}^+$ and  ${\rm N}_2{\rm D}^+$ lines, and in Fig. 9 we show the results obtained in ${\rm N}_2{\rm H}^+$ (1-0) lines (the linewidth corrected for optical depth, or the intrinsic linewidth, is plotted). The figure also shows two theoretical predictions which will be discussed in Sect. 4.2. The decrease of the intrinsic line width with impact parameter in L1521F, although not as steep as in L1544, is clear in Fig. 9, where the average value of $\Delta v$ within bins of 15 $^{\prime \prime }$ (see big dots) ranges between 0.3  $\hbox {{\rm km s}}^{-1}$ at the dust peak to 0.25  $\hbox {{\rm km s}}^{-1}$ at a projected distance of 80 $^{\prime \prime }$.

	Figure 10: Isolated component profile along the $\rm PA = 45\hbox {$^{\rm o}$ }$ axis passing through (-10, 0) ( left panel) and  $-45\hbox {$^{\rm o}$ }$ axis passing through (-10, 0) ( right panel). This figure shows how the linewidth increases going from the edge of the core to the center. Fits are Gaussian fits to this component.
Open with DEXTER

One might interpret this line width gradient as being due to increased optical depth towards the dust peak. This however seems unlikely as illustrated in Fig. 10, where the profiles of the isolated hyperfine component of the ${\rm N}_2{\rm H}^+$ (1-0) line along the 45$^{\rm o}$ and -45$^{\rm o}$ axes, passing through the dust peak position, are shown. If the line is self-absorbed toward the core center, the hyperfine components with the largest statistical weight will be broadened compared to the weakest ones, affecting the hfs fit. However, we performed Gaussian fits to all the hyperfines components, finding the same $\Delta v$ values within the uncertainties. In fact, a similar trend is also observed if one plots the linewidth of the $F_1 ~ F = 1~0 \rightarrow 1~1$ (or weakest) component of the ${\rm N}_2{\rm H}^+$ (1-0) line versus the impact parameter. This line, being thin and symmetric across the core, is not affected by self-absorption.

We believe that this is a characteristic of unstable or supercritical (e.g. Mouschovias & Spitzer 1976) cores on the verge of star formation, more briefly called pre-stellar cores, a term which we use to characterize the subset of starless cores undergoing central infall (see also Ward-Thompson et al. 2003).

   
3.6 Velocity field

Assuming that L1521F is in solid body rotation, we determined the magnitude ${\cal G}$ and the direction $\Theta$ of the corresponding velocity gradient following the procedure described in Goodman et al. (1993). The magnitude ${\cal G}$ ranges between 0.4 and 1 km s^-1 pc^-1 depending on the tracer used, and large variations are also obtained for the gradient direction $\Theta$(see Table 4, Cols. 2 and 3). There is no tendency for ${\cal G}$ to increase for higher density tracers, as observed in L1544 (Caselli et al. 2002a).

Table 4: Results of velocity gradient fits.

To investigate in more detail the internal motions of L1521F, we determined "local'' velocity gradients (see Caselli et al. 2002a), where  ${\cal G}_{\rm l}$ and  $\Theta_{\rm l}$ have been calculated in adjacent $3\times 3$-point grids of the maps. The results are shown in Figs. 11 and 12 for ${\rm N}_2{\rm H}^+$ and ${\rm N}_2{\rm D}^+$ maps, respectively. The arrows across the map indicate the magnitude and the direction of local gradients, and they are centered on the 9-point grid used to estimate the corresponding  ${\cal G}_{\rm l}$ and  $\Theta_{\rm l}$ values. From the figures it is clear that L1521F is not undergoing solid body rotation. The velocity structure is quite complex, showing portions of the core where the gradient direction changes rapidly. For example, if we concentrate on the ${\rm N}_2{\rm H}^+$ (1-0) map (see Fig. 11, top panel), which has the highest sensitivity, three converging velocity gradient patterns are clearly visible: (i) toward South-West in the Northern half, (ii) toward North-East in the SE quadrant, and (iii) toward North-West in the SW quadrant of the core. A similar structure is also present in the other maps (see Fig. 11, bottom panel, and Fig. 12).

Interestingly, the mean of the local-gradient magnitudes ( $\langle{\cal G}_{\rm l}\rangle$) increases going from ${\rm N}_2{\rm H}^+$ (1-0) to ${\rm N}_2{\rm D}^+$ (3-2) (see Col. 4 of Table 4), suggesting that internal motions, although complex ( $\Theta_{\rm l}$ widely varies across the cloud), tend to become more prominent toward the core center. Moreover, ${\rm N}_2{\rm D}^+$ gradients appear larger than those derived from ${\rm N}_2{\rm H}^+$. If the observed velocity structure is at least partially due to inward motions in the core center, the larger local-gradient magnitudes detected in ${\rm N}_2{\rm D}^+$ can be explained by ${\rm N}_2{\rm D}^+$ being a better tracer of the core central regions, as found in L1544 (Caselli et al. 2002a,b; see also Sect. 4.1). However, the magnitude of the local gradients is generally within a few km s^-1 pc^-1 (see Table 4), thus they are produced by small velocity variations ( $\delta v \simeq 0.02{-}0.05$  $\hbox {{\rm km s}}^{-1}$) within $\simeq$0.01 pc, the size of the grid where local gradients have been estimated (see also Sect. 4.2 for discussion).

Figure 11: Velocity gradient vectors in L1521F derived from the two  ${\rm N}_2{\rm H}^+$ maps. The integrated intensity of ${\rm N}_2{\rm H}^+$ (1-0) ( top) and ${\rm N}_2{\rm H}^+$ (3-2) ( bottom) is shown as grey scale. The thick contour marks the integrated intensity half-maximum contour. The arrows across the map are "local'' velocity gradients obtained in a $3\times 3$-point grid of positions centered on the corresponding grid, indicating the magnitude and the direction. The arrow in the bottom right of the figure is the "total'' velocity gradient derived from the whole set of available positions with $A/\sigma _A > 10$ (A being the integrated intensity). Note the large variation in magnitude and direction of the local gradients. The white circles in the top figure locate the positions of the four  ${\rm N}_2{\rm H}^+$ clumps detected by Shinnaga et al. (2004).

Open with DEXTER

	Figure 12: Same as Fig. 11 for ${\rm N}_2{\rm D}^+$ data. A similar pattern is present.
Open with DEXTER

We note that the velocity gradient found with the present ${\rm N}_2{\rm H}^+$ (1-0) observations is significantly smaller than that deduced by Shinnaga et al. (2004) using interferometric observations, probably because they resolve out the more extended emission, with smaller or almost opposite (see Shinnaga et al. 2004) velocity gradients, compared to the central regions.

   
4 Discussion

The previous sections described the results found in our analysis of L1521F. The density profile, the CO depletion factor, the deuterium fractionation and the velocity field in the core have been presented. We found several similarities between L1521F and L1544, including the density structure (with $n(\hbox{${\rm H}_2$ }) \simeq 10^6$  $\hbox{{\rm cm}}^{-3}$ at the center), the amount of integrated CO depletion ( $f_{\rm D}\simeq 15$) toward the dust peak, and the decrease of line width with increasing distance from the cloud center. Differences are however present in the amount of deuterium fractionation (factor of $\sim$2 less than in L1544) and in the velocity structure. In this section we will discuss these findings in view of our knowledge of the chemical and physical processes in dense cloud cores. The discussion is thus split in two subsections, one focussing on the chemistry and the other on the kinematics of L1521F.

   
4.1 Chemistry

In cold clouds, the exothermicity of the proton-deuteron exchange reaction

 

$$\hbox{${\rm H}_{3}^{+}$ }\ ~ + ~ {\rm HD} ~ \rightarrow ~ \hbox{${\rm H}_{2}{\rm D}^{+}$ }\ ~ + \hbox{${\rm H}_2$ }+ \Delta E ,$$ (3)

where $\Delta E/k = 230$ K (e.g. Millar et al. 1989), is responsible for the enhancement of the ${\rm H}_{2}{\rm D}^{+}$/ ${\rm H}_{3}^{+}$ abundance ratio well above the local elemental abundance of deuterium ( $1.5 \times 10^{-5}$; Oliveira et al. 2003). As a consequence, given that  ${\rm H}_{2}{\rm D}^{+}$ transfers its deuterium to neutral species such as CO and ${\rm N}_2$, producing ${\rm DCO}^+$ and  ${\rm N}_2{\rm D}^+$, the $N(\hbox{${\rm DCO}^+$ })/N(\hbox{${\rm HCO}^+$ })$ and  $N(\hbox{${\rm N}_2{\rm D}^+$ })/N(\hbox{${\rm N}_2{\rm H}^+$ })$ column density ratios reach the values observed towards dense cloud cores ($\simeq$0.02-0.2; e.g. Butner et al. 1995; Williams et al. 1998; Caselli et al. 2002b; Sect. 3.4).

However, the freeze-out of neutral species, such as CO, O, and ${\rm N}_2$, boosts the deuterium fractionation (e.g. Dalgarno & Lepp 1984). In fact, a decrease in the fractional abundance of gas phase CO leads to a decrease of the ${\rm H}_{3}^{+}$ and ${\rm H}_{2}{\rm D}^{+}$ destruction rates and to an increase (caused by the higher ${\rm H}_{3}^{+}$ abundance) of the ${\rm H}_{2}{\rm D}^{+}$ formation rate (e.g. Roberts & Millar 2000a; see reaction 3). An empirical relation between CO depletion and deuterium fractionation in prestellar cores has been determined by Bacmann et al. (2003) using doubly deuterated formaldehyde, whose abundance is also predicted to largely increase with CO freeze-out (Roberts & Millar 2000b; see also Roberts et al. 2003).

In the previous sections, we found that in L1521F, CO is depleted (with percentages ranging from 30% at the core edge to 93% at the center, deduced from the integrated CO depletion factor) and, similarly to L1544, the deuterium fractionation is large. The present estimates of  $R_{\rm deut}\equiv N(\hbox{${\rm N}_2{\rm D}^+$ })/N(\hbox{${\rm N}_2{\rm H}^+$ })$ as a function of distance from the dust peak allow us to study the relation between deuterium fractionation and CO freeze-out across a dense core for the first time and test chemical models. In Fig. 13, $R_{\rm deut}$ is plotted as a function of $f_{\rm D}$ for all the positions present in Fig. 8. We note a clear tendency for the deuterium fractionation to increase with integrated CO depletion factor.

Figure 13: Deuterium fractionation ( $R_{\rm deut} \equiv N(\hbox {${\rm N}_2{\rm D}^+$ })/N({\rm N_2H^{+}})$) as a function of integrated CO depletion factor. Filled squares are data from this paper and the filled circle with error bar is our result for the dust peak of L1544 (Caselli et al. 2002b). Thin full curves are predictions from a simple chemical model using different rate coefficients for the proton-deuteron exchange reaction (3): the "standard" rate of  $1.5\times 10^{-9}$ cm3 s-1and the "GHR" rate of  $3.5\times 10^{-10}$ cm3 s-1, recently determined by Gerlich et al. (2002). Dashed curves refer to the same models but with larger value of  $a_{\rm min}$, the minimum radius of dust grains in the adopted MRN grain-size distribution (see Sect. 4.1.1). The dotted curve is the best fit to the data, found if $k_{\rm HD} = 7.5\times 10^{-10}$ cm3 s-1. The thick curve is the result of a more comprehensive chemical model which takes into account the density structure of the core (see Sect. 4.1.2). Typical 1-$\sigma$ error bars are shown. To reproduce the $R_{\rm deut}$ value observed in L1544, ${\rm N}_2$ freeze-out and multiply deuterated forms of ${\rm H}_{3}^{+}$ have to be included in the model (dash-dotted curve).

Open with DEXTER

In the following, we will analyse the observational results with simple chemical models, to better understand the observed  $f_{\rm D}{-}N(\hbox{${\rm H}_2$ })$ and  $R_{\rm deut}{-}f_{\rm D}$ trends shown in Figs. 7 and 13, and the differences between L1521F and L1544. It will be interesting to see whether in L1521F there is evidence of the so-called "molecular hole'', or the region where all species heavier than helium are heavily ($\ga$98%) depleted from the gas phase, as deduced in the case of L1544. In fact, the recent detection of  ${\rm H}_{2}{\rm D}^{+}$ toward the L1544 dust peak is consistent with the presence of a molecular hole in the central $\sim$2800 AU (Caselli et al. 2003; Walmsley et al. 2004). Then, the ${\rm N}_2{\rm H}^+$ and  ${\rm N}_2{\rm D}^+$ emission maps should show a central hole or emission plateau of similar size, but this has not been observed perhaps because of the poor spatial coverage of the central region and the limited spatial resolution ($\sim$2000 AU; Caselli et al. 2002a). Such a "molecular hole'' was predicted by Caselli et al. (1999) and Caselli et al. (2002a) in an attempt to interpret the origin of double peaked optically thin lines for the case of L1544, including the thinnest ${\rm N}_2{\rm H}^+$ (1-0) hyperfine component. In the case of L1521F, molecular lines, although asymmetric, do not clearly show two separated peaks. Thus, if the line asymmetry is due to the presence of a molecular hole in a contracting core, the hole should have a smaller size than in the case of L1544. Evidence for  ${\rm N}_2{\rm H}^+$ depletion toward core centers has also been claimed in B68 (Bergin et al. 2002) and L1512 (Lee et al. 2003), two starless cores with central densities $\sim$10⁵  $\hbox{{\rm cm}}^{-3}$ and close to hydrostatic equilibrium, thus in a different chemical and dynamical phase compared to L1544 and L1521F.

   
4.1.1 A simple analytical model

The four thin curves in Fig. 13 refer to the outputs of a simple steady state analytical chemical model including ${\rm H}_{3}^{+}$, ${\rm H}_{2}{\rm D}^{+}$, ${\rm H}_2$, HD, CO, ${\rm N}_2$, ${\rm HCO}^+$, ${\rm DCO}^+$, ${\rm N}_2{\rm H}^+$, ${\rm N}_2{\rm D}^+$, electrons and negatively charged grains. The recombination on grain surfaces uses the rates from Draine & Sutin (1987), assuming that the grains are bare and that their abundance by number is given by the MRN (Mathis et al. 1977) distribution with upper cutoff radius of 0.25 $\mu$m and lower cutoff radius $a_{\rm min} = 50$ Å (standard case, solid curves), and that all grains are negatively charged (a more realistic value for the fraction of negatively charged grains may be $\sim$0.5; see Flower & Pineau des Forêts 2003). We also considered a larger  $a_{\rm min}$ (=500 Å, dashed curves) to roughly take into account the process of grain coagulation in dense cores (e.g. Ossenkopf & Henning 1994). As expected, larger grains cause a (slightly) higher deuterium fractionation, given that the number density of dust grains decreases and so does the recombination rate on grains of ${\rm H}_{2}{\rm D}^{+}$ molecular ions (further increasing the grain size does not significantly change the result, given that dissociative recombination becomes more important). This simple model also assumes that the electron fractional abundance x(e) varies with density as  $1.3 \times 10^{-5}~n(\hbox{${\rm H}_2$ })^{-0.5}$ (McKee 1989), and that $n(\hbox{${\rm H}_2$ }) = 9.9 \times 10^3 ~ f_{\rm D}^{1.4}$ (in the density range between  $3\times 10^4$  $\hbox{{\rm cm}}^{-3}$ and 10⁶  $\hbox{{\rm cm}}^{-3}$), as found from a linear least square fit to the observed $f_{\rm D}$ data in Fig. 6 and the $n(\hbox{${\rm H}_2$ })$ data derived in Sect. 3.2 from the 1.2 mm continuum observations. This allows us to find a simple relation for the deuterium fractionation as a function of $f_{\rm D}$ in "L1544-like'' cores, if one neglects multiply deuterated forms of  ${\rm H}_{3}^{+}$ (see below):

 

$$R_{\rm deut} \equiv \frac{[\hbox{${\rm N}_2{\rm D}^+$ }]}{[\hbox{... ...}$ }]}{1 + 2/3[\hbox{${\rm H}_{2}{\rm D}^{+}$ }]/[\hbox{${\rm H}_{3}^{+}$ }]} ,$$ (4)

and

 

$$\frac{[\hbox{${\rm H}_{2}{\rm D}^{+}$ }]}{[\hbox{${\rm H}_{3}^{+}... ...rm N}_2$ }} [\hbox{${\rm N}_2$ }] + k_{\rm e} [{\rm e}] + k_{\rm gr} [\rm gr]},$$ (5)

where $k_{\rm HD}$ is the rate coefficient of reaction 3 (see below), [HD]/[${\rm H}_2$ $] \equiv 2\times[$D]/[H $] = 3\times 10^{-5}$ (Oliveira et al. 2003), $k_{\rm CO} = 3.3\times10^{-9}$  $(T/10~{\rm K})^{-0.5}$ cm³ s^-1 is the rate coefficient of the reaction ${\rm H}_{2}{\rm D}^{+}$ + CO  $\rightarrow{\it products}$, $k_{{\rm N_2}} = 1.7\times 10^{-9}$ cm³ s^-1 is the rate coefficient for the reaction ${\rm H}_{2}{\rm D}^{+}$ + ${\rm N}_2$  $\rightarrow$ products, [${\rm N}_2$]/[${\rm H}_2$ $] = 2\times 10^{-5}$ (Lee et al. 1996), $k_{\rm e} = 5.5\times 10^{-7}$  $(T/10\ {\rm K})^{-0.65}$ cm³ s^-1 is the dissociative recombination rate of  ${\rm H}_{2}{\rm D}^{+}$ (see Caselli et al. 1998), and  $k_{\rm gr}$is the rate of recombination of ${\rm H}_{2}{\rm D}^{+}$ on negatively charged dust grains, which, following Draine & Sutin (1987), is given by:

$$k_{\rm gr} ({\rm cm^3 s^{-1}}) 1.6 \times 10^{-7} \frac{a_{\rm min}}{10^{-8}~{\rm cm}} \left( \frac{T}{10 ~ {\rm K}} \right)^{-0.5}$$ =

$$\times \left( 1 + 3.6 \times 10^{-4} \frac{T}{10~ {\rm K}} \frac{a_{\rm min}}{10^{-8}\ {\rm cm}} \right) \cdot$$ (6)

Note that $k_{\rm gr} \times [\rm gr]$ is equivalent to  $\alpha_{\rm gr}$ in Draine & Sutin (1987), where $[{\rm gr}] = n_{\rm gr}/n(\hbox{${\rm H}_2$ })$ is the fractional abundance of dust grains and  $n_{\rm gr} = \int n_{\rm gr}(a) {\rm d}a$. Introducing numerical values into Eq. (5), we obtain an expression for [ ${\rm H}_{2}{\rm D}^{+}$]/[ ${\rm H}_{3}^{+}$] which only depends on $f_{\rm D}$ and  $a_{\rm min}$:

 

$$\frac{[\hbox{${\rm H}_{2}{\rm D}^{+}$ }]}{[\hbox{${\rm H}_{3}^{+}... ...3~s^{-1}})} {7/f_{\rm D} + 2.8 + 1.6/f_{\rm D}^{0.7} + 37.8 ~ f(a_{\rm min})} ,$$ (7)

where

$$f(a_{\rm min}) = \left( \frac{k_{\rm gr}} {1.6 \times 10^{-7}~{\r... ...-1}}} \right) \left( \frac{a_{\rm min}} {10^{-8}~{\rm cm}} \right)^{-2.5} \cdot$$ (8)

The last term in the denominator of Eq. (7) can be neglected as soon as $a_{\rm min} = 0.05~\mu$m or larger, so that, for molecular ions, dissociative recombination becomes more important than recombination on grains. Hence, it can be neglected in dense clouds, where small grains are likely to be depleted on the surface of large grains or coagulated in larger fluffy structures (Ossenkopf 1993). Depletion itself causes  $a_{\rm min} \approx 0.01~\mu$m (Walsmley et al. 2004).

The value of $k_{\rm HD}$ typically used in chemical codes is  $1.5\times 10^{-9}$ cm³ s^-1, which we call the "standard rate'' (see e.g. Roberts et al. 2003). However, Gerlich et al. (2002) have recently measured a lower value for this rate ( $3.5\times 10^{-10}$ cm³ s^-1, the "GHR rate'') and in Fig. 13 results obtained with the "standard'' and "GHR'' rate are shown. The observed data points lie between the two curves, and the best fit (dotted curve in bottom figure) is obtained with $k_{\rm HD} = 7.5\times 10^{-10}$ cm³ s^-1, which may suggest that a rate slightly (factor or $\sim$2) larger than the one recently measured is probably needed to explain observations. However, one should note that in the case of L1544 ( $R_{\rm deut} \sim 0.24$ at  $f_{\rm D}\sim 10$), even the "standard rate'' fails to reproduce the large deuterium fractionation observed toward the dust peak. The only way to reach  $R_{\rm deut} \sim 0.24$ with this analytical model (dash-dotted curve in Fig. 13) is to allow a drop in the abundance of ${\rm N}_2$ in the central regions where  $f_{\rm D} \ge 10$ and include in the chemical scheme D₂H⁺(), which becomes abundant in heavily (CO, ${\rm N}_2$, and O) depleted regions (Roberts et al. 2003; Walmsley et al. 2004). Therefore, the difference in deuterium fractionation between L1521F and L1544 is likely to be due to a different evolutionary stage, with L1521F being less evolved than L1544 (see also Aikawa et al. 2003). The inclusion of the reaction  ${\rm H}_{2}{\rm D}^{+}$ + HD  $\rightarrow$ D₂H⁺ + ${\rm H}_2$ (Gerlich et al. 2002) in this model leads to an increase of the deuterium fractionation by a factor of 2-3. However, in this simple chemical scheme we did not include the so-called "back'' reactions due to ortho-${\rm H}_2$ (see Gerlich et al. 2002), which have the effect of reducing the deuteration degree (see Walmsley et al. 2004); this analysis is beyond the scope of the present paper.

   
4.1.2 A simple chemical model in a centrally concentrated core

The analytic calculation outlined in the previous section has at least two major defects. One is the assumption of no abundance variations within the dense core, which should be computed first in order to estimate molecular abundances as a function of radius and then determine the column densities via integration along the line of sight. The second is our neglect of reactions with species such as O which also act to limit deuterium fractionation. Here, we present a small toy model aimed at overcoming these problems, and already described in Caselli et al. (2002b). This model follows the (steady-state) chemistry of a spherically symmetric cloud with a density profile deduced from the 1.2 mm dust continuum emission (see Sect. 3.2), and constant temperature T= 10 K. Neutral species in the model are ${\rm H}_2$, CO, ${\rm N}_2$, and atomic oxygen. We follow depletion of these species onto dust grains and their desorption due to the cosmic ray impulsive heating of the dust, following the procedure by Hasegawa & Herbst (1993). The abundance of molecular ions such as  ${\rm HCO}^+$, ${\rm N}_2{\rm H}^+$ and corresponding deuterated isotopomers are given by the steady state chemical equations using the istantaneous abundances of neutral species. Analogously, the electron fraction  $x({\rm e})$ can be computed in terms of global estimates for the molecular and metallic ions and using the same istantaneous abundances of CO, ${\rm N}_2$, and O in the gas phase. The approach we have adopted is a simplified version of the reaction scheme of Umebayashi & Nakano (1990), which includes molecular ion recombination on grain surfaces using rates from Draine & Sutin (1987) (see Caselli et al. 2002b for details).

This model furnishes the abundances of gaseous species as a function of radius, and the corresponding column densities as a function of impact parameter are calculated taking into account the appropriate beam convolution to simulate the observations. We used the new value of the dissociative recombination rate for  ${\rm H}_{3}^{+}$ ( $4 \times10^{-7}$ cm³ s^-1 at 10 K) determined by McCall et al. (2003), assumed $k_{\rm HD} = 3.5 \times 10^{-10}$ cm³ s^-1 (see previous section) and varied other parameters (essentially the binding energies on grain surfaces of ${\rm N}_2$ and O). The model has been forced to reproduce within 10% the observed ${\rm C}^{17}{\rm O}$, ${\rm N}_2{\rm H}^+$, and  ${\rm N}_2{\rm D}^+$ column densities toward the dust peak position ( $5.5\times 10^{14}$, $1.6\times 10^{13}$, and  $1.5\times 10^{12}$ cm^-2, respectively), and to reproduce within a factor of 2 the observed column density profiles in the above molecules. The best fit binding energies for ${\rm N}_2$ and O (800 K and 750 K, respectively) are quite close to the values deduced from theoretical calculations and laboratory measurements (750 K and 800 K, for ${\rm N}_2$ and O, respectively; see discussion in Hasegawa et al. 1992 and Bergin & Langer 1997).

Figure 14: Fractional abundances, $n(i)/n(\hbox {${\rm H}_2$ })$, as a function of radius in L1521F. The symbol M+ refers to metals (see text). The adopted density profile is the one described in Sect. 3.2. This model reproduces the observed column densities of CO, ${\rm N}_2{\rm H}^+$, and  ${\rm N}_2{\rm D}^+$. The ${\rm N}_2{\rm H}^+$ abundance decreases towards the center by a factor of about 2, whereas ${\rm N}_2{\rm D}^+$ increases by a factor of $\sim$7 from core edge to core center. Note that the ${\rm H}_{2}{\rm D}^{+}$ abundance is predicted to be $\sim$ $3\times 10^{-10}$at the core center, a factor of 3 smaller than observed in L1544. In analogy with L1544, the ${\rm HCO}^+$ and ${\rm DCO}^+$ abundance profiles present a steep drop within about 5000 AU. The right axis refers to the CO depletion factor (see the curve labelled "$F_{\rm D}$''). The set of parameters used to obtain this model are: $E_{\rm D}$(CO) = 1210 K, $E_{\rm D}$(${\rm N}_2$) = 800 K, $E_{\rm D}$(O) = 750 K, $\zeta = 1.3\times 10^{-17}$ s-1, $a_{\rm min} = 5\times 10^{-6}$ cm, and  $x(\hbox {${\rm N}_2$ }) = 2\times 10^{-5}$, where  $E_{\rm D}(i)$ is the binding energy of species i and $\zeta$ is the cosmic-ray ionization rate (see Caselli et al. 2002b for details of the model).

Open with DEXTER

As illustrated in Fig. 14, the fractional abundance of  ${\rm N}_2{\rm H}^+$ decreases from about  $2\times 10^{-10}$ at the edge of the cloud to about  $1\times 10^{-10}$ at the center. The increase of the ${\rm H}_{2}{\rm D}^{+}$/ ${\rm H}_{3}^{+}$ abundance ratio toward the center boosts the formation of  ${\rm N}_2{\rm D}^+$, which presents an abundance increase of about one order of magnitude from the edge (where $x(\hbox{${\rm N}_2{\rm D}^+$ }) = 3\times10^{-12}$) to the central 4000 AU ( $x(\hbox{${\rm N}_2{\rm D}^+$ }) = 2\times10^{-11}$) of the core. Therefore, in the case of L1521F, we do not need to invoke any "molecular hole'' as in the case of L1544 (see Caselli et al. 2003), although the present data do not allow us to distinguish between the models with different amounts of ${\rm N}_2$ freeze out (and consequent  ${\rm N}_2{\rm H}^+$ depletion) within few thousands AU. On the other hand, recent observations of ortho- ${\rm H}_{2}{\rm D}^{+}$ toward this object (Caselli, van der Tak, Ceccarelli et al., in preparation) strongly suggest that the molecular hole in L1521F must be smaller than in L1544, given that the line is about two times less intense than in L1544. Indeed, we predict here an ${\rm H}_{2}{\rm D}^{+}$ abundance of 3$\times$10^-10 at radii less than 3000 AU, a factor of about 3 lower than in L1544 (see Caselli et al. 2003). This is another indication that L1544 is likely to be more evolved, and more centrally concentrated, than L1521F. This seems to contradict the observational evidence that the central densities in the two cores are quite similar. However, one should keep in mind that the central density values (and the density structure) in the two cores have been obtained assuming isothermal conditions. Allowing the temperature to decrease toward the center, as predicted by models of dense clouds heated by the external radiation field (e.g. Galli et al. 2002) one finds that lower values of the central temperature implies larger central densities (e.g. Evans et al. 2001; Zucconi et al. 2001). One possible solution to the L1521F/L1544 puzzle is that L1544 is colder and more centrally concentrated than L1521F and hence is more depleted in the nucleus. This is consistent with L1544 being more evolved and closer to the "critical'' state than L1521F.

We note that the chemical composition shown in Fig. 14 reproduces the observed $R_{\rm deut}{-}f_{\rm D}$ (see the thick curve in Fig. 13) and  $f_{\rm D}{-}N(\hbox{${\rm H}_2$ })$ (see the thick curve in Fig. 7) relations, without any need to change the value of  $k_{\rm HD}$ (see previous section). This demonstrates the importance of taking into account core structure and differential molecular freeze-out in chemical models. Moreover, note that the depletion factor within the cloud, $F_{\rm D}$, is significantly larger (more than two orders of magnitude) than the integrated CO depletion factor $f_{\rm D}$, which is "diluted'' along the line of sight (compare $F_{\rm D}$ in Fig. 14 with $f_{\rm D}$ in Fig. 7). Finally, the central value of the electron fraction is $\sim$ $5\times 10^{-9}$, implying an ambipolar diffusion time scale (see e.g. Shu et al. 1987) of $\sim$ $2.5 \times10^{13} x({\rm e}) = 1.2\times 10^{5}$ yr, only slightly (factor of $\sim$3) larger than the free-fall time scale, once again suggesting that the core is close to dynamical collapse (although not as close as L1544). As in L1544, the major ion in the core is H₃O⁺, which is due to the presence of a significant fraction of gaseous atomic oxygen in the model (see also Aikawa et al. 2001 for similar results). We should however point out that the more recent models of Aikawa et al. (2003), where surface chemistry is included, predict a much lower abundance of O in the gas phase, given that in this model, an O-atom sticking to a grain is converted to water, which remains on the surface (assuming desorption due to the cosmic-ray impulsive heating of dust grains; see Hasegawa & Herbst 1993). Low ionization potential elements, essentially S⁺, Mg⁺, Fe⁺, Si⁺, Na⁺, generically called "metals'' (M⁺ in Fig. 14), are assumed to freeze out onto dust grains at the same rate as CO molecules. For this reason, they are negligible carriers of positive charges within the core, in our model.

   
4.2 Kinematics

In order to facilitate the comparison between L1521F and L1544, we analyzed the kinematical properties of L1521F using the same models as in Caselli et al. (2002a). In particular, starting from the velocity field predicted by the ambipolar diffusion models of Ciolek & Basu (2000; hereafter CB), we have derived the linewidth gradient and the line profile in a disk-like geometry and compare the results with our observations.

From the analysis of ${\rm N}_2{\rm H}^+$ (1-0) data, we found that the line width, similarly to L1544, decreases with distance from the dust peak (see Fig. 9). As seen in Caselli et al. (2002a), this observational evidence is consistent with the predictions of the CB model. Here, we applied the kinematic analysis of L1544 (Caselli et al. 2002a) to the case of L1521F, assuming that the cloud has a disk-like shape and is centrally concentrated, with the center coincident with the 1.2 mm dust continuum map peak. From the core axial ratio, and following Eq. (1) of CB, the angle between the vertical axis of the model and the plane of the sky is found to be 18$^\circ$. The disk is contracting via ambipolar diffusion and the resultant velocity field is used as input in a model which computes synthetic profiles of optically thin lines for all lines of sight across the model disk (for details, see Caselli et al. 2002a). As in the case of L1544, we also assumed that the density profile and the radial velocity field is given by the CB model at time t = 2.66 Myr ($\equiv$t₃in CB), which best reproduces the continuum observations of both cores. Line broadening is both due to thermal motions ( $\sigma_{\rm th} = 0.05$  $\hbox {{\rm km s}}^{-1}$ at 10 K) and microturbulence described by a turbulent velocity  $\sigma_{\rm tu}$ independent of position.

Figure 15: Observed line profiles of the F1 F = 0 1 $\rightarrow$1 2 hyperfine component of ${\rm N}_2{\rm H}^+$ (1-0) (empty histograms) along the minor ( left panels) and major ( right panels) axis of L1521F compared with predictions based on the Ciolek & Basu (2000) velocity field at time t=t3 (full curves), convolved with a $\sigma = 0.11$  $\hbox {{\rm km s}}^{-1}$ Gaussian to reproduce the observed intrinsic linewidth of the central spectrum (see Fig. 9). The filled histograms are model line profiles computed assuming no turbulent or thermal broadening. The core center is defined by the dust peak, located at offset (-8, +2).

Open with DEXTER

Figure 15 shows the comparison between observed (hystogram) and synthetic (curves) profiles of the ${\rm N}_2{\rm H}^+$ (1-0) isolated hyperfine component ( $F_1F = 01\rightarrow$12) along the minor and major axes of the L1521F core. The (optically thin) weakest component (F₁F = 10 $\rightarrow$11) shows very similar profiles; so that we decided to present the higher sensitivity observations of the isolated component. We also considered two different conditions in the model: (1) the ${\rm N}_2{\rm H}^+$ abundance is constant throughout the core, so that the ${\rm N}_2{\rm H}^+$ column density simply follows the dust, and (2) there is a "hole'' (2000 AU in size) in the ${\rm N}_2{\rm H}^+$ distribution. The difference between the two synthetic profiles is not significant (only a 4% increase of the linewidth toward the center, in the presence of the molecular hole, see Fig. 9), after the convolution with a Gaussian with $\sigma =\sqrt{\sigma_{\rm th} + \sigma_{\rm tu}} = 0.11$  $\hbox {{\rm km s}}^{-1}$, needed to match the intrinsic linewidth toward the center, and thus only one profile is shown in Fig. 15. Also shown in the Fig. 15 are model profiles for the limiting case of no turbulent or thermal broadening (filled histograms). The thing to note is that the agreement with the data is mixed, in the sense that the predicted velocity gradient along the minor axis is observed but it is restricted to the south-west half of the axis. Along the major axis there is no clear gradient in the north-west half of the axis, as expected in absence of disk rotation, but a (0.06  $\hbox {{\rm km s}}^{-1}$) blue shift of the line is present toward south-east.

The synthetic profiles become narrower as we move away from the central 40 $^{\prime \prime }$ of the disk-like cloud, given that the inward velocity in the adopted t₃ model reaches its maximum (0.12  $\hbox {{\rm km s}}^{-1}$) at radius 0.025 pc (or 37 $^{\prime \prime }$) before rapidly dropping to 0 at the cloud edge (see Fig. 2 of CB). We have analysed this line width variation to check its consistency with the observed trend shown in Fig. 9. The solid curve in Fig. 9 is the CB prediction in the case of no central "hole'', whereas the dashed curve illustrates the effects on the line width of a molecular "hole'' in the central 2000 AU (i.e. a region where all heavy elements have condensed onto dust grains). The line width of model lines has been calculated as the second moment of the velocity profile and it has been "normalized'' to the value of the line width observed in the central position by convolving the purely-kinematic profiles with a Gaussian with $\sigma = 0.11$ km s^-1, which can be interpreted as the combination of thermal broadening plus a constant turbulent field across the cloud. The comparison between the curves and the big dots (the binned data, see Sect. 3.5) suggests that our data are consistent with the CB model. In particular, the correspondence between the solid curve and the data indicates that the molecular hole is not present in L1521F, in agreement with the chemical analysis (see previous section). On the other hand, the steeper $\Delta v$-b relation found in L1544 (see Fig. 5 of Caselli et al. 2002a) suggests that the molecular hole is likely to be present in that source, again in agreement with the chemical analysis (see e.g. Caselli et al. 2003).

As shown in Sect. 3.6, the kinematics of L1521F is also characterized by complex motions which may be the result of turbulence (see e.g. Burkert & Bodenheimer 2000) or accretion of material onto the core or a combination of both. Local gradients have been determined in order to quantify these motions, and we found that the magnitude of local gradient vectors tends to increase for higher density tracers. This suggests that turbulence (expected to dissipate more rapidly in denser environments), is probably not the driving source of the observed velocity field. Moreover, unresolved substructure may further complicate the velocity field, as suggested by the ${\rm N}_2{\rm H}^+$ clumps observed by Shinnaga et al. (2004) (see Fig. 11): with the exception of their clump "N4'', the kinematics of the other clumps (N1-N3) is consistent with the velocity field inferred from our local gradient maps (see also Sect. 3.1) and the (marginal, see Fig. 4) evidence for two velocity components in ${\rm N}_2{\rm H}^+$ (3-2), resembling clumps N1 and N2.

Figure 16: Intensity (thin lines) and velocity (thick lines) of  ${\rm N}_2{\rm H}^+$ (1-0) along four different cuts; in the top panel continuous lines refer to the minor axis (PA 115$^{\rm o}$) while the dashed line refers to a cut at PA 70$^{\rm o}$; in the bottom panel, the continuous lines show the behaviour along the major axis (PA 25$^{\rm o}$) and the dashed lines represent the other bisector (PA 160$^{\rm o}$).

Open with DEXTER

It is interesting to compare our data with predictions from the non-magnetic turbulent models of Ballesteros-Paredes et al. (2003), who derived velocity profiles for dense cloud cores. As an example, Fig. 16 shows the velocity cuts observed across L1521F. Within the half maximum contour of  ${\rm N}_2{\rm H}^+$ (1-0), the largest velocity variation observed is $\la$0.04  $\hbox {{\rm km s}}^{-1}$ on a scale of 0.027 pc, corresponding to a velocity gradient of $\la$1.5 km s^-1 pc^-1. On the other hand, Ballesteros-Paredes et al. (2003) found that the smallest velocity variation for their clump 13 at time t1 (see their Fig. 9) is $\ga$0.3  $\hbox {{\rm km s}}^{-1}$ within 0.15 pc or $\ga$2 km s^-1 pc^-1. Thus, current supersonic turbulent models predict velocity gradients which are somewhat too large. We also note that the reversal in the velocity gradient direction observed in L1521F (see bottom panel of Fig. 16) is not present in the model examples shown in Ballesteros-Paredes et al. (2003), but it is predicted by the turbulent models of Burkert & Bodenheimer (2000).

In the two starless cores L1498 and L1517B, Tafalla et al. (2004) found a good correlation between the distribution of CS and the distribution of the high velocity ${\rm N}_2{\rm H}^+$ (1-0). The authors argue that the high velocity wing of the ${\rm N}_2{\rm H}^+$ (1-0) lines comes from a gas shell that is being accreted by the starless core; therefore it has not experienced the high density in the core nucleus and hence its depletion factor is low. We repeated the same experiment in L1521F producing channel maps of the core in ${\rm N}_2{\rm H}^+$ (1-0) but we did not find any strong deviation from the total intensity distribution. However, we did find that the ${\rm N}_2{\rm H}^+$ (1-0) velocity pattern across the core is similar to the ${\rm C}^{18}{\rm O}$ distribution. In fact, as shown in Fig. 17, ${\rm N}_2{\rm H}^+$ (1-0) presents red-shifted velocities where CO is bright. The difference in magnitude of this effect here and in L1517B and L1498 could be due to a brighter core that saturates the emission from the high velocity wing.

Figure 17: ${\rm C}^{18}{\rm O}$ (1-0) integrated intensity map (greyscale; wedge units are K  $\hbox {{\rm km s}}^{-1}$) overlaid by the ${\rm N}_2{\rm H}^+$ (1-0) velocity map (contours). ${\rm N}_2{\rm H}^+$ velocity levels are from 6.37  $\hbox {{\rm km s}}^{-1}$ to 6.55  $\hbox {{\rm km s}}^{-1}$ spaced by 0.02  $\hbox {{\rm km s}}^{-1}$; higher velocity contours are represented as solid lines.

Open with DEXTER

A final thing to note is that Fig. 17, together with Fig. 11 (top figure), suggests the presence of a coherent velocity field resembling low-frequency spatial oscillations of the outer cloud layers around some equilibrium dynamical state as recently proposed by Lada et al. (2003) in the case of the starless core B68. However, L1521F is more massive than B68 and is approaching the critical state for the onset of collapse, so that a situation of near-equilibrium for L1521F may be due to a balance between gravity and a combination of magnetic and thermal forces. If this is the case, the normal modes for L1521F would be more complicated than in the purely thermal case of B68, since the spherical symmetry that is a fair approximation for a thermal-pressure supported cloud such as B68 will no longer be valid for a magnetically supported object (G. E. Ciolek, priv. comm.). In fact, we observe L1521F not to be spherically symmetric (see Table 3).

4.3 Other possible interpretation

Finally, in this section we discuss some more speculative interpretations of the evolutionary state of L1521F. Firstly, it seems that L1521F is likely to be in an earlier stage of evolution than L1544, which is probably colder ($\sim$7 K in the central $\sim$2000 AU, see e.g. Zucconi et al. 2001) and more centrally concentrated (central densities $\sim$10⁷  $\hbox{{\rm cm}}^{-3}$; e.g. Evans et al. 2001) than L1521F, and hence more depleted in the nucleus. For example, in the CB model, this implies that whereas the density profile of L1521F is consistent with the cloud model at time t₃=2.66 Myr, L1544 is closer to t₄=2.68 Myr, and thus the estimated age difference is roughly 20 000 yr. However, these age estimates should be taken with caution given that the two cores may have had a different formation and contraction history and that, in general, cores, like stars, can differ in properties such as size, mass, and temperature regardless of their relative age.

An alternative view can be based on the clumpy structure observed in L1521F (Shinnaga et al. 2004), but not in L1544 (e.g. Williams et al. 1999), using BIMA. If the more complex kinematics in L1521F is due to the unresolved clumpy substructure, one may speculate that L1521F is close to the formation of a group of low mass protostars, whereas L1544 is likely to form one or two stellar objects. Thus, the two cores may show chemical and physical properties characteristic of the initial conditions of different modes of star formation in low mass cores. One should also note that L1521F resides in the middle of the main filament of the Taurus Cloud, whereas L1544 is at the edges of the complex and hence the different environments of the two cores may be responsible for the different kinematics and chemical properties of apparently similar (in the dust continuum and ${\rm C}^{18}{\rm O}$ emission) dense cores.

5 Conclusions

We have analysed physical and chemical properties of L1521F, a starless core in the Taurus Molecular Cloud, with characteristics similar to the pre-stellar core L1544. The main similarities and differences between the two cores are listed below:

1. the dust emission distributions are similar, implying a fairly closely matched density structure, with central densities of  $n_0 \sim 10^6$  $\hbox{{\rm cm}}^{-3}$, the radius of the "flat'' region r₀ = 20 $^{\prime \prime }$, and similar asymptotic power index $\alpha$ (see Sect. 3.2 and Tafalla et al. 2002). In particular, the aspect ratio is quite similar: 1.6 and 1.8 in L1521F and L1544, respectively.

2. The line width decreases with distance from the cloud center ($\sim$0.3  $\hbox {{\rm km s}}^{-1}$) to 80 $^{\prime \prime }$ away ($\sim$0.25  $\hbox {{\rm km s}}^{-1}$; see Fig. 9), in analogy with L1544, and consistent with the predictions of ambipolar diffusion models, although any gravity-driven contraction in 2-3D is expected to get localized line broadening. The particular model which best fits the data will be investigated in the future.

3. The amount of CO freeze-out (integrated CO depletion factor  $f_{\rm D} = 15$) is also comparable to L1544, as is the column density of  ${\rm N}_2{\rm H}^+$ toward the dust peak ($\simeq$ $1.5\times 10^{12}$ cm^-2).

4. The deuterium fractionation toward the L1521F dust peak ( $R_{\rm deut}\sim 0.1$) is however a factor $\sim$2 smaller than in L1544, due to the (factor of 2) smaller column density of ${\rm N}_2{\rm D}^+$ toward L1521F. This can be understood if L1521F is less chemically evolved than L1544, with a smaller (r < 2000 AU) central molecular hole.

5. Unlike in L1544, the velocity field in L1521F maintains a complex structure even at the small scales traced by ${\rm N}_2{\rm D}^+$ and ${\rm N}_2{\rm H}^+$ (3-2) (see Figs. 11 and 12). This may be due to the presence of clumps in the central few thousand AU (as deduced by the interferometric observations of Shinnaga et al. 2004), but could also be caused by normal mode pulsations, as in the case of B68 studied by Lada et al. (2003). The ambipolar diffusion model with infall of Ciolek & Basu (2000) has problems in reproducing the whole velocity field observed across the core. This may be due to the fact that part of the observed bulk motions result from residual core contraction, as suggested by Tafalla et al. (2004) in their study of L1517B and L1498, thus preventing a clear determination of the velocity field within the core nucleus.

6. The line profiles in L1521F show asymmetric structure, although the two peaks are not well separated as in L1544. This is consistent with smaller (factor of 1.5) infall velocities in the central few thousand AU of L1521F.

In any case, the large central density ($\sim$10⁶  $\hbox{{\rm cm}}^{-3}$) and the evidence of central infall (broader line widths toward the center and central infall asymmetry in ${\rm N}_2{\rm H}^+$ (1-0)) indicate that L1521F is another starless core on the verge of star formation, or a pre-stellar core. Although a study of a more complete sample is needed, assuming that L1544 and L1521F are the only two cores in Taurus close to dynamical collapse, and given that the total number of starless cores in Taurus is 44 (Onishi et al. 2002), we can argue that this process of contraction towards the "critical'' stage, or the "L1544-phase'', may last about five percent of the core lifetime.

A more detailed analysis of ${\rm N}_2{\rm H}^+$ line profiles will be presented in a future paper, where the observed chemical abundances will be introduced in a non spherically symmetric Monte Carlo radiative transfer code. This is needed to understand the nature of the ${\rm N}_2{\rm H}^+$ (1-0) line asymmetry, which may be caused by self-absorption plus infall, or to infall plus a molecular hole, or to the relative motion of different clumps, or to a mixture of the above possibilities.

Acknowledgements

The authors are grateful to the staff of the IRAM 30 m antenna, for their support during observations, and to Richard Crutcher and Daniele Galli for discussion. We also thank the referee, Glenn Ciolek, for clarifying several statements in the paper and Floris van der Tak for a critical reading of the submitted manuscript. P.C. and C.M.W. aknowledge support from the MIUR project "Dust and Molecules in Astrophysical Environments''. C.W.L. acknowledges support from the Basic Research Program (grant KOSEF R01-2003-000-10513-0) of the Korea Science and Engineering Foundation.

References

• Aikawa, Y., Ohashi, N., & Herbst, E. 2003, ApJ, 593, 906 [NASA ADS] [CrossRef] (In the text)
• Aikawa, Y., Ohashi, N., Inutsuka, S., Herbst, E., & Takakuwa, S. 2001, ApJ, 552, 639 [NASA ADS] [CrossRef] (In the text)
• Alves, J., Lada, C. J., & Lada, E. A. 1999, ApJ, 515, 265 [NASA ADS] [CrossRef] (In the text)
• André, P., Ward-Thompson, D., & Motte, F. 1996, A&A, 314, 625 [NASA ADS] (In the text)
• Bacmann, A., Lefloch, B., Ceccarelli, C., et al. 2002, A&A, 389, L6 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Bacmann, A., Lefloch, B., Ceccarelli, C., et al. 2003, ApJ, 585, L55 [NASA ADS] [CrossRef] (In the text)
• Ballesteros-Paredes, J., Klessen, R. S., & Vázquez-Semadeni, E. 2003, ApJ, 592, 188 [NASA ADS] [CrossRef] (In the text)
• Bergin, E. A., Alves, J., Huard, T., & Lada, C. J. 2002, ApJ, 570, L101 [NASA ADS] [CrossRef] (In the text)
• Bergin, E. A., & Langer, W. D. 1997, ApJ, 486, 316 [NASA ADS] [CrossRef] (In the text)
• Bianchi, S., Gonçalves, J., Albrecht, M., et al. 2003, A&A, 399, L43 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Burkert, A., & Bodenheimer, P. 2000, ApJ, 543, 822 [NASA ADS] [CrossRef] (In the text)
• Butner, H. M., Lada, E. A., & Loren, R. B. 1995, ApJ, 448, 207 [NASA ADS] [CrossRef] (In the text)
• Caselli, P., Myers, P. C., & Thaddeus, P. 1995, ApJ, 455, L77 [NASA ADS] (In the text)
• Caselli, P., Walmsley, C. M., Terzieva, R., & Herbst, E. 1998, ApJ, 499, 234 [NASA ADS] [CrossRef] (In the text)
• Caselli, P., Walmsley, C. M., Tafalla, M., Dore, L., & Myers, P. C. 1999, ApJ, 523, L165 [NASA ADS] [CrossRef] (In the text)
• Caselli, P., Walmsley, C. M., Zucconi, A., et al. 2002a, ApJ, 565, 331 [NASA ADS] [CrossRef] (In the text)
• Caselli, P., Walmsley, C. M., Zucconi, A., et al. 2002b, ApJ, 565, 344 [NASA ADS] [CrossRef] (In the text)
• Caselli, P., Benson, P. J., Myers, P. C., & Tafalla, M. 2002c, ApJ, 569, 815 [NASA ADS] [CrossRef] (In the text)
• Caselli, P., van der Tak, F. F. S., Ceccarelli, C., & Bacmann, A. 2003, A&A, 403, L37 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Chiar, J. E., Adamson, A. J., Kerr, T. H., & Whittet, D. C. B. 1995, ApJ, 455, 234 [NASA ADS] [CrossRef] (In the text)
• Ciolek, G. E., & Basu, S. 2000, ApJ, 529, 925 [NASA ADS] [CrossRef] (CB) (In the text)
• Codella, C., Welser, R., Henkel, C., Benson, P. J., & Myers, P. C. 1997, A&A, 324, 203 [NASA ADS] (In the text)
• Dalgarno, A., & Lepp, S. 1984, ApJ, 287, L47 [NASA ADS] [CrossRef] (In the text)
• Dore, L., Caselli, P., Beninati, S., et al. 2004, A&A, 413, 1177 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Draine, B. T., & Sutin, B. 1987, ApJ, 320, 803 [NASA ADS] [CrossRef] (In the text)
• Evans, N. J., II, Rawlings, J. C. M., Shirley, Y., & Mundy, L. G. 2001, ApJ, 557, 193 [CrossRef] (In the text)
• Flower, D. R., & Pineau des Forêts, G. 2003, MNRAS, 343, 390 [NASA ADS] [CrossRef] (In the text)
• Frerking, M. A., Langer, W. D., & Wilson, R. W. 1982, ApJ, 262, 590 [NASA ADS] [CrossRef] (In the text)
• Galli, D., Walmsley, C. M., & Gonçalves, J. 2002, A&A, 394, 275 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Gerlich, D., Herbst, E., & Roueff, E. 2002, Planet. Space Sci., 50, 1275 [NASA ADS] [CrossRef] (In the text)
• Gonçalves, J., Galli, D., & Walmsley, C. M. 2004, A&A, 415, 617 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Goodman, A. A., Barranco, J. A., Wilner, D. J., & Heyer, M. H. 1998, ApJ, 504, 223 [NASA ADS] [CrossRef] (In the text)
• Goodman, A. A., Benson, P. J., Fuller, G. A., & Myers, P. C. 1993, ApJ, 406, 528 [NASA ADS] [CrossRef] (In the text)
• Goodwin, S. P., Whitworth, A. P., & Ward-Thompson, D. 2004, A&A, in press [arXiv:astro-ph/0403028] (In the text)
• Green, S. 1975, ApJ, 201, 366 [NASA ADS] [CrossRef] (In the text)
• Hasegawa, T. I., & Herbst, E. 1993, MNRAS, 263, 589 [NASA ADS] (In the text)
• Hasegawa, T. I., Herbst, E., & Leung, C. 1992, ApJS, 82, 167 [NASA ADS] [CrossRef] (In the text)
• Jørgensen, J. K., Schöier, F. L., & van Dishoeck, E. F. 2002, A&A, 389, 908 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Kramer, C., Alves, J., Lada, C. J., et al. 1999, A&A, 342, 257 [NASA ADS] (In the text)
• Lacy, J. H., Knacke, R., Geballe, T. R., & Tokunaga, A. T. 1994, ApJ, 428, L69 [NASA ADS] [CrossRef] (In the text)
• Lada, C. J., Bergin, E. A., Alves, J. F., & Huard, T. L. 2003, ApJ, 586, 286 [NASA ADS] [CrossRef] (In the text)
• Lee, H.-H., Bettens, R. P. A., & Herbst, E. 1996, A&AS, 119, 111 [NASA ADS] (In the text)
• Lee, J.-E., Evans, Neal, J., II, Shirley, Y. L., & Tatematsu, K. 2003, ApJ, 583, 789 [NASA ADS] [CrossRef] (In the text)
• Lee, C. W., Myers, P. C., & Tafalla, M. 1999b, ApJ, 526, 788 [NASA ADS] [CrossRef] (In the text)
• Lee, C. W., Myers, P. C., & Tafalla, M. 2001, ApJS, 136, 703 [NASA ADS] [CrossRef] (In the text)
• Lee, S.-W., Ohashi, N., Moriarty-Schieven, G., & Butner, H. M. 1999a, Proc. of Star Formation 1999, held in Nagoya, Japan, June 21-25, 1999, ed. T. Nakamoto, Nobeyama Radio Observatory, 179 (In the text)
• Mathis, J. S., Rumpl, W., & Nordsieck, K. H. 1977, ApJ, 217, 425 [NASA ADS] [CrossRef] (In the text)
• McCall, B. J., Huneycutt, A. J., Saykally, R. J., et al. 2003, Nature, 422, 500 [NASA ADS] [CrossRef] (In the text)
• McKee, C. F. 1989, ApJ, 345, 782 [NASA ADS] [CrossRef] (In the text)
• Millar, T. J., Bennet, A., & Herbst, E. 1989, ApJ, 340, 906 [NASA ADS] [CrossRef] (In the text)
• Mizuno, A., Onishi, T., Hayashi, M., et al. 1994, Nature, 368, 719 [NASA ADS] [CrossRef] (In the text)
• Mouschovias, T. Ch., & Spitzer, L., Jr. 1976, ApJ, 210, 326 [NASA ADS] [CrossRef] (In the text)
• Oliveira, C. M., Hébrard, G., Howk, J. C., et al. 2003, ApJ, 587, 235 [NASA ADS] [CrossRef] (In the text)
• Onishi, T., Mizuno, A., & Fukui, Y. 1999, PASJ, 51, 257 [NASA ADS] (In the text)
• Onishi, T., Mizuno, A., Kawamura, A., Ogawa, H., & Fukui, Y. 1996, ApJ, 465, 815 [NASA ADS] [CrossRef] (In the text)
• Onishi, T., Mizuno, A., Kawamura, K., Tachihara, K., & Fukui, Y. 2002, ApJ, 575, 950 [NASA ADS] [CrossRef] (In the text)
• Ossenkopf, V. 1993, A&A, 280, 617 [NASA ADS] (In the text)
• Ossenkopf, V., & Henning, Th. 1994, A&A, 291, 943 [NASA ADS] (In the text)
• Penzias, A. A. 1981, ApJ, 249, 518 [NASA ADS] [CrossRef] (In the text)
• Roberts, H., & Millar, T. J. 2000a, A&A, 361, 388 [NASA ADS] (In the text)
• Roberts, H., & Millar, T. J. 2000b, A&A, 364, 780 [NASA ADS] (In the text)
• Roberts, H., Herbst, E., & Millar, T. J. 2003, ApJ, 591, L41 [NASA ADS] [CrossRef] (In the text)
• Shinnaga, H., Ohashi, N., Lee, S.-W., & Moriarty-Schieven, G. H. 2004, ApJ, 601, 962 [NASA ADS] [CrossRef] (In the text)
• Shu, F. H., Adams, F. C., & Lizano, S. 1987, ARA&A, 25, 23 [NASA ADS] (In the text)
• Sofia, U. J., Cardelli, J. A., Guerin, K. P., & Meyer, D. M. 1997, ApJ, 482, L105 [NASA ADS] [CrossRef] (In the text)
• Tafalla, M., Mardones, D., Myers, P. C., et al. 1998, ApJ, 504, 900 [NASA ADS] [CrossRef] (In the text)
• Tafalla, M., Myers, P. C., Caselli, P., Walmsley, C. M., & Comito, C. 2002, ApJ, 569, 815 [NASA ADS] [CrossRef] (In the text)
• Tafalla, M., Myers, P. C., Caselli, P., & Walmsley, C. M. 2004, A&A, 416, 191 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Umebayashi, T., & Nakano, T. 1990, MNRAS, 243, 103 [NASA ADS] (In the text)
• van der Tak, F. F. S., Caselli, P., Walmsley, C. M., et al. 2003, in The Dense Interstellar Medium in Galaxies, Proc. of the fourth Cologne-Bonn-Zermatt Symposium, Sept. 22-26, 2003, in press [arXiv:astro-ph/0311125]
• Williams, J. P., Bergin, E. A., Caselli, P., Myers, P. C., & Plume, R. 1998, ApJ, 503, 689 [NASA ADS] [CrossRef] (In the text)
• Williams, J. P., Myers, P. C., Wilner, D. J., & di Francesco, J. 1999, ApJ, 513, L61 [NASA ADS] [CrossRef] (In the text)
• Zucconi, A., Walmsley, C. M., & Galli, D. 2001, A&A, 376, 650 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Walmsley, C. M., Flower, D., & Pineau des Forêts, G. 2004, A&A, 418, 1035 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Ward-Thompson, D., Kirk, J., & André, P. 2003, in Star Formation at High Angular Resolution, International Astronomical Union. Symposium No. 221, held 22-25 July, 2003 in Sydney, Australia, in press [arXiv:astro-ph/0310066] (In the text)
• Wilson, T. L., & Rood, R. T. 1994, ARA&A, 32, 191 [NASA ADS] (In the text)