text:
# New prospects for polarized hyperon fragmentation functions

## 1 Introduction

The usual parton $`(f=q,\overline{q},g)`$ distribution functions $`f_H(x,Q^2)`$ of a hadron $`H`$, which are extracted from Deep Inelastic Scattering (DIS) for a space-like energy scale $`(Q^20)`$, have counterparts in the time-like region $`(Q^20)`$, which are the hadron fragmentation functions $`D_f^H(z,Q^2)`$. They represent the probability to find the hadron $`H`$ with a fraction $`z`$ of the momentum of the parent parton $`f`$, at a given value of $`Q^2`$. Like distribution functions, fragmentation functions are universal , that is process-independent, and likewise, their $`Q^2`$ dependence is fully predicted by the QCD evolution equations. They can be measured in various high energy collision processes namely, $`e^+e^{}`$ annihilation, $`e^+e^{}(\gamma ^{},Z)HX`$, semi-inclusive DIS, $`lpl^{}HX`$ and single-inclusive $`pp`$ collisions, $`ppHX`$. So far, the best known fragmentation funtions are $`D_f^\pi `$ and $`D_f^K`$, corresponding to the most copiously produced light mesons $`\pi `$ and $`K`$. Our purpose here is to consider the hyperon fragmentation functions and more specifically to examine the case of $`\mathrm{\Lambda }`$ ($`\overline{\mathrm{\Lambda }}`$). The self-analyzing properties of $`\mathrm{\Lambda }`$ ($`\overline{\mathrm{\Lambda }}`$) make this spin-half baryon particularly appealing, to study spin transfer mechanisms. We first recall the results of a QCD analysis of the data for inclusive ($`\mathrm{\Lambda }+\overline{\mathrm{\Lambda }}`$) production in $`e^+e^{}`$ collisions, which yields the first simple and reliable parametrization of the unpolarized fragmentation functions $`D_f^{\mathrm{\Lambda },\overline{\mathrm{\Lambda }}}(z,Q^2)`$. The observed longitudinal polarization of the $`\mathrm{\Lambda }`$’s produced at LEP on the $`Z`$-resonance, leads to some inaccurate information on the spin-dependent fragmentation functions $`\mathrm{\Delta }D_f^\mathrm{\Lambda }(z,Q^2)`$. As we will see, several theoretical models have been proposed for these polarized fragmentation functions which are, so far, badly constrained by existing data. Some predictions can be made for the spin transfer in polarized DIS, but one gets no definite conclusion by comparing them with the present very poor data from HERMES at DESY and E665 at FNAL. We also stress the importance of the $`\mathrm{\Lambda }`$ ($`\overline{\mathrm{\Lambda }}`$) production in neutrino (antineutrino) DIS, which allows a clean flavor and spin separation. We will show some new data from NOMAD at CERN and comment on future experimental possibilities.
We will also give the prospects from $`pp`$ collisions with polarized protons at RHIC-BNL, because there are recent interesting suggestions for measuring the helicity (and transversity) transfer asymmetry in the process $`p\stackrel{}{p}\stackrel{}{\mathrm{\Lambda }}X`$. From its dependence on the rapidity of the $`\mathrm{\Lambda }`$, it is possible to discriminate easily between the various theoretical models, thanks to the high luminosity and the small statistical errors. Of course, this can be easely extended to the production of the other hyperons ($`\mathrm{\Sigma }^{\pm ,0},\mathrm{\Xi }^{,0}`$), which will be also briefly discussed.

## 2  $`e^+e^{}`$ collisions

From the production of $`\mathrm{\Lambda }`$ ($`\overline{\mathrm{\Lambda }}`$) in $`e^+e^{}`$ annihilation one can determine $`D_f^\mathrm{\Lambda }(z,Q^2)`$ by using a standard parametrization, at initial scale $`\mu ^2`$

$$D_f^\mathrm{\Lambda }(z,\mu ^2)=N_fz^{\alpha _f}(1z)^{\beta _f}.$$
(1)
To reduce the number of free parameters occuring above, one must make some simplifying assumptions based on, e.g. $`SU_f(3)`$ arguments. Ten free parameters, needed to include light quarks, heavy quarks and gluons, were obtained by fitting the world data for 14GeV$`\sqrt{s}`$91.2GeV, in a leading and next-to-leading order (LO and NLO) QCD analysis . Although, it leads to a good determination of $`D_f^\mathrm{\Lambda }`$, there is no flavor separation for the dominant contributions $`D_q^\mathrm{\Lambda }(q=u,d,s)`$ and no separation between $`D_q^\mathrm{\Lambda }`$ and $`D_{\overline{q}}^\mathrm{\Lambda }`$, since the data does not separate $`\mathrm{\Lambda }`$ and ($`\overline{\mathrm{\Lambda }}`$). Other sets of quark fragmentation functions into $`\mathrm{\Lambda }`$ ($`\overline{\mathrm{\Lambda }}`$), based on a bag model calculation and a fit to $`e^+e^{}`$ data, have been published recently .

The fragmentation of a longitudinally polarized parton into a longitudinally polarized $`\mathrm{\Lambda }`$ will be described in terms of

$$\mathrm{\Delta }D_f^\mathrm{\Lambda }(z,Q^2)=D_{f(+)}^{\mathrm{\Lambda }(+)}(z,Q^2)D_{f(+)}^{\mathrm{\Lambda }()}(z,Q^2),$$
(2)
where $`D_{f(+)}^{\mathrm{\Lambda }(+)}(z,Q^2)`$ \[$`D_{f(+)}^{\mathrm{\Lambda }()}(z,Q^2)`$\] is the probability to find a $`\mathrm{\Lambda }`$ with positive \[negative\] helicity in a parton $`f`$ with positive helicity. Clearly one obtains the unpolarized fragmentation function $`D_f^\mathrm{\Lambda }`$ by taking in (2), the sum instead of the difference. For the fragmentation of a transversely polarized parton into a transversely polarized $`\mathrm{\Lambda }`$, one uses (2) to define $`\mathrm{\Delta }_TD_f^\mathrm{\Lambda }(z,Q^2)`$ in complete analogy with $`\mathrm{\Delta }D_f^\mathrm{\Lambda }(z,Q^2)`$.

In order to construct $`\mathrm{\Delta }D_f^\mathrm{\Lambda }`$, one usually makes the assumption

$$\mathrm{\Delta }D_f^\mathrm{\Lambda }(z)=C_f^\mathrm{\Lambda }(z)D_f^\mathrm{\Lambda }(z),$$
(3)
at a given scale, where $`C_f^\mathrm{\Lambda }(z)`$ are the spin transfer coefficients, for which several models can be found in the literature. In Ref., $`\mathrm{\Delta }D_f^\mathrm{\Lambda }`$ was neglected for heavy quarks and gluon and they take the following simple parametrization

$$\mathrm{\Delta }D_u^\mathrm{\Lambda }=\mathrm{\Delta }D_d^\mathrm{\Lambda }=N_u\mathrm{\Delta }D_s^\mathrm{\Lambda }\text{and}\mathrm{\Delta }D_s^\mathrm{\Lambda }=z^\alpha D_s,$$
(4)
in terms of two new parameters $`\alpha `$ and $`N_u`$ at a given scale. To cover a wide range of plausible models, they consider three different scenarios:
1) In the non-relativistic quark model, only strange quarks can contribute, therefore $`N_u=0`$.
2) By transposing the estimate of Ref. for the polarized $`u`$ and $`d`$ quark distributions in a $`\mathrm{\Lambda }`$, to the polarized fragmentation functions, one gets $`N_u=0.2`$.
3) By assuming that all polarized fragmentation functions are equal, which is an extreme case, one gets $`N_u=1`$.
In Ref., only two scenarios similar to 1) and 2), were considered and for other alternatives along the same lines, see also Ref.. A different way to construct $`\mathrm{\Delta }D_q^\mathrm{\Lambda }(z)`$ was proposed in Ref., based on our knowledge of the quark distributions and by means of the simple reciprocity relation , the so called Gribov-Lipatov relation

$$D_q^H(z)q_H(x),$$
(5)
where $`q_H(x)`$ is the quark distribution for finding the quark $`q`$ carrying a momentum fraction $`x=z`$ inside the hadron $`H`$.
The longitudinal $`\mathrm{\Lambda }`$ polarization $`P_\mathrm{\Lambda }`$, which can be directly expressed in terms of the $`\mathrm{\Delta }D_f^\mathrm{\Lambda }{}_{}{}^{}s`$ (see ref.), has been measured at LEP by several experiments and is shown in Fig.1. It is large and negative and consistent with the SU(6) quark-diquark model proposed in Ref., as well as scenario 3) from Ref., mentioned above. In a very recent work , one finds various predictions for $`P_B`$, where $`B`$ stands for any hyperon of the octet baryon, obtained by using SU(3) symmetry relations.

## 3 Charged lepton DIS process

These results can be used to make predictions for $`\mathrm{\Lambda }`$ production in semi-inclusive DIS, with only the electron (muon) beam longitudinally polarized. The $`\mathrm{\Lambda }`$ polarization $`P_\mathrm{\Lambda }`$ along its own momentum is related ( for the exact expression see Ref.) to the longitudinal spin transfer $`A^\mathrm{\Lambda }(x,z)`$ defined as

$$A^\mathrm{\Lambda }(x,z)=\frac{\underset{q}{}e_q^2q^N(x)\mathrm{\Delta }D_q^\mathrm{\Lambda }(z)+(q\overline{q})}{_qe_q^2q^N(x)D_q^\mathrm{\Lambda }(z)+(q\overline{q})}.$$
(6)
From the above values of $`N_u`$, one anticipates that $`A^\mathrm{\Lambda }`$ will be small for 1), negative for 2) and large and positive for 3) and the results can be found in Ref.. Another set of predictions is given in Ref. and it will be interesting to compare them with data from future polarized DIS experiments, in particular from HERMES at DESY and COMPASS at CERN.

For the time being, we only have the few data points displayed in Fig.2 for the spin transfer for $`\mathrm{\Lambda }`$ production and in Fig.3 for $`\overline{\mathrm{\Lambda }}`$ production. The three curves represent the predictions from three possible scenarios proposed in Ref. in the framework of the pQCD based model and the use of the reciprocity relation Eq.(5). The dotted curves correspond to pure valence quark contributions and for the dashed curves one has introduced a symmetric quark-antiquark sea, whereas it is asymmetric for the solid curves. Although this last case seems slightly favored, the simultaneous description of $`\mathrm{\Lambda }`$ and $`\overline{\mathrm{\Lambda }}`$ spin transfers remains hard to achieve and therefore more accurate data are badly needed (for a detailed discussion on this point, see Ref. ).

## 4 Neutrino DIS process

The main observation is based on the fact that neutrinos (antineutrinos) can be regarded as a purely polarized lepton beam since neutrinos are left-handed (antineutrinos are right-handed). Therefore they only interact with quarks of specific helicities and flavors. The scattering of a neutrino (antineutrino) beam on a hadronic target provides a source of polarized quarks with specific flavor. This important property makes neutrino (antineutrino) DIS, an ideal laboratory to study the flavor dependence of quark fragmentation.
If one measures, $`\mathrm{\Lambda }`$ and $`\overline{\mathrm{\Lambda }}`$ production in neutrino and antineutrino DIS, one can show  that the four longitudinal polarizations $`P_\nu ^\mathrm{\Lambda }`$, $`P_{\overline{\nu }}^\mathrm{\Lambda }`$, $`P_\nu ^{\overline{\mathrm{\Lambda }}}`$ and $`P_{\overline{\nu }}^{\overline{\mathrm{\Lambda }}}`$ can be expressed in terms of four fragmentation functions, allowing a clean separation between quarks and antiquarks. In Fig.4 we show the NOMAD data for $`P_\nu ^\mathrm{\Lambda }`$, compared with some model calculations from Ref., as well as the predictions for the remaining three processes. Again more data are requiered to test these models but fortunately, one can expect future experimental possibilities from charged currents $`e^\pm p`$ DIS, with HERA at DESY, CEBAF at Jefferson Lab. and eRHIC at BNL.

