text:
# On the Melting of Bosonic Stripes

## Abstract

We use quantum Monte Carlo simulations to determine the finite temperature phase diagram and to investigate the thermal and quantum melting of stripe phases in a two-dimensional hard-core boson model. At half filling and low temperatures the stripes melt at a first order transition. In the doped system, the melting transitions of the smectic phase at high temperatures and the superfluid smectic (supersolid) phase at low temperatures are either very weakly first order, or of second order with no clear indications for an intermediate nematic phase.

Stripe phases of lattice models with broken rotational and translational symmetry can melt in two qualitatively different ways. One scenario is that both symmetries can be restored at a single first or second order transition, and the stripes melt directly into a normal fluid or superfluid phase. The other scenario is that first the translational symmetry is restored when the striped solid melts into a nematic (liquid crystal) phase with broken rotational symmetry. The rotational symmetry is then restored in a second melting transition of the nematic phase.

This quantum lattice problem shows similarities to the long standing problem of the melting of a two dimensional crystal into a continuum model, where there is also either a first order melting or the Kosterlitz-Thouless-Halperin-Nelson-Young (KTHNY) scenario of two Kosterlitz-Thouless transitions with an intervening hexactic phase KTHNY . Clear results for the classical version of this continuum problem were obtained only recently in a simple model of hard disks HardDisk .

Current interest in quantum mechanical stripe phases and their melting stems from the experimental observation of stripes in some high-$`T_c`$ superconductors Tranquada  and from the question whether and how they are related to the occurrence of high temperature superconductivity. Numerically, stripe phases have been found to be competitive ground states of $`t`$-$`J`$-like models tsune ; tjstripes . Analytically some theories of high temperature superconductivity are closely linked to the existence of stripe and nematic phases Kivelson . While it is hard to study stripe phases in strongly correlated fermionic models, because of the negative sign problem of quantum Monte Carlo, we can more accurately investigate bosonic stripes using modern quantum Monte Carlo algorithms loop ; Sandvik . Such bosonic models can appear as effective low energy models neglecting nodal quasiparticles tsune ; bosonmodels . Like the cuprates, these bosonic models show competition and in some models coexistence of superfluidity and charge order.

In this Letter we focus on the simplest bosonic quantum model exhibiting stripe order and determine its finite temperature phase diagram. We carefully investigate its thermal and quantum melting transitions to address the question of how stripe melting occurs in a quantum model of bosonic stripes. The particular hard-core boson Hubbard Hamiltonian we study on a two-dimensional square lattice is

$`H=t{\displaystyle \underset{𝐢,𝐣}{}}(a_𝐢^{}a_𝐣+a_𝐣^{}a_𝐢)+V_2{\displaystyle \underset{𝐢,𝐣}{}}n_𝐢n_𝐣\mu {\displaystyle \underset{𝐢}{}}n_𝐢,`$ (1)
where $`a_𝐢^{}`$ $`(a_𝐢)`$ is the creation (annihilation) operator for hard-core bosons, $`n_𝐢=a_𝐢^{}a_𝐢`$ the number operator, $`V_20`$ the next nearest neighbor Coulomb repulsions and $`\mu `$ the chemical potential. The Hamiltonian (1) is equivalent to an anisotropic spin-1/2 model with a nearest neighbor $`XY`$-like interaction $`J_{xy}=2t`$ and next nearest neighbor Ising interaction $`J_z^{}=V_2`$ in a magnetic field $`h=2V_2\mu `$.

In Fig. 1 we review the ground-state phase diagram of this model which was previously studied by mean-field, renormalization group and local-update quantum Monte Carlo simulations pich ; frey ; batrouni3 . At half filling (density $`\rho =1/2`$ at $`\mu =2V_2`$) the ground state is a smectic for $`t/V_2<0.446\pm 0.006`$ and the quantum melting transition at low temperatures was found to be of first order (translation symmetry is broken only in one dimension, perpendicular to the stripes, hence this phase is a smectic) .

Doping this smectic stripe crystal a stable “supersolid” phase with coexisting stripe order and superfluidity was found, in which the vacancies doped into the stripes form a superfluid.

Before going into details we summarize our key results by presenting the finite temperature phase diagrams along the lines indicated in the ground state phase diagram Fig. 1. At half filling (see Fig. 2a) we find that the transition between smectic and superfluid is of first order at all temperatures. At higher temperatures the transition between smectic and normal fluid is of first order around $`t/V_20.5`$, but then becomes very weakly first order, before turning second order as $`t/V_20`$. In the doped system (see Fig. 2b) coexistence between the smectic phase and superfluidity can be observed at low temperatures around half filling. Here the phase transitions are all either very weakly first order or of second order. Rotational and translational symmetry breaking occur, within the accuracy of our simulations, at the same point.

To determine the finite-temperature phase diagram we have used quantum Monte Carlo simulations, employing a directed loop quantum Monte Carlo algorithm in the in the stochastic series expansion (SSE) representation Sandvik . The simulations were performed in the grand canonical ensemble and do not suffer from any systematic errors apart from finite size effects. In contrast to the loop algorithm loop  the directed loop algorithm remains effective away from half filling where the loop algorithm slows down exponentially.

We determine stripe order by measuring the order parameter

$$O_S=S_n(\pi ,0)+S_n(0,\pi )$$
(2)
where $`S_n(k_x,k_y)`$ is the charge structure factor at the wave vector $`(k_x,k_y)`$. To investigate a nematic phase, we have to look for rotational symmetry breaking in the kinetic energy or in the local charge correlations, using as order parameters

$`O_k={\displaystyle \frac{1}{V}}{\displaystyle \underset{(x,y)}{}}a_{(x,y)}^{}a_{(x+1,y)}^{}a_{(x,y)}^{}a_{(x,y+1)}^{}+H.c.`$ (3)
or alternatively

$`O_N={\displaystyle \underset{(x,y)}{}}n_{(x,y)}n_{(x+1,y)}n_{(x,y)}n_{(x,y+1)}`$ (4)
where $`n_{(x,y)}=a_{(x,y)}^{}a_{(x,y)}^{}`$ is the boson number operator at lattice site $`(x,y)`$ and $`H.c.`$ denotes the Hermitian conjugate.

The extent of the superfluid phase is determined by measuring the superfluid (number) density

$`\rho _s=mTW^2`$ (5)
where $`W`$ is the winding number in one of the directions and $`m`$ the mass of a boson. $`\rho _s`$ is finite in the superfluid phase and jumps to zero from a universal value $`\frac{2}{\pi }mT_c`$ at the Kosterlitz-Thouless transition.

We now discuss the phase diagrams and the nature of the phase transitions in more detail, starting with the half filled system at $`\mu =2V_2`$, where a first order quantum phase transition was found at $`V_2/t=2.24\pm 0.03`$, improving the previous estimate batrouni3 . When raising the temperature we find that the transition remains of first order as the histograms for both the stripe order parameter $`O_S`$ and the rotational symmetry breaking order parameter $`O_N`$ show two clearly separated peaks corresponding to the two coexisting phases (see Fig. 3). Careful finite size scaling of the peak distances on system sizes $`L=8`$ up to $`L=32`$ show that the two peaks stay separate as expected for a first order transition and do not merge as in a second order one. For both order parameters we find the same phase transition line, confirming a melting in a single strong first order transition up to at least $`T/t=2/3`$ \[from point (D) to (B) in Fig. 2a)\].

To investigate whether the stripes melt into a superfluid or normal fluid we measure the superfluid density $`\rho _s`$ at the coexistence line for configurations in the fluid phase (determined by their value of $`O_S`$). At a temperature $`T/t=2/3`$ \[point (B)\], $`\rho _s`$ drops below the universal jump value $`\frac{2}{\pi }mT`$ as a function of system size $`L`$ and will thus scale to zero in the infinite system. Hence the transition at this temperature is of first order to a normal fluid phase in contrast to $`T=0.5t`$ \[point (C)\] where $`\rho _s`$ is well above the universal value for all system sizes indicating a direct first order transition between the smectic and a superfluid phase nelson . There is thus a tricritical point between points $`B`$ and $`C`$ where the Kosterlitz-Thouless phase transition of the superfluid becomes first order.

At higher temperatures and larger repulsion $`V_2`$ it becomes harder to determine the order of the phase transition. At $`V_2=3t`$ no double-peak structure could be seen in the histograms for lattice sizes up to $`L=32`$ and the transition is thus either second order or a very weak first order transition. In the limit of infinite repulsion $`V_2`$ the lattice decouples into two sub-lattices, each of which is equivalent to an Ising antiferromagnet with a second order melting transition. We thus explore the possibility of a second order transition.

For a second order transition the $`4`$-th order cumulant ratios $`C_4=1O^4/3O^2^2`$ have a size independent value at the transition point, which can be found as the crossing points of the cumulants for different system sizes. Our results (see Fig. 4) show that we indeed observe such crossings, which is and indication for second order transition. Both the $`4`$-th order cumulant ratios for both $`O_N`$ and $`O_S`$ cross at the same temperature within the accuracy of our results, indicating that translational symmetry breaking (measured by $`O_S`$) and rotational symmetry breaking (measured by $`O_N`$) happen at the same or at very close temperatures.

The phase diagram of the doped system away from half filling shows an additional “supersolid” phase where the vacancies doped into the stripes form a superfluid smectic with broken translational and rotational symmetry as well as a finite superfluid density. Since the rotational symmetry is spontaneously broken in the smectic phase, the superfluid becomes anisotropic and Eq. (5) for the superfluid density needs to be modified. The geometric mean of the superfluid densities perpendicular and parallel to the stripes $`\rho _s=\sqrt{\rho _{s,||}\rho _{s,}}`$ is the quantity exhibiting the universal jump at the Kosterlitz-Thouless transition and was used in the determination of the phase transition. Care was taken that the anisotropy did not become too large to introduce systematic errors as pointed out by Prokof’ev and Svistunov rhos . The superfluid densities parallel ($`\rho _{s,||}`$) and perpendicular ($`\rho _{s,}`$) to the stripe order were measured for each configuration from the winding numbers in either $`x`$ or $`y`$ direction depending on whether $`S_n(0,\pi )`$ was larger or smaller than $`S_n(\pi ,0)`$.

Although superfluidity coexists with smectic order, it is strongly suppressed and $`T_c`$ and $`\rho _s`$ vanish linearly as half filling is approached. Using $`m_{||}=1/2t`$ and $`m_{}=\frac{1}{2t^2/(4V_2)}`$ as the boson masses parallel and perpendicular to the stripes, respectively, we obtain $`\rho _s=(0.78\pm 0.07)|\rho \frac{1}{2}|`$ consistent with the value obtained previously for a dilute gas of bosons bernardet .

In order to investigate the nature of the doping-driven melting of the smectic phase we repeat a similar procedure as at half filling. In contrast to the half filled case, with a clear jump in the values of $`O_S`$ and $`O_N`$ at the phase transition, or in the doped case of the nearest neighbor model with a jump in density $`\rho `$ at the doping-driven melting transition schmid , no such jump or double-peak structure was observed here. We checked the order parameters $`O_S`$, $`O_N`$, $`O_k`$, and the density $`\rho `$ in systems with size up to $`L=48`$. Repeating the simulations of Ref. batrouni3 , where a local update algorithm was used instead of our nonlocal algorithm we cannot reproduce the jump in the superfluid density seen there. Instead we find the smooth behavior shown in Fig. 5 for the same system, which is also in agreement with analytical considerations using mean-field and spin-wave analysis pich  or renormalization group calculations frey . Since our results are well converged and do not suffer from any systematic error, we believe the differences to be due to the method used to calculate the superfluid densities in Ref. batrouni3 .

