text:
# Production of Super Heavy Elements at GANIL: present status and perspectives

## 1 Introduction

An experimental program dedicated to the study of Super Heavy Elements (SHE) started at GANIL in order to take advantages of the performances of the accelerator and of the powerful velocity filter LISE3.

Indeed because of the very low cross sections for the production of SHE (from 1 $`\mu `$barn for Z=100 to 1 pbarn for Z=112) , very high intensity ECR ion sources delivering about one particle$`\mu `$A (6.27x1012 projectiles/sec) or more should be used. The second important point is the separation between the projectiles and the complete fusion evaporation residues (ER) emitted at forward angle after the production target. Originally, the LISE3 Wien filter  was designed to improve the separation of exotic nuclei produced by fragmentation of the projectile and selected by the LISE spectrometer. Then, an important work has been done to adapt it to fusion experiments improving its ability to work with very high primary beam intensities.

An experimental setup composed of a dedicated rotating target, beam diagnostics, time detectors and a set of Si-detectors (for the implantation and decay products) has been developed. Finally a digital electronics is under development to reach extremely short decay times in the order of few $`\mu `$seconds.

The response of the whole setup was checked via fusion reactions with known cross sections and $`\alpha `$ decay chains: 204,206Fr87 nuclei formed via 86Kr+natSb with cross sections of 10 to 300 $`\mu `$barn . These equipments were then used for an experiment dedicated to the system 86Kr on 208Pb and the transmission and the rejection factor of the Wien filter were precisely measured in an experiment dedicated to the production of the Sg element (Z=106) through the system 54Cr+208Pb. Comparison with known values for the excitation energy function  allowed to check the excellent transmission together with an high rejection factor.

An experimental program has been established for the production of new elements and/or new isotopes with two principal axes: $`i`$ the use of the inverse kinematics with a beam of Lead and $`\mathrm{𝑖𝑖}`$ the use of the neutron-rich 76Ge as projectile. Another program dedicated to the study of the properties of SHE has also started recently.

## 2 Experimental setup

A general overview of the experimental setup, from the reaction chamber to the detection setup is presented in Figure 1.

### 2.1 Reaction chamber

The reaction chamber contains two 67 cm diameter wheels bearing 35 targets (5.0\*1.5 cm2 each) and rotating up to 2000 RPM, allowing targets with low melting temperature (like Pb or Bi) to sustain very intense primary beams. The second wheel support carbon stripper foils in order to re-equilibrate the charge state of evaporation residues. The time structure of the beam is synchronized with the rotation. Two smaller wheels allow us to run with targets having higher melting temperatures. A Si detector continuously monitors the status of each target via elastic scattering reactions.

### 2.2 Wien filter

The LISE3 Wien filter has crossed magnetic and electric fields in a rather compact geometry and its angular acceptance is $`\pm `$36$`\pi `$mrad. It is divided into two identical halves and is followed by a dipole magnet. Several modifications were found to be necessary for fusion experiments: -$`i`$\- the upper plate in the first half of the filter was moved up by 2 cm in order to deflect the beam without hitting the electrode. The beam is now stopped on a new water-cooled collimator, -$`\mathrm{𝑖𝑖}`$\- two pairs of independently movable slits and a beam profiler were installed at mid-filter. An additional background suppression can be performed using the last dipole magnet located between the Wien filter and the detection chamber. The transmission of fusion nuclei from the target position through the quadrupoles, filter and final dipole was studied via a simulation code . The optics are checked and the Wien filter is calibrated using low energy ions (0.35 MeV/u to have the same speed than ER’s after the target).

### 2.3 Detectors

The detection setup (see Figure 2) is composed of two Galottes (carbon foils + MCP’s time detectors) and 13 Si detectors. The time detectors are used in order to tag all the events (ER’s or background) coming from the beam. Since the efficiency of Galotte for light particles might not be 100%, two of them are used. Then, the ”implantation events” are those with an implantation energy signal in coincidence with both Galottes whereas the ”$`\alpha `$” or fission events are those in anti-coincidence. The implantation detector is a 300$`\mu `$m thick Si 48x48 strips. The distance between the strips is reduced ($``$45$`\mu `$m) to have a dead zone inferior to 5%. Two of them are mounted on a movable arm. When emission of $`\alpha `$-particles or fission fragments is identified via electronics, the beam is immediately stopped for a short while during which the first implantation detector is moved out of the beam and the second implantation detector comes in, then the beam is sent again. This technique allows us to wait for very long decay times in the first implantation detector, being completely free of background coming from the beam. The ER’s are implanted at a depth of approximately 10 $`\mu `$m which is smaller than the range in Si for $`\alpha `$ particles with energy above 2.5 MeV. Then, in order to measure the full energy of the $`\alpha `$ particles emitted at backward angles which escape the implantation detector, a set of 8 Si detectors is mounted in a ’tunnel’ geometry. With this configuration, the efficiency for full $`\alpha `$ energy goes up from 55 to 95%. A ’face’ Si detector is placed in front of the ’out of beam’ implantation detector in order to play the same role, the efficiency being then 100%. A veto Si detector is placed after the implantation detector to eliminate the background of light particles coming from the beam. All the Si detectors are mounted on copper supports and can be cooled to reduce the electronic noise.

### 2.4 Data acquisition and Electronics

Specific electronics and data acquisition systems were developed. We use for the data acquisition a complete VXI system for which the standard dead-time is around 150$`\mu `$s. In order to reduce it, a multi-trigger logic with 2 sets of converters (ADC’s) is working in the following way: the hard trigger provided by the electronics enters into two GMT (Ganil Master Trigger) modules, T1 and T2. In case of usual counting rate, all the triggers are treated as ’T1’ and the events are registered into the first set of converters. This way is exactly the same as the usual one used in many experiments. The same signals are also duplicated into the second set of converters but are not yet considered by the acquisition. In case of a new trigger comes during the dead time of GMT1, it is now treated by GMT2 (as a trigger ’T2’) while GMT1 continues to deal with the preceding event. The minimum time between two events is then about 30$`\mu `$s. This procedure was tested with a strong $`\alpha `$-source before the experiment. The used converters are new 14-bits ADC’s which allow to have a good resolution for $`\alpha `$-particles (less than 50 keV) together with a large range of energy, between 300 keV and 250 MeV (necessary for spontaneous fission events).

A fast analysis program allows us to identify ’on line’ $`\alpha `$-decay chains or spontaneous fission and then to stop the beam and move out the implantation detector (see sect. 2.3)

A complete read out digital electronics is actually under development to be able to reach extremely short decay times, in the order of few $`\mu `$sec keeping good energy resolution (less than 50 keV) on a large energy range from 1 MeV to 1 GeV. This new digital electronics will be placed after the analog front-end electronics (preamplifiers) and the main idea of such an approach is to adapt a dedicated processing to each event on a very large dynamical range and without limitation due to the recovery time of the preamplifier.

## 3 Experimental results

### 3.1 Experiment on 86Kr + 208Pb system

In July 1999, the BGS-Berkeley group reported the observation of 3 events of the new element Z=118 identified by $`\alpha `$-decay with chains of 6 $`\alpha `$’s . The reported cross section was much larger than expected ($``$2 pb) and the same system studied with the velocity filter SHIP at Darmstadt did not lead to the observation of this element , justifying to repeat this experiment at GANIL in december 99. Since GANIL is able to deliver an high intensity beam of 86Kr (15$`\mu `$A with charge 10+, i.e. nearly 1013projectiles/sec), our idea was to confirm (or not) Berkeley’s result and, in the case of a positive answer, to obtain additional information: -$`i`$\- be able to reach shorter decay times (by the use of fast electronics, see sect. 2.4.) for possible $`\alpha `$ decay before the minimum time of 120 $`\mu `$sec in the Berkeley experiment, -$`\mathrm{𝑖𝑖}`$\- wait for very long decay times (by the use of movable implantation detectors, see sect. 2.3.), -$`\mathrm{𝑖𝑖𝑖}`$\- since ER’s may be created in an isomeric state which decays via electron capture followed by an electron cascade with an unknown half-live modifying the ionic charge and strongly reducing the transmission in the case of the GSI experiment, the carbon stripper foils were in our experiment located at a distance from the target 3 times larger than in the SHIP experiment.

With a total accumulated dose of 1.1x1018 ions at 5.27 MeV/u on targets of 300 $`\mu `$g/cm2, we didn’t observe any event corresponding to the decay of element 118.

The result of Berkeley was recently retracted .

### 3.2 Experiment on 54Cr + 208Pb system

In order to measure precisely the transmission and the rejection factor of the LISE3 Wien filter for SHE’s, an experiment dedicated to the production of the Sg106 element through the system 54Cr on 208Pb was performed at GANIL in december 2000. With a primary beam intensity of 40 nAp and 2 incident energies (4.698 and 4.756 MeV/u), we observed 10 events of 261Sg106 (1 neutron evaporation channel) and 2 events of 260Sg106 (2 neutrons evaporation channel). The decay chains were identified by position correlation between the implantation of the ER and the decay events in the implantation detector. Comparison with the known excitation energy functions  makes possible to deduce a transmission efficiency above 60$`\%`$ together with a power of suppression of the primary beam and scattered projectiles of $``$2x1010, corresponding to a counting rate of 5-10 Hz in the implantation detector.

## 4 Short and Middle range Perspectives

### 4.1 Inverse kinematics

#### 4.1.1 Method

The use of the inverse kinematics (heavy projectile on light target) has been tried at GANIL in June 1999 and will be used in spring 2002 on the system 208Pb+54Cr to produce 261,260Sg106 isotopes . This technique presents several assets in comparison with direct kinematics

-$`i`$\- the target thickness is no longer limited by multiple scattering of the ER’s and it becomes possible to use thicker targets covering an excitation energy range of typically 10 MeV instead of 3-4 MeV. The counting rate is then maximized and the gain in time can be very important. Indeed, in the search for new isotopes or new elements, the incident energy depends strongly on the mass defects which differ in different mass tables by several MeV whereas the excitation function is rather narrow.

