text:
# hep-ph/0102071RIKEN-BNL preprint Sudakov suppression in azimuthal spin asymmetries

## I Introduction

In this article we will study the effects of Sudakov factors in transverse momentum dependent azimuthal spin asymmetries, like the Collins effect asymmetry of Ref. . We will demonstrate explicitly that such asymmetries suffer from suppression due to these Sudakov factors, in the region where the transverse momentum is much smaller than the large energy scale $`Q^2`$. This Sudakov suppression stems from soft gluon radiation and increases with energy. It implies that tree level estimates of such asymmetries tend to overestimate magnitudes and increasingly so with rising energy. In this paper, the $`Q^2`$ dependence of the suppression due to Sudakov factors will be studied numerically in two examples which are relevant for present-day experimental studies. In addition, the asymptotic $`Q^2`$ dependence of an important subset of asymmetries is calculated analytically, providing an upper bound for the decrease with energy of other asymmetries.

The first example we consider is a Collins effect driven $`\mathrm{cos}(2\varphi )`$ asymmetry in electron-positron annihilation into two almost back-to-back pions , which in principle can be determined from existing LEP data. This asymmetry allows for a determination of the Collins effect fragmentation function, which in turn would be useful for the extraction of the transversity distribution function $`h_1`$. The latter can be done via the single spin asymmetry in semi-inclusive Deep Inelastic Scattering (DIS) which was originally proposed by Collins  –our second example–. In Ref.  Collins already remarks that the Sudakov factors will have the effect of diluting this single spin asymmetry due to broadening of the transverse momentum distribution by soft gluon emission. We will study this effect in a quantitative way to gain insight into the parametric dependences of such Sudakov suppression.

For this we will follow the recent analysis  of a helicity non-flip double transverse spin asymmetry in vector boson production in hadron-hadron scattering. That asymmetry is for instance relevant for the polarized proton-proton collisions to be performed at BNL-RHIC. In Ref.  the effect of the Sudakov factors compared to the tree level asymmetry expression was numerically estimated. It was shown that the inclusion of Sudakov factors cause suppression by at least an order of magnitude compared to the tree level result at scales around $`M_W`$ or $`M_Z`$. Moreover, the suppression increases with energy approximately as $`Q^{0.6}`$ (in the studied range of roughly $`10100`$ GeV). The conclusion was that the Sudakov suppression together with a kinematic suppression (due to explicit lightcone momentum fractions that appear in the prefactor) imply that the asymmetry will be negligible for $`Z`$ or $`W`$ production and is interesting only at much lower energies.

Similarly, Sudakov suppression will turn out to be an important issue for the actual determination of the above mentioned $`\mathrm{cos}(2\varphi )`$ asymmetry from LEP data (or in general at high energy $`e^+e^{}`$ colliders). Due to the lack of knowledge of the nonperturbative Sudakov factor in the case of electron-positron annihilation into almost back-to-back hadrons, solid numerical conclusions about the size and $`Q^2`$ dependence of the suppression cannot be drawn. Nevertheless, one can draw conclusions about what determines the $`Q^2`$ dependence of the transverse momentum distribution of the asymmetries and about the size of the suppression for generic nonperturbative Sudakov factors. For that purpose comparisons to tree level are most instructive, since the evolution of the often unknown distribution and fragmentation functions themselves does not affect the relative Sudakov suppression compared to tree level.

The two examples to be investigated here –the above mentioned $`\mathrm{cos}(2\varphi )`$ asymmetry and the original Collins $`\mathrm{sin}(\varphi )`$ asymmetry in semi-inclusive DIS–, have different transverse momentum weights and angular dependences, producing different Sudakov suppression effects. General conclusions about how these properties affect the $`Q^2`$ dependence can be drawn. In general, a larger power of transverse momentum in the weight of an asymmetry implies larger suppression. Both the numerical calculations at realistic collider energies and the asymptotic behavior of asymmetries exhibit this feature.

Before going into the details of the specific examples, we will first give an overview of the issues regarding factorization theorems for transverse momentum dependent cross sections, in the region where the transverse momentum is very small compared to the hard scale in the process . This will provide the theoretical justification of the factorized asymmetry expressions that will be derived.

The outline of this paper is as follows. In section II we will discuss the essentials of factorization theorems for the Drell-Yan process  (relevant for the double transverse spin asymmetry studied in Ref. ), which are in fact completely analogous to the case of back-to-back jets in electron-positron annihilation  and semi-inclusive DIS (SIDIS). Following the approach of Ref. , we will then study the $`\mathrm{cos}(2\varphi )`$ asymmetry in electron-positron annihilation (section III) and the Collins $`\mathrm{sin}(\varphi )`$ asymmetry in semi-inclusive DIS (section IV). In section V we will discuss the asymptotic behavior of the Sudakov suppression for certain subsets of asymmetries.

## II Factorization and transverse momentum

For definiteness, we will consider here the Drell-Yan process ($`H_1+H_2\mathrm{}+\overline{\mathrm{}}+X`$) and discuss what is known about the factorization of this process. First of all, it is well-known that a local operator product expansion (OPE) cannot be applied, since the process receives relevant contributions off the light cone , i.e. from matrix elements like

$$\underset{X}{}P_1,P_2|J_\mu (\xi )|XX|J_\nu (0)|P_1,P_2\text{with}\xi ^20.$$
(1)
Therefore, one usually considers nonlocal operators and what is called a “working redefinition of twist” , which means the lowest value of $`t`$ in $`1/Q^{t2}`$ at which a function can contribute to the cross section ($`Q`$ is the invariant mass of the lepton pair). The nonlocal operators appear in the leading twist factorization theorem

$$\frac{d\sigma }{d\mathrm{\Omega }dx_1dx_2}=\underset{a,b}{}_{x_1}^1dx_{x_2}^1d\overline{x}\mathrm{\Phi }^a(x)H^{ab}(x,\overline{x};Q^2)\overline{\mathrm{\Phi }}{}_{}{}^{b}(\overline{x}),$$
(2)
where $`a,b`$ are flavor indices and $`H^{ab}(x,\overline{x};Q^2)`$ is the hard (partonic) part of the scattering, which is a function of the hard scale $`Q^2`$ and the lightcone momentum fractions $`x,\overline{x}`$ only. The correlation functions $`\mathrm{\Phi }`$ ($`\overline{\mathrm{\Phi }}{}_{}{}^{a}=\mathrm{\Phi }^{\overline{a}}`$) describe the soft (nonperturbative) physics and are given by nonlocal operator matrix elements

$`\mathrm{\Phi }_{ij}(x)`$ $``$ $`{\displaystyle \frac{d\lambda }{2\pi }e^{i\lambda x}P,S|\overline{\psi }_j(0)^+[0,\lambda n_{}]\psi _i(\lambda n_{})|P,S},`$ (3)
where the path-ordered exponential,

$$^+[0,\lambda n_{}]=𝒫\mathrm{exp}\left(ig_0^\lambda 𝑑\eta A^+(\eta n_{})\right),$$
(4)
renders the matrix element color gauge invariant. The hadronic state $`|P,S`$ is determined by the hadron momentum $`P`$ and spin vector $`S`$. Also, we have suppressed the factorization scale dependence in Eq. (2).

Secondly, if one considers the cross section differential in the transverse momentum $`𝒒_T`$ of the lepton pair, then for the case of $`𝒒_T^2Q_T^2Q^2`$ collinear factorization, like in Eq. (2), does not hold. Instead, the leading twist factorization takes the form

$`{\displaystyle \frac{d\sigma }{d\mathrm{\Omega }dx_1dx_2d^2𝒒_T}}`$ $`=`$ $`{\displaystyle \underset{a,b}{}}{\displaystyle _{x_1}^1}𝑑x{\displaystyle _{x_2}^1}𝑑\overline{x}{\displaystyle d^2𝒌_Td^2𝒑_T\frac{d^2𝒃}{(2\pi )^2}e^{i𝒃(𝒑_T+𝒌_T𝒒_T)}}`$ (6)
$`\times \mathrm{\Phi }^a(x,𝒌_T)H^{ab}(x,\overline{x};Q)e^{S(𝒃,Q)}\overline{\mathrm{\Phi }}{}_{}{}^{b}(\overline{x},𝒑_T)+Y(x_1,x_2,Q,Q_T).`$
The correlation function $`\mathrm{\Phi }^a(x,𝒌_T)`$ is now also a function of transverse momentum, the factor $`e^{S(𝒃,Q)}`$ is the so-called Sudakov form factor and the factor $`Y(x_1,x_2,Q,Q_T)`$ becomes important only when $`Q_TQ`$. Here we use a different notation than in e.g. Refs. , but to make a connection with for instance the notation of Ref.  we note that $`\mathrm{\Phi }(x,𝒌_T)`$ corresponds to the Fourier transform of $`\widehat{𝒫}_{A/i}(x_A,𝒃)\overline{)}P_A/2`$ and the $`e^S`$ terms correspond to each other.

The factor $`Y(x_1,x_2,Q,Q_T)`$ is present to retrieve collinear factorization when $`Q_T`$ is a hard scale itself, i.e. $`Q_TQ`$. We will neglect the $`Y`$ term from here on, since we will only be interested in the region of $`Q_T^2Q^2`$. We will now comment on the other terms in the above factorization theorem and explain the additional restrictions we will impose in the limit $`Q_T^2Q^2`$.

The factor $`e^{S(𝒃,Q)}`$ is the so-called Sudakov form factor, which arises due to exponentiation of soft gluon contributions. This is in contrast to inclusive cross sections, like Eq. (2), in which there is a cancellation of soft gluon contributions. At values $`b^2=𝒃^21/\mathrm{\Lambda }^2`$, the Sudakov form factor is perturbatively calculable and is of the form

$$S(b,Q)=_{b_0^2/b^2}^{Q^2}\frac{d\mu ^2}{\mu ^2}\left[A(\alpha _s(\mu ))\mathrm{ln}\frac{Q^2}{\mu ^2}+B(\alpha _s(\mu ))\right],$$
(7)
where $`b_0=2\mathrm{exp}(\gamma _E)1.123`$ (we choose the usual constants $`C_1=b_0,C_2=1`$). One can expand the functions $`A`$ and $`B`$ in $`\alpha _s`$ and the first few coefficients are known, e.g. . In order to obtain a first estimate of the effect of including the Sudakov factor we will consider only the leading contribution, i.e. take into account only the first term in the expansion of $`A`$: $`A^{(1)}=C_F/\pi `$, which is the same for unpolarized as well as polarized scattering. This leads to the expression

$$S(b,Q)=\frac{16}{332n_f}\left[\mathrm{log}\left(\frac{b^2Q^2}{b_0^2}\right)+\mathrm{log}\left(\frac{Q^2}{\mathrm{\Lambda }^2}\right)\mathrm{log}\left[1\frac{\mathrm{log}\left(b^2Q^2/b_0^2\right)}{\mathrm{log}\left(Q^2/\mathrm{\Lambda }^2\right)}\right]\right].$$
(8)
We will take for the number of flavors $`n_f=5`$ and also $`\mathrm{\Lambda }_{QCD}=200\text{MeV}`$.

The correlation function $`\mathrm{\Phi }^a(x,𝒌_T)`$ is now defined as

$$\mathrm{\Phi }(x,𝒌_T)𝑑k^{}\left[\frac{d^4x}{(2\pi )^4}e^{ikx}P,S|\overline{\psi }(0)^+[0,\mathrm{}]^+[\mathrm{},x]\psi (x)|P,S\right]|_{k^+=xP^+,𝒌_T},$$
(9)
with for each colored field a path-ordered exponential like

$$^+[\mathrm{},x]=𝒫\mathrm{exp}\left(ig_{\mathrm{}}^x^{}𝑑y^{}A^+(y)\right)|_{y^+=x^+,𝒚_T=𝒙_T}.$$
(10)
The correlation function $`\mathrm{\Phi }^a(x,𝒌_T)`$ also has a factorization scale dependence ($`\mu =b_0/b`$), which is linked to that of the hard scattering part $`H`$.

The correlation function $`\mathrm{\Phi }^a(x,𝒌_T)`$ contains the nonperturbative dependence on the transverse momentum, which cannot be calculated perturbatively and should be fitted to experiment. In for example Refs.  it is explained how these functions can be replaced by the ordinary parton distribution functions $`\mathrm{\Phi }^a(x)`$, by introduction of a $`b`$-regulator, e.g. the usual cut-off $`b_{\mathrm{max}}`$, via $`b_{}=b/\sqrt{1+b^2/b_{\mathrm{max}}^2}`$, and by introducing a nonperturbative Sudakov factor $`S_{NP}`$. Schematically, this proceeds as follows. The first term in Eq. (6) can be written as an overall $`𝒃`$ integration of an integrand called $`\stackrel{~}{W}(b)`$,

$$\frac{d\sigma }{d\mathrm{\Omega }dx_1dx_2d^2𝒒_T}=_{x_1}^1𝑑x_{x_2}^1𝑑\overline{x}\frac{d^2𝒃}{(2\pi )^2}e^{i𝒃𝒒_T}\stackrel{~}{W}(b)+Y.$$
(11)
This function $`\stackrel{~}{W}(b)`$ can be trivially rewritten as $`\stackrel{~}{W}(b_{})\times (\stackrel{~}{W}(b)/\stackrel{~}{W}(b_{}))`$. The first term ($`\stackrel{~}{W}(b_{})`$) can be calculated within perturbation theory since $`\alpha _s(b_0/b_{})`$ is always small (we will take the usual $`b_{\mathrm{max}}=0.5\text{GeV}^1`$, hence $`\alpha _s(b_0/b_{}\stackrel{>}{}2)\stackrel{<}{}0.3`$) and the second term can be shown to be of the form

