text:
# Anomalous Phase Transition in Strained SrTiO3 Thin Films

Present Address: \]HASYLAB at DESY, Notkestr. 85, 22603 Hamburg, Germany

## Abstract

We have studied the cubic to tetragonal phase transition in epitaxial SrTiO3 films under various biaxial strain conditions using synchrotron X-ray diffraction. Measuring the superlattice peak associated with TiO6 octahedra rotation in the low temperature tetragonal phase indicates the presence of a phase transition whose critical temperature is a strong function of strain, with TC as much as 50K above the corresponding bulk temperature. Surprisingly, the lattice constants evolve smoothly through the transition with no indication of a phase change. This signals an important change in the nature of the phase transition due to the epitaxy strain and substrate clamping effect. The internal degrees of freedom (TiO6 rotations) have become uncoupled from the overall lattice shape.

preprint: L03-0901

Perovskite films have received a great deal of interest lately due to the potential for creating working technologies based on a variety of interesting properties such as high-TC superconductivity, colossal magneto-resistivity, ferro-electricity, and variable dielectric constants. These properties can be quite different in thin films versus nominally similar bulk samples. The primary reasons for the changes in properties are believed to be strain and defects. There is also an emerging theoretical effort to treat the effects of strain. To understand the mechanisms for these changes we must be able to conduct detailed microscopic measurements of the atomic and electronic structure.

SrTiO3 (STO) is a nearly ferroelectric material with a large dielectric nonlinearity and low dielectric loss at cryogenic temperature, making it ideal for tunable microwave devicesXi et al. (2000). STO is also a good model system for studying structural phase transitions. Bulk STO crystals are cubic at room temperature, with space group Pm3m(O$`{}_{}{}^{1}{}_{h}{}^{}`$), but become tetragonal, space group I4/mcm(D$`{}_{}{}^{18}{}_{4h}{}^{}`$), below about 105K. This phase transition involves the rotation of TiO6 octahedra and has been featured historically in the study of structural phase transitions as the classic example of a soft-mode phase transitionUnoki and Sakudo (1967). The unit cell of the tetragonal phase has a volume four times that of the cubic unit cell, with approximate unit cell of $`\sqrt{2}`$a$`\times `$$`\sqrt{2}`$a$`\times `$2a, where a is the lattice parameter of cubic unit cell. In this paper we describe the tetragonal phase as pseudo-cubic, in order to compare the structure before and after the phase transition. In this pseudo-cubic frame, the tetragonal phase has additional superlattice peaks at half integer index positions. These superlattice peaks disappear as the temperature is raised through the tetragonal-cubic phase transition.

In the bulk, the deviation from cubic symmetry is directly related to the rotation angle of the TiO6 octahedra. This rotation angle has been identified as the order parameter for this phase transitionM ller et al. (1968). Correspondingly, the intensities of the superlattice peaks are proportional to the square of the order parameter and can be used to track the phase transition.

For a second-order structural phase transition, usually the volume of the unit cell changes smoothly through the transition temperature, but lattice parameter versus temperature curves have a sudden change in slope. This is the case for bulk STOLytle (1964)Alefeld (1969). For epitaxial films the in-plane lattice parameters are subject to lateral restraint from the substrate and therefore do not have the freedom to change as in bulk. In this paper we show that for epitaxial films of STO the lattice parameters are indeed constrained though this does not completely inhibit the occurrence of the soft-mode phase transition. The strains do change both the transition temperature and the nature of the transition.

The films studied were STO films grown on (001) LaAlO3 (LAO) single crystal substrate with SrRuO3 (SRO) buffer layers. The samples were grown using the Pulsed Laser Deposition technique with details described elsewhereLi et al. (1998a). X-ray diffraction implies excellent epitaxy with average mosaics around 0.2 degrees and no detectable misoriented regions. Several different measurements of films grown in the same system under the same conditions have been reported elsewhereSirenko et al. (1999)Li et al. (1998b). All of the samples were of the type STO/SRO/LAO with variations of either the STO or SRO thickness. The main study consisted of STO films with thickness varying from 50 to 1000nm on a SRO-buffer layer 350nm thick. A secondary study investigated STO films 200nm thick on SRO buffer layers with thickness from 0 to 350nm. X-ray diffraction measurements were carried out at beamline X22A at the National Synchrotron Lights Source at Brookhaven National Laboratory. X22A has a bent Si(111) monochromator, giving a small beam spot and fixed incident photon energy of 10 keV. The longitudinal resolution with a Si(111) analyzer was at least 0.001Å-1 (HWHM) for an (0,0,2) peak, as measured from the LAO substrate. Below room temperature the sample was cooled in a closed cycle refrigerator with a temperature control of $`\pm `$0.5°.

We choose the coordinate system such that the axis normal to the film surface is the c axis. Like bulk STO, all of our STO films have extra half-integer superlattice peaks at low temperature. The intensity of these peaks decreases as the temperature rises and the peaks vanish at 130$``$170K for different samples. Thus the STO films also have a phase transition with similar symmetry change as seen in the bulk. Fig. 1 shows an example of such a superlattice peak for a \[1000nm STO / 350nm SRO\] sample. In the first panel are raw data scans through the (1/2, 1/2, 7/2) peak at some representative temperatures. The presence of a superlattice peak for SRO is a result of its structural phase transition occuring at $``$800K. The STO and SRO peaks are well resolved and the STO peak intensity decreases with increasing temperature. Fig. 1 panel (b) is a summary plot of the peak intensity versus temperature. Compared to the bulk, the intensity versus temperature curve is shifted towards higher temperature. The intensity does not fall abruptly to zero but includes a tail that extends to higher temperature. This indicates that the phase transition is rounded and occurs over a range of temperature. This rounded transition is most likely due to the presence of a range of local strains within the film. Due to the rounding of the transition, we need to derive a method for extracting the average TC. We extrapolate the linear portion of the curve between 10% to 40% maximum intensity and define the transition temperature as where this line crosses the zero intensity. The use of a linear function is justified by the data itself as well as previous work that has found a power law coefficient b very close to one for this phase transition in bulk samplesDoi and Takahashi (2000). While this method may produce TC too high by a few degrees ($`<`$5K), it allows for an accurate comparison of TC between different films. We believe we can determine $`\mathrm{\Delta }`$TC to about 0.5K.

For samples with a 350nm SRO buffer layer and varying STO thickness, TC increases monotonically with in-plane strain as shown in Fig. 2(a), and in more detail in Fig. 2(b). The difference, $`\mathrm{\Delta }`$TC, is about 14K (or $`\mathrm{\Delta }`$TC/TC$``$10%) between the 1000nm STO sample and 50nm sample although the in-plane strain at room temperature measured by glancing incidence x-ray diffraction (not shown) only changes from 0.01% to 0.23%. The thinner film, with greater strain, has the higher TC, while thicker film has TC close to the bulk value. Thus the structural phase transition is very sensitive to the biaxial strain caused by the lattice mismatch at the interface, with the ratio ($`\mathrm{\Delta }`$TC/TC)/$`\epsilon `$ $``$45, where $`\epsilon `$ is the in-plane strain.

We also studied a series of samples with the same 200nm STO layer but varying the thickness of SRO buffer layer. Here we observed even larger effect of strain on TC, as shown in Fig. 2(c). In this case the largest strains and the largest shifts in TC appear for SRO buffer layers of intermediate thickness. We believe this is due to the complexity of the SRO domain structure. The highest TC is in the sample with a 20nm SRO buffer layer. In this case, TC is shifted about 30K above the STO film on a 350nm buffer and about 40K above the most bulk like of our films \[1000nm STO / 350nm SRO\]. Clearly the overall rise in TC can be quite large. Further detail will be given in a follow up paperHe et al. (a), that also includes a discussion of the temperature-strain phase diagram which has been theoretically predictedPertsev et al. (2000).

Although from the presence of the superlattice peaks we can see that there is a phase transition near 130K, the out-of-plane lattice constants evolve smoothly through the transition temperature with no indication of the transition. This is illustrated in Fig. 3 for the 200nm STO / 350nm SRO film. The variation in lattice constant is merely the thermal expansion with no indication of a phase transition. This is obviously different from the bulk STO where a sudden slope change is clearly seen at 105KAlefeld (1969)Hirota et al. (1995). We ascribe this effect to the clamping effect of the substrate upon the film. The film in-plane lattice parameters are not free to move to their bulk equilibrium values because they are tied to the substrate lattice. Our experiments show that the overall film is tetragonal at all temperatures measured. However this doesn’t prevent the TiO6 octahedra from rotating and ordering leading to a reduced symmetry phase. In fact the tetragonal cell shape seems to favor the rotation of TiO6 octahedra, stabilizing the lower symmetry phase and enhancing TC.

It is a remarkable that the phase transition in the STO films involves a lowering of the internal symmetry due to TiO6 octahedra rotations without a change in the unit cell size and shape. This unique phase transition manifests a decoupling of the internal degree of freedom of TiO6 rotation and the external degree of freedom of the unit cell dimensions. Such a difference between the nature of structural phase transitions in epitaxial films and bulk materials may occur in many situations and should be considered in analyzing a variety of results.

Raman scattering experiments on similar films have studied the temperature dependence of the structural R modes, which are Raman active only below TCAkimov et al. (2000). In bulk STO these phonons disappear at 105K, while in thin films their intensity vanishes at about 120$`\pm `$5K. The mode is detectable in Raman upon the doubling of the unit cell thus folding this mode back to the zone center. This change is dependent upon the rotation of the octahedra and not the presence of a tetragonal unit cell. In contrast, dielectric measurements show obvious thickness dependence and the temperature range over which the dielectric constant is tunable by an electric field is vastly increased versus the bulkLi et al. (1998b). The empirical evidence suggests that this behavior is connected to the presence of a tetragonal unit cell.

A review of the literature reveals that a similar effect may exist in BaTiO3 (BTO) films. Bulk BTO has a series of phase transitions at 183K, 278K and 403K, at each of which the lattice constants change drastically. All of the lower temperature phases are ferroelectric. In contrast, lattice parameter measurements of epitaxial films show no abrupt changes in either the a or c constant, which both increase roughly linearly with temperature from 15K to 800KTerauchi et al. (1992). However, it appears from electrical measurements that the ferroelectric phase transition still occurs in similar BTO films on a Pt/MgO substrateYoneda et al. (1993)Yoneda et al. (2001). Ferroelectricity is due to the shift of position of the Ti ion with respect to the O octahedron surrounding it, not the cubic to tetragonal change in the unit cell. Thus ferroelectricity is due to the internal degree of freedom of the phase transition. The ferroelectric transition in BTO films seems to be another example of an epitaxy strain and substrate clamping effect that permits a phase transition in the internal atomic arrangements (ferroelectricity) but does not allow external dimension changes.