At $`T=t/6`$ we checked whether rotational and translational symmetry breaking occur at the same critical density by comparing estimates obtained using Binder cumulant ratios for $`O_S`$ and $`O_N`$ as well as a determination of the critical point from the maxima of the susceptibility associated with $`O_k`$. The estimates agree within the error bars: $`\rho _c=0.2468\pm 0.0025`$ at $`T=t/6`$.

Finally we want to discuss our results in view of the three possibilities for the nature of the melting transitions of i) a first order transition, ii) a single second order transition or iii) two separate second order transitions. Close to half filling and at low temperatures we have very clear evidence for a single first order transition. The results for thermal melting at higher temperatures and for the doping-driven melting at low temperatures are not as clear-cut. While a very weak first order transition is always possible, it is more likely that we have second order transitions, especially in view of the fact that in the limit $`t/V_2\mathrm{}`$ there will be a second order transition in the Ising universality class. One way to distinguish between the scenarios (ii) of one single transition or (iii) of two second order transitions with an intervening nematic phase is that in the latter case both transitions should be in the Ising universality class since they each break a $`Z_2`$ symmetry. At the transition we would thus expect that the Binder cumulant ratios at $`T_c`$ takes on the value $`C_4=0.6106900(1)`$ kamieniarz . The fact that our results for the cumulant ratio of $`O_S`$ are close but not identical to this universal value for an Ising transition rather indicates secnario (ii), a single second order transition in a different universality class. In the doped system our data both for the cumulant ratios and the scaling of the susceptibility of $`O_k`$ are actually closer to the three dimensioonal $`XY`$ universality class predicted by Ref. frey  than the Ising universality class. More extensive simulations on larger systems, using newly developed flat histogram methods for quantum systems QWL  will be needed to more accurately determine the universality class of this phase transition.

While the exact order and universality class of the melting transition are hard to determine – which is not surprising given the difficulties known from the similar classical problem in the continuum – our results are interesting also with respect to the existence of a nematic phase. Such a phase must be restricted to a very narrow temperature and doping regime, smaller than the resolution of our simulations. It might be possible to stabilize a nematic phase with additional terms in the Hamiltonoan. One suggestion kivelson  is to add a term proportional to the square of the nermatic order parameter $`VO_k^2`$, to stabilize a nematic phase. This term gives two contributions: a nearest neighbor repulsion $`2V_{𝐢,𝐣}n_𝐢n_𝐣`$ and an additional next nearest neighbor hopping term $`V_{i,j}(a_𝐢^{}a_𝐣+a_𝐣^{}a_𝐢)`$. These frustrating terms, the second of which unfortunately causes a negative sign problem for quantum Monte Carlo simulations, might be important for the stability of an extended nematic phase. It is interesting to compare this model with the related frustrated Heisenberg model on a square lattice which however exhibits translational and rotational symmetry breaking only in the ground state at $`T=0`$ singh .

We are grateful for stimulating discussions with G. Batrouni, S. Kivelson and S. Todo and acknowledge support of the Swiss National Science Foundation. The simulations were performed on the Asgard Beowulf cluster at ETH Zürich, using a parallelizing C++ library for Monte Carlo simulations alea .


format: text

text:
# One-dimensional Bosons in Three-dimensional Traps

(August 29, 2003)

## Abstract

Recent experimental and theoretical work has indicated conditions in which a trapped, low-density Bose gas ought to behave like the 1D delta-function Bose gas solved by Lieb and Liniger. Up to now the theoretical arguments have been based on variational - perturbative ideas or numerical investigations. There are 4 parameters: density, transverse and longitudinal dimensions, and scattering length. In this paper we explicate 5 parameter regions in which various types of 1D or 3D behavior occur in the ground state. Our treatment is based on a rigorous analysis of the many-body Schrödinger equation.

This version of the article is a combination of the earlier 4-page arXiv version and the revised 4-page Physical Review Letters version.

It appears to be possible to do experiments in highly elongated traps on ultra-cold Bose gases that are effectively 1D. More precisely, the 1D Bose gas with delta function 2-body interaction, analyzed long ago LL ; L , should be visible, with its quasi-fermionic behavior gir , the absence of Bose-Einstein condensation (BEC) in a dilute limit Lenard ; pita ; girardeau , and an excitation spectrum different from that predicted by Bogoliubov’s theory L ; jackson ; komineas . Several theoretical investigations on the transitions from 3D to an effective 1D behavior were triggered by olshanii . See, e.g., dunjko ; petrov ; kolomeisky ; tanatar ; girardeau2 ; das ; das2 ; menotti . Systems showing indications of such a transition have recently been prepared experimentally bongs ; goerlitz ; greiner ; schreck .

The theoretical work on the dimensional cross-over for the ground state in elongated traps has so far been based either on variational calculations, starting from a 3D delta-potential olshanii ; das2 ; girardeau2 , or on numerical Quantum Monte Carlo studies blume ; astra  with more realistic, genuine 3D potentials, but particle numbers limited to the order of 100. This work is important and has led to valuable insights, in particular about different parameter regions petrov ; dunjko , but a more thorough theoretical understanding is clearly desirable since this is not a simple problem. In fact, it is evident that for a potential with a hard core the true 3D wave functions do not approximately factorize in the longitudinal and transverse variables (otherwise the energy would be infinite) and the effective 1D potential can not be obtained by simply integrating out the transverse variables of the 3D potential (that would immediately create an impenetrable barrier in 1D). It is important to be able to demonstrate rigorously, and therefore unambiguously, that the 1D behavior really follows from the fundamental Schrödinger equation. It is also important to delineate, as we do here, precisely what can be seen in the different parameter regions. The full proofs of our assertions are long and will be given elsewhere LSY , but we emphasize that everything can be rigorously derived from first principles. In this paper we state our main results and outline the basic ideas for the proofs.

We shall always be concerned with the ground state and with large particle number, $`N1`$, which is appropriate for the consideration of actual experiments. In order to make precise statements we take the limit $`N\mathrm{}`$ but the reader can confidently apply these limiting statement to finite numbers like $`N=100`$. Other parameters of the problem are the scattering length, $`a`$, of the two-body interaction potential, $`v`$, and two lengths, $`r`$ and $`L`$, describing the transverse and the longitudinal extension of the trap potential, respectively.

It is convenient to write the Hamiltonian in the following way (in units where $`\mathrm{}=2m=1`$):

$`H_{N,L,r,a}`$ $`={\displaystyle \underset{j=1}{\overset{N}{}}}\left(_j^2+V_r^{}(𝐱_j^{})+V_L(z_j)\right)`$
$`+{\displaystyle \underset{1i<jN}{}}v_a(|𝐱_i𝐱_j|)`$ (1)
with $`𝐱=(x,y,z)=(𝐱^{},z)`$ and with

$`V_r^{}(𝐱^{})`$ $`={\displaystyle \frac{1}{r^2}}V^{}(𝐱^{}/r),`$
$`V_L(z)={\displaystyle \frac{1}{L^2}}V(z/L)`$ $`,v_a(|𝐱|)={\displaystyle \frac{1}{a^2}}v(|𝐱|/a).`$ (2)
Here, $`r,L,a`$ are variable scaling parameters while $`V^{}`$, $`V`$ and $`v`$ are fixed. The interaction potential $`v`$ is supposed to be nonnegative, of finite range and have scattering length 1; the scaled potential $`v_a`$ then has scattering length $`a`$. The external trap potentials $`V`$ and $`V^{}`$ confine the motion in the longitudinal ($`z`$) and the transversal ($`𝐱^{}`$) directions, respectively, and are assumed to be continuous and tend to $`\mathrm{}`$ as $`|z|`$ and $`|𝐱^{}|`$ tend to $`\mathrm{}`$. To simplify the discussion we find it also convenient to assume that $`V`$ is homogeneous of some order $`s>0`$, namely $`V(z)=|z|^s`$, but weaker assumptions, e.g. asymptotic homogeneity lsy2 , would in fact suffice. The case of a simple box with hard walls is realized by taking $`s=\mathrm{}`$, while the usual harmonic approximation is $`s=2`$. It is understood that the lengths associated with the ground states of $`d^2/dz^2+V(z)`$ and $`(^{})^2+V^{}(𝐱^{})`$ are both of the order $`1`$ so that $`L`$ and $`r`$ measure, respectively, the longitudinal and the transverse extensions of the trap. We denote the ground state energy of (One-dimensional Bosons in Three-dimensional Traps) by $`E^{\mathrm{QM}}(N,L,r,a)`$ and the ground state particle density by $`\rho _{N,L,r,a}^{\mathrm{QM}}(𝐱)`$.

While the 3D density is always assumed to be low (in the sense that distance between particles is large compared to the 3D scattering length) the 1D density can be either high or low. In contrast to 3D gases, high density in 1D corresponds to weak interactions and vice versa LL .

To keep the discussion simple let us first think of the case that the particles are confined in a box with dimensions $`r`$ and $`L`$. The 3D particle density is then $`\rho ^{3\mathrm{D}}=N/(r^2L)`$ and the 1D density $`\rho ^{1\mathrm{D}}=N/L`$. The case of quadratic or more general homogeneous trapping potentials will be considered later. We begin by describing the division of the space of parameters into two basic regions. This decomposition will eventually be refined into five regions, but for the moment let us concentrate on the basic dichotomy.

In earlier work lsy1 ; lsy2  we proved that the 3D Gross-Pitaevskii formula for the energy (including its limiting ‘Thomas-Fermi’ case) is correct to leading order in situations in which $`a`$ is small and $`N`$ is large. This energy has two parts: The energy necessary to confine the particles in the trap, which is roughly $`N(r^2+L^2)`$, plus the internal energy of interaction, which is $`N4\pi a\rho ^{3\mathrm{D}}`$. The trouble is that while this formula is correct for a fixed confining potential in the limit $`N\mathrm{}`$ with $`a^3\rho ^{3\mathrm{D}}0`$, it does not hold uniformly if $`r/L`$ gets small as $`N`$ gets large. In other words, new physics can come into play as $`r/L0`$ and it turns out that this depends on the ratio of $`a/r^2`$ to $`\rho ^{1\mathrm{D}}=N/L`$ . There are two basic regimes to consider in highly elongated traps, i.e., when $`rL`$. They are are

* The 1D limit of the 3D Gross-Pitaevskii/‘Thomas-Fermi’ regime
* The ‘true’ 1D regime.

The former turns out to be characterized by $`aL/r^2N0`$, while in the latter regime $`aL/r^2N`$ is of the order one or even tends to infinity. (This is usually referred to as the Girardeau-Tonks region, but it was Girardeau gir , not Tonks, who understood how to calculate the states of the 1D hard-core gas.) These two situations correspond to high 1D density (weak interaction) and low 1D density (strong interaction), respectively. Physically, the main difference is that in the strong interaction regime the motion of the particles in the longitudinal direction is highly correlated, while in the weak interaction regime it is not. Mathematically, this distinction also shows up in our proofs.

In both regions the internal energy of the gas is small compared to the energy of confinement which is of the order $`N/r^2`$. However, this in itself does not imply a specifically 1D behavior. (If $`a`$ is sufficiently small it is satisfied in a trap of any shape.) 1D behavior, when it occurs, manifests itself by the fact that the transverse motion of the atoms is uncorrelated while the longitudinal motion is correlated (very roughly speaking) in the same way as pearls on a necklace. Thus, the true criterion for 1D behavior is that $`aL/r^2N`$ is of the order unity or larger, and not merely the condition that the energy of confinement dominates the internal energy.

In parallel with the 3D Hamiltonian we consider the Hamiltonian for $`n`$ Bosons in 1D with delta interaction and coupling constant $`g0`$ , i.e.,

