text:
# 1 Introduction

## 1 Introduction

The anomalous weak dipole moments of the $`\tau `$ lepton are the tensorial couplings of the $`\mathrm{Z}\tau ^+\tau ^{}`$ vertex. They are zero to first order in the Standard Model (SM). Two types of anomalous weak dipole moments can be distinguished: the magnetic term $`\mu _\tau `$ and the CP-violating electric term $`d_\tau `$. Here, both the real and the imaginary components of each anomalous weak dipole moment are explored, i.e. $`\mathrm{Re}\mu _\tau `$, $`\mathrm{Im}\mu _\tau `$, $`\mathrm{Re}d_\tau `$ and $`\mathrm{Im}d_\tau `$. Radiative corrections in the SM provide nonzero predictions for $`\mu _\tau `$ and $`d_\tau `$  which are below the present experimental sensitivity. This opens the possibility to look for deviations from the SM.

There have been many searches for the CP-violating anomalous weak electric dipole moment of the $`\tau `$ since the beginning of LEP \[3-5\]. In addition, limits on the anomalous weak magnetic dipole moment were obtained more recently .

In this analysis, the previous ALEPH result on $`\mathrm{Re}d_\tau `$  is updated, and $`\mathrm{Re}\mu _\tau `$, $`\mathrm{Im}\mu _\tau `$ and $`\mathrm{Im}d_\tau `$ are determined for the first time in ALEPH. The data sample was collected with the ALEPH detector from 1990 to 1995 at energies around the $`\mathrm{Z}`$ resonance and corresponds to an integrated luminosity of 155 pb-1. Tau leptons are generated in the reaction $`e^+e^{}\tau ^+\tau ^{}`$ at LEP. The method to extract the anomalous weak dipole moments is based on a maximum likelihood fit to the data taking into account all $`\tau `$ spin terms explicitly, including correlations. This is the first time that the complete differential cross section for the production and decay of the $`\tau `$ leptons is considered to estimate the $`\tau `$ anomalous weak dipole moments. The most important semileptonic decays are used: $`\pi `$, $`\rho `$, $`a_1\pi 2\pi ^0`$ and $`a_13\pi `$. The $`\tau `$ spin information is recovered using optimal polarimeters which are different for each decay. The selection and particle identification make use of tools already developed in previous analyses \[6-9\].

The text is organized as follows. The theoretical framework is introduced in Section 2. The most important ALEPH subdetectors for this analysis are covered in Section 3. The data analysis procedure is explained in Section 4, emphasizing the new features of the analysis. The more relevant systematic uncertainties are then discussed in Section 5. The results and conclusions are presented in Section 6.

## 2 Theoretical framework

### 2.1 Production cross section

The currents assumed for photon and $`\mathrm{Z}`$ exchange in $`e^+e^{}\tau ^+\tau ^{}`$ production are

$`\mathrm{\Gamma }_f^{\mu (\gamma )}`$ $`=`$ $`iQ_f\mathrm{e}\gamma ^\mu ,\text{with }f\text{ = }e\text{}\tau ,`$
$`\mathrm{\Gamma }_e^{\mu (\mathrm{Z})}`$ $`=`$ $`i\mathrm{e}\left[v_e\gamma ^\mu a_e\gamma ^\mu \gamma _5\right],`$
$`\mathrm{\Gamma }_\tau ^{\mu (\mathrm{Z})}`$ $`=`$ $`i\mathrm{e}\left[v_\tau \gamma ^\mu a_\tau \gamma ^\mu \gamma _5+i{\displaystyle \frac{\mu _\tau }{2m_\tau }}\sigma ^{\mu \nu }q_\nu +{\displaystyle \frac{d_\tau }{2m_\tau }}\gamma _5\sigma ^{\mu \nu }q_\nu \right],`$ (1)
where $`Q_f\mathrm{e}`$ is the fermion charge; $`a_e`$, $`a_\tau `$, $`v_e`$ and $`v_\tau `$ are the axial vector and vector couplings of the SM; $`\mu _\tau `$ and $`d_\tau `$ are the anomalous weak magnetic and anomalous weak electric dipole moments of the $`\tau `$. In the previous expression both anomalous weak dipole moments are dimensionless quantities. However, the anomalous weak electric dipole moment is often quoted in the literature in units of $`\mathrm{e}\mathrm{cm}`$, by defining the contribution of this dipole moment to the current as $`id_\tau \gamma _5\sigma ^{\mu \nu }q_\nu `$. These different notations are related by the conversion factor e/2$`m_\tau `$ = $`5.552\times 10^{15}\mathrm{e}\mathrm{cm}`$.

Using the currents in Eq. 1, the differential cross section can be expressed as

$`{\displaystyle \frac{d\sigma }{d\mathrm{cos}\theta _\tau }}(\stackrel{}{s}_1,\stackrel{}{s}_2)`$ $`=`$ $`R_{00}+{\displaystyle \underset{\mu =1,3}{}}R_{\mu 0}s_1^\mu +{\displaystyle \underset{\nu =1,3}{}}R_{0\nu }s_2^\nu +{\displaystyle \underset{\mu ,\nu =1,3}{}}R_{\mu \nu }s_1^\mu s_2^\nu .`$ (2)
The $`R_{\mu \nu }`$ terms are functions of the fermion couplings and of the $`\tau `$ production angle $`\theta _\tau `$; $`\stackrel{}{s}_1`$ and $`\stackrel{}{s}_2`$ are unit vectors chosen as the quantisation axes for the spin measurement of the $`\tau ^+`$ and the $`\tau ^{}`$, respectively, in their corresponding rest frames.

The following reference frame has been chosen: the $`z`$ axis is in the outgoing $`\tau ^+`$ direction and the incoming $`e^+`$ is in the $`yz`$ plane. The $`x`$ component is therefore normal to the production plane. The $`y`$ component is called transverse.

Several $`R_{\mu \nu }`$ terms have been already measured by ALEPH. Defining $`\left(R_{\mu \nu }\right)_\pm (R_{\mu \nu }\pm R_{\nu \mu })`$, these terms are the following:

* $`R_{00}`$ = $`d\sigma `$/ $`d\mathrm{cos}\theta _\tau `$ ,
* $`(R_{03})_+/R_{00}=P_\tau (\mathrm{cos}\theta _\tau )`$ is the longitudinal polarisation of the $`\tau `$ ,
* $`R_{22}/R_{00}=R_{11}/R_{00}`$ are the transverse-transverse and normal-normal spin correlations ,
* $`(R_{21})_+/R_{00}`$ are the transverse-normal spin correlations .

The $`R_{\mu \nu }`$ terms most sensitive to $`\mathrm{Re}\mu _\tau `$, $`\mathrm{Im}\mu _\tau `$, $`\mathrm{Re}d_\tau `$ and $`\mathrm{Im}d_\tau `$ are presented below. These terms are obtained from Ref.  after some algebra. The SM contributions are separated from the anomalous (anm) contributions.

$`\mathrm{Re}\mu _\tau `$:

$`\left(R_{02}\right)_+|_{\mathrm{SM}}`$ $``$ $`{\displaystyle \frac{2}{\gamma _\tau }}\mathrm{sin}\theta _\tau |v_\tau |^2Re(v_ea_e^{})+{\displaystyle \frac{1}{\gamma _\tau }}\mathrm{sin}\theta _\tau \mathrm{cos}\theta _\tau (|a_e|^2+|v_e|^2)Re(v_\tau a_\tau ^{})`$
$`\left(R_{02}\right)_+|_{\mathrm{anm}}`$ $``$ $`\gamma _\tau \mathrm{sin}\theta _\tau \mathrm{cos}\theta _\tau (|a_e|^2+|v_e|^2)Re(a_\tau \mu _\tau ^{})`$ (3)
$`+{\displaystyle \frac{2(\gamma _\tau ^2+1)}{\gamma _\tau }}\mathrm{sin}\theta _\tau Re(v_ea_e^{})Re(v_\tau \mu _\tau ^{})+2\gamma _\tau \mathrm{sin}\theta _\tau Re(v_ea_e^{})|\mu _\tau |^2`$
$`\left(R_{32}\right)_+|_{\mathrm{SM}}`$ $``$ $`{\displaystyle \frac{2}{\gamma _\tau }}\mathrm{sin}\theta _\tau Re(v_ea_e^{})Re(v_\tau a_\tau ^{})+{\displaystyle \frac{1}{\gamma _\tau }}\mathrm{sin}\theta _\tau \mathrm{cos}\theta _\tau (|a_e|^2+|v_e|^2)|v_\tau |^2`$
$`\left(R_{32}\right)_+|_{\mathrm{anm}}`$ $``$ $`{\displaystyle \frac{\gamma _\tau ^2+1}{\gamma _\tau }}\mathrm{sin}\theta _\tau \mathrm{cos}\theta _\tau (|a_e|^2+|v_e|^2)Re(v_\tau \mu _\tau ^{})`$ (4)
$`+2\gamma _\tau \mathrm{sin}\theta _\tau Re(v_ea_e^{})Re(a_\tau \mu _\tau ^{})+\gamma _\tau \mathrm{sin}\theta _\tau \mathrm{cos}\theta _\tau (|a_e|^2+|v_e|^2)|\mu _\tau |^2`$
$`\mathrm{Im}\mu _\tau `$:

$`\left(R_{31}\right)_+|_{\mathrm{SM}}`$ $``$ $`{\displaystyle \frac{1}{\gamma _\tau }}\mathrm{sin}\theta _\tau \mathrm{cos}\theta _\tau (|a_e|^2+|v_e|^2)Im(v_\tau ^{}a_\tau )`$
$`\left(R_{31}\right)_+|_{\mathrm{anm}}`$ $``$ $`\gamma _\tau \mathrm{sin}\theta _\tau \mathrm{cos}\theta _\tau (|a_e|^2+|v_e|^2)Im(a_\tau \mu _\tau ^{})`$ (5)
$`+{\displaystyle \frac{2(\gamma _\tau ^21)}{\gamma _\tau }}\mathrm{sin}\theta _\tau Re(v_ea_e^{})Im(v_\tau \mu _\tau ^{})`$
$`\left(R_{01}\right)_+|_{\mathrm{SM}}`$ $``$ $`{\displaystyle \frac{2}{\gamma _\tau }}\mathrm{sin}\theta _\tau Re(v_ea_e^{})Im(v_\tau ^{}a_\tau )`$
$`\left(R_{01}\right)_+|_{\mathrm{anm}}`$ $``$ $`{\displaystyle \frac{\gamma _\tau ^21}{\gamma _\tau }}\mathrm{sin}\theta _\tau \mathrm{cos}\theta _\tau (|a_e|^2+|v_e|^2)Im(v_\tau \mu _\tau ^{})`$ (6)
$`+2\gamma _\tau \mathrm{sin}\theta _\tau Re(v_ea_e^{})Im(a_\tau \mu _\tau ^{})`$
$`\mathrm{Re}d_\tau `$:

$`\left(R_{01}\right)_{}|_{\mathrm{SM}}`$ $``$ $`0`$
$`\left(R_{01}\right)_{}|_{\mathrm{anm}}`$ $``$ $`\gamma _\tau \mathrm{sin}\theta _\tau \mathrm{cos}\theta _\tau (|a_e|^2+|v_e|^2)Re(a_\tau d_\tau ^{})`$ (7)
$`2\gamma _\tau \mathrm{sin}\theta _\tau Re(v_ea_e^{})\left[Re(v_\tau d_\tau ^{})+Re(\mu _\tau d_\tau ^{})\right]`$
$`\left(R_{31}\right)_{}|_{\mathrm{SM}}`$ $``$ $`0`$
$`\left(R_{31}\right)_{}|_{\mathrm{anm}}`$ $``$ $`\gamma _\tau \mathrm{sin}\theta _\tau \mathrm{cos}\theta _\tau (|a_e|^2+|v_e|^2)\left[Re(v_\tau d_\tau ^{})+Re(\mu _\tau d_\tau ^{})\right]`$ (8)
$`2\gamma _\tau \mathrm{sin}\theta _\tau Re(v_ea_e^{})Re(a_\tau d_\tau ^{})`$
$`\mathrm{Im}d_\tau `$:

$`\left(R_{32}\right)_{}|_{\mathrm{SM}}`$ $``$ $`0`$
$`\left(R_{32}\right)_{}|_{\mathrm{anm}}`$ $``$ $`\gamma _\tau \mathrm{sin}\theta _\tau \mathrm{cos}\theta _\tau (|a_e|^2+|v_e|^2)Im(a_\tau d_\tau ^{})`$ (9)
$`+2\gamma _\tau \mathrm{sin}\theta _\tau Re(v_ea_e^{})\left[Im(v_\tau d_\tau ^{})+Im(\mu _\tau d_\tau ^{})\right]`$
$`\left(R_{02}\right)_{}|_{\mathrm{SM}}`$ $``$ $`0`$
$`\left(R_{02}\right)_{}|_{\mathrm{anm}}`$ $``$ $`\gamma _\tau \mathrm{sin}\theta _\tau \mathrm{cos}\theta _\tau (|a_e|^2+|v_e|^2)\left[Im(v_\tau d_\tau ^{})+Im(\mu _\tau d_\tau ^{})\right]`$ (10)
$`+2\gamma _\tau \mathrm{sin}\theta _\tau Re(v_ea_e^{})Im(a_\tau d_\tau ^{})`$
Taking into account that $`a_lv_l`$, the terms can be ordered in sensitivity, and the most sensitive term for each anomalous weak dipole moment is presented first. The quantity $`\gamma _\tau `$ is computed as $`\sqrt{s}/2m_\tau `$. The photon exchange terms are omitted from these expressions for simplicity, although they are taken into account in the final results.

The anomalous weak dipole moments are extracted including all $`R_{\mu \nu }`$ terms in a maximum likelihood fit. In this analysis $`\left(R_{31}\right)_+`$, the most sensitive term to $`\mathrm{Im}\mu _\tau `$, is used for the first time as proposed in Ref. . The terms $`\left(R_{02}\right)_+`$, $`\left(R_{01}\right)_{}`$ and $`\left(R_{32}\right)_{}`$ were previously used in other measurements of the anomalous weak dipole moments.

### 2.2 Tau decay

For each $`\tau `$ decay mode, the differential partial width of a polarised $`\tau `$ is written as

$$d\mathrm{\Gamma }(\stackrel{}{s})=W(1+\stackrel{}{h}\stackrel{}{s})dX,$$
(11)
using the expressions for $`W`$ and $`\stackrel{}{h}`$ from the TAUOLA Monte Carlo program ; $`W`$ is the differential partial width of an unpolarised $`\tau `$, and the $`\stackrel{}{h}`$ vector is the polarimeter of the particular decay mode considered. Both $`W`$ and $`\stackrel{}{h}`$ depend on the four-momenta of the final state particles in the $`\tau `$ rest frame, and they are different for each decay topology. The simplest expressions are those of the $`\tau `$ decay into $`\pi `$. In this case, $`\stackrel{}{h}_\pi `$ is proportional to the $`\pi `$ momentum in the $`\tau `$ rest frame and $`W_\pi `$ is a constant. In the above equation, $`X`$ is a set of independent variables describing the full decay configuration. The number of elements of the set depends on the number of particles in the final state. The set $`X_\pi `$ denotes the set of variables expressing the $`\pi `$ direction in the $`\tau `$ rest frame.

In this analysis, the expressions for $`W`$ and $`\stackrel{}{h}`$ allow the spin information for all the $`\tau `$ decays to be recovered optimally.

### 2.3 The full differential cross section

Once the production cross section and the partial decays are introduced, the full differential cross section of $`e^+e^{}\tau ^+\tau ^{}x_1^+x_2^{}\overline{\nu }_\tau \nu _\tau `$ is built following Refs. , namely

$$\frac{d\sigma }{d\mathrm{cos}\theta _\tau dX_1dX_2}=4\frac{W_1}{\mathrm{\Gamma }_\tau }\frac{W_2}{\mathrm{\Gamma }_\tau }\left[R_{00}+\underset{\mu =1,3}{}R_{\mu 0}h_1^\mu +\underset{\nu =1,3}{}R_{0\nu }h_2^\nu +\underset{\mu ,\nu =1,3}{}R_{\mu \nu }h_1^\mu h_2^\nu \right].$$
(12)
In this equation $`\mathrm{\Gamma }_\tau `$ is the total $`\tau `$ width, $`X_1`$ and $`X_2`$ are the sets of independent variables, and $`\stackrel{}{h}_1`$ and $`\stackrel{}{h}_2`$ are the polarimeters for the decay of the $`\tau ^+`$ and the $`\tau ^{}`$, respectively.

With the definitions $`\overline{R}_{\mu \nu }=R_{\mu \nu }/R_{00}`$, $`H^\mu =Wh^\mu /\mathrm{\Gamma }_\tau `$ ($`\mu ,\nu =0,\mathrm{},3`$), and $`h^0=1`$, the likelihood of an event with the final state topology $`ij`$ is written as

$`L_{ij}(\mu _\tau ,d_\tau |\theta _\tau ,W_1,\mathrm{cos}\theta _{h_1},\varphi _{h_1},W_2,\mathrm{cos}\theta _{h_2},\varphi _{h_2})=`$
$`\text{ }{\displaystyle \underset{\mu ,\nu =0,\mathrm{},3}{}}\overline{R}_{\mu \nu }(\mu _\tau ,d_\tau ,\theta _\tau )H_i^\mu (W_1,\mathrm{cos}\theta _{h_1},\varphi _{h_1})H_j^\nu (W_2,\mathrm{cos}\theta _{h_2},\varphi _{h_2}).`$
The indices ($`i`$, $`j`$) refer to the decay mode of each $`\tau `$, with $`i,j=\pi ,\rho ,\pi 2\pi ^0,3\pi `$, and the quantities ($`W_1,\mathrm{cos}\theta _{h_1},\varphi _{h_1}`$) and ($`W_2,\mathrm{cos}\theta _{h_2},\varphi _{h_2}`$) are the observables related to the decay of the $`\tau ^+`$ and the $`\tau ^{}`$, respectively. The angles ($`\theta _{h_1}`$, $`\varphi _{h_1}`$, $`\theta _{h_2}`$, $`\varphi _{h_2}`$) are the polar and azimuthal angles of the polarimeters of each $`\tau `$, in the reference frame introduced in Section 2.1. The above likelihood is also a function of the centre-of-mass energy. The distributions of the hemisphere observables $`W`$ and $`\mathrm{cos}\theta _h`$ are presented in Figs. 1 and 2.

This likelihood fulfills the normalisation condition

$$\underset{ij}{}L_{ij}(\mu _\tau ,d_\tau |\theta _\tau ,W_1,\mathrm{cos}\theta _{h_1},\varphi _{h_1},W_2,\mathrm{cos}\theta _{h_2},\varphi _{h_2})𝑑X_1𝑑X_2=1.$$
(14)
This is the integral over all possible decay parameters and all possible decay topologies $`ij`$ for a given $`e^+e^{}\tau ^+\tau ^{}`$ event. The normalisation is such that the likelihood depends only upon the net spin polarisation of the produced $`\tau `$ pairs, and not upon $`R_{00}`$ = $`d\sigma `$/$`d\mathrm{cos}\theta _\tau `$.

## 3 Apparatus

The ALEPH detector is described in detail in  and its performance in .

Charged particles are measured with a high resolution silicon vertex detector (VDET), a cylindrical drift chamber (ITC), and a large time projection chamber (TPC). The momentum resolution in the axial magnetic field of 1.5 T provided by a superconducting solenoid is $`\mathrm{\Delta }p/p^2=0.6\times 10^3`$(GeV/$`c`$)-1 for high momentum tracks. The impact parameter resolutions for high momentum tracks with hits in all three subdetectors are $`\sigma _{r\varphi }=23\mu `$m and $`\sigma _z=28\mu `$m.

The tracking devices are surrounded by the electromagnetic calorimeter (ECAL), which is a highly segmented lead/proportional-wire-chamber calorimeter. The calorimeter is read out via cathode pads arranged in projective towers covering $`0.9^{}\times 0.9^{}`$ in solid angle and summing the deposited energy in three sections of depth. A second readout is provided by the signals from the anode wires. The energy resolution is $`\sigma /E=0.009+0.18/\sqrt{E(\mathrm{GeV})}`$.

The ECAL is inside the solenoid, which is followed by the hadron calorimeter (HCAL). Hadronic showers are sampled by 23 planes of streamer tubes giving a digital hit pattern and an analog signal on pads, which are also arranged in projective towers. This calorimeter is used in this analysis to discriminate between pions and muons. Outside the HCAL there are two layers of muon chambers providing additional information for $`\mu `$ identification.

## 4 Data analysis

### 4.1 Selection and decay classification

Events from $`\mathrm{Z}\tau ^+\tau ^{}`$ are retained using a global selection in which each event is divided into two hemispheres along the thrust axis. The selection is that used in the ALEPH measurement of $`P_\tau (\mathrm{cos}\theta _\tau )`$ with the $`\tau `$ direction method . Additional information can be found in  and references therein.

The charged particle identification is based on a likelihood method which assigns a set of probabilities to each particle. A detailed description of the method can be found in Refs. . The probability set for each particle is obtained from ($`i`$) the specific ionisation dE/dx in the TPC, ($`ii`$) the longitudinal and transverse shower profiles in ECAL near the extrapolated track and ($`iii`$) the energy and average shower width in HCAL, together with the number out of the last ten of HCAL planes that fired and the number of hits in the muon chambers.

The photon and $`\pi ^0`$ reconstruction is performed with a likelihood method which first distinguishes between genuine and fake photons produced by hadronic interactions in ECAL or by electromagnetic shower fluctuations . All photon pairs in each hemisphere are then assigned a probability of being generated by a $`\pi ^0`$. High energy $`\pi ^0`$ with overlapping showers are reconstructed through an analysis of the spatial energy deposition in the ECAL towers. All the remaining single photons are considered and those with a high probability of being a genuine photon are selected as $`\pi ^0`$ candidates. Finally, photon conversions are identified following the procedure described in . They are added to the list of good photons and are included in the $`\pi ^0`$ reconstruction.