It appears that the decoupling of the internal and external aspects of a phase transition in films is not universal. One example is NdNiO3. Both bulk and films of this material have a metal-insulator transition around 200K. Recent x-ray diffraction measurements on epitaxial films of NdNiO3 on LaAlO3 have shown that the out of plane lattice constant exhibits a sudden change at the transitionHe et al. (b). Further experiments are ongoing to try to determine what makes the two transitions fundamentally different with respect to the clamping effect of the substrate and the relaxation of the overall unit cell dimensions.

We acknowledge S. P. Alpay for helpful discussion. BW thanks the A.P. Sloan Foundation Research Fellow program for partial support of this work. Work at Brookhaven is supported by Division of Material Sciences, U.S. Department of Energy under contract DE-AC02-98CH10886. The work at Penn State is supported by NSF under grant No. DMR-9702632, and work at UConn is supported by the University of Connecticut Research Foundation.


format: text

text:
# 1 Introduction

## 1 Introduction

We present in this note an update of the predictions for the top quark production cross-section at the Tevatron. These predictions are based on two complementary ingredients:

1. the evaluation of the parton-level cross-sections, carried out in perturbative QCD with the inclusion of the full next-to-leading-order (NLO) matrix elements , possibly improved with the resummation to all orders of perturbation theory of classes of large soft logarithms
2. the proton parton densities (PDFs), which are typically extracted comparing existing data with NLO calculations available for the relevant processes, and extrapolated to the relevant region of $`Q^2`$ using the NLO evolution equations (more recently, accurate estimates of the exact NNLO results have also become available , based on partial evaluations of the three-loop splitting functions).

The numbers we present here are based on the theoretical framework introduced in  and , where the complete NLO calculation of the $`t\overline{t}`$ cross-section was improved with the resummation of leading  and next-to-leading  soft logarithms appearing at all orders of perturbation theory. The introduction of resummation turns out to have only a mild impact on the overall rates (the effects at NLL are typically of the order of a few percent), but improves the stability of the predictions with respect to changes of the renormalization scales. While no progress has occurred since 1998 in the calculation itself, significant development has taken place in the determination of the PDFs. In addition to much improved data from HERA, from fixed-target DIS experiments at FNAL, and to the implementation of Tevatron jet and $`W`$ production data in the fits, progress has occurred in the assessment of the true uncertainties associated with the global fits to these data. This work, which recently received considerable attention (Giele, Keller and Kosower , CTEQ , MRST , Botje , Alekhin ), has led to sets of PDF parameterizations which should provide a meaningful estimate of the “1-$`\sigma `$” uncertainty deriving from PDFs to be associated to any calculations of hard processes in hadronic collisions.

The introduction of these PDF sets “with uncertainties” relaxes the much constrained predictions which used to be anchored to predefined functional parametrizations, and it is natural to anticipate that the range of predictions for a given hard cross-section will be increased.

## 2 Outline of the uncertainty estimate

We shortly outline here the details of our calculation, before presenting the numerical results. Unless explicitly denoted as $`\sigma _{NLO}`$, all of our results are obtained using the NLL-improved formalism of ref. .

### 2.1 Scale uncertainty

The evaluation of the purely theoretical uncertainty is based on the standard exploration of the cross-section dependence on the renormalization ($`\mu _R`$) and factorization ($`\mu _F`$) scales used in the perturbative calculation. In this work, we follow the standard convention of considering the range $`m_{top}/2<\mu <2m_{top}`$, setting $`\mu _R=\mu _F\mu `$. A justification for this choice can be found in , where it was shown that $`\mu m_{top}/2`$ corresponds to a point of minimal sensitivity, providing a maximum of the cross-section in the range $`0.1<\mu /m_{top}<10`$. In the range of mass consistent with the current data, and for the two CM energy values of run I and run II (1.8 and 1.96 TeV, respectively), the relative scale uncertainty at NLO is of the order of $`\pm 10`$%, independent to good approximation of $`\sqrt{S}`$, $`m_{top}`$ and PDF sets. In this region of parameters, the maximum value is obtained for $`\mu m_{top}/2`$, and the minimum for $`\mu =2m_{top}`$. The inclusion of NLL resummation corrections reduces the uncertainty to the level of approximately $`\pm 5`$222This number, as well as all numerical estimates presented in this document, correspond to the choice $`A=2`$, where $`A`$ is the parameter introduced in  to parameterize the uncertainty about subleading higher order terms. In that paper, it was found that $`A=2`$ gives a better estimate of the higher order uncertainties. $`A=0`$, for example, would reduce the scale dependence to only $`\pm 2.5`$%, without changing significantly the central value of the resummed cross-section. This is the effect of very small NLL corrections to the NLO result for small values of $`\mu `$, where the NLO rate is largest, and bigger corrections for large $`\mu `$.

For completeness, we also considered the possibility of varying independently the value of renormalization and factorization scale. These were chosen in the range $`0.5<\mu _R/\mu _F<2`$, with $`0.5<\mu _{R,F}/m_{top}<2`$. We verified (see later) that within this range the results obtained using the choice $`\mu _R=\mu _F`$ are not altered significantly, leading only to a small increase of the upper estimate.

### 2.2 PDF uncertainty

In the framework of , PDFs with uncertainties come in sets of $`n_{PDF}`$ pairs, where $`n_{PDF}`$ is the number of parameters used in the fits. Each pair corresponds to the fit obtained by varying of $`\pm 1\sigma `$ the value of the fit parameter eigenvalues, after diagonalization of the correlation matrix. By construction, the systematic uncertainty obtained for the observable $`𝒪`$ is given by:

$$\mathrm{\Delta }𝒪=\frac{1}{2}\sqrt{\underset{i=1,n_{PDF}}{}(𝒪_{i+}𝒪_i)^2}$$
(1)
where $`𝒪_{i\pm }`$ is the value obtained using the PDF set corresponding to the variation of the $`i`$th eigenvalue within its error range. The central value of the prediciton is obtained using a reference PDF set, typically labelled with $`i=0`$. We explore in this work the sets in the CTEQ6  parameterizations ($`n_{CTEQ}=20`$, corresponding to 40 sets, plus 1 reference set) and in the MRST 2002  compilation ($`n_{MRST}=15`$, corresponding to 30 sets, plus 1 reference set). All sets in the CTEQ compilation have $`\alpha _s(M_Z)`$ =0.118, while those in the MRST one have $`\alpha _s(M_Z)`$ =0.119. The CTEQ sets are labeled as follows: 6M for the default set, and 101-140 for the 20 $`\pm 1\sigma `$ variations. The MRST sets are labeled as 0 for the reference set, and 1-30 for the 15 $`\pm 1\sigma `$ variations. In both cases, CTEQ and MRST, we use the default values of tolerances chosen by the two groups to best represent the uncertainty. In particular, CTEQ selects $`\mathrm{\Delta }\chi ^2=100`$, while MRST selects $`\mathrm{\Delta }\chi ^2=50`$.

In addition, we shall also consider three sets obtained by the MRST group in 2001 , where the values of $`\alpha _s`$ was frozen to $`\pm 1\sigma `$ from the central world average. We shall label these sets as A01L for the low-$`\alpha _s`$ ($`\alpha _s(M_Z)`$ =0.117) fit , A01H for the high-$`\alpha _s`$ ($`\alpha _s(M_Z)`$ =0.121) fit , J01 for a fit based on Tevatron jet data ($`\alpha _s(M_Z)`$ =0.121).

## 3 Results

Table 1 summarizes the results obtained with the PDF sets used in 1998, when the work in ref.  appeared. The numbers agree with what appears in Table 1 of that document.

Table 2 gives the central value and error for the CTEQ sets, for three values of the top mass (170, 175 and 180 GeV) and the two CM energies of interest ($`\sqrt{S}=1800`$ and 1960 GeV). We list the results obtained at the three reference values of the mass scale $`r_\mu =\mu /m_{top}=0.5,\mathrm{\hspace{0.17em}1},\mathrm{\hspace{0.17em}2}`$. Table 3 provides the same information for the MRST sets.

Figure 1 shows the contour plots of the NLL cross-section when $`\mu _R`$ and $`\mu _F`$ are varied independently. The region defined by the oblique solid lines corresponds to $`0.5<\mu _R/\mu _F<2`$. It shows that within this domain the range of NLL rates is compatible with the range obtained using $`\mu _R=\mu _F`$.

In principle one should combine the uncertainty due to PDFs and that due to the scale choice in quadrature. We prefer to add them linearly, since the scale uncertainty is not really a systematic error in the strict sense. We therefore quote our range for the top cross-section as

$$\sigma (r_\mu =2)\mathrm{\Delta }\sigma _{PDF}(r_\mu =2)<\sigma <\sigma (r_\mu =1/2)+\mathrm{\Delta }\sigma _{PDF}(r_\mu =1/2)$$
(2)
The corresponding values are given in Table 4. The similar results for the MRST compilation are provided in Table 5.

Three comments are in order:

1. the uncertainty ranges obtained using the CTEQ sets, for a fixed choice of scale, are almost twice as large as those for the MRST sets. We understand this is the result of the different tolerance criteria used by the two groups (see Appendix B4 of  and Section 6 of  for some discussion). The MRST range increases however if we include in the analsys the 2001 sets with varying $`\alpha _s`$ . This is shown in Table 6. In this case the lowest predictions are obtained from the 2001 A01L fit, with the low value of $`\alpha _s`$ , while the highest prediction comes from the 2001 jet-based J01 fit. After the $`\alpha _s`$ variation is included, the MRST range becomes compatible with that of CTEQ’s.
2. the central values obtained today for the top cross-section are about 3% larger than those obtained in 1998. At $`\sqrt{S}=1.8`$ TeV and $`\mu =m_{top}=175`$ GeV we had 5.06 pb with the set MRSR2 ($`\alpha _s(M_Z)`$ =0.119). We now have 5.19 pb with CTEQ6M, and 5.21 pb with MRST0.
3. the contribution of the PDF systematics to the uncertainty range is large. In the case of the CTEQ sets, it is of the order of 6-7%, larger than that due to the choice of scale. This is a result of the large sensitivity of the top cross-section to the large-$`x`$ gluon content of the proton, which is still poorly known. For CTEQ the largest contribution to the error comes from the two sets 129 and 130333This is consistent with what found in a recent study of jet produciton at the Tevatron .. For these two sets, we find the contribution of the $`gg`$ channel to be respectively 11% and 21% of the total rate. For comparison, the contributions of the $`q\overline{q}`$ production channel for sets 129 and 130 are the same to within 1%. In other words, the PDF uncertainty on the top rate is mostly driven by the poorly known gluon density, whose luminosity in this kinematic range varies by up to a factor of 2 within the 1-$`\sigma `$ PDF range.

While the overall production rate has a large relative uncertainty of approximately $`\pm 15`$%, it is important to point out that the ratio of cross-sections at $`\sqrt{S}=1.96`$ TeV and $`\sqrt{S}=1.8`$ TeV is extremely stable. In the case of the CTEQ sets, for example, we found $`\sigma (1.96)/\sigma (1.8)=1.295\pm 0.015`$ after scanning over the set of scale choices and for $`170<m_{top}<180`$ GeV. The error is about 1%. We therefore consider the prediction of the relative cross-section at the two energies to be a very stable one.