$$H_{n,g}^{1\mathrm{D}}=\underset{j=1}{\overset{n}{}}^2/z_j^2+g\underset{1i<jn}{}\delta (z_iz_j).$$
(3)
We consider this Hamiltonian for the $`z_j`$ in an interval of length $`\mathrm{}`$ in the thermodynamic limit, $`\mathrm{}\mathrm{}`$, $`n\mathrm{}`$ with $`\rho =n/\mathrm{}`$ fixed. The ground state energy per particle in this limit is independent of boundary conditions and can, according to LL , be written as

$$e_0^{1\mathrm{D}}(\rho )=\rho ^2e(g/\rho ),$$
(4)
with a function $`e(t)`$ determined by a certain integral equation. Its asymptotic form is $`e(t)\frac{1}{2}t`$ for $`t1`$ and $`e(t)\pi ^2/3`$ for $`t\mathrm{}`$. Thus

$$e_0^{1\mathrm{D}}(\rho )\frac{1}{2}g\rho \text{for}g/\rho 1$$
(5)
and

$$e_0^{1\mathrm{D}}(\rho )(\pi ^2/3)\rho ^2\text{for}g/\rho 1.$$
(6)
Taking $`\rho e_0^{1\mathrm{D}}(\rho )`$ as a local energy density for an inhomogeneous 1D system we can form the energy functional

$$[\rho ]=_{\mathrm{}}^{\mathrm{}}\left(|\sqrt{\rho }(z)|^2+V_L(z)\rho (z)+\rho (z)^3e(g/\rho (z))\right)𝑑z$$
(7)
with ground state energy defined to be

$$E^{1\mathrm{D}}(N,L,g)=inf\{[\rho ]:\rho (z)0,_{\mathrm{}}^{\mathrm{}}\rho (z)𝑑z=N\}.$$
(8)
By standard methods (cf., e.g., lsy1 ) one can show that there is a unique minimizer, i.e., a density $`\rho _{N,L,g}(z)`$ with $`\rho _{N,L,g}(z)𝑑z=N`$ and $`[\rho _{N,L,g}]=E^{1\mathrm{D}}(N,L,g)`$. We define the mean 1D density of this minimizer to be

$$\overline{\rho }=\frac{1}{N}_{\mathrm{}}^{\mathrm{}}\left(\rho _{N,L,g}(z)\right)^2𝑑z.$$
(9)
In a rigid box, i.e., for $`s=\mathrm{}`$, $`\overline{\rho }`$ is simply $`N/L`$ (except for boundary corrections), but in more general traps it depends also on $`g`$ besides $`N`$ and $`L`$. The order of magnitude of $`\overline{\rho }`$ in the various parameter regions will be described below.

Our main result relates the 3D ground state energy of (One-dimensional Bosons in Three-dimensional Traps), $`E^{\mathrm{QM}}(N,L,r,a)`$, to the 1D density functional energy $`E^{1\mathrm{D}}(N,L,g)`$ in the large $`N`$ limit with $`ga/r^2`$ provided $`r/L`$ and $`a/r`$ are sufficiently small. To state this precisely, let $`e^{}`$ and $`b(𝐱^{})`$, respectively, denote the ground state energy and the normalized ground state wave function of $`(^{})^2+V^{}(𝐱^{})`$. The corresponding quantities for $`(^{})^2+V_r^{}(𝐱^{})`$ are $`e^{}/r^2`$ and $`b_r(𝐱^{})=(1/r)b(𝐱^{}/r)`$. In the case that the trap is a cylinder with hard walls $`b`$ is a Bessel function; for a quadratic $`V^{}`$ it is a Gaussian.

Define $`g`$ by

$$g=\frac{8\pi a}{r^2}|b(𝐱^{})|^4d^2𝐱^{}=8\pi a|b_r(𝐱^{})|^4d^2𝐱^{}.$$
(10)
###### Theorem.

Let $`N\mathrm{}`$ and simultaneously $`r/L0`$ and $`a/r0`$ in such a way that $`r^2\overline{\rho }\mathrm{min}\{\overline{\rho },g\}0`$. Then

$$lim\frac{E^{\mathrm{QM}}(N,L,r,a)Ne^{}/r^2}{E^{1\mathrm{D}}(N,L,g)}=1.$$
(11)
Moreover, if we define the 1D quantum mechanical density by averaging over the transverse variables, i.e.,

$$\widehat{\rho }_{N,L,r,a}^{\mathrm{QM}}(z)\rho _{N,L,r,a}^{\mathrm{QM}}(𝐱^{},z)d^2𝐱^{},$$
(12)
then $`\widehat{\rho }_{N,L,r,a}^{\mathrm{QM}}(z)/\rho _{N,L,g}(z)1`$ in a suitable sense.

Note that because of (5) and (6) the condition $`r^2\overline{\rho }\mathrm{min}\{\overline{\rho },g\}0`$ is the same as

$$e_0^{1\mathrm{D}}(\overline{\rho })1/r^2,$$
(13)
i.e., the average energy per particle associated with the longitudinal motion should be much smaller than the energy gap between the ground and first excited state of the confining Hamiltonian in the transverse directions. Thus, the basic physics is highly quantum-mechanical and has no classical counterpart. The system can be described by a 1D functional (7), even though the transverse trap dimension is much larger than the range of the atomic forces.

The two regimes mentioned previously correspond to specific restrictions on the size of the ratio $`g/\overline{\rho }`$ as $`N\mathrm{}`$, namely $`g/\overline{\rho }1`$ for the limit of the 3D Gross-Pitaevskii regime (weak interaction/high density), and $`g/\overline{\rho }>0`$ for the ‘true’ 1D regime (strong interaction/low density). The precise meaning of $``$ is that the ratio of the left side to the right side tends to zero in the limit considered.

We shall now briefly describe the finer division of these two regimes into five regions altogether. Three of them (Regions 1–3) belong to the weak interaction regime and two (Regions 4–5) to the strong interaction regime. In each of these regions the general functional (7) can be replaced by a different, simpler functional, and the energy $`E^{1\mathrm{D}}(N,L,g)`$ in the theorem by the ground state energy of that functional, in analogy with the fact that in 3D the TF functional, i.e., the GP functional without gradient term, is a limit of the full GP functional. We state the results for the ground state energy, but corresponding results hold for the density. We always assume $`N\mathrm{}`$ and $`r/L0`$.

The five regions are

$``$ Region 1, the Ideal Gas case: $`g/\overline{\rho }N^2`$, with $`\overline{\rho }N/L`$, corresponding to the trivial case where the interaction is so weak that it effectively vanishes in the large $`N`$ limit and everything collapses to the ground state of $`d^2/dz^2+V(z)`$ with ground state energy $`e^{}`$. The energy $`E^{1\mathrm{D}}`$ in (11) can be replaced by $`Ne^{}/L^2`$. Note that $`g/\overline{\rho }N^2`$ means that the 3D interaction energy $`\rho ^{3\mathrm{D}}a1/L^2`$.

$``$ Region 2, the 1D GP case: $`g/\overline{\rho }N^2`$, with $`\overline{\rho }N/L`$, described by a 1D Gross-Pitaevskii energy functional

$$^{\mathrm{GP}}[\rho ]=_{\mathrm{}}^{\mathrm{}}\left(|\sqrt{\rho }(z)|^2+V_L(z)\rho (z)+\frac{1}{2}g\rho (z)^2\right)𝑑z$$
(14)
corresponding to the high density approximation (5) of the interaction energy in (7). Its ground state energy is $`E^{\mathrm{GP}}(N,L,g)=NL^2E^{\mathrm{GP}}(1,1,NgL)`$, by scaling.

$``$ Region 3, the 1D TF case: $`N^2g/\overline{\rho }1`$, with $`\overline{\rho }(N/L)(NgL)^{1/(s+1)}`$, where $`s`$ is the degree of homogeneity of the longitudinal confining potential $`V`$. This region is described by a Thomas-Fermi type functional

$$^{\mathrm{TF}}[\rho ]=_{\mathrm{}}^{\mathrm{}}\left(V_L(z)\rho (z)+\frac{1}{2}g\rho (z)^2\right)𝑑z.$$
(15)
It is a limiting case of Region 2 in the sense that $`NgLNaL/r^2\mathrm{}`$, but $`a/r`$ is sufficiently small so that $`g/\overline{\rho }(aL/Nr^2)(NaL/r^2)^{1/(s+1)}0`$, i.e., the high density approximation in (5) is still valid. In this limit the gradient term in (14) becomes vanishingly small compared to the other terms. The explanation of the factor $`(NaL/r^2)^{1/(s+1)}(NgL)^{1/(s+1)}`$ is as follows: The linear extension $`\overline{L}`$ of the minimizing density $`\rho _{N,L,g}^{\mathrm{GP}}`$ is for large values of $`NgL`$ determined by $`V_L(\overline{L})g(N/\overline{L})`$, which gives $`\overline{L}(NgL)^{1/(s+1)}L`$. In addition condition (13) requires $`g\overline{\rho }r^2`$, which means that $`Na/L(NaL/r^2)^{1/(s+1)}0`$. The minimum energy of (15) has the scaling property $`E^{\mathrm{TF}}(N,L,g)=NL^2(NgL)^{s/(s+1)}E^{\mathrm{TF}}(1,1,1)`$.

$``$ Region 4, the LL case: $`g/\overline{\rho }1`$ , with $`\overline{\rho }(N/L)N^{2/(s+2)}`$, described by an energy functional

$$^{\mathrm{LL}}[\rho ]=_{\mathrm{}}^{\mathrm{}}\left(V_L(z)\rho (z)+\rho (z)^3e(g/\rho (z))\right)𝑑z.$$
(16)
This region corresponds to the case $`g/\overline{\rho }1`$, so that neither the high density (5) nor the low density approximation (6) is valid and the full LL energy (4) has to be used. The extension $`\overline{L}`$ of the system is now determined by $`V_L(\overline{L})(N/\overline{L})^2`$ which leads to $`\overline{L}_{\mathrm{LL}}=LN^{2/(s+2)}`$, in contrast to $`\overline{L}_{\mathrm{TF}}=L(NgL)^{1/(s+1)}`$ in the TF case. Condition (13) means in this region that $`Nr/\overline{L}_{\mathrm{LL}}N^{s/(s+2)}r/L0`$. Since $`Nr/\overline{L}_{\mathrm{LL}}(\overline{\rho }/g)(a/r)`$, this condition is automatically fulfilled if $`g/\overline{\rho }`$ is bounded away from zero and $`a/r0`$. The ground state energy of (16), $`E^{\mathrm{LL}}(N,L,g)`$, is equal to $`N\gamma ^2E^{\mathrm{LL}}(1,1,g/\gamma )`$, where we introduced the density parameter

$$\gamma N/\overline{L}_{\mathrm{LL}}=(N/L)N^{2/(s+2)}.$$
(17)
$``$ Region 5, the GT case: $`g/\overline{\rho }1`$, with $`\overline{\rho }(N/L)N^{2/(s+2)}`$, described by a functional with energy density $`\rho ^3`$, corresponding to the Girardaeu-Tonks limit of the the LL energy density. It corresponds to impenetrable particles, i.e, the limiting case $`g/\overline{\rho }\mathrm{}`$ and hence formula (6) for the energy density. As in Region 4, the mean density is here $`\overline{\rho }\gamma `$. The energy functional is

$$^{\mathrm{GT}}[\rho ]=_{\mathrm{}}^{\mathrm{}}\left(V_L(z)\rho (z)+(\pi ^2/3)\rho (z)^3\right)𝑑z,$$
(18)
with minimum energy $`E^{\mathrm{GT}}(N,L)=N\gamma ^2E^{\mathrm{GT}}(1,1)`$.

We note that the condition $`g/\overline{\rho }1`$ means that Region 4 requires the gas cloud to have aspect ratio $`r/\overline{L}`$ of the order $`N^1(a/r)`$ or smaller, where $`\overline{L}LN^{2/(s+2)}`$ is the length of the cloud. Experimentally, such small aspect ratios are quite a challenge and the situations described in bongs ; goerlitz ; greiner ; schreck  are still rather far from this regime. It may not be completely out of reach, however.