The $`\tau `$ decay classification depends on the number of charged tracks and their identification, and on the number of reconstructed $`\pi ^0`$. It follows the classification for the measurement of $`P_\tau (\mathrm{cos}\theta _\tau )`$ with the $`\tau `$ direction method, described in  and the references therein.

The $`\tau `$ selection efficiencies and the background fractions for the data, presented in Table 1, are estimated from the Monte Carlo simulation. Only statistical errors are quoted. In this data sample the only relevant contamination arises from $`\tau ^+\tau ^{}`$ events with misidentified decay modes ($`\tau `$ background).

### 4.2 Tau direction of flight

The reconstruction of the $`\tau `$ flight direction is mandatory in this analysis to access the event observables, which are functions of the four-momenta of the final state particles in the $`\tau `$ rest frame. This can be achieved in the semileptonic decays, for which the $`\tau `$ direction lies on a cone around the total hadron momentum. For events with both taus decaying semileptonically the $`\tau `$ direction lies along one of the intersection lines of the two reconstructed cones in the case that the taus are produced back-to-back, with equal energies given by $`\sqrt{s}/2`$, and $`m_{\nu _\tau }=m_{\overline{\nu }_\tau }=0`$. However, the two cones may not intersect due to detector effects or radiation. If the cones intersect, the event is considered twice using either solution. If the cones do not intersect, the particle momenta are fluctuated within their measurement errors and the event is accepted if the cones intersect in a minimum number of trials; the average direction is then used .

The effect of using both $`\tau `$ directions has been studied with a Monte Carlo sample having approximately the same size as the data. Table 2 presents the statistical errors obtained in this analysis and using the correct $`\tau `$ direction from the information at the generator level. Not distinguishing between the two $`\tau `$ directions induces some degradation in the overall sensitivity for the four anomalous weak dipole moments.

### 4.3 Candidates and efficiency matrix

The final selection for this analysis requires the $`\tau `$ direction to be successfully reconstructed, as described in Section 4.2, and the event observables ($`W_1`$, $`\mathrm{cos}\theta _{h_1}`$, $`\varphi _{h_1}`$, $`W_2`$, $`\mathrm{cos}\theta _{h_2}`$, $`\varphi _{h_2}`$) to lie in their domains of validity. These requirements decrease the number of candidates by 21%, the main reason being the inability to reconstruct the $`\tau `$ direction for some events.

The final number of candidates in each decay topology is given in Table 3. The efficiency matrix $`ϵ_{ij}`$ for $`i,j=\pi `$, $`\rho `$, $`\pi 2\pi ^0`$, $`3\pi `$ is calculated as a function of the generated polar angle $`\mathrm{cos}\theta _h^{(0)}`$ separately in the barrel and endcaps; the dependence of $`ϵ_{ij}`$ on $`\varphi _h^{(0)}`$ and $`W^{(0)}`$ is quite uniform and has been integrated out. Figure 3 shows the efficiencies for the barrel. The diagonal elements of this figure represent the identification of each decay mode, while the off-diagonal elements represent its misidentification.

### 4.4 Detector effects

The correct approach to introduce the detector effects in the likelihood formula would be via a smearing function $`T_{ij}`$ depending on 12 variables: the set of event observables and the corresponding generated values. The indices $`i`$ and $`j`$ indicate the generated and the reconstructed channel, respectively. With the notation used here, $`T_{ij}=T_{ij}(W_1,\mathrm{cos}\theta _{h_1},\varphi _{h_1},W_2,\mathrm{cos}\theta _{h_2},\varphi _{h_2},W_1^{(0)},\mathrm{cos}\theta _{h_1}^{(0)},\varphi _{h_1}^{(0)},W_2^{(0)},\mathrm{cos}\theta _{h_2}^{(0)},\varphi _{h_2}^{(0)})`$. Because this function cannot be easily calculated, the detector effects are parametrised by the factorised smearing functions $`D_{ij}(x,x^{(0)})`$ ($`x`$ = $`W`$, $`\mathrm{cos}\theta _h`$, $`\varphi _h`$ and $`i,j`$ = $`\pi `$, $`\rho `$, $`\pi 2\pi ^0`$, $`3\pi `$). For $`x`$ = $`W`$, $`\mathrm{cos}\theta _h`$, correlations are neglected, while for $`x`$ = $`\mathrm{cos}\theta _h`$, $`\varphi _h`$, the correlations are taken into account.

The functions $`D_{ij}(x,x^{(0)})`$ give the probability that for generated $`i`$ and reconstructed $`j`$ the smearing introduced by the detector is $`(xx^{(0)})`$ for a certain generated $`x^{(0)}`$ with reconstructed $`x`$. From the definition it follows that

$$D_{ij}(x,x^{(0)})𝑑x=1.$$
(15)
These functions are obtained with the SM Monte Carlo simulation by binning the $`(x,x^{(0)})`$ plane. The binning has been chosen small enough to correctly convolve detector effects with the generated distribution.

In the likelihood expression, Eq. 13, the detector effects are included by replacing the functions $`H_i^\mu `$ by

$`\stackrel{~}{H}_j^\mu (W,\mathrm{cos}\theta _h,\varphi _h)={\displaystyle \underset{i}{}}{\displaystyle H_i^\mu (W^{(0)},\mathrm{cos}\theta _h^{(0)},\varphi _h^{(0)})D_{ij}(W,W^{(0)})}`$
$`\text{ }\times D_{ij}(\mathrm{cos}\theta _h,\mathrm{cos}\theta _h^{(0)})D_{ij}(\varphi _h,\varphi _h^{(0)},\mathrm{cos}\theta _h)ϵ_{ij}(\mathrm{cos}\theta _h^{(0)})dW^{(0)}d\mathrm{cos}\theta _h^{(0)}d\varphi _h^{(0)}.`$
The sum runs over all modes $`i`$ which have been reconstructed as one of the modes $`j`$ used in the analysis, whereby all possible $`\tau `$ decay modes are included in $`i`$. The $`\tau `$ branching fractions are taken into account implicitly in Eq. 16 because the full differential cross section (Eq. 12) contains the probability of generating a $`\tau ^+\tau ^{}`$ pair decaying into specific decay modes with certain final state topologies.

In terms of the effective functions $`\stackrel{~}{H}_i^\mu `$ the likelihood for each event reads

$$L_{ij}=\underset{\mu ,\nu =0,\mathrm{},3}{}\overline{R}_{\mu \nu }(\mu _\tau ,d_\tau |\theta _\tau )\stackrel{~}{H}_i^\mu (W_1,\mathrm{cos}\theta _{h_1},\varphi _{h_1})\stackrel{~}{H}_j^\nu (W_2,\mathrm{cos}\theta _{h_2},\varphi _{h_2}).$$
(17)
### 4.5 Calibration curves

There are two sources for possible bias in the fitting procedure: ($`i`$) the detector effects are handled by the factorised smearing functions $`D_{ij}`$ above which take correlations into account only partially, and ($`ii`$) radiative corrections are not included in the likelihood. It is thus necessary to evaluate the adequacy of the fitting process. This is done with the SCOT Monte Carlo program , interfaced with TAUOLA for the $`\tau `$ decays and with the full detector simulation. The program SCOT describes $`e^+e^{}\tau ^+\tau ^{}`$ production at an energy around the $`Z`$ peak at tree level. It includes the anomalous weak dipole moments and all $`\tau `$ spin effects. The initial state radiation is included by adding a simple radiator function .

The checks are performed by generating various Monte Carlo samples with different values of the anomalous weak dipole moments. The couplings $`a_e`$, $`v_e`$, $`a_\tau `$, $`v_\tau `$ are set to their SM values. The anomalous weak dipole moments $`\mathrm{Re}\mu _\tau `$, $`\mathrm{Im}\mu _\tau `$, $`\mathrm{Re}d_\tau `$, $`\mathrm{Im}d_\tau `$ are varied one by one in an adequate region around zero. The dependence of the reconstructed values on the generated parameters is taken as linear. Significant deviations of the slopes from unity are found for certain decay topologies, and the offset for $`\mathrm{Re}\mu _\tau `$ is not consistent with zero for certain channels. These effects have been studied and are mostly related to not using the correct $`\tau `$ direction and to background effects. In this analysis, this calibration for each anomalous weak dipole moment and decay topology is taken into account to obtain the corresponding individual measurements. The slopes, offsets and $`\chi ^2`$ of the linear fits are presented in Table 4.

The offsets and slopes were derived with a Monte Carlo program which includes only a first order radiator for the initial state bremsstrahlung. The offsets have then been verified with the KORALZ Monte Carlo program  interfaced with the full detector simulation. KORALZ describes $`e^+e^{}\tau ^+\tau ^{}`$ production at an energy around the $`Z`$ peak and with initial state bremsstrahlung corrections up to $`O(\alpha ^2)`$, final state $`O(\alpha )`$ bremsstrahlung and $`O(\alpha )`$ electroweak corrections. However, this program contains only longitudinal spin effects (i.e., the production terms used are $`R_{00}`$, $`R_{03}`$ and $`R_{33}`$) and the anomalous weak dipole moments are set to zero. To check the effect of this approximation, a maximum likelihood fit was built to include the complete $`R_{00}`$, $`R_{03}`$ and $`R_{33}`$ terms and the anomalous parts of the other $`R_{\mu \nu }`$ terms. The values obtained for the four anomalous weak dipole moments for each decay topology are consistent with the corresponding offsets computed with SCOT and the first order radiator. Table 5 shows the global differences between the offsets of the central values of the four anomalous weak dipole moments when using either KORALZ or SCOT. The final results of this analysis are obtained by correcting the individual measurements according to the KORALZ offsets and the SCOT slopes and by including the statistical error of the correction in the systematic uncertainty. The statistical error of the SCOT offsets is also included in the systematic uncertainty as these offsets have contributions from all $`R_{\mu \nu }`$ terms.

## 5 Systematic uncertainties

The systematic uncertainty is calculated for each anomalous weak dipole moment and final state decay topology. The estimates are shown in Tables 6 to 9. The last row of these tables shows the combined systematic uncertainty from each source taking into account the correlations between channels.

The ECAL effects are related to the uncertainty in the global energy scale and the nonlinearity of the response. The global scale is known at the level of 0.25%, through the calibration with Bhabha events. Global variations have been applied to each ECAL module and the effect is propagated to the fitted parameters. The nonlinearity of the response is related to the wire saturation constants of ECAL, which are different in the barrel and endcaps. These saturation constants have been fluctuated within their nominal errors while keeping the measured energy fixed at $`M_\mathrm{Z}/2`$.

The TPC systematic errors are related to the momentum measurement: $`(i)`$ an effect due to the magnetic field acting similarly on positive and negative charged tracks, and $`(ii)`$ a sagitta effect affecting oppositely positive and negative tracks. These two effects are calibrated with dimuon events and the corresponding corrections are applied to the $`\tau `$ data. The systematic errors are then estimated by varying the corrections within their errors for each year.

The systematic errors in the column labeled “Align.” are due to a possible azimuthal tilt between the different parts of the detector.

Variations in the $`\tau `$ branching fractions are considered in the $`\tau `$BF column. This systematic uncertainty was determined from its effect on the calibration curves (Section 4.5).

The experimental errors on the weak parameters $`\mathrm{sin}^2\theta _W`$ and $`M_\mathrm{Z}`$ are propagated to the fitted values. Other weak parameters have negligible effect on the measurements. These effects are summarised in the column labeled “Wpar”.

