text:
# Chemistry and Dynamics in Pre-Protostellar Cores

## 1 INTRODUCTION

Many studies have been done on low mass-star formation in detail, theoretically and observationally, over the past two decades. These studies were possible because many low mass protostellar cores are close, isolated, and geometrically simple. In addition, they evolve slowly, allowing the process to be observed in distinct stages. The information accumulated by those studies now provides the big picture of the evolutionary stages of low-mass star formation. Lada (1987), Myers & Ladd (1993) and André et al. (1993) developed an evolutionary sequence based on the spectral energy distribution (SED) of dust continuum emission through Class 0 to Class III sources, and Adams et al. (1987) supported the evolutionary sequence by modeling the SED theoretically. In addition, several groups (Larson 1969, Penston 1969, Shu 1977, Shu et al. 1987, Foster & Chevalier 1993, Ciolek & Mouschovias 1994, McLaughlin & Pudritz 1997) have developed dynamical evolution models that theoretically predict different structures of density and velocity in protostellar cores, based on differing initial conditions.

In spite of the big picture of low-mass star formation, the initial collapse is poorly understood. The collapse dynamics and the timescale of evolution are crucially dependent on initial conditions. Preprotostellar cores (hereafter PPCs) are believed to be gravitationally bound, but they have no central hydrostatic protostar (Ward-Thompson et al. 1994, 1999). Therefore, PPCs can be considered as the potential sites of future star formation and give us chances to probe the initial conditions of star formation. The series of papers of Ward-Thompson et al. (1994, 1999, 2000, and 2001) and the paper of André et al. (1996) argued that the density of PPCs is relatively constant in the inner cores and falls off at larger radii. This is very different from the singular isothermal sphere that leads to the “inside-out” collapse (Shu 1977). However, they assumed a uniform temperature in PPCs to calculate density structures. Evans et al. (2001) and Zucconi et al. (2001) show that the temperature of a dust grain decreases toward the centers of PPCs, assuming that dust grains are heated only by the attenuated external radiation field and that the radiation from the core itself is optically thin. Based on their results, Evans et al. (2001) argue that the continuum emission is significantly dependent on the temperature structure as well as the density structure, and they show that Bonnor-Ebert spheres can fit well the submillimeter dust continuum emission in PPCs. In agreement, the recent paper of Ward-Thompson et al. (2002) has shown that the color temperatures of PPCs that are calculated from the 170 and 200 $`\mu \mathrm{m}`$ dust emission observed by ISOPHOT decrease toward their centers. As a result of the lower temperature at the center, the density at the center must be more peaked than is the intensity of dust continuum emission. Evans et al. (2001) showed that even a power law could fit the density structure of L1544.

Ciolek & Mouschovias (1994) modeled ambipolar diffusion for the quasi-static evolution of a uniformly magnetized isothermal molecular cloud, finding relatively constant density toward the center. However, there has been debate about the importance of ambipolar diffusion because of differences between inward velocity observed (Tafalla et al. 1998, Williams et al. 1999) and calculated from the magnetic field (Goodman and Heiles 1994, André et al. 1996, Ward-Thompson et al. 1999, Crutcher and Troland 2000, Ward-Thompson et al. 2000). Those differences might be resolved by calculating more accurate density structure and, in turn, the correct temperature structure from the observed data. Ciolek and Basu (2000) showed if the density at the center of a PPC were higher, then their model could explain the observed magnetic field and inward velocity.

In addition to density and temperature distributions, we need to understand the ionization structure because magnetic relaxation by ambipolar diffusion is tightly related to the ionization structure (Rawlings 2000). Ionization fraction can affect the process of gravitational collapse because the coupling between ions and magnetic field plays a key role in preventing free-fall collapse of a dense core and slowing the drift of neutral gas through the magnetic field (ambipolar diffusion). In addition, the gas-phase chemical reactions depend significantly on ionization fraction, and the ionization fraction is highly sensitive to the depletion of molecules because some molecules play a crucial role in destroying ions. The degree of ionization in cold dense cores can be measured by the $`[\mathrm{DCO}^+]/[\mathrm{HCO}^+]`$ abundance ratio (Watson 1977), which is directly related to the electron density and depletion fraction. Recently, many molecular line observations show evidence for chemical differentiation and the depletion of some molecules such as $`\mathrm{C}^{17}\mathrm{O}`$, $`\mathrm{C}^{18}\mathrm{O}`$, and CCS in dense cores (Lada et al. 1994, Kuiper et al. 1996, Black et al. 1995, McMullin et al. 1994, Benson and Myers 1998, Willacy et al. 1998, Kramer et al. 1999, Ohashi et al. 1999, Caselli et al. 1999, Jessop and Ward-Thompson 2001, Tafalla et al. 2002). They suggested that this trend is due to the substantial depletion of molecules onto dust grain surfaces. In addition, several groups (Bergin and Langer 1997, Aikawa et al. 2001, Caselli et al. 2002, Li et al. 2002) modeled chemical evolution, including gas-dust interaction as well as gas-gas interaction, to show the changing depletions of some molecules and ionization fraction with time or density, and most of them compared their results to real observations. Therefore, both gas phase chemical processes and dust-gas interactions should be considered in studying star formation.

Those theoretical and observational results suggest that the relative distributions of several molecules can be used to probe the timescale of the dynamical evolution of a PPC. However, Rawlings et al. (1992) and Rawlings and Yates (2001) modeled several molecular lines combining self-consistently dynamical, chemical, and radiative transfer models, and showed that the abundance variation and the line profiles are very sensitive to the free parameters in the chemical model. Consequently, we need to use multiple transitions and maps in order to constrain these models.

In this paper, we explore the chemical and dynamical differences in three PPCs (L1512, L1544, and L1689B) comparing molecular depletion and ionization fraction by simple analysis of integrated intensity and column density and by modeling several molecular line profiles in detail. In this study, we cover a very early-time molecule (CCS), an early-time molecule (CO), a middle-time molecule ($`\mathrm{HCO}^+`$), and a very late-time molecule ($`\mathrm{N}_2\mathrm{H}^+`$) in chemical evolution. According to the chemical models of Bergin and Langer (1997) and Aikawa et al. (2001) that used an initial density of $`10^4`$ cm-3, the CCS abundance decreases rapidly after a few times $`10^5`$ years since CCS is a polar molecule, and it is easily frozen out onto the surface of dust grains. In addition, CCS is destroyed by HCO+ or H3O+, whose abundances increase until about $`10^6`$ years. Even though CO is a non-polar molecule, it has very high binding energy (1740 K onto water ice mantle grains, Bergin and Langer 1997). As a result, CO starts to be depleted onto the surface of dust grains after about $`10^6`$ years causing the depletion of HCO+ in turn. HCO+ is formed by the interaction between H$`{}_{}{}^{+}{}_{3}{}^{}`$ and CO or H2D+ and CO. However, N2, which is the precursor of N2H+, starts to be abundant much later than other molecules and reaches its maximum around $`10^7`$ years. N2 has low binding energy (750 K, Aikawa et al. 2001) onto the surface of dust grains, so it is not frozen out from the gas phase easily. In addition, H$`{}_{}{}^{+}{}_{3}{}^{}`$, which is the other precursor of N2H+, increases as CO is depleted. Another cause of the enhancement of N2H+ at late times is the decline of the destruction of N2H+ by CO, which is significantly depleted at late times.

We observed multiple positions in each core to study the chemical distribution within a core in addition to studying the chemical differences among the three PPCs. The most important aspect of this study is that we start with reasonable physical models calculated from dust analysis. The density and temperature structures found by Evans et al. (2001) are adopted for modeling molecular line profiles.

In section 2, we summarize observational parameters, and section 3 shows the results of observations. In section 4, we analyze column densities of several molecules. Section 5 compares modeled molecular line profiles and observed line profiles. Finally, we discuss our results regarding depletion and ionization in the view of chemical and dynamical evolution in section 6.

## 2 OBSERVATION

We observed the three pre-protostellar cores in $`\mathrm{C}^{18}\mathrm{O}`$ ($`\mathrm{J}=21`$ and $`\mathrm{J}=32`$), $`\mathrm{C}^{17}\mathrm{O}`$ ($`\mathrm{J}=21`$), $`\mathrm{DCO}^+`$($`\mathrm{J}=32`$), $`\mathrm{HCO}^+`$ ($`\mathrm{J}=32`$), and $`\mathrm{H}^{13}\mathrm{CO}^+`$ ($`\mathrm{J}=32`$) with the 10.4 m telescope of the Caltech Submillimeter Observatory (CSO) at Mauna Kea, Hawaii from 1995 to 2002. We used an SIS receiver with an acousto-optic spectrometer (AOS) with a 50 MHz bandwidth and 1024 channels, and the frequency resolution was about 2 to 2.5 channels, 0.13 to 0.16 $`\mathrm{km}\mathrm{s}^1`$ at around 220 GHz. The pointing uncertainty was about 4$`\mathrm{}`$ on average. We also observed these cores in $`\mathrm{H}^{13}\mathrm{CO}^+`$ ($`\mathrm{J}=10`$), $`\mathrm{N}_2\mathrm{H}^+`$ ($`\mathrm{J}=10`$), and CCS ($`\mathrm{N}_\mathrm{J}=4_33_2`$) with the 45 m telescope of the Nobeyama Radio Observatory in Japan on January 2002. We used the S40, S80, and S100 receivers with an AOS with 40 MHz bandwidth and 2048 channels. The frequency resolution was about 2 channels (about 40 kHz) so that the velocity resolution was about 0.27, 0.14, and 0.13 km s-1 at 43, 86, and 100 GHz, respectively. The observed sources are listed in Table 1; the observed lines and the main beam efficiency ($`\eta _{mb}`$) of each observation are summarized in Table 2. In Table 1, the last column shows the offsets of the dust peaks from the given coordinates; all offsets in figures other than Figure 2 and 3 are relative to the dust peak, assumed to be the center of the core in all modeling.

## 3 RESULTS

### 3.1 Spectra

Figure 1 shows the spectra of molecular lines observed toward the centers of the three PPCs. The centers are the positions of the centroid of the dust emission (Shirley et al. 2000). Most of the molecular lines have not been mapped completely but observed through two orthogonal cuts in $`\alpha \delta `$, crossing at the centers of dust emission in the cores. We do not see any consistent trend in the intensity of each spectrum from core to core. For example, $`\mathrm{C}^{18}\mathrm{O}`$ and $`\mathrm{C}^{17}\mathrm{O}`$ lines are the strongest in L1689B, but $`\mathrm{N}_2\mathrm{H}^+`$ and CCS are the strongest in L1544. All the lines except for $`\mathrm{N}_2\mathrm{H}^+`$ and CCS are much weaker in L1512 than in L1689B. The line width (FWHM) of the $`\mathrm{C}^{18}\mathrm{O}`$ $`\mathrm{J}=21`$ line is about 0.3, 0.5, and 0.6 km s-1 in L1512, L1544, and L1689B, respectively.

### 3.2 Maps

Figure 2 shows maps of three PPCs, where the grey scale indicates 850 $`\mu \mathrm{m}`$ dust continuum intensity and the contours are for integrated intensities of $`\mathrm{C}^{18}\mathrm{O}`$ $`\mathrm{J}=21`$ and $`\mathrm{DCO}^+`$ $`\mathrm{J}=32`$ lines. We can see that the peaks of dust emission and $`\mathrm{C}^{18}\mathrm{O}`$ $`\mathrm{J}=21`$ emission are not consistent. In the map for L1512, especially, the centroid of the dust emission coincides with a hole in the $`\mathrm{C}^{18}\mathrm{O}`$ $`\mathrm{J}=21`$ emission. In contrast, the distribution of $`\mathrm{DCO}^+`$ $`\mathrm{J}=32`$ emission appears to peak toward the center of dust emission in L1544 and L1689B, even though our maps of $`\mathrm{DCO}^+`$ $`\mathrm{J}=32`$ emission are not complete.

## 4 SIMPLE ANALYSIS