Regions 1–3 can be reached as limiting cases of a 3D Gross-Pitaevskii theory. In this sense, the behavior in these regions contains remnants of the 3D theory, which also shows up in the the fact that Bose-Einstein condensation (BEC) prevails in Regions 1 and 2 LSY . Heuristically, these traces of 3D can be understood from the fact that in Regions 1–3 the 1D formula for energy per particle, $`g\rho aN/(r^2L)`$, gives the same result as the 3D formula LY1998 , i.e., scattering length times 3D density. This is no longer so in Regions 4 and 5 and different methods are required.

BEC probably occurs in part of Region 3, besides Regions 1 and 2, but we cannot prove this and it remains an open problem. BEC means LS  that the one-body density matrix $`\gamma (𝐱,𝐱^{})`$ factorizes as $`N\varphi (𝐱)\varphi (𝐱^{})`$ for some normalized $`\varphi `$. This, in fact, is 100% condensation and this is what we prove occurs (in the $`N\mathrm{}`$ limit, of course). The function $`\varphi `$ is the square-root of the minimizer of the 1D GP functional (14) times the transverse function $`b_r(𝐱^{})`$. The proof is similar to the work in LS . BEC is not expected in Regions 4 and 5. Lenard Lenard  showed that the largest eigenvalue of $`\gamma `$ grows only as $`N^{1/2}`$ for a homogeneous gas of 1D impenetrable bosons and, according to papenbrock  and forrester , this holds also for a GT gas in a harmonic trap. (The exponent 0.59 in girardeau  can probably be ascribed to the small number of particles ($`N=10`$) considered.)

We now outline the main steps in the proof of the Theorem, referring to LSY  for full details. The different parameter regions have to be treated by different methods, a watershed lying between Regions 1–3 on the one hand and Regions 4–5 on the other. In Regions 1–3, similar methods as in the proof of the 3D Gross-Pitaevskii limit theorem in lsy1  can be used. This 3D proof needs considerable modifications, however, because in lsy1  the external potential is fixed and the estimates are not uniform in the ratio $`r/L`$.

To prove (11) one has to establish upper and lower bounds, with controlled errors, on the QM many-body energy in terms of the energies obtained by minimizing the energy functionals appropriate for the various regions. The limit theorem for the densities can be derived from the energy estimates in a standard way by variation with respect to the external potential $`V_L`$. As usual, the upper bounds for the energy are easier than the lower bounds, but nevertheless not simple, in particular not for ‘hard’ potentials $`v`$.

The upper bound in Regions 1–3 is obtained from a variational ansatz of the form $`\mathrm{\Psi }(𝐱_1,\mathrm{},𝐱_N)=F(𝐱_1,\mathrm{},𝐱_N)_{k=1}^Nb_r(𝐱_k^{})\sqrt{\rho ^{\mathrm{GP}}}(z_k),`$ with $`F(𝐱_1,\mathrm{},𝐱_N)`$ $`=_{k=1}^Nf(𝐱_k𝐱_{j(k)})`$ where $`𝐱_{j(k)}`$ is the nearest neighbor of $`𝐱_k`$ among the points $`𝐱_j`$, $`j<k`$, $`b_r(𝐱^{})`$ is the lowest eigenfunction of $`(^{})^2+V_r^{}`$, and $`\rho ^{\mathrm{GP}}(z)`$ the minimizer of the 1D GP functional (14). The function $`f`$ is, up to a cut-off length that has to be chosen optimally, the zero energy scattering solution for the two-body Hamiltonian with interaction $`v_a`$. The form of $`F`$, inspired by dyson , is chosen rather than a Jastrow ansatz $`_{i<j}f(𝐱_i𝐱_j)`$ because it is computationally simpler for the purpose of obtaining rigorous estimates.

For an upper bound in Regions 4–5 a natural variational ansatz would appear to be $`\mathrm{\Psi }(𝐱_1,\mathrm{},𝐱_N)=F(𝐱_1,\mathrm{},𝐱_N)_{k=1}^Nb_r(𝐱_k^{})\psi (z_1,\mathrm{},z_N)`$ where $`\psi `$ is the ground state of $`H_{N,g}^{1\mathrm{D}}`$ with the external potential $`V_L`$ added. However, in order to make a link with the exact solution (4) for a homogeneous gas, but also to control the norm of the trial function, it turns out to be necessary to localize the particles by dividing the trap into finite ‘boxes’ (finite in $`z`$-direction), with a finite particle number in each box and making the ansatz with the boundary condition $`\mathrm{\Psi }=0`$ for each box individually. The particles are then distributed optimally among the boxes to minimize the energy. This box method, but with the boundary condition $`\mathrm{\Psi }=0`$, is also used for the lower bounds to the energy. Another essential device for the lower bounds is Dyson’s Lemma that was also used in LS ; lsy1 ; LY1998 . This lemma, which goes back to Dyson’s seminal paper dyson  on the hard core Bose gas, estimates the kinetic and potential energy for a Hamiltonian with a ‘hard’ potential $`v`$ of finite range from below by the potential energy of a ‘soft’ potential $`U`$ of larger range but essentially the same scattering length as $`v`$. Borrowing a tiny part of the kinetic energy it is then possible to do perturbation theory with the ‘soft’ potential $`U`$ and use Temple’s inequality TE  to bound the errors. A direct application of perturbation theory to the original potential $`v`$, on the other hand, is in general not possible.

A core lemma for Regions 4–5 is a lower bound on the 3D ground state energy in a finite box in terms of the 1D energy of the Hamiltonian (3), both with the boundary condition $`\mathrm{\Psi }=0`$. Denoting the former energy by $`E_{\mathrm{box}}^{3\mathrm{D}}`$ and the latter by $`E_{\mathrm{box}}^{1\mathrm{D}}`$, this bound for $`n`$ particles in a box of length $`\mathrm{}`$ in the $`z`$-direction reads

$$E_{\mathrm{box}}^{3\mathrm{D}}\frac{ne^{}}{r^2}E_{\mathrm{box}}^{1\mathrm{D}}\left(1Cn\left(\frac{a}{r}\right)^{1/8}\left[1+\frac{nr}{\mathrm{}}\left(\frac{a}{r}\right)^{1/8}\right]\right)$$

with a constant $`C`$. To prove this bound the ground state wave function is first written as a product of $`_kb_r(𝐱_k^{})`$ and a function $`G(𝐱_1,\mathrm{},𝐱_n)`$. This subtracts $`ne^{}/r^2`$ from $`E_{\mathrm{box}}^{3\mathrm{D}}`$ but the resulting minimization problem for $`G`$ involves the weighted measure $`_kb_r(𝐱_k^{})^2d^3𝐱_k`$ in place of $`_kd^3𝐱_k`$. Nevertheless, Dyson’s Lemma can be used, and after the ‘hard’ potential $`v`$ has been replaced by the ‘soft’ potential $`U`$, it is possible to integrate the transverse variables away and obtain a minimization problem for a 1D many body Hamiltonian with interaction $`d(z_iz_j)=ab_r(𝐱_i^{})^2b_r(𝐱_j^{})^2U(𝐱_i𝐱_j)d^2𝐱_i^{}d^2𝐱_j^{}`$. In the limit considered this converges to a delta interaction with the coupling constant (10). The error terms in the estimate for $`E_{\mathrm{box}}^{3\mathrm{D}}`$ arise both from Temple’s inequality and the replacement of $`d`$ by a delta function, among other things. When the particles are distributed optimally among the boxes to obtain a global lower bound, superadditivity of the energy and convexity of the energy density $`\rho ^3e(g/\rho )`$ are used, generalizing corresponding arguments in LY1998 .

In conclusion, we have reported a rigorous analysis of the parameter regions in which a Bose gas in an elongated trap may or may not be expected to display 1D behavior in its ground state. This takes the form of theorems about the ground state energy and density. We always consider $`N`$ to be large and $`r/L`$ to be small, but this is not enough to distinguish 1D from 3D behavior. The distinction occurs when, additionally, the 1D density is not too large, e.g., $`aL/r^2N^{3/2}`$ in a harmonic trap. The 1D behavior, when it occurs, is described by the delta-function gas solved in LL . We also present a 1D energy functional, analogous to the Gross-Pitaevskii functional, that correctly describes the energy and density in all the 5 parameter regions considered here.


text:
# Single-spin physics with weak bosons at RHIC

## I Introduction

The measurement of longitudinal spin asymmetries in production of $`W^\pm `$ bosons at the Relativistic Heavy Ion Collider (RHIC) will provide an essential probe for spin-dependent quark distributions at high momentum scales $`Q^2`$ Bunce et al. (2000). At a $`pp`$ center-of-mass energy $`\sqrt{s}=500`$ GeV, about $`1.3\times 10^6`$ $`W^+`$and $`W^{}`$ bosons will be produced by the time the integrated luminosity reaches $`800\text{ pb}^1`$. Due to the parity violation in the $`q\overline{q}W`$ coupling, this process permits non-vanishing single-spin asymmetries $`A_L(\xi )`$, defined for any kinematical variable $`\xi `$ as

$$A_L(\xi )\frac{\frac{d\sigma (p^{}pWX)}{d\xi }\frac{d\sigma (p^{}pWX)}{d\xi }}{\frac{d\sigma (p^{}pWX)}{d\xi }+\frac{d\sigma (p^{}pWX)}{d\xi }}.$$
(1)
The Born-level expression for the asymmetry $`A_L(y_W)`$ with respect to the rapidity $`y_W`$ of the $`W`$ boson is particularly simple if the absolute value of $`y_W`$ is large. In that case, the Born-level $`A_L(y_W)`$ reduces to the ratio $`\mathrm{\Delta }q(x)/q(x)`$ of the polarized and unpolarized parton distribution functions Bourrely and Soffer (1994, 1995); Nadolsky (1995). Furthermore, $`A_L(y_W)`$ tests the flavor dependence of quark polarizations.

The original method for extracting the spin-dependent quark distributions out of the RHIC $`W`$ boson data is based on the direct reconstruction of the asymmetry $`A_L(y_W)`$ Bunce et al. (2000); Bland (2000). Unfortunately, such reconstruction is obstructed by specifics of the detection of $`W^\pm `$ bosons at RHIC. First, neither Phenix Morrison et al. (1998) nor Star Harris (1994) detector at RHIC is hermetic. Therefore, it is not possible at those detectors to monitor the energy balance in particle reactions, even in the transverse direction with respect to the beams. Due to the lack of information about the missing energy carried by the neutrino, the determination of $`y_W`$ is in general ambiguous and depends on the assumptions about the dynamics of the process. Second, due to the correlation between the spins of the initial-state quarks and final-state leptons, the measured value of $`A_L(y_W)`$ is strongly sensitive to the experimental cuts imposed on the observed charged lepton. Hence, unless a theory calculation exists to reliably predict the distributions of the leptons from the $`W`$ boson decay, the spin-dependent quark distribution functions cannot be determined with acceptable accuracy.