The finite Monte Carlo statistics also causes systematic uncertainties. The most relevant statistical uncertainty is for the KORALZ offsets. The statistical error of the offsets and slopes obtained with SCOT and the first order radiator (Table 4) are also taken into account. Finally, the statistical error in the calculation of the efficiency matrix is also considered. These effects are shown under the column “MC st.”.

The $`a_1`$ decay dynamics are not well described theoretically. The impact of this was evaluated in the past  by implementing several models in the analysis . The implementation of those models is much more difficult in the present analysis. The uncertainty is estimated by means of three models: the Kühn & Santamaria (KS) model  (used in the fitting formula), the Feindt model  and the Isgur, Morningstar and Reader (simplified) model . The effects of the Feindt and IMR models on $`W`$ and $`h^3`$ are calculated. The corresponding ratios with the $`W`$ and the $`h^3`$ of the KS model are then used to scale the error.

Another source of systematic error is due to fake photons generated by hadron interactions in the ECAL or by electromagnetic fluctuations. This quantity of photon candidates is underestimated in the Monte Carlo simulation compared to the data. This deficit was originally observed in a substantial disagreement between the data and the Monte Carlo simulation for the $`W`$ distribution in the $`\pi 2\pi ^0`$ channel. This discrepancy has notably decreased after weighting the events of the Monte Carlo simulation according to the number of fake photons. This weighting was optimised for other $`\tau `$ analyses . Figure 1 compares the $`W`$ distributions for the data and the Monte Carlo after the approximate weighting. In the end, the effect of fake photons is taken into account by removing fake photons in the simulation and using the difference in the fitted parameters as systematic uncertainty.

The total systematic error and the statistical error for each channel and parameter are shown in the last two columns of the tables.

## 6 Results and conclusions

The final individual measurements of the four anomalous weak dipole moments are obtained applying the offsets and slopes described in Section 4.5. The results for the different decay topologies are presented in Figs. 4 and 5, including both the systematic and the statistical errors. All these measurements are consistent with the SM prediction. Figure 6 shows the relative weights of the different decay channels for the four measured anomalous weak dipole moments.

The final combined results on the four anomalous weak dipole moments are listed in Table 10, showing the statistical, systematic and total errors. The statistical correlations are given in Table 11. The final 95% CL upper limits derived from these measurements are presented in Table 12.

These results supersede the previous ALEPH measurement of $`\mathrm{Re}d_\tau `$ ; the measurement of $`\mathrm{Re}\mu _\tau `$, $`\mathrm{Im}\mu _\tau `$ and $`\mathrm{Im}d_\tau `$, presented in this paper gives the most stringent limits on these quantities to date.

## Acknowledgments

We wish to thank our colleagues from the accelerator divisions for the successful operation of LEP. We are indebted to the engineers and technicians in all our institutions for their contribution to the good performance of ALEPH. Those of us from non-member states thank CERN for its hospitality.


format: text

text:
# Perturbations to Stellar Structure in 2D: Stellar Rotation and Heating in X-ray Binaries.

## 1. Introduction.

Numerical study of stellar structure and evolution is now a well established subject with a long history. Today there is a plethora of 1D stellar evolution codes, from which astrophysicists can choose. All these models rely on the assumption that the star is spherically symmetric.

This assumption has served well for many years. However it is insufficient in cases that require stars to be non-symmetric. Stellar rotation, accretion and binary systems all add asymmetry to a spherical star. The required computational power is now available to look at these perturbations.

We have been working to develop a code that calculates 2D stellar structures in a non-computationally demanding way. So far we have studied perturbations to structure only but we plan to move into 2D stellar evolution. The two areas we have studied are rotation of solar-like stars and heating of a star by a companion star. The results are presented in this paper.

## 2. Numerical Method.

The method we use was developed by Garaud (2002). A 1D model in hydrostatic equilibrium is taken from a standard 1D evolution code and used as the underlying state for perturbation to the structure in 2D. The perturbative equations are linear with respect to thermodynamic quantities, with full non-linearity kept for the flow perturbations. The equations are

$$\rho \text{u}\text{u}=\stackrel{~}{p}\rho \stackrel{~}{\mathrm{\Phi }}\stackrel{~}{\rho }\mathrm{\Phi }+\nu ^2\text{u}+\frac{1}{3}\nu (\text{u}),$$
(1)
$$^2\stackrel{~}{\mathrm{\Phi }}=4\pi G\stackrel{~}{\rho },$$
(2)
$$\frac{\stackrel{~}{p}}{p}=\frac{\stackrel{~}{\rho }}{\rho }+\frac{\stackrel{~}{T}}{T},$$
(3)
$$(\rho \text{u})=0,$$
(4)
and

$$\rho T\text{u}s=k^2\stackrel{~}{T},$$
(5)
where u is fluid velocity, $`p`$ the pressure, $`\rho `$ the density, $`\mathrm{\Phi }`$ the gravitational potential, $`\nu `$ the kinematic viscosity, $`T`$ the temperature, $`k`$ the thermal conductivity and $`s`$ the entropy. The variables with a tilde are the perturbed variables and those without are taken from the 1D hydrostatic equilibrium model. We look for solutions to these equations that are stationary and axisymmetric. Note that we do not perturb the energy generation or opacities. The solution method is to use finite differences in the radial direction and expansion upon Chebyshev polynomials (essentially a Fourier decomposition) in the latitudinal direction.

Once this system is setup all that remains is to impose the boundary conditions appropriate to the system considered.

## 3. Stellar Rotation.

Stellar Rotation is an old problem (Eddington 1925; Sweet 1950). Commonly used models examine meridional currents in stars that are rotating slowly, that is more slowly than the Sun. There are two main effects of rotation on stellar evolution, the increase or decrease of the amount of contraction or expansion of a star during its evolution and the induction of extra mixing of material within the star, especially within regions that are convectively stable. This has been recently reviewed by Pinsonneault (1997).

Garaud (2002) studied the meridional currents in the radiative zone of solar-like stars set up by the differential rotation of an outer convective zone. To investigate the problem a radiative zone model was taken from a 1D standard solar model (Christensen-Dalsgaard et. al. 1991) and perturbed by forcing rotation on the outer boundary and the equations were solved within the radiative zone of the star, with the assumption that the convection is not affected by the radiative zone because the driving forces are larger than any effect from the radiative-zone flow.

The boundary condition imposed differential rotation at the surface of the radiative zone of the form

$$\mathrm{\Omega }_{\text{surface}}(\theta )=\underset{n}{}\mathrm{\Omega }_n\mathrm{cos}^n\theta $$
(6)
along with stress continuity across the surface and $`^2T=0`$ outside. For further details see Garaud (2002).

## 4. Heating of Main-Sequence Stars in X-ray Binaries.

It is possible to study different problems with the same method. Essentially only the boundary conditions require alteration. Therefore we adapted the code to study heating of a main-sequence star in a binary system by X-rays emitted by an accreting neutron-star companion. Many such X-ray binaries are known to exist.

Observational evidence that something non-1D is occurring comes from HZ Herculis where there are variations in the light curve, especially at the minimum when the X-ray source is eclipsed, that cannot be explained by a main-sequence star with a homogeneous surface temperature. It is generally assumed that the external radiation that falls on the star’s surface is thermalized and re-emitted, with some of the energy emitted from the unilluminated side. To achieve this there must be some form of circulation current in the star’s envelope.

There have been previous attempts to model this situation, Kippenhahn & Thomas (1978) assumed that the flux was thermalized very near to the surface and found that the induced flow penetrated to a depth of 20% of the star’s radius, with a 10% increase in the illuminated hemisphere temperature and 2% of the extra flux emitted from the unilluminated surface. The velocities of the vertical flows are subsonic, while the flows parallel to the surface are likely to be supersonic.

Later, Kırbıyık (1982) produced a similar model but with a much shallower penetration of only a few percent of the star’s radius. These analyses of the reflection effect have been used to explain many light curves in close binary systems, with the solution of Kippenhahn & Thomas (1978) being the most commonly citied.

### 4.1. Boundary Conditions

Careful thought must be given to formulation of the boundary conditions because there are many possibilities for boundary conditions that add extra energy to the surface of the star. The X-ray flux is absorbed in a thin layer near the surface. The column density required is of the order of $`0.11\mathrm{g}\mathrm{cm}^2`$. This is reached at a very shallow depth in the stellar surface, under the assumption that there is no absorption in the stellar atmosphere.

To a first approximation we can say that, at the surface, we are adding extra energy with a boundary condition

$$k\frac{T}{r}=\sigma _\mathrm{B}T^4\alpha f(\theta )\mathrm{cos}\theta .$$
(7)
The function $`f`$ is defined as