For reference, we collect the full set of cross-sections (at $`\sqrt{S}=1.96`$ TeV and $`m_{top}=175`$ GeV) for all CTEQ sets and scale choices in Table 8. Here, for the sake of documentation, we provide the NLO rates and the NLL-improved ones separately.

## 4 Conclusions

We reiterate here the main findings of this study. The inclusion of the full PDF systematics, made possible by the recent works of several groups, leads to a more realistic estimate of the top cross-section uncertainty. The latest MRST and CTEQ sets give rise to cross-sections which are typically 3% larger than what obtained with sets available at the time of Run I. In addition to the increase in rate, the size of the uncertainty range has also increased, to a value of the order of $`\pm 15`$%, dominated by the PDF and $`\alpha _s`$ uncertainties. The leading source of PDF uncertainty comes from the (lack of) knowledge of the gluon luminosity at large values of $`x`$. The $`gg`$ contribution can in fact change through the PDF sets by up to a factor of 2 (from 10% to 20% of the total rate at 1.96 TeV). We find that the MRST sets give rise to a smaller PDF uncertainty, a result we ascribe to the tighter tolerances required by MRST in defining the range of the eigenvalues. The MRST uncertainty increases however to values consistent with CTEQ’s once the sets obtained from a $`\pm 1\sigma `$ change of $`\alpha _s(M_Z)`$ are included. This underscores the importance of including the $`\alpha _s`$ uncertainty into the PDF fits in a more systematic fashion. On the same footing, the impact of higher order corrections, as well as of the treatment of higher twist effects in the fitting of low-$`Q^2`$ data, may need some more study before a final tabulation of the PDF uncertainties is achieved .

We collect in Table 7 our final results. This summary table includes the CTEQ6M set and $`\mu =m_{top}`$ as central values, and the most extreme rates extracted from Tables 4, 5 and 6 as lower (with $`\mu =2m_{top}`$) and upper values (with $`\mu =m_{top}/2`$).

In spite of the overall large uncertainty, the ratio of cross-sections at 1.96 and 1.8 TeV is extremely stable, being equal to $`1.295\pm 0.015`$ over the mass range $`170<m_{top}<180`$ GeV.

## Acknowledgements

We thank S. Catani and J. Huston for useful comments and discussions on the topic of this work.


text:
# Baryon Number Transfer in Hadron+Nucleus and Nucleus+Nucleus Collisions: a Link Between Elementary and Complex Interactions

## 1 Introduction

The NA49 experiment at the CERN SPS is a large acceptance hadron spectrometer, based on Time Projection Chambers. . The detector provides precision tracking and particle identification via the measurement of the ionization energy loss in the TPCs. The system is completed by Time of Flight detectors and calorimeters, the latter providing neutron detection above $`x_F=0.2`$ in proton induced reactions. A dedicated counter is designed to control the centrality in p+A interactions, which measures the multiplicity of the recoil (grey) protons.

## 2 Net proton distribution in the elementary collisions

The net proton yield in a hadronic interaction is defined by $`p\overline{p}`$, so antiprotons are subtracted from the protons to cancel the pairproduced component. The net proton distribution (see Fig.1) in p+p interaction is forward-backward symmetric, to first order flat with a diffractive peak. Above $`x_F>0.2`$ the Feynman-scaling holds, which predicts the invariant cross-section as a function of $`x_F`$ being independent of $`\sqrt{s}`$ at high energies. Intuitively, this symmetric distribution has two components: the forward region associated with the projectile, and the backward region with the target. In the following sections, the experimental study of this intuitive picture will be discussed.

## 3 Two component picture of the elementary collisions

The basic idea of the two component picture  is shown on a cartoon (Fig.2): the net proton distribution in p+p is built up from two symmetric components. Assuming, that *the target component is independent of the projectile type of the interaction*, the two components can be experimentally separated. For the measurement, a baryon free projectile was used ($`\pi `$+p) which allows the measurement of the target component alone (see Fig.2).

As $`\pi ^+`$ and $`\pi ^{}`$ has non-zero isospin, the produced net proton on the projectile side will also be non-zero (only $`p\overline{p}+n\overline{n}`$ would vanish). To cancel this effect, the average of $`\pi ^++p`$ and $`\pi ^{}+p`$ collisions were used for the measurement. The net proton distribution in this averaged $`\pi +p`$ is presented in Fig.1

One simple prediction of the two component picture is that the yield of the net protons in $`\pi +p`$ at $`x_F=0`$, so at the center of mass of the colliding system, has to be half of the value measured in $`p+p`$. This prediction is reasonably fulfilled. In the forward region, subtracting the $`\pi +p`$ measurement from the $`p+p`$ values, the *projectile component of the net proton spectrum in p+p* can be obtained, shown in Fig.3.

## 4 Net proton distribution in hadron+nucleus interactions

The net proton distribution in p+Pb at different centralities is shown in Fig 4, compared to net protons in p+p. The centrality of the p+Pb collision is measured using a dedicated detector which determines the multiplicity of the knockout protons from the nucleus. From this multiplicity, applying Glauber-model calculations and detailed detector simulation, the number of collisions ($`\nu `$) is estimated. The distributions steepen up at higher centrality due to the baryon stopping phenomenon .

## 5 Two component picture of hadron+nucleus interactions

The proposed two component picture is more complicated for p+A than for p+p. The following assumption can be made: to first order, on the target side each collision of the projectile produces the target component of the net proton distribution in p+p, i.e. *the target component in pA piles up in proportion with the number of collisions*. The projectile however, as it suffers multiple collisions, develops a *stopped projectile component*. A cartoon describing the picture is shown on Fig.5.

In order to determine the projectile component of the net proton distribution in p+A, the idea of baryon free projectile can be used again. According to the prediction of the two component picture, the target component measured by $`\pi +A`$ scales up from $`\pi +p`$ with the number of collisions; Fig.6 demonstrates this effect.

Subtracting the net proton distribution measured in $`\pi `$+Pb from the net protons in p+Pb at the same centrality, the projectile component of the net proton distribution can be determined in p+Pb (see Fig.7). A smooth progression of baryon number transfer is observed as a function of the number of collisions.

## 6 Two component picture of nucleus+nucleus collisions

In case of AA, the baryon-free projectile concept is not applicable. The only constraint which provides an approximation of the projectile component of the net proton spectrum is the following: at $`x_F=0`$, the value of the projectile component is half of the measured spectra, and the target component dies out fast as we move forward. In Fig.8 the measured net proton distribution and the estimated projectile component is presented.

The conclusion of the studies above is that with increasing centrality, the projectile component of the net proton distribution have a similar behaviour in p+A and A+A interactions, which is a smooth transition of baryon number transfer from the elementary to more complex interactions.

## 7 Reference (isospin corrected) net proton for Pb+Pb collisions

In the Pb nucleus, the n to p ratio is about 60 to 40, so the net proton distribution in a reference N+N collision has to be determined to take care of the neutron content; here N=0.4p+0.6n. Within the framework of the two component picture the isospin symmetry holds separately for the two hemispheres, so the net forward proton in n+p or n+n is equal to the net neutron in p+p. Relying on this assumption, the net proton distribution in an isospin averaged N+N collision can be obtained; this is shown in Fig.9.

## 8 Correlation of the pion multiplicity with the final state baryon in p+p interaction

In case of the p+A and A+A collisions, the baryon number transfer is characterised by the centrality. In the elementary p+p collision the centrality can not be reliably defined, so the study must be done in the other way round: Fig.10 shows that the position of the final state baryon characterises the pion multiplicity p+p events.

The measurement confirms the intuitive expectation: in those events, where the final state baryon is fast (especially in the diffractive region), the multiplicity is low. If the final state baryon is more stopped, more pions are produced. The correlation is independent of the baryon isospin.

## 9 Prediction of the pion multiplicity in Pb+Pb collisions from correlation measurements in p+p

Based on the phenomenon of baryon stopping in central p+A and A+A collisions, and the strong correlation of p+p observables with the final state baryon, the following assumption can be made: *the baryon number transfer characterises the inelasticity (centrality) in p+p, p+A and A+A interactions.* This assumption provides a method to compare minimum bias p+p collisions with central A+A collisions, where the trigger bias is to be taken into account.

In a p+p collision, any quantity, e.g. the pion multiplicity measured as a function of $`x_F`$ of the final state proton defines a correlation function $`C(x)`$. The pion multiplicity in minimum bias p+p then can be calculated by weighing $`C(x)`$ with the net projectile proton distribution; this trivial statement can be expressed in the following way:

$$<\pi >_{pp}=\frac{C(x)P_{pp}𝑑x}{P_{pp}}$$

where $`P_{pp}`$ is the net projectile proton distribution in p+p as a function of $`x_F`$. In A+A, the stopped net projectile proton distribution can be viewed as protons from a superposition of p+p (or N+N) interactions, where events with more stopped final state protons have larger weight; this means that the pion multiplicity in A+A (per participating nucleon) is predicted to be

$$<\pi >_{PbPb}=\frac{C(x)P_{PbPb}𝑑x}{P_{PbPb}𝑑x}$$

In this case, the net proton distribution in Pb+Pb ($`P_{PbPb}`$) is to be applied. The result of this prediction is presented on Fig.11 in comparison with pion multiplicity measurements from PbPb interactions at different centralities.

The prediction gives a very good qualitative descriprion of the measurements. If the assumptions are correct, this means that the full increase of pion multiplicity in central Pb+Pb interaction can be explained with a superposition of elementary collisions, taking proper care of the centrality selection.

## 10 Conclusions

An experimental study of the two component picture of hadronic interactions at $`\sqrt{s}`$=17.2 GeV was discussed. In this picture, the net proton distribution is built up by a target and a projectile component; these components can be separated experimentally with a baryon free (pion) projectile. The projectile component of the net proton distribution in p+A and A+A collisions is found to behave similarly at increasing centrality, namely that the baryon number moves towards the center. Based on the assumption that the baryon transfer characterises the inelasticity of p+p, p+A and A+A collisions, a model independent procedure was developed to compare minimum bias p+p and central Pb+Pb collisions, which showed that the increase of the pion multiplicity in A+A compared to p+p interactions can be fully predicted only by taking care of the sizeable baryon stopping in A+A.

## Acknowledgement

This work was supported by the Director, Office of Energy Research, Division of Nuclear Physics of the Office of High Energy and Nuclear Physics of the US Department of Energy (DE-ACO3-76SFOOO98 and DE-FG02-91ER40609), the US National Science Foundation, the Bundesministerium fur Bildung und Forschung, Germany, the Alexander von Humboldt Foundation, the UK Engineering and Physical Sciences Research Council, the Polish State Committee for Scientific Research (5 P03B 13820 and 2 P03B 02418), the Hungarian Scientific Research Foundation (T032648, T14920, T32293 and F034707), the EC Marie Curie Foundation, and the Polish-German Foundation.

## Notes

* Permanent address: Eötvös Loránd University, Department of Atomic Physics, Budapest, Hungary;
E-mail: dezso.varga@cern.ch