In weak boson production at the Fermilab Tevatron $`p\overline{p}`$ collider, the mass $`M_W`$ of the $`W`$ boson has been measured to a great precision Affolder et al. (2001); Abazov et al. (2002). The information on $`M_W`$ is extracted out of the lepton-level distributions $`d\sigma /m_T`$, $`d\sigma /dp_{Te}`$, or $`d\sigma /dp_{T\nu }`$, where $`p_{Te}|𝐩_{Te}|`$ and $`p_{T\nu }|𝐩_{T\nu }|`$ are the transverse momenta of the electron and neutrino, and $`m_T\sqrt{2p_{Te}p_{T\nu }2𝐩_{Te}𝐩_{T\nu }}`$ is the transverse mass of the lepton pair. The value of $`M_W`$ can be determined from the analysis of the position and shape of the Jacobian peaks in each of these three distributions. For most of the events contributing near the Jacobian peak, the transverse momentum $`q_T`$ of the $`W`$ boson is very small, of order of a few GeV. Therefore, at each order of the perturbative QCD calculation, there are large logarithmic contributions of the form $`\alpha _S^nq_T^2\mathrm{ln}^m(M_W^2/q_T^2)`$, where $`n=1,2,\mathrm{}`$ and $`m=0,1,\mathrm{},2n1`$. These large logarithmic contributions can be summed to all orders in $`\alpha _S`$ by applying the transverse momentum resummation formalism formulated in the impact parameter space by Collins, Soper, and Sterman (CSS) Collins et al. (1985). Without such a resummation calculation performed at the decay lepton level Balazs and Yuan (1997), it would not be possible to determine $`M_W`$ from the experimental data. This is because at hadron colliders the longitudinal momentum of the neutrino from the $`W`$ boson decay is not observed by the detectors; hence, the invariant mass of the $`W`$ boson cannot be measured directly. Instead, the determination of the $`W`$ boson mass relies on the analysis of the shape of the Jacobian peaks in the lepton-level distributions sensitive to $`q_T`$. Again, without a resummation calculation at the secondary lepton level, it would not be possible to predict the shape of these distributions accurately enough to determine $`M_W`$. We refer the reader to Ref. Balazs and Yuan (1997) for a detailed discussion of the phenomenology of the $`W`$ boson physics at the Tevatron collider.

Similarly, at RHIC, a resummation calculation is needed to reliably predict the distribution of the secondary leptons from the $`W`$ boson decay. As we already mentioned, the reason is the necessity to deduce spin-dependent parton luminosities based on the observation of the secondary charged lepton only. For RHIC purposes, the resummation calculation at the lepton level has to be modified to include the spin dependence, which has been performed in Ref. Nadolsky and Yuan (2003). In this paper, we will apply the calculation in Ref. Nadolsky and Yuan (2003) to model single-spin asymmetries of $`W`$ boson production in realistic RHIC conditions. We put the emphasis on the discussion of the asymmetries $`A_L(y_{\mathrm{}})`$ and $`A_L(p_T\mathrm{})`$ for the distributions in the rapidity $`y_{\mathrm{}}`$ and transverse momentum $`p_T\mathrm{}`$ of the observed charged lepton. We will argue that these directly observed asymmetries provide a viable alternative to the commonly discussed asymmetry $`A_L(y_W)`$.

The paper is organized as follows. In Section II, we briefly summarize the resummation formalism used for this study. A next-to-leading order (NLO) calculation is also presented for comparison to the resummation calculation. Since good understanding of spin-averaged cross sections is needed to derive polarized PDFs from the spin asymmetries, Section III discusses the uncertainties due to the imperfect knowledge of unpolarized parton distributions, as well as the potential for RHIC experiments to reduce these uncertainties. Section IV demonstrates that a correct model of multiple parton radiation is important to describe the differential distributions of the final-state leptons. We then employ the resummation calculation to realistically evaluate the feasibility of the reconstruction of the asymmetry $`A_L(y_W)`$ from the observed charged lepton data. Section V shows that the spin asymmetries for the distributions of the observed secondary leptons provide an attractive alternative to the asymmetry $`A_L(y_W)`$: they can be measured directly and discriminate efficiently between the different PDF sets. We also show that the RHIC experiments can explore the spin dependence of the nonperturbative contributions from the multiple parton radiation. The main findings of this paper are recapped in the conclusion.

## II  Theoretical essentials

The primary goal of this paper is to discuss production of $`W^\pm `$ bosons, which will be used at RHIC to obtain information about the longitudinally polarized parton distribution functions (PDFs) in the quark sector. Due to the unique feature of the maximal parity violation in the $`Wq\overline{q}`$ coupling, the parton-level cross sections for spin-dependent $`W`$ boson production have pronounced single-spin asymmetries. This feature gives $`W`$ boson production an edge on parity-conserving processes, in which non-trivial hadronic dynamics can be probed only through the measurement of more complex double-spin asymmetries.

Denote the unpolarized distribution of a parton $`a`$ in a nucleon $`A`$ as $`f_{a/A}(\xi ,\mu _F)`$ and the longitudinally polarized distribution of $`a`$ in the nucleon $`A`$ as $`\mathrm{\Delta }f_{a/A}(\xi ,\mu _F)`$. Then,

$`f_{a/A}(\xi ,\mu _F)`$ $``$ $`f_{+/+}(\xi ,\mu _F)+f_{/+}(\xi ,\mu _F),`$ (2)
$`\mathrm{\Delta }f_{a/A}(\xi ,\mu _F)`$ $``$ $`f_{+/+}(\xi ,\mu _F)f_{/+}(\xi ,\mu _F),`$ (3)
where $`f_{h_a/h_A}(\xi _a,\mu _F)`$ is a helicity-dependent parton distribution function, *i.e.*, a probability of finding a parton $`a`$ with the momentum $`p_a^\mu =\xi _ap_A^\mu `$ and helicity $`h_a`$ in a hadron $`A`$ with the momentum $`p_A^\mu `$ and helicity $`h_A.`$ The PDFs depend on the factorization scale $`\mu _F`$, which in our calculation is assumed to coincide with the QCD renormalization scale. To find the polarized PDFs from an experiment with one longitudinally polarized beam, one measures the unpolarized cross section,

$$\frac{d\sigma }{d\xi }=\frac{1}{2}\left(\frac{d\sigma (p^{}pWX)}{d\xi }+\frac{d\sigma (p^{}pWX)}{d\xi }\right),$$
(4)
as well as the single-spin cross section,

$$\frac{d\mathrm{\Delta }_L\sigma }{d\xi }=\frac{1}{2}\left(\frac{d\sigma (p^{}pWX)}{d\xi }\frac{d\sigma (p^{}pWX)}{d\xi }\right).$$
(5)
Here $`\xi `$ denotes any kinematic variable, such as the rapidity $`y_W`$ of the $`W`$ boson. The single-spin asymmetry is defined as

$$A_L(\xi )\frac{d\mathrm{\Delta }_L\sigma /d\xi }{d\sigma /d\xi }.$$
(6)
In Ref. Nadolsky and Yuan (2003), we presented fully differential unpolarized cross sections $`d\sigma /d\xi `$, single-spin cross sections $`d\mathrm{\Delta }_L\sigma /d\xi ,`$ and double-spin cross sections $`d\mathrm{\Delta }_{LL}\sigma /d\xi `$ for production and decay of $`\gamma ^{},`$$`W^\pm ,`$ and $`Z^0`$ bosons. To reliably predict the rate at any point of the momentum phase space, these cross sections include the all-order sum of large soft and collinear logarithms, which dominate at small transverse momentum $`q_T`$ of the vector boson. This all-order sum is combined, without double counting, with the finite-order $`𝒪(\alpha _S)`$ cross section, which dominates at large $`q_T`$. The single-spin resummed cross sections in the narrow width approximation were first obtained in Ref. Weber (1993). The main result of Ref. Nadolsky and Yuan (2003) is the (more involved) derivation of the resummed cross sections that also accounts for the decay of the vector bosons and angular distributions of the final-state leptons. Our result was derived in the $`\overline{MS}`$ factorization scheme Mertig and van Neerven (1996); Vogelsang (1996a, b); therefore, it is fully compatible with the $`\overline{MS}`$ parton distributions.

The resummed cross sections from Ref. Nadolsky and Yuan (2003) are implemented in the numerical resummation program Legacy and Monte-Carlo integration program RhicBos,111The Fortran code and input grids for RhicBos can be downloaded from http://hep.pa.msu.edu/~nadolsky/RhicBos/.  which were used to produce the results discussed in this paper. Legacy is a C++/Fortran program that quickly and accurately generates the resummed cross sections on a grid of points in the kinematical phase space. RhicBos reads in the cross section grids, performs their multi-dimensional integration, and produces the output that is optimized for RHIC specifics. Since the finite-order cross section is characterized by a smaller theoretical uncertainty at $`q_TQ`$, RhicBos switches from the resummed cross section ($`W+Y`$) to the finite-order cross section at the point where the resummation cross section becomes smaller than the finite-order cross section at $`q_T`$ above the position of the maximum in $`d\sigma /dq_T`$. To compare the resummation calculation to the finite-order approach, we have also calculated the finite-order $`𝒪(\alpha _S)`$ cross section (next-to-leading order, or NLO, cross section) with the help of a phase space slicing method. This method introduces a separation transverse momentum scale $`q_T^{sep}`$, below which the integral of $`d\sigma /dq_T^2`$ over the region $`0q_T^2(q_T^{sep})^2`$ is calculated analytically, using the small-$`q_T`$ approximation in Eq. (78) of Ref. Nadolsky and Yuan (2003). The integral over the region $`q_T^2(q_T^{sep})^2`$ is integrated numerically using the exact $`𝒪(\alpha _S)`$ cross section. As we will show later, in the NLO approach, predictions for several important distributions depend strongly on the value of the auxiliary parameter $`q_T^{sep}`$, so that the resummed cross sections must be used in order to model those distributions. Some other distributions, such as $`d\sigma /dQ^2`$, are not sensitive to $`q_T^{sep}`$ (up to higher-order corrections), so that the resummation and NLO calculations give close results for such distributions.

The electroweak parameters in this paper were obtained in the on-shell scheme according to the procedure described in Section II.B of Ref. Nadolsky and Yuan (2003) for the Fermi constant $`G_F=1.16639\times 10^5\text{ GeV}^2`$, $`W`$ boson mass $`M_W=80.419`$ GeV, and $`Z^0`$ boson mass $`M_Z=91.187`$ GeV. The width $`\mathrm{\Gamma }_W`$ of the $`W`$ boson was evaluated as $`\mathrm{\Gamma }_W=3G_FM_W^3/(2\pi \sqrt{2})`$. The finite-order and asymptotic cross sections were calculated for the QCD factorization scale $`\mu _F=Q`$. The resummed cross section (see Eqs. (79) and (80) in Ref. Nadolsky and Yuan (2003)) was calculated for the canonical values of the scale parameters, $`C_1=C_3=2e^{\gamma _E}`$, where $`\gamma _E0.577`$ is the Euler constant, and $`C_2=1`$. Unless stated otherwise, the nonperturbative Sudakov factor in the resummed cross sections was evaluated by using the recent Gaussian parametrization Landry et al. (2002) derived from the global analysis of transverse momentum distributions in unpolarized vector boson production. The projected statistical errors $`\delta A_L`$ in the measurement of the single-spin asymmetries were estimated according to Eq. (13) in Ref. Bunce et al. (2000):

$$\delta A_L=\sqrt{\frac{1}{NP^2}\frac{1}{N}A_L^2},$$
(7)
where $`N=\sigma `$ is the number of the events for the given luminosity $``$ and cross section $`\sigma `$. Unless stated otherwise, the assumed integrated luminosity is $`=800\text{ pb}^1`$, and polarization $`P`$ of the beam $`A`$ is $`70\%`$.

## III  Physics potential in unpolarized $`pp`$ collisions

To extract the single-spin cross section $`d\mathrm{\Delta }_L\sigma /d\xi `$ from the single-spin asymmetry $`A_L(\xi )`$, the spin-averaged cross sections $`d\sigma /d\xi `$ must be known well. Obviously, the uncertainty in the knowledge of $`d\sigma /d\xi `$ must be small enough as compared to the targeted uncertainty in the measurement of $`d\mathrm{\Delta }_L\sigma /d\xi `$. Therefore, an interesting question is: how well do we currently know the unpolarized cross sections at RHIC energies, and how important is their measurement in the upcoming RHIC experiments? The goal of this section is to argue that at present the uncertainties in the unpolarized cross sections at RHIC are sizeable; to reliably predict the unpolarized rate, these uncertainties must be reduced by dedicated RHIC measurements.