In this section, we use a simple analysis with standard techniques for comparison to similar work done by others, and, in the next section, we will deal with detailed models. We use dust emission as the tracer of column density least affected by chemistry in this section, and we assume that the dust emissivities do not change with radius. If temperature decreases toward the center, the column density calculated from dust emission with the assumption of an isothermal case could be underestimated toward the center. Even though the absolute column density calculated from the dust emission can change with the assumed dust opacity by a factor of 2 or 3, its relative distribution with radius can be traced relatively well by dust emission.

### 4.1 Integrated Intensities

We can compare the distribution of a molecular component with that of the dust component within a core using the integrated intensity of an optically thin molecular line and the dust continuum intensity observed at a submillimeter wavelength. In the upper panels of Figure 3, we plot the logarithmic integrated intensities of CCS $`\mathrm{N}_\mathrm{J}=4_33_2`$, H13CO+ $`\mathrm{J}=10`$, and $`\mathrm{N}_2\mathrm{H}^+`$ $`\mathrm{J}=10`$ lines in order to compare to the 850 $`\mu \mathrm{m}`$ dust continuum intensity through the cuts with constant declination shown in Figure 2. The dust continuum intensity and the H13CO+ intensity have been shifted upward to compare the shapes of intensity distributions with radius more effectively, and the intensity of dust continuum has been measured into apertures of 20$`\mathrm{}`$, spaced by 20$`\mathrm{}`$ for consistency with the H13CO+ and $`\mathrm{N}_2\mathrm{H}^+`$ observations.

The dust emission of L1512 is about half of that of L1544, and the distribution is not as peaked as that of L1544 and L1689B. L1544 and L1689B have similar dust emission even though it is more peaked in L1544. However, the CCS intensity of L1689B is more similar to that of L1512 rather than to that of L1544, and the CCS intensity of L1544 is about twice those in L1512 and L1689B. The CCS intensity in all the three PPCs is not peaked as much toward the center as the dust emission. However, it has a bump toward the centers of L1512 and L1544 (see §6.1). In L1689B, the bump is not certain because of the low signal-to-noise ratio.

If we check the relative variations of the integrated intensity of the $`\mathrm{H}^{13}\mathrm{CO}^+`$ $`\mathrm{J}=10`$ line with radius in each core, in L1689B, it is skewed to the west but follows well the dust continuum emission, decreasing at the edge of the core. On the other hand, the distribution of the integrated intensity of the $`\mathrm{H}^{13}\mathrm{CO}^+`$ $`\mathrm{J}=10`$ line in L1512 is not well matched with dust continuum, and the integrated intensity is even depressed toward the center. We also see the depression of $`\mathrm{H}^{13}\mathrm{CO}^+`$ emission at the center of L1544. This depression was predicted by Caselli et al. (2002) in their chemical models, but not observed, probably because their angular resolution (28$`\mathrm{}`$) was worse than ours (18$`\mathrm{}`$). We suspect that the optical depth of the $`\mathrm{H}^{13}\mathrm{CO}^+`$ $`\mathrm{J}=10`$ line in L1544 and L1689B is not negligible. Therefore, the depression in those cores could be partially caused by optical depth. However, L1689B does not show a depression even though it has a density structure similar to L1544 suggesting that depletion is playing a major role in L1544.

The $`\mathrm{N}_2\mathrm{H}^+`$ intensity is well peaked in L1544 and shows almost the same shape as the dust emission. L1689B also shows similar shapes in dust emission and $`\mathrm{N}_2\mathrm{H}^+`$ intensity. However, the $`\mathrm{N}_2\mathrm{H}^+`$ intensity of L1689B is only two thirds of that of L1544, even though these two cores have very similar dust emission and density structures, with central densities of $`10^6`$ cm-3. In L1689B, the peak of the intensity of $`\mathrm{N}_2\mathrm{H}^+`$ is shifted from that of the intensity of the dust emission similarly to the intensity of $`\mathrm{H}^{13}\mathrm{CO}^+`$. However, the ratio of the intensity of $`\mathrm{H}^{13}\mathrm{CO}^+`$ to the intensity of $`\mathrm{N}_2\mathrm{H}^+`$ at the center of L1689B, is greater than that in the other two cores.

In L1512, with a central density of $`10^5`$ cm-3, $`\mathrm{N}_2\mathrm{H}^+`$ emission is more extended than the dust, unlike in the other cores. Since the optical depth of the $`\mathrm{N}_2\mathrm{H}^+`$ $`\mathrm{J}=10`$ line in L1512 is smaller than in L1544 (Benson et al. 1998), this might indicate the possible presence of depletion at the dust peak.

In the following section, we compare the distributions of CO isotopes to that of dust emission by using molecular hydrogen column densities calculated from them instead of using the intensity distribution because the optical depth effect in the C17O line can be calculated, and the standard abundance of CO is known.

### 4.2 Column Density

We compare the $`\mathrm{H}_2`$ column densities calculated from observed molecular line intensities with that calculated from submillimeter continuum emission. This simple analysis assumes that 1) dust continuum emission traces all material along a line of sight in a core, and dust and gas are well mixed, 2) the core is isothermal, and dust and gas are well coupled so that they have same temperature, 3) molecular abundances are constant through the core, and 4) all levels are in LTE. Therefore, the excitation temperatures in all levels are the same as the kinetic temperature and the same as the dust temperature.

#### 4.2.1 The $`\mathrm{H}_2`$ column density from molecular emission lines

If we assume a line is very optically thin so that absorption can be completely neglected, the equation of radiative transfer is given by

$$\frac{dI_\nu }{ds}=\frac{h\nu _{ul}A_{ul}}{4\pi }n_u,$$
(1)
where $`n_u`$ is the number density of upper level, and $`A_{ul}`$ is the Einstein coefficient. For the transition $`\mathrm{J}\mathrm{J}1`$ in a linear molecule, the resulting relation between $`N(x)`$, the column density of molecule $`x`$, and the integrated intensity of the line, is

$$N(x)=\frac{3kQe^{E_u/kT_{ex}}}{8\pi ^3\nu \mu ^2J}T_R𝑑v,$$
(2)
where $`E_u=hBJ(J+1)`$, $`B`$ is the rotational constant, $`T_{ex}`$ is the excitation temperature above the ground state, $`\mu `$ is the dipole moment, and $`Q`$ is the partition function. Equation 2 is valid only in the limit that $`\tau 0`$. For finite optical depth, an optical depth correction can be applied:

$$N_{thick}=N_{thin}\frac{\tau _\nu }{1e^{\tau _\nu }}.$$
(3)
In principle, one should also correct for the presence of the cosmic background radiation, $`T_{cmb}`$, but for the lines considered here, at $`\lambda <1.4`$ mm, the Planck function in temperature units is less than 0.23 K, completely negligible compared to likely values of $`T_{ex}`$. The effect on $`N(x)`$ of neglecting $`T_{cmb}`$ is less than 5% for $`T_{ex}10`$ K and $`\tau 1`$, typical of the conditions relevant here. These effects are less than any plausible calibration uncertainty.

The rotational constant and the dipole moment of each molecule used for this calculation are listed in Table 3. For the excitation temperature, we used the kinetic temperature of 10 K, a typical value in dark molecular clouds, under the assumption that all lines are thermalized, that is, $`T_i=T_{ex}=T_k=10`$ K. If we use 20 K for this calculation, the result is changed only by a factor of 1.2 for the $`\mathrm{J}=21`$ line of $`\mathrm{C}^{18}\mathrm{O}`$ and $`\mathrm{C}^{17}\mathrm{O}`$.

We calculated the column densities using $`\mathrm{C}^{18}\mathrm{O}`$ $`\mathrm{J}=21`$ and $`\mathrm{C}^{17}\mathrm{O}`$ $`\mathrm{J}=21`$ lines. First, the lines were assumed to be optically thin. However, we can calculate the optical depth of the $`\mathrm{C}^{17}\mathrm{O}`$ $`\mathrm{J}=21`$ line directly by comparing the relative intensities of its 9 hyperfine components (Ladd et al. 1998). The optical depth of the $`\mathrm{C}^{17}\mathrm{O}`$ $`\mathrm{J}=21`$ line was obtained by a fitting method in the program CLASS. The panels in the third row of Figure 1 shows the results of fitting the hyperfine components of the $`\mathrm{C}^{17}\mathrm{O}`$ $`\mathrm{J}=21`$ line in the center of each core. The total optical depths for L1544 and L1689B are not small enough to be ignored (Table 4), even though each individual hyperfine component is optically thin. In order to correct for this optical depth effect, we applied equation (3) to the above calculation for an optically thin line. For the $`\mathrm{C}^{17}\mathrm{O}`$ $`J=21`$ line, $`\tau _\nu `$ is the optical depth of the line center that is calculated by multiplying the relative strength (0.693) of the main group, which has four close components, by the total optical depth. The correction factors at the centers of L1544 and L1689B are about 1.4. Note that we do not use the excitation temperature of $`\mathrm{C}^{17}\mathrm{O}`$ $`\mathrm{J}=21`$ obtained by fitting the hyperfine structure to calculate the column density of $`\mathrm{C}^{17}\mathrm{O}`$ because this method assumes that $`T_{ex}`$ is constant along the line of sight. Our detailed models (see §5.3.1) shows that this assumption leads to incorrect results. Also we want to use the same method for $`\mathrm{C}^{18}\mathrm{O}`$ and $`\mathrm{C}^{17}\mathrm{O}`$ and compare the results.

Finally, we convert the column density of each molecule into the $`\mathrm{H}_2`$ column density using the abundance ($`X`$) of molecule $`x`$;

$$N(\mathrm{H}_2)=\frac{N(x)}{X(x)}\mathrm{cm}^2.$$
(4)
We used a CO abundance of $`2.7\times 10^4`$ (Lacy et al. 1994), about three times bigger than the abundance used in other studies such as Bacmann et al. (2002) and Jorgensen et al. (2002). The abundance ratios \[CO/$`\mathrm{C}^{18}\mathrm{O}`$\] and \[$`\mathrm{C}^{18}\mathrm{O}`$/$`\mathrm{C}^{17}\mathrm{O}`$\] are about 560 ($`\pm `$ 25) and 3.2 ($`\pm `$ 0.2), respectively (Wilson and Rood, 1994). The abundance of each molecule used in calculating the $`\mathrm{H}_2`$ column density is listed in Table 3.

#### 4.2.2 The $`\mathrm{H}_2`$ column density from dust continuum emission

If the emission from dust is optically thin, the observed flux density ($`S_\nu `$) can be related to the column density of gas by

$$N(\mathrm{H}_2)=\frac{S_\nu }{\mu m_H\kappa _\nu B_\nu (T_d)\mathrm{\Omega }},$$
(5)
where $`\mu `$ is the mean molecular weight, $`m_H`$ is the atomic mass unit, $`\kappa _\nu `$ is the mass opacity of dust per gram of gas, $`B_\nu `$ is the Planck function, and $`\mathrm{\Omega }`$ is the aperture solid angle. $`\mathrm{\Omega }=(\pi \theta ^2)/(4ln2)`$ for a circular Gaussian aperture.

We use column 5 of the Ossenkopf-Henning dust opacity table (Ossenkopf & Henning 1994), which represents agglomerated dust grains with thin ice mantles, and assume a gas-to-dust mass ratio of 100 to get $`\kappa _{850}=0.018\mathrm{cm}^2\mathrm{g}^1`$. We use the 850 µm data from Shirley et al. (2000) but measure the flux in an aperture of 33$`\mathrm{}`$, the size of main beam in the observations of C17O. If the dust temperature, $`T_d`$ is 10 K, the $`\mathrm{H}_2`$ column density can be calculated from

$$N(H_2)=3.36\times 10^{22}\times S_{850}(\mathrm{Jy})\mathrm{cm}^2.$$
(6)
If the temperature were 20 K, the calculated $`\mathrm{N}(\mathrm{H}_2)`$ would be lower by a factor of 3.3.

#### 4.2.3 Results