$$f(\theta )=\{\begin{array}{cc}1\hfill & \frac{\pi }{2}<\theta <\frac{\pi }{2},\hfill \\ 0\hfill & \text{otherwise},\hfill \end{array}$$
(8)
where $`\theta `$ is the latitude, $`\alpha `$ the strength of the X-ray flux, $`\alpha =L_X/(4\pi d^2)`$. Because we are dealing with small perturbations we expand $`TT+\stackrel{~}{T}`$ and assuming $`\stackrel{~}{T}`$ is small we obtain

$$k\frac{\stackrel{~}{T}}{r}=4\sigma _\mathrm{B}T^3\stackrel{~}{T}\alpha f(\theta )\mathrm{cos}\theta .$$
(9)
This boundary condition can also be used if the illuminating star is a main-sequence star.

We also have conditions on the flow such that

$$u_r(r=1)=0\mathrm{and}u_\theta (\mu =\pm 1)=0,$$
(10)
along with a stress free boundary on the flow at the surface. We have defined $`\mu =\mathrm{cos}\theta `$, r the dimensionless radial coordinate, $`r=0`$ is at the centre and $`r=1`$ is at the surface.

### 4.2. Numerical Results

Only one model, the solar radiative zone, has been studied successfully numerically. By this we mean that we have removed the upper part of the Sun so we are directly heating the region at $`0.6R_{}`$. This was done to test the numerics on a well-known background state and helped us develop analytical tools that can be generalised to any stellar radiative zone. The viscosity has also been boosted to resolve the thin thermo-viscous boundary layers.

We find that circulation is induced but limited to a thin boundary layer near the surface. Similarly to previous authors we find that there are two circulation cells. Figure 1 shows the temperature fluctuations versus depth. The perturbations only penetrate to a shallow depth. Figure 2 plots the value of $`\psi `$ which represents the nature of the induced flows. Radial velocities are subsonic but latitudinal velocities do reach supersonic values for illumination equivalent to a source of $`10^{35}\mathrm{erg}\mathrm{s}^1`$ at $`1`$ au. The flow velocities are linearly dependent on the illuminating flux.

Because the effect of the heating is confined to a small boundary layer we realised that we could look at the problem analytically. This provides us with some generic analytical results on the heating that we can use for any star and any form of heating.

### 4.3. Analytic Modelling

We have investigated the situation and solved the equations analytically using boundary layer approximations to keep only the highest order radial derivatives. We use the same method of Fourier expansion of the variables. We simplify the equations by restricting them to the surface layer and using only the dominant terms in that region. We start by solving the mass conservation equation for incompressible flow

$$(\rho 𝐮)=0.$$
(11)
We assume the density variation is negligible and define a circulation variable $`\psi `$ such that

$$\rho u_r=\frac{1}{r^2}\frac{\psi }{\mu }\mathrm{and}\rho u_\theta =\frac{1}{r\sqrt{1\mu ^2}}\frac{\psi }{r}.$$
(12)
We now make the variables dimensionless by scaling distance to the star’s radius, $`R`$, and velocities using the viscous diffusion timescale so that $`\psi ^{}=\psi R\nu `$ becomes dimensionless. Thus the energy equation (5) and vorticity equation (1) become

$$(1\mu ^2)\frac{g}{T}\frac{\stackrel{~}{T}}{\mu }=\frac{\nu ^2}{R^3}\frac{^4\psi ^{}}{r^4}$$
(13)
and

$$\frac{\nu R\rho c_pN^2T}{kg}\frac{\psi ^{}}{\mu }=^2\stackrel{~}{T}.$$
(14)
where $`g`$ is the acceleration due to gravity, $`c_p`$ the specific heat capacity at constant pressure and $`N^2`$ the buoyancy frequency.

To solve the equations, we expand the variables $`\stackrel{~}{T}`$ and $`\psi `$ as

$$\stackrel{~}{T}T_0(r)+\mu T_1(r)+(2\mu ^21)T_2(r)$$
(15)
and

$$\frac{\psi ^{}}{\mu }=u_0(r)+\mu u_1(r)+(2\mu ^21)u_2(r).$$
(16)
To zeroth order the system reduces to Laplace’s equation for the temperature perturbation. For the higher orders we use equations (12)–(16) to obtain equations of the form

$$T_i=\delta _i^6\frac{^6T_i}{r^6},$$
(17)
when $`i>0`$. These have a solution of the form

$$T_i=𝒯_i\mathrm{exp}\left(\delta _i(r1)\right),$$
(18)
with similar equations for the flow velocities. The value of $`\delta `$ is of most interest because it determines the depth of penetration below the surface. With this analysis we obtain

$$\delta _1=\left(\frac{2\rho c_pN^2R^4}{\nu k}\right)^{\frac{1}{6}}\mathrm{and}\delta _2=\left(\frac{6\rho c_pN^2R^4}{\nu k}\right)^{\frac{1}{6}}.$$
(19)
This gives us six solutions for $`\delta _i`$ in the complex plane and we use those that decay with depth. We find the full solution becomes

$$T_i(r)=𝒯_i\left(\mathrm{exp}\left(\delta _i(r1)\right)+2\mathrm{exp}\left(\frac{1}{2}\delta _i(r1)\right)\mathrm{cos}\left(\frac{\sqrt{3}}{2}\delta _i(r1)\right)\right).$$
(20)
Then by fitting to the boundary condition (9) we obtain

$$𝒯_i=A_i\alpha \left(12\sigma _\mathrm{B}T^3+\frac{2k\delta _i}{R}\right)^1.$$
(21)
Orders only differ by a numerical factor, for expansion to second order $`A_1=\frac{1}{2}`$ and $`A_2=\frac{2}{3\pi }`$. The solutions for the flow components $`u_r`$ and $`u_\theta `$ can be deduced from equations (13) and (21). Table 1 shows values of $`\delta _1`$ for various stellar models, computed with the Eggleton code (see Pols et al. 1996). The stars are above $`2M_{}`$ to ensure fully radiative envelopes. In the numerical study we had to boost the viscosity by a few orders of magnitude to resolve the boundary layer. But in the analytic work we use real values for the viscosity (Spitzer 1962) and so the penetration is shallower, varying as $`\delta \nu ^{\frac{1}{6}}`$. However if we use the same values as for the numerical calculations the results agree closely. See figure 1 for the comparison of results by this method and those from numerical results.

The temperature perturbations are also listed in table 1, with example flow velocities in table 2. The sources are placed at a distance of $`1`$au from the secondary. This is much larger than the typical separation in an X-ray binary (a few $`R_{}`$ but we cannot yet deal with such intense irradiation – see Discussion). For one of the stellar models we show the variation of the values with different illuminating fluxes.

### 4.4. Discussion

From the tables it can be seen that a circulation is present even if the irradiating flux is low, although the perturbation is too small to be observed. With larger fluxes a problem arises because the latitudinal flow becomes supersonic. This agrees with the results of Kippenhahn and Thomas but would produce shocks in the flow for which we have not allowed in this work.

The amount of extra flux emitted from the unilluminated side of the star can be estimated by comparison of the strength of the dipole mode to the other terms. We find that about $`60\%`$ of the flux is emitted from the irradiated surface with $`40\%`$ from the unilluminated side. This is a smaller difference than found by previous authors and is due to the large zeroth order temperature perturbation. The value is large because we have not allowed the star to expand.

There are two interesting deviations from the work of Kippenhahn and Thomas. First we have much shallower penetration similar to the results of Kırbıyık and second there is more flux emitted from the opposite side. In our next step we shall test these analytic results for various stars numerically.

In these results it should also be noted that the boundary conditions are simple. While they are a reasonable estimate for illumination by black-body-like stars, if the bulk of the incident radiation is in the form of X-rays then energy will not all be deposited at the surface of the star but at a greater depth.

## 5. Future Work

The code calculates 2D perturbations to the hydrostatic equilibrium of any stellar radiative region. However it can only deal with linear perturbations to the thermodynamic variables. Therefore the next phase is to develop the code to model non-linear effects and solve for the full structure, not just perturbations to an assumed hydrostatic equilibrium. In particular in the above study, we have not allowed the star to expand or shrink when heated and we shall need to allow deformation, deviations from spherical shape. We shall also improve upon the boundary conditions to include a model atmosphere as in Tout et al. (1989) and to include external radiation pressure from the irradiating flux, along with the possible deeper penetration of the incident radiation. To study less massive stars we shall also be looking at expanding the code to deal with convection zones by some parameterisation that does not require us to follow the full dynamics of convection.

## References

Christensen-Dalsgaard, J., Gough, D. O. & Thompson M. T. 1991, ApJ, 378, 866

Eddington, A. S. 1925, Observatory, 48, 73

Garaud, P. 2002, MNRAS, 335, 707

Kippenhahn, R., & Thomas, H. C. 1978, A&A, 63, 265

Kirbiyik, H. 1982, MNRAS, 200, 907

Pinsonneault, M. H. 1997, ARA&A, 38, 143

Spitzer, L. 1962, Physics of Fully Ionized Gases (New York, London : Interscience Publishers)

Sweet, E. 1950, MNRAS, 110, 548

Tout, C. A., Eggleton, P. P., Fabian, A. C., & Pringle, J. E. 1989, MNRAS, 238, 427

Tout, C. A., Pols, O. R., Eggleton, P. P. & Han, Z. 1996, MNRAS, 281, 257


text:
# Determining the BDMPS transport coefficient via medium-modified fragmentation functions

## Abstract

In nucleus-nucleus collisions at RHIC and LHC, partons produced at high transverse momentum can undergo multiple scattering within the collision region prior to fragmenting into hadrons. We have studied the resulting medium-modified fragmentation function based on a calculation of the BDMPS-Z medium-induced gluon radiation for a dense, expanding medium of small finite extension. Here we explain how the BDMPS transport coefficient $`\widehat{q}`$ which measures the energy density attained in the collision, can be extracted from the observed modification of high-$`p_{}`$ hadroproduction. We also comment on the significant remaining uncertainties in extracting $`\widehat{q}`$ from data.

To leading order in perturbative QCD, high-$`p_{}`$ hadroproduction in proton–proton collisions is described by the factorization formula

$$\frac{d\sigma ^h}{dp_t^2dy}=K(\sqrt{s})\frac{dz}{z^2}𝑑y_2\underset{i,j}{}x_1f_i^A(x_1,Q^2)x_2f_j^B(x_2,Q^2)\frac{d\sigma ^{ijk}}{d\widehat{t}}D_{kh}(z,Q^2).$$
(1)
Here $`f_i^A`$, $`f_j^B`$ are the proton parton distribution functions (PDF), $`d\sigma ^{ijk}/d\widehat{t}`$ is the partonic cross section and $`D_{kh}(z,Q^2)`$ describes the fragmentation of a parton $`k`$ into the hadron $`h`$ carrying a fraction $`z`$ of the momentum. Eq. (1) leads to a fair description of the shape of high-$`p_{}`$ hadronic spectra while the normalisation has to be adjusted by an energy-dependent $`K`$-factor (see e.g. ). Also to NLO, the disagreement between theory and experiment lies essentially in an albeit reduced normalisation factor .

For nucleus–nucleus collisions, additional corrections arise due to parton density effects. At sufficiently high $`p_{}`$ where the factorised form (1) is expected to hold, two corrections are known: i) Density effects in the initial state can be parametrised by nuclear PDFs for which global fits are available . At RHIC energies, however, this is only a small correction  . ii) Density effects in the final state affect the fragmentation of the high $`p_{}`$ parton produced in the hard collision. Several studies indicate that the dominant final state correction is due to medium-induced gluon radiation and the ensuing parton energy loss. If the high $`p_{}`$ parton looses with probability $`P_E(ϵ)`$ a fraction $`ϵ`$ of its energy while propagating through the medium, then its medium-modified fragmentation function can be written as

$`D_{kh}^{(\mathrm{med})}(z,Q^2)={\displaystyle _0^1}𝑑ϵP_E(ϵ){\displaystyle \frac{1}{1ϵ}}D_{kh}({\displaystyle \frac{z}{1ϵ}},Q^2).`$ (2)
We have calculated the quenching weights $`P_E(ϵ)`$ for quarks and gluons starting from the BDMPS-Z  gluon radiation spectrum $`dI/d\omega `$ generalised to the case of an expanding medium of small finite size . For the case of a static medium, the resulting finite size corrections to the BDMPS-expression of $`dI/d\omega `$ are shown in Fig. 1. They become important for gluon energies much smaller than the characteristic gluon frequency $`\omega _c=\widehat{q}L`$. The latter depends on the in-medium path-length $`L`$ and the BDMPS transport coefficient $`\widehat{q}`$ which measures the amount of transverse momentum squared transferred from the medium to the hard parton per unit path-length. For a static medium, the radiation spectrum depends only on $`\omega _c`$ and the dimensionless parameter combination $`R=\frac{1}{2}\widehat{q}L^3`$. As seen in Fig. 1, the BDMPS case is the limit $`L\mathrm{}`$ of our calculation keeping $`\omega _c`$ finite.

Quenching weights are evaluated in the independent gluon emission approximation  and show both a discrete and a continuous contribution (see Fig. 2),

$`P_E(ϵ)`$ $`=`$ $`{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{1}{n!}}\left[{\displaystyle \underset{i=1}{\overset{n}{}}}{\displaystyle 𝑑\omega _i\frac{dI(\omega _i)}{d\omega }}\right]\delta \left(ϵ{\displaystyle \underset{i=1}{\overset{n}{}}}{\displaystyle \frac{\omega _i}{E}}\right)\mathrm{exp}\left[{\displaystyle 𝑑\omega \frac{dI}{d\omega }}\right]`$ (3)
$`=`$ $`p_0\delta (ϵ)+p(ϵ).`$ (4)
Here, $`P(ϵ)`$ is a generalised probability of norm 1. For a medium of finite length $`L`$, there is a finite probability $`p_0`$ that the fragmentation is not disturbed by the medium. This discrete contribution decreases with increasing in medium path-length or increasing density of the medium.