Table 1 lists the unpolarized cross sections and event rates for the processes $`pp(W^\pm \mathrm{}\nu _{\mathrm{}})X`$ and $`pp(Z^0\mathrm{}\overline{\mathrm{}})X`$ at RHIC for various values of center-of-mass energies $`\sqrt{s}`$ and integrated luminosities $``$. The cross sections were derived at NLO using the updated CTEQ6 parton distribution functions Pumplin et al. (2002); Stump et al. (2003). The cross sections $`\sigma `$ are calculated using the best-fit PDF parametrization CTEQ6M. The uncertainty $`\delta \sigma _{PDF}`$ in the cross section $`\sigma `$ is generated by the uncertainty in the unpolarized parton luminosities, which are themselves phenomenological functions known with a finite precision. This uncertainty is evaluated with the help of the Hessian matrix analysis employed in Ref. Pumplin et al. (2002). Specifically, $`\delta \sigma _{PDF}`$ corresponds to the maximal variation of $`\sigma `$ for all possible PDF parametrizations lying in the acceptable range of $`\chi ^2`$. According to Table 1, the relative PDF uncertainties are of order $`25\%`$ at $`\sqrt{s}=200`$ GeV and $`710\%`$ at $`\sqrt{s}=500`$ GeV. These are the errors that exist prior to the measurements at RHIC.

The above PDF uncertainties are large because $`W`$ and $`Z`$ boson production at the $`pp`$ collider RHIC probes the sea quark distributions (primarily $`\overline{u}(x)`$ and $`\overline{d}(x)`$) at $`x0.1`$. For instance, the major part of $`W`$ boson events occurs in the central-rapidity region $`y_W0`$ and corresponds to the Born-level momentum fractions $`x=Q/\sqrt{s}`$ close to $`0.4`$ and $`0.16`$ at $`\sqrt{s}=200`$ GeV and $`\sqrt{s}=500`$ GeV, respectively. At such $`x`$ the sea quark PDFs are not as constrained by the existing (mostly DIS) data, which better determine the valence-dominated distributions $`u(x)`$ and $`d(x).`$ Therefore, RHIC provides the information that is complementary to the DIS experiments at HERA and collider experiments at the $`p\overline{p}`$ collider Tevatron.

This point is illustrated by Fig. 1, which shows the relative PDF uncertainty $`\delta \sigma _{PDF}/\sigma `$ in $`W^\pm `$ and $`Z^0`$ boson production at various existing colliders. To compare the potential of $`pp`$ and $`p\overline{p}`$ colliders of the same center-of-mass energy, this figure also shows the PDF uncertainties at fictitious $`p\overline{p}`$ colliders with $`\sqrt{s}=200`$ and $`500`$ GeV. According to the figure, the PDF uncertainties at RHIC substantially exceed the uncertainties at the Tevatron and LHC. They are also up to 2.5 times larger than the PDF uncertainties at the fictitious $`p\overline{p}`$ colliders of RHIC energies. The PDF uncertainty is the largest in $`W^{}`$ boson production at RHIC, which is sensitive primarily to the $`d\overline{u}`$ combination of the parton distributions.

On the other hand, for the planned integrated luminosities, RHIC will measure the spin-averaged rates with better accuracy, particularly in $`W^+`$ boson production. At $`\sqrt{s}=500`$ GeV and $`=800\text{ pb}^1`$, the $`2\sigma `$ relative statistical errors ($`2N^{1/2}`$) in the rate are $`0.6\%,`$ $`1.2\%`$, and $`2\%`$ in $`W^+`$, $`W^{}`$, and $`Z^0`$ boson production, respectively. At $`\sqrt{s}=200`$ GeV and $`=320\text{ pb}^1`$, the respective errors are $`10\%`$, $`18\%,`$ and $`44\%`$. We see that the projected statistical errors in $`W`$ boson production are smaller than the current PDF errors. Therefore, we conclude that the measurement of the spin-averaged cross sections at RHIC will drastically reduce the uncertainty in the unpolarized parton distributions in the probed region of $`x`$. We emphasize that the measurement of the spin-averaged cross sections should be an important part of the RHIC physics program, since it complements information about the unpolarized PDFs from the other colliders and reduces systematic uncertainties in the measured polarized cross sections.

## IV Reconstruction of spin asymmetries from RHIC data

### IV.1 Impact of multiple parton radiation on differential distributions

To extract useful information from the raw data in $`W`$ boson production, it is necessary to unfold the kinematical acceptance of the signal events, which can only be reliably predicted by a resummation calculation, similar to what has been done for the $`W`$ boson physics at the Tevatron Balazs and Yuan (1997). The correct lepton-level calculation is even more indispensable given that RHIC cannot register the transverse energy carried by the decay neutrino.

For definiteness, let us concentrate on the production of $`W^+`$ bosons. Of course, all conclusions in this section hold for $`W^{}`$ boson production as well. Fig. 2 shows the transverse momentum distributions of the $`W^+`$ boson and charged lepton from its decay, as predicted by the resummation package Legacy-RhicBos. Here and after, the numerical results correspond to the c.m. energy $`\sqrt{s}=500`$ GeV and integrated luminosity $`=800\text{ pb}^1`$. It is evident that most of the $`W`$ bosons are produced with small, but non-zero transverse momenta. Such non-zero $`q_T`$ is acquired through radiation of soft and collinear partons, which cannot be approximated by finite-order perturbative calculations. In order to obtain reliable predictions for differential cross sections, dominant logarithmic terms $`\alpha _S^n\mathrm{ln}^m\left(q_T^2/Q^2\right)`$ (where $`0m2n1`$) associated with such radiation should be summed through all orders of the perturbative series.

Furthermore, consider the distribution $`d\sigma /dp_T\mathrm{}`$ in Fig. 2(b). If the transverse momentum $`q_T`$ and width $`\mathrm{\Gamma }_W`$ of the $`W`$ boson were negligible (as in the Born-level calculation), $`d\sigma /dp_T\mathrm{}`$ would show a singular Jacobian peak exactly at $`p_T\mathrm{}=M_W/2`$. The smearing of the Jacobian peak in Fig. 2(b) is due to the multiple soft gluon radiation and non-zero width of the $`W`$ boson.

The example of $`d\sigma /dp_T\mathrm{}`$ illustrates that multiple parton radiation and heavy boson decay may drastically alter the Born-level predictions. The resummation effects are important if the distribution is sensitive to $`q_T`$ directly (as $`d\sigma /dp_T\mathrm{}`$), or if the phase space is constrained in a region with slower convergence of perturbation theory (as $`d\sigma /dy_{\mathrm{}}`$ at large $`|y_{\mathrm{}}|`$). In the remainder of this section, we utilize this accurate calculation to evaluate feasibility of the measurement of the boson-level asymmetry $`A_L(y_W)`$, which is commonly discussed in literature as the most convenient observable to probe polarized sea quark distributions. Many aspects of our discussion are based on the analysis of two-variable distributions $`d^2\sigma /(dp_T\mathrm{}dy_{\mathrm{}})`$ in the transverse momentum and rapidity of the observed charged lepton. These distributions for the unpolarized and single-spin $`W^+`$ and $`W^{}`$ boson production are shown in Fig. 3.

### IV.2 $`W`$ boson rapidity at the Born level

Consider the Born-level expression for the following asymmetry in $`W^+`$ boson production induced by the up quark and down antiquark via $`u\overline{d}W^+\mathrm{}^+\nu _{\mathrm{}}`$:

$`A_L(y_W,\theta ^{})`$ $``$ $`{\displaystyle \frac{d\mathrm{\Delta }_L\sigma /(dQ^2dy_Wd\mathrm{cos}\theta ^{})}{d\sigma /(dQ^2dy_Wd\mathrm{cos}\theta ^{})}}`$ (8)
$`=`$ $`{\displaystyle \frac{\mathrm{\Delta }u(x_A)\overline{d}(x_B)(1+\mathrm{cos}\theta ^{})^2+\mathrm{\Delta }\overline{d}(x_A)u(x_B)(1\mathrm{cos}\theta ^{})^2}{u(x_A)\overline{d}(x_B)(1+\mathrm{cos}\theta ^{})^2+\overline{d}(x_A)u(x_B)(1\mathrm{cos}\theta ^{})^2}}.`$
Here $`x_{A,B}(Q/\sqrt{s})e^{\pm y_W},`$ and $`\theta ^{}`$ is the polar angle of the charged lepton in the rest frame of the $`W`$ boson, with the $`z`$ axis pointing in the moving direction of the polarized proton beam. The single-spin asymmetry $`A_L(y_W)`$ in the rapidity of the $`W`$ boson can be obtained from Eq. (8) by integrating $`\mathrm{cos}\theta ^{}`$ out, i.e., integrating over the moving direction of $`\mathrm{}^+`$ and $`\nu _{\mathrm{}}`$:

$`A_L(y_W)`$ $`=`$ $`{\displaystyle \frac{_1^1\frac{d\mathrm{\Delta }_L\sigma }{dQ^2dy_Wd\mathrm{cos}\theta ^{}}d\mathrm{cos}\theta ^{}}{_1^1\frac{d\sigma }{dQ^2dy_Wd\mathrm{cos}\theta ^{}}d\mathrm{cos}\theta ^{}}}={\displaystyle \frac{\mathrm{\Delta }u(x_A)\overline{d}(x_B)+\mathrm{\Delta }\overline{d}(x_A)u(x_B)}{u(x_A)\overline{d}(x_B)+\overline{d}(x_A)u(x_B)}}.`$ (9)
The resulting expression for $`A_L(y_W)`$ is commonly discussed in literature as an example of a simple observable, which is straightforwardly sensitive to $`\mathrm{\Delta }u(x)`$ and $`\mathrm{\Delta }\overline{d}(x).`$ Similarly, the Born-level asymmetry $`A_L(y_W)`$ in $`W^{}`$ boson production is straightforwardly sensitive to $`\mathrm{\Delta }\overline{u}(x)`$ and $`\mathrm{\Delta }d(x)`$. It has been also argued in Ref. Gehrmann (1998) that $`𝒪(\alpha _S)`$ corrections cancel to a good degree between the numerator and denominator of $`A_L(y_W).`$ Unfortunately, as we will now demonstrate, the measurement of $`A_L(y_W)`$ is obstructed by the effects of $`W`$ boson decay and the limited detector acceptance, which seriously complicate both the reconstruction of $`y_W`$ and determination of $`A_L(y_W)`$ from the observed data.

Let us first discuss the determination of the rapidity of the $`W`$ boson. Given that neither of the two RHIC detectors measures the missing energy, the four-momentum of the $`W`$ boson cannot be deduced from the momenta of its decay products. Therefore the rapidity of the $`W`$ boson cannot be measured directly.

Despite the impossibility to measure $`y_W`$ in general, it has been proposed Bland (2000); Bunce et al. (2000) to statistically choose the correct $`y_W`$ *in a certain kinematical region* based on several assumptions about the QCD dynamics of the process. Let us outline the main idea of this method.

Denote the rapidity and transverse momentum of the lepton in the lab frame as $`y_{\mathrm{}}`$ and $`p_T\mathrm{}`$, respectively. Similarly, denote the corresponding variables in the rest frame of the $`W`$ boson as $`y_{\mathrm{}}^{}`$ and $`p_T\mathrm{}^{}`$, respectively. In the naive Born approximation, the $`W`$ bosons have zero width ($`Q=M_W`$) and zero transverse momentum ($`q_T=0`$). Therefore, the lepton variables in the different reference frames are related as

