text:
# Vacuum diamagnetism of mesoscopic metallic samples

\[

## Abstract

Zero-point energy of surface plasmon modes of mesoscopic metal samples in an external magnetic field has been considered. The magnetic response of plasmon vacuum has shown to be diamagnetic. In thin films this surface vacuum diamagnetism is at least of the same order of magnitude as the magnetism of the bulk electrons. Thus, a novel type of magnetism shown by thin metallic samples, which is complementary to the well known Pauli paramagnetic and Landau diamagnetic contributions to the magnetism of the electron gas has been demonstrated.

PACS no.: 78.67.-n, 74.78.Na

\]

Transport and magnetic properties of mesoscopic systems are the topics of very active current research, due to the strong drive towards development of novel nanodevices and quantum computing. However, there still exist numerous current contradictions between mesoscopic theories and experiment . While a lot of attention is concentrated on the potential role of sample imperfections such as magnetic impurities, at least some of these contradictions may be explained by previously unaccounted intrinsic mechanisms. One of the missing intrinsic elements in many current experimental and theoretical studies is the effect of surface plasmons  on the transport and magnetic properties of mesoscopic metallic samples, such as wires, rings, and various other shapes. Importance of surface plasmons in mesoscopic phenomena is clear from the recent paper  and from the following qualitative consideration.

The usual justification for not considering surface plasmons in low temperature measurements is that surface plasmons are not excited if the sample size is of the order of a few micrometers and the temperature is low. However, even if there are no real plasmon quanta in the system, the zero-point fluctuations of the electromagnetic field of all the possible plasmon modes in the system have to be considered. The importance of zero-point fluctuations in mesoscopic systems clearly manifests itself in the observations of negative vacuum energy density between metal plates separated by submicrometer distances (the Casimir effect, see for example ). The energy density in such a mesoscopic cavity depends on the dielectric constant of the material between the metal plates. A mesoscopic metallic sample constitutes a similar mesoscopic resonator for surface plasmons. As has been shown in numerous papers, magnetic field produces substantial modifications of the surface plasmon dispersion law via modification of the dielectric tensor of the metal , and via the Aharonov-Bohm (AB) effect (the frequencies of surface plasmon modes in nanotubes and mesoscopic rings change periodically by about 10 percent due to the AB effect ). Thus, magnetic field applied to a mesoscopic sample effectively changes its ”dielectric constant” as seen by the surface plasmons, and hence, changes the zero-point energy of the plasmon field. This fact may be interpreted as an additional vacuum contribution to the magnetic moment $`\mu _{vac}=E_0/H`$ of the metal sample. According to the results presented in , this vacuum contribution may be quite large, because of the rather large magnetic moments $`\mathrm{}\omega /H\mu _B`$ of the individual plasmon modes.

In this Letter I am going to calculate the plasmon zero-point energy of some metallic mesoscopic samples with simple geometries in an external magnetic field. The magnetic response of these plasmon vacuums appears to be rather large and diamagnetic. In thin films this predicted novel surface vacuum diamagnetism appears to be at least of the same order of magnitude as the magnetism of the bulk electrons. Thus, a novel type of magnetism shown by thin metallic samples, which is complementary to the well known Pauli paramagnetic and Landau diamagnetic contributions to the magnetism of electron gas will be demonstrated.

Let us consider the magnetic response of the surface plasmon vacuum of a square $`a\times a`$ region of a thin metal film with thickness $`d<<a`$ to an applied perpendicular magnetic field $`H`$ (Fig.1). Let us start by considering the dispersion law of a surface plasmon (SP), which propagates along the metal-dielectric interface in zero magnetic field. The SP field decays exponentially both inside the metal and the dielectric. Let us assume that both metal and dielectric completely fill the respective $`z<0`$ and $`z>0`$ half-spaces. In such a case the dispersion law can be written as

$$k^2=\frac{\omega ^2}{c^2}\frac{ϵ_dϵ_m(\omega )}{ϵ_d+ϵ_m(\omega )},$$
(1)
where we will assume that $`ϵ_m=1\omega _p^2/\omega ^2`$ according to the Drude model, and $`\omega _p`$ is the plasma frequency of the metal. This dispersion law is shown in Fig.2 for the cases of metal-vacuum and metal-dielectric interfaces. It starts as a ”light line” in the respective dielectric at low frequencies and approaches asymptotically $`\omega =\omega _p/(1+ϵ_d)^{1/2}`$ at very large wave vectors. The latter frequency corresponds to the so-called surface plasmon resonance. Under the surface plasmon resonance conditions both phase and group velocity of the SPs is zero, and the surface charge and the normal component of the electric field diverge. Since at every wavevector the SP dispersion law is located to the right of the ”light line”, the SPs of the plane metal-dielectric interface are decoupled from the free-space photons due to the momentum conservation law.

We are mostly interested in the region of the dispersion law near the plasmon resonance since our samples of interest are ”mesoscopic”, and hence $`a`$ is much smaller than $`2\pi c/\omega _{sp}`$, where $`\omega _{sp}`$ is the surface plasmon resonance frequency. The corrections to surface plasmon dispersion law in an applied magnetic field were calculated by Chiu and Quinn . Let us initially consider the most simple case of magnetic field $`H`$ applied perpendicular to the metal film bounded by vacuum. At large plasmon wavevectors $`k`$ the surface plasmon frequency

$$\omega _{sp}((\omega _p^2+\omega _c^2)/2)^{1/2}$$
(2)
does not depend on $`k`$ (Fig.2), where $`\omega _p`$ and $`\omega _c=eH/mc`$ are the plasma and cyclotron frequencies, respectively . As a result, all the surface plasmon modes have the same magnetic moment $`\mathrm{}\omega /H`$, and the total magnetic moment of the plasmon vacuum can be written as

$$\mu _{vac}\mu _B/2^{1/2}\mathrm{\Sigma }_k\omega _c/\omega _p,$$
(3)
where $`\mu _B`$ is the Bohr magneton, and summation has to be done over all the plasmon modes of the square region of the thin metal film under consideration. The surface plasmon eigenmodes of this square region are defined by the two-component wave vector $`(k_x,k_y)=\pi /a\times (n_x,n_y)`$, where $`n_x`$ and $`n_y`$ are integer. Due to Landau damping  the summation over all possible surface plasmon wave vectors has to be cut off at $`k_{max}k_F`$ (the electron Fermi momentum). Thus, the total number of plasmon modes on the top and bottom interfaces of the metal film is roughly $`8(k_Fa/\pi )^2`$, and

$$\mu _{vac}8a^2n_e^{2/3}(\omega _c/\omega _p)\mu _B$$
(4)
For thin metal films this surface vacuum diamagnetism is at least of the same order of magnitude as the contribution of the bulk electrons. Detailed description of various contributions to the magnetism of metallic samples can be found, for example in . The dominating mechanism of the magnetism of conductivity electrons is the paramagnetic contribution first obtained by Pauli  as

$$\chi _e^{pm}=\frac{3^{1/3}m\mu _B^2}{\pi ^{4/3}\mathrm{}^2}n_e^{1/3},$$
(5)
where $`m`$ is the electron mass. The diamagnetic contribution first obtained by Landau  is usually smaller by approximately a factor of 3 . Thus, we only need to compare the magnitudes of the Pauli paramagnetic contribution and the surface plasmon vacuum contribution to the electron magnetism of our thin metal film sample. Let us assume the free-electron model value for the plasma frequency of electron gas in the metal :

$$\omega _p^2=\frac{4\pi e^2n_e}{m},$$
(6)
and compare these two contributions. The Pauli paramagnetic moment of our sample can be written as follows:

$$\mu _{pm}a^2n_e^{2/3}\omega _c\frac{3^{1/3}md}{2\pi ^{4/3}\mathrm{}n_e^{1/3}}\mu _B=8a^2n_e^{2/3}(\omega _c/\omega _d)\mu _B,$$
(7)
where $`d`$ is the film thickness, and some characteristic frequency $`\omega _d=16\pi ^{4/3}\mathrm{}n_e^{1/3}/(3^{1/3}md)`$ is introduced. From this expression we immediately see that at small film thicknesses $`d`$, such that $`\omega _p\omega _d`$, Pauli paramagnetic and surface plasmon vacuum diamagnetic contributions have the same order of magnitude. The characteristic film thickness necessary for this situation to occur can be written as

$$d\frac{16\pi ^{4/3}\mathrm{}n_e^{1/3}}{3^{1/3}m\omega _p}\frac{8\pi ^{5/6}\mathrm{}}{3^{1/3}m^{1/2}en_e^{1/6}}2nm$$
(8)
Thus, according to this simple estimate Pauli paramagnetic and surface plasmon vacuum diamagnetic contributions to the magnetism of a thin film sample would be approximately equal at $`d2`$nm. Moreover, the contributions of Landau diamagnetism and plasmon vacuum diamagnetism should be about the same at $`d6`$nm. As a result, experimental measurements of the magnetic response of mesoscopic thin film samples performed as a function of film thickness in the 2-20 nm range have reasonable chance of detecting the surface plasmon vacuum diamagnetism. I should also point out that in reality the vacuum diamagnetic moment may be at least an order of magnitude larger than the value determined by equation (4), since the cut-off wave vector is known only by an order of magnitude. In addition, magnetic response measurements may be performed as a function of the dielectric constant of the substrate and/or absorbed layer on the interfaces of the metal film. According to Chiu and Quinn , in the presence of a dielectric layer on the metal surface and at large plasmon wavevectors $`k`$ the surface plasmon frequency in perpendicular magnetic field is given by

$$\omega _{sp}((ϵ_d^1\omega _p^2+\omega _c^2)/2)^{1/2},$$
(9)
where $`ϵ_d`$ is the dielectric constant of the layer. As a result, the magnetic moment of individual plasmon modes $`\mathrm{}\omega /H`$, and the total magnetic moment of the plasmon vacuum (defined by equations (3) and (4)) are multiplied by the refractive index of the dielectric $`n_d=ϵ_d^{1/2}`$. Since atomic monolayer quantities of the adsorbents are sufficient to shift the plasmon resonance , such measurements may be very useful in separating the relative contribution of the vacuum diamagnetism into the overall magnetic response of the sample (since there is no reason why a thin absorbed layer would alter either Pauli or Landau contributions). In addition, increase of the overall effect by another factor of 2 or 3 means that even thicker metal films (few tens of nanometers) may be used to make a mesoscopic sample of interest. This additional increase in the scale of the film thickness also means that calculations of plasmon eigenfrequencies based on the macroscopic Maxwell equations are much more reliable.