The transport coefficient is proportional to the density of scattering centres, $`\widehat{q}(\tau )n(\tau )`$. The dynamical expansion of the collision region can be parametrised by a power law decrease in proper time, $`n(\tau )=n_0\left(\frac{\tau _0}{\tau }\right)^\alpha `$, where $`\alpha =0`$ for static, $`\alpha =1`$ for Bjorken expansion. Fig. 3 shows the gluon radiation spectrum for different values of $`\alpha `$ with parameters $`\omega _c`$ and $`R`$ written in terms of the linearly weighted transport coefficient

$`\overline{\widehat{q}}={\displaystyle \frac{2}{L^2}}{\displaystyle _{\tau _0}^{\tau _0+L}}𝑑\tau (\tau \tau _0)\widehat{q}(\tau ).`$ (5)
After rescaling according to (5), the radiation spectra $`dI/d\omega `$ calculated for different values $`\alpha `$ coincide (see Fig. 3). This establishes a dynamical scaling law  which simplifies calculations considerably since results for a dynamically expanding medium can be obtained by scaling the result obtained for a static medium. It extends earlier findings of dynamical scaling of the average energy loss  to the $`\omega `$-differential radiation spectrum .

Based on the quenching weights (3) for static and expanding media, we have calculated medium-modified fragmentation functions according to (2), using the KKP parametrisation  for parton fragmentation in the vacuum (see Fig. 4). The virtuality $`Q`$ of the parent parton is identified with its initial transverse energy. For a medium showing one-dimensional Bjorken expansion ($`\alpha =1`$), the transport coefficient $`\overline{\widehat{q}}`$ is related to the gluon rapidity density as

$`R={\displaystyle \frac{1}{2}}\overline{\widehat{q}}L^3={\displaystyle \frac{L^2}{R_A^2}}{\displaystyle \frac{dN^g}{dy}}.`$ (6)
The calculation presented in Fig. 4 allows to determine medium modifications to the transverse momentum spectrum (1) at sufficiently high $`p_{}`$. While this last step is not yet performed, one can estimate the effect based on the observation  that the hard partonic cross section in (1) essentially weighs the fragmentation function with the power $`z^6`$. Inspecting the corresponding unintegrated moment $`z^6D_{\pi /u}(z,Q^2p_{}^2)`$ in Fig. 4, one estimates that a factor 2 suppression of the hadronic pion spectrum for $`p_{}`$ 6–7 GeV corresponds to an initial gluon density of 500–1000. For an in-medium path-length of $`L=7`$ fm, this amounts to a line-averaged transport coefficient $`\overline{\widehat{q}}=`$ (300 MeV)2 – (600 MeV)2/fm which corresponds to $`p_{}^2/L_{\mathrm{med}}`$ of several GeV2/fm within the first two fm/c.

For a reliable determination of $`\widehat{q}`$ from data, the evaluation of (1) is not sufficient. One also has to discriminate the medium modification discussed here from other corrections to (1), such as “soft” hadronic contributions at low $`p_{}`$ and higher order perturbative contributions at high $`p_{}`$. Much of this can be done by studying the dependence of spectra on the in-medium path-length e.g. via the centrality dependence of high-$`p_{}`$ hadroproduction and monojet/dijet production. Another discriminatory tool results from the fact that medium-induced gluon radiation is a medium-enhanced power correction which shows a power law decrease with increasing $`p_{}`$. This is in contrast e.g. to perturbative NLO corrections which mainly affect the normalisation but not the shape of spectra. This will also provide an important additional cross check at LHC: according to (6), if the rapidity density increases by a factor 5 from RHIC to LHC, and if the suppression factor 2 for RHIC $`\pi ^0`$ spectra  at $`p_{}`$ 6–7 GeV can be attributed solely to medium-induced gluon radiation, then a factor 2 suppression of the corresponding spectrum at LHC is expected in the range $`p_{}`$ 40–50 GeV. At even higher $`p_{}`$ still accessible at LHC, this suppression would decrease according to a power law, thus leading to a modified $`p_{}`$-shape of hadronic spectra.

C.A.S. is supported by a Marie Curie Fellowship of the European Community programme TMR (Training and Mobility of Researchers), under the contract number HPMF-CT-2000-01025.


text:
# Geometric Random Inner Products: A New Family of Tests for Random Number Generators

## I Introduction

Monte Carlo methods are among the most widely used numerical algorithms in computational science and engineering Dongarra and Sullivan (2000). The key element in a Monte Carlo calculation is the generation of random numbers. Although a truly random number sequence produced by either a physical process such as nuclear decay, an electronic device etc., or by a computer algorithm, may not actually exist, a new and computationally easy-to-implement scheme to investigate random number generators is always highly desirable.

There have been many proposed schemes for the quality measure of random number generators Knuth (1998); W. Press, S. Teukolsky, W. Vetterling, and B. Flannery (1992); Marsaglia (1985); Garcia (2000); Giordano (1997); Gentle (1998); A. M. Ferrenberg, D. P. Landau, and Y. J. Wong (1992); I. Vattulainen, T. Ala-Nissila, and K. Kankaala (1994). These computational tests are based either on probability theory and statistical methods (for example: the $`\chi ^2`$ test, the Smirnov-Kolmogorov test, the correlation test, the spectral test, and the DieHard battery of randomness tests), or on mathematical modeling and simulation for physical systems (for example: random walks and Ising model simulations). These methods also open the door to studying the properties of random number sequences such as randomness and complexity Kolmogorov and Uspenskii (1987). Some important attempts at an operational definition of randomness were previously developed by Kolmogorov and Chaitin (algorithmic informational theory) Kolmogorov (1998); Chaitin (1987, 1990, 2001) and by Pincus (approximate entropy) Pincus (1991).

In this paper, we study a new method to measure $`n`$-dimensional randomness which we denote by GRIP (Geometric Random Inner Products). The GRIP family of tests is based on the observation that the average scalar products of random vectors produced in geometric objects (e.g., circles and spheres), define geometric constants which can be used to evaluate the quality of random number generators. After presenting the simplest example of a GRIP test, we exhibit a computational method for implementing GRIP, which is then used to analyze a number of random number generators. We then discuss the GRIP formalism in detail and show how a random number sequence, when converted to random points in a space defined by a geometric object, can produce a series of known geometric constants. Later we introduce additional members and include them within the GRIP family. We then present the computational results for configurations of four, six, and eight random points, along with a consideration of some key issues. Finally, we conclude by discussing how the GRIP test measures the quality of random number generators by explicitly adding a new quantitative property to random number sequences along with the three known qualitative properties summarized in Ref. Kolmogorov and Uspenskii (1987).

## II General Description of the GRIP Formalism

The GRIP scheme is derived from the theory of random distance distribution for spherical objects, and can be generalized to other geometric objects with arbitrary densities Tu (2001); Tu and Fischbach (2002). First, three random points ($`\stackrel{}{r}_1`$, $`\stackrel{}{r}_2`$, and $`\stackrel{}{r}_3`$) are independently produced from the sample space defined by a geometric object. We then evaluate the average inner product of $`\stackrel{}{r}_{12}\stackrel{}{r}_{23}`$ from two associated random vectors, $`\stackrel{}{r}_{12}=\stackrel{}{r}_2\stackrel{}{r}_1`$ and $`\stackrel{}{r}_{23}=\stackrel{}{r}_3\stackrel{}{r}_2`$. For a geometric object such as an $`n`$-ball of uniform density with a radius $`R`$, the analytical result is a geometric constant which can be expressed in terms of the dimensionality $`n`$ of the space Tu (2001); Tu and Fischbach (2002):

$$\stackrel{}{r}_{12}\stackrel{}{r}_{23}_n=\frac{n}{n+2}R^2.$$
(1)
A simple derivation of Eq. (1) can be found in the Appendix.

The following procedures are the numerical implementation of our testing programs. A random number sequence produced from a random number generator is used to generate a series of three random points $`\stackrel{}{r}_1`$, $`\stackrel{}{r}_2`$, and $`\stackrel{}{r}_3`$ such that these random points are uniformly distributed in an $`n`$-dimensional spherical ball $`𝐁`$ of radius $`R`$, where

$$𝐁=\{(x_1,x_2,\mathrm{},x_n):x_1^2+x_2^2+\mathrm{}+x_n^2R^2\}.$$
(2)
We then compute a series of values for $`\stackrel{}{r}_{12}\stackrel{}{r}_{23}`$. If $`\stackrel{}{r}_{12}\stackrel{}{r}_{23}`$ is evaluated $`N`$ times, then statistically we expect

$$\underset{N\mathrm{}}{lim}\frac{1}{N}\underset{i=1}{\overset{N}{}}\left(\stackrel{}{r}_{12}\stackrel{}{r}_{23}\right)_i=\frac{n}{n+2}R^2,$$
(3)
as predicted by Eq. (1).

## III Random Number Generators

We now apply the GRIP test to the following random number generators frequently used in Monte Carlo simulations.

1. RAN0 - a linear congruential generator Knuth (1998); W. Press, S. Teukolsky, W. Vetterling, and B. Flannery (1992):

$$x_n=16807\times x_{n1}\mathrm{mod}2147483647.$$
(4)
2. RAN3 - a lagged Fibonacci generator Knuth (1998); W. Press, S. Teukolsky, W. Vetterling, and B. Flannery (1992):

$$x_n=\left(x_{n55}x_{n24}\right)\mathrm{mod}2^{31}.$$
(5)
3. R31 - a generalized feedback shift register (GFSR) generator Knuth (1998); A. M. Ferrenberg, D. P. Landau, and Y. J. Wong (1992); I. Vattulainen, T. Ala-Nissila, and K. Kankaala (1994):

$$x_n=x_{n31}x_{n3},$$
(6)
where $``$ is the bitwise exclusive OR operator.
4. durxor - a generator selected from IBM ESSL (Engineering and Scientific Subroutine Library) ess .
5. durand - a generator selected from IBM ESSL (Engineering and Scientific Subroutine Library) and the sequence period of durand is shorter than durxor ess .
6. ran\_gen - one of the subroutines in IMSL libraries from Visual Numeric ims .
7. Random - a Fortran 90/95 standard intrinsic random number generator Metcalf and Reid (1999).
8. Weyl - a Weyl sequence generator B. L. Holian, O. E. Percus, T. T. Warnock, and P. A. Whitlock (1994); Tretiakov and Wojciechowski (1999):

$$x_n=\left\{n\alpha \right\},$$
(7)
where $`\left\{x\right\}`$ is the fractional part of $`x`$, and $`\alpha `$ is an irrational number such as $`\sqrt{2}`$.
9. NWS - a nested Weyl sequence generator B. L. Holian, O. E. Percus, T. T. Warnock, and P. A. Whitlock (1994); Tretiakov and Wojciechowski (1999):