-$`\mathrm{𝑖𝑖}`$\- the ER’s are more strongly focussed at forward angles resulting in a better transmission through the LISE3 filter.

-$`\mathrm{𝑖𝑖𝑖}`$\- because of the higher energy after the target, the ionic charge distribution of ER’s is better estimated and its relative width is smaller. Here again the transmission is favoured.

-$`\mathrm{𝑖𝑣}`$\- finally, the quality of the data is better: due to the deeper implantation depth, the energy deposited by the ER’s is higher, both fission fragments are stopped and escaping $`\alpha `$ particles lose more than 2 MeV and then cannot be missed.

The main drawback of this approach concern the velocity difference between the beam and the ER’s which is much smaller than in direct kinematics. Then, higher electric and magnetic fields are necessary in the Wien filter. With the new power supplies installed this year, the LISE3 filter will be able to get the same spatial separation at mid-filter as in direct kinematics.

#### 4.1.2 Use of 208Pb beam

Several systems can be studied for different purposes (summarized in Table 1):

$`i\mathrm{𝑁𝑒𝑤}\mathrm{𝑝𝑎𝑡ℎ𝑠}\mathrm{𝑡𝑜}\mathrm{𝑘𝑛𝑜𝑤𝑛}\mathrm{𝑖𝑠𝑜𝑡𝑜𝑝𝑒𝑠}\mathrm{𝑤𝑖𝑡ℎ}{}_{}{}^{\mathit{208}}\mathrm{𝑃𝑏}:`$ Produce the same isotopes of the odd Z compound nuclei than ER’s formed at GSI with a 209Bi target . Here, two points are of interest: first, measure and then compare the cross sections which will give a clue to the respective roles of macroscopic and structure effects: charge asymmetry of the entrance channel and closed shells effects. Secondly, produce ER’s in different levels, especially isomeric states.

$`\mathrm{𝑖𝑖}\mathrm{𝑁𝑒𝑤}\mathrm{𝑖𝑠𝑜𝑡𝑜𝑝𝑒𝑠}\mathrm{𝑤𝑖𝑡ℎ}{}_{}{}^{\mathit{207}}\mathrm{𝑃𝑏}:`$ Produce isotopes lighter by 1 neutron than the ones known for several elements: 260Bh107, 265Mt109, 271111 and 276112. For the last two cases, a chain of 3 $`\alpha `$’s will be observed before reaching a known isotope, i.e. 3 new isotopes will be observed in one event.

$`\mathrm{𝑖𝑖𝑖}\mathrm{𝑁𝑒𝑤}\mathrm{𝑒𝑙𝑒𝑚𝑒𝑛𝑡𝑠}:`$ Above Z=109 the counting rates are very low. For 277112, measurements were made at one beam energy only, not necessarily at the maximum of the excitation function. The broader energy coverage provided by inverse kinematics may reveal a larger cross section. Moreover, the 1n channel cross section decreasing quickly with Z, it may be that the 0n cross section becomes larger even through the fusion probability drops very much at the corresponding incident energy. Added to the uncertainty on the compound nucleus mass defect, the exploration of a broad range of incident energies corresponding to 0n+1n channels leads to lengthy measurements at several incident energy in direct kinematics: 6 beam energies were tried for Se+Pb at GSI in the search of element 116  whereas an excitation energy range of 12 MeV or more can be covered with a single incident energy in inverse kinematics.

### 4.2 Use of 76Ge32 beam

The main interest of using 76Ge beam is the direct cold fusion production of the element 273110 using a 198Pt target . Indeed, this isotope has been observed in 3 events only through the decay of element 277112 and possibly 3 events by direct production in hot fusion reaction (34S+244Pu=273110+5n). A very large range of $`\alpha `$-decay energies, from 9.18 to 11.17 MeV with lifetimes from 0.11 to 170 msec have been observed. What we propose is to make a precise $`\alpha `$-spectroscopy of this isotope since the population of different levels (especially isomeric states) may be due to overreaching the deformed shell N=162. This would require a ”large” number of events.

A more general interest for the 76Ge beam is provided by the accessible region of production. Indeed, it makes possible, using targets ranging from 192Os76 to 205Tl81, to produce elements 108 to 113 heavier by two neutrons (except for element 111) than the isotopes produced using Pb/Bi targets with beams ranging from 48Ca to 70Zn . Clearly, comparing with Pb/Bi based reactions, the entrance channel is more symmetric and there are no closed shells, the cross sections will be then reduced. But on the other hand, the accessible compound nuclei are richer by 2 neutrons (except for element 111) and this will be in favour of the cross sections.

An estimation of the reduction factors (see Figure 3) can be made using the compound nucleus 220Th produced via different entrance channels having similar Bass barriers : Ar18+Hf72 and Zn30+Nd60 systems have no closed shells and the more symmetric system has a cross section smaller by a factor 14. Assuming a steady variation with the charge asymmetry, the expected cross sections for Ca20+Yt70 and Zr40+Sn50 are estimated. The actual cross sections are much higher, by a factor 6 and 23 respectively, due to the role of the closed shells . These factors may now be used to estimate the reduction factor which applies when the same CN is formed on 76Ge based reactions in comparison with Pb/Bi based reactions. The loss of shell effects is estimated to reduce the cross section by one order of magnitude (represented by a solid arrow in Figure 3). At Z=106, the decrease of asymmetry reduces the cross section by a factor about 12 (represented by a dashed arrow in Figure 3) whereas this asymmetry disappears for 76Ge+208Pb reaction producing element 114.

An estimation of the gain in the cross section due to the neutron richness can be made also considering the Th isotopes. Detailed measurements show that 2 more neutrons in the CN increase the residue cross section by a factor 9 . As seen in Figure 3, the gain for Z=110 using 64Ni beam instead of 62Ni is in the order of 5. In our estimation for 76Ge beam, we consider an average factor of 7 (represented by a solid arrow on Figure 3). Then, the estimated cross sections are represented by solid crosses. They are smaller than the Pb/Bi based cross sections.

### 4.3 Study of the structure of the Super Heavy Elements

A new experimental program is proposed at GANIL  on the study of transfermium isotopes by the use of spectroscopic information in order to have a better understanding of their shell structure. This program is going to start with the $`\alpha `$, $`\gamma `$ and electron spectroscopy of 251Md and 251Fm. 251Md is populated by the $`\alpha `$-decay of 255Lr produced in the 2 neutrons channel of the 48Ca+209Bi fusion evaporation reaction. The 251Md should decay with less than 10% by $`\alpha `$-decay on 247Es and with more than 90% by EC on 251Fm. Since the ground state parity and spin of the first level of 251Fm are known, ground state of 251Fm and possibly of 255Lr could be determined. The experimental setup will be adapted with a tunnel of Si detectors to detect conversion-electrons and a set of Ge clover detectors from Exogam will be added in a close geometry.

## Acknowledgements

The authors would like to acknowledge all the technicians and engineers from mechanics, detector and electronic groups who made possible this experimental program, especially at LPC Caen, GANIL and Saclay. We are also indebted to the target laboratories of LNS Catania, IFU Sao Paulo, GSI Darmstadt and IPN Orsay for the excellent quality of their products. The users support group (Service des aires) and the accelerator staff of GANIL should be mentioned for their continuous support and efficient performances at all stages of this work. Finally, a great help was provided by the data acquisition group GIP (GANIL) for the electronics and by the CHARISSA collaboration for lending us their VXI ADC’s.


format: text

text:
# Do Young Neutron Stars Which Show Themselves As AXPs, SGRs and Radio Pulsars Accrete?

## 1 Introduction

In the last few years soft gamma ray repeaters (SGRs) and anomalous X-ray pulsars (AXPs), which are considered as a separate class of X-ray pulsars, have attracted much attention (see Mereghetti 2001 and references therein). The only important difference between SGRs and AXPs is that SGRs have active periods showing gamma ray bursts. A puzzling property of these objects is that their X-ray luminosities (Lx) are 1-3 orders of magnitude higher than their rate of rotational energy loss (Ė). Besides, repeating gamma ray bursts of SGRs is the most important property of this class of objects. In order to explain the main physical peculiarities of these objects, not only a theory overstepping the limits of the standard NS/pulsar physics is needed, but also examining another type of NS, known as radio quiet NS, with low X-ray luminosity is required.

There is no radio radiation observed from AXPs/SGRs, only the upper limits for radio fluxes (and luminosities) are known (Gaensler et al. 2001). It does not necessarily mean that these objects do not radiate in the radio band only because of failing to observe radio pulses, because most of the pulsars are born with low luminosities (Lyne et al. 1998; Allahverdiev et al. 1997) and also because of the beaming factor. However, if somehow it is shown that AXPs and SGRs do not radiate in the radio band, then there will be a sharp drop in the radio luminosity of pulsars including AXPs and SGRs as we go to higher magnetic fields. It is difficult to check this possibility because it is not easy to increase the number of such observational data as the number of AXPs and SGRs are small and they are located at large distances.

The so called magnetar model is proposed by Thompson and Duncan (1995) as a potential explanation for the nature of AXPs and SGRs. According to this model these sources are isolated rotating NSs with surface dipole magnetic fields of 1014-1015 G, extracting their luminosity from the decay of the magnetic field.

The alternative model is that the X-ray luminosity is due to accretion from a fall-back disk that is left over from the supernova that formed the NS (Chatterjee et al. 2000; Alpar 2001). In this work we will show the incompetence of the $`existing`$ accretion models in explaining the AXPs, SGRs, and also radio pulsar properties. In section 2, the similarities and differences between SGR bursts and flares of flare stars are shown qualitatively and quantitatively. In section 3, we investigated the possibility of a fossil disk remaining around the NS after the SN explosion, and if such a disk exists, the effects of this disk on the observed properties of the NS. In section 4, we compared the specific angular momentum transfer in AXPs and SGRs with the specific angular momentum transfer in X-ray binaries. In section 5, after a brief review of history of accretion theory, we showed that the proposed propeller mechanism for radio pulsars can not account for the evolution of radio pulsars on the P-Ṗ diagram. The $`existing`$ fall-back disk models proposed to explain the nature of AXPs and SGRs are not well developed that some predictions of these models contradict some of the observational data. These will be shown below.