Except for the extreme sensitivity to the absorbed layers, the vacuum diamagnetism described above looks similar to the Pauli and Landau contributions with respect to temperature changes. Since surface plasmon eigenfrequencies involved are of the order of a few electron-Volts, no considerable changes in the surface plasmon magnetic response may be expected from absolute zero up to the room temperature. The best candidates for vacuum diamagnetism observations may be the mesoscopic samples made of gold, silver, copper or aluminum, since these metals exhibit very pronounced plasmon resonances . Even in the presence of absorbed layers the vacuum diamagnetic contribution defined by equation (4) is rather small due to the small $`(\omega _c/\omega _p)`$ factor. For a gold sample at $`H=1`$T this ratio is $`(\omega _c/\omega _p)10^4`$. Thus, a 1 micrometer by 1 micrometer square sample would have vacuum diamagnetic moment of $`10^4\mu _B`$ at $`H=1`$T , or a few tens of Bohr magnetons at $`H=10`$G. However, similar sensitivity has been achieved recently by Deblock $`et`$ $`al.`$  who reported measurements of the magnetic response of individual mesoscopic silver rings and obtained that the measured response of an individual silver ring is diamagnetic in the limit of zero magnetic field, which was not consistent with the available theoretical predictions. According to these measurements performed on an ensemble of $`1.5\times 10^5`$ 1 micrometer by 1 micrometer square silver rings with the thickness of 70 nm, the magnetic moment of an individual ring oscillates with the number of magnetic flux quanta with an amplitude of approximately $`30\mu _B`$ and the period of $`20`$G. Regardless of the nature of these oscillations, which we will briefly discuss below, the result of Deblock $`et`$ $`al.`$ for the unexpected diamagnetic response of an individual ring seams consistent with both the sign and the magnitude (according to Fig.4 from , $`\mu _{ring}10\mu _B`$) of the vacuum diamagnetic contribution described above.

Nevertheless, the results of  should be taken with some degree of caution, because of the chosen experimental geometry. Looking at the geometry of the experiment presented in the Fig.1 of , we immediately observe that the plasmon modes of individual rings in the sample are strongly coupled to each other. The surface plasmon wavelengths of the individual rings may be roughly estimated as $`\lambda _pL/n`$ where $`L=4\mu m`$ is the ring perimeter, and $`n`$ is an integer. This corresponds to the plasmon field decay outside a ring with an exponent $`\lambda _p/2\pi 640/n`$ nm. On the other hand, the distance between the individual rings in the periodic array shown in Fig.1 of  appears to be $`550nm`$. As a result, the plasmon modes of the individual rings are coupled to each other, and the calculation of surface plasmon spectrum of an array of such coupled cylindrical rings represents rather difficult problem. The fact that the rings are coupled indicates that Deblock $`et`$ $`al.`$ may not have achieved proper measurement of the magnetic response of an individual mesoscopic silver ring. The proper measurement can be done if the spacing of the individual rings will be made substantially larger than the ring perimeter, so that individual plasmon modes will be decoupled.

In some geometries the frequencies of individual plasmon modes may exhibit linear dependence on the applied magnetic field . For example, if magnetic field is applied parallel to the metal film (in the $`x`$-direction), than according to Chiu and Quinn  at large plasmon wavevectors $`k`$ the surface plasmon frequency is given by

$$\omega _{sp}((\omega _p^2+\frac{1}{2}\omega _c^2)/2)^{1/2}\pm \frac{1}{2}\omega _c,$$
(10)
where the two signs correspond to plasmon propagation in the $`y`$-direction. Such plasmon modes have magnetic moments $`\mathrm{}\omega /H\pm \mu _B`$. In a similar way, cylindrical surface plasmon modes of nanoholes and nanotubes  with nonzero angular momenta also have magnetic moments $`\pm \mu _B`$ . In rotationally symmetric samples vacuum magnetic moments of these “left “and “right “modes compensate each other, so that only the contributions quadratic in the magnetic field are left in the expression for the total zero-point energy of the plasmon vacuum, and we would come up with an expression more or less similar to equation (4) for the vacuum magnetic moment. However, this may not be the case for an asymmetric sample. While surface plasmon modes with very large wave vectors may not be affected much by the asymmetry of the sample shape, small number of low $`k`$ modes will be very sensitive to it. A small number of such “asymmetric “plasmon modes may be responsible for the diamagnetic moment $`10\mu _B`$ observed by Deblock $`et`$ $`al.`$ in zero magnetic field if the rings used in their experiment are not perfectly symmetric. In addition, asymmetry between the “left “and “right “plasmon modes may be caused by natural or magnetic field induced optical activity  of the mesoscopic sample itself and/or the absorbed layer on its surface. Optical activity of materials and samples is usually described by the difference in the refractive indices $`n_+`$ and $`n_{}`$ for the left and right circular polarizations (in other words, modes with opposite angular momenta) of light. Because of the different refractive index, an effective optical length of an optically active cylindrical sample is different for the left and right plasmon modes. As a result, the total number of the left and right plasmon modes may be slightly different, and the cylinder may posses an additional vacuum magnetic moment

$$\mu _{vac}(n_+n_{})han_e^{2/3}\mu _B,$$
(11)
where $`h`$ is the height of the cylinder and $`a`$ is its radius. Here we should also point out that all metals exhibit magnetic field induced optical activity , so that this effect may be also quite relevant for the mesoscopic metallic samples in an external magnetic field.

Another potentially important (although separate) question is how the Aharonov-Bohm effect may affect the vacuum diamagnetism of mesoscopic samples. Very recently Chaplik $`et`$ $`al.`$ demonstrated that the frequencies of surface plasmon modes in nanotubes and mesoscopic rings change periodically by about 10 percent due to the AB effect . The reason for this effect is the fact that both Fermi energy and the polarization operator of the nanotubes and nanorings changes periodically with magnetic flux due to the periodic dependence of the single-particle spectrum on magnetic flux. These periodic changes lead to the periodic changes of the surface plasmon eigenfrequencies, which can be found by solving the Poisson equation. Thus, magnetic moments of these modes $`\mathrm{}\omega /H`$ experience periodic oscillations too. According to the results of numerical calculations by Chaplik $`et`$ $`al.`$  the amplitude of these oscillations is of the order of $`\mu _B`$ for the zero-angular momentum ($`m=0`$) plasmon mode of a carbon nanotube at $`ka=1`$ (Fig.2(b) from ), while this mode does not have any angular momentum in zero magnetic field (for this mode $`d\omega /dH=0`$ at $`H=0`$). Similar to the classical result for cylindrical surface plasmons with nonzero angular momentum, plasmons of the nanotube with $`m0`$ have nonzero AB-induced magnetic moments (Fig.4 from ). While exact numerical calculations for the total AB-induced vacuum magnetic moment would be rather cumbersome, because of the complicated dependencies of the frequencies of the plasmon modes with arbitrary $`m`$ and $`k`$ on the magnetic flux , there is absolutely no reason for the total magnetic moment not to experience periodic oscillations with the amplitude of a few $`\mu _B`$ due to the AB effect. Thus, the periodic small oscillations of the magnetic response of the silver rings observed by Deblock $`et`$ $`al.`$  may also be at least partially attributed to the oscillations of vacuum magnetism.

In conclusion, zero-point energy of surface plasmon modes of mesoscopic metal samples in an external magnetic field has been considered. The magnetic response of plasmon vacuum has shown to be diamagnetic. In thin films this surface vacuum diamagnetism is at least of the same order of magnitude as the magnetism of the bulk electrons. Thus, a novel type of magnetism shown by thin metallic samples, which is complementary to the well known Pauli paramagnetic and Landau diamagnetic contributions to the magnetism of the electron gas has been demonstrated.

This work has been supported in part by the NSF grants ECS-0210438 and ECS-0304046.

Figure captions.

Fig.1 Model geometry for the calculations of the magnetic response of surface plasmon vacuum in the case of a square $`a\times a`$ region of a thin metal film with thickness $`d<<a`$ in an applied perpendicular magnetic field $`H`$.

Fig.2 Surface plasmon dispersion law for the cases of metal-vacuum and metal-dielectric interfaces in zero magnetic field (solid lines) and in the magnetic field applied perpendicular to the film surface (dashed lines). At large plasmon wavevectors $`k`$ the changes in surface plasmon frequency in the applied magnetic field do not depend on $`k`$.


format: text

text:
# Grain Growth in the Dark Cloud L1251

## 1 Introduction

The ratio of total to selective extinction, $`R_V`$=$`A_V`$/$`E_{BV}`$, is a measure of wavelength dependence of the interstellar extinction, representing the shape of the extinction curve. $`R_V`$ is an important parameter to characterize dust properties at optical wavelengths because it is closely related to the size and composition of dust grains (e.g., Mathis, Rumpl, & Nordsieck 1977; Hong & Greenberg 1978; Kim, Martin, & Hendry 1994). It is well established that $`R_V`$ has a uniform value of $``$3.1 in the diffuse interstellar medium (e.g., Whittet 1992), but it tends to be much higher in dense molecular clouds ($`R_V`$=4$``$6; e.g., Vrba, Coyne, & Tapia 1993; Strafella et al. 2001). Many theoretical predictions attribute the higher $`R_V`$ in the dense cloud interior to a significant change of dust properties, particularly in size (e.g., Mathis 1990; Ossenkopf & Henning 1994). Since the timescale of grain growth due to coagulation of smaller grains and/or accretion of molecules in the gas phase onto the larger grain surface depends greatly on density (e.g., Tielens 1989), the higher $`R_V`$-values observed in dense molecular clouds are most likely to represent the growth of dust grains.

There have been a number of observational reports on $`R_V`$ measurements to date (e.g., Larson et al. 2000; Whittet et al. 2001). However, they are mostly based on the spectrophotometry of a limited number of stars lying in the background of the clouds, and the $`R_V`$-values have been measured only toward the selected stars. A large-scale distribution of $`R_V`$ over an entire cloud, as well as its local variation inside, are therefore still poorly known.

In order to investigate how much $`R_V`$ can vary in a single cloud, we have carried out imaging observations of L1251 at multiple wavelengths, $`B`$, $`V`$, $`R`$, and $`I`$, using the 105 cm Schmidt telescope at the Kiso Observatory, Japan. This cloud is located in the Cepheus Flare ($`\alpha `$$``$$`22^h35^m`$, $`\delta `$$``$$`75^{}00^{}`$; J2000.0) and has an elongated shape. Since L1251 is an isolated dark cloud located at a relatively high Galactic latitude ($``$$`15^{}`$), it is expected that the cloud is free from contamination due to unrelated clouds in the same line of sight. Kun & Prusti (1993) estimated its distance to be 300$`\pm `$50 pc with a color excess measurement and spectral classification of stars. They also found 12 Hα emission-line stars at the eastern part of the cloud, five of which are detected as $`IRAS`$ point sources. Global molecular distribution of the cloud was first revealed by Sato & Fukui (1989) and Sato et al. (1994) through 13CO and C18O ($`J`$=1$``$0) observations carried out with the 4 m radio telescope at the Nagoya University (HPBW=2.7). Their observations clearly revealed the cometary “head-tail” morphology of L1251. They also found that two of the $`IRAS`$ point sources located in the eastern “head” part of the cloud are associated with molecular outflows, while no sign of ongoing star formation was found in the western “tail.” These characteristics of L1251 make this cloud a suitable site to investigate the possible difference in the dust properties (e.g., $`R_V`$) between star-forming and non-star-forming regions in a single cloud.