$$x_n=\left\{n\left\{n\alpha \right\}\right\}.$$
(8)
10. SNWS - a shuffled nested Weyl sequence generator B. L. Holian, O. E. Percus, T. T. Warnock, and P. A. Whitlock (1994); Tretiakov and Wojciechowski (1999):

$`s_n`$ $`=`$ $`M\left\{n\left\{n\alpha \right\}\right\}+{\displaystyle \frac{1}{2}},`$ (9)
$`x_n`$ $`=`$ $`\left\{s_n\left\{s_n\alpha \right\}\right\},`$ (10)
where $`M`$ is a large positive integer.

The computational results obtained from Eq. (3) when $`n=3`$ and $`n=9`$ are presented in Table 1. Results for random number generators based on other algorithms such as the Data Encryption Standard (DES) Knuth (1998); W. Press, S. Teukolsky, W. Vetterling, and B. Flannery (1992) can be found in Ref. Tu and Fischbach  along with the computed results obtained from other geometric objects. We note that both the ran\_gen and RAN0 generators perform better overall, while the NWS and Weyl generators (which are based on the Weyl sequence method) are ranked lowest compared to the other generators. The reasons why this is the case will be discussed later.

## IV GRIP Analysis

In the following, we analyze the relationship between GRIP and a random number sequence, and show how a good random number sequence, when converted to random points in a a space defined by a geometric object, can produce a series of known $`n`$-dimensional geometric constants. A random number sequence generated from a random number generator can be written as,

$$a_1a_2a_3a_4a_5a_6a_7a_8a_9a_{10}\mathrm{}\mathrm{}.$$
(11)
When the sequence is converted to represent random points in a $`2`$-dimensional geometric object, the random numbers in Eq. (11) can then be grouped in pairs as

$$\left(a_1a_2\right)\left(a_3a_4\right)\left(a_5a_6\right)a_7a_8a_9a_{10}\mathrm{}\mathrm{},$$
(12)
where Cartesian coordinates are used. The first set of random points $`\{\stackrel{}{r}_1,\stackrel{}{r}_2,\stackrel{}{r}_3\}`$ can thus be identified as

$$\stackrel{}{r}_1=(a_1,a_2),\stackrel{}{r}_2=(a_3,a_4),\stackrel{}{r}_3=(a_5,a_6).$$
(13)
GRIP then uses $`\stackrel{}{r}_1`$, $`\stackrel{}{r}_2`$, and $`\stackrel{}{r}_3`$ to evaluate the average scalar product which can be computed by rewriting,

$$\stackrel{}{r}_{12}\stackrel{}{r}_{23}=\frac{1}{N}\underset{i=1}{\overset{N}{}}\underset{j=1}{\overset{2}{}}\left(a_{6i4+j}a_{6i6+j}\right)\left(a_{6i2+j}a_{6i4+j}\right),$$
(14)
where $`N`$ is a large positive integer. When the geometric object is a circle of radius $`R`$ and uniform density, we expect $`\stackrel{}{r}_{12}\stackrel{}{r}_{23}0.5R^2`$ as predicted by Eq. (1).

The analysis for $`2`$-dimensional GRIP can be immediately generalized to the $`n`$-dimensional case. When the sequence in Eq. (11) is used to generate random points in a $`n`$-dimensional spherical object, we can regroup Eq. (11) as follows:

$$\left(a_1\mathrm{}a_k\right)\left(a_{k+1}\mathrm{}a_{2k}\right)\left(a_{2k+1}\mathrm{}a_{3k}\right)\left(\mathrm{}\right)\left(\mathrm{}\right)\left(\mathrm{}\right)\mathrm{}\mathrm{}.$$
(15)
The average scalar product of $`\stackrel{}{r}_{12}\stackrel{}{r}_{23}`$ can then be expressed as

$$\stackrel{}{r}_{12}\stackrel{}{r}_{23}=\frac{1}{N}\underset{i=1}{\overset{N}{}}\underset{j=1}{\overset{n}{}}\left(a_{3in2n+j}a_{3in3n+j}\right)\left(a_{3inn+j}a_{3in2n+j}\right).$$
(16)
When the geometric object is an $`n`$-ball with a radius $`R=1`$ and a uniform density, we expect from Eq. (1) that the result of Eq. (16) should be a geometric constant, $`n/(n+2)`$.

## V GRIP Members

For practical computational purposes, we may wish to transform a random number sequence from a uniform density distribution to one which is non-uniform. One of the most important non-uniform density distributions is the Gaussian (normal) distribution $`P(r)`$ with mean zero and standard deviation $`\sigma `$,

$$P(r)=\frac{1}{(2\pi )^{n/2}\sigma ^n}e^{(1/2)\left(r^2/\sigma ^2\right)}.$$
(17)
Here $`_0^{\mathrm{}}P(r)𝑑r=1`$, $`r=\left(x_1^2+\mathrm{}+x_n^2\right)^{1/2}`$, and $`n`$ is the space dimensionality. One can use either the Box-Muller transformation method to generate a random number sequence with a Gaussian density distribution, or use available subroutines from major computational scientific libraries such as IBM ESSL and IMSL ess ; ims . By applying the probability density function of the random distance distribution as discussed in Ref. Tu and Fischbach (2002), one can add a new GRIP member to investigate the quality of a Gaussian random number generator, and this new GRIP test can be expressed as:

$$\stackrel{}{r}_{12}\stackrel{}{r}_{23}_n=n\sigma ^2.$$
(18)
A very common situation arises when one has to produce random points uniformly distributed on the surface of an $`n`$-sphere of radius $`R`$. Some general computational techniques for doing this are summarized in Refs. Knuth (1998); Tu (2001). We can then use

$$\stackrel{}{r}_{12}\stackrel{}{r}_{23}_n=R^2,$$
(19)
to examine the quality of such transformed random number generators as discussed in Ref. Tu and Fischbach .

Another application of the GRIP formalism is in stochastic geometry. We can design a test scheme for a configuration utilizing any number of random points Tu and Fischbach , and these tests can be included in the GRIP family. Among the tests are:

1. Four uniform random points configuration for an $`n`$-ball of radius $`R`$

$`\left(\stackrel{}{r}_{12}\stackrel{}{r}_{23}\right)\left(\stackrel{}{r}_{34}\stackrel{}{r}_{41}\right)_n`$ $`=`$ $`{\displaystyle \frac{n\left(n+1\right)}{\left(n+2\right)^2}}R^4,`$ (20)
$`\left(\stackrel{}{r}_{12}\stackrel{}{r}_{34}\right)\left(\stackrel{}{r}_{23}\stackrel{}{r}_{41}\right)_n`$ $`=`$ $`{\displaystyle \frac{2n}{\left(n+2\right)^2}}R^4,`$ (21)
$`\stackrel{}{r}_{13}\stackrel{}{r}_{24}_n`$ $`=`$ $`0.`$ (22)
2. $`2m`$ uniform random points configuration for an $`n`$-ball of radius $`R`$

$$\left(\stackrel{}{r}_{12}\stackrel{}{r}_{23}\right)\mathrm{}\mathrm{}\left(\stackrel{}{r}_{2m\mathrm{1\hspace{0.17em}\hspace{0.17em}2}m}\stackrel{}{r}_{2m\mathrm{\hspace{0.17em}\hspace{0.17em}1}}\right)_n=(1)^m\frac{n\left(n^{m1}+1\right)}{\left(n+2\right)^m}R^{2m},$$
(23)
where $`2m`$ ($`m=2,`$ $`3`$, $`4`$ etc.) is a positive even number.

A derivation of Eq. (20) can be found in the Appendix. We summarize the computational results for Eq. (23) when $`m=2,`$ $`3,`$ $`4`$ in in Tables 2, 3, and 4. A discussion of other results, such as Eqs. (18) and (19), can be found in Ref. Tu and Fischbach .

We observe that all of the generators except NWS and Weyl perform significantly better in $`n=3`$ than in $`n=9`$ using the GRIP test based on $`\left(\stackrel{}{r}_{12}\stackrel{}{r}_{23}\right)\left(\stackrel{}{r}_{34}\stackrel{}{r}_{41}\right)_n`$. We also note from Table 2, and the $`n=9`$ results (from R31 to RAN0), that these results are clearly biased to larger numbers compared to the expected value. One interpretation may be that $`\left(\stackrel{}{r}_{12}\stackrel{}{r}_{23}\right)\left(\stackrel{}{r}_{34}\stackrel{}{r}_{41}\right)_9`$ is a more sensitive and dedicated computational test for investigating random number generators than other GRIP tests. We also note that the results for $`n=9`$ are overall worse than $`n=3`$, and that the results for $`\left(\stackrel{}{r}_{12}\stackrel{}{r}_{23}\right)\left(\stackrel{}{r}_{34}\stackrel{}{r}_{41}\right)_9`$ reveal a more significant bias than in any of the other cases. These results suggest that the GRIP test either in higher dimensions (large $`n`$), or using a configuration of four random points, can serve as a more computationally sensitive test to detect non-random patterns hidden in random number sequences. Finally we note that it is not surprising that the NWS and Weyl generators are ranked worst among all cases in our GRIP test. As reported previously in Tretiakov and Wojciechowski (1999), these two show unacceptable non-random behavior and strong correlations.

## VI Conclusions

We have presented a new computational paradigm for evaluating the quality of random number generators. We demonstrate how GRIP helps to understand complexity and randomness by adding a new property, besides three known properties (typical, chaotic, and the stability of frequencies) Kolmogorov and Uspenskii (1987), for random number sequences. This quantitative feature shows how a random number sequence, when converted to random points in a space defined by a geometric object, can produce a series of known geometric constants. Ten random number generators were selected to run our GRIP tests, and they are ranked based on the errors between the numerical and analytical results. Finally we note that one implication of our work is that computational scientists should test the random number generators they use in their simulations, and verify that their random number generators pass as many proposed tests as possible.

*

## Appendix A Derivation of $`\left(\stackrel{}{r}_{12}\stackrel{}{r}_{23}\right)_n`$ and $`\left(\stackrel{}{r}_{12}\stackrel{}{r}_{23}\right)\left(\stackrel{}{r}_{34}\stackrel{}{r}_{41}\right)_n`$

We derive the analytical result of Eq. (1) for a circle ($`n=2`$) of radius $`R`$ and uniform density. The same derivation can be applied to the case of $`n`$ dimensions where $`n3`$. We label three independent random points as $`1`$, $`2`$, and $`3`$ in Fig. 1, and then calculate

$$\stackrel{}{r}_{12}\stackrel{}{r}_{23}=r_{12}r_{23}\mathrm{cos}\theta =r_{12}r_{23}\mathrm{cos}\alpha ,$$
(24)
where $`\alpha +\theta =\pi `$. From the triangle formed by the random points, we then have

$$r_{31}^2=r_{12}^2+r_{23}^22r_{12}r_{23}\mathrm{cos}\alpha .$$
(25)
Extending this $`2`$dimensional case to the $`n`$-dimensional case, and combining Eqs (24) and (25), we then evaluate