## 5 Hyperons in $`pp`$ collisions

Let us now move to another area where we can test these polarized fragmentation functions, by studying the spin transfer parameters in $`pp`$ collisions. A new polarized $`pp`$ collider at RHIC-BNL, which starts operating now, will allow to undertake a vast spin physics programme at high center-of-mass energies, up to $`\sqrt{s}`$ = 500GeV. First, we are interested in the reaction $`\stackrel{}{p}p\stackrel{}{\mathrm{\Lambda }}X`$, where both the initial proton and the produced $`\mathrm{\Lambda }`$ are longitudinaly polarized and we consider the helicity transfer parameter $`D_{LL}^\mathrm{\Lambda }`$ defined as follows

$$D_{LL}^\mathrm{\Lambda }(\eta ,p_T)=(d\sigma _{++}d\sigma _+)/(d\sigma _{++}+d\sigma _+),$$
(7)
where $`\eta `$ and $`p_T`$ are the rapidity and the transverse momentum of the outgoing $`\mathrm{\Lambda }`$. $`d\sigma _{++}`$ \[$`d\sigma _+`$\] denotes the cross section where the proton and the $`\mathrm{\Lambda }`$ helicities have the same \[opposite\] sign. When $`\eta `$ is positive the $`\mathrm{\Lambda }`$ has the direction of $`\stackrel{}{p}`$ and $`p_T`$ is always assumed to be large. $`D_{LL}^\mathrm{\Lambda }`$ is directly related to the $`\mathrm{\Delta }D_f^\mathrm{\Lambda }{}_{}{}^{}s`$, whose exact expression can be found in Ref., and therefore it is possible to make predictions, which are shown in Fig.5. Again one can see that the sign and magnitude of $`D_{LL}^\mathrm{\Lambda }`$ are strongly correlated to the values of $`N_u`$ and given the smallness of the expected statistical errors, indicated in the figure along the horizontal axis, it will be now rather easy to pin down the right underlying mechanism.
One notices that in all cases $`D_{LL}^\mathrm{\Lambda }`$ is small for negative rapidity. This is due to the fact that the dominant subprocess is $`gqgq`$ and we have assumed that $`\mathrm{\Delta }D_g^\mathrm{\Lambda }=0`$. This feature is rather general and it was discussed in more details in Ref.. It characterizes also other helicity transfers as shown in Fig.6, where we display the predicted helicity transfers for all the hyperons $`\mathrm{\Lambda },\mathrm{\Sigma }^{\pm ,0},\mathrm{\Xi }^{,0}`$, resulting from two models of the fragmentation functions.

Finally, we can also consider the case where both the proton and the $`\mathrm{\Lambda }`$ are transversely polarized and the corresponding spin transfer parameter $`D_{NN}^\mathrm{\Lambda }`$, is defined as in Eq.(7), where helicity is replaced by transverse polarization. This spin observable is related to the quark transversity distributions, so called $`h_1^q(x,Q^2)`$ and to the transversity fragmentation functions $`\mathrm{\Delta }_TD_f^\mathrm{\Lambda }(z,Q^2)`$ mentioned in Section 2. Needless to say that we have no experimental information on these quantities. However we can derive some bounds on $`D_{NN}^\mathrm{\Lambda }`$, using positivity arguments. We recall that positivity leads to the inequality

$$2|h_1^q(x,Q^2)|q(x,Q^2)+\mathrm{\Delta }q(x,Q^2),$$
(8)
whose validity has been established up to NLO QCD corrections. A similar constraint holds for hadron fragmentation functions ,namely

$$2|\mathrm{\Delta }_TD_q^H(z,Q^2)|D_q^H(z,Q^2)+\mathrm{\Delta }D_q^H(z,Q^2).$$
(9)
By saturating the above inequalities, we get an estimate for an upper and lower bounds for $`D_{NN}^\mathrm{\Lambda }`$ which are displayed in Fig.7. The larger allowed region corresponds to scenario 3), mentioned earlier, due to the value of $`N_u`$ and from the expected statistical accuracy, we get rather small error bars.

To summarize, we have now some knowledge of the unpolarized fragmentation functions $`D_f^\mathrm{\Lambda }(z,Q^2)`$, but we miss a better flavor separation, which can be properly achieved by means of neutrino (antineutrino) DIS or charged currents in $`e^\pm p`$ DIS. The corresponding polarized fragmentation functions $`\mathrm{\Delta }D_f^\mathrm{\Lambda }(z,Q^2)`$ are poorly known and more data is needed for testing several theoretical models proposed in the literature. Future measurements of $`\mathrm{\Lambda }`$ and $`\overline{\mathrm{\Lambda }}`$ production in semi-inclusive charged lepton DIS will improve the situation and this can be generalized to the case of the production of the other hyperons $`\mathrm{\Sigma }^{\pm ,0},\mathrm{\Xi }^{,0}`$. The new polarized $`pp`$ collider at RHIC-BNL, with high luminosity and high energy, is also a very powerful facility to study unpolarized and polarized fragmentation functions for all octet baryon members.

Acknowledgements

It is my pleasure to thank the orginazers S.D.Bass, A. De Roeck and A. Deshpande for their invitation and for setting up this excellent workshop in such pleasant and stimulating atmosphere.


format: text

text:
# Detection of Galaxy Spin Alignments in the PSCz Shear Field

## 1 INTRODUCTION

In the standard gravitational instability paradigm, the linear tidal torque theory explains the origin of the angular momentum of a galaxy as arising from its tidal coupling with the surrounding matter at the early proto-galactic stage (Hoyle, 1949; Peebles, 1969; Doroshkevich, 1970; White, 1984), and predicts a correlation of its spin axis with the intermediate principal axis of local shear (Dubinski, 1992; Catelan & Theuns, 1996; Lee & Pen, 2000). In this model, each component of the proto-galaxy angular momentum vector in the shear principal axis frame is represented by a difference of a pair of the three eigenvalues of the shear tensor, and thus the intermediate component proportional to the difference between the largest and the smallest eigenvalues of the shear tensor dominates the direction of the angular momentum. This correlation of proto-galaxy spin axis with the initial tidal shear field is intrinsic and local, and distinct from the weak gravitational lensing effect Crittenden et al. (2001).

An important question is whether the intrinsic alignments of proto-galaxy spins with the tidal shear fields remain till present or not. If so, it will have a considerable impact on weak lensing searches since the intrinsic shear effect on galaxy spins will play a role of systematic errors. Furthermore recently Lee & Pen (2000,2001) have shown that the relatively easily measurable galaxy spin field can be used to reconstruct the initial tidal shear and mass density fields provided that the effect of intrinsic shear on galaxy spins persist. Several numerical tests (Heavens et al., 2000; Porciani et al., 2001, 2001) and observational searches (Pen et al., 2000; Brown et al., 2000) have so far found some evidences for the existence of intrinsic spin alignments, but early attempts were all indirect.

The detection of intrinsic correlations between the tidal shear and galaxy spin fields from the observed universe requires independent measurement of the two fields. Here we attempt to measure the tidal shear and galaxy spin fields from two independent all-sky catalogs and measure the local directional correlations of the two fields directly.

## 2 FIELD MEASUREMENTS

One would like to measure the intrinsic shear in the initial Lagrangian field, which is, however, practically impossible to derive from present observational data. What one can do at most is to measure the evolved final Eulerian shear and approximate the Lagrangian shear to it. This approximation can be improved by applying a Wiener filter to the observational data.

For the measurement of shear, we have used the PSCz (a complete redshift survey from the IRAS Point Source Catalog) data, which provides a flux-limited all-sky catalog of 15500 galaxies detected by IRAS. It covers $`84\%`$ of the sky down to a flux limit of $`0.6`$ Jy at $`60\mu `$m. The remaining $`16\%`$ of the sky is a gap region containing the galactic plane and survey gaps excluded in the PSCz Saunders et al. (2000).

The flux-limited catalog is not a homogeneous sample since the fraction of galaxies observed decreases with distance. For the measurement of tidal shear, it is statistically simpler to use a volume-limited sample where the fraction of galaxies is constant with distance. We converted the flux limited PSCz catalog into volume-limited samples by dividing the whole volume into the 7 spherical shells of distance range $`[0,50],[50,75],[75,100],[100,125]`$,$`[125,150]`$,
$`[150,175]`$ and $`[175,200]`$ (in unit of $`h^1`$Mpc) as in the previous literature (Yahil et al., 1001; Szapudi et al., 2000) and including only galaxies luminous enough to be above the flux limit if placed at the outer edge of each volume. We ended up with 7 volume-limited samples containing $`1730,1541,1356`$, $`1001,701,586`$ and $`442`$ galaxies respectively.

We first measure the raw density field $`\delta (𝐱)`$ of a chosen PSCz sample (say, $`V`$), assuming that the mass density is proportional to the number density of galaxies: $`\delta (𝐱)=\frac{1}{\rho _0}_{𝐱_iV}\delta _D(𝐱𝐱_i)1`$, where $`\delta _D`$ is the Dirac-delta function, and $`\rho _0`$ is the mean number density of galaxies of the sample $`V`$. We calculate the raw density field of each sample separately as if there were no mass source outside the sample under consideration.

To minimize the shot noise due to the finite number of mass sources and optimize the approximation of the Lagrangian shear to the Eulerian shear, we convolve this density field with the Wiener filter, based on the standard power law correlation function for the galaxies. We find the Fourier space Wiener filter to be $`W(k)=\frac{Ak^{1.2}}{Ak^{1.2}+\rho _0^1}`$ with $`A=0.903\times 2\pi ^2\times (4.43)^{1.8}`$. Here $`\rho _0^1`$ is the noise power spectrum derived from the Poissonian distribution. Note that each sample has a different Wiener filter distinguished by a different value of $`\rho _0`$.

Since the Wiener filtered gravitational potential in Fourier space is $`\varphi (𝐤)=k^2\delta (𝐤)W(k)`$, its real counter-part is written as $`\varphi (𝐱)=_V\delta (𝐱^{})\stackrel{~}{W}(𝐱𝐱^{})d^3𝐱^{}`$ where $`\stackrel{~}{W}`$ is a modified Wiener filter defined in Fourier space as $`\stackrel{~}{W}(k)=k^2W(k)`$.

Now we find the gravitational potential due to the mass density of the chosen sample $`V`$ at some field point $`𝐱`$:

$`\varphi (𝐱)`$ $`=`$ $`{\displaystyle \frac{1}{2\pi ^2\rho _0}}{\displaystyle \underset{𝐱_iV}{}}{\displaystyle k\stackrel{~}{W}(k)\frac{\mathrm{sin}(k|𝐱𝐱_i|)}{|𝐱𝐱_i|}𝑑k},`$ (1)
$`+{\displaystyle \frac{2}{\pi }}{\displaystyle \stackrel{~}{W}(k)\left\{\frac{R_1}{k}\mathrm{cos}(R_1k)\frac{R_2}{k}\mathrm{cos}(R_2k)+\frac{1}{k^2}\mathrm{sin}(R_2k)\frac{1}{k^2}\mathrm{sin}(R_1k)\right\}𝑑k},`$
where $`R_1`$ and $`R_2`$ are the inner and outer radii of $`V`$. The total gravitational potential at $`𝐱`$ will be the super-positions of all the gravitational potentials calculated from the seven PSCz samples separately. The tidal shear field can be finally determined as the second derivative of the total gravitational potential such that $`T_{ij}(𝐱)=_{\mathrm{samples}}_i_j\varphi (𝐱)`$.