text:
# Angular Momentum Evolution of Young Stars: Toward a Synthesis of Observations, Theory, and Modeling1footnote 11footnote 1The conference was held in Albuquerque, New Mexico, 2002 June 5, as a Topical Session at the 200th Meeting of the American Astronomical Society.

## 1 Overview of Angular Momentum Evolution of Young Stars—Donald Terndrup (Ohio State University)

During the last decade we have witnessed a tremendous growth in the amount and quality of data on stellar rotation, which now includes samples in ages from T Tauri stars to clusters as old as the Hyades. Along with this has been an increase in the sophistication of theoretical models of the spindown of stars before and on the main sequence. Many of us in the field now believe we have a good handle on the three principal ingredients of these models:

$``$ The distribution of initial rotation rates in the pre-main sequence, coupled with the (possible) influence of accretion disks.

$``$ Angular momentum loss from a magnetized solar-like wind, which operates on time scales of tens to hundreds of Myr.

$``$ Prescriptions for the angular momentum distribution in convective regions of the star and angular momentum transport in radiative regions.

The initial conditions for rotation are set during the star formation process. T Tauri stars are logical starting points, as they represent the earliest easily-observable phases of evolution at which stars have reached their final mass. These stars are at the beginning of the hydrostatic phase of evolution prior to the ignition of nuclear reactions, and typically have modest accretion disks of order 0.01 M. There is extensive data on rotation rates of T Tauri stars, beginning with the work of Vogel & Kuhi (1981), which shows that there is already a significant range of rotation rates among protostars with ages of 1 Myr or less. Recent surveys in the Orion Nebula Cluster (ONC), which has an age of 1–3 Myr, have provided hundreds of rotation periods as will be reviewed below.

The observation that the youngest stars rotate at a small fraction (10–20%) of the breakup velocity indicates that they must have had an efficient method for shedding angular momentum. The most common approach to modeling this is to invoke “disk-locking,” as introduced by Königl (1991). In the presence of an accretion disk, a star could maintain a nearly constant surface rotation rate as it contracts (Shu et al. 1994); this corotation between the star and the accretion disk would govern the stellar rotation rate. Theoretical models of angular momentum evolution have therefore included the accretion disk lifetime as a parameter describing the initial conditions (Sills, Pinsonneault, & Terndrup 2000 and references therein).

Later on, when stars have arrived on the main sequence, the stars experience continuing angular momentum loss, envisioned as resulting from a magnetized stellar wind (Weber & Davis 1967). In all young clusters, most stars are rotating slowly, with $`v\mathrm{sin}i<20`$ km s-1, but a minority are rotating must faster. The loss mechanism efficiently removes stellar angular momentum: in clusters such as Alpha Persei (age 50 – 80 Myr), the fastest rotators have equatorial velocities near 200 km s-1, but these have fallen to below 10 km s-1 by the age of the Hyades (Radick et al. 1987). The existence of rapid rotation at ages of $`100`$ Myr suggests that the loss rate is saturated or suppressed at high rotation rates. Observations of stellar activity (e.g., Krishnamurthi et al. 1998) and a variety of theoretical arguments lead to the idea of a saturation threshold (Kawaler 1988). In this picture the loss rate scales as $`\omega ^3`$ for slow rotators (Skumanich 1972), where $`\omega `$ is the angular velocity, but scales as $`\omega _{\mathrm{crit}}^2\omega `$ for $`\omega >\omega _{\mathrm{crit}}`$. The saturation threshold depends on stellar mass, in the sense that it has a lower value for stars of lower mass (Krishnamurthi et al. 1997).

The internal distribution of angular momentum within a star must also be understood in a complete theory of angular momentum evolution. Low-mass stars are initially fully convective. Their surface rotation rates can be used to infer an initial angular momentum if the angular momentum distribution in convective regions is specified, since the convective overturn timescale is much shorter than the evolutionary timescale. Once on the main sequence, however, stars with masses above 0.35 solar masses develop radiative cores. Stars with higher masses have much thinner convective zones: at 1 solar mass, for example, only 10% of the moment of inertia is in the convective envelope. When a radiative core develops, the internal transport of angular momentum in the radiative regions must be accounted for, since the timescales for evolution and for angular momentum transport are comparable. All this means that the models have to begin with assumptions about the (radial) angular momentum profile; obvious limiting cases are to assume initial solid body rotation or constant specific angular momentum. The models also need to have an adjustable efficiency of internal angular momentum transport.

## 2 Rotation of Solar-Mass Stars at Different Stages of Evolution: Observations

### 2.1 Rotation of Protostars—Thomas Greene (NASA Ames) & Charles Lada (CfA)

It is important to determine the angular momenta of embedded protostars to understand better the star-formation accretion process and to establish the angular momentum “zero point” from which stars evolve. Unfortunately, the angular momenta of protostars are largely unknown. Very high extinctions ($`A_V40`$ mag) make them impossible to observe in visible light with high-resolution spectrographs. Further, these objects are not very bright in the near-IR ($`K10`$), and they have very large continuum veilings ($`r3`$) at these wavelengths, greatly diluting the strength of their photospheric absorption lines and star spots.

Despite these difficulties, we have recently obtained high-resolution absorption spectra of some protostars using sensitive IR spectrographs on large telescopes. These data show that protostars with flat-spectrum and Class I spectral energy distributions (SEDs) have spectral types, surface gravities, and radii similar to more evolved T Tauri stars (TTS). However, the embedded protostars observed to date have considerably different $`v\mathrm{sin}i`$ and angular momenta. These objects appear to be rotating significantly faster than most TTS; typical protostars have $`v\mathrm{sin}i40`$ km s-1 (Greene & Lada 1997, 2002). Recent ASCA and Chandra X-ray observations also suggest that these protostars have short rotation periods.

Our Keck NIRSPEC observations of YLW 15 have provided the first measurement of the angular momentum and detailed astrophysics of a Class I protostar (Greene & Lada 2002). With a projected rotation velocity $`v\mathrm{sin}i=50`$ km s-1, its co-rotation radius is about 2 $`R_{}`$. This is about equal to the magnetic coupling radius predicted by the Shu et al. (1994) X-wind/magnetospheric-accretion model for the observed mass accretion rate of $`\dot{M}=1.6\times 10^6M_{}`$ yr-1. Like other Class I protostars, YLW 15 will evolve into a TTS in $`10^5`$ yr. Since TTS in the same $`\rho `$ Oph cloud rotate only about one-third as fast (e.g., Greene & Lada 1997), YLW 15 (and other similar protostars) must shed approximately two-thirds of its angular momentum in about this amount of time. Within the framework of magnetic star-disk interaction, this is possible provided that its mass accretion rate drops quickly to $`10^7M_{}`$ yr-1.

These observations are in their infancy, but they may prove essential for understanding initial stellar angular momentum and star-disk interactions during the protostellar evolutionary phase. Soon we expect to extract statistically significant results on these issues from our recently-completed Keck NIRSPEC observations of a few dozen embedded protostars in nearby star-forming clouds.

### 2.2 Rotation of Pre–Main-Sequence Stars—William Herbst (Wesleyan University)

Low-mass stars often have large, stable, non-axisymmetric, cool spots on their surfaces which makes it possible to measure their rotation periods by photometric monitoring. There is nothing we can determine with as much accuracy for PMS stars as their rotation rates.

Angular velocities accurate to about 1% are currently available for $`1500`$ PMS stars (see Fig. 2). Approximately 40% of the stars with $`12.5<I<15.5`$ in the Orion Nebula Cluster (ONC), the best studied example, now have rotation periods measured for them. There is no reason to believe that this sample is significantly biased with respect to rotation properties (Rhode et al. 2001; Herbst et al. 2002) or other characteristics, so it provides an excellent means of studying rotation among very young stars. We focus our discussion on angular velocity instead of angular momentum, because the latter requires knowledge of the radii of stars, which are particularly uncertain at PMS ages.

The first surprising fact about PMS rotation rates is the breadth of the distribution—it extends over at least a factor of 30 from $`2/3`$ day to $`20`$ days. The second surprising fact is that it is strongly mass dependent. Stars in the range of $`1/4`$ to 1 $`M_{}`$ have a bimodal distribution (Fig. 3) with peaks near 2 d and 8 d, as first shown by Attridge & Herbst (1992) and Choi & Herbst (1996), while lower mass stars have a unimodal distribution and generally spin faster (Herbst et al. 2001). This behavior is not limited to the ONC, but is seen also in NGC 2264 (Lamm et al. 2002). As NGC 2264 is a somewhat older cluster, the peaks in the higher mass distribution are at 1 d and 4 d, instead of 2 d and 8 d as in the ONC. This is evidence for spin-up of a significant percentage of the PMS stars in NGC 2264. Assuming that it is three times the age of the ONC (i.e. 3 My as opposed to 1 My), the amount of spin-up (i.e. a factor of 2) is consistent with expectation based on PMS contraction models and conservation of angular momentum.

It is interesting to speculate about how these (mass dependent) distributions could have arisen during the early PMS or proto-star stages. The leading hypothesis is disk-locking. One possibly relevant piece of information is that, in the ONC, we find a distinct correlation between rotation rate and strength of near-IR excess (indicative of the presence of an inner circumstellar disk). In particular, the $`IK`$ excess emission is $`0.55\pm 0.05`$ for slowly rotating stars (i.e. those with angular velocities less than 1, corresponding to rotation periods longer than 6.3 days), while it is $`0.17\pm 0.05`$ for rapid rotators (with $`\omega >2`$ corresponding to $`P<3.1`$ days). The first order interpretation of this is that disks tend to brake stars, as the disk-locking theory predicts. However, under some scenarios one expects a correlation in the opposite sense, as Stassun et al. (2001) have discussed. A deeper understanding of the meaning of the near-IR correlation and its relation to disk-locking theory will require further study, particularly aimed at establishing the presence and nature of disks around PMS stars with known rotation rates.

### 2.3 Rotation of solar-type stars at zero age and beyond—Sydney A. Barnes (University of Wisconsin, Madison)

The most basic and most solidly known result in stellar rotation is that of the solar rotation rate, on the surface, and more recently, in the interior. We know that upper main sequence (F and earlier) stars spin fast, and that lower main sequence (G, K & M) stars spin slowly. This is related respectively to the absence and presence of surface convection zones and winds on these stars (Schatzman 1962).

For lower main sequence stars, the spindown is described by the Skumanich (1972) relationship: $`vt^{1/2}`$, where $`v`$ and $`t`$ are the stellar rotation velocity and the age of the star respectively. This suggests that on the main sequence, the rate of angular momentum ($`J`$) loss is related to the angular velocity ($`\omega `$) by $`dJ/dt\omega ^3`$. Apparently, this relationship breaks down among very fast rotators, a phenomenon referred to as magnetic “saturation”.

But there are complications when one deals with young stars. It appears that in addition to contracting on the PMS, a star can also interact with its circumstellar disk, and this can have profound effects on the rotational evolution. Once on the main sequence, winds dominate the rotational evolution, except that internal transport of angular momentum can also play an important role.