## 2 Bursts and Flares Due to Magnetic Activities

In general bursts, and in particular gamma-ray bursts, are common astronomical phenomena. In this section we compare SGR bursts with the bursts of well known objects in astrophysics. Bursts are observed, for instance, for the Sun, flare stars (UV-Cet type), and T-Tauri stars. Flares have fast rise-times of luminosity and longer decay times. For all of the flaring objects, the flares are related to the magnetic activities and have short characteristic times ($``$minutes for UV-Cet stars). Intense flares are followed by less-intensive ones. Durations of SGR bursts (10 ms-10 s) bear resemblance to (taking into account the very small size of the NS) the durations of flares; roughly, duration of a flare is directly proportional to the size of the star. Flares occur on a large portion of the surface of flare stars. This must be also true for SGRs since the bursts, particularly the giant bursts, can not be produced on a small portion of the neutron star’s surface. Radius of a typical UV-Cet type star is about Rflare=1010 cm, and of a neutron star is about RNS=106 cm. Taking typical flare duration ($`\mathrm{\Delta }t_{flare}`$) to be 1 minute, the duration of the SGR burst can be estimated to be 60 ms. This value is consistent with the observed $`\mathrm{\Delta }t_{SGR}`$ values. The difference between flares and SGR bursts is the power output, and time intervals between the strong flares (for UV Cet stars it is, on the average, about 10-100 times shorter compared to SGR bursts). The ratio of the energy of outburst to the persistent radiation energy for flare stars (UV-Cet type) can be as large as 104 and for SGRs 108. Since AXPs and SGRs are considered as one class of objects, we must enquire if AXPs also have outbursts, however, not frequently and/or with low intensity.

SGR bursts have very large energies in the range $`10^{4044}`$ erg with characteristic timescales of 10 ms-10 s. Such bursts can not be brought into existence by the $`existing`$ accretion (onto a NS) models and such bursts have never been observed in accreting binary systems even though X-ray luminosities of some binaries are about two orders of magnitude higher than the persistent X-ray radiation of AXPs and SGRs. It must be noted that, the X-ray bursts in low mass X-ray binaries (LMXBs) are due to accretion powered thermonuclear reactions on the surfaces of old NSs (Lipunov 1992). For these bursts to occur the magnetic field of the burster should be small. The burster characteristic has not been observed for high mass X-ray binaries (HMXBs). For transient X-ray binaries the bursts are due to change in the accretion rate. The source of these bursts observed from transients and bursters, which have very different physical mechanisms, is not the magnetic field of the NS, but accretion. Moreover, the spectra of the bursts observed from X-ray binaries and SGRs are different.

As we see above, characteristics of the SGR bursts are similar to the characteristics of flares. This is an evidence that the origin of SGR bursts is the magnetic field of the NS. If this is the case, magnetic field of the NS must be $``$1014-1015 G to power the SGR bursts with energies of $``$1040-1044 erg (Thompson & Duncan 1995). The $`existing`$ accretion models avoid to explain the SGR bursts and this is the most important handicap of such models.

## 3 The Fall-Back Disk Model

### 3.1 Can a Disk Remain Around a Neutron Star After Supernova Explosion?

The X-ray luminosity (Lx) of AXPs and SGRs can be explained by accretion, however, the compactness of the values of Lx is achieved due to the fact that the parameters Ṁ and B are adjustable in a wide range in the accretion models.

In the fall-back disk models, it is very important to know how much mass of gas may remain close to the NS after the supernova explosion. Gravitational energy of a NS with one M and a radius about 1.5$`\times 10^6`$ cm can be estimated as (1-2)$`\times 10^{53}`$ ergs which is about 5-10$`\%`$ of Mc2. The gravitational energy is converted into heat and rotational energy, and this heating energy, about 1053 ergs, transforms into the energy of neutrinos and antineutrinos during supernova explosion (Zeldovich and Guseinov 1965, Guseinov 1966). From the study of SNRs, it is known that their kinetic energies lie in the interval $``$3$`\times 10^{49}10^{51}`$ ergs. From the observations of SNe it is known that the explosion energy, on the average, lies in the interval 1050-1051 ergs and the thrown out mass of gas is $``$(0.5-4)M. This shows that only a small part of the energy of the neutrinos transform into explosion energy, but, in principle, this is enough to sweep out all the mass which are present near the NS. On the other hand, whether some fall-back matter remain near the NS (and also its quantity) depends on the velocity distribution of the thrown out matter.

The part of the gravitational energy which is transformed into rotational energy of the NS mainly depends on initial value of angular momentum (and its distribution) of the collapsed star. Even if we know these values, in order to calculate the rotational energy of the NS we must also know the value of the magnetic viscosity as well as the dynamics of the collapse and the explosion. At this point, there are many uncertainties and therefore angular momentum and rotational energy of the NS should be estimated from the initial period of pulsars. Since, the initial period is about 10 ms, initial rotational energy of the NS must be about 5$`\times 10^{51}`$ ergs. The collapse transforms the angular momentum from the NS (which is about to be born) to external parts and this can sweep out the surrounding matter (Bisnovatyi-Kogan 1971; Amnuel et al. 1973). Actually, observational data do not show presence of gaseous disks or planets near single pulsars with a magnetic field in excess of 1010 G. Optical observations given in Kaplan et al. (2001) showed that there is no optical counterpart of SGR 0526-66 and no accretion disk is present around the NS. Moreover, due to precession of the disc around the NS and/or inhomogeneities in the disc, variations in the dispersion measure values of PSRs must have been observed at least for some of the nearby young PSRs. Therefore, there is no basis for the presence of any matter of fall-back around AXPs and SGRs, but we can not completely exclude the possibility of a little bit of gas, which can not be observed, remaining near the neutron star.

### 3.2 Are AXPs and SGRs Accreting Systems?

As we indicated above, an unobservable amount of matter might remain around the NS. Below, we examine if such a small amount of matter spins up or spins down the NS. When the collapse begins to slow down, the SN explosion occurs and the shell around the NS is thrown away. Due to conservation of angular momentum, the rotational speed of the shell will decrease as it moves away from the NS. At the beginning, the magnetic field of the NS is frozen to the shell. This slows down the NS, and speeds up the shell (Bisnovatyi-Kogan 1971 and Amnuel et al. 1973). Angular speed of the shell falling back increases. Consequently, the NS will be spun-up by the fall-back matter if the propeller mechanism is not working simultaneously, i.e. if the fall-back matter does not lose angular momentum. However, the $`existing`$ accretion from a fall-back disk models do not explain how the fall-back disk loses angular momentum instead of gaining it. Models based on fall-back disk should include propeller mechanism as well as accretion mechanism working together simultaneously.

In the $`existing`$ fall-back disk models the parameters of the disk (related to the magnetic field and rotational speed of NS) is chosen such that the X-ray luminosity is about $`10^{3436}`$ erg/s and Ṗ $`10^{13}10^{11}`$ s/s. So, the $`existing`$ models, in fact, are not related to the properties and history of the fall-back matter. On the other hand, the only accretion model (onto a 1012 G NS) which takes into account the small ages of AXPs and SGRs is the model of accretion from a fall-back disk.

## 4 Comparisons Between AXPs/SGRs and X-ray Pulsars in Binaries

We examined the differences between Lx/$`|`$Ė$`|`$ values of AXPs, SGRs, and X-ray pulsars in binary systems. Here, Lx is the X-ray luminosity which depends only on the rate of accreted mass, Ṁ${}\_{x}$, and Ė is the rate of rotational energy change which depends on both Ṁ${}\_{x}$ and specific angular momentum, Lspc. Therefore, Lx/$`|`$Ė$`|`$ will be inversely proportional to Lspc. In Table 1, Lx and $`|`$Ė$`|`$ values (and their ratio) of X-ray pulsars in high mass X-ray binaries (HMXBs), a LMXB, and AXPs/SGRs with known Ṗ values are displayed (we do not include transient X-ray binaries in this table, because it is difficult to determine the correct value of Lx/$`|`$Ė$`|`$ for such systems). For HMXBs, Lx/$`|`$Ė$`|`$ values range from $`3\times 10^3`$ to $`3\times 10^6`$, with one exception, namely H0115-737, which has a lower Lx/$`|`$Ė$`|`$ value of 65. This is because of the small spin period value of this pulsar. There is only 1 LMXB with reliable Ṗ and Lx values and its Lx/$`|`$Ė$`|`$ value is (5-7.5)$`\times 10^2`$, less than all of the Lx/$`|`$Ė$`|`$ values of HMXBs which have spin periods close to spin periods of AXPs and SGRs.

Accretion onto a NS can spin it up or down. However, for binary systems with the same parameters, we expect that as Ṁ${}\_{x}$ increases Ṗ, $`|`$Ė$`|`$ and Lx increase. There is disk accretion in LMXBs and unit mass of the accreted matter has more angular momentum compared to the unit mass of the accreted matter in HMXBs since, on the average, they have higher orbital velocity. Because of this reason, for LMXBs Lx/$`|`$Ė$`|`$ is smaller than for HMXBs, on the average. Without the propeller effect, value of Lx/$`|`$Ė$`|`$ for fall-back matter accreted onto the NS must not be smaller than for X-ray pulsars in HMXBs if their spin period values are close to each other. As seen from Table 1, 4 of 7 AXPs/SGRs seem to have lower values than LMXBs that accretion from fall-back matter without propeller effect contradicts the observational data. On the other hand, the propeller effect must be weak enough not to be observed in the plerionic parts of the SNRs, because the SNRs which have genetic connections with AXPs have pure shell type structures (Gaensler et al. 2001; Tagieva & Ankay 2002).

## 5 The Other Aspects of Accretion: Single Stars