There are, however, two sources of bias in the measurement of tidal shear: the existence of a gap region excluded in the PSCz and the redshift distortion effect on the orientation of the shear principal axes along the line of sight. To minimize these errors, first we compensate for the existence of the gap by adding the gap contribution to the total mass density, assuming a uniform mass distribution in the gap. Second, we cancel out the redshift distortion bias in the determination of shear by readjusting the shear principal axes to have a statistically isotropic distribution. Basically we force the angle between the shear principal axes of each shear and the line-of-sight direction to have a uniform distribution in the range of $`[0,2\pi ]`$ without changing its original ranking order. The most important intermediate axis has been adjusted first, the major axis as next, and finally the minor axis has been automatically determined from the orthogonality condition.

Note that these procedures minimize the systematic bias in the measurement of shear orientation but also lower the inferred strength of intrinsic shear effects since the shear axes may deviate further from the true orientations. Finally we renormalize the resulting shears into the unit traceless shear since our eventual interest is not in the magnitude but the axis orientation of the tidal shear.

Regarding the galaxy spin field, we consider the unit spins of 12122 Tully disk galaxies as used by Pen et al. (2001). Let us briefly review the work of Pen et al.. The unit spin vector of each disk galaxy is initially estimated from the major-minor axis ratio and the position angle, assuming that a disk galaxy is an infinitely thin disk with its spin axis perpendicular to the disk plane. The measured axis ratio and position angle of a disk galaxy determine the inclination and orientation angles of the disk. The inclination angle is only determined up to a sign ambiguity. We take into account this sign ambiguity in the determination of line-of-sight component of a unit spin, and construct two unit spin vectors at each Tully galaxy position. Therefore we end up with twice as many unit spins as Tully galaxies.

Real disk galaxies, however, are not infinitely thin disks. Their intrinsic thickness and deviation in shape from a perfect ellipse create biases in the statistical distribution of inclination angles. Although this statistical bias has a different source from the redshift distortions, its effect on the galaxy spin axis can be removed by the same procedure that we use for the adjustment of shear principal axis. We adjust $`\alpha `$ to have a uniform distribution in range of $`[0,\pi ]`$, as would be the unbiased case.

This statistical bias might depend on both the galaxy distance and Hubble type (intrinsic morphology). We divide the Tully galaxies into the same distance bins, resulting in five non-zero samples. The number of the included galaxies, the distance range of each Tully sample are listed in the 1st and 2nd columns of Table 1. For each Tully sample, we carry out the uniform redistribution of galaxy inclination angles within each Hubble type separately.

Now that we have independently measured tidal shear and galaxy spin fields, we can estimate the strength of intrinsic correlations between the tidal shear and galaxy spin fields.

## 3 MEASUREMENTS OF INTRINSIC SHEAR EFFECTS

Lee & Pen (2000, 2001) have shown that the intrinsic shear effect on galaxy spins can be quantified by one dimensionless parameter, $`a`$ in range of $`[0,0.6]`$, and derived its optimal estimation formula: $`a=26_{i=1}^3\widehat{\lambda }_i^2\widehat{L}_i^2`$ where $`\{\widehat{\lambda }_i\}`$ are the three eigenvalues of shear, and $`(\widehat{L}_i)`$ is the unit spin vector measured in the shear principal axis frame. The associated statistical uncertainty, $`\sigma `$, can be calculated as the standard deviation of $`a`$ assuming the null hypothesis that the orientations of galaxy spins are completely random with no intrinsic alignment: $`\sigma =2/\sqrt{5N}`$ where $`N`$ is the total number of independent unit spin vectors Lee & Pen (2001). Each galaxy is used twice, once for each degenerate orientation. We have verified numerically that using $`N`$ given by twice the number of galaxies correctedly reproduces the standard deviation for an uncorrelated sample.

We calculate the values of $`a`$ at the positions of all Tully galaxies using this formula, and averaged them within each sample separately. Table 1 lists the observed signal $`a^{obs}`$ of each sample in the 3rd column. We have also checked the effect of measurement errors in the axis ratios and position angles of Tully galaxies on the observed signal of intrinsic shear effect. Given the average mean errors of axis ratio and position angle de Vaucouleurs et al. (1991), we have calculated the propagation error due to measurement inaccuracy and found it to be small, less than $`10\%`$ of the signal for every sample. Fig. 1 illustrates a statistical sample of intrinsic shear effect on the Tully unit spins in the distance range $`[98100]`$ $`h^1`$ Mpc.

The observed signal does not give the full strength of intrinsic shear effect. In our analysis the shear field has been convolved by the adaptive Wiener-filter to minimize noise, but the full strength of intrinsic shear effect should be measured from the shear field smoothed with a top-hat galactic radius. The observed signal, $`a^{obs}`$, can be rescaled to the full signal, $`a`$, by $`a=a^{obs}[\mathrm{cov}(\delta _w,\delta _{th})]^2`$. Here $`\mathrm{cov}(\delta _w,\delta _{th})`$ is the covariance between wiener filtered and top-hat filtered density fields ($`\delta _w`$,$`\delta _{th}`$ respectively). Approximating the density power spectrum on galactic scales as a power law $`P(k)=k^2`$, and using a top-hat radius of $`0.55`$$`h^1`$ Mpc which corresponds to the typical size of a galaxy halo Bardeen et al. (1986), we rescale the observed signals from each sample separately and find the full signal of intrinsic shear effect, which are listed in the 4th column of Table 1. This results from each Tully sample except for the 1st one are consistent with one another, detecting a nonzero signal of a few percent. Therefore, the null hypothesis that the orientations of galaxy spin axes in the local shear frame are completely random with no intrinsic alignment can be safely abandoned.

In the null hypothesis, the distribution of the correlation parameter for each sample is Gaussian with zero mean. The non-zero value of $`a_i`$ ($`i=1,\mathrm{},5`$) of the Tully sample would just be statistical fluctuations with standard deviation $`\sigma _i`$. Treating each sample as independent, we define five independent variables, $`A_1,\mathrm{},A_5`$ as $`A_i_{a_i}^{\mathrm{}}𝑑a\frac{1}{\sqrt{2\pi }\sigma _i}\mathrm{exp}(\frac{a^2}{2\sigma _i^2})`$. In the null hypothesis each $`A_i`$ is a random variable, uniformly distributed in the range $`[0,1]`$. The probability of the null hypothesis, that is, the probability that the Tully sample doe not detect intrinsic alignments of galaxy spins is given as the cumulative probability, $`P(<W)`$ with $`W\mathrm{\Pi }_{i=1}^5A_i`$. One can show that this probability is given as $`P(W)=[\mathrm{log}(W)]^4/4!`$. From the results of Table 1, we find $`W=3.25\times 10^8`$ and $`P(<W)=1.53\times 10^4`$, which rules out the null hypothesis at the confidence level higher than $`99.98\%`$.

Finally we compute the weighted average of the signals by $`\overline{a}=\frac{_i(a_i/\sigma _i^2)}{_i(1/\sigma _i^2)}`$ and $`\overline{\sigma }=[_i(1/\sigma _i^2)]^{1/2}`$. Excluding the inconsistent result from the 1st sample, and we find $`\overline{a}=0.17\pm 0.04`$, which is a 25% accurate measurement of the intrinsic shear alignment effect.

## 4 DISCUSSIONS AND CONCLUSIONS

We have presented an observational evidence for the existence of intrinsic shear effect on present-day galaxy spin orientations. We have used current observational data (the PSCz and the Tully galaxy catalogs) to measure the tidal shear and the galaxy spin fields independently, and measured the strength of intrinsic shear effect directly. We have used five samples in different distant range and measured the signal of each sample separately. All the samples except for the innermost one have produced non-negligible signal of intrinsic shear effects on galaxy spins, and thus we have been able to rule out the null hypothesis that galaxy spins are randomly aligned in the tidal shear field at a confidence level higher than $`99.98\%`$. The mean strength of the signal averaged over the four samples except for the innermost one has been shown to be $`\overline{a}=0.17\pm 0.04`$. The reconstruction procedure is likely to systematically underestimate $`\overline{a}`$.

Yet, there are two statistical problems in the measurement of $`\overline{a}`$, which we would like to address before concluding. The first concern is about excluding the 1st sample. When combining the samples that have high signal-to-noise ratio, we should give each sample a weight proportional to its volume. In such a weighting, the 1st Tully sample has the smallest weight. Furthermore, structures are coherent on scales of many Mpc, so that the “cosmic variance”is expected to be large in the small sample. Our estimator for $`a`$ assumed that the shear is known without any errors, i.e. we neglected“cosmic variance”. This first sample has a small statistical error, but it is not clear how to include this sample when we in fact know that the shear field is not known to nearly such good accuracy. Cosmic shear has a correlation length comparable to the size of this bin, so the coherent errors in this bin are expected to be large. Hence, given the potentially large contributions of cosmic variance, we exclude it in the calculation of $`\overline{a}`$.

Excluding the 1st sample, however, is the simplest way to treat the effects of cosmic variance on $`\overline{a}`$. For a more realistic treatment, one should include the cosmic variance error, assuming a particular relation between the sample size and the value of $`a`$, and use a maximum likelihood analysis for $`\overline{a}`$. The second concern is about the assumed statistical model for $`\overline{a}`$. Here we have assumed a null hypothesis of $`\overline{a}`$ equal zero, and have shown that alignments are not random at the $`99.98\%`$ level. If cosmic variance is indeed the origin for the small apparent signal in the first bin, the estimate of $`\overline{a}`$ has an additional uncertainty from cosmic variance. The quantitative estimate of the cosmic scatter, however, is beyond the scope of this letter.

We finally conclude that the null hypothesis of no intrinsic alignment is very safely ruled out, and the resulting signal represents the average strength of intrinsic shear effect from the galaxy samples that have consistent nonzero signals. Our result is consistent with the linear tidal torque model based on the gravitational instability. It also hints that the present galaxy spin field might be used as a probe to measure the initial tidal shear and mass density fields, and eventually map the three dimensional matter distribution of the universe (Lee & Pen, 2001). Applying the results of (Crittenden et al., 2001), we conclude that weak lensing surveys are strongly contaminated by this effect at low redshift, but only negligibly at high redshift.

We thank B. Tully and E. Branchini for their catalogs and useful comments. We also thank the referee for helping us improve the manuscript. This work has been supported by the Academia Sinica and NSERC.


text:
# Gas flows around two young stellar clusters in NGC2264

## 1 Introduction

The mass of the molecular ISM is dominated by the largest clouds which tend to form stars in groups, rather than in isolation (Lada et al., 1994). Star formation is inefficient, however, and as the binding mass of a molecular cloud is lost through the action of stellar winds, radiation, and outflows, most groups evaporate (Adams, 2000) and become isolated stellar systems. Thus an increased understanding of clustered star formation is essential for learning about the origin of most Galactic disk stars, perhaps including the Sun.

In this paper, we present millimeter wavelength observations of two star formation sites, IRS1 and IRS2, in the NGC2264 region of northern Monoceros. Lying at a distance of 760 pc (Sung, Bessell, & Lee, 1997), this region is associated with a moderate mass cloud, $`M3\times 10^4M_{}`$ (Oliver, Masheder & Thaddeus, 1996), that contains 30 IRAS sources, most classified as Class I protostars (Margulis, Lada, & Young, 1989), and $`360`$ near-infrared sources (Lada et al., 1993).

IRS1, also known as Allen’s source (Allen, 1972), is an early B star situated just north of the Cone nebula. Thompson et al. (1998) postulate that it has triggered the formation of several low mass stars around it and Ward-Thompson et al. (2000) map a number of nearby massive dust condensations suggesting that high mass star formation may also be occuring. Molecular spectral lines show signatures of both outflow (Margulis, Lada, & Snell, 1988; Schreyer et al., 1997) and infall (Wolf-Chase & Gregersen, 1997; Williams & Myers, 1999).