Thus, the present framework for interpreting the observations involves: pre-main sequence evolution, magnetic saturation, disk interaction, winds, and internal transport. The complexity of these phenomena requires observational guidance to decode. On the PMS, a wealth of data has been provided by Herbst and collaborators, Stassun et al. (1999) and Rebull (2001). On the main sequence, two $`v\mathrm{sin}i`$ datasets have emerged recently, for the open clusters M 34 (Soderblom et al. 2001) and NGC 2516 (Terndrup et al. 2002). With ages of approximately 250 Myr and 150 Myr respectively, these observations are helping to fill in the gap in the observations between the young clusters ($``$ 100 Myr) and the Hyades (600 Myr).

Today, significant numbers of rotation periods among main sequence stars are available222We reference here only the largest datasets. in IC 2391 (Patten & Simon 1996), IC 2602 (Barnes et al. 1999), IC 4665 (Allain et al. 1996), Alpha Per (Prosser & Grankin 1997; Prosser et al. 1993), the Pleiades (van Leeuwen et al. 1987; Krishnamurthi et al. 1998), and Hyades/Coma (Radick et al. 1987; Radick et al. 1990), spanning an age range from 30 Myr to 600 Myr. In the immediate future, we can expect the release of datasets in NGC 1039 (M 34), NGC 2516, and NGC 3532 (Barnes and collaborators), two of which were displayed in preliminary form at the June 2002 AAS meeting. The author is aware of on-going studies in M 34 and M 35 (S. Meibom, personal communication). These and other observations will ensure that it is possible to build a time sequence for a range of stellar masses, thus completely removing the ambiguities inherent in $`v\mathrm{sin}i`$ data.

It is already possible to identify trends in the complete dataset. These were pointed out during the meeting, and will be published in the immediate future. The color-magnitude diagrams of the newly studied clusters, often generated in conjunction with membership studies of one sort or another, appear to be very promising, and enable the identification even of very low mass cluster members. Trends in the amplitudes of rotational variability are slowly emerging, and we are beginning to decode its dependence on stellar mass, age, and rotation rate.

Lest implications on the PMS and during the early main sequence evolution distract us, it is worth noting that the Skumanich (1972) relationship is still an excellent description of main sequence stellar spindown, once the stars have been binned by stellar mass. This is explained at length in Barnes (2001), which considered whether and how the rotation rates of the extremely well-studied sample of Mt. Wilson stars can be connected to those in open clusters, and whether the rotation rates of the host stars of extra-solar planets could be distinguished from the others. The result is that, with a few minor exceptions that can be explained, a Skumanich-type spindown is a perfectly adequate description of the behavior of both the Mt. Wilson stars and of the planet host stars. The planet host stars are not particularly slow rotators, in contrast to expectations from disk-locking scenarios of stellar rotation.

## 3 Rotational Evolution of Solar-Mass Stars: Theory

### 3.1 Internal Angular Momentum Transport Processes—Marc Pinsonneault (Ohio State University)

Hydrodynamic mechanisms (meridional circulation, shear, and GSF instabilities) can cause mixing and they are effective in rapidly rotating stars. However, these mechanisms occur over a long timescale in slowly rotating stars. Turbulence in convective regions can excite waves where the restoring force is variations in the potential ($`g`$-modes); such waves can effectively transport angular momentum throughout the radiative cores of stars. Significant mixing would only occur near the edges of convection zones. However, the impact of this mechanism upon stellar models has been controversial and the overall impact on the behavior of stars is not clear. Finally, if there exists a magnetic field in the radiative core of a star it will tend to enforce corotation along field lines without inducing significant mixing. The impact upon stellar properties, however, depends significantly upon the overall morphology of the field.

We have two important constraints on the timescale for angular momentum transport in stars. Helioseismic data provide a measurement of the internal solar rotation as a function of radius and latitude. There is a shear layer located beneath the convection zone with a characteristic width of less than 0.05 solar radii, and the radiative core below the shear layer appears to rotate as a solid body. This indicates that angular momentum transport in the solar interior was highly effective, but it does not directly provide an indication of the relevant timescale. The second constraint is the spindown of slowly rotating stars. Because only a fraction of the total moment of inertia is contained in the outer convection zone, there will be a difference in the time evolution of the surface rotation rates between models where the timescale for internal angular momentum transport is short (e.g. the whole star must be spun down by the wind) and models where the timescale for internal angular momentum transport is long (e.g. initially only the convection zone is spun down, and the core is coupled at a later epoch.)

Models which include only angular momentum transport from hydrodynamic mechanisms predict that as stars lose angular momentum from a magnetized wind there will be an initial phase where the envelope spins down with respect to the core. At later ages the models will reach a rough equilibrium where the core is rotating significantly faster than the surface. The observed spindown of the slow rotator population in young open clusters is consistent with the early core-envelope decoupling predicted by such models, but the flat solar rotation profile is inconsistent with the strong internal angular velocity gradients predicted by such models. This combination of evidence suggests that there is an additional angular momentum transport mechanism that operates over timescales longer than 100 Myr but shorter than the lifetime of the Sun. Possible methods for distinguishing between waves and magnetic fields were also briefly discussed.

### 3.2 Magnetic Star-Disk Coupling Mech-anisms—Frank Shu (Tsing Hua Univ.)

333Written by K. Stassun from notes provided by F. Shu

The idea of magnetic “disk-locking” in some sense originates with the theory developed by Ghosh & Lamb (1979) in the context of neutron stars. In their model, a stellar dipole field, $`B`$, everywhere threads a circumstellar disk, which is truncated at an inner radius, $`R_t`$, where magnetic and ram pressures balance. In this picture, if $`R_t`$ is within the co-rotation radius, $`R_c`$ (the radius in the disk that rotates with Keplerian angular velocity equal to the stellar angular velocity), then gas interior to $`R_c`$ spins up the star while gas exterior to $`R_c`$ spins it down.

But the Ghosh & Lamb picture presents some problems in the context of young, low-mass stars. By construction, the magnetic field is strong enough marginally to compete with centrifugal forces in the funnel flow. Thus, since the disk must be much more dense than the funnel flow, the field threading the disk cannot resist wrapping except at the single radius $`R_c`$. Continuous wrapping of the field would lead to magnetic reconnection, but the observed X-ray activity of classical T Tauri stars (those with disks) is similar to that of weak-lined T Tauri stars.

Our own X-wind model avoids these difficulties by bringing the dipole field in the disk together at the X-point, $`R_xR_c`$. All of the unperturbed dipole exterior to $`R_x`$ is brought in by disk accretion to $`R_x`$, and half of the interior field is brought to $`R_x`$. This pinch towards $`R_x`$ is caused by back-reaction to the funnel flow and X-wind. This pinch distorts the dipole configuration near $`R_xR_c`$ so that differential rotation of the disk no longer wraps the field, even though the disk is much denser than the gas in the flow/X-wind, making a steady state configuration possible.

The “dipole” assumption has problems when confronted by observation, however. For example, the results of polarization studies are inconsistent with well-organized fields at the stellar surface, and the hot spots that have been observed on classical T Tauri stars have much smaller areas than implied by the dipole assumption. Furthermore, completely convective stars are not expected to have dipole field geometries near the surface.

Fortunately, X-wind theory is not dependent on the (unperturbed) magnetic field being dipolar. The critical invariant in the theory is the assumption of “trapped magnetic flux” at the X-point, $`R_x`$. Relaxing the dipole assumption, we find that “disk-locking” in the X-wind model is only very slightly modified in the “multipole” case, but the model now allows for the smaller hot spots and polarization levels that are observed.

We make two critical assumptions:

(1) Disk locking: $`\mathrm{\Omega }_{}=\mathrm{\Omega }_x`$, and

(2) Flux trapping: $`\mathrm{\Phi }_t=\alpha (GM_{}\dot{M}_D/\mathrm{\Omega }_{})^{1/2}`$,

where $`\mathrm{\Phi }_t`$ is the magnetic flux at the truncation radius, and $`\alpha `$ is an order unity coefficient that depends on details of X-wind theory.

In the X-wind model, the trapped flux is 1/3 in the outflowing X-wind, 1/3 in the “dead zone” (static coronal gas), and 1/3 in the funnel flow. Since the 1/3 of the trapped flux in the funnel flow is anchored in the stellar surface, we have the prediction that:

$`B(f4\pi R_{}^2)=1/3\mathrm{\Phi }_t`$,

where $`B`$ is the field strength in the footpoint on the stellar surface, and $`f`$ is the surface covering fraction of that footpoint. Assuming that $`B`$ does not vary much from star to star, this leads to the testable prediction that:

$`fR_{}^2(M_{}\dot{M}_DP_{\mathrm{rot}})^{1/2}`$,

where $`P_{\mathrm{rot}}=2\pi /\mathrm{\Omega }_{}`$.

Recently, Johns-Krull & Gafford (2002) compared observations of T Tauri stars in Taurus-Auriga to the “disk-locking” prediction of this generalized, “multipole” X-wind model. Their results, shown in Fig. 4, seem to verify this prediction.

We close by commenting on the eventual breakdown of disk locking. At some point, as $`\dot{M}_D`$ decreases in time, disk locking must break down, because the star has too much inertia for the then-meager funnel flow to force the star to rotate at $`\mathrm{\Omega }_{}=\mathrm{\Omega }_x`$. However, unlike the situation where the stellar field is a fixed dipole, disk locking may have some flexibility built into it from the ability of corotating multipole fields to accommodate even fairly low disk accretion rates. Because of this complication, and because we do not know how stellar dynamos work in the presence of interactions with a surrounding accretion disk, we have no first-principles, quantitative theory of disk locking.

## 4 Panel Discussion: Observations and Theory vis-a-vis Modeling

The session included a panel discussion which focused on a set of guiding questions posed by Keivan Stassun (University of Wisconsin), who moderated the discussion. The aim of the panel discussion was for the observers, theorists, and modelers to begin synthesizing work from one another’s domains, to confront areas of uncertainty, and to identify future research directions. The specific questions discussed by the panel were:

Disks:

1. How can we understand the observed relationship between rotation and infrared excesses \[among PMS stars in Orion\] in the context of star-disk coupling?

Observations in Orion show a positive correlation between rotation period and magnitude of near-IR excess. This would seem to suggest that the inner edges of the circumstellar disks in these stars know something about the rotation rate of the star. This is puzzling because near-IR emission probes the dust in the disk, not the gas that is thought to magnetically couple to the star. And, if any correlation might be expected, it is in the opposite sense of what is seen. The panel was not able to put closure on this issue, and it remains unclear what this observational result might imply for the physics of magnetic star-disk interaction.
2. Does star-disk coupling always act to remove angular momentum from a star?

Models of PMS angular momentum evolution typically invoke “disk-locking” as a way of extracting angular momentum from the star. However, strictly speaking, accretion torques can act positively on the star as well as negatively, depending on the balance between stellar magnetic field strength and the accretion rate in the disk. Panelist Frank Shu confirmed that in some cases stellar spin-up can result from star-disk interaction.