Classical accretion (onto a single star) theory is a very developed and well known subject in astrophysics. To explain the radiation from stars, in early times, accretion from interstellar medium and contraction of stars were proposed. Later it was understood that thermonuclear reactions in the cores of stars are the source of radiation of stars, then there was no need for accretion anymore. On the contrary, it was found that stars have winds, sometimes very strong winds.

Salpeter (1964) tried to explain the X-ray sources (when they were first found by Giacconi et al. 1962) based on accretion from interstellar medium onto a single neutron star. So, accretion onto single stars was called to mind. It is well-developed for systems including NSs and black holes. However, accretion from interstellar medium is found to be wrong, because it is found that if there is accretion onto a single NS, then the NS becomes hotter and the speed of sound in the surrounding matter increases, so that the accretion rate decreases by 6 (Schwartsman 1970) or 8 (if the magnetic field frozen in the interstellar medium is also taken into account, Amnuel and Guseinov 1972) orders of magnitude. X-ray sources were explained as accreting binary sytems (Zeldovich and Guseinov 1966).

After PSRs were discovered, it was understood that rapidly rotating single NSs also produce winds (Pacini 1967). Such winds are much more efficient than radiation pressure in sweeping out the surrounding matter around PSRs (Lipunov 1992). Then accretion theory could only work for slowly rotating single NSs with low magnetic fields. For 30 years, this has been searched for accreting old single NSs in optical and soft X-ray bands without considerable success (Danner 1998). Also there was no success in finding fluctuating (accreting from interstellar medium) single black holes (Schwartsman 1970). Moreover, there is no sign of fall-back matter in the central parts, particularly in the near-environments, of point sources of historical (i.e. very young) SNRs Crab (PSR J0534+2200), 3C58 (RX J0201.8+6435, Torii et al. 2000; Bocchino et al. 2001), Cas A (CXO J2323+5848, McLaughlin et al. 2000; Kaplan et .al 2001).

### 5.1 Propeller Mechanism for Radio Pulsars

Despite the facts against accretion onto single NSs given above, Alpar et al. (2001) following Menou et al. (2001) have proposed that the distribution of radio pulsars in P-Ṗ diagram can be explained by the combination of magnetic dipole torque and propeller torque of a fall-back disk, instead of pure magnetic dipole radiation torque. Alpar et al. (2001) predicted the evolutionary tracks (represented with the 2 curves in Figure 1) of pulsars on the P-Ṗ diagram. Constant B-field lines are also shown in Figure 1 to make a comparison between them. According to Peng et al. (1982) and Huang et al. (1982) neutrino emission from pulsars and magnetic dipole radiation of superfluid neutrons, respectively, also yield tracks similar to those of Alpar et al. (2001). As seen in Figure 1, down to a minimum period derivative value, pulsars follow dipole-dominant radiation tracks and after that point pure dipole and combined dipole+propeller tracks diverge from each other. Period derivative becomes Ṗ $`P^3`$ in the propeller-dominant phase. This diversion also reflects that, after the dipole-dominant phase is over, the ages predicted by the model of Alpar et al. (2001) and the characteristic ages (which are determined by the effect of pure magnetic dipole radiation torque) follow different paths as can be seen in Figure 1.

The real ages of pulsars are the kinetic ages which are valid for all models; the kinetic age is proportional to the distance of the PSR from the Galactic plane (see e.g. Lyne & Graham-Smith 1998). So, the dipole+propeller model can be tested by comparing the model’s age predictions with the kinetic ages. According to the dipole+propeller model, older pulsars must be located in the upper right part of the P-Ṗ diagram, and these pulsars must be far away from the Galactic plane, i.e. their kinetic ages must be larger. In order to check this, we constructed the P-Ṗ diagram by representing the pulsars with $`|z|`$ $`<`$ 200 pc with + symbol and the pulsars with $`|z|`$ $`>`$ 400 pc by open circles in Figure 1, where $`|z|`$ is the distance of pulsar from the Galactic plane.

In general, the birth places of pulsars are very close to the Galactic plane. The average scale height at the time of birth is about 60 pc, similar to OB stars’ average scale height. However, in the outer parts of the Galaxy, the star formation regions might deviate from the geometric plane of the Galaxy. Optical observations of cepheids with high luminosities and of red supergiants located at distances about 5-10 kpc from the Sun, in the directions l$`200^o330^o`$, showed that the star formation regions are located below the Galactic plane by about 300pc, and the star formation regions in the directions l$`70^o100^o`$ are located above the plane by about 400pc. The star formation regions at about 3-5 kpc from the Sun in the directions l$`270^o320^0`$ are located about 150 pc below the geometric plane of the Galaxy (Berdnikov 1987). These deviations of the locations of star formation regions from the Galactic plane have strong influence on the kinetic ages of young pulsars. We did not include the pulsars located in the deviated parts of the star formation regions indicated above. We also did not include the pulsars beyond 5 kpc due to large uncertainties in the distance measurements. As seen from Figure 1, along the path (from left to right) of the dipole+propeller tracks, kinetic ages first increase and then start to decrease. This is inconsistent with the propeller model. Therefore, dipole+propeller model can not account for the evolution of pulsars on the P-Ṗ diagram. Age predictions due to pure magnetic dipole radiation mechanism are in good agreement with the kinetic ages (see Figure 1). Also, Allakhverdiev & Tagieva (2002) showed that magnetic dipole evolutionary tracks are more reliable; pulsars evolve roughly with constant magnetic field, number of pulsars increase with characteristic time on constant magnetic field lines.

Only a few of the pulsars with characteristic ages less than $`10^6`$ years are more than 400pc away from the plane (see Figure 1) which could be related to their place of birth being high above the Galactic plane (the progenitor can be a runaway star) and their speed being large.

There is also another possibility to test the propeller model of Alpar et al. (2001) which predicts the change of Ṗ (P̈) to be the highest on the propeller-dominant parts of the evolutionary tracks. In order to find the braking index (n=$`\ddot{\mathrm{\Omega }}\mathrm{\Omega }/\dot{\mathrm{\Omega }}^2`$, where $`\mathrm{\Omega }`$ is the angular frequency of the pulsar) the value of P̈ must be large. About 40 years of observations have shown that, in that part of the diagram, P̈ values of pulsars are not large enough to be measured. This also contradicts the propeller model.

Glitch is a common phenomenon for young radio pulsars. According to the model of Alpar et al. (2001) for pulsars with P$`>`$ 1 s and Ṗ$`>`$ 10-14 s/s no glitch must be observed, because according to their model these pulsars must be older than about 107 yrs. If glitches from these pulsars are observed, then this will directly show that they are young.

## 6 Conclusions

In this work, we have focused on if a residual (fossil) disk may remain around single NSs (AXPs, SGRs, and radio pulsars) after the supernova explosion and the effects of such a disk on the NS properties. We approached the problem from various different ways. We tried to clarify the physical nature of these objects and the existing accretion models in the light of observational data. In summary, we have showed that:
1) The characteristics of SGR bursts are similar to the characteristics of flares of flare stars all of which have the same source (origin): the magnetic field of the star, but not accretion onto the star.
2) After the supernova explosion, only a little amount of mass may remain around a single NS. The explosion energy of the thrown out matter is enough to sweep out most of the mass around the NS.
3) The surrounding matter can also be swept out by angular momentum transfer due to magnetic braking and this is supported by observational data which do not show gaseous disks or planets (or the effects of these) near single pulsars with a magnetic field in excess of 1010 Gauss.
4) Even if such a little amount of mass remains around the NS, it spins up the NS (instead of spinning it down) unless the propeller mechanism is working together with accretion simultaneously. We showed this also by Lx/Ė considerations in section 4. However, the propeller effect should be weak enough not to be observed since all the AXPs genetically connected to SNRs are pure shell type remnants. Fall-back disk models should be constructed taking these facts into account.
5) Comparing the kinetic ages of pulsars with the age predictions of the propeller+dipole model for radio pulsars, we showed that there is a strong contradiction between them. So, propeller model does not work for radio pulsars.

$`Acknowledgments`$. It is a pleasure to thank Oktay H. Guseinov for fruitful and stimulating discussions. We thank Gennadi S. Bisnovatyi-Kogan for comments on the manuscript. We also thank TÜBİTAK, the Scientific and Technical Research Council of Turkey, for support through TBAG-ÇG4.

*  Haberl et al. 1995;  Hanson et al. 1989;  Burderi et al. 1998;  Levine et al. 1991;  Woo et al. 1996;  Burderi et al. 2000;  Tsunemi et al. 1996;  Mereghetti & Stella 1995;  Chakrabarty et al. 1997;  Angelini et al. 1995;  Hutchings et al. 1974;  Mavromatakis 1993;  Robba et al. 1996;  Weisskopf 1984;  Li et al. 1978;  Vrtilek et al. 1997;  Tjemkes et al. 1986;  Yokogawa et al. 2000;  Bonnet-Bidaut & van der Klis 1981;  Kelley et al. 1983;  Clark 2000;  Clark et al. 1994;  Rubin et al. 1997;  Robba et al. 1992;  Vasisht & Gotthelf 1997;  Gotthelf et al. 1999;  Sugizaki et al. 1997;  Israel et al. 1999a;  Baykal & Swank 1996;  Fahlman & Gregory 1981;  Kaspi et al. 1999;  Morini et al. 1988;  Corbet & Mihara 1997;  Seward, et al. 1986;  Mereghetti, 1995;  Oosterbroek et al. 1998;  Kouveliotou et al. 1998;  Kouveliotou et al. 1999;  Sonobe et al. 1994;  Marsden et al. 1999;  Israel et al. 1999b;  White et al. 1996;  Israel et al. 1994.


text:
# A Coupled Dynamical and Chemical Model of Starless Cores of Magnetized Molecular Clouds: I. Formulation and Initial Results

## 1 Introduction

Through a close interplay between theory and observation, a “standard” picture for the formation of isolated, low-mass stars has emerged (Shu, Adams & Lizano 1987). At the heart of this picture lie the so-called “dense cores”, studied extensively by Myers and coworkers (e.g., Myers 1999). These dense cores are intimately associated with star formation, with roughly half of them already harboring infrared sources (Beichman et al. 1986). The other half are termed “starless cores”, and they are the focus of our investigation.