$`\mathrm{C}^{18}\mathrm{O}`$ $`\mathrm{J}=21`$ and $`\mathrm{C}^{17}\mathrm{O}`$ $`\mathrm{J}=21`$ lines have been used to calculate the $`\mathrm{H}_2`$ column density assuming constant molecular abundances. Basically, we assume the lines are optically thin. If a molecular line is optically thick, then the calculation with equation (2) and (4) underestimates the $`\mathrm{H}_2`$ column density. As mentioned in the previous section, however, we can calculate the optical depth of the $`\mathrm{C}^{17}\mathrm{O}`$ $`\mathrm{J}=21`$ line by fitting its hyperfine structure. We have corrected the column density of the $`\mathrm{C}^{17}\mathrm{O}`$ $`\mathrm{J}=21`$ line for optical depth in L1544 and L1689B even though the correction is not very big.

The results of fitting the hyperfine components of $`\mathrm{C}^{17}\mathrm{O}`$ $`\mathrm{J}=21`$ show that the total optical depth of this line is about 0.1 in L1512, but it is not negligible at the central positions in L1544 ($`\tau 1.0`$) and L1689B ($`\tau 0.8`$). In L1512, the error of the calculated optical depth is very big, because of the low S/N in the weak components. However, the ratio ($``$ 4.3) between the strongest component and the second strongest component in L1512 is bigger than that ($``$ 3.3) in the two other cores, so the total optical depth would be less than that of L1689B, at least. The abundance of $`\mathrm{C}^{18}\mathrm{O}`$ is about 3.2 times greater than that of $`\mathrm{C}^{17}\mathrm{O}`$. Therefore, the optical depth of $`\mathrm{C}^{18}\mathrm{O}`$ is about 3.2 greater than the total optical depth of $`\mathrm{C}^{17}\mathrm{O}`$, and the optical depth of $`\mathrm{C}^{18}\mathrm{O}`$ $`\mathrm{J}=21`$ must be considered in calculating column densities. We plot the $`\mathrm{H}_2`$ column density calculated from the $`\mathrm{C}^{17}\mathrm{O}`$ $`\mathrm{J}=21`$ line in the lower panes of Figure 3, where the optical depth effect has been corrected by equation (3). However, in the figure, we use the $`\mathrm{H}_2`$ column density calculated from the $`\mathrm{C}^{18}\mathrm{O}`$ $`\mathrm{J}=21`$ line based on the optically thin line approximation without correcting for its optical depth in order to show the significance of the optical depth effect. In fact, $`\mathrm{C}^{17}\mathrm{O}`$ has been observed in fewer positions so that we cannot calculate the optical depth of every $`\mathrm{C}^{18}\mathrm{O}`$ spectrum.

In the figure, we see significant differences between the $`\mathrm{H}_2`$ column densities calculated from the 850 $`\mu \mathrm{m}`$ emission and the $`\mathrm{C}^{18}\mathrm{O}`$ $`\mathrm{J}=21`$ line in all three PPCs. At central positions, the $`\mathrm{H}_2`$ column density calculated from the $`\mathrm{C}^{18}\mathrm{O}`$ $`\mathrm{J}=21`$ line is about 15 times smaller in L1512 and L1544, and 8 times smaller in L1689B than the $`\mathrm{H}_2`$ column densities calculated from the $`\mathrm{S}_{850}`$ $`\mu \mathrm{m}`$ emission. However, after correcting for the optical depth of the $`\mathrm{C}^{17}\mathrm{O}`$ $`\mathrm{J}=21`$ line, the difference between the $`\mathrm{H}_2`$ column densities calculated from $`\mathrm{S}_{850}`$ and from $`\mathrm{C}^{17}\mathrm{O}`$ has been reduced. In L1689B, especially, the difference is just a factor of 2 after the correction.

Another interesting thing is that the distribution of $`\mathrm{C}^{17}\mathrm{O}`$ in L1544 does not show the central depression observed in Caselli et al. (1999). This is probably due to the larger beam (33″) in this observation compared to that (21″) of Caselli et al. (1999), similarly to the case of H13CO+ (see §4.1). This result suggests that inadequate resolution can prevent us for calculating the actual depletion of molecules.

### 4.3 Conclusions

The integrated intensities and the $`\mathrm{H}_2`$ column densities calculated from several molecules have been compared with the intensity of 850 $`\mu \mathrm{m}`$ dust emission and the $`\mathrm{H}_2`$ column densities calculated from the dust emission in our simple analysis of molecular line observations. The first result of this analysis is that the dust emission is more peaked than the molecular emission. The depression of molecular emission toward the centers is partially caused by the optical depth effect on molecular lines. We could see this effect comparing the column densities calculated from $`\mathrm{C}^{17}\mathrm{O}`$ $`\mathrm{J}=21`$ (marginally optically thin and corrected by $`\tau `$) and $`\mathrm{C}^{18}\mathrm{O}`$ $`\mathrm{J}=21`$ (optically thick but not corrected by $`\tau `$) lines. The correction for the optical depth effect by the $`\mathrm{C}^{17}\mathrm{O}`$ $`\mathrm{J}=21`$ line removes some of the difference between the $`\mathrm{H}_2`$ column densities calculated from the $`\mathrm{C}^{18}\mathrm{O}`$ $`\mathrm{J}=21`$ line and dust emission. The remaining considerable discrepancy would be caused by the depletion of CO molecules. This result shows that we can overestimate the depletion of a molecule if we use an optically thick line without correcting for the optical depth. Jorgensen et al. (2002) have compared the H2 column densities calculated from the 1.3 mm dust continuum and the C17O 1$``$0 line and found similar depletion factors of CO to ours in L1544 and L1689B. We will explore the depletion in detail through modeling molecular line profiles using the Monte Carlo method in the next section. We emphasize that the amount of the discrepancy varies from core to core. After the correction for the optical depth effect, the $`\mathrm{H}_2`$ column density calculated from the $`\mathrm{C}^{17}\mathrm{O}`$ $`\mathrm{J}=21`$ line is 11 in L1512, 9 in L1544, and 2 in L1689B times smaller than the $`\mathrm{H}_2`$ column density calculated from dust emission (Table 5).

The simple method used in this analysis clearly reveals evidence for depletion, but it has several limitations for quantitative analysis. First, the beam size of each line observation is not the same. If we had mapped these cores completely in each line, then we would have been able to convolve the maps using the same beam size. Second, the variations of the excitation temperature and the abundance of a molecule along the line of sight are not included in this simple analysis. In addition, the absolute $`\mathrm{H}_2`$ column density calculated from dust emission may be uncertain by a factor of 2 to 3 because of the uncertainty of dust opacity. Due to these limitations, we have to model molecular lines and compare them to the observed lines in order to study depletion more in detail.

## 5 DETAILED MODELS

### 5.1 Method

We used the Monte Carlo (hereafter MC) method to calculate the radiative transfer of molecular lines (Bernes 1979, Choi et al. 1995). The MC code for this work has been developed by Choi et al. (1995). The MC code generates model photons at a random position, in a random direction, and at a random frequency with proper random number distributions. These photons go through a one-dimensional spherically symmetric molecular cloud adjusting the level populations of molecules. The MC code calculates the excitation by the model photons and uses statistical equilibrium to adjust each level population until the criteria of convergence are satisfied. The MC method is more powerful than the LVG method because it can deal with arbitrary distributions of systematic velocity, density, kinetic temperature, microturbulence, and abundance self-consistently. Therefore, we can combine reasonable physical models from analysis of dust emission with chemical models to produce molecular line simulations, and the depletion of molecules with radius can be tested.

For the MC simulation, we needed collision rates at lower gas temperatures than are usually available, so we extrapolated linearly to 5 K the collision rates for CO (Flower and Launay 1985), which we used for $`\mathrm{C}^{18}\mathrm{O}`$ and $`\mathrm{C}^{17}\mathrm{O}`$. Similarly, we extrapolated the downward collision rates for $`\mathrm{HCO}^+`$ (Flower 1999) to use for $`\mathrm{H}^{13}\mathrm{CO}^+`$ and $`\mathrm{DCO}^+`$. Upward rates were calculated for detailed balance for each isotope. The size of each core is fixed at 0.15 pc, and microturbulent velocity is assumed constant through each modeled core.

Once the MC code calculates each level population, we simulate specific molecular line profiles using the virtual telescope (VT) simulation to compare the modeled profiles with observed molecular line profiles. We can simulate a molecule with multiple transitions and multiple positions simultaneously in a specific core with the parameters such as beam sizes and main beam efficiencies that we used in the actual observations. The remaining lack of realism is caused by geometric differences between a 1-dimensional simulation and the actual core. As seen in Figure 2, the three cores are not perfectly spherically symmetric, but instead they are elongated. In order to compare simulated line profiles with real observations, we averaged the lines that have been observed at the same distance from the centers, e.g., the lines in the off-positions of $`\pm 30\mathrm{}`$ in $`\alpha `$ and $`\delta `$. However, some lines have not been averaged in all four directions because the observation did not cover all parts. We need to map cores fully and use a 2-dimensional code to do more accurate comparisons.

We model the transitions of four molecules, $`\mathrm{C}^{18}\mathrm{O}`$ ($`\mathrm{J}=21`$ and $`\mathrm{J}=32`$), $`\mathrm{DCO}^+`$ ($`\mathrm{J}=32`$), $`\mathrm{H}^{13}\mathrm{CO}^+`$ ($`\mathrm{J}=10`$ and $`\mathrm{J}=32`$), and $`\mathrm{HCO}^+`$ ($`\mathrm{J}=32`$). If the gas and dust of a core are well mixed and well coupled, we can apply the distributions of density and temperature of the dust to the gas. However, for a density lower than about $`10^4\mathrm{cm}^3`$, the gas temperature is not coupled well with the dust temperature. The excitation temperatures of the lines that have high critical density such as $`\mathrm{DCO}^+`$ and $`\mathrm{H}^{13}\mathrm{CO}^+`$ are mainly dependent on density so that the inner, denser region mainly contributes to the intensities of those molecular lines, and the outer, less dense region does not change our results much. However, we should in principle correct for the difference between gas and dust temperatures in the outer, less dense parts of cores, especially for the molecular lines such as $`\mathrm{C}^{18}\mathrm{O}`$ lines that have low critical density. We leave this correction for a future work.

### 5.2 Physical Models

We use two different types of physical models: Bonnor-Ebert spheres (Bonnor 1956, Ebert 1955); and Plummer-like models (Whitworth & Ward-Thompson 2001). Both these models allow an initial density distribution that is nearly constant for small radii, but approaches a power law at large radii. The Bonnor-Ebert sphere provides a good fit to the dust continuum emission of the three PPCs that we are modeling (Evans et al. 2001). The Plummer-like model can fit the dust continuum data for L1544 (Whitworth & Ward-Thompson 2001). It has the advantage that it provides a simple solution for the density and velocity field during a free-fall collapse from the initial state, thereby providing a velocity field, which we need for modeling some of our lines. The initial Plummer-like radial density profile (Whitworth and Ward-Thompson 2001) is

$$\rho (r,t=0)=\rho _{flat}\left[\frac{R_{flat}}{(R_{flat}^2+r^2)^{1/2}}\right]^\eta .$$
(7)
Even though this model is very simple, in that it assumes a pressure-free collapse, it may be able to explain the later phase of the evolution of PPCs. Li et al. (2002) show that the evolution of a core driven by ambipolar diffusion becomes closer to the evolution of a non-magnetic, free-falling core at higher density.