IRS2, discovered by Sargent et al. (1984), is a less luminous and less well studied source $`6^{}`$ to the north of IRS1. A small cluster and nebulosity is visible in the I band image of Lada et al. (1993). Margulis, Lada, & Snell (1988) found a CO outflow toward this source, Wolf-Chase, Walker, & Lada (1995) mapped an associated dense molecular clump (which they called MonOB1-D), and Williams & Myers (1999) found preliminary evidence for collapse.

Here, we combine continuum observations of the dusty envelopes around the most embedded protostars with spectral line observations of the gas in these two cluster forming regions. The molecular envelopes around the two clusters have a similar mass but are morphologically distinct and present an interesting comparative study. In both cases, the self-absorbed HCO+ line is used to diagnose relative motions and show large scale collapse onto each stellar group that appears to be disrupted by protostellar outflows on smaller scales. In the bright, compact IRS1 clump, the optically thin H13CO+ line is used to map the molecular abundance, velocity centroid, and dispersion and thereby constrain the modeling of the HCO+ profiles. In IRS2, where individual protostars are resolved, the size and velocity dispersion is used to estimate an upper limit for their formation time and to determine the gravitational stability of the group. The observations are described in §2, and their analysis in §3. We conclude in §4.

## 2 Observations

The observations were carried out at the Heinrich Hertz Telescope111the Heinrich Hertz Telescope is operated by the Submillimeter Telescope Observatory on behalf of Steward Observatory and the Max-Planck-Institut fuer Radioastronomie. on Mt. Graham, AZ over a four day period in April 2001. We used the facility 19 channel bolometer to observe continuum emission at $`870\mu `$m and the 230 GHz receiver to observe the $`J=32`$ lines of HCO+ and H13CO+. The weather was good during the run and the zenith opacity at $`870\mu `$m ranged from 0.3 to 0.55; for the spectral line observations system temperatures ranged from 350 to 450 K.

The bolometer observations were carried out in dual-beam raster mode with an azimuth scanning velocity of $`8^{\prime \prime }`$/sec and a wobbler frequency of 2 Hz. Individual maps had sizes $`8^{}\times 6^{}`$, sampled at $`10^{\prime \prime }`$ in azimuth and elevation, and took 40 minutes to complete. Skydips were made before and after each map to determine the sky opacity and observations of Jupiter and Uranus were used to calibrate the data. Two maps were made of each source at different hour angles and therefore different scan directions across the source. The data were reduced and co-added using the NIC software package resulting in an rms noise in the final maps of 0.14 Jy per $`21^{\prime \prime }`$ beam.

The spectral line data were also taken in on-the-fly mode. Maps were made by scanning alternately in right ascension and declination at $`3^{\prime \prime }`$/sec with the reference $`12^{}`$ west of the map center. Individual map sizes were $`3^{}\times 2^{}`$, sampled at $`10^{\prime \prime }`$, and took $`25`$ minutes each. System temperatures varied from 350 K to 425 K. Maps were repeated several times and made at different offsets, then co-added and stitched together in the CLASS software package so as to achieve an rms noise temperature of 0.8 K (HCO+) and 0.4 K (H13CO+) per 250 kHz ($`0.28`$ $`\mathrm{km}\mathrm{s}^1`$) channel. The temperature scale was calibrated through the chopper wheel method (Rohlfs & Wilson, 1996), and we then applied main beam and forward efficiencies $`\eta _{\mathrm{mb}}=0.46,\eta _{\mathrm{fss}}=0.92`$. Observations of the standard source Orion IRc2 were made daily as an overall check of calibration, pointing, and setup. For both the bolometer and line observations, pointing was checked through observations of Jupiter every 2–3 hours.

## 3 Results

### 3.1 Bolometer data

The $`870\mu `$m continuum emission toward IRS1 and IRS2 is shown in Figure 1. IRS1 is the bright, compact source to the south. Its elongated shape shows signs of substructure; Ward-Thompson et al. (2000) resolve five distinct condensations with higher resolution IRAM and JCMT observations. IRS2 lies to the north and is more fragmented with several low intensity peaks, suggestive of a slightly more evolved cluster in which the individual protostars are beginning to shed their circumstellar envelopes and migrate away from their birthplace.

The integrated fluxes, for flux densities greater than 0.5 Jy beam-1, of the two clusters are similar, 94 Jy for IRS1 (consistent with Ward-Thompson et al.) and 84 Jy for IRS2 indicating that they have similar dust masses and make a good comparative study. The fluxes convert to masses through

$$M=\frac{F_\nu d^2}{\kappa B_\nu (T_d)}=8.4\left(\frac{F_\nu }{\mathrm{Jy}}\right)M_{}$$

where we have taken the mass opacity $`\kappa =0.009`$ cm2 g-1, dust temperature $`T_d=17`$ K, as in Ward-Thompson et al. (2000) and distance $`d=760`$ pc. Thus the total interstellar mass around IRS1 is $`790M_{}`$ and $`710M_{}`$ around IRS2. From the projected area of each group, 0.45 pc2 for IRS1 and 0.62 pc2 for IRS2, we estimate average column and volume densities $`N_{\mathrm{H}_2}=8\times 10^{22}`$ cm-2, $`n_{\mathrm{H}_2}=5\times 10^4`$ $`\mathrm{cm}^3`$ for IRS1 and $`N_{\mathrm{H}_2}=5\times 10^{22}`$ cm-2, $`n_{\mathrm{H}_2}=3\times 10^4`$ $`\mathrm{cm}^3`$ for IRS2.

### 3.2 Line data

The average density of the material traced by the continuum observations is similar to that of HCO+ emitting gas. Indeed, the morphology of the integrated line emission follows the dust (see Figure 2) and we therefore use the line data to investigate the dynamics of the circumcluster material.

The optically thick HCO+ spectra provide information on relative motions, i.e. infall and outflow, in the gas (e.g., Walker et al. (1986)) and the optically thin H13CO+ show the intensity weighted mean velocity and linewidth along each line of sight. We also use the latter to measure the molecular abundance in the IRS1 clump.

The individual protostars in the compact IRS1 clump merge together in the $`29^{\prime \prime }`$ beam of the line observations and we examine average, group properties. Objects in the more extended IRS2 clump are resolved, however, and we can determine mass infall rates onto both the cluster and a single protostar. In addition, by comparing the cluster’s size and velocity dispersion, we estimate its age.

#### 3.2.1 IRS1

Line and continuum data in IRS1 are compared in Figure 2. The left panel overlays contours of integrated H13CO+ intensity over a grayscale of the dust emission. Their similar appearance demonstrates that we can use the line data to examine motions of the circumcluster material. The right panel plots the two quantities on a point by point basis and shows an approximately linear trend over a decade in scale with the implication that the H13CO+ is both optically thin and its abundance constant across the clump. By converting each quantity to a column density (assuming an excitation temperature of 20 K for the gas and dust temperature of 17 K) the best fit H13CO+ abundance relative to H2 is found to be $`x(\text{H}\text{13}\text{CO}\text{+})=2.0\times 10^{11}`$. For an isotopic ratio, $`{}_{}{}^{12}\mathrm{C}/^{13}\mathrm{C}=50100`$, the implied HCO+ abundance lies in the range, $`x(\text{HCO}\text{+})=12\times 10^9`$, similar to values in other OB star forming regions (Bergin et al., 1997; Ungerechts et al., 1997). The lack of variation across the clump is demonstrated by the two dashed lines bracketing the data at $`x(\text{H}\text{13}\text{CO}\text{+})=[1.0,4.0]\times 10^{11}`$.

Higher moments of the H13CO+ emission show the variation of mean velocity and dispersion across the cluster (Figure 3) The former increases from 7.0 $`\mathrm{km}\mathrm{s}^1`$ toward the East to 8.4 $`\mathrm{km}\mathrm{s}^1`$ at center, but then returns to a lower value, $`7.8`$ $`\mathrm{km}\mathrm{s}^1`$ in the West. There is no clearly defined rotation across the clump and toward the central region, $`\mathrm{\Delta }\alpha =10^{\prime \prime }`$ to $`50^{\prime \prime }`$, the velocity difference, $`0.3\mathrm{km}\mathrm{s}^1`$, is much less than the dispersion. The dispersion map itself also shows little variation in this central region and has an average value, $`\sigma =1.0`$ $`\mathrm{km}\mathrm{s}^1`$.

The HCO+ spectra are self-absorbed across the clump: they generally show two peaks and a central dip at the mean velocity of H13CO+. Changes in the velocity of the absorbing gas relative to the emitting gas alter the observed profile and are a diagnostic of relative motions in the clump (Leung & Brown, 1977). The relatively poor resolution of the line data, $`29^{\prime \prime }0.1`$ pc, does not permit a detailed analysis of the dynamics around individual protostars within this cluster and we therefore focus on global properties, as illustrated in Figure 4. Here, average spectra are shown, computed over different regions of the clump corresponding to different $`870\mu `$m flux densities. This is similar to radial averaging but takes into account the non-circular shape of the clump. At low flux densities, corresponding to the outer regions, the average spectra are relatively weak and the optical depths are low. Nevertheless, the HCO+ spectra are clearly asymmetric and show the red-shifted self-absorption characteristic of inward motions (Walker et al., 1986). At higher flux densities, toward the clump center, the line intensities and optical depths increase and the HCO+ profiles show two peaks with a central dip. However, the spectra become more symmetric with increasing flux density until at the highest values at the clump’s center, the two HCO+ peaks are of similar intensity suggesting zero average inflow in this region.

Infall spectra are notoriously difficult to interpret (Menten et al., 1987), and the change in profile shapes observed in IRS1 require detailed modeling to understand the flows of material in this complex region. Since the clump has a single, central peak and no clear velocity gradient, we use the 1-d radiative transfer code, RATRAN, of Hogerheijde & van der Tak (2000) to model, simultaneously, the dust and line data. The continuum and H13CO+ observations fix the size, temperature, density profile, H13CO+ abundance, mean velocity, and velocity dispersion and tightly constrain the model, leaving only the velocity profile and isotopic ratio as the remaining parameters in the model output HCO+ line profiles.

The model consists of 10 logarithmically spaced, radially symmetric, shells from $`r=10^{16}`$ cm to $`r=r_{\mathrm{max}}=10^{18}`$ cm (0.32 pc). Ward-Thompson et al. (2000) show that the spectral energy distribution of IRS1 is well fit by the sum of two greybodies consisting of a compact component at 38 K and a larger, cooler component at 17 K that dominates the mass. Based on this, we set the temperature of the outer layer at 17 K and the inner layers, $`r<6\times 10^{17}`$ cm at 38 K. The gas temperature is fixed to be the same as the dust as would be expected at the high densities in the model. Figure 2 shows that the clump is moderately, but not strongly, peaked and we find that the density profile $`n_{\mathrm{H}_2}(r)=10^5(r_{\mathrm{max}}/r)^{3/2}`$ $`\mathrm{cm}^3`$, when convolved with a $`29^{\prime \prime }`$ beam, matches the dust and velocity integrated H13CO+ maps well. The abundance, $`x(\text{H}\text{13}\text{CO}\text{+})=2.0\times 10^{11}`$, systemic velocity, $`v=8.0`$ $`\mathrm{km}\mathrm{s}^1`$, and velocity dispersion, $`\sigma =0.8`$ $`\mathrm{km}\mathrm{s}^1`$, are set constant in all layers as indicated by the H13CO+ data. The high relative motions between the layers contribute to broaden the output H13CO+ spectra close to the observed 1.0 $`\mathrm{km}\mathrm{s}^1`$. With these parameters, we find that we can reproduce the (radially averaged) continuum fluxes and H13CO+ spectra to within 20%.