The starless cores are thought to be condensed out of strongly magnetized, turbulent background clouds. Direct Zeeman measurements to date, as compiled by Crutcher (1999), suggest that the field strength is close to the critical value required for the molecular cloud support, after likely geometric corrections (Shu et al. 1999). Core formation can be driven by either turbulence decay (Myers 1999) or ambipolar diffusion. We shall focus on the latter, which has been studied quantitatively by many authors (e.g., Nakano 1979; Lizano & Shu 1989; Basu & Mouschovias 1994; Ciolek & Mouschovias 1994). A general conclusion is that dense cores are formed on a time scale several times the free fall time of the background. This relatively long formation time should leave a strong imprint on the chemistry of the cores.

The goal of our investigation is to predict the spatial distributions of various molecular species at different stages of cloud evolution. For this initial study, we will concentrate on CO, CS, CCS, NH3, N2H+ and HCO+, which are often used to probe the physical conditions of star-forming clouds. We take into account the feedback of chemistry on the cloud dynamics, through ion abundances, which regulate the ambipolar diffusion and thus cloud contraction rate. Through a detailed comparison of the predicted and observed chemical abundances of starless cores, we seek to provide further support for the standard picture of isolated star formation involving ambipolar diffusion. A starless core well suited for such a purpose is the dense core of L1544, for which an extensive data set is becoming available (e.g., Caselli et al. 2001a,b).

The L1544 starless core is located in the eastern part of the Taurus molecular cloud complex, at an estimated distance of 140 pc (Elias 1978). It is an elongated core that shows substantial infall motions (up to $`0.1`$ km s-1) on both large ($`0.1`$ pc) and small ($`0.01`$ pc) scales, based on single-disk observations of CO, 13CO, C18O, CS, C34S, HCO+, and N2H+ (Tafalla et al. 1998) and interferometric observations of N2H+ (Williams et al. 1999). Interferometric observations of CCS by Ohashi et al. (1999) reveal a ring-like structure, which shows evidence for collapse as well as rotation. The low central emission is most likely due to a depletion of CCS, judging from the fact that dust continuum peaks inside the ring (Ward-Thompson, Motte & Andre 1999). Indeed, the column density distribution of the L1544 core is rather well determined111See, however, Evans et al. (2001) who find that a singular isothermal sphere cannot be ruled out from the modeling of SCUBA data alone, when the lower temperature at the core center relative to edge is taken into account., not only from the dust continuum emission, but also from absorption against a mid-infrared background (Bacmann et al. 2000). The core has a small, high column (and volume) density central plateau (of radius $`2900`$ AU), surrounded by an envelope in which the column (and volume) density decreases rapidly outward, before it joins rather abruptly with the background at a radius of order $`10^4`$ AU (Bacmann et al. 2000). The plateau-envelope structure has been interpreted in terms of the ambipolar diffusion-driven dynamical evolution of a magnetized cloud (Li 1999; Ciolek & Basu 2000), and the large column density contrast between the central plateau and the background ($`20`$; Bacmann et al. 2000) suggests that the L1544 core is rather evolved (Williams et al. 1999), and may be on the verge of forming a star (or stellar system). Further evidence for the core being in an advanced stage of evolution comes from IRAM 30m observations of Caselli et al. (1999), who found that CO is depleted by a factor of order $`10`$ at the dust continuum peak. These authors also noted that the (likely) optically thin lines of D13CO+ and HC18O+ are double-peaked, which may be another indication that these species are centrally depleted and their emission comes mainly from an infalling shell. The wealth of molecular line and dust continuum data, coupled with the recent Zeeman measurement of magnetic field strength (Crutcher & Troland 2000) and the determination of magnetic field direction from submillimeter polarization measurements (Ward-Thompson et al. 2000), makes L1544 a Rosette Stone for theories of dense core condensation leading to isolated low-mass star formation.

Many of the dynamical and chemical features of L1544 are shared by L1498, another quiescent, elongated starless core also located in the Taurus molecular complex. Lemme et al. (1995) studied the core in detail, and found that the distributions of both C18O and CS are ring-like, indicating a central depletion of these two species. The CO depletion factor is inferred to be of order $`8`$ or higher (Willacy, Langer & Velusamy 1998), comparable to that of L1544, based on ISOPHOT 100 and 200 $`\mu `$m observations. Single-dish and interferometric observations of CCS (Kuiper, Langer & Velusamy 1996) show that its distribution is also ring-like. Moreover, the NH3 emission appears to be centrally peaked, again similar to L1544. Perhaps most strikingly, there is evidence for infall motions on a large ($`0.1`$ pc) scale in the double-peaked profiles of both optically thick CS lines and the likely optically thin C34S(2-1) transition (Lemme et al. 1995); the latter could come from an infalling shell, as in L1544 (Caselli et al. 1999). The optically thin tracer N2H+ shows, however, a single peaked spectrum (Lee, Myers & Tafalla 2001), indicating that it is centrally peaked, in contrast to C34S. These similarities to L1544 point to a relatively advanced stage of pre-protostellar evolution for L1498, making it another valuable testing ground for chemical and dynamical models of (low-mass) star formation.

Common to both sources is the ring (or shell)-like structure of CCS, a well-known “early time” species (Bergin & Langer 1997) tracing chemically “young” material of order $`10^5`$ yrs or less. This prominent feature should provide a strong constraint on models. Additional constraints should come from the inferred large CO depletion near the core center as well as the spatial extent and speed of the infall motions present in both L1544 and L1498. The question we want to address is: can the standard picture of dense core formation involving ambipolar diffusion satisfy these chemical and kinematic constraints simultaneously? A partial answer to this question has been provided by Bergin & Langer (1997), who considered the chemical evolution of a parcel of cloud material whose density increases with time in a prescribed way (i.e., a standard time dependent but depth independent model). They found that differential depletion of various molecular species onto dust grains must be a key ingredient in understanding the differentiated chemical structure observed in L1498.

Aikawa et al. (2001) went one step further, and followed the evolution of molecular abundances of a dynamically collapsing starless core as a whole, with an emphasis on the spatial distribution. The study was based on the self-similar solution of Larson (1969) and Penston (1969) and its artificially delayed analogs. Among the dynamical models adopted, they found that the undelayed Larson-Penston model provides the best overall fit to the observational data on CO, CCS and N2H+ for L1544. If this result is robust, its implications would be far reaching. It is well known that the Larson-Penston solution describes the collapse of a cloud essentially on a free-fall time scale. The fact that it matches observations better than its delayed analogs appears to pose a serious challenge to the standard picture of star formation involving strongly magnetized cloud cores, which are formed through ambipolar diffusion on a time scale much longer than the free-fall time scale. However, the large infall velocity of the Larson-Penston solution, approaching 3.3 times the sound speed, is clearly incompatible with those inferred for L1544 and other starless cores, which are typically less than the sound speed (Lee et al. 2001). Furthermore, the best-fit model of Aikawa et al. under predicts the abundance of N2H+ by a large factor of $`20`$. These discrepancies motivate us to reexamine the molecular evolution of starless cores using a more sophisticated dynamical model, one that is coupled to the cloud chemistry.

The dynamical model we will adopt is that of Li (1999). It is a spherical model of starless cores that takes into account the dynamic effects of magnetic field approximately. As envisioned in standard scenario of low-mass star formation, the cores in the model evolve through ambipolar diffusion, and their dynamics are coupled to chemistry through the abundances of charged species. Previously, we (Shematovich, Shustov, & Wiebe 1997; 1999) have modeled the chemistry of dynamically evolving clouds in detail. The chemical model will be combined with the dynamical model into a coupled dynamical-chemical model of magnetized starless cores. An advantage of the combined model is that it allows for a determination of not only the time evolution of various molecular species but also their spatial distributions and kinematics. We describe the formulation of the combined model in §2. In §3, we present representative model results, focusing on the density distribution, velocity field and abundance distributions of several commonly used molecular species. These results are compared with observations of L1544. We conclude and discuss our main results and future refinements in §4.

## 2 Formulation of Coupled Dynamical-Chemical Model

### 2.1 Dynamical Model

Ambipolar diffusion-driven evolution of strongly magnetized clouds has been modeled by many authors, although a complete understanding is still lacking. The numerical complexity of the problem forced most workers to adopt either a quasi-static approximation (e.g., Nakano 1979; Lizano & Shu 1989) or a simplified disk-like geometry (e.g., Basu & Mouschovias 1994; Ciolek & Mouschovias 1994). Some 2D axisymmetric dynamical models are available (Fiedler & Mouschovias 1993; Desch & Mouschovias 2001), although incorporation of a full-blown chemical network into such models would be prohibitively expensive computationally at present. A simpler approach to the dynamical problem was pioneered by Safier, McKee & Stahler (1997). By retaining only the pressure component of the magnetic forces they were able to retain a spherical geometry, which simplifies the dynamics greatly. Indeed, they were able to obtain analytic solutions to both core formation and collapse assuming further that the thermal pressure is negligible and that the ions (and magnetic field lines tied to them) move much more slowly than neutrals. These last two assumptions have been relaxed in Li (1998) and the resulting set of governing equations can be solved efficiently, particularly in a Lagrangian form (Li 1999). The formulation of the dynamical problem by Li (1999) is employed in our coupled dynamical-chemical model.

For reference, we list below the equations governing the dynamic evolution of a spherical, magnetized cloud (taken from Li 1999) in dimensionless radius $`\zeta `$, mass $`m`$, density $`\widehat{\rho }`$, velocity $`u`$, magnetic field strength $`b`$ and time $`\tau `$:

$$\frac{\zeta }{m}=\frac{1}{\widehat{\rho }\zeta ^2},$$
(1)
$$\frac{u}{\tau }=\frac{m}{\zeta ^2}\zeta ^2\frac{}{m}\left(\frac{\widehat{\rho }}{2\alpha _\mathrm{c}}+\frac{b^2}{2}\right),$$
(2)
$$\frac{}{\tau }\left(\frac{b}{\widehat{\rho }\zeta }\right)=\frac{}{m}\left(\frac{1.4}{\nu _{\mathrm{ff}}}\frac{b^2\zeta ^2}{\widehat{\rho }^{1/2}}\frac{b}{m}\right),$$
(3)
$$u=\frac{\zeta }{\tau }.$$
(4)
The density $`\rho `$ and field strength $`B`$ are scaled by their initial values at the cloud center $`\rho _\mathrm{c}`$ and $`B_\mathrm{c}`$. The time $`t`$ is measured in units of the initial free-fall time

$$t_{\mathrm{ff},\mathrm{c}}=\frac{1}{(4\pi G\rho _c)^{1/2}},$$
(5)
at the cloud center, and the velocity $`V`$ in units of the Alfvén speed

$$V_{\mathrm{A},\mathrm{c}}=\frac{B_\mathrm{c}}{(4\pi \rho _\mathrm{c})^{1/2}}.$$
(6)
Scales for the radius $`r`$ and mass $`M`$ are then $`r_\mathrm{c}=V_{\mathrm{A},\mathrm{c}}t_{\mathrm{ff},\mathrm{c}}`$ and $`M_\mathrm{c}=4\pi \rho _\mathrm{c}r_\mathrm{c}^3`$. There are two dimensionless quantities in the governing equations. The first, $`\alpha _\mathrm{c}`$, is the ratio of magnetic pressure to thermal pressure at the cloud center at $`t=0`$. For simplicity, it will be taken to be unity, consistent with the current data on the magnetic fields in the interstellar clouds (Crutcher 1999). The second is the magnetic coupling parameter $`\nu _{\mathrm{ff}}`$ defined as the ratio of the local free-fall time (at a given density $`\rho `$) to the magnetic field–neutral coupling time scale

$$\nu _{\mathrm{ff}}=\frac{\gamma \rho _\mathrm{i}}{\sqrt{4\pi G\rho }},$$
(7)
where $`\rho _i`$ is the total mass density of ions, and $`\gamma `$ is an ion-neutral drag coefficient. An approximate value for the drag coefficient of $`\gamma =3.5\times 10^{13}`$ cm3 g-1 s-1 is given in Shu (1992), obtained assuming a mean molecular weight of ions of $`30`$. We will adopt the above simple prescription for the coupling parameter in this initial study, although a more elaborate treatment taking into account of contributions from charged dust grains is possible (e.g., Nishi, Nakano & Umebayashi 1991; Li 1999) and will be implemented in the future.

### 2.2 Chemical Model

The chemical model used here has been described in detail by Shematovich et al. (1997, 1999). Here we summarize the main features briefly. We solve the equations of chemical kinetics describing reactions in interstellar clouds of relatively low density ($`10^210^8`$ cm-3) and temperature ($`T10100`$ K). Molecules are formed and destroyed both in gas-phase and on grain surfaces. The dust grains are supposed to be well mixed with gas and contain 1% of the cloud mass. We label different species, both neutral and charged, by $`i=0,1,\mathrm{},M`$, with $`i=0`$ denoting electrons. The evolving chemical composition is described through the number densities $`n_i^\mathrm{g}(r,t)`$ in gas-phase and $`n_i^\mathrm{d}(r,t)`$ on icy mantles. The exchange between gas and dust is driven by accretion and desorption processes. In general, the chemical abundances of all species are computed by integrating the equations

$$\frac{\mathrm{d}}{\mathrm{d}t}n_i^\mathrm{g}(r,t)=\underset{j}{}\underset{l}{}K_{lj}^\mathrm{g}n_l^\mathrm{g}n_j^\mathrm{g}n_i^\mathrm{g}\underset{j}{}K_{ij}^\mathrm{g}n_j^\mathrm{g}k_i^{\mathrm{acc}}n_i^\mathrm{g}+k_i^{\mathrm{des}}n_i^\mathrm{d},$$
(8)
$$\frac{\mathrm{d}}{\mathrm{d}t}n_i^\mathrm{d}(r,t)=\underset{j}{}\underset{l}{}K_{lj}^\mathrm{d}n_l^\mathrm{d}n_j^\mathrm{d}n_i^\mathrm{d}\underset{j}{}K_{ij}^\mathrm{d}n_j^\mathrm{d}+k_i^{\mathrm{acc}}n_i^\mathrm{g}k_i^{\mathrm{des}}n_i^\mathrm{d},$$
(9)
where $`K_{ij}^\mathrm{g}`$ and $`K_{ij}^\mathrm{d}`$ are the gas-phase and grain surface chemical reaction rates. The coefficients $`k_i^{\mathrm{acc}}`$ and $`k_i^{\mathrm{des}}`$ define the accretion and desorption rates for the $`i`$th component. These equations are integrated, along with the dynamical equations, at each point of the flow. The computed abundances are then used to determine the magnetic coupling parameter $`\nu _{\mathrm{ff}}`$, and to obtain column densities for various commonly used molecular tracers for comparison with observations.

#### 2.2.1 Gas-Phase Reactions

We consider gas-phase reactions in the dense, dark regions of molecular clouds that are well shielded from the external UV radiation field. The reactions are initiated by galactic cosmic rays which ionize the most abundant species with an adopted rate of $`1.3\times 10^{17}`$ s-1. The primary ions that are formed react quickly with abundant hydrogen and CO molecules, defining the main paths of chemical evolution.

We use the UMIST Rate95 file (Millar, Farquhar & Willacy 1997). The included elements are H, He, C, N, O, S, Si, Na, Mg, and Fe. Our chemical reaction network is mainly limited to those species containing no more than two atoms of trace elements (i.e. all elements other than H and He), with a few exceptions. One of these exceptions is C2S, which has been mapped interferometrically in L1498 (Kuiper et al. 1996), L1544 (Ohashi et al. 1999), and others. Its abundance and spatial distribution provide a crucial constraint on models, as discussed by Aikawa et al. (2001) and shown below. In total, our gas-phase network contains about 2000 reactions between 59 neutral and 83 ion species.

#### 2.2.2 Grain Surface Reactions

The formation of icy mantles due to accretion and desorption has been taken into account. The accretion rate is given by

$$k_i^{\mathrm{acc}}=\pi a^2<v_i>S_in_{\mathrm{gr}}s^1,$$
(10)
where $`a`$ and $`n_{\mathrm{gr}}`$ are the grain radius (taken to be 10-5 cm) and dust number density, $`<v_i>`$ is the average thermal velocity of the species, and $`S_i`$ the sticking probability. We assume $`S_i=0.3`$ for all neutral species except H2 and He (Willacy, Rawlings & Williams 1994). Two desorption mechanisms are included. They are thermal evaporation and cosmic ray induced impulsive thermal evaporation (Hasegawa & Herbst 1993). The dust temperature is assumed to be 7 K everywhere, close to the value inferred by Evans et al. (2001) and Zucconi, Galli & Walmsley (2001) near the center of L1544. The values of the sticking probability and desorption energies are uncertain. The sensitivity of our model to these uncertainties will be explored elsewhere.

After being accreted onto the grain surface, the species are allowed to react with one another mainly through the hydrogen addition reactions, with rates taken from Hasegawa, Herbst & Leung (1992). The surface products include saturated molecules (such as CH4, H2O and NH3), homonuclear molecules (such as C2, N2, and O2) and their hydrogenated forms, among others. The accretion limit to the surface hydrogen addition chemistry was taken into account through the first modification of rate equations proposed by Caselli, Hasegawa, & Herbst (1998). In addition to 53 surface reactions (including molecular hydrogen formation), we also consider the recombination of ions on grain surfaces. The products of dissociative recombination are assumed to return into the gas phase immediately. The effects of surface reactions on abundance distributions will be fully examined in a subsequent paper.

### 2.3 Initial and Boundary Conditions

We adopt as a starting point for our calculations an isolated molecular clump in mechanical equilibrium, with a magnetic pressure equal to the thermal pressure everywhere. The gas temperature is assumed to be 10 K and is kept constant throughout the computations. The initial molecular hydrogen central density is set to $`n_{\mathrm{H}_2}=10^3`$ cm-3 for all models. Without ambipolar diffusion, the cloud would remain in the equilibrium forever. At $`t=0`$ ambipolar diffusion is switched on. The cloud evolution is then followed numerically with reflection boundary conditions (i.e., zero spatial gradient) at the cloud center and with a free pressure boundary at the cloud outer edge.

The initial chemical composition is assumed to be mostly atomic. Species with non-zero initial abundances are listed in Table 1. Their abundances are given with respect to the number density of hydrogen nuclei. All hydrogen is assumed to be in a molecular form initially.

### 2.4 Numerical Method

We divide the computational domain into spatial zones of equal mass. At each time step, we first advance the chemical model to evaluate ion density $`\rho _\mathrm{i}`$ in each spatial zone. Equations (1), (2), and (4) are then solved to update the velocity, radius of zone boundary and density. A new distribution of magnetic field is obtained by solving the diffusion equation (3), with the magnetic coupling parameter calculated from the ion densities at the beginning of the time step.

The value of time step is computed from the usual magnetic Courant-Friedrich-Lewy condition to ensure dynamic stability. The chemical model is advanced over this time interval independently, with its own algorithm and smaller time substeps that are determined from the required abundance accuracy. We solve the chemical kinetics equations with the standard LSODE package, and terminate the computation when the central density reaches $`n_{\mathrm{H}_2}=10^6`$ cm-3. This density is close to the value inferred at the center of L1544. At higher densities the effects of charged dust grains on $`\nu _{\mathrm{ff}}`$ may become dominant (Nishi et al. 1991), and such effects have not been taken into account in the current version of the code.

## 3 Representative Model Results

The initial conditions for forming starless cores are not well constrained, either observationally or theoretically. There is in principle a large parameter space to be explored. In this initial study, we will focus on two models that highlight the effects of magnetic field on the dynamics and chemistry. In the first model, we consider a standard, magnetic cloud with a mass of 20 M (StM+B). The second model is identical to the first except that the magnetic field strength is set to zero (i.e., non-magnetic; StM–B). For reference, we also consider the dynamics of a model that is the same as the first except that the magnetic coupling parameter is fixed at the canonical value of $`\nu _{\mathrm{ff}}=10`$ instead of being computed self-consistently. These models are summarized in Table 2.