We use the Bonnor-Ebert sphere model that best fit the dust emission (Evans et al. 2001), and we assume a purely microturbulent velocity field, except for the $`\mathrm{HCO}^+`$ $`\mathrm{J}=32`$ lines in L1544 and L1689B and the $`\mathrm{H}^{13}\mathrm{CO}^+`$ lines in L1544. Because those line profiles show deep self-absorption features indicative of systematic radial velocities, we use a Plummer-like model for them. To constrain the various parameters in the Plummer-like model, we modeled the dust emission, as was done by Evans et al. (2001). The resulting values that fit best the data for both L1689B and L1544 are as follows: $`\eta =3`$, the initial $`\rho _{flat}`$ (t=0) of $`5\times 10^5\mathrm{cm}^3`$, and $`R_{flat}`$ of 2800 AU. Whitworth and Ward-Thompson found a steeper density structure ($`\eta =4`$) and a bigger $`R_{flat}`$ (5350 AU) because they used an isothermal sphere for their model. When the temperature decline toward the center is included, a smaller $`R_{flat}`$ is needed. The evolutionary timescale that fits the dust emission well is half the timescale at which the central point-mass is formed (the start of the Class 0 stage). The maximum velocity in this model is about 0.08 km/s, which is consistent with the results of Tafalla et al. (1998) and Williams et al. (1999) in L1544. The resulting density structure is compared to that of the best fitting Bonnor-Ebert sphere in Figure 4a; they are very similar. Likewise, the resulting dust temperature structure is very similar (Figure 4b). The velocity field from the Plummer-like model is shown in Figure 4c.

While the density and temperature structure of the best-fitting Plummer-like model is very similar to that of the best fitting Bonnor-Ebert sphere, the former does contain a systematic velocity field. To check whether that field affected the results for lines modeled with the Bonnor-Ebert sphere, we modeled all other lines using the Plummer-like model. The velocity structure does not affect the line profiles of $`\mathrm{C}^{18}\mathrm{O}`$, $`\mathrm{DCO}^+`$, and $`\mathrm{H}^{13}\mathrm{CO}^+`$; the two models have given the same results.

### 5.3 Results of Models

The main purpose of this modeling is to find how much a specific molecule is depleted, and how the depletion is distributed with radius by simulating molecular line profiles at multiple positions within a core and by comparing the simulated line profiles with observed line profiles. We tried several functional forms, such as exponential, power laws, Plummer-like, and step functions, for the distribution of a molecular abundance. In addition, we applied the results of Li et al.’s dynamical and chemical model (2002) to our modeling, but their model has too high a systematic velocity to fit the observed molecular line profiles. However, the group is improving their chemical and dynamical models while communicating with us. Figure 5 and Table 6 shows the comparison between results of different functional forms of depletion with radius in the C18O 2$``$1 line simulations. We used integrated temperatures to calculate $`\chi ^2`$ in each model. Among those distributions, a step functional distribution fits best the observed line profiles with radius. This seems reasonable because the results of chemical models show a sharp decrease in abundance with density or time. Here, we mainly show the results of the modeling with the abundance distribution of a step function that has three free parameters: the undepleted abundance of a molecule ($`X_0`$), the fractional depletion ($`f_D=X/X_0`$), and the radius inside which a molecule is depleted ($`r_D`$). Each can be different in different species. Therefore, the total depletion of a molecule is dependent on the combination of $`f_D`$ and $`r_D`$. We emphasize that the derived depletion factor from this modeling has an uncertainty of about a factor of 3 because of the uncertainty of the density structure calculated from the dust analysis.

We compared the distributions of observed and modeled integrated temperatures with radius to find the best fits. Figure 6.a shows the distribution of the reduced $`\chi ^2`$ of the models in the C18O 2$``$1 and 3$``$2 lines of L1512 in the space of $`f_D`$ and $`r_D`$. The depletion radius is well constrained around 0.075 pc, but the fractional depletion is constrained only to be greater than about 25. This trend occurs in every molecule that is significantly depleted, so we consider $`f_D`$ to be the lower limit of the depletion factor. In other words, the data are consistent with complete depletion inside some $`r_D`$. However, the distribution of $`\chi ^2`$ constrains only the upper limit of the depletion radius in H13CO+ and DCO+ of L1689B and in DCO+ of L1544, which show small $`r_D`$ and $`f_D`$ (see Table 8 and Table 9). For example, the $`\chi ^2`$ in H13CO+ of L1689B changes within a factor of about 2 (from 3 to 7) in all range of $`f_D`$ within $`r_D<0.011`$ (Figure 6.b). As we mentioned above, the density structure has uncertainty of a factor of about 3, so the calculated depletion in this model could have the same uncertainty.

#### 5.3.1 $`C^{18}O`$ $`J=21`$ and $`J=32`$

The best fits of $`\mathrm{C}^{18}\mathrm{O}`$ $`\mathrm{J}=21`$ and $`\mathrm{J}=32`$ line profiles in the three cores are shown in Figure 7. The parameters of the best fit models are summarized in Table 7. We used the standard $`\mathrm{C}^{18}\mathrm{O}`$ abundance of $`4.82\times 10^7`$ for the undepleted abundance ($`X_0`$) for all cores, which was used to calculate the column density of $`\mathrm{H}_2`$ ($`\mathrm{\S }`$ 4.2.1). According to the results, CO molecules are depleted within 0.075 pc (about 110″) by a factor of 25 in L1512, so CO is depleted everywhere in our C18O map. In the case of L1544, the depletion factor is 25 within 0.045 pc (about 70″). $`\mathrm{C}^{18}\mathrm{O}`$ $`\mathrm{J}=21`$ and $`\mathrm{J}=32`$ line profiles simulated simultaneously at different positions from the center fit the observed line profiles well in L1512 and L1544. On the other hand, we could not find a model to fit well all the observed line profiles from the center to the off-position of $`90\mathrm{}`$ in L1689B. The antenna temperatures of observed $`\mathrm{J}=21`$ and $`\mathrm{J}=32`$ lines are almost the same or increase outward from the center in L1689B, even though the outermost lines are not well sampled. This might indicate that the large molecular cloud surrounding this source (Loren 1989) contributes to the observed $`\mathrm{C}^{18}\mathrm{O}`$ lines. We tested a core that has a warm envelope surrounding an inner, denser core. The core has the same density and temperature structures inside 0.15 pc as those calculated from dust emission, and a constant density and temperature of $`10^3`$ cm-3 and 50 K from 0.15 to 0.6 pc, which are based on the 13CO observations (Loren 1989). The temperature of 50 K is the upper limit in the $`\rho `$ Oph molecular complex. The modeled line profiles are shown in Figure 7c, which fit better the observed line profiles than those in the previous model that does not have a warm envelope. In this model, the abundance distribution is the same as the previous one. In order to avoid confusion by the envelope, we need to use a more rare isotope or higher J transition to trace the high density core. Even though L1512 and L1544 are believed to be well isolated, they might also be surrounded by bigger clouds. The modeled results, however, fit well the observed $`\mathrm{C}^{18}\mathrm{O}`$ $`\mathrm{J}=21`$ to the off-position of 90″. The influence of bigger structures in L1512 and L1544 does not seem significant compared to L1689B.

Jessop and Ward-Thompson (2001) modeled the C18O J$`=`$2$``$1 and J$`=`$3$``$2 lines of L1689B using a density structure with the same central density ($`\rho _{flat}`$) and the same $`R_{flat}`$ as calculated by André et al. (1996) and a constant temperature. After comparing to observations they argued that the flat distribution of C18O is caused by the significant (about 95%) depletion of CO. They pointed out that the gas temperature must be above 14 K to have better agreement with observations. However, they compared the model to the observation within 32$`\mathrm{}`$, which cannot give a good constraint on the outer envelope. In addition, the lower density and higher temperature at the inner region than those assumed in this study make the transitions optically thinner than the actual lines, so the depletion factor to fit the observations could be overestimated. We modeled $`\mathrm{C}^{18}\mathrm{O}`$ lines using our density structure and a constant temperature of 14 K. The best-fit model has $`r_D0.035`$ and $`f_D10`$, which are bigger than our best-fit model, where the temperature structure is calculated from the analysis of dust emission. The modeled $`\mathrm{C}^{18}\mathrm{O}`$ J$`=`$2$``$1 lines within 60$`\mathrm{}`$ are stronger than the observed lines. However, the line at 90$`\mathrm{}`$ is still weaker than the observed line even though it is stronger than that in our best-fit model without a warm envelope.

#### 5.3.2 $`H^{13}CO^+`$ $`J=10`$ and $`J=32`$

Except for the C18O lines, $`X_0`$ is not fixed in modeling molecular lines, but it can be easily decided by the outermost line profile. The results of the best models in $`\mathrm{H}^{13}\mathrm{CO}^+`$ $`\mathrm{J}=10`$ and $`\mathrm{J}=32`$ are shown in Figure 8 and Table 8. We have observed the $`\mathrm{H}^{13}\mathrm{CO}^+`$ $`\mathrm{J}=10`$ line with a $`20\mathrm{}`$ grid out to the distance of $`60\mathrm{}`$ from the center with an $`18\mathrm{}`$ beam. The $`\mathrm{J}=32`$ line has been observed with $`30\mathrm{}`$ spacing with a $`26\mathrm{}`$ beam. In L1512, the $`\mathrm{H}^{13}\mathrm{CO}^+`$ $`\mathrm{J}=32`$ line, whose critical density is about $`10^6`$ cm-3, was not detected, so we simulated only the $`\mathrm{J}=10`$ line. The depletion factor in L1512 is 25 within 0.021 pc (38″). In L1689B, the best-fit model of $`\mathrm{J}=10`$ line also fits $`\mathrm{J}=32`$ line very well. The depletion factor and radius in L1689B are much smaller than those in L1512. However, the undepleted abundance of $`\mathrm{H}^{13}\mathrm{CO}^+`$ in L1512 is 1.4 times bigger than in L1689B.

The observed $`\mathrm{H}^{13}\mathrm{CO}^+`$ $`\mathrm{J}=10`$ line, whose critical density is about $`10^5`$ cm-3, shows a blue skewed profile in L1544. Contracting cores show stronger blue peaks and self-absorption dips in optically thick lines (Myers et al. 1996). Evidence of inward motion in L1544 has been observed by Tafalla et al. (1998) and Williams et al. (1999). Therefore, we used the best-fit Plummer-like model in the simulation of the $`\mathrm{H}^{13}\mathrm{CO}^+`$ $`\mathrm{J}=10`$ line. In L1544, models fit well the shapes of the observed $`\mathrm{H}^{13}\mathrm{CO}^+`$ $`\mathrm{J}=10`$ line profiles, but the modeled $`\mathrm{J}=32`$ line is stronger than the observed line by a factor of 2. The rms noise of this line is about 0.015 K (5$`\sigma `$ detection), and the baseline is poor, so the discrepancy may not be very significant. There are several possibilities. If there is molecular gas surrounding L1544, this less dense gas would increase the $`\mathrm{J}=10`$ line more than the $`\mathrm{J}=32`$ line. Therefore, our models would overestimate the abundance of the inner, denser core to fit the $`\mathrm{J}=10`$ line. As a result, the modeled $`\mathrm{J}=32`$ line can be stronger than the observed line. Another possible explanation is that the depletion fraction increases toward the denser, inner region. Since the $`\mathrm{J}=32`$ line traces denser gas than the $`\mathrm{J}=10`$ line, the $`\mathrm{J}=32`$ line in our models where we used a step function of the abundance would be stronger than the observed line if $`\mathrm{H}^{13}\mathrm{CO}^+`$ is depleted more at the center. The uncertainties in the observed line or in the density and temperature structures could explain the difference between the modeled and observed lines.

#### 5.3.3 $`DCO^+`$ $`J=32`$

The modeling of $`\mathrm{DCO}^+`$ $`\mathrm{J}=32`$ does not constrain free parameters well because the observations in the line do not cover positions far from the center. The results are shown in Figure 9 and in Table 9. In L1512, only the center position was observed, so we just used a constant abundance to fit the line profile. The abundance of $`\mathrm{DCO}^+`$ that fits the line profile well is $`2.8\times 10^{10}`$. However, the line profiles of two positions in L1544 and three positions in L1689B were compared with the results of models. The factor and the radius of the depletion of $`\mathrm{DCO}^+`$ molecule are much smaller than those of CO molecules calculated by simulating $`\mathrm{C}^{18}\mathrm{O}`$ $`\mathrm{J}=21`$ and $`\mathrm{J}=32`$ lines. This result is consistent with the results of the chemical model of Caselli et al. (2002). The radius of the depletion is less (20″) than the beam size (35″) in L1689B and similar to the beam size in L1544.