Under these constraints, it is not possible to fit the observed HCO+ spectra as closely since only the isotopic abundance and relative velocities between the layers can be varied. Through experimentation with the models, it was found that, given the observed density, temperature, abundance, and velocity profiles, the increasing symmetry of the HCO+ spectra toward the cluster center can only be explained by an outer layer collapsing onto an expanding inner region. Figure 4 plots radially averaged model HCO+ and H13CO+ spectra as dotted lines in comparison with the observed data for the case where the isotopic ratio, $`{}_{}{}^{12}\mathrm{C}/^{13}\mathrm{C}=75`$ (implying $`x(\text{HCO}\text{+})=1.5\times 10^9`$), and the outer layer ($`r=6\times 10^{17}`$ to $`1\times 10^{18}`$ cm) collapses at a speed $`v_{\mathrm{in}}=0.3`$ $`\mathrm{km}\mathrm{s}^1`$ onto the inner layers ($`r=10^{16}`$ cm to $`6\times 10^{17}`$ cm) expanding at a speed $`v_{\mathrm{exp}}=1.0`$ $`\mathrm{km}\mathrm{s}^1`$.

The combination of large scale collapse and central expansion conspire to match the observed spectral variations from core edge to center because the optical depth is very high in the low density, outer layers and produces a broad (blue and redshifted), self-absorption dip along lines of sight toward the center but only redshifted self-absorption along lines of sight near the edges that do not intercept the inner expanding layer. The velocity of the absorption dip is not observed to change from core center to edge and is reproduced in the model via the increased projection of the radial flow along lines of sight toward the edges. The model spectra have broader self-absorption dips than the data because of the small number of layers in the model and the abrupt change in temperature from 17 K in the outer layer to 38 K in inner layers. Applying a smoother temperature transition and more layers produces more rounded self-absorption dips but at the expense of adding additional parameters to the model.

Because of the radial dependence, $`M(<r)r^{3/2}`$, most of the mass of the core is in the outer layers and altering the velocities of the innermost layers, $`r10^{17}`$ cm, does not change the output profiles noticeably. That is, the model resembles the two layer model described by Myers et al. (1996) but with a spherical geometry. With that caveat, the output HCO+ profiles are very sensitive to the optical depth and infall speed of the outer layer and expansion speed of the layer interior to that (which together account for about three quarters of the total mass) and through experimentation with the models we found that the column density and velocities of these layers are constrained to within $`20\%`$ by the requirement of matching the observed redshifted absorption at large clump radii and near-symmetry at small radii.

The outer, infalling, layer contains 1/2 the total mass, or about $`M_{\mathrm{in}}=400M_{}`$, and the mass infall rate onto the cluster is therefore

$$\dot{M}_{\mathrm{in}}=M_{\mathrm{in}}v_{\mathrm{in}}/r_{\mathrm{max}}=4\times 10^4M_{}\mathrm{yr}^1.$$

This is one to two orders of magnitude greater than infall rates onto solar and intermediate mass protostars (Zhou, 1995) as expected for a global infall onto a group of such objects and is similar to the value expected for a gravitational flow, $`\dot{M}=\sigma ^3/G`$, where $`\sigma =1.0`$ $`\mathrm{km}\mathrm{s}^1`$ is the average velocity dispersion of the core (Figure 3), and $`G`$ is the gravitational constant (Shu, 1977).

The model requires an expanding shell to fit the increasing symmetry toward the core center. Physically, this is probably due to the protostellar outflows in the region, one of which appears to be directed along our line of sight (Schreyer et al., 1997). Shu, Adams, & Liazno (1987) hypothesize that (low mass) stars self-determine their final mass through the removal of a collapsing envelope by their outflows (see also Velusamy & Langer (1998)). By comparing clusters in different evolutionary states and examining the gas flows around them, it will be possible to see if such a process is relevant in a clustered environment where most stars form, and therefore assess its role in the origin of the stellar IMF (Adams & Fatuzzo, 1996).

#### 3.2.2 IRS2

Continuum and line emission were lower in IRS2 than in IRS1. The proportion between them, however, was approximately the same and the H13CO+ abundance was found to be consistent with that measured in IRS1, though with a greater uncertainty due to the lower signal-to-noise ratio. Due to the relative weakness of the emission, only the central part of the cluster was mapped in HCO+ and H13CO+. Figure 5 plots the boundary of the line data on a map of continuum emission. Slight red-shifted self-absorption is found in the average HCO+ spectrum shown in the left panel. The self-absorption dip is small due to the low optical depth when averaged over the entire mapped region, but is more pronounced in individual spectra, particularly around the dust peaks where the column density is greatest. Both red- and blue-shifted asymmetries are found (see below) but the dominant tendency is for red-shifted asymmetry suggesting, as in IRS1, large scale collapse onto the cluster.

Given the complex resolved structure in IRS2, similar fits to the spectra at different size-scales as in §3.2.1, are not feasible with a 1-d model. However, both simple plane-parallel and spherical two layer model fits to the average spectrum imply an infall speed $`v_{\mathrm{in}}0.3\pm 0.1`$ $`\mathrm{km}\mathrm{s}^1`$. This is similar to the speed of the collapsing layer around IRS1 and implies a similar mass infall rate, $`\dot{M}_{\mathrm{in}}4\times 10^4M_{}\mathrm{yr}^1`$.

The greater separation of the protostars in IRS2 allows an investigation of the flows around each one individually. Seven cores, labeled a-g in Figure 5 are apparent as peaks in the dust emission and interpolated spectra (at a resolution of $`29^{\prime \prime }`$) toward the six that lie within the boundary of the line data are shown in Figure 6. A complex mixture of dynamical states are found. Core c has the classic infall profile and core d shows a strong reversed signature. Cores e and f show possible signs of infall and reversed flow respectively but the optical depth is too low to be sure. The HCO+ spectra toward cores a and b are not readily interpreted in the same manner, due to multiple peaks possibly from absorption by overlapping outflows.

Gaussian fits to the H13CO+ spectra show that the non-thermal velocity dispersion is typically in the range 0.7 to 1.0 $`\mathrm{km}\mathrm{s}^1`$ but has a pronounced local minimum of 0.4 $`\mathrm{km}\mathrm{s}^1`$ around core b. The resolution of these data is too low to test collapse models in detail but we do not find a correlation between the HCO+ asymmetry and the H13CO+ dispersion, analogous to that seen in the Serpens cluster (Williams & Myers, 2000) and as might be expected from motions driven by the dissipation of turbulence.

To estimate the dynamical influence of each protostar, we have examined average spectra around each one with different smoothing lengths. Although there is confusion from merging at angular scales $`45^{\prime \prime }`$, we find that the overall infall dynamics of the cluster eventually dominate and that the influence of any single protostar does not exceed a radius $`60^{\prime \prime }0.2`$ pc.

Core c shows the clearest infall signature with a deep, red-shifted, absorption dip. The infall spectra extend for $`60^{\prime \prime }`$ away from the core before merging with the reversed profiles of spectra associated with core d. The core also has the largest envelope with a total continuum flux of 12.3 Jy and integrated H13CO+ intensity of 2.7 K $`\mathrm{km}\mathrm{s}^1`$ within a radius $`r=0.17`$ pc, implying a mass $`M_{\mathrm{env}}=100M_{}`$. The high column density through this envelope is responsible for the deep self-absorption dip, and does not imply a particularly high infall speed. A model fit using the spherically symmetric RATRAN code implies an infall speed $`v_{\mathrm{in}}=0.3`$ $`\mathrm{km}\mathrm{s}^1`$, and mass infall rate $`\dot{M}_{\mathrm{in}}=2\times 10^5M_{}\mathrm{yr}^1`$. The speed is the same as the collapse around the entire cluster but the mass infall rate is considerably less because the fit requires only an outer $`10M_{}`$ layer to be collapsing. Dividing the overall mass infall rate onto the cluster by the number of protostars within it gives an average mass infall rate per protostar about 3 times higher than measured around core c, the one with the clearest individual collapse signature: clearly not all the material falling onto the cluster makes it all the way to the scale of protostars.

The one other core, d, with a deep self-absorption dip shows the opposite asymmetry, indicative of an expanding, rather than infalling, outer layer. It also has the largest H13CO+ linewidth, 1.0 $`\mathrm{km}\mathrm{s}^1`$, in the region. Using a similar model as for core c, the expansion speed is estimated to be 0.5 $`\mathrm{km}\mathrm{s}^1`$. This is likely due to an embedded star although there is no IRAS or MSX point source at its center: higher sensitivity mid-infrared observations are required to determine the stellar content of each core.

Because the envelopes around individual protostars in the cluster are resolved, the optically thin H13CO+ spectra give a measure of the core-core velocity dispersion in the cluster. The six cores, a-f, have a dispersion $`\sigma _{\mathrm{c}\mathrm{c}}=0.90`$ $`\mathrm{km}\mathrm{s}^1`$ and extend over an area with radius $`r=0.44`$ pc. Together, these imply a virial mass, $`M_{\mathrm{virial}}=3r\sigma _{\mathrm{c}\mathrm{c}}^2/G=250M_{}`$. This is less than the total interstellar mass estimated from the continuum observations but is likely much greater than the final stellar mass of the system (Lada et al., 1993). Thus, although the system is currently bound, as the surrounding gas is dispersed through, e.g., stellar outflows, it is likely to become unbound (Lada, Margulis, & Dearborn, 1984).

An unbound cluster expands and its age may be estimated from its size and expansion rate. This effect is seen in optically visible clusters in Orion over a range of ages $`312`$ Myr (Blaauw, 1991; Brown, Dekker, & de Zeeuw, 1997) but the morphology of the two groups IRS1,2 in Figure 1 suggests that it may apply to considerably earlier times. In IRS2, if $`\sigma _{\mathrm{c}\mathrm{c}}`$ is due to the expansion of the protostellar group, the formation time can be estimated from its crossing time, $`t_{\mathrm{cross}}=r/\sigma _{\mathrm{c}\mathrm{c}}=5\times 10^5`$ yr, characteristic of late Class I sources (Mundy, Looney, & Welch, 2000). The gravitational attraction of the gas around the stellar group should slow their expansion, and the crossing time is therefore an upper limit to the cluster age. Given the continued infall onto the cluster and the large envelope mass, it is not clear if the cluster is actually expanding and whether the crossing time is a useful measure of the cluster age at these early stages in its formation but the suggestive morphology and flux densities of the continuum maps are compelling. Additional observations of other embedded stellar groups in this and other regions will show if there is a consistent pattern. In this regard, we note that a very short integration toward a more extended (and therefore possibly more evolved) group of IRAS sources $`11^{}`$ north of IRS2 did not detect $`870\mu `$m emission to an rms flux level of 0.1 Jy beam-1.

## 4 Conclusion

We have mapped two embedded stellar groups, IRS1 and IRS2, in the NGC2264 region in the continuum at $`870\mu `$m and the $`J=32`$ lines of HCO+ and H13CO+. The former shows the thermal emission from the dusty envelopes around the clusters and their protostellar members and the latter show the flows around them; the optically thick HCO+ diagnosing relative motions (infall and outflow) and the optically thin H13CO+ revealing the bulk motion of the clusters and the protostellar velocity dispersion within IRS2.