### 3.1 Dynamics

The main dynamic quantities of the models (density distribution and velocity field) are displayed in Figs. 13, along with the ionization fraction and magnetic coupling parameter for the StM+B model. As mentioned earlier, all models start with a central density of $`n_{\mathrm{H}_2}=10^3`$ cm-3, typical of clouds in transition from being diffuse to dark. Note that when a certain center density is reached, the difference in density profile is modest, with that in the non-magnetic model somewhat peakier than those in the magnetic models. The difference in the time it takes to reach a given density is more pronounced. In Table 3 we list the time intervals taken to increase the central density by successive factors of 10. The time for the first factor-of-10 increase, to a value of 104 cm-3 characteristic of dense molecular cores, is $`1`$ Myr for the non-magnetic model. It is substantially shorter than those for the magnetic models ($`45`$ Myrs). The times for subsequent 10-fold increases become shorter and shorter, and the differences in time between the models also get smaller. The time for the second (and third) factor-of-10 increase is longer for the magnetic models than for the non-magnetic one by only a factor of $`2.5`$ (and $`2`$). In other words, after formation a magnetic starless dense core typically lasts for only a few free-fall times before stellar birth. That the evolutionary time scale becomes closer to that of a non-magnetic, free-falling model at higher densities is an important feature of the ambipolar diffusion-driven cloud evolution. This characteristic cannot be captured by simply adding a constant delay factor to a non-magnetic model, as done for example in Aikawa et al. (2001). The differences in time scale among the models will have a decisive influence on the chemistry, with potentially observable consequences.

Another potentially observable difference among the models is the velocity field. From Figs. 13 we find that the infall speed is typically a few tenths of 1 km s-1, comparable to the isothermal sound speed of the cloud ($`a=0.188`$ km s-1 at 10 K). As expected, the non-magnetic model collapses right away, approaching a maximum speed $`2.3`$ times larger than the sound speed $`a`$ by the time of a $`10^3`$-fold increase in the central density. Even this maximum speed is substantially lower than that in a Larson-Penston flow, which is 3.3 times the sound speed. The infall speed in magnetized models increases more gradually, reaching a maximum value of $`0.35`$ km s-1 (or about twice the sound speed) for the StM+B and reference model when the central density reaches $`n_{\mathrm{H}_2}=10^6`$ cm-3, roughly the density inferred for the central region of L1544 (Ward-Thompson et al. 1999). Note that the density distribution flattens near the center for all three models, more so in the magnetic models than in the non-magnetic one. The radius of the central flat region is of order $`23\times 10^3`$ AU for the StM+B model, consistent with those deduced by Ward-Thompson et al. (1999) and Bacmann et al. (2000). The predicted peak H2 column density of $`10^{23}`$ cm-2 also agrees with the value inferred by Ward-Thompson et al. (1999) within a radius of $`900`$ AU. However, the observationally inferred infall velocity for L1544 of order 0.1 km s-1 (Tafalla et al. 1998; Williams et al. 1999; Ohashi et al. 1999) is lower than predicted. The discrepancy, although substantial, is less severe in our model than in those of Aikawa et al. (2001). It could be further reduced by taking into account the projection effects in a more realistic flattened geometry and the effects of magnetic tension on cloud dynamics (Ciolek & Basu 2000).

### 3.2 Ionization Fraction and Magnetic Coupling Parameter

Before showing results on chemical abundances, we comment on the ionization fraction $`x_\mathrm{i}`$ (defined as the ratio of the number density of all ions to that of hydrogen nuclei) and magnetic coupling parameter $`\nu _{\mathrm{ff}}`$, both of which are shown in Fig. 2 for the StM+B model. The fraction $`x_\mathrm{i}`$ decreases with density, as expected. Its lowest value of $`3\times 10^9`$ at the center of the cloud is somewhat higher than that inferred for L1544 by Caselli et al. (2001b), although the inference is model dependent. The ionization fraction yields a magnetic coupling parameter $`\nu _{\mathrm{ff}}`$ within a factor of two of the canonical value of 10.

### 3.3 Chemical Abundances and Comparison with L1544

#### 3.3.1 Observationally Inferred Abundances in L1544

We concentrate on the abundances of six species, including CO, CS, CCS, HCO+, NH3 and N2H+. Other species will be considered in a subsequent paper of the series. For CO, CCS, HCO+, and N2H+, the column density distributions have been inferred from observations in L1544. Caselli et al. (2001b) found that CO is heavily depleted near the center of L1544, by a factor up to $`10`$. Its column density distribution is flat or slightly depressed near the center, with an average central value of $`1.5\times 10^{18}`$ cm-2, although large scatter exists in the data (see their Fig. 5). Both HCO+ and N2H+ are centrally peaked, with a peak column density of $`10^{14}`$ and $`2\times 10^{13}`$ cm-2 respectively. Ohashi et al. (1999) mapped L1544 interferometrically in CCS. They found that CCS peaks in a ring of radius $`7,500`$ AU. The peak CCS column density is $`4\times 10^{13}`$ cm-2, which is $`1.4`$ times higher than the value at the center (Aikawa et al. 2001). Observational data on NH3 and CS appear to be less complete. Using the Nobeyama 45 m telescope, Suzuki et al. (1992) observed L1544 in NH3 and derived an average column density of $`2.7\times 10^{14}`$cm-2 within a region of radius $`6,000`$ AU, corresponding to a beamwidth of $`80^{\prime \prime }`$. L1544 has been observed in CS (Tafalla et al. 1998), although the peak value and spatial distribution of its column density are not yet available.

#### 3.3.2 Standard Magnetic Model

The column density distributions of the six species of interest for the StM+B model, obtained by integrating the number densities along the line of sight at various distances from line of sight through the cloud center, are shown in Fig. 4. We will concentrate on the last time (t=5.69 Myrs) when the density distribution fits approximately that observed (Ward-Thompson et al. 1999). It is immediately evident from the figure that a central hole exists in the distribution of CCS, with a radius ($`7,000`$ AU) and depth ($`1.3`$) that match almost exactly those observed (Ohashi et al. 1999). Note the differentiated chemical structure, with the centrally peaked NH3 and N2H+ surrounded by a CCS ring, which is in turn embedded within a slightly broader CS ring222This is not exactly what is observed in L1498, where CS appears to concentrate near the inner edge of the CCS ring (Kuiper et al. 1996). The situation with L1544 is unclear, although there is some indication that a central hole exists in CS as well (Tafalla 2001; priv. comm.).. The predicted peak column density for N2H+ of $`1.3\times 10^{13}`$ cm-2 is within a factor of 2 of the observed value ($`2\times 10^{13}`$ cm-2). The predicted peak column density for NH3 of $`5.5\times 10^{14}`$ cm-2 is about twice the value $`2.7\times 10^{14}`$ cm-2 given by Suzuki et al. (1992). Averaging the column density over a radius of 6,000 AU (the beamwidth) yields a value $`2\times 10^{14}`$ cm-2, which is in a closer agreement with Suzuki et al.’s value. From Fig. 4 we see that HCO+ is also centrally peaked, although less sharply than N2H+, again in agreement with observation (see Fig. 5 of Caselli et al. 2001b). Its predicted peak column density of $`4.5\times 10^{13}`$ cm-2 is about a factor of two lower than observed ($`10^{14}`$ cm-2). The CO column density in this model shows a relatively flat distribution near the center, with a maximum value of $`1.6\times 10^{18}`$ cm-2, almost identical to the observed value. Its spatial distribution appears to be consistent with that inferred by Caselli et al. (2001b) within errors. The model predicts that the column density distribution of CS should have a shallow depression near the center. Data for testing this prediction are not yet available. We conclude that, within a factor of two or so, the predictions of the StM+B model on the abundances of CO, N2H+, HCO+ and NH3 are consistent with the data available on L1544.

The abundance and spatial distribution of CCS deserve special attention. We have mentioned the close fit to the observed spatial distribution and will now concentrate on abundance. The CCS molecule is a well-known “early time” species. Its abundance should decrease, according to Bergin & Langer (1997), rapidly after a few times $`10^5`$ years, as a result of neutral atomic carbon being converted into CO (see also Suzuki et al. 1992 and Aikawa et al. 2001). One might expect to find little CCS in the magnetic cloud of the StM+B model, which reaches the observed state after nearly 6 million years of evolution. However, this turns out not to be the case, as shown in Fig. 4. The predicted peak CCS column density is $`1.3\times 10^{13}`$ cm-2, about a third of the observed value ($`4\times 10^{13}`$ cm-2). The agreement would be closer, if one takes into account of likely projection effects. The CCS ring is seen almost edge-on, judging from its large aspect ratio ($`3`$; Ohashi et al. 1999). The nearly edge-on geometry implies a smaller peak column density perpendicular to the ring, by a factor of $`3`$. It turns out that CCS plays a crucial role in discriminating models. We will postpone a more detailed discussion of this molecule to § 4.1. Here, we simply note that the StM+B model appears to reproduce all abundance data on L1544 reasonably well. It is the “standard” model against which other models will be compared.

#### 3.3.3 Free-Falling Non-Magnetic Model

Without magnetic support, the cloud considered in the previous subsection collapses promptly, reaching the observed state in 1.56 million years. This faster dynamic evolution affects the chemical abundances, as shown in Fig. 5. The effects on CO, NH3, N2H+ and HCO+ are relatively modest. The peak column density of CO differs little from that of the standard model. Its spatial distribution appears more centrally peaked, however, which is more difficult to reconcile with observation. The distributions of N2H+, NH3 and HCO+ are centrally peaked, as in the standard model. Compared with the standard model, the peak value for N2H+ is down by a factor of $`2`$, making it significantly below the observed value (by a factor of $`4`$). The column density of NH3 is also down by a similar factor of $`2`$, making the fit to the value inferred by Suzuki et al. (1992) worse. The column density of HCO+ is lower by a modest factor of $`20\%`$.