$$\frac{\stackrel{~}{W}(b)}{\stackrel{~}{W}(b_{})}=\mathrm{exp}\left(\mathrm{ln}(Q^2/Q_0^2)g_Q(b)g_A(x,b)g_B(\overline{x},b)\right)\mathrm{exp}\left(S_{NP}(b)\right).$$
(12)
The functions $`g`$ are not calculable in perturbation theory and need to be fitted to experiment. The functions $`g_A(x,b)`$ and $`g_B(\overline{x},b)`$ parameterize the intrinsic transverse momentum of the functions $`\stackrel{~}{\mathrm{\Phi }}{}_{}{}^{a}(x,b)`$ and $`\stackrel{~}{\overline{\mathrm{\Phi }}}{}_{}{}^{b}(\overline{x},b)`$ (in the polarized case $`\stackrel{~}{\mathrm{\Phi }}^a`$ is actually a function of $`𝒃`$, see further comments below). It will be important later on that the functions $`g_A`$ and $`g_B`$ are independent of the scale $`Q`$. Finally, a lightcone expansion allows to express $`\stackrel{~}{\mathrm{\Phi }}^a(x,b_{})`$ in terms of $`\mathrm{\Phi }^a(x)`$.

We will consider $`\stackrel{~}{\mathrm{\Phi }}^a(x,𝒃_{})\stackrel{~}{\overline{\mathrm{\Phi }}}{}_{}{}^{b}(\overline{x},𝒃_{})H(x,\overline{x};Q)`$ only to lowest order in $`\alpha _s`$. The reason is that the perturbative tail of $`\stackrel{~}{\mathrm{\Phi }}^a(x,𝒃_{})`$ is not known for the Collins function. Because of the unknown size of the distribution and/or fragmentation functions appearing in the asymmetries, the magnitude of the asymmetries cannot be estimated. Therefore, here we will be interested in the effect of soft gluons on the transverse momentum dependence of the asymmetries, which determines the $`Q^2`$ dependence of the $`Q_T`$ dependence. Although the hard scattering part is restricted to tree level, one can draw conclusions on the size of the effect of Sudakov factors by comparison to the tree level result, because both results will receive the same corrections to the hard part. Corrections to the hard scattering part can be easily inserted with hindsight (for a one-loop example cf. Ref. , pages 432 and 433). As said we are interested in the $`Q^2`$ evolution of the $`Q_T`$ dependence of the asymmetries. This $`Q^2`$ dependence arises in a nontrivial way via Bessel functions, such that a numerical study is warranted.

In Refs.  unpolarized scattering is considered, where only dependence on the length of $`𝒃`$ appears. Here we will consider polarized scattering and the above sketched procedure of Refs.  to factor a nonperturbative part from $`\stackrel{~}{W}(b)`$, can also be applied after $`\mathrm{\Phi }^a(x,𝒌_T)`$ is parameterized in terms of functions of $`𝒌_T^2`$. Explicit factors of $`𝒌_T𝑺_T`$ or $`\overline{)}𝒌_T`$ can be included in the hard scattering part $`H`$.

The above factorization theorem for the Drell-Yan process provides the justification of the factorized expression for the double transverse spin asymmetry studied in Ref. . In Ref.  an analogous factorization theorem had been discussed earlier for back-to-back jets in electron-positron annihilation, which comes down to a replacement of $`\mathrm{\Phi }(x,𝒌_T)\mathrm{\Delta }(z,𝒌_T)`$ –the fragmentation correlation function–. In Ref.  it is argued that a similar factorization theorem holds for semi-inclusive DIS. Therefore, factorization theorems analogous to the one discussed above (differing only by obvious replacements), provide the theoretical foundation for the examples to be considered below.

It is interesting to see what happens if one takes tree level everywhere in the above factorized expression Eq. (6) (ignoring the $`Y`$ term). At tree level one finds the reduction

$$\frac{d^2𝒃}{(2\pi )^2}e^{i𝒃(𝒑_T+𝒌_T𝒒_T)}e^{S(𝒃)}\delta ^2(𝒑_T+𝒌_T𝒒_T),$$
(13)
such that one arrives at the expression used by Ralston and Soper for the hadron tensor

$$𝒲^{\mu \nu }=\frac{1}{3}d^2𝒑_Td^2𝒌_T\delta ^2(𝒑_T+𝒌_T𝒒_T)\text{Tr}\left(\mathrm{\Phi }(x_1,𝒑_T)V_1^\mu \overline{\mathrm{\Phi }}(x_2,𝒌_T)V_2^\nu \right)|_{p^+,k^{}}+\left(\begin{array}{c}qq\\ \mu \nu \end{array}\right).$$
(14)
The vertex $`V_i^\mu `$ can be either the photon, $`Z`$ or $`W`$ boson vertex. In the above unpolarized case $`\mathrm{\Phi }(x,𝒑_T)`$ is parameterized as $`\frac{1}{2}f_1(x,𝒑_T^2)\overline{)}P`$, where $`f_1`$ is the parton momentum distribution function. Collins  has shown that polarization does not affect the factorization theorems and therefore, we conclude that the tree level formalism of Ralston and Soper Eq. (14) is in accordance with the factorization theorem Eq. (6) also in the case when polarization is taken into account. Clearly, Eq. (14) is only applicable for $`Q_T`$ values of the order of the intrinsic transverse momentum. To extend the formalism of Ralston and Soper to larger values of $`Q_T`$, but still under the restriction of $`Q_T^2Q^2`$, one can use the above factorized expression Eq. (6) (without the $`Y`$ term) beyond tree level.

## III Unpolarized asymmetry in electron-positron annihilation

In this section we investigate a $`\mathrm{cos}(2\varphi )`$ asymmetry in the process of electron-positron annihilation into two almost back-to-back pions $`e^+e^{}\pi ^+\pi ^{}X`$. This process is very similar to the Drell-Yan process, so we can employ the same expression as given in Eq. (6), but with the correlation functions $`\mathrm{\Phi }`$, replaced by the fragmentation correlation function

$$\mathrm{\Delta }_{ij}(k)=\underset{X}{}\frac{d^4z}{(2\pi )^4}e^{ikz}0|\psi _i(x)^{}[x,\mathrm{}]|P,S;XP,S;X|^{}[\mathrm{},0]\overline{\psi }_j(0)|0,$$
(15)
which is used to define transverse momentum dependent fragmentation functions.

The asymmetry requires observing the transverse momentum of the vector boson compared to the two pions. The tree level asymmetry expression for the $`\mathrm{cos}(2\varphi )`$ asymmetry in the process $`e^+e^{}\pi ^+\pi ^{}X`$ was discussed in Refs. , to which we refer for details. Here we shortly repeat the essentials. The asymmetry depends on the fragmentation function $`H_1^{}`$ associated with the Collins effect . The definition of the Collins function is given by

$$\mathrm{\Delta }(z,𝒌_T)=\frac{M}{4P^{}}\left\{D_1(z,𝒌_T^2)\frac{\overline{)}P}{M}+H_1^{}(z,𝒌_T^2)\frac{\sigma _{\mu \nu }k_T^\mu P^\nu }{M^2}\right\},$$
(16)
where we have only displayed the fragmentation functions that can be present for unpolarized hadron production ($`D_1`$ is the ordinary unpolarized fragmentation function). The $`\mathrm{cos}(2\varphi )`$ asymmetry is an azimuthal spin asymmetry in the sense that the asymmetry arises due to the correlation of the transverse spin states of the quark-antiquark pair. On average the quark and antiquark will not be transversely polarized, but for each particular event the transverse spin need not be zero and moreover, the spin states of the quark and antiquark are exactly correlated. Subsequently, the directions of the produced pions are correlated due to the Collins effect and this correlation does not average out after summing over all quark polarization states. The Collins effect correlates the azimuthal angle of the transverse spin of a fragmenting quark with that of the transverse momentum of the produced hadron (both taken around the quark momentum), via a $`\mathrm{sin}(\varphi )`$ distribution of their difference angle $`\varphi `$.

One finds for the leading order unpolarized cross section, taking into account both photon and $`Z`$-boson contributions ,

$`{\displaystyle \frac{d\sigma (e^+e^{}h_1h_2X)}{d\mathrm{\Omega }dz_1dz_2d^2𝒒_T}}={\displaystyle \frac{3\alpha ^2}{Q^2}}z_1^2z_2^2{\displaystyle \underset{a,\overline{a}}{}}\{K_1^a(y)\left[D_1\overline{D}_1\right]`$ (18)
$`+[K_3^a(y)\mathrm{cos}(2\varphi _1)+K_4^a(y)\mathrm{sin}(2\varphi _1)]\left[(2\widehat{𝒉}𝒑_T\widehat{𝒉}𝒌_T𝒑_T𝒌_T){\displaystyle \frac{H_1^{}\overline{H}_1^{}}{M_1M_2}}\right]\}.`$
The convolutions $``$ will be dicussed below (Eqs. (37) and (40)), but first we will give the definition of the functions $`K_i^a(y)`$ appearing in this expression and comment on the frame in which this cross section is expressed. As before, $`a`$ is the flavor index; the functions $`K_i^a(y)`$ are defined as

$`K_1^a(y)`$ $`=`$ $`A(y)\left[e_a^2+2g_V^le_ag_V^a\chi _1+c_1^lc_1^a\chi _2\right]{\displaystyle \frac{C(y)}{2}}\left[2g_A^le_ag_A^a\chi _1+c_3^lc_3^a\chi _2\right],`$ (19)
$`K_3^a(y)`$ $`=`$ $`B(y)\left[e_a^2+2g_V^le_ag_V^a\chi _1+c_1^lc_2^a\chi _2\right],`$ (20)
$`K_4^a(y)`$ $`=`$ $`B(y)\left[2g_V^le_ag_A^a\chi _3\right],`$ (21)
which contain the combinations of the couplings

$`c_1^j`$ $`=`$ $`(g_V^j{}_{}{}^{2}+g_A^j{}_{}{}^{2}),`$ (22)
$`c_2^j`$ $`=`$ $`(g_V^j{}_{}{}^{2}g_A^j{}_{}{}^{2}),`$ (23)
$`c_3^j`$ $`=`$ $`2g_V^jg_A^j,`$ (24)
where $`g_V`$ and $`g_A`$ are the vector and axial-vector couplings to the $`Z`$ boson. The propagator factors are given by