In this paper, we report results of the $`BVRI`$ observations toward L1251. The observational procedures are described in $`\mathrm{\S }`$ 2. The distribution of extinction ($`A_B`$ and $`A_V`$) and color excess at each band ($`E_{BV}`$, $`E_{VR}`$, $`E_{VI}`$) are derived from the obtained data, which reveal the global dust distribution of the cloud. These are shown in $`\mathrm{\S }\mathrm{\S }`$ 3 and 4. We also derive $`R_V`$ distribution over the cloud on the basis of color excess measurements and an empirical relation between $`R_V`$ and $`A_\lambda `$/$`A_V`$ (Cardelli, Clayton, & Mathis 1989). A significant difference has been found in $`R_V`$ distribution between the dense head and the relatively diffuse tail of the cloud, suggesting a larger dust size in the dense head than in the tail. We introduce our procedure to create $`R_V`$ map and compare the obtained $`R_V`$ distribution with the other data in $`\mathrm{\S }`$ 5. Our conclusions are summarized in $`\mathrm{\S }`$ 6.

## 2 Observations and Data Reduction

We carried out imaging observations of L1251 from 2000 November 27 to December 2 using the 105 cm Schmidt telescope equipped with the 2kCCD camera at the Kiso Observatory. The camera has a backside-illuminated 2048 $`\times `$ 2048 pixel chip covering a wide field of view of $``$$`50^{}`$$`\times `$$`50^{}`$ with a scale of 1.′′5 pixel-1 at the focal plane of the telescope. Using the camera, we mapped L1251 with the $`B`$, $`V`$, $`R`$, and $`I`$ filters at a $`40^{}`$ grid spacing along the equatorial coordinates, so that adjacent frames could have an overlap of $``$10 on the sky. We performed integrations of 200, 100, 100, and 50 s for the images obtained with the $`B`$, $`V`$, $`R`$, and $`I`$ filters, respectively. The typical seeing was 3.′′0$``$4.′′5 during the observations. In total, we obtained 10 frames per filter, covering the entire extent of the cloud.

To calibrate the images, we obtained dome flat frames in the beginning of the observations and frequently obtained bias frames during the observations. In addition, we observed some sets of standard stars listed by Landolt (1992) to correct the atmospheric extinction and calibration.

We processed the calibrated images in a standard way using the IRAF packages.111IRAF is distributed by the National Astronomy Observatory, which is operated by the Association of Universities for Research in Astronomy, Inc., under coorperative agreement with the National Science Foundation. After replacing values at the defected pixels in the detector (i.e., dead pixels) by those interpolated using surrounding pixels, we applied the stellar photometry package DAOPHOT (Stetson 1987) to the images to detect stars and then to make photometric measurements. To improve the photometry we further measured the instrumental magnitudes of the detected stars by fitting a point-spread function (PSF) and then calibrated them for the atmospheric extinction. Among all of the detected stars we selected the ones having an intensity greater than $`4.0`$ $`\sigma `$ noise level of the sky background. Some false detection, however, could not be completely ruled out, which we later removed by eye inspection.

For calculating the plate solution, we compared the pixel coordinates of a number of the detected stars with the celestial coordinates of their counterparts in the Digitized Sky Survey (DSS), which resulted in a positional error of less than 1′′ (rms). All of the stars detected in our observations were summarized separately with respect to the filters we used ($`BVRI`$), which we will adopt to generate an extinction map using a star count method. The limiting magnitudes are 20.0, 19.7, 19.0, and 18.3 mag for the $`B`$, $`V`$, $`R`$, and $`I`$ bands, respectively. Among all of the identified stars $``$33,000 were detected in all of the four bands. These stars are useful for probing the distribution of the color excess in the cloud. We converted the instrumental magnitudes of these stars into the magnitudes in the standard Johnson-Cousins photometric system (Johnson & Morgan 1953; Cousins 1976), and we also measured their colors ($`B`$$``$$`V`$, $`V`$$``$$`R`$, and $`V`$$``$$`I`$). For the color excess measurements, we constructed a $`BVRI`$-matched star list for $``$30,600 stars whose photometric errors in the color measurements, $`\sigma _{V\lambda }`$, are less than 0.15 mag. One should note that we excluded all of the known young stellar objects (Sato & Fukui 1989; Sato et al.1994; Kun & Prusti 1993) cross-identified in our photometries in order to avoid errors in the following analyses since they are unrelated to the global extinction in L1251.

## 3 Distribution of Extinction

On the basis of the stars detected in the $`B`$ and $`V`$ band, we derived extinction (i.e., $`A_B`$ and $`A_V`$) maps of L1251 using a star count method. Given a cumulative stellar density in the extinction-free reference field $`N_{\mathrm{ref}}`$ as a function of the $`\lambda `$-band magnitude $`m_\lambda `$ (i.e., the Wolf diagram; Wolf 1923), a logarithmic cumulative stellar density log $`N`$ measured at the threshold magnitude $`m_{\lambda ,0}`$ can be converted into the extinction $`A_\lambda `$ as

$`A_\lambda (l,b,m_{\lambda ,0})=m_{\lambda ,0}f^1\left[\mathrm{log}N(l,b,m_{\lambda ,0})\right],`$ (1)
where $`l`$ and $`b`$ are the Galactic coordinates, and $`f^1`$ is the inverse function of $`f`$ defined as log $`N_{\mathrm{ref}}`$ = $`f(m_\lambda )`$. The function $`f`$ is often assumed to be linear as $`f=a+bm_\lambda `$, where $`b=\left[d\mathrm{log}N_{\mathrm{ref}}/dm_\lambda \right]_{m_\lambda =m_{\lambda ,0}}`$ (e.g., Dickman 1978).

To derive $`N(l,b,m_{\lambda ,0})`$ in equation (1) from the obtained data, we first set circular cells with a 4 diameter spaced by 2 along the Galactic coordinates, and we counted the number of stars in the cells, setting the threshold magnitude to be $`m_{B,0}`$ = 19.5 mag and $`m_{V,0}`$ = 19.2 mag, which is well above the detection limit in our observations. We then smoothed the stellar density map to a 6 resolution with Gaussian filtering for a better signal-to-noise ratio to produce the final stellar density map $`N(l,b,m_{\lambda ,0})`$, which we adopted to derive the $`A_B`$ and $`A_V`$ distribution. After smoothing typical stellar density in the reference field and in the most opaque region was $``$20 cell-1 and $``$1 cell-1 for B band and $``$30 cell-1 and $``$2 cell-1 for V band, respectively.

We composed the stellar luminosity function $`f_{\mathrm{ref}}(m_\lambda )`$ in the reference field, i.e., the area expected to be free from extinction in the observed region (see Fig. 1). We found that the function is actually not linear, but has a curved shape. We also found that stellar densities at the threshold magnitude $`m_{\lambda ,0}`$ systematically decrease along with the Galactic latitudes $`b`$. To deal with this problem, we fitted the average luminosity function $`f_{\mathrm{ref}}(m_\lambda )`$ by a 4th-order polynomial and scaled it to the background stellar density $`N_{\mathrm{back}}`$ at $`m_\lambda `$$`=`$$`m_{\lambda ,0}`$ as

$`f(l,b,m_\lambda )={\displaystyle \frac{\mathrm{log}N_{\mathrm{back}}(l,b,m_{\lambda ,0})}{f_{\mathrm{ref}}\left(m_{\lambda ,0}\right)}}{\displaystyle \underset{n=1}{\overset{4}{}}}c_nm_\lambda ^n,`$ (2)
where $`c_n`$ represents the best-fitting coefficients. The value $`f_{\mathrm{ref}}(m_{\lambda ,0})`$ corresponds to the mean stellar density at $`m_{\lambda ,0}`$ in the reference field. We determined $`N_{\mathrm{back}}(l,b,m_{\lambda ,0})`$ by assuming an exponential dependence on $`b`$ as

$`N_{\mathrm{back}}(l,b,m_{\lambda ,0})=N_0\mathrm{exp}\left(\alpha \left|b\right|\right),`$ (3)
where $`N_0`$ and $`\alpha `$ are constants determined using values of $`N(l,b,m_{\lambda ,0})`$ in the reference field. We obtained the coefficients ($`N_0`$, $`\alpha `$) to be (61.5, 0.07) for $`B`$ band and (106.2, 0.09) for $`V`$ band.

We numerically calculated $`f^1`$, the inverse function of $`f`$ in equation (2), and used it in equation (1) to derive $`A_B`$ and $`A_V`$ after substituting the actually measured stellar density $`N(l,b,m_{\lambda ,0})`$ into $`f^1`$. The resulting $`A_V`$ map is shown in Figure 1. The distribution of $`A_B`$ (not shown) is very similar to that of $`A_V`$. Total uncertainty in $`A_\lambda `$ varies from region to region, depending on the number of stars falling in the cells. Typical uncertainties in the reference field ($`A_\lambda `$ $``$ 0 mag) and in the most opaque region are estimated to be $``$0.5 and $``$0.8 mag for the $`B`$ band, and $``$0.35 and $``$0.7 mag for the $`V`$ band, respectively.

The $`A_V`$ map reveals an interesting dust distribution in L1251 with a cometary head-tail morphology, which is similar to the molecular distribution traced in 13CO (Sato et al. 1994).

## 4 Distribution of Color Excess

We derived the color excess distribution in L1251 using the near-infrared color excess (NICE) method originally developed by Lada et al. (1994). The method compares the color index of stars in the reference field with those in the dark clouds. If we assume that the population of the stars across the observed field is invariable, the mean stellar color in the reference field can be used to approximate the mean intrinsic color of stars (i.e., $`<`$$`V`$$``$$`\lambda `$$`>`$int$`=`$$`<`$$`V`$$``$$`\lambda `$$`>`$ref), and the color excess can be derived by subtracting $`<`$$`V`$$``$$`\lambda `$$`>`$ref from the observed color of stars. To obtain the color excess maps, we adopted an “adaptive grid” technique (Cambrésy 1998, 1999) to the $``$30,600 stars detected in all of the four bands ($`B`$, $`V`$, $`R`$, and $`I`$). This technique consists in fixing the number of counted stars $`N`$ per variable cell to achieve a uniform noise level in the resulting map. To estimate the mean color excess $`<`$$`E`$V-λ$`>`$ ($`\lambda `$$`=`$$`B`$, $`R`$, or $`I`$) in circular cells spaced by 2 along the Galactic coordinates with undefined sizes, we calculated the mean color of the closest $`N`$ stars to the center positions of the cells by setting $`N`$ to be 10 and derived their shifts from the mean color in the reference field $`<`$$`V`$$``$$`\lambda `$$`>`$ref as