The difference in CS and CCS is more pronounced. Unlike in the standard model, the column density distributions of CS and CCS in the non-magnetic model are centrally peaked instead of showing a central depression. While the observational situation with CS is still unclear at present, a central peak in CCS is definitively ruled out. One potentially desirable feature of the non-magnetic model is that the peak column density of CCS is about twice the value in the standard model, making it closer to the observationally inferred value without projection corrections. However, its overall fit to all available data is significantly worse than that of the standard model. This result lends some support the standard picture of low-mass star formation involving core formation through ambipolar diffusion over a time scale substantially longer than free fall. The peak value and spatial distribution of the column density of CS would further constrain the models but are not yet available. They are urgently needed.

We note that our non-magnetic model is similar to the undelayed model of Aikawa et al. (2001), except that we have included some surface reactions and adopted a smaller sticking probability.

## 4 Discussion and Conclusions

### 4.1 Magnetic Cloud Support and CCS Abundance

The success of the standard, magnetic model in fitting the data on L1544 hinges to a large extent on the relatively high abundance predicted for CCS. This result is somewhat surprising in view of the fact that CCS is a well-known “early time” species, whose abundance should decline rapidly after a few times $`10^5`$ years according to Suzuki et al. (1992), Bergin & Langer (1997), and Aikawa et al. 2001). Part of the reason for the persistently high CCS abundance in the standard model even after some 5.69 Myrs of evolution (the time it takes the cloud to reach the observed state) is that we adopted an initial cloud density of $`n_{\mathrm{H}_2}=10^3`$ cm-3, which is lower than those of Bergin & Langer ($`10^{3.5}`$ cm-3) or Suzuki et al. and Aikawa et al. ($`10^4`$ cm-3). As a result, the bulk of the $`5.69`$ Myrs is spent in the “pre-dense core” phase when the central density remains below $`10^4`$ cm-3 (see Table 3). Only a small fraction ($`10\%`$) of that time is spent at densities above $`10^4`$ cm-3. That a cloud spends a relatively short time (a few times the free-fall time) at high densities after forming a magnetically supercritical core is a general feature of ambipolar diffusion-driven evolution. It helps to keep the CCS abundance higher than one would expect based on the total time of cloud evolution.

The fractional abundances of CCS are shown in Fig. 6 as a function of radius for both the magnetic and non-magnetic models. Note the “humps” on the distributions of CCS abundance in the magnetic model, which are not apparent in the non-magnetic model. The exact origin of the humps is unclear. We suspect that they are related to the depletion of CO in the central high density region, similar to the “depletion” peak proposed by Ruffle et al. (1997) for HC3N. Ruffle et al. (1997, 1999) showed that C2H and HC3N in dense, cold cores are characterized by a secondary “late-time” maximum in their fractional abundances. The second peak is thought to be mainly caused by the increased depletion of CO, which allows C+ ion to react more readily with H2 than with oxygen bearing species. The net result is an increase in the rates of production of CH and other carbon-bearing molecules (without oxygen). It was found that high late-time C2H and HC3N maxima are achieved only when the freeze-out time scale is long compared to the chemical time scale. Because C2H and other (neutral and ionized) late-time hydrocarbons are the precursors of CCS formation in the reactions with neutral and ionized sulfur through

$$\mathrm{C}_2\mathrm{H},\mathrm{C}_2\mathrm{H}_2,\mathrm{C}_2\mathrm{H}_3,\mathrm{}+\mathrm{S}^+\mathrm{HCCS}^++\mathrm{}.,$$
(11)
$$\mathrm{C}_2\mathrm{H}_2^+,\mathrm{C}_2\mathrm{H}_3^+,\mathrm{C}_2\mathrm{H}_4^+,\mathrm{}+\mathrm{S}\mathrm{HCCS}^++\mathrm{}.,$$
(12)
$$\mathrm{HCCS}^++\mathrm{e}\mathrm{CCS}+\mathrm{H},$$
(13)
CCS should behave in a way similar to the C2H and HC3N. The competition between depletion and late-time hydrocarbon chemistry enhanced by CO-depletion appears to be the most likely cause of the humps on the radial profiles of CCS abundance.

In any case, the humps have apparently kept the CCS abundance of the magnetic model close to, and in some regions exceeding, that of the non-magnetic model. They are largely responsible for the higher-than-expected value of the CCS column density in the magnetic model. In addition, the decline of CCS abundance outside the hump, coupled with a steep decrease in the hydrogen number density with radius, may explain the observed rapid decline of the CCS column density outside the ring (Ohashi 2001; priv. comm.).

### 4.2 Depletion of Molecules onto Dust Grains

#### 4.2.1 Depletion of CO and Other Species

In the standard model, CO is heavily depleted in the central high density region. The depletion is most clearly seen in the first panel of Fig. 7, where the fractional abundance is plotted as a function of radius. The abundance is significantly below the canonical value of $`4.7\times 10^5`$ inside a radius of $`10^4`$ AU, by a factor up to $`30`$. The size of the depletion region is roughly consistent with that deduced by Caselli et al. (2001b; their Fig. 2). The depletion is less evident in the CO column density distribution (see Fig. 4), which is flat near the center. The ratio of CO and hydrogen in peak column density is a factor of $`5`$ below the canonical value. This depletion factor is about half of the value ($`9`$) inferred by Caselli et al. (2001b). Their inference is based on the observations of C17O, C18O, and dust emission (Ward-Thompson et al. 1999). In particular, an dust opacity of 0.005 cm2g-1 was adopted to convert the dust continuum flux into a hydrogen column density. At densities as high as $`10^6`$ cm-3, an opacity of $`0.01`$ cm2g-1 may be more appropriate (Ossenkopf & Henning 1994), which would lower the estimate of the hydrogen column density by half and bring a closer agreement between the model prediction and observation. Alternatively, a flattened geometry, as observed for L1544, and/or a somewhat higher CO adsorption energy (see Aikawa et al. 2001) may enhance the CO depletion in the column density distribution.

Besides CO, there are other species that are strongly depleted in the high density central region. These include CS, CCS and to a lesser extent HCO+, as can be seen from Fig. 7. The nitrogen bearing species, N2H+ and NH3, are on the other hand hardly depleted; if anything, their abundances are slightly enhanced near the center. This differential depletion pattern is in agreement with previous findings (e.g., Bergin & Langer 1997; Aikawa et al. 2001). It may have profound effects on the molecular line profiles used to probe cloud kinematics.

#### 4.2.2 Effects of Depletion on Line Profiles

An interesting observational fact about L1544 is that some optically thin lines in L1544, such as CCS (Ohashi et al. 1999) and HC18O+ (Caselli et al. 2001b), are double-peaked. As noted by Caselli et al. (1999, 2001b), the double-peaked optically thin lines can be explained in a collapsing cloud provided that the line-emitting molecules are strongly depleted near the center. We have seen from Fig. 7 that CCS and, to a lesser extent, HCO+ (and thus HC18O+) are indeed centrally depleted in our standard model. A unique strength of the coupled dynamical-chemical model is that it allows for a simultaneous determination of the spatial distribution of a given species and the velocity field, the two ingredients for line profile modeling. In a subsequent paper, we will produce synthetic line profiles for both optically thin and thick lines (such as CS and HCO+; Tafalla et al. 1998), which should enable us to constrain the models of core formation leading to star formation using both the spatial and kinematic information of molecular species.

### 4.3 Conclusions

We have developed a coupled dynamical and chemical model for the starless cores of strongly magnetized molecular clouds. The coupling is achieved through ions, the abundances of which are determined self-consistently from a chemical network. The ionic abundances control the time scale of cloud dynamic evolution via ambipolar diffusion. We have concentrated on a representative model in which an isolated magnetic cloud increases its central number density by a factor of $`10^3`$ to $`n_{\mathrm{H}_2}=10^6`$ cm-3 in 5.69 million years, with an eye on explaining the observational data on the chemical abundances of the well-studied starless core L1544. We find that the predicted peak values of CO, CCS, N2H+ and HCO+ column densities are within a factor of two or so of those inferred for L1544 from observations. The spatial distributions of these species are also consistent with those observed. An alternative, non-magnetic model was also considered for comparison. With the same initial conditions as in the magnetic model, the non-magnetic cloud collapses promptly, reaching the observed state in merely 1.56 million years. It produces a worse overall fit to the available data on L1544. In particular, the column density distribution of CCS predicted by the model is centrally peaked, which is not observed.

We conclude that our initial results of modeling lend some support to the standard picture of dense core formation out of strongly magnetized molecular clouds involving ambipolar diffusion over several dynamic times. There are, however, a number of uncertainties to which the results may be sensitive to, including the sticking probability and adsorption energies for various species. A parameter survey is needed to firm up the conclusion.

### 4.4 Future Refinements

There are several aspects of the coupled model that we wish to improve upon in the future. These include both dynamics and chemistry. Even though our spherical model captures the essence of the ambipolar diffusion-driven cloud evolution, the effects of magnetic tension (which tend to flatten a cloud) cannot be treated. Indeed, the opposite extreme, a disk-like geometry, may describe the observed elongated mass distribution in L1544 better. It should be straightforward to extend the model to this geometry, which would allow for the inclusion of rotation whose presence has been inferred for L1544 by Ohashi et al. (1999) based on CCS data. On the chemistry side, a major uncertainty is the adsorption energies for CO and other species. Different sets of adsorption energies have been considered by Aikawa et al. (2001), and a similar parameter study is needed for our model. Also uncertain are the probability of molecules sticking onto dust grains and grain surface processes, and an exploration of different cases is desirable. Another area of improvement would be the treatment of magnetic coupling coefficient, by including the effects of charged dust grains (especially small ones) and external UV radiation field.

The authors are grateful to P. Murphy and B. Turner for access to NRAO computing facilities, and to P. Caselli and the referee for helpful comments. This research is supported in part by RFBR Grant 01-02-16206.