#### 5.3.4 $`HCO^+`$ $`J=32`$

Infall asymmetry in PPCs has been observed (Tafalla et al. 1998, Lee et al. 1999, Gregersen and Evans 2000) in optically thick lines. If the gas of the outermost envelope of a core that has inward motion has sufficient optical depth and has a lower excitation temperature than the gas of the inner envelope in a given molecular line, the molecular line shows infall asymmetry, in which the intensity ratio of blue peak to red peak increases with time (Gregersen et al. 1997). We model one of the optically thick lines, $`\mathrm{HCO}^+`$ $`\mathrm{J}=32`$ in the three PPCs and compare the results with the observational results of Gregersen and Evans (2000). In order to model an $`\mathrm{HCO}^+`$ $`\mathrm{J}=32`$ line that has infall asymmetry, we used the best-fit Plummer-like model. Figure 10 shows the difference in MC modeling of the $`\mathrm{HCO}^+`$ $`\mathrm{J}=32`$ line for two different physical models: Bonnor-Ebert sphere and Plummer-like model. The latter includes an infall velocity structure and produces a blue asymmetry. The $`\mathrm{HCO}^+`$ $`\mathrm{J}=32`$ line in L1512 does not show a self-absorption dip and asymmetric profile so that we model the line in L1512 using the best-fit Bonnor-Ebert sphere.

The results of modeling $`\mathrm{HCO}^+`$ $`\mathrm{J}=32`$ in the three PPCs are shown in Figure 11, and in Table 10. Because L1512 has only one spectrum at the position 30$`\mathrm{}`$ away from the center, we used a constant abundance of $`\mathrm{HCO}^+`$ for this core. The abundance of the best fit is about $`1\times 10^9`$. In L1544, the intensity of the $`\mathrm{HCO}^+`$ $`\mathrm{J}=32`$ line at 30$`\mathrm{}`$ is half that at the center, unlike the $`\mathrm{H}^{13}\mathrm{CO}^+`$ 1-0 line, whose intensity is almost flat to 40$`\mathrm{}`$. As a result, an exponential abundance distribution as well as a step function also fits well the line profiles. In L1689B, the $`\mathrm{HCO}^+`$ $`\mathrm{J}=32`$ line has a stronger red-peak asymmetry at the center and 15$`\mathrm{}`$ away, but it has a stronger blue-peak asymmetry at 30$`\mathrm{}`$ and 45$`\mathrm{}`$. This complication of line profiles causes difficulty in fitting the exact line profiles in this core. Therefore, we tried to fit the integrated temperatures of the lines with radius. The exponential function of the abundance shows similar results to the step function. In order to test the functional form of the abundance, we need a more extended map. The interesting result is that the $`\mathrm{HCO}^+`$ $`\mathrm{J}=32`$ line of L1689B has a similar strength to that of L1544, unlike $`\mathrm{H}^{13}\mathrm{CO}^+`$ $`10`$, which is stronger in L1689B than in L1544. This might be due to the fact that L1689B is embedded in a larger molecular cloud that causes more absorption of the $`\mathrm{HCO}^+`$ 3$``$2 line.

On the whole, the abundance of $`\mathrm{HCO}^+`$ in this result is very small compared to the abundance that is expected from $`\mathrm{H}^{13}\mathrm{CO}^+`$ and the typical ratio of C and $`{}_{}{}^{13}\mathrm{C}`$ in the local interstellar medium of about 77 (Wilson and Rood 1994). However, our results show the abundance ratio of $`\mathrm{HCO}^+`$ and $`\mathrm{H}^{13}\mathrm{CO}^+`$ is approximately from 3 to 15 in the three PPCs. (We do not think that the ratio in L1512 is correct because we used only one spectrum in each line to calculate the abundance.) In addition, the depletion factors and the depletion radii calculated from the $`\mathrm{HCO}^+`$ line are different from those calculated from the $`\mathrm{H}^{13}\mathrm{CO}^+`$ lines. These differences and the low abundance ratio may result from the gas components outside the cores that we modeled in this study because the $`\mathrm{HCO}^+`$ line can trace less dense gas than $`\mathrm{H}^{13}\mathrm{CO}^+`$. As a result, the $`\mathrm{HCO}^+`$ line can be absorbed by the surrounding material significantly, but the $`\mathrm{H}^{13}\mathrm{CO}^+`$ line may be not absorbed as much as the $`\mathrm{HCO}^+`$ line. In our models, we limited the size of every core to 0.15 pc so that our models cannot account for the gas outside the core. Then, we would have underestimated the abundance trying to fit a line that has been weakened by the absorption by the surrounding gas component. However, this explanation is just one possibility. We cannot ignore the uncertainties of many other parameters, so we leave this as future work.

## 6 DISCUSSION

### 6.1 Depletion

We concluded that CO is significantly depleted in L1512 and L1544. It is difficult to estimate the depletion of CO in L1689B because of the gas surrounding this core and the optical depths of the lines used, but the depletion appears to be much less. The $`\mathrm{HCO}^+`$ molecule is depleted in L1512 and L1544, but we do not see significant depletion of this molecule in L1689B, based on modeling the $`\mathrm{H}^{13}\mathrm{CO}^+`$ $`\mathrm{J}=10`$ line. The optical depth and possible surrounding material also affect the calculation of the depletion of $`\mathrm{HCO}^+`$ in L1544. $`\mathrm{DCO}^+`$ is not significantly depleted in any of the PPCs. This result is consistent with results in many other starless cores where the distribution of DCO+ fits well the distribution of dust continuum (Myers 2002, personal communication). However, we need to use better resolution to do a more precise study because the calculated depletion radius is not bigger than the actual beam size (35$`\mathrm{}`$). Caselli et al. (2002) show that the depletion radius of DCO+ is about 3000 AU in L1544.

The analysis of the integrated intensity shows that the earliest-time molecule CCS is not peaked toward the centers of any of the PPCs, unlike the dust emission. Thus, CCS is depleted significantly in all three. As we saw, however, the CCS intensity in L1544 is twice stronger than that in L1689B even though the dust emission of L1689B and L1544 is similar. This result suggests that the depletion of CCS is connected to the depletion of CO. Li et al. (2002) modeled chemical abundance changes in PPCs using a dynamical model with a magnetic field, which drives the ambipolar diffusion, and they got a relatively high abundance of CCS. They indicated that this abundance could result from the interplay between significant depletion of CO molecules and late-time hydrocarbon chemistry. Ruffle et al. (1997, 1999) also showed the same result in their models of the time-dependence of several molecules. According to their models, the abundance of CCH and CCS have later secondary maxima only when the freeze-out timescale is long compared to the chemical timescale or the collapse timescale. CCS increases quickly to the first maximum of its abundance in time and is depleted very early as the density increases because CCS is destroyed by increasing HCO+ or H3O+, and this polar molecule can be tightly bound in grain mantles. However, in the region where CO molecules are depleted significantly, the $`\mathrm{C}^+`$ ion reacts more with molecules that are not oxygen-bearing species so that the precursors of CCS such as CCH increase again. This results in a second increase of the CCS abundance. The trend has not shown up in the models of Bergin & Langer (1997) and Aikawa et al. (2001) because the models used dense initial conditions in which a core evolves very quickly dynamically, and CCS did not show the second maxima. We think that the CCS molecule is possibly in the second enhancement in the late stage of chemical evolution resulting from the significant depletion of CO in L1512 and L1544. The integrated intensity of the CCS line shows bumps in the center regions of L1512 and L1544 (Figure 3). If CCS is depleted with time without the second enhancement, we should see holes rather than bumps. The other possible explanation of the enhancement of CCS is mixing of some atomic carbon from a translucent envelope into the core (van Dishoeck 2002, personal communication).

According to the comparison between the intensity of $`\mathrm{N}_2\mathrm{H}^+`$ lines and the dust emission intensity, this molecule is possibly depleted in L1512, but not in L1544 and L1689B, where its intensity follows well the dust emission intensity. $`\mathrm{N}_2\mathrm{H}^+`$ is known as a very late-time molecule, but Bergin et al (2002) have shown that $`\mathrm{N}_2\mathrm{H}^+`$ is depleted in a cold dark cloud, B68. The weak intensity of $`\mathrm{N}_2\mathrm{H}^+`$ line in L1689B, compared to in L1544, might indicate that this molecule has not yet reached its maximum because of less abundant precursors (N2 and H$`{}_{}{}^{+}{}_{3}{}^{}`$) but a more abundant destroyer (CO).

Bergin and Langer (1997) modeled chemical evolution in PPCs including gas-grain interactions as well as gas-phase reactions and compared two kinds of grain mantle properties: a weakly bound CO mantle and tightly bound $`\mathrm{H}_2\mathrm{O}`$ mantle. Their results show that sulfur-bearing molecules such as CS and SO are very sensitive to density and depleted seriously as the density of a core increases, regardless of the properties of dust grains. On the other hand, CO and $`\mathrm{HCO}^+`$ molecules are depleted in the model with $`\mathrm{H}_2\mathrm{O}`$ grain mantles but remain in gas phase in the model with CO grain mantles. In contrast, $`\mathrm{NH}_3`$ and $`\mathrm{N}_2\mathrm{H}^+`$ do not show depletion in any model because of the low binding energy of $`\mathrm{N}_2`$ to $`\mathrm{H}_2\mathrm{O}`$ and CO mantles. Therefore, the comparison of the results of their chemical models to our results indicates that the dust grains in the three PPCs have tightly bound $`\mathrm{H}_2\mathrm{O}`$ mantles. In addition, Aikawa et al. (2001) predicted the abundance distributions of molecules by numerical chemical models simply by applying a constant delay factor to a Larson-Penston model and considering different sets of adsorption energies, and they compared their results with the observed molecular abundances in L1544. According to their results, CCS and CO molecules are depleted, but $`\mathrm{N}_2\mathrm{H}^+`$ and $`\mathrm{NH}_3`$ molecules are more abundant in more slowly collapsing cores at the density peak.

To summarize our results on depletion, CCS, CO, and $`\mathrm{HCO}^+`$ are significantly depleted in L1512 and in L1544. On the other hand, the depletion of $`\mathrm{DCO}^+`$ in the two cores is not significant. $`\mathrm{N}_2\mathrm{H}^+`$ is not depleted in L1544, but it is possibly depleted in L1512. In contrast, only CCS is significantly depleted in L1689B. It is hard for us to estimate the depletion of CO in L1689B because of the optical depth effect and the material surrounding the cores.

### 6.2 Ionization

The coupling between ions and the magnetic field is crucial in regulating star formation rates if cores are magnetically subcritical. The timescale of ambipolar diffusion depends on the ionization fraction ($`\mathrm{t}_{\mathrm{AD}}2.5\times 10^{13}x(e)`$, Shu et al. 1987). In addition, the gas-phase chemical reactions are also dependent on the ionization fraction. Therefore, calculating the ionization fraction accurately is significant for understanding the chemistry as well as the timescale for ambipolar diffusion. However, it is not easy to calculate the ionization fraction in protostellar cores because it depends on the depletions of molecules as well as density structures and cosmic-ray ionization rate.

We used equation (16) in Caselli (2002) to calculate the ionization fractions:

$$x(e)=\frac{2.7\times 10^8}{[\mathrm{DCO}^+]/[\mathrm{HCO}^+]}1.9\times 10^7\left[\frac{1}{\mathrm{f}_\mathrm{D}}+\frac{\mathrm{x}(\mathrm{O})/\mathrm{f}_\mathrm{D}^{^{}}}{10^4}\right],$$
(8)
where $`\mathrm{f}_\mathrm{D}`$ is the depletion factor of CO, x(O) is the abundance of atomic oxygen ($`1.5\times 10^4`$, Caselli et al. 1998), $`\mathrm{f}_\mathrm{D}^{^{}}`$ is the depletion factor of O, and $`[\mathrm{DCO}^+]/[\mathrm{HCO}^+]`$ is the abundance ratio of two molecules, which can be calculated by comparing the column densities of $`\mathrm{DCO}^+`$ and $`\mathrm{H}^{13}\mathrm{CO}^+`$ and using 77 $`(\pm 7)`$ for the $`{}_{}{}^{12}\mathrm{C}/^{13}\mathrm{C}`$ ratio (Wilson & Rood 1994). Caselli et al. (1998) and Caselli (2002) derived this equation assuming that (i) the main source of $`\mathrm{HCO}^+`$ is the reaction between $`\mathrm{H}_3^+`$ and CO, (ii) the deuterium fractionation is due to the reaction between $`\mathrm{H}_3^+`$ and HD, and (iii) molecular ions are destroyed mainly by electrons and neutral species such as CO and O. In fact, Caselli et al. (2002) have showed that N2H+ and N2D+ are better tracers of ionization fraction because HCO+ and DCO+ are possibly depleted in the central region. However, we do not have N2D+ data, so we use $`\mathrm{DCO}^+`$ and $`\mathrm{HCO}^+`$, which can give correct information in the region with $`r>3000`$ AU.

We calculate, for each line of sight, the average $`x(e)`$ along the line of sight. This approach is not completely self-consistent because the method is only valid for homogeneous clouds. We use our model to simulate observation with the same beam size. The column density of each molecule has been calculated from the simulated line based on the assumption of an optical thin line. Figure 12 shows the ionization fraction with radius in the three PPCs.

There are some caveats in this calculation. First, we could not assess the depletion factor of CO in L1689B correctly because of the optical depth problem and the possible material surrounding the core ($`\mathrm{\S }`$ 5.2.1). Second, the model of $`\mathrm{H}^{13}\mathrm{CO}^+`$ lines in L1544 does not match the observations very well so that the calculated column density from the simulated $`\mathrm{H}^{13}\mathrm{CO}^+`$ $`\mathrm{J}=32`$ has more uncertainty. Third, since we do not have enough observed data to assess the depletion radius of $`\mathrm{DCO}^+`$ in L1512, the column density of $`\mathrm{DCO}^+`$ is uncertain in this core too. In addition to the caveats related to models, the equation that we used to calculate ionization fraction includes the abundance and the depletion fraction of O that are not well known. We supposed that the depletion of atomic oxygen might be proportional to the depletion factors of other molecules. We used $`\mathrm{f}_\mathrm{D}^{^{}}`$ of 3, 2, and 1 for L1512, L1544, and L1689B, respectively.

In spite of these caveats, the three cores have similar average ionization fractions (about $`5\times 10^7`$) along the line of sight toward the centers. These results are similar to the ionization fractions calculated in other protostellar cores (Caselli et al. 1998) but much bigger than the result (about $`10^9`$ at the center) of Caselli et al. (2002). Caselli et al. (1998) have used the integrated temperatures of the actual line profiles, but Caselli et al. (2002) have used the results of their “L1544 best-fit” chemical model, which takes into account the cloud density structure, molecular freeze out, and the recombination of molecular ions on grain surfaces, in order to infer the electron fraction as a function of cloud radius. Our data do not rule out still smaller values of $`x(e)`$ at small radii. The ionization fraction in L1512 and L1544 decreases toward the center, which trend has been shown in Caselli et al. (2002) and Caselli (2002). In contrast, in L1689B, it increases toward the center. As a result, the ionization fraction of L1689B around 0.06 pc is smaller than those of L1512 and L1544 by a factor of 10.

### 6.3 Timescales

The timescale of ambipolar diffusion ($`\mathrm{t}_{\mathrm{AD}}2.5\times 10^{13}x(e)`$, Shu et al. 1987) calculated from the above ionization fraction toward the centers is about $`1.25\times 10^7`$ years. This timescale applies to cores that are magnetically subcritical. This ambipolar timescale represents the duration of contraction from the present stage to the stage of Class 0. These are much bigger than the free fall timescale ($`10^6`$ years). However, $`x(e)`$ could be much smaller at the center. Even though these three cores have similar ambipolar diffusion timescale toward the centers, our results suggest that they have different chemical status from each other.

L1512 seems to have evolved chemically more than L1544 showing that CCS, CO, HCO+, and possiblely N2H+ are depleted. Caselli (2002) modeled a core that is less dense and less centrally concentrated than L1544 to show the chemistry in an early stage in the dynamical evolution. According to their result, molecules in the model core are not significantly depleted. The model core is similar to L1512 in density structure, but surprisingly, L1512 shows a bigger depletion radius of CO than L1544. In addition, N2H+ at the center of L1512 might start to be depleted. This result could result from the slower contraction of L1512 than L1544 so that L1512 could have a longer time to evolve chemically. Another possibility is a more efficient molecular freeze out in L1512 because the dust grains in this object are smaller than in the other two cores, so that the depletion timescale in L1512 is shorter (van Dishoeck et al. 1993; Caselli et al., in prep.). This could be possible because L1512 is less dense than the other two cores and has not have enough time for grain coagulation (Evans et al. 2001).

On the other hand, L1689B has not evolved much chemically even though the density structure is very similar to L1544 so that it is in a later stage of dynamical evolution than is L1512. This tells us that the chemical evolution is not a simple function of density, unlike the results of other chemical models. According to this result, we should not ignore the possibility that L1689B could be magnetically supercritical, and the core has evolved faster than the other two sources dynamically so that L1689B has not had enough time to evolve chemically. The broader lines in L1689B (Figure 1 and $`\mathrm{\S }`$ 3.1) indicate more active kinematics than the other two sources. The total masses of L1544 and L1689B are similar according to the dust modeling. Therefore, L1689B has to have a weaker magnetic field or a smaller ionization fraction in order to be magnetically supercritical. Crutcher & Troland (2000) measured the magnetic field along the line of sight in L1544 ($`B_{\mathrm{LOS}}11`$ $`\mu `$G) by using the Zeeman effect of the OH lines. There is no specific measurement of magnetic fields in L1689B. However, Troland et al. (1996) measured a OH Zeeman effect in the $`\rho `$ Oph cloud indicating $`B<10`$ $`\mu `$G. In addition, our calculation shows the ionization fraction in L1689B is smaller than in the other two cores in the outer part, where the calculation is more credible. The other possibility of the faster dynamical evolution of L1689B than L1544 is a larger external pressure (Galli et al. 2002) because L1689B is in the $`\rho `$ Oph complex (see §5.3.1).

Based on our results, we suggest that the stage of dynamical evolution of a core cannot be simply probed by its chemical status because the chemical evolution depends on the size of dust grains or the absolute timescale, during which the core has been in a given environment, as well as its density structure that shows the relative stage of dynamical evolution. We summarize the dynamical and chemical evolutionary stages of each core in Table 11.

## 7 SUMMARY

The summary of our results is as follows:

1) The difference between the column densities inferred from dust emission and the emission of a given molecular line can partially result from the optical depth of the line. However, the depletion of a molecule must explain the remaining difference. In L1689B, the main cause of the difference between the dust emission and the CO emission is the optical depth in the lines. However, in L1512, the depletion of CO is the major source of the difference between the dust emission and the CO emission.

2) In our Monte Carlo simulation of molecular lines, a step function with depletion in the inner region produces the best fit to the abundance variation of a molecule in PPCs.

3) The CO and $`\mathrm{HCO}^+`$ molecules are depleted significantly in L1512 and L1544, but not in L1689B. The gas around L1689B makes it difficult to calculate the actual depletion of CO.

4) The depletion of $`\mathrm{DCO}^+`$ is not substantial in the three PPCs.

5) The CCS molecule is substantially depleted in the three PPCs. However, CCS might be in the second enhancement in L1512 and L1544 because of the significant depletion of CO.

6) The distribution of the integrated intensity of $`\mathrm{N}_2\mathrm{H}^+`$ $`\mathrm{J}=10`$ suggests the possible depletion of $`\mathrm{N}_2\mathrm{H}^+`$ in L1512. On the other hand, $`\mathrm{N}_2\mathrm{H}^+`$ follows dust emission very well in L1544. In L1689B, $`\mathrm{N}_2\mathrm{H}^+`$ has not had enough time to reach its maximum.

7) The velocity structures calculated by a simple dynamical model that uses a Plummer-like density profile and a free-fall collapse model can produce the asymmetry of $`\mathrm{HCO}^+`$ $`\mathrm{J}=32`$ lines in L1544 and L1689B. However, the abundance calculated by modeling the lines is much smaller than what is expected from modeling $`\mathrm{H}^{13}\mathrm{CO}^+`$. The surrounding gas could possibly account for the low $`\mathrm{HCO}^+`$ abundance of the models.

8) The ionization fraction is similar toward the centers of the three cores even though L1512 is less centrally condensed than L1544 and L1689B. However, L1512 has a bigger depletion radius of CO than L1544 and possible depletion of N2H+ at the center. In contrast, L1689B does not show the maturity in chemical evolution compared to the other two cores. This difference suggests that the chemical evolution depends on the size of dust grains or the absolute timescale during which a core stays in a given environment as well as the density structure of the core that shows the relative dynamical evolutionary stage.

9) The evolved density structure and the young chemistry of L1689B suggest that this core is possibly experiencing a free-fall collapse rather than ambipolar diffusion. As a result, L1689B evolves dynamically too fast to evolve chemically. This would be possible if the magnetic field in L1689B is weak enough, or the external pressure is big enough to make L1689B magnetically supercritical.

## 8 Acknowledgments

We are grateful to Paola Caselli, the referee of this paper for helpful comments. We thank the staff of the Caltech Submillimeter Observatory and the Nobeyama Radio Observatory for assistance in using 10.4 m and 45 m telescopes. We also thank E. van Dishoeck, P. Myers, and J. Rawlings for their useful comments. We are very grateful to Z.-Y. Li and V.I. Shematovich for providing the data from their models. We thank State of the Texas and the NSF (Grant AST-9988230) for support. J.-E. Lee thanks the Rotary Foundation for support through an Ambassadorial Scholarship.


format: text

text:
# Ξ⁻ and Ξ̄⁺ baryon production in Au+Au collisions at √s__NN=130 GeV

## 1 Introduction

The study of multi-strange baryons in ultra-relativistic heavy ion collisions was first motivated by the prediction of an enhancement of strange $`q\overline{q}`$ pair production due to the formation of a quark gluon plasma (QGP) compared to that found from a hadron gas state . The enhancement is thought to be more pronounced for multi-strange baryons and anti-baryons since their production is strongly suppressed in hadronic interactions due to their high mass thresholds . Also, by comparing the multi-strange baryon yields and ratios to other particles with those predicted from particle production models, one can test whether the system has reached chemical equilibrium or not. The study of transverse momentum spectra will allow us to test the thermal equilibration of the system.

The Relativistic Heavy Ion Collider (RHIC) allows for the study of nuclear matter in extreme conditions such as in the case of the theoretical QGP. We present here preliminary measurements of multi-strange baryon production at mid-rapidity as a function of the collision centrality in $`\sqrt{\mathrm{s}_{_{\mathrm{NN}}}}=130`$ GeV Au+Au collisions from the STAR experiment.

## 2 Experimental setup and data analysis

The data presented here were taken with the STAR detector during the year 2000 Au+Au at $`\sqrt{\mathrm{s}_{_{\mathrm{NN}}}}=130`$ GeV run. The experimental setup for this period consisted mainly of a large cylindrical Time Projection Chamber (TPC) used for charged particle tracking and identification. A Central Trigger Barrel (CTB) and two Zero Degree Calorimeters (ZDC) were used for triggering. More details of the experimental setup can be found in . The data were divided into three centrality classes corresponding to $`010\%`$, $`1025\%`$ and $`2575\%`$ of the total hadronic cross section as described in .