$`E_{V\lambda }(l,b)=\left[{\displaystyle \underset{i=1}{\overset{N}{}}}{\displaystyle \frac{\left(V\lambda \right)_i}{N}}\right](l,b)V\lambda _{\mathrm{ref}},`$ (4)
where ($`V`$$``$$`\lambda `$)i is the color index of the $`i`$th star in a cell. The distance to the $`N`$th star from the cell center gives the radius of the cell, corresponding to one half of the angular resolution. $`<`$$`V`$$``$$`\lambda `$$`>`$ref is the mean color of the $``$2000 cells in the reference field, and we determined $`<`$$`B`$$``$$`V`$$`>`$ref, $`<`$$`V`$$``$$`R`$$`>`$ref, and $`<`$$`V`$$``$$`I`$$`>`$ref to be 1.02, 0.55, and 1.23 mag, respectively.

In the reference field, we found that $`<`$$`V`$$``$$`\lambda `$$`>`$ slightly but systematically decreases with increasing Galactic latitudes $`b`$. This may be caused by the presence of diffuse dust layer of the Galaxy and/or gradual change in the stellar population toward the Galactic plane. We calibrated the effect by substituting $`<`$$`V`$$``$$`\lambda `$$`>`$ref in equation (4) for the following equation assuming an exponential dependence on $`b`$ as

$`V\lambda (l,b)=A\mathrm{exp}\left(\alpha \left|b\right|\right),`$ (5)
where $`A`$ and $`\alpha `$ are the fitting coefficients. We determined the coefficients ($`A`$, $`\alpha `$) to be (1.60, 0.03), (2.60, 0.05), and (0.84, 0.03) for $`<`$$`B`$$``$$`V`$$`>`$ref, $`<`$$`V`$$``$R$`>`$ref, and $`<`$$`V`$$``$$`I`$$`>`$ref, respectively.

We then smoothed the color excess maps with a Gaussian filter (FWHM$`=`$4) to reduce noise. The resulting resolution of the maps is typically $``$5 in the reference field and $``$10 in the most opaque region. An example of the obtained color excess maps is shown in Figure 2, where we show the spatial distribution of $`E_{BV}`$. The distributions of $`E_{VR}`$ and $`E_{VI}`$ (not shown) are similar to that of $`E_{BV}`$. L1251 is also traced in the all-sky $`E_{BV}`$ map with 6 resolution (Schlegel, Finkbeiner, & Davis 1998), which is derived from far-infrared emission data ($`COBE`$/$`DIRBE`$, $`IRAS`$/$`ISSA`$). The cloud appears similar in both the all-sky map and our color excess maps, indicating that the optical color excess and the far-infrared thermal emission of dust trace the same region.

Uncertainty in deriving $`<`$$`E_{V\lambda }`$$`>`$ is given by following equation (e.g., Arce & Goodman 1999):

$`\sigma _{E\left(V\lambda \right)}=\left[\sigma _{\mathrm{cell}}^2+{\displaystyle \underset{i=1}{\overset{N}{}}}\left({\displaystyle \frac{\sigma _{\left(V\lambda \right)_i}}{N}}\right)^2\right]^{1/2}`$ (6)
where $`\sigma _{(V\lambda )_i}`$ is the photometric error in $`V`$$``$$`\lambda `$ of the $`i`$th star in each cell. The value $`\sigma _{\mathrm{cell}}`$ is the standard deviation of the $`N`$-sampling distribution of the mean, where the original distribution is the $`V`$$``$$`\lambda `$ color of the stars in the reference field. It is equivalent to the standard deviation of the $`<`$$`E_{V\lambda }`$$`>`$ distribution in the reference field. We estimated $`\sigma _{\mathrm{cell}}`$ to be 0.04, 0.03, and 0.05 mag for the color of $`B`$$``$$`V`$, $`V`$$``$$`R`$, and $`V`$$``$$`I`$, respectively. In the $`E_{BV}`$, $`E_{VR}`$, and $`E_{VI}`$ maps, a typical $`\sigma _{<E_{V\lambda }>}`$-value calculated using equation (6) is 0.05, 0.03, and 0.05 mag, respectively.

We derived the color excess diagram $`E_{BV}`$ versus $`E_{VI}`$ to see the global trend of reddening by dusts in L1251, which is shown in Figure 3. We investigated the relation separately in the star-forming head ($`l`$$`>`$114.2) and in the non-star-forming tail ($`l`$$``$114.2). For comparison, we also show $`E_{BV}`$ versus $`E_{VI}`$ relations in the cases of $`R_V`$$`=`$3.1 and 6.5, adopting an empirical extinction law reported by Cardelli et al. (1989).

It is clear that there is a significant difference in the $`E_{BV}`$ versus $`E_{VI}`$ relations obtained in the two parts of the cloud. In the non-star-forming tail the data are distributed along the line of $`R_V`$$`=`$3.1, a typical value in the diffuse ISM, and the linear least-squares fit to the data actually yields a slope of $`E_{VI}`$$`/`$$`E_{BV}`$$`=`$1.21, which corresponds to $`R_V`$$`=`$3.2. On the other hand, the data in the star-forming head are distributed above the line of $`R_V`$$`=`$3.1 in the range of $`E_{BV}`$$``$0.4 mag, and they gradually get close to the line of $`R_V`$$`=`$6.5 for larger $`E_{BV}`$. It is evident that optical properties of dust grains significantly vary even in a single cloud. The different extinction laws are well separated into two groups only when we divide them into star-forming and non-star-forming parts. Since $`R_V`$ increases with increasing size of dust grains, larger $`R_V`$-value in the star-forming head is indicative of the evolution (i.e., growth in size) of dust grains. We note that $`E_{VI}`$ increases in the same way (along with $`R_V`$$``$3.1) up to $`E_{BV}`$$``$0.4 mag both in the tail and in the head part, which indicates dust optical properties remain normal at the diffuse boundary of the cloud.

## 5 Derivation of $`R_V`$ Distribution and Evidence for Grain Growth

In order to investigate the local variation of the dust properties in L1251, we derived the $`R_V`$ distribution using a relationship between $`R_V`$ and $`E_{V\lambda }`$/$`E_{BV}`$ that we produced from the $`R_V`$-dependent extinction law derived by Cardelli et al. (1989). They found that observed $`A_\lambda `$/$`A_V`$ versus $`R_V^1`$ relations toward a number of line-of-sight stars are well fitted over the wavelength range from 0.125 $`\mu `$m to 3.5 $`\mu `$m as

$`A_\lambda /A_V=a(x)+b(x)R_V^1,`$ (7)
where $`x`$ is $`\lambda ^1`$ $`\mu `$m-1. The wavelength-dependent coefficients $`a(x)`$ and $`b(x)`$ are shown in equations (2a-2b) and (3a-3b) of Cardelli et al. (1989). By using equation (7) we reproduced $`R_V`$ and $`E_{V\lambda }`$/$`E_{BV}`$ relations with the color excess ratios as $`E_{V\lambda }`$/$`E_{BV}`$=($`A_V`$-$`A_\lambda `$)$`/`$($`A_B`$-$`A_V`$).

Color excess ratios of $`E_{VR}`$/$`E_{BV}`$ and $`E_{VI}`$/$`E_{BV}`$ are used to obtain the $`R_V`$-values at the cells that we used to derive the color excess maps. In each cell we numerically calculated the $`R_V`$-value, which minimizes the sum of residual variance; i.e., $`[`$($`E_{VR}`$/$`E_{BV}`$)obs$``$($`E_{VR}`$/
$`E_{BV}`$)model$`\}`$2 $`+`$ $`\{`$($`E_{VI}`$/$`E_{BV}`$)obs$``$($`E_{VI}`$/$`E_{BV}`$)model$`]`$2. We carried out the $`R_V`$ fitting only in the region with $`E_{BV}`$$``$0.4 mag in order to avoid a large error in $`R_V`$ arising from too small $`E_{BV}`$ values. The resulting $`R_V`$ map is shown in Figure 4. As seen in the figure, $`R_V`$ in the star-forming head is systematically higher than in the tail with the maximum value of $`R_V`$$``$6. The mean value of $`R_V`$ in the head and tail are $``$4.6 and $``$3.2, respectively. The feature in the $`R_V`$ distribution is consistent with what we found in the color excess diagrams derived in the previous section (Fig. 3).

Uncertainties in the $`R_V`$ measurement, $`\sigma `$$`_{R_V}`$, are estimated using Monte Carlo simulations. We first generated noise that follows a Gaussian distribution with the standard deviation equal to the 1 $`\sigma `$ error of $`E_{V\lambda }`$$`/`$$`E_{BV}`$ measurements, and then we performed the $`R_V`$ fitting to the noise-added $`E_{V\lambda }`$$`/`$$`E_{BV}`$ data. We repeated this calculation 1000 times in each cell, and regarded the standard deviation of the resulting $`R_V`$-values as $`\sigma `$$`_{R_V}`$. The value of $`\sigma `$$`_{R_V}`$ varies from $``$0.5 to $``$0.9 (0.69 on the average) and tends to be relatively higher in cells with small $`E_{BV}`$. In addition to the random noise our $`R_V`$ map may suffer from a systematic error due to a number of unidentified pre-main-sequence (PMS) stars formed in the dense head. A certain fraction of PMS stars are known to have a blue color excess (e.g., Kenyon & Hartmann 1995), which may result in an overestimate in $`R_V`$. However, this possible error should be negligible in our case, as we discuss in the Appendix.

We further checked the validity of the color-excess-based $`R_V`$ measurements by using $`A_B`$ and $`A_V`$ derived with star count method. Though large uncertainties (typically $``$0.5 mag in our data) are included in $`A_B`$ and $`A_V`$ data, $`R_V`$ can be directly computed as $`R_V`$$`=`$$`A_V`$$`/`$$`(A_BA_V)`$ without using an empirical model. We made $`A_V`$ versus $`A_B`$ diagram for the same cells used in $`R_V`$ measurements and separated them into the star-forming head ($`l`$$`>`$114.2) and non-star-forming tail ($`l`$$``$114.2). The slope of the $`A_V`$ versus $`A_B`$ relations yields mean $`R_V`$-value of data points as $`<`$$`R_V`$$`>`$$`=`$1$`/`$($`\alpha `$$``$1), where $`\alpha `$$`=`$$`A_B/A_V`$. The linear least-squares fit to each dataset resulted in $`\alpha `$head$`=`$1.23$`\pm `$0.09 and $`\alpha `$tail$`=`$1.34$`\pm `$0.07, which correspond to $`<`$$`R_V`$$`>`$head$``$4.3 and $`<`$$`R_V`$$`>`$tail$``$3.0, respectively. These values are consistent with the results of color-excess-based $`R_V`$ measurements.