$`y_{\mathrm{}}`$ $`=`$ $`y_{\mathrm{}}^{}+y_W,`$ (10)
$`p_T\mathrm{}`$ $`=`$ $`q_T+p_T\mathrm{}^{}\begin{array}{c}=\\ \stackrel{}{q_T=0}\end{array}p_T\mathrm{}^{},`$ (13)
where $`y_{\mathrm{}}^{}`$ and $`p_T\mathrm{}^{}`$ are the functions of the polar angle $`\theta ^{}`$ of the lepton in the $`W`$ boson rest frame:

$`y_{\mathrm{}}^{}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\mathrm{ln}{\displaystyle \frac{1+\mathrm{cos}\theta ^{}}{1\mathrm{cos}\theta ^{}}},`$ (14)
$`p_T\mathrm{}^{}`$ $`=`$ $`{\displaystyle \frac{Q}{2}}\mathrm{sin}\theta ^{}\begin{array}{c}\begin{array}{c}=\\ \stackrel{}{Q=M_W}\end{array}\end{array}{\displaystyle \frac{M_W}{2}}\mathrm{sin}\theta ^{}.`$ (18)
The comments below the equations show the assumptions involved in obtaining the rightmost expressions in Eqs. (13) and (18).

For a given set of $`y_{\mathrm{}}`$ and $`p_T\mathrm{}`$, Eqs. (10)-(18) imply two solutions $`y_+`$ and $`y_{}`$ (corresponding to $`\theta ^{}`$ and $`\pi \theta ^{}`$) for $`y_W`$. One of them can be correctly chosen if the magnitude of the rapidity $`y_{\mathrm{}}`$ of the charged lepton is large (i.e., if the lepton is observed in the forward or backward rapidity region). In these special regions, one of the two solutions for $`y_W`$ can always be discarded, because it takes a value outside of the allowed rapidity range of the collider. In other words, the chosen solution is the one with the smaller magnitude. However, in reality, $`q_T`$ never vanishes, as predicted by the resummation calculation (see Fig. 2(a)). Similarly, $`\mathrm{\Gamma }_W`$ is not zero. Therefore, the above approximation never holds exactly, and the size of the error is determined by the difference between the Born-level and exact dynamics.

By comparing the leading order, NLO, and resummation calculations, we have found that the approximation holds reasonably well when $`q_TM_W`$, $`p_T\mathrm{}^{}p_T\mathrm{}`$, and $`y_{\mathrm{}}`$ is large. This point is illustrated in Fig. 4, which shows the statistically averaged differences $`(|y_\pm y_W|)\sigma /\sigma `$ between the exact $`y_W`$ and each of the two solutions $`y_\pm `$ in the two-dimensional plane of $`y_{\mathrm{}}`$ and $`p_T\mathrm{}`$, predicted by the resummation calculation. It can be seen that the difference is small for one of the solutions in the respective region $`y_{\mathrm{}}1`$ or $`y_{\mathrm{}}1`$, and $`p_T\mathrm{}2530`$ GeV. At the same time, the difference is unacceptably large for the other solution. Correspondingly, in each special region the distribution in the correct solution $`y_+`$ or $`y_{}`$ closely reproduces the distribution in true $`y_W`$. For example, Fig. 5 shows that the distributions $`d\sigma /dy_{}`$ agrees well with the distribution $`d\sigma /dy_W`$ in the region $`10\text{ GeV}p_T\mathrm{}30\text{ GeV}`$ and $`1.2y_{\mathrm{}}2.4`$.

Except for the special regions in Fig. 4, the approximate solutions $`y_\pm `$ do not agree well with the exact $`y_W`$. An important feature to note is that in the special regions the event rate is small. As shown in Fig. 3, most of the event rate comes from $`p_T\mathrm{}M_W/2`$ and $`y_{\mathrm{}}0`$. Furthermore, the magnitude of the error in $`y_W`$ depends on the shape of the cross section, which, of course, is unknown *a priori* in the polarized case.

### IV.3 Distortion effects on $`A_L(y_W`$)

Let us now turn to another problem in the measurement of $`A_L(y_W)`$: distortions in the deduced shape of $`A_L(y_W)`$ due the experimental cuts. To obtain $`A_L(y_W)`$, the experimental cross sections have to be integrated over all possible directions of the lepton motion, i.e., over $`\mathrm{cos}\theta ^{}`$ in the interval $`1\mathrm{cos}\theta ^{}1`$ (see Eq. (9) for an example of such integration in the Born-level analysis). Unfortunately, the acceptance of RHIC detectors does not permit such integration, so that the resulting asymmetry has the residual dependence on the acceptance range in $`\mathrm{cos}\theta _{}`$ and does not agree with the true asymmetry $`A_L(y_W)`$.

Indeed, $`\mathrm{cos}\theta ^{}`$ is related to the rapidity of the lepton $`y_{\mathrm{}}`$ via Eqs. (10) and (14):

$$y_{\mathrm{}}=y_W+\frac{1}{2}\mathrm{ln}\frac{1+\mathrm{cos}\theta ^{}}{1\mathrm{cos}\theta ^{}}.$$

Therefore, the restrictions on the range of $`y_{\mathrm{}}`$ reduce the range of integration over $`\mathrm{cos}\theta ^{}.`$ Specifically, the Phenix muon detectors are able to detect a muon if the muon rapidity $`y_\mu `$ and azimuthal angle $`\varphi _\mu `$ satisfy $`1.2|y_\mu |2.4`$ and $`0\varphi _\mu 2\pi `$. The Phenix electromagnetic calorimeter can register electrons with $`|y_e|0.35`$ and $`0\varphi _e\pi `$. The Star detector registers electrons in the barrel electromagnetic calorimeter that covers $`|y_e|1.0`$ and $`0\varphi _e2\pi `$. In addition, a selection cut $`p_T\mathrm{}1020`$ GeV is expected to be imposed on the charged leptons to suppress background contributions, which also constrains $`\mathrm{sin}\theta ^{}`$ through Eqs. (13) and (18).

The effect of the above cuts is demonstrated in Fig. 6, which shows the single-spin asymmetry in $`W^+`$ boson production calculated without constraints on $`y_{\mathrm{}}`$ and $`p_T\mathrm{}`$ (solid line), as well as with the constraints $`1.2<|y_{\mathrm{}}|<2.4,`$ $`p_T\mathrm{}>20\text{ GeV}`$ (circles) and $`|y_{\mathrm{}}|<1,`$ $`p_T\mathrm{}>20`$ GeV (boxes). The results shown here are derived using the resummation calculation.222For the asymmetries $`A_L(y_W)`$, the resummation and NLO calculations give close predictions (within one estimated statistical error $`\delta A_L`$), even though these calculations predict quite different shapes for the unpolarized and polarized differential cross sections.  According to the figure, there is a substantial difference between the asymmetries calculated with and without experimental cuts. This difference arises due to the different dependence of the unpolarized and polarized cross sections on angular distributions of the charged leptons, which affects the asymmetry because $`\mathrm{cos}\theta ^{}`$ is not integrated out completely.

The findings in this section can be summarized as follows. The extraction of $`A_L(y_W)`$ from RHIC data is not an easy task, because (a) the direct measurement of $`y_W`$ is impossible, and (b) the reconstructed $`A_L(y_W)`$ is distorted by the experimental acceptance cuts. In particular, the two-fold ambiguity in the determination of $`y_W`$ is unavoidable at $`y_{\mathrm{}}0`$ and $`p_T\mathrm{}M_W/2,`$ i.e., in the region with the largest cross section. To take the full advantage of the large event rate at RHIC, we ought to seek an alternative to $`A_L(y_W)`$, which is free of the complications discussed above. In the next section, we turn our attention to the spin asymmetries of the lepton-level cross sections $`d\sigma /dy_{\mathrm{}}`$ and $`d\sigma /dp_T\mathrm{}`$, which may serve as such alternative observables. Based on the resummation calculation, we analyze the sensitivity of the lepton-level asymmetries to the choice of the experimental cuts. We also present theoretical predictions for various PDF sets.

## V Single-spin asymmetries at the lepton level

### V.1 Impact of the cuts and soft radiation on lepton-level cross sections

Let us first separately discuss the unpolarized and polarized cross sections. To study the detector effects on the shape of various distributions, we will compare the cross sections for the kinematical acceptance of Phenix and Star detectors, as well as cross sections without any kinematical cuts. We also compare the resummation calculation to the NLO calculation, with the latter performed in the phase space slicing method for the separation parameter $`q_T^{sep}=1.6`$ GeV. We concentrate on the distributions $`d\sigma /dp_T\mathrm{}`$, where the differences between the resummed and NLO predictions are the most visible.

Fig. 7 presents the distribution $`d\sigma /p_T\mathrm{}`$ for the positively charged leptons from the decay of $`W^+`$ bosons. The shape of the NLO cross sections (dashed curves) for $`p_T\mathrm{}`$ around and above $`M_W/2`$ (i.e., the Jacobian peak) is strongly affected by the arbitrary parameter $`q_T^{sep}`$ employed in this calculation. Hence, the only reliable predictions for this part of the phase space are given by the resummed cross sections (shown as the solid curves). For smaller $`p_T\mathrm{}`$, the resummation calculation agrees with the NLO calculation, because it is formulated to match the NLO calculation in the region where the multiple gluon emission is not important. Furthermore, without any kinematic cuts, the integrated rate of the resummation calculation agrees with that of the NLO calculation to a percent level. As shown in Fig. 7(b), the selection of the leptons in the large rapidity region distorts the shape of the $`p_T\mathrm{}`$ distribution by suppressing the contributions near the Jacobian peak. On the other hand, the selection of the central-rapidity region $`|y_e|<1.0`$ (cf. Fig. 7(c)) suppresses the contributions with $`p_T\mathrm{}10`$ GeV. Both features can be understood from the shape of the differential distribution $`d^2\sigma /(dp_T\mathrm{}dy_{\mathrm{}})`$ given in Fig. 3(a). According to this figure, the small (large) rapidity $`|y_{\mathrm{}}|`$ statistically corresponds to the large (small) transverse momentum $`p_T\mathrm{}`$.

As a $`pp`$ collider, RHIC produces different numbers of $`W^{}`$ and $`W^+`$ bosons because of different parton luminosities. Due to the different rapidity distributions for $`W^+`$ and $`W^{}`$ bosons, the rapidity cuts have a different effect on $`d\sigma /dp_T\mathrm{}`$ in these processes. Fig. 8 shows the distributions $`d\sigma /dp_T\mathrm{}`$ in the process $`pp(W^{}\mathrm{}^{}\overline{\nu }_{\mathrm{}})X`$. We find that the shape of the distributions without cuts shown in Fig. 8(a) is very close to the shape of the corresponding distribution in Fig. 7(a), even though the overall normalizations are obviously different. On the other hand, the shapes do not coincide in the presence of the cuts, in particular in the case of the Phenix selection cuts shown in Figs. 7(b) and 8(b).

The $`p_T\mathrm{}`$ distributions in the single-spin $`W^+`$ boson production are shown in Fig. 9.

The cross sections $`d\mathrm{\Delta }_L\sigma /dp_T\mathrm{}`$ in $`W^+`$ boson production are negative, because the charged weak coupling is purely left-handed, and the dominant polarized PDF $`\mathrm{\Delta }u(x)`$ is positive in the RHIC range of $`x`$. Again, we see that the NLO calculation cannot reliably predict the distribution of $`p_T\mathrm{}`$ around and above $`M_W/2`$. For the Phenix kinematics, cf. Fig. 9(b), a sharper peak is developed in the distribution of $`p_T\mathrm{}`$ around $`M_W/2`$ as compared to the unpolarized case, cf. Fig. 7(b). The different shape is caused by the different rapidity dependence in the unpolarized and single-spin cases, as shown in Figs. 3(a) and 3(c). Finally, the distribution $`d\sigma /dp_T\mathrm{}`$ for the charged leptons $`\mathrm{}^{}`$ in the $`W^{}`$ events is shown in Fig. 10.