Winds:

1. What do we know (or wish we knew) about solar-analog winds from young stars and their role in early angular momentum evolution?

Panelist Marc Pinsonneault suggested that more quantitative observations of stellar magnetic field strengths and of mass outflow rates would be extremely helpful in more rigorously constraining models. Indeed, there is some uncertainty about the basic relationship between stellar rotation, dynamo-generated field strength, and mass outflow rates. In addition, observations that probe the internal rotation profiles of young stars are needed, specifically to test the prediction that young stars are strong differential rotators internally.

With respect to future research directions, the following were presented as key areas of observational focus:

$``$ rotation measurements prior to 1 Myr

$``$ rotation measurements in the “age gaps” (primarily between 1 Myr and 100 Myr); connect rotation period distributions at various ages

$``$ better assessment of biases in measured rotation rates at various ages

$``$ ascertain the role of stellar multiplicity on the angular momentum evolution of single PMS stars

## 5 New Frontiers at Other Masses: Inter-mediate-Mass Stars and Brown Dwarfs

### 5.1 Intermediate Mass Stars: From the Birthline to the Main Sequence—Sidney C. Wolff (NOAO), S. Strom (NOAO), L. Hillenbrand (Caltech)

We have measured projected rotational velocities ($`v\mathrm{sin}i`$) for a sample of 145 stars located in the Orion star-forming complex with masses ($`M`$) between 0.41 and 15 M. These measurements establish, as nearly as direct observations can, the initial values of specific angular momentum ($`jJ/M`$) for stars in this mass range. The goal of this work is to try to account for the initial values of $`j`$ in terms of theories of star formation and angular momentum regulation. We then map the initial values of $`j`$ onto the ZAMS and discuss the processes that affect the evolution of angular momentum from the birthline to the ZAMS.

Our data (see Fig. 5) show that there is a continuous power law relationship between $`j`$ and $`M`$ for stars on convective tracks with masses in the range $`0.5`$ (and probably less) to 2–3 M and that this power law merges smoothly with the relationship for more massive stars, which are already on the ZAMS. This power law can be compared with the predictions of models that posit that stars are “locked” to circumstellar accretion disks until they are released at the birthline. If we assume that the accretion rate is $`10^5`$ M yr-1 for 1 M stars, and that mass accretion rates are proportional to stellar mass, we can account for the observed slope and zero point of the upper envelope of the observed power law relationship between $`j`$ and $`M`$.

If we compare the initial values of $`j`$ with those of stars on PMS radiative tracks or on the main sequence, we find the primary difference to be a break in the power law for stars with $`M<2`$ M that have completed the convective phase of evolution. This sharp decrease in $`j`$ has been known for a long time for field stars, but the decrease is seen already in the Orion stars, which are less than 1 Myr old and are still on PMS radiative tracks. We show that the break in the power law is a consequence of core-envelope decoupling at the time of the transition from the convective to radiative tracks. The difference in rotation rates seen for stars on either side of this transition are consistent with conservation of angular momentum in shells and inconsistent with solid body rotation.

We also argue that the systematic variations seen in $`v\mathrm{sin}i`$ along the ZAMS for stars with masses greater than 2–3 M, i.e. a maximum in $`v\mathrm{sin}i`$ in the late B-type stars, lower $`v\mathrm{sin}i`$ in the early B-type stars, and a large fraction of B2–B4 stars rotating at less than 100 km/s, can be accounted for by variations in the predicted mass-radius relationship along the birthline, primarily associated with deuterium shell burning, with the amount of subsequent spin-up being determined by the difference between the radius on the birthline and the radius on the ZAMS.

### 5.2 Rotation and Activity of Ultra-Cool Dwarfs—Adam Burgasser (UCLA)

The study of angular momentum evolution in ultra-cool dwarf stars and brown dwarfs is just beginning, as these faint objects have until quite recently eluded detection (for a recent review on brown dwarf searches, see Basri 2000). However, early indications show a connection between rotation and stellar activity quite opposite that seen in hotter stars (Noyes et al. 1984).

High-resolution $`v\mathrm{sin}i`$ measurements by Basri et al. (2000) and Reid et al. (2002) have shown that in the field, cool L dwarfs (Kirkpatrick et al. 1999; Martín et al. 1999) are systematically faster rotators than M dwarfs, with no L dwarfs observed to date rotating slower than $`v\mathrm{sin}i=10`$ km/s. At the same time, magnetic H$`\alpha `$ emission shows a steep and steady decline in both frequency and strength beyond spectral type M7 (Gizis et al. 2000; Burgasser et al. 2002c), possibly due to the low ionization fractions in the upper atmospheres of the cooler dwarfs (Mohanty et al. 2002). This implies that rapidly rotating L dwarfs, which can approach $`v\mathrm{sin}i`$ $``$ 60 km/s (Basri et al. 2000), do not have strong emission, in contrast to hotter rapidly-rotating stars which are typically the most active.

When they are young, however, brown dwarfs are apparently not as extreme in their rotation velocities. Joergens & Guenther (2001) have shown that substellar objects in Cha I rotate as rapidly as their T-Tauri counterparts. Basri (2002) has proposed that the absence of magnetic activity, and hence magnetic braking, allows brown dwarfs to retain their angular momentum as they contract. However, this appealing hypothesis is complicated by the detection of disks around young brown dwarf candidates (Muench et al. 2001; Testi et al. 2002) and the presence of magnetic flares even in the L dwarfs (Hall 2001). Investigations of rotation velocities in young cluster brown dwarfs are currently underway to clarify this picture (G. Basri, priv. comm.). Theoretical studies of disk braking in these cool, low-mass systems are desperately needed.

Photometric variability studies have also begun to measure rotation periods in ultracool dwarfs, but results have been somewhat confusing. Work by Bailer-Jones & Mundt (1999, 2001), Martín et al. (2001), Clarke et al. (2002), and Gelino et al. (2002) have found photometric variations in late-M and L dwarfs on order 0.01–0.1 mag with aperiodic or periodic (timescales of a few hours) modulations. However, periods have been difficult to characterize in general, as some objects show multiple and equally significant periods (Martín et al. 2001), while others exhibit evolution in the photometric periods (Bailer-Jones & Mundt 2001; Gelino et al. 2002).

The difficulty in obtaining robust rotational periods for these objects is likely related to the presence of condensate clouds in their atmospheres (Ackerman & Marley 2001), which may produce rotationally modulated features that evolve over a similar timescale (Gelino et al. 2002). Indeed, the transition from L dwarfs to T dwarfs (Burgasser et al. 2002a; Geballe et al. 2002) may be driven by the sudden disruption of cloud layers by convective cells (Burgasser et al. 2002b), obscuring any rotational modulation. Various studies of photometric variability are currently underway, however, because of the feasibility of such projects on small (1-4m) telescopes, in contrast to the high sensitivity required to measure $`v\mathrm{sin}i`$ in these faint dwarfs.

Keivan G. Stassun
University of Wisconsin—Madison

Donald Terndrup
Ohio State University


text:
# SCUBA Observations of High Redshift Radio Galaxies

## 1 Introduction

The discovery of a tight relation between the masses of the stellar bulges of galaxies and the masses of their central black holes (Magorrian98aj, ) has reemphasized the importance of high redshift active galactic nuclei (AGNs) as probes of galaxy formation. They are therefore the target of many studies of starformation in their host galaxies. Since a large proportion of the optical and UV emission from a starburst is absorbed and reprocessed by dust to restframe far-IR wavelengths, a good way to investigate the roles of dust and star formation in these distant galaxies is to measure their (sub-)millimeter continuum emission. Another advantage of observing at those wavelengths is that dimming due to cosmological distance is compensated for as the peak in the dust emission is redshifted into the observed passband (e.g., BlainLongair93mnras, )

The formation of massive elliptical galaxies may be best studied by tracing the evolution of powerful radio galaxies. They require the presence of an actively-accreting super massive black hole ($`10^9\mathrm{M}_{}`$; Lacy01apj, ), are very extended and luminous at many wavelenghts allowing detailed studies, and their selection at radio wavelengths is not affected by the presence of dust.

Archibald et al. (2001) have conducted a submillimeter survey of 47 radio galaxies over an interval of $`0.7<z<4.4`$. They found a substantial increase in 850 $`\mu \mathrm{m}`$ detection rate with redshift and that the average 850 $`\mu \mathrm{m}`$ luminosity $`L_{850}`$ rises at a rate $`(1+z)^{34}`$ out to $`z4`$. While it is true that a buried AGN could be responsible for heating the dust, there are strong indications that starbursts (SFR $`1000\mathrm{M}_{}\mathrm{yr}^1`$) are the dominating mechanism: (i) there is no strong correlation between radio power and submillimeter luminosity (Archibald et al., 2001; Reuland03prepb, ), (ii) for all hyperluminous infrared galaxies the starburst outshines the AGN at rest-frame wavelengths longwards of 50 $`\mu \mathrm{m}`$ (Rowan-Robinson, 2000; Farrah et al., 2002), and (iii) the extent (up to 30 kpc) of the CO and dust-emission found in 3 HzRGs (Papadopoulos00apj, ; DeBreuck02prep, ) is not readily explained by a central process.

Until recently only few HzRGs at $`z>3`$ were known, hampering statistical studies. However, using radio and near-IR selection we have significantly increased the sample of $`z>3`$ HzRGs at Keck (e.g., DeBreuck01aj, ). Here we present submillimeter observations of these. Additionally, we have observed an intriguing class of HzRG-candidates for which no spectroscopic redshift could be determined. One of these candidates had been studied before in detail (WN J0305+3525; Reuland03apj, ) and was found to be very submillimeter luminous.

## 2 Sample Selection and Observations

We have selected all known HzRGs with redshifts $`z>3`$ without prior SCUBA observations. We aimed to observe a significant sample of HzRGs to complement the observations of Archibald et al. (2001) and obtain better statistical info at the highest redshifts. The sample is the result of an ongoing effort by our group and others (e.g., DeBreuck01aj, ) to find distant radio galaxies based on Ultra Steep Radio Spectrum (USS) and near-IR identification selection criteria. Also 10 HzRG-candidates have been observed that did not show optical emission lines in deep spectroscopic observations at Keck.

The observations have been carried out over 4 semesters in 4 years (August 1999 $``$ January 2002) with SCUBA (Holland98spie, ) at the James Clerk Maxwell Telescope. We observed at 450 $`\mu \mathrm{m}`$ and 850 $`\mu \mathrm{m}`$ wavelengths in photometry mode. The optical depth $`\tau _{850}`$ varied between 0.14 and 0.38 with an average value of 0.26. The data were reduced following standard procedures outlined in the SCUBA Photometry Cookbook.

## 3 Results and Discussion

### 3.1 HzRGs with $`z>3`$