In Figure 4 it is particularly noteworthy that there is a good coincidence between the local peaks of the $`R_V`$ distribution and the locations of YSOs. There are three strong $`R_V`$ peaks in the head at ($`l`$,$`b`$)$`=`$(114.60, 14.47; $`R_V`$=6.1), (114.50, 14.67; $`R_V`$=5.8), and (114.27, 14.77; $`R_V`$=5.6). The first two peaks show a good agreement with the locations of the two outflow sources reported by Sato & Fukui (1989) and Sato et al. (1994), though the third peak has no apparent corresponding source. In addition, the three $`R_V`$ peaks also coincide with the positions of molecular cloud cores found in C18O (named “B”, “C”, and “E” by Sato et al. 1994; see their Table 1 and Fig. 2), indicating that the region with high $`R_V`$-values is rich in molecular gas.

To investigate the dependence of $`R_V`$ on the dust column density, we derived a diagram of $`A_V`$ versus $`R_V`$. For a direct comparison we remeasured the $`A_V`$-values with the same cells used to produce the $`R_V`$ map. In Figure 5 we show the diagram measured both in the dense head and in the diffuse tail. It is clear that the $`R_V`$ increases along with $`A_V`$ in the head, while that in the tail remains constant regardless of $`A_V`$.

All of the above results obtained from the comparison of the $`R_V`$ map with the other data, i.e., the coincidence of the $`R_V`$ peaks and the outflow sources, as well as the C18O cores, and the difference in the $`A_V`$ versus $`R_V`$ relation between the head and the tail, strongly indicate that $`R_V`$ is higher in the dense star-forming head than those in the diffuse tail without any signs of star formation. Because a higher value of $`R_V`$ is expected for dusts with a larger size, the dusts in the head of L1251 may actually be larger than in the diffuse tail, indicating the growth of dust grains through coagulation and accretion of smaller dust particles as well as molecules in gas phase onto larger dust grains in the dense molecular cloud cores experiencing star formation. The $`R_V`$ distribution derived in this study may not be sufficient to probe the densest region of L1251 ($`A_V`$$`>`$4$``$5 mag) in detail because of a large optical depth at the observed wavelengths. However, a similar change of dust properties most likely representing a large dust size in molecular clouds was recently reported through submillimeter observations of the dust emission (e.g., MCLD123.5$`+`$24.9: Bernard et al. 1999; Taurus: Stepnik et al. 2003). The idea of grain growth is also supported by their study based on the optically much thinner submillimeter observation, implying that dust grains with a larger size than in the diffuse ISM may be common in dense molecular clouds.

Our procedure to derive the $`R_V`$ distribution on the basis of the average color excess measurements is useful to investigate a large-scale variation of the dust optical properties. It may be applicable also to the near infrared data in the future, with which we should be able to measure the $`R_V`$ distribution in the densest region of L1251 to confirm the significant change of $`R_V`$ therein.

## 6 Summary

We have carried out optical observations toward L1251 at multiple bands, $`B`$, $`V`$, $`R`$, and $`I`$, using the 105 cm Schmidt telescope equipped with the 2kCCD camera at the Kiso Observatory in Japan. Extinction and color excess maps of the cloud were derived using a star count method and a NICE method. From the observed $`E_{V\lambda }`$$`/`$$`E_{BV}`$ ($`\lambda `$$`=`$$`R`$,$`I`$), we also derived the $`R_V`$ distribution using an empirical $`R_V`$ versus $`A_\lambda `$$`/`$$`A_V`$ relation (Cardelli et al. 1989).

(1) We investigated the color excess distribution of the cloud on the basis of the $``$30,600 stars detected in our observations. Color excess diagram, $`E_{BV}`$ versus $`E_{VI}`$, indicates a significant difference in $`R_V`$ between the star-forming dense head ($`l`$$`>`$114.2) and the non-star-forming diffuse tail ($`l`$$``$114.2).

(2) We derived the $`R_V`$ distribution over the cloud on the basis of the color excess measurements and an empirical relation between $`R_V`$ and $`A_\lambda `$$`/`$$`A_V`$ (Cardelli et al. 1989). From the resulting $`R_V`$ map we found that $`R_V`$ is systematically higher in the dense head ($`R_V`$$`=`$4$``$6) than in the diffuse tail ($`R_V`$$``$3.2), which is consistent with what we found in the $`E_{BV}`$ versus $`E_{VI}`$ diagram. We confirmed the validity of color-excess-based $`R_V`$ measurements by comparing them with independently measured $`R_V`$-values; i.e., $`A_V`$$`/`$($`A_B`$$``$$`A_V`$), which were directly computed using $`A_B`$ and $`A_V`$ data from star count analysis.

Local peaks of the $`R_V`$ distribution in the head of L1251 coincide with the positions of the young stars accompanied by molecular outflows and the molecular cloud cores identified in C18O. We also found that $`R_V`$ increases with increasing $`A_V`$ in the dense star-forming head, while $`R_V`$ is likely to be independent of $`A_V`$ in the diffuse tail. Since $`R_V`$ is closely related to the size of dusts, we suggest that these results are indicative of the growth of dust grains in the dense head due to coagulation and accretion of smaller dust particles and molecules onto the larger dust grains.

We are very grateful to M. Tamura, K. Tatematsu, and T. Mizuno for their helpful comments and suggestions. Thanks are due to the staff at Kiso Observatory including K. Nakada, T. Miyata, T. Soyano, K. Tarusawa, Y. Tanaka, H. Mito, T. Aoki, N. Ito, and S. Nishiura for their kind support during the observations. This work was financially supported by the Grant-in-Aid for Scientific Research by the Japanese Ministry of Education, Science, Sports and Culture (Nos. 13640233 and 14022214).

## Appendix A A Possible Systematic Uncertainty in $`R_V`$

Here we consider the degree of uncertainties caused by the presence of unknown faint PMS stars in the star-forming head of L1251, which escaped detection in the previous Hα survey (Kun & Prusti 1993). It is well known that a certain fraction of PMS stars show the optical blue excess, which is prominent in $`U`$$``$$`B`$ (e.g., Kenyon & Hartmann 1995) and possibly in $`B`$$``$$`V`$. Their blue intrinsic color should result in an underestimate of $`E_{BV}`$, for which we may overestimate $`R_V`$.

In order to examine the fraction of blue excess PMS stars at $`B`$$``$$`V`$, we made a $`B`$$``$$`V`$ versus $`V`$$``$$`I`$ diagram of 96 PMS stars in the Taurus-Auriga region (see Table A1 of Kenyon & Hartmann 1995). We found 23 ($``$25$`\%`$) sources show a blue excess of greater than 0.2 mag compared with the main-sequence stars, while the rest of the sources have a similar color to the main-sequence stars with some reddening.

If we assume that the ratio of the blue excess PMS stars measured in the Taurus-Auriga region is also valid in L1251, only two to three stars (25$`\%`$ of $`N`$$`=`$10) in each counting cell are expected to be blue excess PMS stars even in an extreme case that all of the stars detected in our observations are PMS stars. However, at the maximum estimate of the blue excess PMS stars (about two to three stars per cell), there might be an influence on the value of $`E_{BV}`$ calculated as the average color of 10 stars (eq. ). We therefore recalculated $`E_{BV}`$ taking the median color of 10 stars instead of the average, because it should be much less affected by a few blue excess PMS stars. We compared the two values of $`E_{BV}`$, the median and the average colors of 10 stars, and found that they are quite consistent with each other even in the densest head of L1251 with a typical difference of only $``$0.02 mag. This means that the possible error in $`R_V`$ due to the blue excess PMS stars is negligible, and the high $`R_V`$-value derived in the star-forming dense head is most likely to be real.


text:
# The Impact of Pollution on Stellar Evolution Models

## 1 Introduction

The discoveries of dozens of giant planets in orbits with semi-major axes less than 1 AU (for a complete listing see www.exoplanets.org), which began in 1995 (Mayor & Queloz, 1995), posed a problem for planet formation theories. Such theories indicate that giant planets ought to form in orbits with semi-major axes greater than 1 AU. However, through interactions with other planets and the protoplanetary disk it is possible for these planets to migrate inward until they reach a stable orbit (Lin & Shigeru, 1997). The inward migration can sweep material in the disk onto the surface of the star. Gonzalez (1997) pointed out that if the accreted material is metal-rich then stars with planets (SWPs) will exhibit enhanced metallicity, a concept referred to as “stellar pollution”.

The case has been made in recent years that SWPs are metal rich, for example Laws & Gonzalez (2002). The two main explanations for these stars to be metal-rich are a high intrinsic metallicity and the accretion of $``$5 $`\mathrm{M}_{}`$ of iron from the protoplanetary disk, for a more complete discussion see Murray & Chaboyer (2002).

Based upon the apparent deficiency of heavy elements in the asteroid belt, and considerations of orbital dynamics in the early solar system, Murray et al. (2001) found that the formation of Jupiter likely led to the accretion of $``$0.5 $`\mathrm{M}_{}`$ of iron onto the Sun. Their analysis of the correlations between \[Fe/H\]–age–mass in local field stars lead Murray et al. (2001) to suggest that stellar pollution effects in the solar neighborhood appear to be quite common at a similar level. Using a Monte Carlo approach they estimate that field stars have accreted an average of $``$0.5 $`\mathrm{M}_{}`$ of iron with standard deviation equal to half the mean. Hence, Murray et al. (2001) suggest that most protoplanetary disks are similar to that of the Sun and not the stars in the SWPs sample.

In a paper on modeling the evolution of SWPs, Ford et al. (1999) discuss the need for stellar evolution models which have enhanced surface metallicities relative to the interiors. The authors briefly explain their approach to the problem which involves increasing the metallicity of the surface convection zone and/or decreasing the interior metallicity, and discuss the impact of such models on the Sun and a few SWPs. The main thrust of the paper is, however, modeling SWPs with standard stellar evolution models.

In this paper, stellar evolution models are presented which allow for the addition of metal-rich material to the surface. The behavior of these models is compared standard stellar evolution models. The modifications made to standard stellar evolution models are explained in $`\mathrm{\S }`$2. The basic impact of these changes on the evolution of a star is presented $`\mathrm{\S }`$3. In $`\mathrm{\S }`$4 the metal-rich parent star $`\tau `$ Bootis is modeled with both polluted and unpolluted models and the results are compared. In $`\mathrm{\S }`$5 polluted model isochrones are compared to observational data on the Hyades; the likelihood that the Hyades stars have experienced low-level pollution is discussed. The attempt to create polluted stellar evolution models is summarized in $`\mathrm{\S }`$6.

## 2 Stellar Models