The two groups have a similar total mass, $`M10^3M_{}`$, but IRS1 has a much higher peak flux density and is more compact than the more dispersed, lower surface brightness IRS2. In IRS1, comparing the continuum flux density with the integrated line intensity shows that the H13CO+ abundance is constant across the clump at the $`29^{\prime \prime }`$ resolution of these data. Average HCO+ spectra are fit with a radially symmetric model constrained by the dust and H13CO+ data and show a large scale collapse with an infall speed $`v_{\mathrm{in}}=0.3`$ $`\mathrm{km}\mathrm{s}^1`$ and mass infall rate $`\dot{M}_{\mathrm{in}}=4\times 10^4M_{}\mathrm{yr}^1`$, falling onto an inner core expanding at speed $`v_{\mathrm{exp}}=1.0`$ $`\mathrm{km}\mathrm{s}^1`$. More detailed observations of the competing processes of infall and outflow will show its role in determining stellar masses.

Large scale collapse with a similar infall speed and mass infall rate are found in IRS2 but detailed modeling is precluded by the complex structure that is resolved within. Flows around individual protostars show a mixture of dynamical states ranging from localized collapse with the same speed as the large scale collapse and mass infall rate, $`\dot{M}_{\mathrm{in}}=2\times 10^5M_{}\mathrm{yr}^1`$, to expansion at speeds $`v_{\mathrm{out}}0.5`$ $`\mathrm{km}\mathrm{s}^1`$. By averaging spectra over increasingly larger size scales, the influence of any one protostar was found to be limited to within a radius $`0.2`$ pc.

The H13CO+ spectra in IRS2 was used to measure the protostellar velocity dispersion and thereby infer the virial mass of the system. This was found to be less than the total mass enveloping the cluster but more than the stellar mass indicating that the group is destined to become unbound as the surrounding gas and dust are dispersed by protostellar winds and outflows. The expansion of an unbound stellar group may already have started in IRS2 where the greater separation of individual protostars and the low column density toward them, relative to IRS1, both suggest a more evolved system. In this case, an upper limit to its age, $`5\times 10^5`$ yr, can be determined kinematically by assuming a constant expansion speed equal to the observed protostellar velocity dispersion. Additional observations of other clusters are necessary to see if there is a pattern between cluster sizes, velocity dispersion, and evolutionary state, but this technique – a well known tool for optical clusters – may be a useful new approach for following protocluster and protostellar evolution.

The dominant mode of stellar birth is in groups and these observations provide new insight into their formation. To examine cluster formation at earlier times will require higher resolution observations of dense, compact cores such as IRS1, while large scale, high sensitivity imaging of the dust and gas around more evolved systems than IRS2 can follow the evolution at later times than discussed here.

We thank the scientific staff at the HHT for their assistance during the observation run, Michiel Hogerheijde for his advice in setting up and using the RATRAN radiative transfer program, and Joe McMullin, the referee, for his helpful comments and suggestions.


text:
# Tracing protons through the Galactic magnetic field: a clue for charge composition of ultra-high energy cosmic rays.

## I Introduction

Ultra-high energy cosmic rays (UHECR) are a subject of an active debate for over more than 20 years. The data accumulated during this time experiments ; takeda  provide a compelling evidence that the GZK cutoff GZK  predicted in the UHECR spectrum at energies $`E10^{20}`$ eV may be absent. Many models explaining this puzzle have been proposed (for reviews see Refs. reviews ). The actual focus of the debate is the question of whether an unconventional astrophysical model can be constructed which explains the observed super-GZK events, or a new physics is required. This key question still remains open.

With the accumulation of the high-energy events important signatures emerge which discriminate efficiently between different models. It has been known for some time that UHECR events form clusters clusters1 ; takeda ; Uchihori:2000gu . Recent analysis shows that clustering is statistically significant Tinyakov:2001ic . The small angular size of clusters of order $`2.5^{}`$, consistent with the experimental angular resolution, suggests that they are due to point sources. The models which do not reproduce this feature are therefore disfavored.

Observation of clusters implies that some sources of UHECR may be identified. Recently, significant correlations of arrival directions of UHECR with positions of BL Lacertae (BL Lacs) were found Tinyakov:2001nr . BL Lacs are blazars (ANGs with relativistic jet directed along the line of sight) characterized in particular by the (near) absence of the emission lines. The correlations between UHECR and BL Lacs are most significant at angles of order $`2.5^{}`$ and are present at relatively low energies $`E>2.4\times 10^{19}`$. Such tightness suggests strongly the existence of neutral primary particles. Association with distant BL Lacs combined with neutrality of primaries rules out most of the models of UHECR, leaving the models based on neutrino Zbur  (the Z-burst models), models with hypothetical “immune messengers” hm , or models involving violation of the Lorentz invariance VLI . It should be noted that ultra-high energy photons also cannot be excluded at present Kalashev:2001qp .

Regardless of whether there exist neutral primary particles, correlations with BL Lacs imply that the acceleration mechanism of UHECR production actually works. Therefore, protons are involved and should be present in the UHECR flux at energies around or below the GZK cutoff. If identified, such protons provide independent evidence that BL Lacs are indeed sources of UHECR. The purpose of this paper is to address this issue.

Knowing actual sources of UHECR is an extremely powerful tool even when statistics is limited. As we will see shortly, this tool can be used to study charge composition of UHECR. The idea is to use the bending of charged primary particles in the Galactic (GMF) GMF\_bend  and extragalactic (EGMF) EGMF\_bend  magnetic fields. While deflections in EGMF are random and therefore unpredictable, deflections in GMF are regular. If charges, energies of particles and GMF are known, the original directions (before entering GMF) can be restored. If the effect of EGMF is small, the actual positions of sources can be reconstructed.

In practice, one has to solve the inverse problem: to reconstruct charges assuming set of potential sources. We show that this is possible at least in a statistical sense: assuming that a substantial fraction of UHECR are protons significantly improves correlations with BL Lacs. This implies  simultaneously that GMF model is roughly correct, and that the effect of EGMF is small.

The paper is organized as follows. In Sect. II we briefly review the present knowledge about the Galactic magnetic field and fix the GMF parameters for further calculations. In Sect. III we analyse correlations between BL Lacs and UHECR when the latter are assigned non-zero charges. Sect. IV contains discussion and concluding remarks.

## II Galactic Magnetic Field

The Galactic magnetic field can be divided in two parts: the disc and the halo. Each one has regular and turbulent components. While the strength of turbulent component is larger, it is the regular field which gives dominant contribution into CR propagation. Most of the information on the regular component of the disc is obtained from the Faraday rotation measurements of pulsars and extragalactic radio sources. The latter are used also for the reconstruction of the halo field. Magnetic field in the disc resembles the spiral structure of our Galaxy. It may either reverse direction between different spiral arms (bisymmetric, or BSS model), or there may be no reversals (axisymmetric, or ASS model). Several field reversals were detected RL ; HQ ; ID ; HMQ ; Frick  which are consistent with BSS model (note however that discrimination between the two models is complicated by small-scale irregularities in the magnetic field). In our calculations we adopt the BSS model. Simple analytical representation of the spiral field structure (see Refs. SF ; HQ ) contains the following parameters:

* Distance from the Sun to the Galactic center, $`R=8.5`$ kpc.
* Local (at the Sun position in the Galaxy) field strength, $`B_0`$.
* Pitch angle $`p`$ which determines the direction of the local magnetic field (the field points in the direction $`l=90+p`$ in the Galactic coordinates).
* Distance $`d`$ to the first field reversal. Negative $`d`$ means that nearest reversal occurs in the direction to the Galactic center, positive corresponds to the opposite direction.

In terms of these parameters the field in the disk is written as

$$B_\theta =B\mathrm{cos}(p),B_r=B\mathrm{sin}(p).$$

The magnitude $`B=B(r,\theta )`$ has the spiral structure,

$$B(r,\theta )=B(r)\mathrm{cos}\left(\theta \beta \mathrm{ln}\left(\frac{r}{R}\right)+\varphi \right),$$
(1)
where $`\beta 1/\mathrm{tan}(p)`$, a constant phase $`\varphi `$ is given by

$$\varphi =\beta \mathrm{ln}\left(1+\frac{d}{R}\right)\frac{\pi }{2},$$
(2)
and

$$B(r)=B_0\frac{R}{r\mathrm{cos}(\varphi )}.$$
(3)
In the last expression the standard assumption that magnitude of the field decreases as $`r^1`$ in radial direction is made. It is also assumed that $`B(r)=\mathrm{const}`$ at $`r<4`$ kpc. Note that precise dependence of the disc field far away from the Sun position is not important for our study as only a small fraction of the observed UHECR passes through this region. Note also that the constant phase $`\varphi `$ can be absorbed in another parameter $`r_0`$, so that Eq. (1) becomes $`B(r,\theta )=B(r)\mathrm{cos}(\theta \beta \mathrm{ln}(r/r_0))`$, which is the parametrization used in Refs. SF ; HQ . We find it more convenient to work with parameters directly related to the local field.

To proceed, we need to fix the parameters $`B_0`$, $`p`$ and $`d`$. All studies of GMF based on pulsar and extragalactic rotation measures converge on $`B_0=1.4\mu \mathrm{G}`$, see RL ; HQ ; ID ; Frick  and recent reviews beck ; han . We adopt this value.

Pitch angle was found to be close to zero or positive in early publications, $`p=+5^{}`$ RK , $`p=2^{}`$ RL , but decreased in more recent studies: $`p=8^{}`$ was obtained in Refs. HQ ; Heiles ; ID ; HMQ , while in Ref. Frick  $`p=15^{}`$ was found as an average pitch angle in nearby spiral arms. Following reviews beck ; han  we take $`p=8^{}`$.

Field reversal was found to be at $`d=0.20.3`$ kpc in Ref. HQ  (however, the best fit value of $`r_0`$ in the same paper corresponds to $`d=0.48`$), at $`d=0.4`$ kpc in Ref. RL , at $`d=0.6`$ kpc in Refs. RK ; Vale ; Frick . Finally, in the review beck  $`d=0.6`$ kpc was cited, while the review han  follows Ref. HQ . We take $`d=0.5`$ kpc.

The simplest approximation for the halo field is obtained by taking the disk field and extending it outside of the disk with exponentially decreasing amplitude,

$$B(r,\theta ,z)=\mathrm{exp}\left(\frac{|z|}{h}\right)B(r,\theta ).$$
(4)
This introduces one more parameter, the height $`h`$. Small values of $`B_z=\mathrm{\hspace{0.17em}0.2}\mu \mathrm{G}`$ found in HQ  for the halo field support this approximation. The disc has a height of $`h=1.5`$ kpc according to reviews beck ; han  and we adopt this value.

In addition, the halo fields above and below the disk (more precisely, its parts parallel to the disk) may be parallel or anti-parallel. In the latter case the halo field may be approximated by

$$B(r,\theta ,z)=\mathrm{sign}(z)\mathrm{exp}\left(\frac{|z|}{h}\right)B(r,\theta ).$$
(5)
The first case, Eq. (4), corresponds to the quadrupole-type model which we denote $`\mathrm{BSS}_\mathrm{Q}`$, while the second case, Eq. (5), corresponds to the dipole-type model denoted as $`\mathrm{BSS}_\mathrm{D}`$ in what follows. There are indications in favour of the $`\mathrm{BSS}_\mathrm{D}`$ global structure HMBB ; beck ; han  of the halo field.

In calculations of UHECR propagation in GMF, the parameters of the latter are often chosen following early work Stanev:1996qj . We do not use the conventions of Ref. Stanev:1996qj  because they are ambiguous: the parameter $`r_0`$ used there does not correspond to the cited local field strength, the value of pitch angle does not correspond to cited direction of the local field and references for the assumed value of the halo height are not given.

To summarize, we adopt Eqs. (1-5) for the Galactic magnetic field with the following set of parameters

$`B_0=1.4\mu \mathrm{G}p=8^{}`$
$`d=0.5\mathrm{kpc}h=1.5\mathrm{kpc}`$ (6)
and assume that this field extends to $`R_{\mathrm{max}}=20`$ kpc in all directions. It should be stressed that these parameters are chosen on the basis of rotation measurements and their choice has nothing to do with the propagation of UHECR.