$`\chi _1`$ $`=`$ $`{\displaystyle \frac{1}{\mathrm{sin}^2(2\theta _W)}}{\displaystyle \frac{Q^2(Q^2M_Z^2)}{(Q^2M_Z^2)^2+\mathrm{\Gamma }_Z^2M_Z^2}},`$ (25)
$`\chi _2`$ $`=`$ $`{\displaystyle \frac{1}{\mathrm{sin}^2(2\theta _W)}}{\displaystyle \frac{Q^2}{Q^2M_Z^2}}\chi _1,`$ (26)
$`\chi _3`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Gamma }_ZM_Z}{Q^2M_Z^2}}\chi _1.`$ (27)
The above cross section is expressed in the following frame. A normalized timelike vector $`\widehat{t}`$ is defined by $`q`$ (the vector boson momentum) and a normalized spacelike vector $`\widehat{z}`$ is defined by $`\stackrel{~}{P}^\mu `$ = $`P^\mu (Pq/q^2)q^\mu `$ for one of the outgoing momenta, for which we choose $`P_2`$ ($`P_i`$ are the hadron momenta),

$`\widehat{t}^\mu `$ $``$ $`{\displaystyle \frac{q^\mu }{Q}},`$ (28)
$`\widehat{z}^\mu `$ $``$ $`{\displaystyle \frac{Q}{P_2q}}\stackrel{~}{P}_2^\mu =2{\displaystyle \frac{P_2^\mu }{z_2Q}}{\displaystyle \frac{q^\mu }{Q}},`$ (29)
The azimuthal angles lie inside the plane orthogonal to $`\widehat{t}`$ and $`\widehat{z}`$. In particular, $`\varphi ^{\mathrm{}}`$ gives the orientation of $`\widehat{l}_{}^\mu \left(g^{\mu \nu }\widehat{t}^\mu \widehat{t}^\nu +\widehat{z}^\mu \widehat{z}^\nu \right)l_\nu `$, the perpendicular part of the lepton momentum $`l`$.

In the cross sections we also encounter the following functions of $`y=l^{}/q^{}`$, which in the lepton center of mass frame equals $`y=(1+\mathrm{cos}\theta _2)/2`$, where $`\theta _2`$ is the angle of $`\widehat{z}`$ with respect to the momentum of the incoming lepton $`l`$:

$`A(y)`$ $`=`$ $`\left({\displaystyle \frac{1}{2}}y+y^2\right)\stackrel{cm}{=}{\displaystyle \frac{1}{4}}\left(1+\mathrm{cos}^2\theta _2\right),`$ (30)
$`B(y)`$ $`=`$ $`y(1y)\stackrel{cm}{=}{\displaystyle \frac{1}{4}}\mathrm{sin}^2\theta _2,`$ (31)
$`C(y)`$ $`=`$ $`(12y)\stackrel{cm}{=}\mathrm{cos}\theta _2.`$ (32)
Since we have chosen $`P_2`$ to define the longitudinal direction, the momentum $`P_1`$ can be used to express the directions orthogonal to $`\widehat{t}`$ and $`\widehat{z}`$. One obtains $`P_1^\mu `$ = $`g_{}^{\mu \nu }P_{1\nu }`$ (see Fig. 1),

$$P_1^\mu =z_1q_T^\mu =z_1Q_T\widehat{h}^\mu ,$$
(33)
where we define the normalized vector $`\widehat{h}^\mu `$ = $`P_1^\mu /|𝑷_1|`$. The angle $`\varphi _1`$ is between $`\widehat{𝒉}`$ and $`\widehat{𝒍}_{}`$.

The asymmetry $`A(𝒒_T)`$ at the $`Z`$ mass (ignoring the photon and interference contributions, which can be easily included) is now defined by:

$$\frac{d\sigma (e^+e^{}h_1h_2X)}{d\mathrm{\Omega }dz_1dz_2d^2𝒒_T}\left\{1+\mathrm{cos}(2\varphi _1)A(𝒒_T)\right\},$$
(34)
with

$`A(𝒒_T)`$ $``$ $`{\displaystyle \frac{_ac_1^ec_2^aB(y)[(2𝒒_T𝒑_T𝒒_T𝒌_T𝒒_T^2𝒑_T𝒌_T)H_1^{}\overline{H}{}_{1}{}^{}]}{Q_T^2M_1M_2_b\left(c_1^ec_1^bA(y){\scriptscriptstyle \frac{1}{2}}c_3^ec_3^bC(y)\right)\left[D_1\overline{D}_1\right]}}.`$ (35)
A confirmation of this asymmetry would confirm the Collins effect without the need of polarization. To avoid repeating irrelevant factors, we will first focus on the numerator and denominator of the term

$$\kappa =\frac{\left[\left(2𝒒_T𝒑_T𝒒_T𝒌_T𝒒_T^2𝒑_T𝒌_T\right)H_1^{}\overline{H}_1^{}\right]}{Q_T^2M_1M_2\left[D_1\overline{D}_1\right]}.$$
(36)
In Refs.  tree level was considered, which means the expressions are valid in the region where the observed transverse momentum is very small compared to the hard scale(s), applicable only in the region of intrinsic transverse momentum. Eq. (36) is a tree level expression if one uses the convolution notation

$$\left[D\overline{D}\right]d^2𝒑_Td^2𝒌_T\delta ^2(𝒑_T+𝒌_T𝒒_T)D^a(z_1,z_1^2𝒑_T^2)\overline{D}{}_{}{}^{a}(z_2,z_2^2𝒌_T^2),$$
(37)
which only involves intrinsic transverse momenta. In order to go beyond this region, we include the Sudakov factor arising from resummed perturbative corrections to the transverse momentum distribution, by considering the approach of Ref. . This will extend the range of applicability from the region of intrinsic transverse momentum to the region of moderate $`𝒒_T`$ values (still under the restriction that $`𝒒_T^2Q^2`$). This is outlined for the general case in the previous section; here we will apply it to the $`\mathrm{cos}(2\varphi )`$ asymmetry.

Resummation of soft gluons into Sudakov form factors (see e.g. ) results in a replacement in Eq. (37) of

$$\delta ^2(𝒑_T+𝒌_T𝒒_T)\frac{d^2𝒃}{(2\pi )^2}e^{i𝒃(𝒑_T+𝒌_T𝒒_T)}e^{S(b)},$$
(38)
leading to (suppressing the flavor indices)

$`\left[D\overline{D}\right]`$ $``$ $`{\displaystyle \frac{d^2𝒃}{(2\pi )^2}e^{i𝒃𝒒_T}e^{S(b)}\stackrel{~}{D}(z_1,b)\stackrel{~}{\overline{D}}(z_2,b)}`$ (39)
$`=`$ $`{\displaystyle \frac{1}{2\pi }}{\displaystyle _0^{\mathrm{}}}𝑑bbJ_0(bQ_T)e^{S(b)}\stackrel{~}{D}(z_1,b)\stackrel{~}{\overline{D}}(z_2,b).`$ (40)
The function $`\stackrel{~}{D}`$ denotes the Fourier transform of $`D`$. In order to compute the above expression, we will assume a Gaussian transverse momentum dependence for $`D_1(z,z^2𝒌_T^2)`$:

$$D_1(z,z^2𝒌_T^2)=D_1(z)R_u^2\mathrm{exp}(R_u^2𝒌_T^2)/\pi z^2D_1(z)𝒢(|𝒌_T|;R_u)/z^2.$$
(41)
For details see Ref. . Taking the Fourier transform of Eq. (41) yields

$$\stackrel{~}{D}{}_{1}{}^{}(z,b^2)=D_1(z)\mathrm{exp}(\frac{b^2}{4R_u^2})/z^2.$$
(42)
The numerator in Eq. (36),

$`\left[\left(2𝒒_T𝒑_T𝒒_T𝒌_T𝒒_T^2𝒑_T𝒌_T\right)D\overline{D}\right]{\displaystyle \frac{d^2𝒃}{(2\pi )^2}e^{i𝒃𝒒_T}e^{S(b)}}`$ (44)
$`\times {\displaystyle }d^2𝒑_Td^2𝒌_T(2𝒒_T𝒑_T𝒒_T𝒌_T𝒒_T^2𝒑_T𝒌_T)e^{i𝒃(𝒑_T+𝒌_T)}D(z_1,z_1^2𝒑_T^2)\overline{D}(z_2,z_2^2𝒌_T^2),`$
cannot be treated exactly like the denominator. A model for the transverse momentum dependence of the function $`H_1^{}`$ is needed. In order to get a first estimate, or rather an upper bound for the asymmetry, one might be inclined to assume the maximally allowed function, by saturating the bound satisfied by $`H_1^{}`$

$$|𝒌_T|\left|H_1^{}(z,|𝒌_T|)\right|zM_hD_1(z,|𝒌_T|),$$
(45)
producing a $`1/|𝒌_T|`$ behavior of $`H_1^{}(z,𝒌_T)`$. However, this is not consistent with the fact that the Collins effect should vanish in the limit $`𝒌_T0`$.

A model by Collins  suggests the following transverse momentum dependence. Collins’ parameterization for the fragmentation function $`H_1^{}`$ is (note that Collins uses the function $`\mathrm{\Delta }\widehat{D}_{H/a}ϵ_T^{ij}s_{Ti}k_{Tj}H_1^{}`$, where $`s_T`$ is the transverse spin of the fragmenting quark)

$$\frac{H_1^{}(z,𝒌_T^2)}{D_1(z,𝒌_T^2)}=\frac{2M_CM_h}{𝒌_T^2+M_C^2}\text{Im}\left[A^{}(k^2)B(k^2)\right]\frac{(1z)}{z},$$
(46)
where $`M_h`$ is the mass of the produced hadron and $`M_C`$ is the quark mass that appears in a dressed fermion propagator $`i(A(k^2)\overline{)}k+B(k^2)M_C)/(k^2M_C^2)`$, the functions $`A`$ and $`B`$ are unity at $`k^2=M_C^2`$. But for the present purpose, the additional fall-off with $`1/𝒌_T^2`$ on top of the Gaussian fall-off is not needed. We will restrict to a Gaussian fall-off (with unknown magnitude) and assume the simple form for $`H_1^a(z,z^2𝒌_T^2)=H_1^a(z)𝒢(|𝒌_T|;R)/z^2`$. Here one should take the radius $`R`$ to be larger than $`R_u`$ of the unpolarized function $`D_1`$, such as to satisfy the bound Eq. (45) for all $`|𝒌_T|`$. Also, we will assume that the fragmentation functions for both hadrons are Gaussians of equal width, i.e. we take $`R_1=R_2=R`$ and $`R_{u1}=R_{u2}=R_u`$ and also $`M_1=M_2=M`$.

One finds

$`{\displaystyle d^2𝒑_Td^2𝒌_T\left(2𝒒_T𝒑_T𝒒_T𝒌_T𝒒_T^2𝒑_T𝒌_T\right)e^{i𝒃(𝒑_T+𝒌_T)}𝒢(𝒑_T^2;R)𝒢(𝒌_T^2;R)}=`$ (48)
$`{\displaystyle \frac{1}{4R^4}}\left[2(𝒒_T𝒃)^2𝒒_T^2𝒃^2\right]\mathrm{exp}\left({\displaystyle \frac{b^2}{2R^2}}\right),`$
which after application to Eq. (44) yields

$$\left[\left(2𝒒_T𝒑_T𝒒_T𝒌_T𝒒_T^2𝒑_T𝒌_T\right)D\overline{D}\right]=\frac{d^2𝒃}{(2\pi )^2}e^{i𝒃𝒒_T}(\frac{1}{4R^4})\left[2(𝒒_T𝒃)^2𝒒_T^2𝒃^2\right]e^{S(b)}\stackrel{~}{D}(z_1,b)\stackrel{~}{\overline{D}}(z_2,b).$$
(49)
It is important to note that the factor $`1/R^4`$ stems from the intrinsic transverse momentum of the functions $`H_1^{}(z,𝒌_T)`$ and therefore, is not dependent on the scale $`Q`$ (see discussion after Eq. (12)).

In analogy to $`A_{TT}^V`$ , Eq. (36) can be transformed into (keeping in mind that in the end the summation over $`a,b`$ in numerator and denominator should be performed separately)

$$\kappa =\frac{H_1^a(z_1)\overline{H}{}_{1}{}^{a}(z_2)}{4M^4R^4D_1^b(z_1)\overline{D}{}_{1}{}^{b}(z_2)}𝒜(Q_T)=\frac{H_1^{(1)a}(z_1)\overline{H}{}_{1}{}^{(1)a}(z_2)}{D_1^b(z_1)\overline{D}{}_{1}{}^{b}(z_2)}𝒜(Q_T),$$
(50)
where for future reference we have also expressed the asymmetry in terms of the function $`H_1^{(1)}(z)=z^2d^2𝒌_T𝒌_T^2/(2M^2)H_1^{}(z,z^2𝒌_T^2)`$ which occurs frequently in $`Q_T`$-weighted cross sections (e.g. ); furthermore,

$$𝒜(Q_T)M^2\frac{_0^{\mathrm{}}𝑑bb^3J_2(bQ_T)\mathrm{exp}\left(S(b_{})S_{NP}(b)\right)}{_0^{\mathrm{}}𝑑bbJ_0(bQ_T)\mathrm{exp}\left(S(b_{})S_{NP}(b)\right)}.$$
(51)
This is similar to the factor $`𝒜(Q_T)`$ of Ref. , but with the replacement of $`J_0(bQ_T)J_2(bQ_T)`$ in the numerator. We note that unlike the case of $`A_{TT}^V`$ investigated in Ref. , the present asymmetry does not need to oscillate as a function of $`Q_T`$. Rather the vanishing of the numerator after integration over $`d^2𝒒_T`$ is due to the angular integration. Note that $`\kappa `$ only depends on $`Q_T`$ and not on $`\varphi _1`$. Also, we note that this asymmetry has a kinematic zero at $`Q_T=0`$, since $`\widehat{𝒉}`$ cannot be defined in that case (and indeed $`J_2(Q_T=0)=0`$). This is also seen from Ref. , where a similar asymmetry factor has been investigated as a function of $`Q_T`$.

Here the main focus will be on the factor $`𝒜(Q_T)`$, which is a measure for the effect of the Sudakov factors on the asymmetry compared to the tree level result $`\mathrm{exp}\left[(R^2R_u^2)Q_T^2/2\right]M^2Q_T^2R^6/R_u^2`$ (cf. Eq. (53)), which is valid only for values of $`Q_T`$ of the order of the intrinsic transverse momentum.

In the above expression we have introduced the usual cut-off $`b_{\mathrm{max}}`$, via $`b_{}=b/\sqrt{1+b^2/b_{\mathrm{max}}^2}`$ and replaced $`\frac{1}{2}b^2/R^2S_{NP}(b)`$, for which we take one of the standard nonperturbative smearing functions, needed to describe the low $`𝒒_T`$ region properly. It is important to realize that $`S_{NP}(b)`$ is introduced only in part to take care of the smearing due to the intrinsic transverse momentum (cf. Eq. (12)), hence one cannot simply equate $`\frac{1}{2}b^2/R^2`$ with $`S_{NP}(b)`$. But taking into account the term $`\mathrm{exp}\left(\frac{1}{2}b^2/R^2\right)`$ in addition to $`S_{NP}(b)`$ will just produce a change in the coefficient of the $`b^2`$ term in $`S_{NP}(b)`$. To keep the unpolarized cross section unaffected, we will therefore introduce as nonperturbative term $`\mathrm{exp}\left[S_{NP}(b)+\frac{1}{2}b^2/R^2\right]`$, in order not to count the contribution from intrinsic transverse momentum twice. It is also worth mentioning again that since $`RR_u`$ (important at tree level), $`S_{NP}(b)`$ need not be the same in numerator and denominator (much less relevant however, since it affects only $`g_A,g_B`$ and not $`g_Q`$, cf. Eq. (12)) and that in principle, it can depend on $`z_1`$ and $`z_2`$, but we will not take into account these refinements.

Here we will take for the nonperturbative Sudakov factor the parameterization of Ladinsky-Yuan (Ref.  with $`x_1x_2=10^2`$),

$$S_{NP}(b)=g_1b^2+g_2b^2\mathrm{ln}\left(\frac{Q}{2Q_0}\right),$$
(52)
with $`g_1=0.11\text{GeV}^2,g_2=0.58\text{GeV}^2,Q_0=1.6\text{GeV}`$ and $`b_{\mathrm{max}}=0.5\text{GeV}^1`$. Part of the results will depend considerably on this choice and this issue will be addressed in detail below. The value $`g_1=0.11\text{GeV}^2`$ can be viewed as an intrinsic transverse momentum of $`𝒑_T^2=1/R^2=(220\text{MeV})^2`$.

The reason we have chosen the parameterization of Ladinsky-Yuan , which is fitted to the transverse momentum distribution of $`W/Z`$ production in $`pp`$ $`(p\overline{p})`$ scattering, is that unfortunately there is no nonperturbative Sudakov factor available for $`e^+e^{}A+B+X`$, except for the energy-energy correlation function obtain from low energy data . This is surprising considering the wealth of data from LEP experiments. Moreover, the nonperturbative Sudakov factor in the energy-energy correlation function ($`_{A,B}𝑑z_1z_1𝑑z_2z_2Q^2𝑑\sigma /𝑑Q_T^2`$) as fitted in Ref.  is not the same as the one in the differential cross section as a function of the lightcone momentum fractions. For related discussions, also relevant for SIDIS, see Refs. .

Although the nonperturbative Sudakov factors for Drell-Yan and $`e^+e^{}`$ need not be related, the $`S_{NP}`$ of Ref.  can be viewed as a generic one and allows us to study the general features of the Sudakov suppression. We will also investigate the dependence on the nonperturbative Sudakov factor by varying the parameters. At a later stage one can always insert a more appropriate (phenomenologically determined) nonperturbative Sudakov factor into the asymmetry expression Eq. (51). Our results also underline the importance of a good determination of $`S_{NP}`$.

In Fig. 2 the asymmetry factor $`𝒜(Q_T)`$ is given at the scales $`Q=30\text{GeV}`$, $`Q=60\text{GeV}`$ and $`Q=M_Z`$.

The maximum of the asymmetry factor $`𝒜(Q_T)`$ at $`Q=30\text{GeV}`$ ($`S_{NP}(b)=1.41b^2`$), $`Q=60\text{GeV}`$ ($`S_{NP}(b)=1.81b^2`$) and $`Q=M_Z`$ ($`S_{NP}(b)=2.05b^2`$) is seen to be $`0.57`$ (at $`Q_T3.6\text{GeV}`$), $`0.31`$ (at $`Q_T3.8\text{GeV}`$) and $`0.22`$ (at $`Q_T4\text{GeV}`$), respectively. One observes that the magnitude of the asymmetry factor goes down with increasing energy and the position of the maximum –and also the average $`Q_T`$– shifts to higher values of $`Q_T`$.

Now we will discuss the dependence of these results on the choice of $`b_{\mathrm{max}}`$ and $`S_{NP}`$. Taking a higher value of $`b_{\mathrm{max}}`$ increases the Gaussian width. The above choice of $`b_{\mathrm{max}}=0.5\text{GeV}^1`$ can be considered as optimistic ($`1/b_{\mathrm{max}}`$ is the scale down to which one trusts perturbation theory), enhancing the asymmetry factor somewhat.

The asymmetry factor decreases with increasing Gaussian smearing width in $`S_{NP}`$, to which it has a considerable sensitivity. Empirically, we find that if the Gaussian width is reduced by a factor $`\alpha `$, then the maximum of the asymmetry increases roughly by $`\alpha `$ and the corresponding value of $`Q_T`$ by $`\sqrt{\alpha }`$. In Fig. 3 this is illustrated for the asymmetry factor $`𝒜(Q_T)`$ at $`Q=M_Z`$. The solid curve is for $`S_{NP}(b)=2.05b^2`$ and the dashed curve for $`S_{NP}(b)=1.37b^2`$, where the latter width is taken from a recent two parameter fit .

The decrease with energy of the maximum of the asymmetry was found not to be very sensitive to changes to the Gaussian width in $`S_{NP}`$. We find that the decrease goes as $`Q^{0.9}Q^{1.0}`$.

So far we have only considered the asymmetry as a function of $`Q_T`$, which means the cross section needs to be kept differential in the angle $`\varphi _1`$ and the magnitude $`Q_T`$ of the transverse momentum $`𝒒_T`$. Since the asymmetry does not vanish after integration over $`Q_T^2`$ and since the $`Q_T^2`$-integrated cross section has been studied using LEP data , we will now consider that case. Note that the integration has to be done in numerator and denominator separately. Since the denominator has no dependence on $`\varphi _1`$ (the $`\mathrm{cos}(2\varphi _1)`$ dependence belongs to the numerator) one can in fact integrate over $`𝒒_T`$ completely. This corresponds to the inclusive cross section in which soft gluon contributions cancel, hence this integration over the denominator of the asymmetry factor $`𝒜(Q_T)`$ yields 1 and has no dependence on Sudakov factors.

The $`Q_T^2`$-integrated result should be compared to tree level, therefore we will first give the tree level expression. The tree level asymmetry is (neglecting the small $`c_3^a`$ dependent term in the denominator)

$$A^{(0)}(Q_T)=\frac{Q_T^2R^2\mathrm{exp}(R^2Q_T^2/2)}{4M^2R_u^2\mathrm{exp}(R_u^2Q_T^2/2)}\frac{\mathrm{sin}^2\theta _2}{1+\mathrm{cos}^2\theta _2}\frac{_ac_2^aH_1^a(z_1)\overline{H}{}_{1}{}^{a}(z_2)}{_bc_1^bD_1^b(z_1)\overline{D}{}_{1}{}^{b}(z_2)}.$$
(53)
Note that if one integrates the cross section over $`Q_T^2`$, one should be careful to retain all dependences on $`Q_T`$ in both numerator and denominator separately. From the above expression we also infer that $`𝒜(Q_T)`$ itself should be compared with the tree level quantity $`\mathrm{exp}\left[(R^2R_u^2)Q_T^2/2\right]M^2Q_T^2R^6/R_u^2`$ as mentioned before.

In Fig. 4 we have displayed the comparison of $`𝒜(Q_T)`$ at $`Q=90\text{GeV}`$ and the tree level quantity using the values $`R_u^2=1\text{GeV}^2`$ and $`R^2/R_u^2=3/2`$, which were chosen such as to minimize the magnitude. The value $`R_u^2=1\text{GeV}^2`$ can be regarded as too small already. We conclude that inclusion of Sudakov factors has the effect of suppressing the tree level result by at least an order of magnitude. This is important to keep in mind when making predictions of transverse momentum dependent azimuthal spin asymmetries based on tree level expressions.

For the particular case of a produced $`\pi ^+`$ and $`\pi ^{}`$, we will make the following simplifying assumptions. We assume $`D_1^{u\pi ^+}(z)=D_1^{\overline{d}\pi ^+}(z)`$, $`D_1^{d\pi ^{}}(z)=D_1^{\overline{u}\pi ^{}}(z)`$ and neglect unfavored fragmentation functions like $`D^{d\pi ^+}(z)`$; and similarly for the Collins functions. As a consequence of these assumptions the fragmentation functions can be taken outside the flavor summation and appear as a square (at the average $`z_1,z_2`$). After $`Q_T^2`$ integration one arrives at a cross section differential in the angle $`\varphi _1`$ of $`𝒒_T`$:

$$\frac{d\sigma }{d\mathrm{\Omega }d\varphi _1}\left(1+A\left[\frac{H_1^{}}{D_1}\right]^2F(y)\mathrm{cos}(2\varphi _1)\right),$$
(54)
where

$$F(y)=\frac{_{a=u,\overline{d}}c_1^ec_2^aB(y)}{_{b=u,\overline{d}}\left(c_1^ec_1^bA(y){\scriptscriptstyle \frac{1}{2}}c_3^ec_3^bC(y)\right)}\frac{\mathrm{sin}^2\theta _2}{1+\mathrm{cos}^2\theta _2}\frac{_a(g_V^a{}_{}{}^{2}g_A^a{}_{}{}^{2})}{_b(g_V^b{}_{}{}^{2}+g_A^b{}_{}{}^{2})},$$
(55)
which is largest at $`\theta _2=90^{}`$: $`F_{\mathrm{max}}0.5`$ (for a plot of the full factor $`F(\theta _2)`$ see Ref. ).

For the prefactor $`A`$ one finds at tree level $`A^{(0)}=1/(2M^2R^2)=𝒑_T^2/(2M^2)`$ and one should assume a typical intrinsic transverse momentum squared value, in the range of $`𝒑_T^2(200700\text{MeV})^2`$. For pions this means $`A^{(0)}112`$. This can be compared to $`A^{(0)}=6/\pi `$ of Ref. , which yields $`p_T^2(270\text{MeV})^2`$ (or $`R^213\text{GeV}^2`$), which seems to be a reasonable value (remember that $`R^2`$ must be larger than $`R_u^2`$).

After including the Sudakov factors our numerical calculation yields $`A=0.07`$ ($`S_{NP}(b,Q=M_Z)=2.05b^2`$), much smaller than the tree level values discussed above. The Gaussian width in $`S_{NP}`$ is not crucial for this conclusion; if one takes a much smaller width of for instance $`1\text{GeV}^2`$, then one obtains $`A=0.11`$. It has to be emphasized that the width used in $`S_{NP}(b,Q)`$ should increase with increasing energy.

Our result shows that upon including Sudakov factors one retrieves parton model characteristics (also noted in Ref. ), but with transverse momentum spreads that are significantly larger than would be expected from intrinsic transverse momentum (this is supported by the presently available parameterizations of $`S_{NP}`$ in various processes, see e.g. ). This means that a large tree level value for $`A^{(0)}`$ (e.g. $`6/\pi `$), would lead to an extracted Collins function that is considerably smaller than if Sudakov factors would have been included. After the nonperturbative Sudakov factor has been obtained from LEP data, one can estimate exactly how much. For our choice of $`S_{NP}(b,Q=M_Z)=2.05b^2`$, the effect on the magnitude of the Collins function is an additional factor $`\sqrt{6/(0.07\pi )}5`$ compared to the extraction using the tree level expression.

We conclude that the Sudakov factors produce a strong suppression compared to the tree level result. A tree level analysis applied at $`Q=M_Z`$ is expected to overestimate the asymmetry at least by an order of magnitude and therefore, it will underestimate the Collins function significantly. Hence, Sudakov factors should be taken into account when extracting the Collins function from LEP data or in general, from $`e^+e^{}`$ data obtained at high values of $`\sqrt{s}`$.

## IV Collins asymmetry in SIDIS

The Collins function $`H_1^{}`$ originally was shown to lead to a single spin azimuthal asymmetry in semi-inclusive DIS . This asymmetry has received much attention, since it would provide an additional way of accessing the transversity function $`h_1`$. A preliminary determination of the Collins asymmetry from SMC data has been performed , yielding an asymmetry for $`\pi ^+`$ production of $`11\%\pm 6\%`$. Also, the $`\mathrm{sin}\varphi `$ azimuthal asymmetry as recently measured by the HERMES Collaboration  might be related to the Collins asymmetry, providing further indication that the Collins function is nonzero. Measurements by HERMES and COMPASS are expected to provide more conclusive information on the Collins effect and its magnitude. Since these experiments are performed at different energies, it is important to know the $`Q^2`$ dependence of the asymmetry. Here we will investigate the $`Q^2`$ dependence of the transverse momentum ($`Q_T`$) distribution of the asymmetry.

From Ref.  we extract the expression for this single spin asymmetry (extended to include contributions from $`Z`$-boson exchange; the expressions for $`W`$-boson exchange can be found in Ref. ):

$$\frac{d\sigma (\mathrm{}H\mathrm{}^{}hX)}{dxdzdyd\varphi _{\mathrm{}}d^{\mathrm{\hspace{0.17em}2}}𝒒_T}=\frac{2\alpha ^2xz^2s}{Q^4}\underset{a,\overline{a}}{}\left\{K_1^a(y)\left[f_1D_1\right]\right|𝑺_T|K_3^a(y)\mathrm{sin}(\varphi _h^{\mathrm{}}+\varphi _S^{\mathrm{}})\left[\widehat{𝒉}𝒌_T\frac{h_1H_1^{}}{M_h}\right]+\mathrm{}\},$$
(56)
where $`𝑺_T`$ is the transverse spin of the incoming hadron $`H`$. The couplings $`K_1^a(y)`$ and $`K_3^a(y)`$ are of the same form as before, except that now

$`A(y)`$ $`=`$ $`\left(1y+{\displaystyle \frac{1}{2}}y^2\right),`$ (57)
$`B(y)`$ $`=`$ $`(1y),`$ (58)
$`C(y)`$ $`=`$ $`y(2y),`$ (59)
and $`y=(Pq)/(Pl)q^{}/l^{}`$ ($`l`$ is the momentum of the beam lepton). Also, since $`Q^2`$ is now space-like, the width $`\mathrm{\Gamma }_Z`$ can be ignored:

$`\chi _1`$ $`=`$ $`{\displaystyle \frac{1}{\mathrm{sin}^2(2\theta _W)}}{\displaystyle \frac{Q^2}{Q^2+M_Z^2}},`$ (60)
$`\chi _2`$ $`=`$ $`\left(\chi _1\right)^2.`$ (61)
The azimuthal angle $`\varphi _h^{\mathrm{}}`$ ($`\varphi _S^{\mathrm{}}`$) around the three-momentum of the virtual boson is between the lepton’s three-momentum and the outgoing hadron’s three-momentum (transverse spin), cf. Ref. .

If we write the cross section as

$$\frac{d\sigma (\mathrm{}H\mathrm{}^{}hX)}{dxdzdyd\varphi _{\mathrm{}}d^{\mathrm{\hspace{0.17em}2}}𝒒_T}\left\{1\right|𝑺_T|\mathrm{sin}(\varphi _h^{\mathrm{}}+\varphi _S^{\mathrm{}})A(𝒒_T)\},$$
(62)
the asymmetry analyzing power is given by

$$A(𝒒_T)\frac{_aK_3^a(y)\left[𝒒_T𝒌_Th_1H_1^{}\right]}{Q_TM_h_aK_1^a(y)\left[f_1D_1\right]}.$$
(63)
Assuming again Gaussian transverse momentum dependence and including Sudakov factors one arrives at

$$A(𝒒_T)=\frac{_aK_3^a(y)h_1^a(x)H_1^a(z)}{2M_h^2R^2_bK_1^b(y)f_1^b(x)D_1^b(z)}𝒜(Q_T),$$
(64)
where $`R^2`$ is the Gaussian width of $`H_1^{}`$ and the asymmetry factor $`𝒜(Q_T)`$ is defined as

$$𝒜(Q_T)M_h\frac{_0^{\mathrm{}}𝑑bb^2J_1(bQ_T)\mathrm{exp}\left(S(b_{})S_{NP}(b)\right)}{_0^{\mathrm{}}𝑑bbJ_0(bQ_T)\mathrm{exp}\left(S(b_{})S_{NP}(b)\right)}.$$
(65)
This should be compared with Eq. (51) of the previous section. In Fig. 5 the asymmetry factor $`𝒜(Q_T)`$ is given at the scales $`Q=30\text{GeV}`$, $`Q=60\text{GeV}`$ and $`Q=M_Z`$.

The maximum of the asymmetry factor $`𝒜(Q_T)`$ at $`Q=30\text{GeV}`$, $`Q=60\text{GeV}`$ and $`Q=M_Z`$ is seen to be $`0.60`$ (at $`Q_T3.1\text{GeV}`$), $`0.42`$ (at $`Q_T3.4\text{GeV}`$) and $`0.34`$ (at $`Q_T3.5\text{GeV}`$), respectively. Again we note that the magnitude of the asymmetry factor goes down with increasing energy and the position of the maximum shifts to higher values of $`Q_T`$.

As before we have used the nonperturbative Sudakov factor by Ladinsky and Yuan, Eq. (52), in the absence of a well-established $`S_{NP}`$ for SIDIS (cf. Ref. ). But the results obtained here can nevertheless be viewed as generic, only changes in the specific numbers are expected. Like for the previous $`\mathrm{cos}(2\varphi )`$ asymmetry, the asymmetry factor decreases with increasing Gaussian smearing width. Empirically we find that if the Gaussian width in $`S_{NP}`$ is reduced by a factor $`\alpha `$, then both the maximum of the asymmetry factor and the corresponding value of $`Q_T`$ increase roughly by $`\sqrt{\alpha }`$. However, the decrease with energy of the maximum of the asymmetry is not very sensitive to changes in this Gaussian width. We find that the decrease goes as $`Q^{0.5}Q^{0.6}`$, which is considerably slower than for the $`\mathrm{cos}(2\varphi )`$ asymmetry. In the next section we will discuss the comparison between different asymmetries in more detail.

We will now compare these results to the tree level expression for this Collins asymmetry. At tree level the convolution $``$ (where $`w`$ denotes a weight function) is given by

$$\left[w(𝒑_T,𝒌_T)fD\right]d^{\mathrm{\hspace{0.17em}2}}𝒑_Td^{\mathrm{\hspace{0.17em}2}}𝒌_T\delta ^2(𝒑_T+𝒒_T𝒌_T)w(𝒑_T,𝒌_T)f^a(x,𝒑_T^2)D^a(z,z^2𝒌_T^2).$$
(66)
If Gaussian transverse momentum distribution and fragmentation functions are assumed, one obtains at tree level

$$A^{(0)}(Q_T)=\frac{Q_TR^2\mathrm{exp}(R^2Q_T^2/2)}{2M_hR_u^2\mathrm{exp}(R_u^2Q_T^2/2)}\frac{_aK_3^a(y)h_1^a(x)H_1^a(z)}{_bK_1^b(y)f_1^b(x)D_1^b(z)}.$$
(67)
Therefore, one should compare $`𝒜(Q_T)`$ with the tree level quantity $`\mathrm{exp}\left[(R^2R_u^2)Q_T^2/2\right]M_hQ_TR^4/R_u^2`$.

In Fig. 6 we have displayed the comparison of $`𝒜(Q_T)`$ at $`Q=90\text{GeV}`$ and the tree level quantity using again the values $`R_u^2=1\text{GeV}^2`$ and $`R^2/R_u^2=3/2`$, which minimize the maximum asymmetry value. As before we conclude that inclusion of Sudakov factors has the effect of suppressing the tree level result. It is also clear that for the Collins asymmetry (which arises from a single Collins effect) the Sudakov suppression is less severe than for the case of the $`\mathrm{cos}(2\varphi )`$ asymmetry (which depends on the Collins effect squared). One also observes a pronounced increase of the average $`Q_T`$ upon inclusion of Sudakov factors. The difference to tree level becomes even more pronounced as the choice of $`R_u^2`$ is increased to larger values more appropriate for a tree level analysis. Hence, also for transverse momentum dependent azimuthal spin asymmetries in SIDIS one needs to include Sudakov factors in order to extract distribution and fragmentation functions reliably.

We now turn to a more general discussion of the different types of transverse momentum dependent azimuthal spin asymmetries, categorized by different transverse momentum weights and angular dependences. The comparisons are most cleanly done for the asymptotic $`Q^2`$ case.

## V Asymptotic behavior

In the previous sections we have studied asymmetries with different types of weights, namely $`2\widehat{𝒉}𝒑_T\widehat{𝒉}𝒌_T𝒑_T𝒌_T`$ and $`\widehat{𝒉}𝒌_T`$. The former (and also $`𝒑_T𝒌_T`$ in the asymmetry $`A_{TT}^V`$ of Ref. ) is typical of a double transverse spin asymmetry. It has two powers of transverse momentum and is therefore expected to decrease faster than the (single spin) asymmetries which have weights with a single power. The reason for this difference in decrease is that the asymmetries, which are convolutions of transverse momentum distributions, are largest if these distributions are large and have a large overlap. The distributions are largest at small $`𝒑_T`$ and $`𝒌_T`$ and the overlap is largest at small transverse momentum $`𝒒_T`$. However, the explicit powers of transverse momentum in the weights tend to suppress the region of small transverse momentum. More powers of transverse momentum in the weight thus implies more suppression. The power $`n`$ of $`b^n`$ in the numerator compared to the denominator is a measure of this effect. A weight proportional to $`𝒑_T^i𝒌_T^j`$ will lead to $`n=2`$ (cf. Eq. (51)), whereas $`𝒑_T^i`$ or $`𝒌_T^j`$ will lead to $`n=1`$ (cf. Eq. (65)). Larger $`n`$ means, larger suppression and since the transverse momentum distributions broaden and decrease in magnitude with increasing energy, the convolution also decreases. Asymmetries with $`n=1`$ will fall off slower with energy and the Sudakov suppression is less severe. This natural expectation is clearly observed in the numerical analysis of the examples investigated here and in Ref. .

It is interesting to see what is the asymptotic behavior of the Sudakov suppression. As is usually done, we will assume that at high values of $`Q`$ the nonperturbative factor is irrelevant (true only asymptotically) and consider asymmetry factors of the form

$$𝒜_{n,m}(Q_T)M^n\frac{_0^{\mathrm{}}𝑑b^2b^nJ_m(bQ_T)\mathrm{exp}\left(S(b)\right)}{_0^{\mathrm{}}𝑑b^2J_0(bQ_T)\mathrm{exp}\left(S(b)\right)}.$$
(68)
The asymptotic behavior of the denominator at $`Q_T=0`$ (relevant for the unpolarized differential cross section) has been studied before  by means of a saddle point approximation. Unfortunately, such a saddle point approximation cannot be applied to the ratios $`𝒜_{n,m}(Q_T)`$ for $`m0`$ and $`Q_T0`$, since one cannot take $`J_m(bQ_T)J_m(b_{SP}Q_T)`$ to deduce the $`Q^2`$ dependence. For instance, since $`b_{SP}Q_T`$ is in general not small, i.e. $`m`$ (even under the assumption $`Q_T^2Q^2`$), no approximation to the Bessel functions can be used in the region of interest (i.e., around the value of $`Q_T`$ where the asymmetry is maximal).

In Ref.  the asymmetry factor $`𝒜_{2,0}(Q_T)`$ arises, which does allow for a straightforward saddle point approximation at $`Q_T=0`$, where the asymmetry is maximal. Therefore, we simply assume that the weight $`2\widehat{𝒉}𝒑_T\widehat{𝒉}𝒌_T𝒑_T𝒌_T`$ has similar behavior as $`𝒑_T𝒌_T`$. From our numerical studies we conclude that the maximum of $`𝒜_{2,2}(Q_T)`$ actually decreases faster with energy than the maximum of $`𝒜_{2,0}(Q_T)`$ (at $`Q_T=0`$). Moreover, for the asymptotic behavior of $`𝒜_{n,0}(Q_T=0)`$ one can obtain an analytic expression. First we would like to comment on the saddle point approximation of expressions of the form

$$_0^{\mathrm{}}𝑑b^2b^n\mathrm{exp}\left(S(b)\right)=\frac{d\mathrm{ln}(b^2\mathrm{\Lambda }^2)}{(\mathrm{\Lambda }^2)^{1+\frac{n}{2}}}\mathrm{exp}\left\{(1+\frac{n}{2})\mathrm{ln}(b^2\mathrm{\Lambda }^2)S(b)\right\}.$$
(69)
As before we will keep only the leading term in the expansion of $`A(\alpha _s(\mu ))`$ and $`\alpha _s(\mu )/\pi =1/(\beta _1\mathrm{ln}(\mu ^2/\mathrm{\Lambda }^2))`$:

$$S(b,Q)=_{b_0^2/b^2}^{Q^2}\frac{d\mu ^2}{\mu ^2}\frac{C_F}{\beta _1\mathrm{ln}(\mu ^2/\mathrm{\Lambda }^2)}\mathrm{ln}\frac{Q^2}{\mu ^2}.$$
(70)
The integral Eq. (69) has a saddle point at

$$b_{SP}=\frac{b_0}{\mathrm{\Lambda }}\left(\frac{Q}{\mathrm{\Lambda }}\right)^{C_F/[C_F+(1+\frac{n}{2})\beta _1]}.$$
(71)
For $`n=0`$ we retrieve the power $`0.41`$ of Refs. . For general $`n`$ one finds in the saddle point approximation (defining $`\gamma _n=C_F/[C_F+(1+\frac{n}{2})\beta _1]`$)

$$_0^{\mathrm{}}𝑑b^2b^n\mathrm{exp}\left(S(b)\right)\stackrel{SP}{}\left(\frac{b_0^2}{\mathrm{\Lambda }^2}\right)^{(1+\frac{n}{2})}\left(\frac{Q^2}{\mathrm{\Lambda }^2}\right)^{(1+\frac{n}{2})\gamma _n+\frac{C_F}{\beta _1}(1+\gamma _n+\mathrm{ln}(\gamma _n))}.$$
(72)
In this way one finds for the denominator of $`𝒜_{2,0}(Q_T=0)`$ an approximate power behavior of $`(Q^2)^{0.94}`$ and for the numerator $`(Q^2)^{0.62}`$, when one takes $`C_F=4/3`$ and $`\beta _1=23/12`$ (corresponding to 5 flavors, but 6 flavors only makes a few percent difference). The approximated asymmetry factor $`𝒜_{2,0}(Q_T=0)`$ then has an asymptotic power behavior of $`(Q^2)^{0.32}`$, which is a similar decrease as found numerically in Ref.  and apparently holds also for very large values of $`Q`$. At lower energies, the saddle point is not very pronounced and the introduction of the Gaussian smearing and $`b_{}`$ will have an effect on the approximation, hence we view this agreement as a coincidence. Of course, the saddle point approximation should become better with increasing $`Q^2`$.

For the ratios $`𝒜_{n,0}(Q_T=0)`$ one can also derive an analytic expression for the asymptotic $`Q^2`$ dependence:

$$𝒜_{n,0}(Q_T=0)=M^n\frac{_0^{\mathrm{}}𝑑b^2b^n\mathrm{exp}\left(S(b)\right)}{_0^{\mathrm{}}𝑑b^2\mathrm{exp}\left(S(b)\right)}=c_n\left(\frac{M^2b_0^2}{\mathrm{\Lambda }^2}\right)^{\frac{n}{2}}\left(\frac{Q^2}{\mathrm{\Lambda }^2}\right)^{\frac{C_F}{\beta _1}\mathrm{ln}c_n},$$
(73)
where $`c_n=(1+C_F/\beta _1)/(1+n/2+C_F/\beta _1)`$. For $`n=2`$ this yields $`(Q^2)^{0.32}`$ for 5 flavors, which shows that the saddle point approximation is an extremely good approximation (the power of $`Q^2`$ differs by $`𝒪(10^9)`$). One now expects that the maximum value of $`𝒜_{1,1}(Q_T)`$ (appearing in the Collins asymmetry) to fall off at least as fast as $`𝒜_{1,0}(Q_T=0)`$ (for all $`Q`$), which asymptotically goes as $`(Q^2)^{0.18}`$. We note that the introduction of the Gaussian smearing and $`b_{}`$ have a considerable effect on the approximation at lower energies as our numerical studies demonstrate (in the studied range of $`10100\text{GeV}`$ we found approximately $`1/Q`$ for $`\mathrm{max}(𝒜_{2,2}(Q_T))`$ and $`1/\sqrt{Q}`$ for $`\mathrm{max}(𝒜_{1,1}(Q_T))`$). Our main conclusion is that for very high $`Q`$ transverse momentum dependent azimuthal spin asymmetries will fall off as a fractional power of $`1/Q`$ and that the behavior (or an upper bound to it) can approximately be found by looking at the power of $`b`$ in the integrand.

Finally, we want to mention that the increase of the Sudakov suppression with energy provides a solution to the following problem. In general, azimuthal spin asymmetries in collinear configurations, where partons are collinear to the parent hadrons, are suppressed by explicit powers of the hard energy scales. The transverse momentum factors in the weights have to be generated in the hard scattering part and for this one pays a price in terms of inverse powers of the hard scale. The collinear case can be viewed as averaged over the transverse momenta, relevant for $`𝒒_T`$-integrated cross sections, which receives its main contribution from the region of small transverse momentum. On the other hand, in the case of $`𝒒_T`$-dependent differential cross sections at small $`𝒒_T`$ (compared to $`Q`$), small parton transverse momentum needs to be included (one deals with nearly collinear partons therefore). In this case there need not be an explicit power suppression to generate similar azimuthal spin asymmetries (the average transverse momentum is now a scale in the problem, so dimensionless ratios using $`Q_T`$ rather than $`Q`$ can be formed). This seems counterintuitive, since by considering a differential cross section for a nearly collinear configuration, one can obtain information without power suppression, that is power suppressed in the $`𝒒_T`$-integrated cross section, which receives its main contribution from the region of small transverse momentum also (for an example, see Ref. ). But we find that instead of explicit power suppression the presence of Sudakov factors gives rise to partial power suppression. In this sense the power suppression is effectively replaced by a Sudakov suppression. This means that the $`𝒒_T`$-dependent azimuthal spin asymmetries vanish with increasing energy, as do their $`𝒒_T`$-integrated counterparts. Hence, departing from the collinear configuration (that is, by including small transverse momenta) does not allow one to measure azimuthal spin asymmetries in (differential) cross sections at arbitrary high energies.

## VI Conclusions

In this article we have shown by quantitative examples how transverse momentum dependent azimuthal spin asymmetries are suppressed by Sudakov factors, in the region where the transverse momentum is much smaller than the large energy scale $`Q^2`$. Physically, the Sudakov suppression stems from broadening of the transverse momentum distribution due to recoil from soft gluon radiation and the suppression increases with energy. This implies that tree level estimates of transverse momentum dependent azimuthal spin asymmetries tend to overestimate the magnitudes and increasingly so with rising energy.

The size and $`Q^2`$ dependence of the Sudakov suppression have been studied numerically for two such asymmetries, both arising due to the Collins effect. The size of the suppression (compared to tree level) depends considerably on the nonperturbative Sudakov factor that must be determined from experiment; however, the $`Q^2`$ dependence of the suppression turned out to be much less sensitive to the nonperturbative input.

We observe that in general the larger the power $`n`$ of the transverse momentum in the weight of an asymmetry, the larger the suppression. For the Collins effect driven $`\mathrm{cos}(2\varphi )`$ asymmetry in electron-positron annihilation into two almost back-to-back pions ($`n=2`$), the Sudakov suppression was numerically found to be approximately $`1/Q`$ for the maximum of the asymmetry and an upper bound for the asymptotic behavior was found to be $`1/Q^{0.6}`$. For the Collins effect single spin asymmetry in semi-inclusive deep inelastic scattering ($`n=1`$), the Sudakov suppression was numerically found to be approximately $`1/\sqrt{Q}`$ for the maximum of the asymmetry and an upper bound for the asymptotic behavior was found to be $`1/Q^{0.4}`$. For the maximum of the asymmetries of the type $`𝒜_{n,0}(Q_T)`$, the asymptotic $`Q^2`$ dependence could be calculated analytically. This provides upper bounds on the fall-off of $`𝒜_{n,m}(Q_T)`$, with $`m0`$.

Since our results depend on the input for the nonperturbative Sudakov factors $`S_{NP}`$, which is not (well) determined for the above processes, the numerical conclusions about the size and $`Q^2`$ dependence of the suppression should be viewed as generic, not as precise predictions. Therefore, we would like to stress the need for an extraction of the nonperturbative Sudakov factor from the process $`e^+e^{}A+B+X`$ (for any two, almost back-to-back hadrons $`A`$ and $`B`$) and from SIDIS. Considering the wealth of data from the LEP and HERA experiments this should pose no problem.

The Sudakov suppression of the transverse momentum distribution of azimuthal spin asymmetries is significant already for $`Q`$ values in the range of $`10100\text{GeV}`$ as can be concluded from comparison to tree level. In single spin asymmetries the magnitude of the suppression is less severe than in double spin asymmetries and in both cases a pronounced shift of the average $`Q_T`$ to higher values is observed. We conclude that it is essential to take into account Sudakov factors in transverse momentum dependent azimuthal spin asymmetries.

## ACKNOWLEDGMENTS

I would like to thank John Collins, Anatoli Efremov, Zi-wei Lin, George Sterman and Werner Vogelsang for helpful comments and discussions. Also, I thank Chris Dawson for C$`^{_{++}}`$ help. Furthermore, I thank RIKEN, Brookhaven National Laboratory and the U.S. Department of Energy (contract number DE-AC02-98CH10886) for providing the facilities essential for the completion of this work.


format: text

text:
# Searching for fluctuations in the IGM temperature using the Lyman 𝛼 forest

## 1. Introduction

The Lyman $`\alpha `$ forest has become one of the major tools of cosmology. It has been used to set constraints on many of the parameters of the cosmological model, such as the amplitude and slope of the spectrum of primordial density fluctuations and the nature of the dark matter (Croft et al. (1998, 1999); White & Croft (2000); Narayanan et al. (2000); Zaldarriaga, Hui & Tegmark (2000); Croft et al. (2000)). In the popular cosmologies, the forest arises rather naturally and both cosmological simulation (Cen et al. (1994); Hernquist et al. (1995); Zhang et al. (1995); Miralda-Escude et al. (1996); Muecket et al. (1996); Wadsley & Bond (1996); Theuns et al. (1998)) as well as analytical models have been used to understand its properties (Bi et al. (1992); Reisenegger & Miralda-Escude (1995); Bi & Davidsen (1997); Gneidn & Hui (1996); Croft et al. (1997); Hui & Gnedin 1997a ; Hui et al. 1997b ).

The statistical properties of the forest flux are affected by the physical properties of the intergalactic medium (IGM) such as its temperature or the slope of its equation of state affect. In turn these physical properties of the IGM are sensitive to the ionization history of the universe and the cosmological parameters. An intensive effort has been devoted to understanding what different scenarios predict for the forest and in constraining those scenarios with the available observations (eg. Hui & Gnedin 1997a ; Gnedin & Hui (1998); Schaye et al. (1999); McDonald et al. (2000); Croft et al. (2000)).

The temperature of the IGM has been determined in two different ways, by measuring the widths of the lines (Schaye et al. (1999); McDonald et al. (2000)) and by studying the small scale power spectrum of the flux (Zaldarriaga, Hui & Tegmark (2000)). Both methods find values of the temperature around $`T20,000K`$ which are higher than theoretical expectations (ie. Hui & Gnedin 1997a ; Hui et al. 1997b ). Both these type of studies as well as most theoretical models assume that the temperature of the IGM is constant in space. The aim of this paper is to formulate a method for testing this assumption, illustrate it using the spectrum of one quasar and investigate what type of constraints on the temperature fluctuations could be obtained with more data.

A possible explanation for the high temperature of the IGM around redshift $`z3`$ is heat deposited in the gas by the double ionization of He. Depending upon the nature of the source responsible for the reionization of the universe and its clumpiness, HeII reionization can be substantially delayed relative to H reionization and could happen around $`z3`$ (Miralda-Escude & Rees (1994); Madau & Meiksin (1994); Miralda-Escudé et al. (2000)).

A natural consequence of this scenario are fluctuations in the temperature of the IGM around $`z3`$. In fact studies of the HeII Lyman $`\alpha `$ forest show that the mean absorption is increasing rapidly with redshift for redshifts between 2 and 3. Moreover several gaps of size $`1000\mathrm{k}\mathrm{m}/\mathrm{s}`$ where the transmitted flux is high have been observed (Reimers et al. (1997); Anderson et al. (1998); Heap et al. (2000); Smette et al. (2000)).

A uniform IGM temperature has been assumed in most of the theoretical and observational Lyman $`\alpha `$ forest work so far. We live in an inhomogeneous universe, inhomogeneous heating of the IGM due to a late HeII reionization is only one of the possible sources of temperature fluctuations. It is timely to investigate possible ways of detecting temperature fluctuations. In this paper we use the small scale power spectrum of the forest flux as a measure of temperature (Theuns et al. (2000); Zaldarriaga, Hui & Tegmark (2000)) an introduce the tools necessary to search for fluctuations in the small scale power.

The paper is organized as follows, in section §2 we summarize the properties of the small scale power spectrum and its dependence on the IGM temperature, in §3 we introduce the tools we use in the search for temperature fluctuations, in §4 we apply our technique to Q1422+231 (Kim et al. (1997)), a quasar at redshift $`z=3.62`$ and in §5 we estimate how much data is needed to substantially improve our constraints. We conclude in §6

## 2. The power spectrum as a measure of the IGM temperature

The power spectrum of the flux in the Lyman $`\alpha `$ forest can be used as an indication of the IGM temperature (Theuns et al. (2000); Zaldarriaga, Hui & Tegmark (2000)). Figure 1 shows the predictions of a set of models with different temperatures together with the measurements of McDonald et al. (1999). The simulations used in this paper are described in Zaldarriaga, Hui & Tegmark (2000). We use PM simulations to solve for the dark matter and then use analytic scaling relations to construct mock flux spectra. Our simulations have $`128^3`$ particles in a box that at $`z=3`$ is $`3200\mathrm{k}\mathrm{m}/\mathrm{s}`$ in size.

Given that the small scale power is sensitive to the temperature of the IGM it is reasonable to try to detect fluctuations in the IGM temperature by searching for fluctuations in the small scale power. To do so we have to measure the power in different positions along the spectrum and/or on lines of sight towards different quasars. We could use Fourier analysis on pieces of the spectra but we prefer to use wavelets, a natural tool to obtain measures of power that are simultaneously localized in real and frequency space. In the next section we review the properties of the wavelet expansion that are relevant for our work.

## 3. Wavelets

Wavelets have been used in a number of statistical studies of the forest. In Meiksin (2000) wavelets were introduced for data compression and several statistical properties of the wavelet coefficients were presented.

For our purposes the more relevant work is that of Theuns & Zaroubi (2000). Just as we will do in this paper, the authors used wavelet coefficients to measure the small scale power in the spectrum. They presented the cumulative probability distribution of the squares of the wavelet coefficients and focus most of their effort on comparing models with different temperatures and equations of state and thus different small scale power. They conclude that the wavelet coefficients in these different models have different enough distributions to be able to tell the models apart observationally. In Zaldarriaga, Hui & Tegmark (2000) we used Fourier analysis and measurements of the small scale power in the literature to constrain the equation of state of the IGM. Both Theuns & Zaroubi (2000) and Zaldarriaga, Hui & Tegmark (2000) are very similar in this respect, they propose the same “thermometer” for the IGM, the small scale power in the forest.

Furthermore, Theuns & Zaroubi (2000) argue that because wavelet coefficients are localized in real space they can be used to measure fluctuations in the temperature. They showed results for the cumulative distribution of a model with two temperatures and compare them to models with one single temperature. Unfortunately it is unclear how much of the difference is due to the models having a different mean temperatures and how much is due to the fluctuations themselves, or how much data is needed to get tight constraints on fluctuations. Our work in this paper builds on the ideas presented in Theuns & Zaroubi (2000) and Zaldarriaga, Hui & Tegmark (2000) with special focus on statistical tests of temperature fluctuations.

In this section we introduce the concepts necessary to understand the use of wavelets in the context of this paper.

### 3.1. Definitions

The best way to introduce the necessary concepts is by looking at the wavelet decomposition of a piece of spectra from Q1422+231 (Kim et al. (1997)). In figure 2 we show $`3200\mathrm{km}/\mathrm{s}`$ of the spectra together with wavelet coefficients. This is approximately a 10th of the Lyman $`\alpha `$ part of the spectra. We will use all the spectra later, here we chose to show a small portion mainly for clarity but also because this is the length of the simulations we will have to compare with. The wavelet transform of the spectra is obtained by convolving the spectra with a set of test functions which we will call $`\psi _n(xx_i)`$,

$$a_n(x_i)=𝑑xF(x)\psi _n(xx_i),$$
(1)
$`x`$ denotes position along the spectra and $`F(x)`$ is the flux. The test functions are characterized by two numbers, the index $`n`$ and $`x_i`$. Increasing values of $`n`$ correspond to higher spatial frequencies, $`x_i`$ determines the position of the wavelet along the spectrum. The values of $`x_i`$ are determined by the wavelet scheme. For the wavelets we are using there are $`2^n`$ values of $`x_i`$ in each order $`n`$. Figure 2 shows a few examples of these wavelets to illustrate how as $`n`$ increases, the wavelet corresponds to higher spatial frequencies. In this work we use the Daubenchies 4 wavelets implemented in a Numerical Recipes routine (Press et al. (1992)).

For our purposes the wavelets are just a set of filters that we apply to the data. They have the advantage of being localized in both Fourier and spatial domain, of being an orthogonal and complete set, and perhaps most importantly that packages that perform the transform are readily available.

To further de-mystify the wavelets we can look at the power spectrum of the forest from a wavelet perspective. Rather than computing squares of Fourier coefficient we can obtain a measure of power by averaging the squares of the wavelets coefficients of a fixed $`n`$ but different positions along the spectrum. We have

$$\widehat{P}\frac{1}{N}\underset{i=1}{\overset{N}{}}a_n^2(x_i),$$
(2)
the constant of proportionality depends on the normalization convention of $`\psi _n`$. The wavelets are well localized in real space so they must have a broad response in Fourier space. These response or window functions for a few orders are shown in figure 3. By construction window functions of successive $`n`$ differ only by a factor of 2 scaling in the $`x`$ axis. In figure 3 we also show the power spectrum measured in the conventional Fourier way or using wavelts. The $`x`$ coordinate of each wavelet measurement is set to the peak of the window function, but the windows are clearly very broad.

We should point out that the $`n`$ label is arbitrary in the sense that the label for a wavelet that probes a particular range of spatial scales depends on the length of the spectrum being analyzed. For example, if we start with a piece of spectrum that is twice as long ($`6400\mathrm{km}/\mathrm{s}`$) all the labels would be shifted by one. That is, the new $`n=7`$ wavelet probes the same scales as the old $`n=6`$ did. As long as one is consistent with the normalization convention the new wavelet coefficients of order $`n+1`$ will be equal to the old coefficients of order $`n`$ at the same position along the spectrum. To avoid confusion we will use the labels corresponding to $`3200\mathrm{k}\mathrm{m}/\mathrm{s}`$ in all our discussions even when we analyze the full length of Q1422+231 spectrum. Figure 3 indicates exactly what scales we probe when we discuss each order $`n`$.

It is useful to discuss how the power in the wavelet coefficients depends on the temperature. In figure 4 we illustrate the rate of change of power with temperature obtained in our simulations. Clearly the smallest spatial scales are most sensitive. The figure also shows that for $`n=7`$ a factor of $`2`$ change in the temperature, which could be expected in models were He has been doubly ionized in patches of the spectrum, could lead to a factor of $`2.5`$ change in the observed power. For the figure me assumed an isothermal model, the temperature of the gas is independent of its density. In Zaldarriaga, Hui & Tegmark (2000) we showed that the small scale power is sensitive to the temperature over a range of overdensities $`0.9<\rho /\overline{\rho }<1.8`$. The ratio of temperatures in figure 4 should be interpreted as the ratio of temperatures on that range of overdensities. This fact can be particularly important in the case of energy injection by He reionization as the process will not only change the temperature of the gas but also its equation of state, so the relevant parameter is by how much the temperature is increased in the range of overdensities that are probed by the small scale power.

In summary wavelets of different scales $`n`$ measure power on different spatial frequencies. Their frequency response is broad, the trade-off for the wavelets being localized in real space. The wavelets probing the small scales are very sensitive to the temperature of the IGM and we will use these coefficients as probes of the temperature. We will take advantage of the localization in real space of the wavelets to search for temperature fluctuations.

### 3.2. Statistical properties of the wavelet coefficients

Now that we have defined the wavelet coefficients we will turn to some analysis of their statistical properties in the case of the Lyman $`\alpha `$ forest flux.

In the rest of the paper we will be interested in the small scale power, that we will be using as a thermometer. In figure 5 we show a histogram of the amplitude of the $`n=7`$ coefficient obtained in our simulations. The variance of this distribution proportional to the flux power spectrum (equation 2). In the the top panel we show the histograms for two different temperatures, $`T=400(800)(\mathrm{km}/\mathrm{s})^2=2.4(4.8)\times 10^4\mathrm{K}`$. Clearly the model with the largest temperature has the smaller small scale power, the variance of the distribution is smaller.

In the bottom panel of figure 5 we have rescaled the $`x`$ axis by dividing it by the variance. The two histograms now fall on top of each other, although there are some small differences. For comparison we also show the histograms of the coefficients of Q1422+231 which have a distribution which is remarkably similar to the coefficients in the simulation.

In figure 5 we show a Gaussian together with the histograms of the coefficients in the simulations and in Q1422+231. There are clear differences between the Gaussian and the rest. The amplitudes of the wavelets of the flux have a probability of being near zero which is much larger than one would expect in a Gaussian of the same variance. This is clearly noticeable in the example of figure 2, several stretches of the spectrum can be identified where the amplitude of the wavelet coefficients are very small. The distribution of the coefficient on the simulation also has tails, it is much probable to get large coefficients than in a Gaussian distribution.

To lay the foundation for our method to probe for temperature fluctuations we also need to study the spatial correlations of the wavelet coefficients. In figure 6 we show the correlation coefficient of the wavelet amplitudes ($`CC`$) as a function of separation along the spectra ($`\delta v`$) for several scale indices $`n`$,

$$CC(\delta v)=\frac{<a_n(x)a_n(x+\delta v)>}{<a_n^2>}.$$
(3)
The correlation coefficient decays fast with separation and then levels off at a value below $`1\%`$. In our future analysis we will be considering only the smaller scales, corresponding to $`n=6`$ or $`7`$. The correlations in these cases become negligible (below $`1\%`$) for separations larger than $`100\mathrm{k}\mathrm{m}/\mathrm{s}`$.

### 3.3. Analytic tools

In order to detect temperature fluctuations we want to look for variations in the statistical properties of the wavelet coefficients as we look at different parts of the spectra or as we compare different lines of sight. It is clear that any statements about fluctuations will be statistical in nature, the variance in the distribution of the wavelet coefficients is our thermometer, not the coefficients themselves. Thus we cannot measure the temperature on a pixel by pixel basis, we need to average over many pixels.

Our measure of power is the square of the wavelet coefficient, thus we will introduce averages of the square of these coefficients which we will call $`w`$,

$$w=\frac{1}{N_p}\underset{i}{}a_n^2(x_i).$$
(4)
The sum is done over contiguous portions of the spectra of length $`L`$. Different portions do not overlap, and there are $`N_p`$ coefficient in each chuck of size $`L`$. $`N_p`$ depends on the order $`n`$ one is considering, $`N_p=2^n(L/3200\mathrm{k}\mathrm{m}/\mathrm{s})`$. By definition the mean of $`w`$ is the measure of power and thus temperature. The mean only depends on the index $`n`$ and not on $`L`$. In figure 7 we show the histogram or probability distribution (PDF) of $`w`$ for several choices of $`L`$ in our simulation. We show the histograms for two temperatures. It is clear from the figure that the model with the larger temperature has a smaller mean $`w`$.

In figure 8 we have rescaled the $`x`$ axis in figure 7 using the mean of $`w`$ ($`\overline{w}`$) for each temperature. The histograms for the two temperatures are quite similar, we are thus encouraged to find some simple analytic model. Our statistic $`w`$ is an average of squares of many variables which as figure 6 shows quickly become uncorrelated. We could expect that a $`\chi ^2`$ distribution would be a good model for the observed distributions. The effective number of degrees of freedom $`\nu `$ of the $`\chi ^2`$ distributions need to be determined. Furthermore because the different panels correspond to values of $`L`$ that increase by factors of 2, the effective number of degrees of freedom should increase by factors of 2. In figure 8 we plot this simple analytic approximation,

$`P(w|\overline{w})`$ $`=`$ $`{\displaystyle \frac{\nu }{\mathrm{\Gamma }(\nu /2)2\overline{w}}}e^uu^{\nu /21}`$
$`u`$ $`=`$ $`w\nu /2\overline{w},`$ (5)
where for the different panels we choose $`\nu =1,2,4,8`$. Changing slightly the values of $`\nu `$ improves or worsens the agreements depending on the panel. We choose these values doing a simple $`\chi `$ by eye of the above plots. Our constraints on the temperature fluctuations are insensitive to the details of this choice, different values of $`\nu `$ give the same answer, as long as the analytic formulas fit the panels reasonably well. Note also that the distribution for the data of Q1422+231 is also quite similar, there are some differences for the largest smoothing window but in that case there are only a few samples Q1422+231, so the discrepancy is not significant.

We will use this analytic form for the pdf of $`w`$ throughout the rest of our work. Although the approximation is not perfect our analytic formulas give a good description what is observed in the simulations. We are going to be searching for variations in this distribution as we move along the spectra, so it is very useful that we have an analytic form and furthermore that this analytic form depends on only one parameter $`(\overline{w})`$. When formulated this way, the question of fluctuations in the temperature can be very simply stated. Does a model with a single value of $`\overline{w}`$ fit the data better than a model were $`\overline{w}`$ varies in some particular way across the spectrum.

We could take a different approaches to look for fluctuations in the statistical properties of the coefficient. We could directly compare the observed distributions with what is obtained in simulations where a uniform temperature was assumed and search for differences. The problem with this approach is that one must establish that this differences are not a consequence of systematic effects introduced in the simulation process by some missing physics or some numerical artifact. We believe that our method is more robust as it does not use simulations directly, only to motivate our choice of distribution. We are essentially looking for changes in the parameters of the distribution we use to model the data as we go along the spectrum. Any change in this parameters is a signal of fluctuations even if our distribution is not a perfect model of the data.

## 4. Constraints on temperature fluctuations

Our objective is to use the small scale power as a thermometer, and this thermometer to search for fluctuations. In the previous sections we have established a model for the probability distribution of the square of the wavelets coefficients. This model is dependent on one parameter, which we chose to be the mean power on that scale. The objective of this section is to describe the formalism necessary to constrain spatial fluctuations in the mean power, such as those that would be produced by fluctuation in the temperature.

Let as call $`𝐰`$ the vector with the square of the amplitude of the wavelet coefficient as a function of position along the spectra. Let us denote $`N`$ the dimension of that vector, which is given by the ratio of the length of the spectrum to length of the averaging filter in equation (4). For simplicity we will call each of the entries of these vectors pixels, although they are a combination of several pixels in the initial spectrum.

In the previous sections we found a model for $`\mathrm{P}(w_i|\overline{w}_i)`$, the probability of measuring $`w_i`$ given that mean $`\overline{w}_i`$ of $`w`$ at position $`i`$ along the spectrum, (equation 3.3). To describe our model for the fluctuation we need to introduce another vector ($`\overline{𝐰}`$) that specifies the value of $`\overline{w}`$ along the spectrum. The dimension of $`\overline{𝐰}`$ is also N. The probability of measuring the vector $`𝐰`$ is given by,

$`\mathrm{P}(𝐰)`$ $`=`$ $`{\displaystyle 𝑑\overline{𝐰}\mathrm{P}(\overline{𝐰})\mathrm{P}(𝐰|\overline{𝐰})}`$ (6)
$`=`$ $`{\displaystyle 𝑑\overline{w}_1\mathrm{}𝑑\overline{w}_n\mathrm{P}(\overline{𝐰})\underset{i=1}{\overset{N}{}}\mathrm{P}(w_i|\overline{w}_i)}.`$
Different models for the temperature distribution differ in what $`\mathrm{P}(\overline{𝐰})`$ is. For example for a uniform temperature model where for every pixel the mean power is $`\overline{w}`$ we have,

$$\mathrm{P}(\overline{𝐰})=\delta ^D(\overline{w}_1\overline{w})\mathrm{}\delta ^D(\overline{w}_N\overline{w})$$
(7)
By combining equations (3.3), (6) and (7) we obtain the probability of measuring $`𝐰`$ under the assumption that the temperature is uniform ($`\mathrm{P}^u`$),

$`\mathrm{ln}\mathrm{P}^u`$ $`=`$ $`{\displaystyle \underset{i}{}}\mathrm{ln}\mathrm{P}(w_i|\overline{w})`$ (8)
$`=`$ $`{\displaystyle \underset{i}{}}(\nu /21)\mathrm{ln}(u_i)u_i`$
$``$ $`N\mathrm{ln}[2\overline{w}\mathrm{\Gamma }(\nu /2)/\nu ],`$
where $`u_i=w_i\nu /2\overline{w}`$.

We want to constrain fluctuations in the temperature by comparing the probability of having measured a vector $`𝐰`$ in a particular model for the temperature fluctuations ($`\mathrm{P}^f`$) with the same probability in a uniform temperature model ($`\mathrm{P}^u`$). The ratio of $`\mathrm{P}^f`$ and $`\mathrm{P}^u`$ is the likelihood ratio for the two models. In order to compute $`\mathrm{P}^f`$ we need to specify the value of $`\mathrm{P}(\overline{𝐰})`$, our model for the temperature fluctuations. Clearly there are many different choices of $`\mathrm{P}(\overline{𝐰})`$ which describe different physical scenarios and it would be impossible for us to try every possibility. Instead in this section we have chosen some particular examples and found constraints in those scenarios.

In all our models we will assume that there are only two different temperature values possible and that the mean power for each of these two temperatures $`\overline{w}_A`$ and $`\overline{w}_B`$, differ by a factor $`\beta =\overline{w}_B/\overline{w}_A`$. We will also define $`f`$ as the fraction of the spectrum with mean $`\overline{w}_A`$. The mean power averaged over all pixels becomes,

$$\overline{w}=f\overline{w}_A+(1f)\overline{w}_B.$$
(9)
When comparing different models we will keep this mean, $`\overline{w}`$ fixed to the observed value in the spectrum. We are not interested in comparing models with different mean power but models with or without spatial temperature fluctuations in the power.

There are two probability ratios that define the problem,

$`\mathrm{ln}{\displaystyle \frac{\mathrm{P}(w_i|\overline{w}_A)}{\mathrm{P}(w_i|\overline{w})}}`$ $`=`$ $`u_i(1f)(1\beta )`$
$`+`$ $`{\displaystyle \frac{\nu }{2}}\mathrm{ln}(f+\beta (1f))`$
$`\mathrm{ln}{\displaystyle \frac{\mathrm{P}(w_i|\overline{w}_B)}{\mathrm{P}(w_i|\overline{w})}}`$ $`=`$ $`u_if{\displaystyle \frac{(\beta 1)}{\beta }}`$ (10)
$`+`$ $`{\displaystyle \frac{\nu }{2}}\mathrm{ln}({\displaystyle \frac{f+\beta (1f)}{\beta }}),`$
where $`u_i=w_i\nu /2\overline{w}`$. The expressions for $`\mathrm{P}^f/\mathrm{P}^u`$ in the different models for the fluctuations will always be expressed in terms of these two quantities.

### 4.1. Independent pixels

We first consider a model in which every pixel is independent of the rest and can be in one of two temperatures, corresponding to mean powers $`\overline{w}_A`$ and $`\overline{w}_B`$. The prior on the small scale power becomes,

$$\mathrm{P}(\overline{𝐰})=\underset{i}{}[f\delta ^D(\overline{w}_i\overline{w}_A)+(1f)\delta ^D(\overline{w}_i\overline{w}_B)].$$
(11)
The likelihood ratio $`\mathrm{P}^f/\mathrm{P}^u`$ is,

$`\mathrm{ln}[\mathrm{P}^f/\mathrm{P}^u]`$ $`=`$ $`{\displaystyle \underset{i}{}}\mathrm{ln}[f{\displaystyle \frac{P(w_i|\overline{w}_A)}{P(w_i|\overline{w})}}`$ (12)
$`+`$ $`(1f){\displaystyle \frac{P(w_i|\overline{w}_B)}{P(w_i|\overline{w})}}].`$
This expression remains constant if we transform $`\beta 1/\beta `$ and $`f(1f)`$.

We use equation (12) and the full spectra of Q1422+231 to put simultaneous constraints on $`f`$ and $`\beta `$ which are shown in figure 9. We do not need to consider values of $`\beta <1`$ because they can be obtained from $`\beta >1`$ using the symmetry property of the likelihood ratio. We considered models with a likelihood ratio smaller than 0.05 as ruled out. As can be expected the tightest constraints on $`\beta `$ are obtained when $`f0.5`$. With the data we have we can say that $`\beta <3.5`$ ($`f0.5`$).

As could have been guessed the tightest constraints are obtained when both temperatures occupy a significant fraction of the spectrum, $`f0.5`$. However it is interesting to note, that the contour is not symmetric around $`f=0.5`$. Figure 9 shows that it is harder to put constraints when $`f`$ is near zero than when it is near one. The reason for this is that when $`f0`$ and $`\beta >1`$, most of the spectrum has the larger small scale power $`\overline{w}_B`$ (low temperature) and only a small fraction has the lower power $`\overline{w}_A`$ (high temperature). Even when $`\overline{w}`$ is high, figure 8 shows that it is more likely that in any given pixel the value of $`w`$ is low, thus it is hard to detect the small excess of pixels with low $`w`$s that is produced by the small fraction of the spectrum that has $`\overline{w}_A`$.

### 4.2. Sharp transition

We now consider a case in which there is a sharp transition between $`\overline{w}_A`$ and $`\overline{w}_B`$ which occurs at pixel number $`i_{tr}`$. The prior becomes,

$$\mathrm{P}(\overline{𝐰})=\underset{i<i_{tr}}{}\delta ^D(\overline{w}_i\overline{w}_A)\underset{i>i_{tr}}{}\delta ^D(\overline{w}_i\overline{w}_B).$$
(13)
The likelihood ratio $`\mathrm{P}^f/\mathrm{P}^u`$ for this case becomes,

$`\mathrm{ln}[\mathrm{P}^f/\mathrm{P}^u]`$ $`=`$ $`{\displaystyle \underset{i<i_{tr}}{}}\mathrm{ln}[{\displaystyle \frac{P(w_i|\overline{w}_A)}{P(w_i|\overline{w})}}]`$ (14)
$`+`$ $`{\displaystyle \underset{i>i_{tr}}{}}\mathrm{ln}[{\displaystyle \frac{P(w_i|\overline{w}_B)}{P(w_i|\overline{w})}}].`$
For this model we define $`f=i_{tr}/N`$. We compute $`\mathrm{P}^f/\mathrm{P}^u`$ for all possible values of $`i_{tr}`$ and plot the constraints on $`\beta `$ in figure 9. We plot the results in terms of $`f`$ rather than $`i_{tr}`$. For this model we plot constraints from both $`\beta `$ bigger and smaller than one, as both cases are not equivalent. One corresponds to a jump from low to high temperature and the other to the opposite case.

These model is the one where the tightest constraints can be obtained because it is the less random. Given $`i_{tr}`$ the temperature of all pixels is specified. We can can constrain $`0.5\beta 2`$ for almost all values of $`f`$.

### 4.3. Fixed size bubbles

Finally we will consider a case in which there are bubbles of a fixed size with power $`\overline{w}_B`$ and the rest of the spectrum has $`\overline{w}_A`$. To compute the $`\mathrm{P}(𝐰)`$ under this scenario we reinterpret equation (6) as an average over $`\mathrm{P}(\overline{𝐰})`$ of $`\mathrm{P}(𝐰|\overline{𝐰})`$. We can use a Montecarlo technique to perform this average. We create realizations simply by laying down $`N_b`$ bubbles of a fixed size in the spectrum. We choose randomly the centers of the bubbles. For each realization we compute the likelihood using equation (4) and the average that over all realizations to compute the $`\mathrm{P}(𝐰)`$. The model has two free parameters, the size of each individual bubble and the number of bubbles (or equivalently $`f`$, the fraction of the spectrum having $`\overline{w}_A`$ and $`\overline{w}_B`$).

In figure 9 we show the constraints on $`f`$ and $`\beta `$ for two particular bubble sizes, $`1200\mathrm{km}/s`$ and $`2400\mathrm{km}/s`$. In this model $`\beta `$s larger and smaller than one are not equivalent. Our results for independent pixels can be interpreted as a constraint when the typical size of the bubble is $`400\mathrm{km}/s`$. As can be seen in figure 9 our constraints end up being very similar for all bubble sizes. For $`f0.5`$ the ratio of powers have to be in the range $`0.3\beta 3.5`$.

## 5. Future

In the previous sections we have developed a method for constraining temperature fluctuations. We illustrated our method using data from Q1422+231. The constraints on the ratio of the powers obtained are somewhat larger than what one might expect to be present in the universe. In this section we compute the number of spectra needed to substantially improve the constraints.

The statistical question we will ask is what constraints on $`\beta `$ can be obtained from a certain amount of data in the independent pixel model under the assumption that the universe has uniform temperature. Equation (12) implies that the ratio of the probabilities of a uniform model to the fluctuating temperature model ($`\mathrm{P}^f/\mathrm{P}^u`$) is,

$`\mathrm{ln}(\mathrm{P}^f/\mathrm{P}^u)`$ $`=`$ $`{\displaystyle \underset{i}{}}\mathrm{ln}[f{\displaystyle \frac{P(w_i|\overline{w}_A)}{P(w_i|\overline{w})}}`$ (15)
$`+`$ $`(1f){\displaystyle \frac{P(w_i|\overline{w}_B)}{P(w_i|\overline{w})}}]`$
$``$ $`N{\displaystyle 𝑑w\mathrm{P}(w|\overline{w})}`$
$`\mathrm{ln}`$ $`[f{\displaystyle \frac{P(w|\overline{w}_A)}{P(w|\overline{w})}}+(1f){\displaystyle \frac{P(w|\overline{w}_B)}{P(w|\overline{w})}}]`$
In the last line we approximated the sum by a mean over the distribution from which the $`w_i`$ are assumed to be drawn, $`\mathrm{P}(w|\overline{w})`$, and we used their independence. The last line of equation (15) is a measure of distance between two probability distribution functions called Kullback-Leibler entropy (Kullback (1959)).

We can set the ratio in equation (15) to $`0.05`$ and solve for $`N`$ as a function of $`\beta `$ and $`f`$. We show the results in figure 10. We plot only points for $`\beta >1`$, other values of $`\beta `$ can be obtained using the transformation $`\beta 1/\beta `$, $`f(1f)`$. We have changed variables form $`N`$ to $`N_q`$, the number of quasars with the same number of pixels as Q1422+231. The figure shows that ten quasars are needed to obtain constraints on $`\beta 2`$ for $`f0.5`$. The figure also illustrates the fact that it is very difficult to obtain constraints that are much tighter than a factor 2 in $`\beta `$.

It is interesting to understand why it gets so much harder to obtain constraints when $`\beta `$ approches one. Figure 10 shows that the number of quasars scales as $`N_q1/(\beta 1)^4`$. The key point is that we are comparing models that by construction have the same mean $`\overline{w}`$ (the same mean temperature). Thus what we are actually comparing is the variance of the $`w`$ distribution in both models. We can compute how many samples are needed to distinguish two models that have different variance. We define the variance as $`v=<w^2>\overline{w}^2`$, its value depends on the model (fluctuating or uniform temperature), but $`\overline{w}`$ is the same across models. We are assuming all the pixels are independant so we can estimate the variance as,

$$\widehat{v}=\frac{1}{N_p}\underset{i}{}w_i^2\overline{w}^2.$$
(16)
The question becomes how many samples $`(N_p)`$ are needed to distinguish a model with uniform temperature with one characterized by $`\beta `$ and $`f`$. We define the signal to noise ratio as,

$$(S/N)^2=[v^fv^u]^2/\mathrm{var}(v),$$
(17)
the ratio of the difference in variance between the model with fluctuations and the one without to the variance of the variance (which we calculate in the uniform model).

In our example the signal to noise ratio becomes,

$$(S/N)^2=cN_p\frac{(1\beta )^4f^2(1f)^2}{[f+\beta (1f)]^4},$$
(18)
where the constant $`c`$ depends on the shape of the distribution of $`w`$ and is defined as $`c=<w^2>^2/(<w^4><w^2>^2)`$. The scaling of the signal to noise with the number of pixels is $`(1\beta )^4`$ as figure 5 shows. It is the fact that we are looking at the difference in the variance of $`w`$ for models constructed to have the same $`\overline{w}`$, what makes distinguishing the models so difficult as $`\beta 1`$. This fact is what forces $`N_q`$ to scale as $`1/(\beta 1)^4`$ instead of $`1/(\beta 1)^2`$.

Finally with more data we could not only get better constraints on model parameters such as $`\beta `$ but also constrain more ambitious models. For example, rather than have fixed size bubbles we could draw the bubble sizes from some distribution or allow for correlations between their positions. Moreover with more data one has to study in detail certain sources of systematic error in our method like contamination by metal lines or the effect of evolution along the spectra.

## 6. Conclusions

The uniformity of the temperature in the IGM is an untested assumption in most theoretical modelling of the Lyman $`\alpha `$ forest. We have presented a statistical method to test this assumption. Our method uses the small scale power of the forest flux as a thermometer and searches for fluctuations in this power. As a result the method is statistical in nature. An important advantage of the method we suggest is that it does not require a direct comparison with numerical simulations. It searches for changes in the statistical properties in different parts of the spectra.

We illustrated our method with Q1422+231. We assumed that in different regions of the spectrum the IGM could be in one of two temperatures. We found that the ratio of small scale power in these two type regions could not differ by more than a factor of $`3.5`$ if the two type of regions occupy a comparable fraction of the spectrum. This ratio of power translates in a factor of approximately $`2.5`$ in temperature. If the transition between the two regimes is sharp, the constraints are more stringent, the ratio of powers cannot differ by more than a factor of $`2`$.

We have calculated how much data is needed to improve this constraints substantially. We showed that ten quasars are needed to constrain fluctuations of factors of two in the power, or equivalently factors of $`1.7`$ in temperature. Furthermore with more data, more ambitious models for the fluctuations can be considered.

Acknowledgements The author acknowledges very useful comments and suggestions from Joop Schaye and John Bahcall. Support for this work is provided by the Hubble Fellowship HF-01116-01-98A from STScI, operated by AURA, Inc. under NASA contract NAS5-26555.