Stellar evolution models were constructed using the Dartmouth stellar evolution code (Chaboyer et al., 2001). Helium and heavy element diffusion are included using Thoul et al. (1994). All models use a solar calibrated mixing length $`\alpha `$ = 1.75. The models were evolved from pre main sequence, fully convective polytropes to the main sequence turnoff in $``$800 time steps and include overshooting of the surface convection zone.

Murray et al. (2001) suggest that we may model the depletion of lithium from the photosphere by introducing a convective overshoot to the surface convection zone. They present a model for the mass of the surface mixed layer (SML) as a function of stellar mass and \[Fe/H\] which has been implemented in the code for this work. Each time the code calls for the surface convection zone to be mixed we calculate the mass of the SML and mix down to its base. For a plot of SML mass as a function of stellar mass, see Figure 2 of Murray et al. (2001).

In the literature on stellar pollution the amount of accreted material is often referred to in Earth masses ($`\mathrm{M}_{}`$) of iron. We shall follow this convention throughout. Following arguments made by Murray et al. (2001), protoplanetary material is roughly 20% iron by mass. Note that we are essentially changing the mass fraction of heavy elements (Z) in the models and the accretion of 1 $`\mathrm{M}_{}`$ of iron corresponds to the addition of 5 $`\mathrm{M}_{}`$ of Z.

Stars which accret rocky material will experience an increase in heavy element composition in the outer layers. Our approach is to rescale the composition of the SML without altering its mass. Starting with an amount of iron, we take the ratio of the mass of rocky material to the mass of the SML ($`\mathrm{Z}_{\mathrm{rock}}`$) and add this to the initial mass fraction of metals ($`\mathrm{Z}_\mathrm{i}`$). The new mass fraction of metals is

$$Z_f=Z_i+Z_{rock}$$
(1)
A new heavy element mass fraction is calculated and the abundances of hydrogen (X) and helium (Y) are adjusted. We rescale X and Y as follows

$$\frac{X_f}{X_i}=1\frac{Z_{rock}}{1Z_i}$$
(2)
$$\frac{Y_f}{Y_i}=1\frac{Z_{rock}}{1Z_i}$$
(3)
The $`i`$\- and $`f`$-subscripts refer to initial and final compositions, respectively.

In models which accrete more than about $`1\mathrm{M}_{}`$ of iron we found it necessary to increase the metallicity over multiple time steps. Spreading the pollution over several time steps allows us to avoid the problems which occur when the conditions in the model change too quickly for the code to adjust. In addition, it is more in line with current theories which describe the time scales of stellar pollution as taking place over tens of millions of years. Neither altering the time scale nor changing the exact age at which rocky material is accreted (after the zero-age main sequence, ZAMS) has a significant impact on evolution. Unless otherwise noted, we begin to add the rocky material when the model reaches the ZAMS111 We note that if the pollution were to occur significantly before the the star reaches the ZAMS, then accreted material would be mixed throughout the star, as stars in the early pre main sequence phase of evolution are fully convective. For stellar pollution to have observational consequences, it must occur after the star’s convection zone has thinned near the end of the pre main sequence. . Since all stars considered metal rich compared to the Sun we rescale Y and Z (from Solar values) in the ratio $`\mathrm{\Delta }`$Y/$`\mathrm{\Delta }`$Z = 2 (Chaboyer et al., 1999).

Finally, to calculate \[Fe/H\] it is sufficient to use \[Fe/H\]=log($`mathrmZ_{}/Z_{}`$) since the relative change in X is insignificant compared to the relative change in Z. However, this assumes the star in question has the same abundance of iron as the Sun. The accreted material is roughly 20% iron by mass, while the Sun is only about 0.14% iron by mass (iron comprises 7.3% of $`\mathrm{Z}_{}`$ by mass, and $`\mathrm{Z}_{}`$=0.019). We account for this difference by introducing a factor to the \[Fe/H\] formula which accounts for the increased abundance of iron on the surface of the star.

$$\left[\frac{Fe}{H}\right]=log\left(\frac{Z_{}}{Z_{}}\right)+log\left(\frac{M_P(Fe)+M_{SML}(Fe)}{M_{SML}(Fe)}\right)$$
(4)
Where $`\mathrm{M}_\mathrm{P}`$(Fe) is the mass of iron in the polluting material and $`\mathrm{M}_{\mathrm{SML}}`$(Fe) is the mass of iron in the SML.

## 3 Basic Results

The most significant impact of stellar pollution on the evolution of star is the difference in effective temperature. Given two main sequence stars of similar mass, age, and surface metallicity a polluted star will have a higher effective temperature than an unpolluted star (for a specific example see $`\mathrm{\S }`$4).

The altered behavior can be interpreted as follows: the unpolluted model has a higher mass fraction of heavy elements in the deep interior (below the SML). This leads to a higher opacity in the deep interior of the unpolluted model, making it more difficult for energy to flow outward. The unpolluted model is thus physically larger and cooler than the polluted model.

The effect of adding a given amount of iron to a star decreases with SML mass because mixing the polluting material over a larger region dilutes the effect.

## 4 $`\tau `$ Bootis

Recent evidence that SWPs are metal-rich (Laws & Gonzalez, 2002) compared to a sample of field stars lends support to the concept that these stars have accreted metal-rich material before the planetary system reached equilibrium. If the accretion occurred after the star left the pre main sequence convective phase the effects should be noticeable in the form of high metallicity.

$`\tau `$ Bootis has a short-period, giant planet. It also has a high metallicity, \[Fe/H\]=0.32$`\pm `$0.06 (Gonzalez & Laws, 2000). To illustrate the effects stellar pollution may have when stellar models are used to infer properties of real stars, observational data for $`\tau `$ Boo was fitted by standard and polluted models.

A series of polluted models was constructed with masses between 1.30 and 1.40 $`\mathrm{M}_{}`$ adding 1-5 $`\mathrm{M}_{}`$ of iron when the model reaches the ZAMS. Figure 1 shows the evolutionary tracks of two 1.36 $`\mathrm{M}_{}`$ stellar models. The solid line is a model polluted with 2 $`\mathrm{M}_{}`$ of iron. The dashed line is an unpolluted model with bulk metallicity chosen to match the surface metallicity of the polluted model. The polluted model is marked at several points to show how the model evolves with time. The unpolluted model evolves in the same manner.

Table 1 contains a small sample of polluted models with luminosities, temperatures, and ages at closest approach to the observational data. For comparison we include models constructed using the same evolution code but which have not been polluted.

The data in Table 1 show that, for a given mass, the polluted models are slightly older than the unpolluted models: from a few percent at the high mass end (around 1.40 $`\mathrm{M}_{}`$) up to 20% at the low mass end (1.34 $`\mathrm{M}_{}`$). The polluted models have a best fit mass which is 0.02 - 0.03 $`\mathrm{M}_{}`$ lower than the best fit from the unpolluted models. The purpose of this table is only to illustrate the general trends which pollution has on stellar models, and not to constrain the stellar mass of, or amount of pollution, experienced by, $`\tau `$ Boo.

## 5 The Hyades

Murray et al. (2001) estimate that stars in the solar neighborhood have accreted 0.5$`\pm `$0.25 $`\mathrm{M}_{}`$ of iron. The sample used to make this estimate was comprised of field stars which have a broad range of ages and bulk metallicities. The variation in age and metallicity of the sample, along with the fact that younger stars tend to be more metal-rich than older stars, complicated the analysis of Murray et al. (2001).

Stars in the Hyades are likely to have formed at the same time from material of the same composition. As a result, Hyades stars are likely to yield better constraints than field stars on stellar pollution. The relatively young age of stars in the Hyades, combined with a uniform composition, should imply a tight main sequence. However, planet formation and accretion of metal-rich material will vary from one star to another which will cause scatter in the color-magnitude diagram (CMD).

Quillen (2002) proposes that the scatter of Hyades stars about a theoretical isochrone can be used to limit the level of pollution experienced by Hyades stars. Quillen (2002) did not have access to polluted stellar evolution models, and so based her proposal on some plausible assumptions that the effects of stellar pollution would have on the evolution of stars. With the modified evolution code presented in this paper it is possible to directly test Quillen’s proposal.

As a preliminary test, three isochrones were constructed with varying levels of pollution. Each isochrone consists of models ranging in mass from 0.6 to 1.8 $`\mathrm{M}_{}`$ in steps of 0.05 $`\mathrm{M}_{}`$. Each model begins on the pre main sequence, receives the polluting material at about 150 Myr (at which point all models have reached or passed the ZAMS), and is evolved to 650 Myr. Uncertainty in the final age of the isochrones by $`\pm `$50 Myr alters the scatter by an amount roughly an order of magnitude below the scatter due to pollution effects.

The models in each isochrone have a common initial $`Z=0.024`$ leading to $`<`$\[Fe/H\]$`>`$0.13 in accord with Boesgaard & Friel (1990) who find $`<`$\[Fe/H\]$`>`$ = 0.127$`\pm `$0.022 in the range of 6000 - 7000 K. In the “mean” polluted isochrone, models are polluted with 0.5 $`\mathrm{M}_{}`$ of iron. The other two isochrones correspond to pollution amounts of $`2\sigma `$ above and below the mean, respectively 1.0 and 0.0 $`\mathrm{M}_{}`$ of iron.

The top panel in Figure 2 is a CMD showing the main sequence from 0.6 – 1.8 $`\mathrm{M}_{}`$ of the Hipparcos Hyades sample of de Bruijne et al. (2001) (also used by Quillen (2002) ). Buser & Kurucz (1992) color tables were used to convert the models from luminosity and effective temperature to $`\mathrm{M}_\mathrm{V}`$ and B–V. In addition, the CMD plot includes the mean polluted isochrone and a fit to the data made using the LOWESS (Locally Weighted Regression) technique (Cleveland, 1979; Gebhardt et al., 1994, 1995). In the plot the dashed line is the LOWESS fit and the solid line is the mean isochrone. The LOWESS fit is weighted by the the local points, and represents the best fit to the data that is possible with a single line.

The bottom panel of Figure 2 shows the scatter of the Hyades sample about the LOWESS fit. Superimposed are the differences between the $`2\sigma `$ isochrones and the mean isochrone. The difference between the mean and $`+2\sigma `$ isochrone is represented by the dashed line; the difference between the mean and the $`2\sigma `$ isochrone is represented by the dotted line. Note that the dispersion of the data points is smaller than that predicted by the $`\pm 2\sigma `$ models, showing that the dispersion predicted by the stellar pollution models is not ruled out by the data.