In the course of our program we observed 23 HzRGs, of which 12 were detected with a signal-to-noise ratio larger than 3. Adding these data to the sample of Archibald et al. (2001) increases the number of 3$`\sigma `$-detections at $`z>3`$ from 5 to 16. There was no attempt to correct for a contribution to the submillimeter continuum from the non-thermal radio emission, since most sources were selected using an USS criterion. Moreover, extrapolation from the radio regime to submillimeter wavelengths is very uncertain at best, and is likely to result in an overestimate of the non-thermal flux density due to steepening of the radio spectrum. For this reason (and following Willott02mnras, ) we have chosen to use the “un-corrected” values from Archibald et al. (2001). Figure 1 shows two histograms representing the number of HzRGs that were observed at each redshift and the number of HzRGs that were actually detected given a 2$`\sigma `$ or 3$`\sigma `$ detection criterion. Regardless of the actual detection limit used, a marked increase in the fractional detection rate from $`15\%`$ at $`z<2.5`$ to $`50\%`$ at $`z>2.5`$ is apparent, confirming the result found by Archibald et al. (2001).

We have correlated the inferred submillimeter luminosities with radio power to search for signs of significant heating by the AGN. No strong correlation is apparent (Fig. 2). Recently, it was found that the size of the radio source and $`L_{850}`$ are correlated for a sample of radio-loud QSOs at $`z=1.5`$, with smaller (i.e. younger; Blundell99aj, ) sources being more luminous (Willott02mnras, ). This effect was attributed to a possible relation between the jet-triggering event and a short-lived starburst or quasar-heated dust in QSOs. The effect was not found for a matched sample of radio galaxies at $`z1.5`$. We have checked whether there is a correlation between $`L_{850}`$ and radio source size (age) in our larger sample. Figure 2 shows that they are not strongly correlated. However, size of the radio source depends not only on age but also on jet-power and enviroment which is hard to correct for in our heterogeneous sample. Therefore we cannot exclude that $`L_{850}`$ and age of the radio source are correlated.

The (sub-)millimeter emission of starbursts should be largely indepent of viewing angle. Unification schemes (e.g., Barthel89apj, ) therefore suggest that we can combine observations of HzRGs with those of QSOs, allowing us to study the evolution of starformation in their host galaxies over a larger redshift range. Figure 3 shows the average observed (sub-)millimeter fluxdensity as a function of redshift for mostly radio-quiet QSO samples (Carilli01apj, ; Omont02aa, ; Isaak02mnras, ; Priddey02mnras, ) and HzRGs. Reassuringly, the trends mimic each other and they indicate that the average observed flux density rises to redshifts of $`z4`$ after which it declines.

### 3.2 Obscured “no-$`z`$” HzRGs

During our search for HzRGs among USS sources, we found that 33% of the selected canditates were at redshifts $`z>3`$ (DeBreuck01aj, ). However we also found that $`24\%`$, mostly with a compact radio morphology, do not show clear signs of emission lines and that a significant fraction ($`10\%`$ of the entire sample) does not even show continuum emission. Possible explanations for this may be that (i) these are galaxies in the “redshift desert” ($`1.5<z<2.3`$) where only weak emission lines fall in the optical window, (ii) they are pulsars, (iii) they are extremely obscured AGN (iv) or at such high redshift ($`z>7`$) that Ly$`\alpha `$ is shifted out of the optical passband.

Figure 4 shows one of these objects, WN J0305+3525, which was studied in detail (Reuland03apj, ). It is a strong (sub-)millimeter detection with $`S_{850}`$ = $`12.5\pm 1.5`$ mJy and $`S_{1250}`$ = $`4.2\pm 0.6`$ mJy, and appears associated with a number of faint $`K21`$ mag and $`J23`$ mag objects. The radio morphology is compact ($`\theta <1.9`$$`^{\prime \prime }`$), as is typical of the other candidates. Based on circumstantial evidence like the near-IR magnitudes and morphology, the (sub-)millimeter ratios and compact radio source morphology it was argued that WN J0305+3525 is a highly obscured HzRG possibly at a redshift $`z=3\pm 1`$ in a young stage of its evolution (Reuland03apj, ). We have detected 3 more of these HzRG-candidates with flux densities of $`S_{850}5\pm 1`$ mJy resulting in a detection rate of 4 out of 10, very similar to what is found for the confirmed HzRGs. This further supports the idea that they may be young, extremely obscured HzRGs.

## 4 Summary

* Submillimeter detection rate appears to be primarily a function of redshift. If this is interpreted as being due to a change in the intrinsic far-IR luminosity, it would be consistent with a scenario in which the bulk of the stellar population of radio galaxies forms rapidly around redshifts of $`z=35`$ after which they are more passively evolving.
* A significant fraction (10%) of USS sources may be very obscured radio galaxies. The obscuration may be related to their young evolutionary stage.
* In contrast to the submillimeter sources, HzRGs have accurately determined redshifts and host identifications. Even with the advent of SIRTF this is not likely to change. Therefore, HzRGs are a key population for studies of galaxy formation in the early universe, allowing detailed follow-up mm-interferometry observations to study their dust and gas content.

The work of M.R. and W.v.B. was performed under the auspices of the U.S. Department of Energy, National Nuclear Security Administration by the University of California, Lawrence Livermore National Laboratory under contract No. W-7405-Eng-48.


text:
# Dark energy, dissipation and the coincidence problem

## Acknowledgments

This work has been partially supported by the Spanish Ministry of Science and Technology under grant BFM 2000-C-03-01 and 2000-1322, and the University of Buenos Aires under Project X223.


text:
# Dissipative flow and vortex shedding in the Painlevé boundary layer of a Bose Einstein condensate

## Abstract

Raman et al.R  have found experimental evidence for a critical velocity under which there is no dissipation when a detuned laser beam is moved in a Bose-Einstein condensate. We analyze the origin of this critical velocity in the low density region close to the boundary layer of the cloud. In the frame of the laser beam, we do a blow up on this low density region which can be described by a Painlevé equation and write the approximate equation satisfied by the wave function in this region. We find that there is always a drag around the laser beam. Though the beam passes through the surface of the cloud and the sound velocity is small in the Painlevé boundary layer, the shedding of vortices starts only when a threshold velocity is reached. This critical velocity is lower than the critical velocity computed for the corresponding 2D problem at the center of the cloud. At low velocity, there is a stationary solution without vortex and the drag is small. At the onset of vortex shedding, that is above the critical velocity, there is a drastic increase in drag.

Dilute Bose-Einstein condensates have recently been achieved in confined alkali-metal gases and the study of vortices therein is one of the key issues. Raman et al. R ; RO , Onofrio et al. O  have studied dissipation in a Bose Einstein condensate by moving a blue detuned laser beam through the condensate at different velocities. They found experimentally a critical velocity for the onset of dissipation. This critical velocity has been related to the one found by Frisch et al. FPR  for the problem of a 2D superfluid flow around an obstacle in the framework of Nonlinear Schrodinger Equation (NLS): below a critical velocity, the flow is stationary and dissipationless, while beyond this critical velocity, the flow around the disc becomes time dependent and vortices are emitted. Numerical simulations have been done for this type of problem in 2D HB  and 3D JMA ; W . In particular, the direct 3D simulation of W  shows the plot of the drag against the velocity. A critical velocity can be numerically computed when the drag becomes nonzero, but no precise mechanism of vortex nucleation is described by the authors. This critical velocity has been analyzed theoretically for a homogeneous 2D system SZ  and an inhomogeneous 2D system C ; FS .

In this paper, we want to take into account the 3D geometry of the experiment of R ; O ; RO . Our aim is to understand the mechanism of vortex nucleation in the boundary region. Indeed the analysis of FPR  allows to understand what is happening in the interior of the cloud, where the kinetic energy is negligible in front of the interaction energy. In the region where the laser beam crosses the boundary of the cloud, the sound velocity gets small, since the amplitude of the wave function becomes small. There, the kinetic energy term can no longer be neglected in front of the trapping and interaction terms. We blow up this region in such a way that the trapping potential varies linearly with the distance to the boundary and far away from the laser beam, the wave function is then given by a Painlevé equation. We analyze the behavior of the wave function in the frame of the laser beam. The real experiments are quite complex, and in particular here we do not take into account the oscillations and acceleration of the beam but we believe that our analysis allows to understand the mechanism of increase of drag. One of our main results is that there is always a drag around the laser beam and this drag grows continuously. At low velocity, the drag is not a consequence of the shedding of vortices, and finally of a time dependent density and velocity field. The origin of this drag is in the radiation condition for the wavefield: the motion changes continuously the structure of the solution seen in the frame of reference of the ”fluid” at infinity. We study the transition toward a time dependent regime of vortex shedding, which happens at a critical velocity. The critical velocity that we find is lower than the 2D critical velocity at the center of the cloud coming from the computation of FPR . Vortices are nucleated close to the boundary of the cloud and the tubes grow and detach to form rings that move downstream. When tubes are emitted, significantly large drag values are observed. The drag increases smoothly as the velocity increases.

The dynamics can be modeled using the Gross Pitaevskii equation at zero temperature with an external trapping potential $`V_{tr}=m/2(\omega _x^2x^2+\omega _y^2y^2+\omega _z^2z^2)`$.

$$i\mathrm{}_t\mathrm{\Psi }=\frac{\mathrm{}^2}{2m}\mathrm{\Delta }\mathrm{\Psi }+(V_{tr}+Ng|\mathrm{\Psi }|^2)\mathrm{\Psi }.$$

If an object is moved inside the condensate, $`V_{tr}`$ has to be replaced by $`V_{tr}+V_{ob}`$, where $`V_{ob}`$ depends on $`xvt`$. Based on the experimental data of R ; O , we take $`a=mg/4\pi \mathrm{}^2=2.94nm`$, $`N=\mathrm{1.2.10}^7`$, $`\omega _y=\omega _z=377s^1`$, and $`\omega _x=\lambda \omega _z`$, with $`\lambda =0.3`$. We also define the characteristic length $`d=(\mathrm{}/m\omega _z)^{1/2}=2.71\mu m`$ and a small nondimensionnalized parameter $`\epsilon `$ given by

$$\epsilon =(\frac{d}{8\pi Na})^{2/5}.$$

We find that $`\epsilon =\mathrm{6.21\; 10}^3`$ which may be viewed as small parameter and allow rescaling the equation near the edge of the condensate. Re-scaling the distances by $`R=d/\sqrt{\epsilon }=34.4\mu m`$, the time by $`1/(\epsilon \omega _z)`$, we have $`\psi (𝐫,t)=R^{3/2}\varphi (\stackrel{~}{𝐫},\stackrel{~}{t})`$ where $`\stackrel{~}{𝐫}=R𝐫`$. In these new units, the radii of the condensate are $`R_y=R_z=0.65`$ and $`R_x=2.18`$. The laser beam is modeled by an obstacle which is a cylinder $`𝒞`$ of axis $`z`$ and radius $`l=0.19`$ on which $`\psi =0`$. We will work in the frame where the obstacle is stationary. Outside the obstacle, the equation can be rewritten as