For the rapidity distributions $`d\sigma /dy_{\mathrm{}}`$ without imposed cuts on $`p_T\mathrm{}`$, the NLO and resummed cross sections are close to one another, in accordance with the cancellation of soft contributions in sufficiently inclusive observables. However, according to the discussion above, the resummation and NLO predictions for $`d\sigma /dy_{\mathrm{}}`$ can be very different if cuts on $`p_T\mathrm{}`$ are imposed.

### V.2 Single-spin asymmetries

We now turn our attention to the single-spin asymmetry. Since the only straightforward signature of the $`W`$ bosons at RHIC is the observation of the secondary charged leptons, it is important to understand these asymmetries at the lepton level. In Figs. 11 and 12, we show the lepton-level asymmetries $`A_L(y_{\mathrm{}})`$ in the $`W^+`$ and $`W^{}`$ events, respectively, for various cuts on $`p_T\mathrm{}`$. In these figures, the positive rapidity direction is in the moving direction of the polarized proton beam. Similarly, Figs. 13 and 14 show the asymmetries $`A_L(p_T\mathrm{})`$ for various cuts on $`y_{\mathrm{}}`$. All asymmetries are derived using the Gehrmann-Stirling PDF sets A and B Gehrmann and Stirling (1996) and GRSV-2000 valence-like PDF set Gluck et al. (2001) in the resummation calculation. We found by experimenting with other available PDF sets in Refs. Gehrmann and Stirling (1996); Gluck et al. (1996, 2001) that the shown graphs characterize well the possible span in the asymmetry due to the different polarized PDF parametrizations. The projected statistical errors are estimated according to Eq. (7) for the bin sizes shown in the figure.

According to the figures, the lepton-level asymmetries can be accurately measured for both $`W^+`$ and $`W^{}`$ bosons. These directly observed asymmetries can efficiently discriminate between different PDF sets; hence, they provide a viable alternative to the less accessible asymmetry $`A_L(y_W)`$. For both $`W^+`$and $`W^{}`$ boson samples (cf. Figs. 11 and 12), the power to discriminate between different PDF sets is larger in the region of the positive rapidity $`y_{\mathrm{}}`$, especially for the $`W^{}`$ boson data. It is interesting to observe that, in different ranges of $`y_{\mathrm{}},`$ the largest variations due to the choice of the PDFs appear in different ranges of $`p_T\mathrm{}`$. In $`W^+`$ boson production, the largest variation in the asymmetries $`A_L(p_T\mathrm{})`$ is concentrated at $`p_T\mathrm{}3040`$ GeV in the backward rapidity region $`y_{\mathrm{}}<0`$ (Fig. 13(a) and (b)) and at $`p_T\mathrm{}3040`$ GeV in the forward rapidity region $`y_{\mathrm{}}>0`$ (Fig. 13(c) and 13(d)). In $`W^{}`$ boson production the largest variation is observed for very forward leptons, $`1.2<y_{\mathrm{}}<2.4`$ (Fig. 14(d)). If an additional condition $`p_T\mathrm{}>20`$ GeV is imposed, other regions of $`y_{\mathrm{}}`$ (cf. Fig. 14(a)-14(c)) show a smaller sensitivity to different PDF sets.

During the first years of RHIC spin program, it may not be possible to distinguish between the positively and negatively charged leptons from the $`W`$ boson decay, particularly before the solenoid field of the Star detector is fully operational. It is therefore interesting to examine the possibility of measuring the single-spin asymmetries in the combined sample of $`W^^+`$ and $`W^{}`$ bosons, without distinguishing between charged leptons and anti-leptons. Such asymmetries are shown in Figs. 15 and 16. Since there are about three $`W^+`$ bosons produced per each $`W^{}`$ boson, these asymmetries are close to the $`W^+`$ boson asymmetries shown in Figs. 11 and 12. Therefore, these asymmetries will be able to distinguish between different PDF sets for a sufficiently high number of events.333Of course, the projected error bars in Figs. 15 and 16 (evaluated for $`=800\text{ pb}^1`$) will have to be adjusted to agree with the actual integrated luminosity at the time of the measurement. However, they cannot completely replace the single-spin asymmetries for independent $`W^+`$ and $`W^{}`$ boson samples, which have better discriminating power and allow us to distinguish between contributions of different parton flavors.

### V.3 Nonperturbative contribution

As was discussed above, the transverse momentum distributions for vector bosons are affected by the soft and collinear radiation, which can be described only by means of all-order resummation. In addition, the distributions at very small $`q_T`$ are sensitive to nonperturbative contributions characterized by large impact parameters $`b1\text{ GeV}^1`$. As a result, the shape of the lepton-level distribution $`d\sigma /dp_T\mathrm{}`$ around its peak at about $`p_T\mathrm{}=M_W/2`$ is affected by both perturbative and nonperturbative QCD radiation, in addition to the nonzero width of the $`W`$ boson. While the perturbative soft and collinear contributions can be calculated order-by-order in the resummation formalism, the nonperturbative contributions (associated with the nonperturbative transverse motion of partons inside the proton) can only be approximated by phenomenological parametrizations. In particular, based on the factorization properties of the resummed cross sections, one can hypothesize a simple relationship between the nonperturbative Sudakov factors in the unpolarized, single-spin, and double-spin vector boson production:

$$S_{NP}|_{singlespin}=\frac{1}{2}\left(S_{NP}|_{unpolarized}+S_{NP}|_{doublespin}\right).$$
(19)
Furthermore, $`S_{NP}`$ is usually assumed to be independent of the type of the electroweak hard probe, so that it is the same in $`\gamma ^{},`$ $`W,`$ and $`Z`$ boson production. If the relationship (19) for the spin dependence of $`S_{NP}`$ indeed holds in nature, the nonperturbative contributions in single-spin $`W`$ boson production can be determined unambiguously once the parametrizations of $`S_{NP}`$ on the right-hand side of Eq. (19) are available, e.g., from the global analysis of the distributions $`d\sigma /dq_T`$ in the unpolarized and double-spin $`\gamma ^{}`$ or $`Z^0`$ production.

It is well known from the unpolarized studies (see, e.g., Refs. Collins et al. (1985); Davies et al. (1985); Arnold and Kauffman (1991)) that the importance of the nonperturbative contributions in the massive boson production at high-energy colliders is strongly reduced as compared to Drell-Yan pair production in fixed-target experiments. Nonetheless, this dependence is not entirely negligible at small transverse momenta of $`W`$ bosons, roughly for $`q_T5`$ GeV or less. To give an idea about the typical size of the nonperturbative contributions, Fig. 17 shows the event rate for the single-spin cross section $`d\mathrm{\Delta }_L\sigma /dp_T\mathrm{}`$ at $`=800\text{ pb}^1`$, integrated over the Star acceptance range $`1y_{\mathrm{}}1`$. This rate was calculated using two parameterizations Ladinsky and Yuan (1994); Landry et al. (2002) of the nonperturbative part, which were found from the analysis of unpolarized vector boson production. Correspondingly, the resummed cross sections were calculated with the assumption that $`S_{NP}|_{singlespin}=S_{NP}|_{unpolarized}`$. It can be seen that the sensitivity to the nonperturbative input is small, but, nonetheless, visible near the Jacobian peak. For comparison, we also included the NLO cross section calculated using the phase space slicing method. The NLO curve substantially deviates from the resummation curves, and, moreover, its shape can be drastically modified by varying the phase space slicing parameter $`q_T^{sep}`$. In contrast, practically all features of the resummation curve are determined unambiguously by the perturbative contribution, and the remaining small dependence on the nonperturbative contributions is pinned down by utilizing a phenomenological parameterization discussed above. Needless to say, RHIC should explore the spin properties of nonperturbative contributions by studying distributions $`d\sigma /dq_T`$ in the double-spin $`\gamma ^{}`$ and $`Z^0`$ production and $`d\sigma /dp_T\mathrm{}`$ in single-spin and double-spin $`W`$ boson production.

## VI Conclusion

Due to the maximal parity violation in the weak current coupling, the longitudinal single-spin asymmetry of $`W`$ boson production is strongly sensitive to the polarized parton distributions. At the center-of-mass energy of 500 GeV, $`W`$ boson production events at RHIC can clearly discriminate between various polarized quark luminosities. As a $`pp`$ collider, RHIC probes the sea quark PDFs in the region of relatively high $`x,`$ where such PDFs are not adequately constrained even in the unpolarized case. Therefore, RHIC experiments also have a good potential to reduce the uncertainties in the unpolarized sea quark PDFs in a dedicated study of spin-averaged cross sections.

Since Phenix and Star detectors do not determine the missing transverse energy associated with the neutrino from $`W`$ boson decay, and because the charged leptons from the decay can only be observed in a part of the $`4\pi `$ solid angle, the reconstruction of the asymmetries at the level of $`W`$ bosons is not easy to realize. As an alternative to the asymmetry $`A_L(y_W)`$ for the $`W`$ boson rapidity distribution $`d\sigma /dy_W`$, we advocate the measurement of the asymmetries $`A_L(y_{\mathrm{}}`$) and $`A_L(p_T\mathrm{})`$ for the distributions $`d\sigma /dy_{\mathrm{}}`$ and $`d\sigma /dp_T\mathrm{}`$ in the rapidity $`y_{\mathrm{}}`$ and transverse momentum $`p_T\mathrm{}`$ of the secondary charged lepton. These lepton-level asymmetries can be directly observed at RHIC detectors, they are sizeable, and they can be used to determine the polarized PDFs with high precision.

To reliably predict the distributions of the charged leptons in all available phase space, a theoretical calculation is needed for all-order summation of large logarithmic contributions corresponding to $`W`$ bosons with a small transverse momentum. At the lepton level, the largest effect of the logarithm resummation is observed in the distributions $`d\sigma /dp_T\mathrm{}`$ and $`d\mathrm{\Delta }_L\sigma /dp_T\mathrm{}`$ in the proximity of the Jacobian peak ($`p_T\mathrm{}M_W/2`$). The importance of the resummation can be verified by observing a large difference between the resummed and NLO cross sections near the Jacobian peak, which happens only when the effect of soft gluon radiation is important.

In this paper, we apply the state-of-art resummation calculation Nadolsky and Yuan (2003) for spin-dependent vector boson production to predict reliably the distributions of the observed leptons from $`W`$ boson decay. In addition to the spin-averaged and single-spin cross sections, we presented a detailed study of single-spin asymmetries for the decay lepton in the presence of kinematical cuts imposed by Phenix and Star detectors. Furthermore, we propose to study the distribution $`d\mathrm{\Delta }_L\sigma /dp_T\mathrm{}`$ at $`p_T\mathrm{}M_W/2`$ to learn about the spin (in)dependence of nonperturbative contributions in the CSS resummation formalism. The study of such contributions for different beam polarizations and bosons of different types ($`\gamma ^{},`$$`W^\pm `$, $`Z^0`$) will shed light on basic properties of QCD factorization and provide important clues about the nature of the intrinsic transverse motion of partons inside the proton.

## Acknowledgments

Authors would like to thank C. Balazs, D. Boer, G. Bunce, M. Grosse Perdekamp, S. Gupta, J. Kiryluk, N. Saito, M. Stratmann, W. Vogelsang, and members of the CTEQ collaboration for helpful discussions. We are grateful to F. Olness for his comments on the manuscript. We thank the organizers of RHIC Spin workshops, where preliminary results of this work were presented. The work of P. M. N. has been supported by the U.S. Department of Energy and Lightner-Sams Foundation. The research of C. P. Y. has been supported by the National Science Foundation under grant PHY-0100677.