The LOWESS technique is used to fit the Hipparcos Hyades sample because the mean polluted isochrone diverges slightly from the main sequence for $`\mathrm{B}\mathrm{V}`$ greater than $``$0.65 mag. This difference between the models and the data is likely due to problems with the treatment of convection, or in the color calibration of the models. We are not interested in studying this problem in this paper, rather the goal is to compare the scatter in the Hyades stars about a smooth main sequence to the predicted scatter caused by pollution. The LOWESS fit allows a comparison between the observed scatter and the theoretical prediction which is not hampered by large scale errors introduced by the color conversion and/or treatment of convection. In general, the theoretical isochrone is a reasonable fit to the data, with the the only consistent deviation between the LOWESS fit and the mean polluted isochrone occurring for $`\mathrm{B}\mathrm{V}0.65`$ mag.

The fit of the mean polluted isochrone to the observational data has a reduced $`\chi ^2`$ of $``$3.3 and the LOWESS fit has a reduced $`\chi ^2`$ of 2.1. The mean polluted isochrone has the lowest $`\chi ^2`$ of all the isochrones constructed and is therefore the standard. The LOWESS fit is the best is possible single line fit to the data, and the large value of $`\chi ^2`$ implies that their is an inherent dispersion in the data.

As a more stringent test, the method of Gebhardt et al. (1994, 1995) for determining the velocity dispersion as a function of radius within a globular cluster has been adapted to study the scatter in the CMD. In this method, the value of the error bar for each star in the Hyades sample is removed in quadrature from the absolute value of $`\mathrm{\Delta }\mathrm{M}_\mathrm{V}`$ (scatter) value of that star. The square root of the resulting value is the $`\mathrm{M}_\mathrm{V}`$ dispersion, analogous to the velocity dispersion in the globular cluster studies of Gebhardt et al. (1994, 1995). If the size of the error bar is larger than the absolute value of $`\mathrm{\Delta }\mathrm{M}_\mathrm{V}`$ for a given data point, the $`\mathrm{M}_\mathrm{V}`$ dispersion is set to zero for that point. The LOWESS fit is then determined for $`\mathrm{\Delta }\mathrm{M}_\mathrm{V}`$ as a function of B V. This produces “dispersion profile” for the Hyades sample. Note that the dispersion is a non-negative quantity, any information about the scatter being above or below the mean isochrone is lost.

The final step is to estimate the error in the derived dispersion profile. This is done in a manner somewhat analogous to the standard bootstrap re-sampling technique for error analysis. In particular, for each data point one randomly selects a new value for $`\mathrm{M}_\mathrm{V}`$ using a Gaussian distribution, with a mean equal to the original value of $`\mathrm{M}_\mathrm{V}`$ and a one-$`\sigma `$ width equal to the quoted one-sigma error in the data point. The process is performed one thousand times for each star and a mean dispersion profile with confidence bands is calculated from the results.

Having obtained the dispersion profile for the Hyades data, it is necessary to put the theoretical prediction in a similar form for comparison. Using the mean polluted isochrone as a guide, the mass of each star with $`0.2\mathrm{B}\mathrm{V}1.0`$ in the Hyades sample was determined. Stellar models with these particular values of mass were evolved 1000 times, where each of the stellar evolution models had a random amount of stellar pollution added to it at 150 Myr. The amount of the polluting material was drawn from a Gaussian distribution with standard deviation equal to half the mean. This is done for two cases, means equal to 0.5 and 1.0 $`\mathrm{M}_{}`$ of iron.

The result of this simulation is a set of 2000 theoretical CMDs, where each CMD contains the same number of data points as the Hyades data set, but for each each of the stars in the CMD has been polluted by a varying amount. To compare with the dispersion profile in the data, the individual theoretical CMD is treated as a “mock” data set and put through the exact same procedure as the real data set to generate a dispersion profile of $`\mathrm{M}_\mathrm{V}`$ as a function of B V. This was done each of the 2000 theoretical CMDs, first for 1000 CMDs which had a mean amount of pollution equal to 0.5 $`\mathrm{M}_{}`$ of iron and then for the 1000 CMDs with a mean pollution of 1.0 $`\mathrm{M}_{}`$ of iron.

Figure 3 is a plot of the dispersion in the data and in the averaged mock data sets. In the both panels, the thick line is the mean dispersion in the data and the thin line is the mean dispersion in the averaged mock data sets. The 1.0 $`\mathrm{M}_{}`$ case is the top panel; the 0.5 $`\mathrm{M}_{}`$ case is the bottom panel. In each panel, 90% confidence bands are drawn about the mean dispersion in the Hyades data (dotted lines) and the averaged mock data sets (dashed lines). From Figure 3 it is clear that neither the 0.5 nor the 1.0 $`\mathrm{M}_{}`$ case is ruled out by the observed dispersion in the Hyades data at the 90% confidence level. However, looking at the trends apparent in these two cases, it is clear that pollution at the mean level of $``$1.5 $`\mathrm{M}_{}`$ of iron is ruled out by the data.

## 6 Conclusions

Stars with giant planets in tight orbits, such as $`\tau `$ Boo, are likely to have accreted $``$5 $`\mathrm{M}_{}`$ of iron (Murray & Chaboyer, 2002). A method is presented for modifying standard stellar evolution models to realistically account for this stellar pollution effects. Such polluted stars will exhibit a significantly hotter surface than stars of similar mass and age with a correspondingly high bulk metallicity. Polluted stellar evolution models may prove useful in studying the evolution of stars which prove difficult for standard stellar evolution models.

On the other hand, most stars in the solar neighborhood have not been observed to have giant planets, and appear to have accreted an average of 0.5 $`\mathrm{M}_{}`$ of iron (Murray et al., 2001). At such a low level, pollution effects in the CMD are found to be very small. As an example, we look at pollution effects in the Hyades open cluster. Stars in the Hyades all formed from the same material at roughly the same time. At 650 Myr the stars we have chosen are all on the main sequence and so past the era when pollution would have occurred. Monte Carlo simulations show that the dispersion in $`\mathrm{\Delta }\mathrm{M}_\mathrm{V}`$ of the Hyades about a polluted isochrone is consistent with the observed data, assuming that the pollution in the Hyades is at a similar level to the field star sample. However, if the level of pollution in the Hyades were a factor of 3 or more greater than that in the field stars, the analysis performed would be able to rule it out. We are currently investigating the mass-\[Fe/H\] relationship in the Hyades to see if this data can put stronger constraints on stellar pollution.

Research supported in part by a NSF CAREER grant 0094231 to BCC. BCC is a Cottrell Scholar of the Research Corporation. We would like to thank Karl Gebhardt for guidance and computer programs relating to the LOWESS technique.


text:
# 𝐽/𝜓 Polarization in 800-GeV limit-from𝑝-Cu Interactions

Present address: \]University of Illinois, Urbana, IL 61801

Present address: \]University of Cincinnati, Cincinnati, OH 45221

Present address: \]Florida State University, Tallahassee, FL 32306

Present address: \]University of Illinois, Urbana, IL 61801

On leave from \]Kurchatov Institute, Moscow, Russia

Present address: \]Indiana University, Bloomington, IN 61801

FNAL E866/NuSea Collaboration

## Abstract

We present measurements of the polarization of the $`J/\psi `$ produced in 800-GeV proton interactions with a copper target. Polarization of the $`J/\psi `$ is sensitive to the $`c\overline{c}`$ production and hadronization processes. A longitudinal polarization is observed at large $`x_F`$, while at small $`x_F`$ the state is produced essentially unpolarized or slightly transversely polarized. No significant variation of the polarization is observed versus $`p_T`$.

J/psi(3100) polarization, J/psi(3100) hadroproduction, angular distribution

A detailed understanding of the production mechanism of charmonium is important for ongoing research at the Relativistic Heavy Ion Collider, where this process can play an important role in the search for quark-gluon plasma formation Matsui and Satz (1986), as well as for the investigation of the gluon contribution to the proton spin structure Robinett (1991); Teryaev and Tkabladze (1997). In the non-relativistic QCD formalism (NRQCD), charmonium polarization can probe Bodwin et al. (1995) details of the production process Cho and Wise (1995) that are perturbatively calculable.

NRQCD is an effective field theory that approximates the full QCD Lagrangian for large quark masses. Uncalculable matrix elements for the production of various Fock-space components of the quarkonium wavefunction are ranked, using simple scaling rules, according to their order in $`v`$, the relative quark-antiquark velocity in the quarkonium rest frame; for charmonium, this has a value of $`0.5`$. The leading matrix elements can then be extracted from fits to experimental data and used in calculating other processes.

A determination Kniehl and Kramer (1999) of the various matrix elements from available high-energy data defines to a large extent the expected production properties of charmonium, one of which is the polarization. For a spin-1 particle, given the fraction $`\xi `$ of particles produced in the $`j_z=0`$ (“longitudinal”) state, we can define the polarization $`\lambda =(13\xi )/(1+\xi )`$; $`\lambda `$ is positive (negative) for transverse (longitudinal) polarization. While several of the intermediate $`c\overline{c}`$ states are color-octet states and must be followed by multiple gluon emission before a physical charmonium state is produced, heavy-quark symmetry implies that gluon radiation leaves the quark spins unchanged, providing definite predictions for the spin state of the final charmonium Cho and Wise (1995).

The polarization of the $`J/\psi `$ has been measured with relatively high statistical precision only in fixed-target experiments, in pion Biino et al. (1987) and proton Alexopoulos et al. (1997); Gribushin et al. (2000) interactions with solid nuclear targets. No significant polarization has been seen in either, with the exception of an intriguing large longitudinal polarization at the highest $`x_F`$ value of the pion-induced data. Measurements in collider experiments, at energies and transverse momenta where theoretical calculations should be more robust and which are free of potential complications from nuclear effects, suffer from low statistics Affolder et al. (2000). Finally, in a recent study of fixed-target bottomonium production Brown et al. (2001), the $`\mathrm{{\rm Y}}(1S)`$ was found to be largely unpolarized, while the $`\mathrm{{\rm Y}}(2S,3S)`$ states had strong transverse polarization.

In this experiment, we accumulated Chang (1999) a much larger sample — approximately 9 million — of reconstructed $`J/\psi `$ than any previous study, in interactions of an 800-GeV/$`c`$ proton beam with a copper target. $`J/\psi `$ decays were measured using the Fermilab Meson-East dimuon spectrometer Moreno et al. (1991), which consisted of three dipole magnets, SM0, SM12, and SM3, and three stations of drift chambers and trigger scintillator hodoscopes. Data were collected during a month-long dedicated run in which the copper beam dump inside SM12 was used in place of a target; the first dipole was switched off during this run. A copper absorber filtered out all hadrons and electrons produced in the interactions, allowing only muons to enter the spectrometer. Additional hadron rejection was obtained using a muon identifier, consisting of proportional tubes and scintillator hodoscopes behind a thick absorber at the downstream end of the apparatus. Events were recorded when the trigger condition of two oppositely-charged muons was satisfied.

Events were reconstructed offline from the recorded hits in the drift chambers. Tracks were traced back through the magnetic fields to the dump/target, where a vertex was formed, consistent with the beam position. Energy losses and multiple Coulomb scattering in the absorber and the dump were taken into account in the traceback. The momentum of the muons was determined by their bending in SM3 and this was used to calculate their trajectories in SM12.