$$2i_t\psi =\mathrm{\Delta }\psi +\frac{1}{\epsilon ^2}(\rho _{TF}|\psi |^2)\psi ,$$

where $`\rho _{TF}=\rho _0(\lambda ^2x^2+y^2+z^2)`$ is the Thomas Fermi limit density and $`\rho _0=0.42`$ is the rescaled chemical potential. Note that $`|\psi |^2`$ is close to its Thomas Fermi value $`\rho _{TF}`$ except near the obstacle and near the boundary of the cloud. This boundary layer has a thickness of order $`\epsilon ^{2/3}`$ so that we rescale the domain with $`\psi (\stackrel{~}{x},\stackrel{~}{y},\stackrel{~}{z})=\epsilon ^{1/3}u(x,y,z)`$, where $`x=\stackrel{~}{x}/\epsilon ^{2/3}`$, $`y=\stackrel{~}{y}/\epsilon ^{2/3}`$ and $`z=(\sqrt{\rho _0}\stackrel{~}{z})/\epsilon ^{2/3}`$, $`v=\stackrel{~}{v}\epsilon ^{2/3}`$. By blowing up the boundary of the cloud near $`z=0`$, and truncating at $`z=L`$, the rescaled layer thickness, we see that the modulus of the stationary solution in the boundary layer for $`|x|`$ and $`|y|`$ large, that is far away from the obstacle, is given by the solution of the first Painlevé equation FF ; DPS .

$$p^{\prime \prime }+(2z\sqrt{\rho _0}p^2)p=0,p(0)=0,p(L)=\sqrt{2\sqrt{\rho _0}L}.$$
(1)
We choose the size of the boundary layer $`L`$ so that $`\epsilon ^{2/3}L=3\sqrt{\rho _0}/10`$. This is based on the consideration that, on the one hand, $`\epsilon ^{2/3}L`$ should be suitably small so that $`2z\sqrt{\rho _0}`$ is a good approximation for $`\rho _{TF}=\rho _0\stackrel{~}{z}^2`$ in the boundary layer and on the other hand the critical velocity at $`z=L`$ is not too different from the critical velocity at the center of the cloud. The obstacle is now a cylinder of radius $`a=l/\epsilon ^{2/3}=5.6`$.

The obstacle moves at the rescaled velocity $`v=v_{exp}/(\epsilon ^{1/3}\omega _zR)`$, and in the frame of the obstacle, the equation becomes

$$2i_tu=\mathrm{\Delta }u2iv_xu+(2z\sqrt{\rho _0}|u|^2)u.$$
(2)
We want to understand the behaviour of solutions depending on $`v`$. If we restrict (2) to $`z=L`$, we can perform a similar analysis to FPR  and get the value of the critical velocity for the onset of vortex shedding and find $`v_c^2=2\sqrt{\rho _0}L/11=2c_s^2/11`$, where $`c_s`$ is the sound velocity. Of course, we cannot apply this analysis in the low density region, since there the sound velocity gets close to 0. Another mechanism has to be understood. The rescaled drag around the obstacle is

$$drag=\frac{1}{2}_𝒞(u_x\overline{u}_n\overline{u}_xu_n)𝑑l𝑑z.$$
(3)
We first analyze the stationary solution of (2) in the very low density region, where the system is very dilute and one can neglect the nonlinear term. In fact, a precise condition is that $`p^2`$ is less than $`v^2`$, which gives a truncation point $`z_c`$ at which $`p^2(z_c)=v^2`$. It is rather straightforward in classical scattering theory to compute the perturbed wavefield and finally the drag on the obstacle (a related problem, the scattering of sound by a cylinder, is treated in M ). In the low density region, it is reasonable to look for $`u`$ with the following ansatz

$$u(x,y,z)=p(z)\psi (x,y)e^{ivx}.$$
(4)
We can first approximate $`p(z)`$ in this region by an Airy function given by the solution of $`p^{\prime \prime }+2zp\sqrt{\rho _0}=0`$, that is, by defining $`\overline{z}^3=1/(2\sqrt{\rho _0})`$, we have

$$p(z)\sqrt{2}Ai(\frac{z}{\overline{z}})\frac{1}{\sqrt{2\pi }}(\frac{\overline{z}}{z})^{\frac{1}{4}}exp(\frac{2}{3}(\frac{z}{\overline{z}})^{\frac{3}{2}}).$$
(5)
Then, outside the obstacle, $`\psi `$ is a solution of the 2D Helmholtz equation

$$\mathrm{\Delta }\psi +v^2\psi =0$$
(6)
with $`\psi =0`$ for $`r=a`$, the obstacle boundary, and $`\psi e^{ivx}`$ at infinity. This solution can be computed M  in terms of Bessel functions $`J_k`$ and $`N_k`$. One finds that, at leading order for $`v`$ small, the 2D drag of $`\psi `$ is proportional to

$$v^2J_0^2(v)J_1(v)N_2(v)/N_0^2(v)v/\mathrm{log}^2v.$$
(7)
The total drag has to be multiplied by the integral of $`p^2`$ along the $`z`$ axis to the truncation point $`z_c`$ defined by $`p^2(z_c)=v^2`$. Direct calculation gives

$$\frac{v}{\mathrm{log}^2v}_{\mathrm{}}^{z_c}p^2𝑑zC\frac{v^3}{\mathrm{log}^{8/3}v}.$$
(8)
In conclusion, the total scattering drag tends to zero at low speed. It is plotted in Figure 1 (solid line).

To understand how the solutions of (2) behave, we numerically integrate the equation in a computational domain of dimension $`60\times 60\times L`$ with periodic boundary conditions in $`x`$ and $`y`$ and taking $`u=0`$ on the boundary of the obstacle and away from the condensate ($`z=0`$). At the truncated surface $`z=L`$ inside the condensate, we use the condition $`(/z)(u/p)=0`$, where $`p`$ is the solution to the Painlevé equation (1) described before. The numerical solution is computed based on a continuous piecewise quadratic finite element approximation in space and the Runge-Kutta fourth-order in time integration scheme. Using $`p`$ as initial condition, we first compute the solution of (2) for some time by adding a damping coefficient, that is, we replace $`iu_t`$ in (2) by $`iu_t(1+i\gamma )`$. For small velocity, this effectively drives the numerical solution of (2) close to a stationary solution. Then, we continue the integration with a much reduced damping coefficient $`\gamma =0.02`$ or with no damping $`\gamma =0`$.

In what follows, we will divide the velocity by the sound velocity at center $`c_s=\sqrt{2\rho _0}/\epsilon ^{1/3}`$. In Figure 1, we plot the drag vs. the velocity divided by $`c_s`$. For a given velocity value, the drag is the value obtained through time averaging of (3). We have verified that with different small values of $`\gamma `$, the drag calculation remains essentially the same.

For $`v`$ small, we find that the solution is almost stationary. Surface oscillations are present near $`z=0`$, and the drag is small, but not zero. See Figure 2 for plots of the solution. The drag computed in this regime fits very well to the cubic growth given by (8).

When $`v`$ is increased, at a critical velocity $`v_c/c_s0.2`$, the surface oscillations develop into small handles that move up and down the obstacle without detaching; see Figure 3.

There is no stationary solution, but no vortex shedding either: the small handles move up the obstacle to a critical $`z`$ value and down. This instability may be related to the one discussed by Anglin A1 : in our scaling, the critical velocity found in A1  is 0.2. At this stage, the solutions do not produce large drag nor vortex shedding.

It is only for larger velocities ($`v/c_s>0.25`$) that the handles move up to the top, detach from the obstacle and produce significant drag. This is a wholly nonlinear phenomenon and most likely cannot be described by a linear analysis.

Let us describe the solutions for $`v/c_s>0.25`$ illustrated in Figure 4.

The vortex handles seem to first nucleate near $`z=0`$ and are top connected to the obstacle. As time increases, the bottom ends move away from the obstacle in a slightly down stream direction while the top end moves up along the obstacle (Figure 4a). When the top ends of the vortices become close to $`z=L`$, the bottom ends reverse their trend of moving away from obstacle. Instead, they move back to the bottom of the obstacle, as if the handles prefer certain curvature (Figure 4b). Eventually, the top ends of the handle move away from the obstacle and produce a pair of vortex tubes with their bottom ends at the bottom of the obstacle (Figure 4c). The handles merge into a half vortex ring, this half ring moves both upward and downstream (Figure 4d). Near $`z=0`$, the solution can be approximated by the solution (4) and this solution does not have vortices, so the instability creates the vortex but the vortex moves away. Vortex detachment happens only at sufficiently high density, in the region where the nonlinear term in the equation dominates. The direction of the vortex displacement is due to the velocity of the flow and the self interaction of the vortex on itself, which gives a movement along its normal vector. Meanwhile, while the vortex ring starts to detach from the obstacle, another pair of vortex handles is forming near the obstacle. The above process repeats itself. Note that we have truncated the domain close to the boundary of the cloud, so that the half ring we compute would correspond to a closed ring in the experiments.

We have to point out that the critical velocity we have found for the onset of vortex shedding is lower than the critical velocity for the 2D problem at $`z=L`$. In this case $`v_{2D}/c_s=0.35`$. So the inhomogeneity in the condensate lowers the critical velocity from the 2D value. One can check that for different $`L`$, the critical velocity does not change. This is verified by our numerical computation where we have used two boxes with one about 50% higher in $`z`$ than the other, and there is little change in the drag plots, nor there is any significant difference in the dynamic behavior of the solutions.

In the experiments R ; RO ; O , the drag is plotted vs velocity and a critical velocity can be defined when a sharp slope is observed in the drag plot. The critical velocity in R  is very similar to ours, though slightly smaller. This is certainly due to the finite extent of the condensate in the $`x`$, $`y`$ direction. Indeed, our simulations have not taken into account that the cloud is narrower in the $`y`$ direction than along the $`x`$. We can check that for the 2D problem, this geometry lowers the velocity. On the other hand, our computations indicate that the inhomogeneity in the $`z`$ direction and the soft boundary of the laser beam are well accounted for by our problem.

Summary. We have studied the onset of dissipation in the Painlevé boundary layer of a BEC when a detuned laser beam is moved in the condensate. We do a change of frame and blow up the low density region near the boundary of the cloud to write the equation for the wave function in this region: $`z=0`$ is now the boundary of the cloud and $`z`$ large is the center. For small velocity, there is a drag around the obstacle due to radiation, but no vortex is generated: it is a stationary flow, which is supersonic near $`z=0`$, but subsonic for $`z`$ larger. On the other hand, when the critical velocity is reached, the instability propagates towards the top, a vortex handle is nucleated and detaches from the obstacle to form vortex rings that move away. Our aim was to understand the origin of vortex shedding. The critical velocity is lower than for the 2D problem. There is a drag for all velocity, it increases smoothly with the velocity, and there a significant increase at the onset of vortex shedding.

###### Acknowledgements.

The authors would like to acknowledge discussions with Vincent Hakim and Marc Etienne Brachet. Qiang Du is supported in part by a NSF grant DMS-0196522.