The multi-strange $`\mathrm{\Xi }^{}`$ particles were reconstructed via their decay topologies, $`\mathrm{\Xi }^{}\mathrm{\Lambda }\pi ^{}`$ followed by $`\mathrm{\Lambda }p\pi ^{}`$ with a branching ratio of $`100\%`$ and $`64\%`$ respectively. The reconstruction of $`\mathrm{\Xi }^{}`$ particles was handled by first reconstructing $`\mathrm{\Lambda }`$ candidates as described in . These $`\mathrm{\Lambda }`$ candidates were then extrapolated backwards and combined with negatively charged tracks to determine $`\mathrm{\Xi }^{}`$ candidates, and similarly for the charge conjugate $`\overline{\mathrm{\Xi }}^+`$ decay. Particle identification from energy loss of charged tracks in the TPC as well as simple geometric and kinematic cuts were applied at both steps to reduce the large combinatorial background inherent to such a reconstruction process. For $`\mathrm{\Xi }^{}`$ and $`\overline{\mathrm{\Xi }}^+`$ the raw yields were extracted from the invariant mass distribution by counting the entries within $`\pm 15`$ MeV/c2 about the expected mass and then subtracting the background. The background under the peak was estimated by sampling two regions on either side.

## 3 Results and Discussion

The invariant mass distribution of $`\mathrm{\Xi }^{}`$($`\overline{\mathrm{\Xi }}^+`$) candidates was then histogramed in transverse mass $`m_T=\sqrt{p_T^2+m_0^2}`$ bins and the signal extracted in each bin as described above. Each $`m_T`$ bin was then corrected for detector acceptance and reconstruction efficiency by the Monte Carlo technique where simulated $`\mathrm{\Xi }^{}`$ particles were embedded into real events. The data cover $`|y|<0.75`$, simulations showed that the acceptance and efficiency for $`\mathrm{\Xi }`$ is constant, in $`y`$, over this range.

### 3.1 Centrality dependence of invariant yields and slope parameter

Each of the obtained preliminary $`m_T`$ spectra for $`\mathrm{\Xi }^{}`$ and $`\overline{\mathrm{\Xi }}^+`$ for the three centrality classes used in this analysis was then fit to an exponential function ($`\frac{dN}{dy}e^{(m_Tm_0)/T}`$) to determine simultaneously the inverse slope, $`T`$, and the rapidity density, $`\frac{dN}{dy}`$, integrated over all $`m_T`$. The results are summarized in table 1. We observe that the slope parameter is, within error bars, the same for $`\mathrm{\Xi }^{}`$ and $`\overline{\mathrm{\Xi }}^+`$ and is not strongly dependent on centrality. Also the $`\frac{dN}{dy}`$ for $`\mathrm{\Xi }^{}`$ and $`\overline{\mathrm{\Xi }}^+`$ is found to increase linearly with $`dN_h^{}/d\eta `$. Thus, no noticeable change of the production mechanism for $`\mathrm{\Xi }`$ baryons is observed in the covered centrality region. Systematic errors are estimated to be $`20\%`$. The feeddown from weak decay of $`\mathrm{\Omega }^{}`$($`\overline{\mathrm{\Omega }}^+`$) on $`\mathrm{\Xi }^{}`$($`\overline{\mathrm{\Xi }}^+`$) is estimated to be less than $`2\%`$ and thus it has been neglected for this analysis. The $`\mathrm{\Xi }^{}`$ and $`\overline{\mathrm{\Xi }}^+`$ particle yields in the measured $`m_T`$ region correspond to $`75\%`$ of the total yield.

### 3.2 Systematics of the inverse slope parameter

Figure 2 shows the mid-rapidity inverse slope parameters as a function of particle mass, for central Pb+Pb collisions at $`\sqrt{\mathrm{s}_{_{\mathrm{NN}}}}=`$17.2 GeV  and preliminary results from Au+Au collisions at $`\sqrt{\mathrm{s}_{_{\mathrm{NN}}}}=`$130 GeV . Note that both, statistical errors and systematic errors (horizontal lines) are reported on STAR points, while SPS points only include statistical errors. We observe the following : the inverse slope parameters are higher in RHIC collisions than in SPS collisions and the strange baryons seem to deviate from the collective flow pattern driven by lighter particles at both energies. The relatively low inverse slope parameter found for the $`\mathrm{\Xi }`$ baryons may be due to the functional form and the $`m_T`$ range used for the fit or it could indicate, as proposed in , that they freeze-out earlier during the system’s evolution.

### 3.3 Non-identical particle ratios

Figure 3 shows the mid-rapidity ratios (a) $`\mathrm{\Xi }/h^{}`$ and (b) $`\mathrm{\Xi }/\mathrm{\Lambda }`$ from heavy ion collisions at SPS , plus preliminary RHIC results. The $`\overline{\mathrm{\Xi }}^+/h^{}`$ ratio increases with beam energy indicating a higher exitation energy from the collision reached at RHIC energies compared to SPS energies. On the other hand, the $`\mathrm{\Xi }^{}/h^{}`$ ratio remains constant from SPS to RHIC energies due to the competing effect of the droping net-baryon density. Also from SPS to RHIC energies, the $`\mathrm{\Xi }^{}/\mathrm{\Lambda }`$ ratio increases due again to the drop of the net-baryon density which causes a smaller fractional rise on $`\mathrm{\Lambda }`$ than on $`\mathrm{\Xi }^{}`$, while the $`\overline{\mathrm{\Xi }}^+/\overline{\mathrm{\Lambda }}`$ is constant. The behaviour of $`\overline{\mathrm{\Xi }}^+/\overline{\mathrm{\Lambda }}`$ suggests a saturation of the strangeness production from SPS to RHIC. Note that $`\mathrm{\Xi }^{}/\mathrm{\Lambda }`$ and $`\overline{\mathrm{\Xi }}^+/\overline{\mathrm{\Lambda }}`$ ratios have been corrected for feeddown from weak decays of $`\mathrm{\Xi }`$ on $`\mathrm{\Lambda }`$. For SPS values the quoted  upper limits of $`5\%`$ for $`\mathrm{\Lambda }`$ and $`10\%`$ for $`\overline{\mathrm{\Lambda }}`$ were used.

Particle ratios between dissimilar species can be used to test particle production models. Using preliminary STAR data we have obtained $`\mathrm{\Xi }^{}/\pi ^{}=0.0088\pm 0.0004`$. ALCOR, a quark coalescence model, overestimates the $`\mathrm{\Xi }^{}/\pi ^{}`$ ratio by almost a factor of two . A statistical model assuming local chemical equilibrium and strangeness neutrality  is, within systematic errors, in agreement with the experimental value of $`\mathrm{\Xi }^{}/\pi ^{}`$ ratio. Finally the authors of  need to assume chemical non-equilibrium in their QGP model to better address the $`\mathrm{\Xi }^{}/\mathrm{\Lambda }`$ ratio.

## 4 Conclusion

In summary, the $`\mathrm{\Xi }^{}`$ and $`\overline{\mathrm{\Xi }}^+`$ invariant yields at mid rapidity increase linearly with the number of produced negative hadrons. The production of strangeness per produced hadron seems to saturate from SPS to RHIC energies, which is consistent with the system being in chemical equilibrium. The extracted inverse slope parameters are the same for both $`\mathrm{\Xi }^{}`$ and $`\overline{\mathrm{\Xi }}^+`$ and are not strongly dependent with centrality. For central collisions this inverse slope parameter is higher than that from collisions at the SPS energy and shows a deviation from a collective flow pattern that may indicate an earlier freeze-out of multi-strange baryons.


text:
# Novel Microwave Absorption Due to Strong Coupling between Josephson Plasma and the Josephson Vortex Array in Bi2Sr2CaCu2O8+δ

## Abstract

We have investigated the Josephson plasma excitations in magnetic fields parallel to the $`ab`$-plane in Bi2Sr2CaCu2O8+δ single crystals in a wide microwave frequency region (9.8 – 75 GHz). It was found that there are two kind of phase-collective modes: one increases with magnetic fields in higher fields, approaching linear asymptotic dependence with $`ck`$, which lies well above the inherent plasma frequency $`\omega _p`$, while the other does not show considerable field dependence. The higher linear mode is attributed to the Josephson plasma mode propagating along the reciprocal lattice vector of the Josephson vortex, whereas the lower one can be ascribed to the oscillation mode of Josephson vortices.

The interplay between the Josephson plasma (JP) and the Josephson vortex (JV), which is due to the linear and non-linear phase dynamics of the Josephson junction, has attracted much attention since the 1960s. In single junctions, the excitation spectrum has been investigated by Fetter and Stephen Fetter and Stephen (1968). They derived a gapless “vortex” mode attributed to the dc sliding of JVs and a “plasma” mode identical to the transverse JP, in which the propagation vector $`𝒌`$ is equal to the primitive reciprocal lattice vector of the one-dimensional JV array. The frequency of the traveling plasma wave can be controlled by the field parallel to the block (barrier) layer $`H_{}`$ according to the dispersion relation of the transverse plasma $`\omega _p^T=\omega _p\sqrt{1+(ck)^2}`$, where $`\omega _pc/\sqrt{ϵ}\lambda _c`$ is the inherent plasma frequency, with $`ϵ`$ and $`\lambda _c`$ being the high-frequency dielectric constant and the penetration depth along the $`c`$ axis. This plasma mode explains the so-called Ech resonance, which was found in the current–voltage ($`IV`$) characteristics of Pb/PbO/Pb junctions, with the resonance voltage being proportional to $`H_{}`$Eck et al. (1964).

In systems where the atomic-scale of weak Josephson junctions are stacking, namely an intrinsic Josephson junction (IJJ), which is realized in anisotropic high-$`T_c`$ superconductors represented by Bi2Sr2CaCu2O8+δ (BSCCO) Kleiner and Müller (1994), the physics attributed to JP, JV, and their interplay are drastically richerTachiki et al. (1994); Clem and Coffey (1990). Because of the larger charge screening length ($`\mu `$ 10 Å) compared with the thickness of the CuO2 bilayers ($`t=3`$ Å), Josephson coupling between not only adjacent layers but also next-adjacent layers contributes to the charge transport along the $`c`$ axis Machida et al. (1999). This longitudinal coupling yields the longitudinal JP mode propagating along the $`c`$ axis with a dispersion relation of $`\omega _p^L=\omega _p\sqrt{1+ϵ\mu ^2}`$ and rich variety of the dynamical nature of the two-dimensional JV lattice, which has been predicted by theoretical calculationsMachida et al. (2000) and partly found experimentallyHechtfischer et al. (1997).

So far, studies on the JP of high-$`T_c`$ superconductors in finite magnetic fields have concentrated on the case of the $`c`$ axis fieldShibauchi et al. (1999); Kakeya et al. (2001), where vortices perpendicular to the layers can be regarded as static matter to suppress the Josephson coupling, representing fluctuations of pancake vorticesKoshelev (1996). In contrast, in the field parallel to the $`ab`$ plane, vortices should be treated as dynamical matter because of the strong coupling between JP and JV. The microwave absorption associated with JP in the presence of JV is one of the key phenomena in revealing the linear and nonlinear dynamics of IJJ, which may also bring about high-frequency applications of high-$`T_c`$ materials.

In this letter, we describe the field dependence of the plasma frequency $`\omega _p(H_{})`$ in fields parallel to the $`ab`$-plane through the Josephson plasma resonance (JPR) measurements as functions of the microwave frequency, field intensity, and temperature, in order to reveal the entire nature of JP in the presence of JVs. We obtained two resonance modes at higher and lower frequencies in good accordance with recent theoretical results. The higher-frequency mode, which increases in proportion with the field, is attributed to JP mode propagating along the primitive reciprocal lattice vector of the JV lattice, while the lower frequency mode, which depends only weakly on the field is interpreted as the vortex oscillation mode modified by the vortex pinning and stacking effects of IJJ. This result explains the preliminary results of Matsuda et al., who have reported vanishing of the plasma resonance in fields very close to the $`ab`$ planeMatsuda et al. (1997).