The invariant mass of the muon pair was calculated from the muon momenta and opening angle at the vertex. Since the energy loss and multiple scattering are only known on average and not on an event-by-event basis, the mass resolution, typically $`500`$ MeV (FWHM), was not sufficient to separate the $`J/\psi `$ and $`\psi ^{}`$ peaks in this run using this extended target. While the two charmonium states cannot be separated, it is estimated that the $`\psi ^{}`$ contributes only about 1% to the total event count, based on the relative production cross sections and branching ratios into muons. Therefore, the $`\psi ^{}`$ is not considered in the following discussion.

The mass distributions were plotted in bins in the Feynman-$`x`$ variable $`x_F`$, transverse momentum $`p_T`$, and the dimuon polar angle $`\vartheta `$ in the dimuon rest frame. We use the Collins-Soper frame Collins and Soper (1977). This is identical to the Gottfried-Jackson frame Gottfried and Jackson (1964), used in several earlier experiments, for $`p_T=0`$, and to a very good approximation equivalent even at the highest $`p_T`$ values in this experiment.

The invariant mass distribution in each $`(x_F,p_T)`$-bin were fitted to a Gaussian peak plus an exponential or polynomial background. All of the parameters, including the mass, were determined by the fit. The number of events under the peak gave the combined $`J/\psi `$ and $`\psi ^{}`$ triple-differential, unnormalized cross section in $`x_F`$, $`p_T`$, and $`\vartheta `$. Distributions as a function of $`\vartheta `$ were then formed in each $`x_F`$ and $`p_T`$ bin.

The spectrometer and trigger acceptance was calculated with the help of a Monte Carlo simulation of the $`J/\psi `$ production process, which included all the measured magnetic fields and detector efficiencies and geometry. The known properties of $`J/\psi `$ production from previous experiments were used. However, because of the binning in $`x_F`$ and $`p_T`$, exact knowledge of the form of the cross section as a function of these two variables was not crucial. Events were generated with a flat $`\vartheta `$ distribution, corresponding to unpolarized production. Simulated data were passed through the same analysis chain as the real data. The Monte Carlo reproduced quite accurately the main features of the data, including the mass and vertex resolutions.

The azimuthal-angle ($`\varphi `$) dependence of the production cross section was assumed to be flat in the Monte Carlo, since the two interacting hadrons were unpolarized. The corresponding measured acceptance-corrected distribution was essentially flat, consistent with a 2% uncertainty in the direction of the incoming beam. The decay distribution was also assumed to be flat in $`\varphi `$, consistent with previous results Biino et al. (1987). No attempt was made to extract the $`\varphi `$ dependence of the decay in this analysis.

The $`\vartheta `$ distributions of the data in all bins were divided by the corresponding ones of the simulation, resulting in acceptance-corrected $`\vartheta `$ distributions. These were then used to calculate the polarization in $`x_F`$ and $`p_T`$ bins, according to the formula $`d\sigma /d\mathrm{cos}\vartheta =A(1+\lambda \mathrm{cos}^2\vartheta )`$, where the normalization constants $`A`$ were left free. Results were obtained separately for two experimental runs with the current in the SM12 magnet set to 2040 and 2800 A respectively, resulting in substantially different acceptances.

Figure 1 shows the polarization parameter $`\lambda `$ as a function of $`x_F`$ in four $`p_T`$ bins; it can be seen that results from the two magnet settings are in reasonable agreement, giving some confidence that the acceptance is understood and providing an estimate for the magnitude of the relevant systematic uncertainties. The $`x_F`$ dependence of $`\lambda `$ appears to be independent of $`p_T`$.

Systematic errors from various sources were considered. Inexact knowledge of the $`p_T`$ dependence of the production cross section, coupled with a strong $`p_T`$ dependence of the acceptance versus decay angle $`\vartheta `$, led to an uncertainty of $`\pm 0.06`$ in $`\lambda `$, independent of $`x_F`$. Additional contributions included mass-peak fitting errors (0.04–0.08, depending on the $`x_F`$ bin) and uncertainties in the exact position (0.02) and angle (0.02–0.04) of the incoming beam and in the fields in the analyzing magnets (0.01). Uncertainties from various sources are largely uncorrelated; they were added in quadrature for the overall systematic error.

Figure 2 presents the polarization parameter $`\lambda `$ as a function of $`p_T`$ for two $`x_F`$ ranges, where the two data sets have been statistically combined. The two intervals in $`x_F`$ approximately correspond to regions were the gluon-gluon and quark-antiquark processes are dominant. No significant $`p_T`$ dependence is seen in either region after the systematic errors are taken into account. At large $`p_T`$, charmonium production is understood Braaten and Yuan (1993) to be dominated by gluon bremsstrahlung with subsequent fragmentation into a $`c\overline{c}`$ pair. In this case, the charmonium state is expected Cho and Wise (1995); Braaten et al. (2000) to retain to a large degree the transverse polarization of the high-$`p_T`$, on-shell gluon. However, the $`p_T`$ range of the experiment, extending to 4 GeV, may not be sufficient to see clearly such an effect. It must be noted that this effect is also not observed at the high $`p_T`$ values available at collider energies Affolder et al. (2000), where, if anything, polarization appears to be longitudinal. In the following, we assume there is no significant $`p_T`$ dependence and the results are presented integrated over all $`p_T`$.

Figure 3 shows the polarization $`\lambda `$ as a function of $`x_F`$ for the combined data set. The one-sigma systematic uncertainty is shown as a dark band. Also shown are the results previously obtained by the Chicago-Iowa-Princeton collaboration Biino et al. (1987), which used a 252-GeV pion beam and a tungsten target. Our results suggest that the $`J/\psi `$ is produced slightly transversely polarized at small-to-intermediate $`x_F`$. Within the systematic uncertainties the results are consistent with the CIP experiment Biino et al. (1987), which saw no polarization in this range, albeit with much larger statistical uncertainties. It should be noted that in this range the $`q\overline{q}`$ annihilation process is expected to play a more important role with a pion beam than with a proton beam, where $`gg`$ fusion dominates. For $`x_F>0.6`$, the polarization turns to longitudinal, outside the margins of the systematic error. This is similar to a pattern seen in the CIP experiment, although the behavior appears smoother as a function of $`x_F`$ and begins at smaller $`x_F`$ than the sudden turn-over near $`x_F0.85`$ in the latter. However, the results are not incompatible. It is also interesting to note that the clear change in polarization in this experiment roughly coincides with the transition from gluon-fusion dominance to quark-annihilation dominance.

Integrated over the entire $`x_F`$ range, the measured polarization is $`\lambda =0.069\pm 0.004\pm 0.08`$ (statistical and systematic errors). This small value, consistent with no polarization if systematic errors are taken into account, is in agreement with previous proton-beam experiments Alexopoulos et al. (1997); Gribushin et al. (2000), which had insufficient statistics to study the $`x_F`$ dependence.

As a cross-check of the analysis, the polarization of the Drell-Yan continuum was also studied using the same technique. This can be done only for dimuon invariant masses greater than 4 GeV, since the $`J/\psi `$ peak dominates the spectrum at lower masses. An additional complication arises from random coincidences of uncorrelated pairs of muons, which in the case of the $`J/\psi `$ were removed by the fitting procedure. These were estimated and subtracted from the Drell-Yan sample by studying the distributions of same-charge muon pairs. A fit to the formula $`1+\lambda \mathrm{cos}^2\vartheta `$ for dimuon masses from 4 to 7 GeV gave $`\lambda =0.98\pm 0.04`$, in good agreement with the 100% transverse polarization expected and previously observed Brown et al. (2001) in the Drell-Yan process. While the kinematic range is not the same as for the $`J/\psi `$ measurement, this agreement with expectations increases confidence in the soundness of the $`J/\psi `$ analysis. The level of statistics were not adequate to investigate any $`x_F`$ dependence of the Drell-Yan polarization.

The results are not in agreement with published predictions for $`J/\psi `$ polarization based on NRQCD, which are in the range $`0.31<\lambda <0.63`$ Beneke and Rothstein (1996), or with similar predictions Vanttinen et al. (1995) based on an early model that considers only color-singlet intermediate states. The small positive values at $`x_F0.5`$ are very similar to the results obtained for $`\mathrm{{\rm Y}}(1S)`$ by this experiment in a similar $`p_T`$ range, and in sharp contrast to the essentially 100% transverse polarization of the (unresolved) $`\mathrm{{\rm Y}}(2S)`$ and $`\mathrm{{\rm Y}}(3S)`$ states Brown et al. (2001).

To understand the behavior of the polarization, all sources of $`J/\psi `$ production must be considered. While $`b`$-quark production is not a major source at these energies, almost half of the produced $`J/\psi `$’s are the decay products of higher-mass charmonium resonances, mainly the $`{}_{}{}^{3}P_{J}^{}`$ states $`\chi _{cJ}`$. Feed-down from $`\chi _{c2}`$, produced exclusively in the $`J_z=2`$ state from gluon-gluon fusion, results in 100% transverse $`J/\psi `$’s, increasing the observed values of $`\lambda `$. The turn-over at high $`x_F`$ may reflect the transition, at $`x_F0.6`$, to the quark-annihilation graph Kowitt et al. (1994), which produces a mix of $`\chi _{c2}`$ spin states. As for the overall level, it cannot be explained without a substantial contribution from $`\chi _{c1}`$, which in general produces longitudinal $`J/\psi `$’s. Beneke and Rothstein Beneke and Rothstein (1996) originally estimated the $`\chi _{c1}`$ contribution to be about 10 times smaller than that from $`\chi _{c2}`$, but it was later found experimentally Alexopoulos et al. (2000) to be of comparable size. NRQCD calculations which retain higher orders in $`v`$ can also accommodate substantial $`\chi _{c1}`$ contributions Beneke (1997). It would probably require a contribution near the experimental and theoretical upper limits in order to obtain $`\lambda `$ values as small as experimentally measured, within the context of NRQCD.

In addition, it must be noted that no existing calculation takes into account nuclear effects, which strongly affect the production cross section Leitch et al. (2000). The formation length of the $`J/\psi `$ at these energies is generally longer than the nuclear size Farrar et al. (1990) and it is conceivable that color-singlet and -octet components of the wavefunction are absorbed differently while propagating through the nuclear medium, resulting in a different mix of Fock states compared to the free-nucleon production, thus altering the polarization. Furthermore, the polarization of the $`\chi _c`$ states themselves is predicted to have a nuclear dependence Gerland et al. (1998) which would feed down to the $`J/\psi `$ polarization. Finally, at the highest values of $`x_F`$, higher-twist effects may become important Vanttinen et al. (1995).

###### Acknowledgements.

This work was supported in part by the U.S. Department of Energy.