## III Correlation analysis

We start by specifying the sets of BL Lacs and UHECR. Following Ref. Tinyakov:2001nr , we take BL Lacs form the QSO catalog veron  which contains 306 confirmed BL Lacs. The choice of UHECR set is motivated as follows. Apriori, best results should be achieved with the largest set having best angular and energy resolution. In addition, this should be (relatively) high-energy set, because the uncertainties in the deflection angle (which result from uncertainties in GMF and energies of particles) should not dominate over the angular resolution. This requirement is satisfied at $`E>4\times 10^{19}`$ eV. Therefore, we chose all published AGASA events with energy $`E>4\times 10^{19}`$ eV takeda . This set contains 57 events.

For the quantitative measure of correlations we use the probability $`p(\delta )`$ introduced in Refs. Tinyakov:2001ic ; Tinyakov:2001nr , which is a probability to have certain excess of events within the angle $`\delta `$ from any of the sources (BL Lacs). This probability is calculated by counting how often the excess observed in the real data occurs in the Monte-Carlo (MC) simulations. The MC configurations are generated as described in Refs. Tinyakov:2001ic ; Tinyakov:2001nr . The difference is that now charges are assigned and arrival directions are corrected for GMF both in the real data and in each MC set. The energies of MC events are taken from the real data. Exactly the same treatment of the real data and MC sets is important to prevent appearance of artificial correlations.

The charge assignment can be done by any algorithm as long as the same algorithm is used in MC simulations. In the simplest case one assigns equal charges to all particles. Although correlations which are due to neutral particles are destroyed in this case, this effect may be compensated by charged particles which move closer to their sources. One may expect this situation when fractions of neutral and charged particles are comparable.

As explained in Ref. Tinyakov:2001nr , there are two possible strategies to estimate the significance of correlations on the basis of $`p(\delta )`$. One may impose cuts on BL Lac set and adjust them in such a way that correlations are maximum. Likewise, one may look for the ”best” values of the parameters of GMF. In this case the penalty factor should be calculated and included in the probability. Since reconstruction of particle trajectories in GMF is rather time-consuming, this approach is not very practical in the case at hand. Thus, we do not adjust the parameters of GMF, Eq. (6), and impose no cuts on BL Lac set except for a single cut in the apparent magnitude. Instead of penalty calculation we present explicit dependence of the probability $`p(\delta )`$ on this cut. We draw our main conclusions from the behaviour of the curve, rather than from the value of the probability at the minimum.

Correlations found in Ref. Tinyakov:2001nr  were strongest at $`\delta 2.5^{}`$. We therefore fix this value of $`\delta `$ in our calculations. Since $`\delta `$ is not adjusted to minimize the probability, there is no penalty factor associated with that.

Consider first the case of the symmetric (quadrupole) model $`\mathrm{BSS}_\mathrm{Q}`$. Fig. 1 shows the dependence of $`p(2.5^{})`$ on the cut in magnitude (the cut $`\mathrm{mag}<20`$ corresponds to inclusion of practically all BL Lacs). Solid and dotted lines correspond to the charge $`+1`$ or $`1`$ assigned to all particles, respectively.

We see that in both cases $`p(2.5^{})`$ has minima which are comparable in depth (although both not very deep). One is tempted to conclude that both charges may be present. This kind of situation is expected in the Z-burst models. One may notice, however, that these minima are not equally significant. The minimum at $`Q=+1`$ is wider and corresponds to much higher statistics (correlations are present at the level of $`1\%`$ even when all BL Lacs are included). Moreover, event-by-event analysis shows that 12 out of 14 events contributing to the minimum of probability at $`Q=+1`$ are situated in the Northern hemisphere. This strange feature suggests that the field in the Southern hemisphere is wrong, while in the Northern one it is roughly correct. The obvious thing to try is the $`\mathrm{BSS}_\mathrm{D}`$ model in which the Southern field has different sign.

Consider the case of the asymmetric (dipole) model $`\mathrm{BSS}_\mathrm{D}`$. In our problem changing the direction of the GMF in the Southern hemisphere to the opposite is equivalent to flipping sign of charges of the events in the Southern hemisphere. According to previous results, this should increase the correlations. Indeed, the situation changes. Fig. 2 shows the dependence of the probability $`p(2.5^{})`$ on the cut in magnitude in the case of the $`\mathrm{BSS}_\mathrm{D}`$ model.

Three cases $`Q=1,0,+1`$ are shown. The correlations in the case $`Q=+1`$ have improved by almost two orders of magnitude as compared to the $`\mathrm{BSS}_\mathrm{Q}`$ model and are now at the level below $`10^3`$ in a wide range of magnitudes. Even with no cuts on BL Lacs the significance of correlations at $`Q=+1`$ is below $`10^3`$. On the contrary, in the case $`Q=1`$ the correlations are now absent. In the case $`Q=0`$ correlations with BL Lacs satisfying $`\mathrm{mag}<18`$ are observed at the level of $`2\%`$. For completeness, Fig. 3 shows the the dependence of $`p(\delta )`$ on the angle $`\delta `$ at the cut $`\mathrm{mag}<18`$.

## IV Discussion and conclusions

The pairs BL Lac – cosmic ray separated by less than $`2.5^{}`$ are listed in Table I (imposing the cut $`\mathrm{mag}<18`$). BL Lacs appearing in this table are probable sources of UHECR. For each BL Lac the name, Galactic coordinates l and b, and redshift z (when known) are given. Corresponding cosmic ray energy E (in units of $`10^{19}`$ eV) and probable charge Q are listed in the last two columns.

We have assigned $`Q=0`$ to those particles which contribute to correlations at $`\delta =2.5^{}`$ in the case when all particles are assumed to be neutral. Likewise, we have assigned charge $`Q=+1`$ to those events which fall within $`2.5^{}`$ from any of the BL Lacs when they are assumed to have charge $`+1`$ and deflection in GMF is taken into account. Some events satisfy both requirements; corresponding entry in Table I reads “0 or +1”.

Lines 1-8 of Table I contain four BL Lacs which correlate with two cosmic rays each. These are most probable sources. In the lines 9-16 eight BL Lacs are listed which correlate with singlets. Finally, lines 17-22 contain three cosmic rays for which sources are ambiguous (each ray has two neighbouring BL Lacs within $`\delta <2.5^{}`$).

Examining Table I one notices two striking regularities (we do not attempt to assign statistical significance to those in the present paper). First, most objects in this table are X-ray selected radio loud BL Lacs. Second, the fraction of BL Lacs with unknown redshifts in Table I is much larger than in average over the whole BL Lac catalog. In doublet part of the Table this fraction is 6:2, in singlet section this fraction is 7:1.

Absence of emission lines (more precisely, their weakness and narrow width) is a defining feature of BL Lac family within the general blazar class. Therefore, it is not surprising that redshifts of roughly half of confirmed BL Lacs are not known. Increased fraction of such BL Lacs in Table I may mean that the absence of emission lines is important for a blazar to became emitter of UHECR.

As it is generally assumed, BL Lacs with unknown redshift may be far away: the observational bias would then explain that the redshifts of such objects are unknown more often. Then however new physics or extreme astrophysics would be necessary to explain observed correlations.

Within the conventional framework, when propagation is limited by interactions with comic backgrounds, one should expect sources to be relatively nearby. Closest BL Lac with known redshift is at $`z=0.029`$. This is well outside of the GZK sphere. It should be noted, however, that most of the cosmic rays in Table I have energies below the GZK cutoff. As it was found in Ref. Kalashev:2001qp , the flux of protons at energy $`E=8\times 10^{19}`$ is not attenuated substantially if sources are located at $`z<0.03`$, and the flux at $`E=5\times 10^{19}`$ is not attenuated if sources are located at $`z<0.1`$. Assuming this range of redshifts for charged entries in Table I, we see that there is no problem to explain all of them except lines 3 and 5. The ray (21-22) has sufficiently low energy so that BL Lac B2 0804+35 with z=0.082 can be a source without any problem with attenuation. Likewise, the BL Lac 19 can be an actual source of the ray (19-20) if it has redshift $`z0.1`$. Thus, the correlations due to charged particles listed in Table I can be explained within the framework of conventional physics if sources with unknown redshifts are assumed to be within $`z0.1`$.

Consider now the entries corresponding to neutral particles. According to Ref. Kalashev:2001qp , photons can be UHECR primaries at any energy $`E>10^{19}`$ eV provided that EGMF is smaller than $`10^9`$ G, energy spectrum at the source is hard, $`E^\alpha `$ with $`\alpha <2`$ and the maximum energy is high enough. (Note that correlations with charged particles imply that EGMF is in this range anyway.) Therefore, entries 3 and 17 can be explained by photons. The real difficulty is with rays 5,7,8 and 16. Note, however, that according to MC simulations, there should be 7 background events in the Table I on average.

To summarize, the idea of using the Galactic magnetic field as a mass-spectrograph of UHECR seems to work. The correlations between UHECR and BL Lacs substantially improve when arrival directions of cosmic rays are corrected for GMF. If not a statistical fluctuation, this implies the following: i) cosmic rays of highest energies contain a substantial fraction of protons ii) extragalactic magnetic fields have little effect on propagation of UHECR even from cosmological distances iii) the model of the Galactic magnetic field described by Eq. (5) is roughly correct.

Finally, the significance of correlations with charged particles is $`p<10^3`$. This may be considered as yet more evidence that BL Lacs are sources of UHECR. Interestingly, the events and BL Lacs contributing into the correlation with lowest $`p`$ found in Ref. Tinyakov:2001nr  and the charged events of Table I do not overlap, so the two correlations should be considered as independent.

## Acknowledgments

We are grateful to J. L. Han, M. E. Shaposhnikov, D. V. Semikoz and P. Veron for valuable comments and discussions. The work is supported by the Swiss Science Foundation, grant 21-58947.99 and by INTAS grant 99-1065.


text:
# Coherent electron transport in a Si quantum dot dimer

\[

## Abstract

We show that the coherence of charge transfer through a weakly coupled double-dot dimer can be determined by analyzing the statistics of the conductance pattern, and does not require large phase coherence length in the host material. We present an experimental study of the charge transport through a small Si nanostructure, which contains two quantum dots. The transport through the dimer is shown to be coherent. At the same time, one of the dots is strongly coupled to the leads, and the overall transport is dominated by inelastic co-tunneling processes.

\]

The ability to preserve quantum coherence over large distances and during extended period of time plays a key role in the quest for alternative schemes for conventional electronics and nonclassical electronic behavior. For this purpose, low-dimensional structures, in particular quantum dots, have an obvious advantage, because the $`k`$-space for inelastic scattering events is reduced due to the reduced dimensionality. Conventionally, coherence is probed by quantum interference effects, such as weak localization. In closed quantum dots, coherence has been investigated by embedding the dot in one arm of an Aharonov-Bohm (AB) interferometer. This method requires a host material with large phase coherence length $`l_\varphi `$, larger than the total length of both arms of the interferometer.

Is it possible to measure coherence in a mesoscopic device embedded in a material with small $`l_\varphi `$? We encountered this problem during our studies of Si nanostructures, where in the host two-dimensional electron gas $`l_\varphi <1000`$ Å and interferometric methods cannot be used. Coherence in Si nanostructures is of particular interest because Si has intrinsically long spin relaxation time, which is of great importance for future spintronic devices. In this work we present a new method to discriminate between coherent and incoherent transport through a double dot system. The method is based on statistical analysis of the conductance pattern, and does not rely on large $`l_\varphi `$ in the surrounding contact regions.