$$\stackrel{}{r}_{12}\stackrel{}{r}_{23}_n=\frac{1}{2}r_{12}^2+r_{23}^2r_{31}^2_n=\frac{1}{2}r^2_n=\frac{1}{2}_0^{2R}P_n(r)r^2𝑑r=\frac{n}{n+2}R^2,$$
(26)
where we have utilized the fact that $`\stackrel{}{r}_{12}`$, $`\stackrel{}{r}_{23}`$, and $`\stackrel{}{r}_{31}`$ are three independent random vectors. The functions $`P_n(r)`$ in Eq. (26), which can be found in Refs. Tu (2001); Tu and Fischbach (2002); Kendall and Moran (1963); Santaló (1976); Solomon (1978); Ambartzumian (1990); Klain and Rota (1997); Fischbach (1996), are the probability density functions for the random distance $`r`$ between two random points in an $`n`$-dimensional spherical ball of radius $`R`$ and uniform density.

We consider next the analytical result in Eq. (20) for a circle ($`n=2`$) of radius $`R`$ and uniform density. A similar derivation can lead to Eqs. (21), (22), and (23), as well as to the case of $`n`$ dimensions where $`n3`$. We begin by expressing $`4`$ random points $`\stackrel{}{r}_1`$, $`\stackrel{}{r}_2`$, $`\stackrel{}{r}_3`$, and $`\stackrel{}{r}_4`$ in Cartesian coordinates, where $`\stackrel{}{r}_i=(x_i,y_i)`$. The expression in Eq. (20) can then be evaluated by writing

$$\left(\stackrel{}{r}_{12}\stackrel{}{r}_{23}\right)\left(\stackrel{}{r}_{34}\stackrel{}{r}_{41}\right)_2=\frac{_R^R𝑑x_1_{\sqrt{R^2x_1^2}}^{\sqrt{R^2x_1^2}}𝑑y_1\mathrm{}_R^R𝑑x_4_{\sqrt{R^2x_4^2}}^{\sqrt{R^2x_4^2}}f_1\times f_2𝑑y_4}{_R^R𝑑x_1_{\sqrt{R^2x_1^2}}^{\sqrt{R^2x_1^2}}𝑑y_1\mathrm{}_R^R𝑑x_4_{\sqrt{R^2x_4^2}}^{\sqrt{R^2x_4^2}}𝑑y_4}=\frac{3}{8}R^4,$$
(27)
where

$`f_1`$ $`=`$ $`\left(x_2x_1\right)\left(x_3x_2\right)+\left(y_2y_1\right)\left(y_3y_2\right),`$
$`f_2`$ $`=`$ $`\left(x_4x_3\right)\left(x_1x_4\right)+\left(y_4y_3\right)\left(y_1y_4\right).`$
A derivation of the general result using the probability density functions $`P_n(r)`$ in Eq. (26) can be found in Ref. (Tu and Fischbach ).

###### Acknowledgements.

The authors wish to thank T. K. Kuo and Dave Seaman for helpful discussions and the Purdue University Computing Center for computing support. This work was supported in part by the US Department of Energy contract DE-AC02-76ER1428.


text:
# 26 September, 2002 HIP-2002-43/TH Structure Functions are not Parton ProbabilitiesTalk at ICHEP 2002, Amsterdam (July 2002).

## 1 DIS DYNAMICS

Following the pioneering work of Drell, Levy and Yan  in 1969, the deep inelastic scattering ($`\mathrm{}N\mathrm{}+X`$, DIS) cross section has been commonly thought to be given by the probability of finding a parton in the target $`N`$ wave function, i.e., by the parton density. In the Bjorken limit of photon energy $`\nu `$ and virtuality $`Q^2`$ tending to infinity at fixed $`x_B=Q^2/2m_N\nu `$, final state interactions do not affect the DIS cross section – in a scalar theory such as considered in Ref. . In QCD, gluon exchange modifies the relation between the cross section and the parton density , and causes observed effects such as diffraction, shadowing and transverse spin asymmetry. In this talk I give a brief account of this.

The standard DIS ‘handbag’ shown in Fig. 1a may be viewed as a Light-Cone (LC) time-ordered diagram where the ordering of the vertices from left to right corresponds to an ordering in LC time $`y^+=y^0+y^3`$. In a frame where the virtual photon moves along the negative $`z`$-axis $`(q^{}2\nu ,q^+m_Nx_B)`$ there are no $`\gamma ^{}q\overline{q}`$ splittings of the kind shown in Fig. 1b, since the quark lines must have $`p^+>0`$. Due to the large value of $`q^{}`$ the scattering takes place almost instantaneously in LC time, $`y^+1/\nu `$.

The dashed lines in Fig. 1a indicate that the struck quark can (and typically does) rescatter elastically on the color field of the target. Since the exchanged gluon has spin 1 this soft (instantaneous, coulombic) scattering does not decrease with the energy $`\nu `$ of the quark – whereas scalar exchange would bring a factor $`1/\nu `$. Thus rescattering affects the DIS cross section at leading twist if it occurs within the coherence length $`\nu /Q^2=1/2m_Nx_B`$ of the virtual photon. Even at moderately small values of $`x_B`$ the coherence length is comparable to the size (several fm) of a nuclear target. Note that the rescattering occurs in vanishing LC time $`y^+`$ due to the relativistic motion of the struck quark. The DIS amplitudes nevertheless acquire complex dynamical phases from contributions of on-shell intermediate states between the rescatterings. The interference between amplitudes with and without rescattering gives rise to diffraction and shadowing .

## 2 SUBTLETIES OF LC GAUGE

QCD factorization in DIS determines the expression for the parton distribution in a general gauge to be

$`f_{\mathrm{q}/N}(x_B,Q^2)={\displaystyle \frac{1}{8\pi }}{\displaystyle 𝑑y^{}\mathrm{exp}(ix_Bp^+y^{}/2)}`$
$`\times N(p)|\overline{\mathrm{q}}(y^{})\gamma ^+[y^{};0]\mathrm{q}(0)|N(p)`$ (1)
where all fields are evaluated at $`y^+=y_{}=0`$. The path-ordered exponential

$$[y^{};0]\mathrm{P}\mathrm{exp}\left[ig_0^y^{}𝑑w^{}A^+(w^{})\right]$$
(2)
ensures the gauge invariance of the matrix element and contains the rescattering effects discussed above.

The exponential (2) reduces to unity in LC gauge ($`A^+=0`$). The matrix element (2) can then be expressed as an overlap of the target Fock state amplitudes $`\psi _n`$ at $`y^+=0`$ ,

$`f_{\mathrm{q}/N}\stackrel{A^+0}{}𝒫_{\mathrm{q}/N}={\displaystyle \underset{n}{}}{\displaystyle \underset{i}{}dx_id^2k_i}`$
$`\times |\psi _n(x_i,k_i)|^2{\displaystyle \underset{j=q}{}}\delta (x_Bx_j)`$ (3)
This appears like the probability of finding a quark in the initial state and is the origin of the misconception that the DIS cross section is an incoherent sum over the target quark density. However, the rescattering and interference effects discussed above have not vanished, they are just hidden by the peculiarities of LC gauge.

In Ref.  the properties of a perturbative, abelian model of DIS was studied in Feynman and LC gauge. In the $`x_B0`$ limit the dominant amplitudes were found to have on-shell intermediate states between the Coulomb gluon exchanges (dashed lines in Fig. 1). Such states can kinematically occur only after the $`\gamma ^{}`$ has been absorbed, and thus belong to the final state. Each on-shell rescattering amplitude is gauge invariant - hence the on-shell contributions cannot be eliminated by a gauge choice.

The model furthermore showed that in Feynman gauge the DIS cross section is affected only by diagrams where the struck quark rescatters. Indeed, final state interactions (FSI) between spectators cannot affect the DIS cross section at leading twist in this gauge (see, e.g.,  for a proof).

In LC gauge the situation was found to be reversed. Rescattering of the struck quark does not affect the DIS cross section, which is consistent with the fact that the path-ordered exponential $`(\text{2})`$ reduces to unity for $`A^+=0`$. Hence the reinteractions producing on-shell intermediate states in LC gauge occur between spectators. The singularity of the gluon propagator $`𝒟_{LC}^{\mu \nu }(k)`$ at $`k^+=0`$ in LC gauge enhances the $`k^+0`$ exchanges between spectators so that they modify the DIS cross section at leading twist .

The spectator interactions occur within a vanishing LC time $`y^+`$ of $`𝒪\left(1/\nu \right)`$ after the $`\gamma ^{}`$ interacts, even though they are spread over a longitudinal distance of $`𝒪`$(fm) in the target. Hence they can formally be included in the target Fock amplitudes $`\psi _n(y^+=0)`$. The $`k^+=0`$ gluon exchanges are then ‘zero modes’ which propagate along the light-like $`(y^+=0)`$ quantization surface of the Fock states. The Fock state wave functions thus defined are complex, with the rescattering and interference effects contained within their norm. However, the inclusion of the rescattering effects as zero modes of the target LC wave function appears artificial from a physical point of view.

## 3 REMARKS

The insight that DIS dynamics involves rescattering and interference effects has led to interesting further developments. Brodsky, Hwang and Schmidt  showed that the complex phase arising from rescattering causes a transverse spin asymmetry at leading twist in semi-inclusive DIS. Thus the azimuthal distribution of the current jet around the virtual photon momentum depends on the transverse polarization of the target nucleon in $`\gamma ^{}+N_{}jet+X`$. This “Sivers effect”  had previously been thought to be of higher twist . The effect is unrelated to the transverse spin carried by the target quark and thus complicates the measurement of the quark transversity distribution.

In semi-inclusive DIS the quark fields of the parton distribution (2) are at non-zero transverse separation, $`y_{}0`$. The integral in the path-ordered exponential (2) must then cover also this transverse distance (at asymptotic LC time) . The transverse path involves $`A^{}`$ and does not vanish in LC gauge – in fact it is the origin of the transverse spin asymmetry.

Diffractive DIS is described by two-gluon (Pomeron) exchange amplitudes like the one in Fig. 1b. The scattering is most naturally viewed in a frame where the virtual photon moves along the positive $`z`$-axis, i.e., with $`q^+2\nu `$. In this frame there is a $`\gamma ^{}q\overline{q}`$ transition which at low $`x_B`$ typically occurs long before the target. Crossing symmetry requires the amplitude to be dominantly imaginary at high energy.

The frame of Fig. 1b is related to that of Fig. 1a (where $`q^{}2\nu `$) by a rotation of $`180^{}`$ around the $`y`$-axis. Due to the nonperturbative nature of the target this rotation cannot in general be performed. The amplitude of Fig. 1b is often interpreted as corresponding, in the frame of Fig. 1a, to scattering off a ‘Pomeron component’ in the target wave function. The dynamics in the two frames can, however, be unambiguously related within the model of Ref. . There are no miracles – the amplitude remains imaginary after the rotation, and this phase arises from a rescattering in the final state. Thus views of the Pomeron as residing in the nucleon wave function before the $`\gamma ^{}`$ interacts are incorrect.

The above developments motivate comparisons of DIS with other hard processes and renewed checks of the universality of parton distributions. It has become clear that the transverse spin asymmetry has opposite sign in DIS and the Drell-Yan process . Furthermore, the dynamics of the rescattering which gives rise to shadowing is quite different in these two processes .

Acknowledgments. The work I have described was done in collaboration with S. J. Brodsky, N. Marchal, S. Peigné and F. Sannino. I am also grateful for helpful discussions with M. Diehl and X. Ji.