We have measured JPR in three BSCCO single crystals grown by the modified traveling solvent floating zone method. Two of the crystals were under-doped (U1, U2), while the other was optimallydoped (OP), with all having been obtained after appropriate annealing processes in controlled atmospheres. The superconducting transition temperatures $`T_c`$’s of these crystals were determined by magnetization measurements with a SQUID magnetometer as 70.2, 72.5, and 90.5 K for U1, U2, and OP, respectively. All JPR data displayed here were obtained in U1. JPR measurements were made in a microwave frequency range between 9.8 and 75 GHz, with both reflective and transmission-type bridge balance circuits with rectangular TE102 mode cavity resonators. To selectively excite the longitudinal plasma mode only, the sample was placed inside the cavity so that the oscillating electric field parallel to the $`c`$-axis was exerted homogeneously over the $`ab`$-plane of the crystal Kakeya et al. (1998). Frequency-stabilized microwaves were generated by the signal swept generator (HP83650B) or Gunn oscillators, and the resonance was detected by changes in the cavity impedance (Q-factor) by sweeping temperatures. The magnetic field was applied by a split-pair superconducting magnet, and the angle between field direction and the CuO2 plane was adjusted by rotating the cavity resonator with respect to the magnet by a precision rotator. The alignment of the field direction exactly parallel to the $`ab`$-plane was determined from the symmetry of the angular dependence of the resonance within an accuracy of 0.01 degrees.

Figures 1 (a) shows the JPR absorption curves at 25.5 GHz obtained in various constant fields by sweeping temperature. Two clear resonance lines without hysteresis are found at higher and lower temperatures with a temperature gap in between. A sharp and symmetric resonance was observed just below $`T_c`$ in a low field region. With increasing fields, the higher temperature resonance mode shifts to lower temperatures below 1 kOe, then turns back to higher temperatures, and finally disappears above 1.5 kOe. In contrast, the resonance found at lower temperatures appears above 1.2 kOe. The resonance shifts to lower temperatures without a considerable change in its line-shape and intensity as the field is increased. The measurements with sweeping temperature at constant field were ideally reversible with well-defined low-resonance peaks at the higher temperatures and rather broad higher-field resonances at lower temperatures. We refer to the higher-temperature and lower-temperature resonances as HTB (higher temperature branch) and LTB peaks, respectively. That these two resonance branches are observed simulatanously and have different field dependences is peculiar to JPR in $`𝑯ab`$ and is in sharp contrast to JPR in $`𝑯c`$, where one resonance line whose resonance field monotonically decrease with temperature is observed in the configuration of the longitudinal excitationKakeya et al. (1998).

To study the doping dependence of these features, we performed similar measurements in crystals OP and U2. The JPR results at $`\omega /2\pi =`$ 42 GHz in OP were found to be qualitatively similar to the ones at 19 GHz in U1, although the resonance temperatures and fields were different. This similarity can be understood by considering the frequency ratio of $`\omega /\omega _p`$, which is less than one and closer to unity in more under-doped crystals for a given $`\omega `$ because $`\omega _p`$ is lower in higher anisotropic crystals. Provided that 19 GHz in U1 and 42 GHz in OP give the same $`\omega /\omega _p`$ value, we derive $`\omega _p`$ of OP as 130 GHz, which is a typical plasma frequency for optimallydoped BSCCO. Here, we used $`\omega _p`$ of U1 as 56.8 GHz, as described below. This scaling is also valid for results in U2, with slightly higher doping than U1. Therefore, the resonance splitting is observed at $`\omega \omega _p/3`$, and the magnitude of the temperature gap between LTB and HTB is larger at $`\omega `$, being closer to $`\omega _p`$. In contrast, $`\omega `$ is much smaller than $`\omega _p`$, LTB merges to HTB as observed at 9.8 GHz in U1, resulting in a single branch. This is the reason why Matsuda et al. did not observe vanishing of the resonance in an optimallydoped crystal in Ref Matsuda et al. (1997).

The resonance peaks of LTB and HTB as a function of magnetic fields for various frequencies are plotted in Fig. 2 (a) and (b), respectively. With increasing frequency, both LTB and HTB shift to lower temperatures. Because the decrease of LTB is larger than the one of HTB, LTB goes above our experimental temperature range above 30 GHz, whereas HTB could be observed throughout the temperature range below 60 GHz. The JPR can be clearly observed even in zero magnetic field, as seen in Fig. 1(a). This zero-field resonance occurs at a particular temperature $`T_0`$ and at a corresponding microwave frequency. The frequency dependence of $`T_0`$ has been reported previously and has been explained by using a simple two-fluid model with a temperature-independent scattering rate Kadowaki et al. (2000). By extrapolating the temperature dependence of the plasma frequency $`\omega _p(T)`$ to $`T=0`$, the inherent plasma frequency $`\omega _p`$ is estimated to be 56.8 GHz, which yeilds $`\lambda _c=217\mu `$m with $`ϵ=15`$.

The bending-over behavior of HTB becomes greater at higher frequencies, and the JPR above 30 GHz extends to temperatures even higher than $`T_0`$. This surprising behavior is certainly beyond our current understandings of JPR in perpendicular fields, where the plasma frequency should decrease both with application of high magnetic field and with increasing temperature because of the suppression of the interlayer coherence $`\mathrm{cos}\phi _{l,l+1}(H,T)`$ in response to the fluctuation of pancake vortices. Above 1.5 kOe, HTB obtained at higher frequencies lies at higher fields at the given temperatures, suggesting that the HTB resonance can be observed even above $`\omega _p`$. As shown in Fig. 1 (b), at 74.8 GHz where no JPR can be observed in the case of $`𝑯c`$, we detected a clear resonance. This resonance shows a remarkable angular variation and therefore can be found only within $`\pm `$ 3 degrees of the $`ab`$-plane above 3 kOe. The resonance field of this JPR mode shifts to high fields as the temperature is increased and finally disappears above approximately 5 kOe. A similar resonance was obtained at another three frequencies above $`\omega _p`$ (61.7, 65.9, and 74.3 GHz). These are plotted in Fig. 2 (b).

In Fig. 3(a), the frequency-field diagrams extracted from Fig. 2 are shown. It is clear that there are two resonance modes well-separated by a frequency gap. The higher one extracted from HTB shows an upturn at approximately 1 kOe, then monotonically rises to higher frequencies even above $`\omega _p`$, while the lower one extracted from LTB shows a gradual decrease with increasing field. Because LTB can be observed only in a finite field, the existence of JVs is considered to be required for LTB. At low fields below 0.5 kOe, the extrapolation of LTB to low fields tends to approach zero frequency, suggesting that LTB may be in a gapless mode at the zero-field limit. From these considerations we attribute LTB and HTB to the vortex and plasma branches as suggested in Ref. Fetter and Stephen (1968).

Recently, Bulaevskii et al. has proposed a qualitative picture of the plasma resonance in a layered superconductor on the basis of the single junction model Bulaevskii et al. (1997). With increasing parallel field $`H_{}`$, the frequency of the plasma resonance increases while the resonance intensity decreases because of the decrease in the optical weight of the plasma, and the resonance finally disappears in the high field limit $`H_{}H_0\mathrm{\Phi }_0/\gamma s^2`$, where $`\mathrm{\Phi }_0`$ is the flux quantum and $`s`$ is the interlayer spacing. In contrast, the vortex resonance survives even in $`H_{}H_0`$ and lies at a non-zero frequency by introducing periodic pinning of JVs. Thus, we attribute HTB and LTB to the plasma oscillation and the oscillation of JVs.

Quite recently, Koshelev and Machida have derived absorption spectra of the JP mode in a layered superconductor in high parallel fields, with a uniform oscillating electric field being applied parallel to the $`c`$-axis Koshelev and Machida (2001). Assuming that all block layers are fully occupied by JVs forming a distorted triangular lattice according to the anisotropy, they analytically and numerically calculated the dispersion relation and absorption spectra of the JP mode. Because the periodicity of the JV lattice induces a modulation of the interlayer phase difference, the $`𝒌`$-vector of the JP mode corresponds to the primitive reciprocal lattice vector of the JV lattice $`𝒒`$, which consists of the $`c`$-axis component $`q_c`$ fixed on $`\pi /s`$ due to the intrinsic pinning and the $`ab`$-plane component $`q_{ab}`$ proportional to the magnetic field $`H_{}`$ as $`2\pi sH_{}/\mathrm{\Phi }_0`$, as depicted in Fig. 3 (b). Therefore, a tilted JP wave as a mixture of the longitudinal plasma with $`k_z=q_c`$ and the transverse plasma with $`k_{xy}=q_{ab}`$ can be excited. The plasma frequency is determined by the dispersion relation along $`𝒒`$ lying between the longitudinal and transverse plasma modes and rises with magnetic field because of the linear increase in $`q_{ab}`$. It should be noted that the dispersion is very close to the longitudinal dispersion for $`q_{ab}q_c`$. As a result, the peak frequency of the dissipation spectrum, corresponding to the resonance peak, is proportional to $`H_{}`$ in high fields as

$$\omega _p(H_{})=\omega _p\frac{\pi H_{}\gamma s^2}{\mathrm{\Phi }_0}.$$
(1)
Assuming that $`\gamma =1070`$, a value that is quite reasonable for under-doped BSCCO samples, this equation provides a good agreement with the experimental data above 3 kOe at all temperatures, as shown in the inset of Fig. 3(a), where we used $`s=15`$ Å. It should be noted that this equation has been derived in Ref. Bulaevskii et al. (1997). In this situation, all block layers are expected to be occupied by JVs above $`H^{}=\mathrm{\Phi }_0/4.6\gamma s^2`$ Ichioka (1995), which is estimated to be 2.1 kOe for $`\gamma =1070`$. This value is consistent with the fact that the linear field dependence of HTB is violated below 3 kOe, as shown in Fig. 3(a). We therefore conclude that the linear increase in $`\omega _p(H_{})`$ of HTB is due to the increase in $`q_{ab}`$ of the JV lattice, while $`q_c`$ is fixed on $`\pi /s`$ because of the triangular lattice. This result suggests that the excited plasma frequency can be easily controlled by adjusting the magnetic fields parallel to the layers even above $`\omega _p`$.

The absorption intensity of the HTB resonance is thought to be diminished as frequency is increased because the intensity is weaker at higher temperature and higher frequency. This diminishing is predicted due to the increase in the quasiparticle damping along the $`c`$-axis, which is directly connected to the tunneling mechanism of high-$`T_c`$ superconductors. As discussed previously Tachiki et al. (1994), the higher quasiparticle conductivity along the $`c`$-axis broadens the resonance line and weakens the resonance intensity.

The origin of the intensity of the vortex resonance (LTB), which is comparable to the plasma resonance in our experiments, has not been conclusive, although a much weaker absorption below $`\omega _p`$ was obtained by numerical calculations in a single junction model in which the randomness of the critical Josephson current along the layer was introduced Kakeya et al. (2002). The frequency obtained in the calculations is much lower than the experimental results. This discrepancy can be attributed to the pancake pinning yielded by misalignments of magnetic fields and thermal fluctuations as predicted by Sonin Sonin (2001).

The initial decrease in HTB frequency with increasing field below 3 kOe cannot be explained by their approach, which may in part be due to the thermal fluctuation of JVs because the slope is larger at higher temperatures. Furthermore, it may be important to note that in the strongly coupled region JV lattice and JP soliton-like wave excitations as well as continuous waves are expected as described above. Further study is needed to obtain better understanding of these complicated dynamical behaviors.

In summary, we have for the first time found two microwave excitation modes with a temperature-dependent gap in parallel magnetic fields of BSCCO. Strong coupling between JP and the JV lattice are responsible for these two modes. We have concluded that in high fields where all block layers are occupied with JVs, one lying at higher frequencies comes from the JP wave propagating in a tilted direction with respect to the JV lattice structure, while the other seems to originate from the oscillation of the Josephson vortices.

###### Acknowledgements.

We thank M. Machida, A. Koshelev, T. Koyama, and M. Tachiki for their fruitful discussions and critical comments. This work has been supported by CREST, the Japan Science and Technology Cooperation.