Statistical properties of single quantum dots have been extensively investigated over the past ten years. In a weakly coupled double-dot dimer the statistical properties of each dot, such as the peak heights distribution, are almost uncorrelated. However, the way the individual conductances are combined into the total conductance of the dimer depends on whether the transport is coherent or sequential. Thus, a proper deconvolution of the total conductance can identify the type of the transport through the whole structure. This method directly probes the coherence during the charge transfer through several nanostructures, which is of a paramount importance for any practical applications.

In the following, we analyze charge transport through a Si double-dot device. We show that electrons are transferred coherently through the dimer, even though the transport is dominated by inelastic co-tunneling processes. This result was not anticipated a priori, since conventional wisdom associates inelastic processes with decoherence.

The sample is a Si quantum dot fabricated from a silicon-on-insulator wafer. A narrow bridge with a lithographically defined dot is formed from the top Si layer; the bridge is connected to wide source and drain regions via two constrictions. Subsequently, a 50 nm thick layer of SiO2 is thermally grown around the dot, followed by a poly-Si gate. For a detailed description of sample preparation see Ref. . In most cases, more than one dot is naturally formed in such Si nanostructures, and the origin of these additional dots is a subject of ongoing research. For this particular study we have chosen a device that exhibits two distinct periods of conductance oscillations at low temperatures. As will be shown below, the device consist of two dots, both participating in the charge transport.

The conductance $`G`$ through the sample is plotted in Fig. 1a as a function of gate voltage $`0<V_g<9`$ V. The temperature dependence is shown in the inset. At high temperatures, $`15`$ K $`<T<60`$ K, the data is consistent with the conventional theory of Coulomb blockade (CB) in a single dot. This dot will be called dot 1 throughout the paper. At $`T<15`$ K the behavior of the conductance changes qualitatively – oscillations with a much smaller period are superimposed on the main dot CB oscillations. A remarkably large number of these fast oscillations – more than $`500`$ – can be resolved in a single scan. This coexistence of two periods suggests that two dots are involved in the transport. In fact, we can rule out interference effects as an origin of the fast oscillations by showing that there is an electrostatic coupling between the dots.

It is well known that quantum dots can be used as sensitive electrometers, and, in our case, each of the dots can be potentially used to measure the charge on the other dot. In Fig. 2 a normalized peak position of fast oscillations, $`V_p^{^{}}(N)=[V_g^p(N)/\mathrm{\Delta }V_{g2}N]`$, is plotted as a function of $`V_g`$, where $`V_g^p(N)`$ is the position of the $`N`$-th peak and $`\mathrm{\Delta }V_{g2}=14`$ mV is the average peak spacing. For periodic oscillations $`V_p^{^{}}`$ should be a constant, independent of $`N`$ (and $`V_g`$). Vertical lines mark positions of the CB peaks in dot 1. Slips of $`V_p^{^{}}`$ appear every time an electron is added into dot 1. These slips should be expected for capacitively coupled dots: each electron, added into dot 1, increases the electrostatic potential of dot 2 by $`\mathrm{\Delta }\varphi =\frac{eC_c}{C_{\mathrm{\Sigma }1}C_{\mathrm{\Sigma }2}}`$, resulting in a $`\mathrm{\Delta }\varphi \frac{C_{\mathrm{\Sigma }2}}{C_{g2}}`$ shift of the CB peaks. Here, $`C_c`$, $`C_{g2}`$, $`C_{\mathrm{\Sigma }1}`$ and $`C_{\mathrm{\Sigma }2}`$ are the cross capacitance between the dots, the gate capacitance of dot 2, and the total capacitance of dots 1 and 2. The slips are extended over a few periods of the fast oscillations, due to the finite broadening of the CB peaks in dot 1. We also observed two large slips extended over $`20`$ periods (see inset). These probably reflect charging of some other traps, which do not participate in the transport. In principle, slips in peak positions can occur for an interferrometer geometry as well, as a result of the $`\pi `$ phase shift accumulated each time an electron enters the dot. This effect would cause a slip of half a period per an added electron, which is much larger than experimentally observed. Therefore, the existence of the slips, correlated with CB peaks from dot 1, is an unambiguous experimental evidence for the presence of the second dot.

How are the dots coupled? If the dots are strongly coupled electrostatically, $`C_cC_{\mathrm{\Sigma }1},C_{\mathrm{\Sigma }2}`$, the dimer will behave as a single dot and the conductance should demonstrate single period oscillations. This is clearly not the case in our sample. In the opposite regime, $`C_cC_{\mathrm{\Sigma }1},C_{\mathrm{\Sigma }2}`$, one can distinguish between the two possibilities: i) the dots are connected in parallel, and ii) the dots are connected in series. At low gate voltages $`V_g<3`$ V the fast oscillations are suppressed in the valleys of the CB in dot 1 (see left inset of Fig 1a), implying that the dots are connected in series.

In the regime of sequential tunneling through two weakly coupled dot connected in series the total conductance $`G_{\mathrm{seq}}^1G_1^1+G_2^1`$, where $`G_{1,2}`$ are the conductances of each dot. Sequential tunneling is in qualitative disagreement with the data. In such regime amplitude of the fast oscillations $`\delta R`$ of the total resistance $`R_{\mathrm{seq}}=G_{\mathrm{seq}}^1`$ should originate from the CB in the extra dot, $`\delta RG_2^1`$. As shown in Fig. 1b, $`\delta R`$ is strongly correlated with $`R_1`$ and changes by two orders of magnitude within a few periods of the fast oscillations (there are $`8`$ periods of the fast oscillations per slow period). Such strong modulations of the peak height with such small ($`8`$ periods) correlation length is not expected for either weakly or strongly coupled dots.

As Fig 1a shows, the amplitude of the fast oscillations $`\delta G`$ is correlated with the envelope of the total conductance $`G_{\mathrm{env}}`$: the amplitude is larger at the peaks and smaller at the valleys. Indeed, the ratio $`\delta G/G_{\mathrm{env}}`$, plotted in Fig. 1c, is practically $`V_g`$-independent up to 6 V and gradually decreases with further increase of $`V_g`$. Moreover, this ratio, rather than $`\delta G`$ or $`\delta R`$, is almost the same near the minima and the maxima of the main dot CB oscillations. This observation hints that the total conductance should be treated quantum mechanically as a transmission problem. In this case, the total conductance is proportional to the product of the transmission through each dot, $`G_{QM}\mathrm{\Gamma }_{\mathrm{total}}=\mathrm{\Gamma }_1\mathrm{\Gamma }_2`$, and the transport through the whole dimer is coherent.

In the second part of the paper we use our knowledge of the transport through the dimer to analyze some intriguing features in the temperature dependence of the total conductance, and to show that the transport is dominated by inelastic co-tunneling. Parameters of dot 1 can be easily extracted in the usual way. The obtained gate capacitance $`12`$aF is consistent with the geometrical estimates for the lithographically defined dot. The charging energy $`E_{c1}`$, and the mean level spacing $`\mathrm{\Delta }_1`$ are extracted from the statistics of the peak spacings, see left panel of Fig. 3. The spacings form a broad distribution, consistent with $`\mathrm{\Delta }_1E_{c1}4`$meV. The mean level broadening is $`\mathrm{}\mathrm{\Gamma }_10.8`$ meV, less than $`\mathrm{\Delta }_1`$. The peak widths strongly fluctuate at low temperatures $`k_BT<\mathrm{}\mathrm{\Gamma }_1`$, where the width is determined by the coupling to the leads. As expected from the Random Matrix Theory, both distributions are asymmetric.

In Fig. 4 and the right panel of Fig. 3 we present the results of the $`T`$-scaling of peak width and the distributions of peak spacing and width for dot 2, performed similarly to the analysis of the main dot. The standard CB analysis does not work for this dot: i) the peak widths do not scale linearly with temperature; ii) the peak shapes are not Lorentzian, although their width, $`\mathrm{\Gamma }_2`$, saturates at low temperatures at 4 mV, which translates to $`2`$ meV $`k_BT`$; and iii) there is no appreciable fluctuations of the peak width even at low temperatures, where the $`T`$-dependence of each peak has already saturated. Nevertheless, one can analyze the distribution of the peak spacing. The sharpness of the distribution requires the mean level spacing to be much smaller than the charging energy, $`\mathrm{\Delta }_2E_{c2}`$ ($`E_{c2}`$ does not fluctuate). Thus, the following set of inequalities is satisfied $`\mathrm{\Delta }_2,k_BT\mathrm{}\mathrm{\Gamma }_2<E_{c2}`$. This means that dot 2 is in the strong coupling regime, where the standard CB theory is not applicable. The above features of the second dot can be understood within the framework the co-tunneling theory in the strong coupling regime.

There are two contributions to the conductance of a strongly coupled asymmetric dot. One is from elastic co-tunneling (the dot always remains in its ground state) and the other involves inelastic processes, which create particle-hole excitations in the dot. Elastic co-tunneling dominates at low temperatures, resulting in a strongly fluctuating $`G`$ and, correspondingly, peak width. At higher temperatures the leading mechanism of electron transport is inelastic co-tunneling and the conductance shows regular non-fluctuating periodic modulations as a function of $`V_g`$, and gradually evolves with temperature from CB peaks into smooth oscillations. The relevant energy scale is $`E_cr^2\mathrm{cos}^2(\pi 𝒩)`$, where $`r`$ is the smallest of the reflection coefficients at the barriers, and $`𝒩=V_g/eC_{g2}`$ measures the charge that minimizes the dot electrostatic energy, regardless of charge quantization. The conductance at the maxima and minima of the oscillations are estimated to be $`G_{\mathrm{max}}(e^2/\pi \mathrm{})k_BTr^2/(E_cr^2)`$, and $`G_{\mathrm{min}}(e^2/\pi \mathrm{})(k_BT)^2r^2/(E_cr^2)^2`$, respectively. The ratio between the maxima and minima is expected to have a linear $`T`$-dependence

$`G_{\mathrm{min}}/G_{\mathrm{max}}k_BT/(E_cr^2).`$

Indeed, as shown in Fig. 4b, the ratio $`G_{\mathrm{min}}/G_{\mathrm{max}}`$ for the fast oscillations has linear $`T`$-dependence regardless of whether it is measured near the peaks (solid symbols) or the valleys (open symbols) of dot 1 CB oscillations. The fast oscillations are observed up to $`k_BT=E_cr^2`$, thus strong coupling effectively renormalizes the charging energy.

The anomalous temperature dependence of the conductance, observed in the sample, can be naturally understood for the coherent transport through the dimer. Experimentally, the conductance near the peaks of dot 1 increases with temperature (see inset in Fig. 1a), contrary to the prediction of the CB theory that the peak height should be temperature independent for $`k_BT<\mathrm{}\mathrm{\Gamma }_1`$ and $`1/T`$ for $`\mathrm{}\mathrm{\Gamma }_1<k_BT<\mathrm{\Delta }_1<E_{C1}`$ . Because transport is coherent, the total conductance is $`GG_1G_2/(e^2/\pi \mathrm{})`$, and the anomalous temperature dependence is a result of the $`T`$-dependence of $`G_2`$. As we have shown above, $`G_2`$ increases with $`T`$ due to the inelastic co-tunneling processes. Thus, the observed anomalous temperature dependence is an additional argument in favor of both coherent transport through the dimer and inelastic co-tunneling in dot 2.

In conclusion, we proposed a new method to identify coherent transport using two quantum dots, which does not require large phase coherence length in the host material. We studied the electron transport through a Si double-dot structure, and established that the electrons are transferred coherently through the whole dimer. Due to strong coupling to the leads of one of the dots, the transport is dominated by the inelastic co-tunneling processes, altering the conventional temperature dependence of the CB oscillations. In the future, this method can be used to measure $`l_\varphi `$ in the interconnect region by varying the distance between the dots.

This work was supported by the ARO, ONR and DARPA.