text:
# 1 Introduction

## 1 Introduction

The detection of Dark Matter and the studies of their properties  are of fundamental importance in particle physics and cosmology. The Weakly Interacting Massive Particles (WIMPs) are good candidates for “Cold” Dark Matter, and their experimental searches have gathered a lot of interests in recent years. The most promising avenue is to detect the nuclear recoil signatures due to elastic scatterings of WIMPs on the target isotopes. The typical energy depositions are only of the order of 10 keV, imposing big experimental challenges in terms of the detection of weak signals as well as background control at low energy close to detection threshold. A wide spectrum of experimental techniques is being pursued .

The DAMA experiment observed an annual modulation of nuclear recoil events  with NaI(Tl) scintillating crystal detectors, which can be interpreted as positive evidence of WIMPs due to the difference of the relative velocities of the Earth from the Halo sea within the year. This, however, contradicts the limit from the CDMS experiment  based on cryogenic technique.

There is still much room for new detector concept to thoroughly test the DAMA parameter space, and to push the sensitivities further. It would be of great interest if the sensitivities of WIMP searches can probe the level predicted by the various Super-Symmetry models.

There are potential merits of using CsI(Tl) scintillating crystals  for WIMP search and other low-energy low-background experiments . An experiment towards 200 kg of CsI(Tl) crystal scintillators to study low energy neutrino interactions at the Kuo-Sheng power reactor is being pursued , while the adaptation of the crystal for Dark Matter searches are the focus of several on-going projects .

The high-A content of the CsI enhances the sensitivities for the spin-independent interactions (which depends on $`\mathrm{A}^2`$) between the WIMPs and the target, relative to most other candidate target isotopes. The high-Z composition allows a compact design and provides large suppression of background due to ambient radioactivity if a three dimensional fiducial volume definition can be realized. Both 133Cs and 127I are 100% in their respective isotopic abundance. Being close in their mass numbers, the response to nuclear recoil from the interactions with WIMPs would be similar, allowing simpler interpretation of experimental signatures.

As a detector, the crystal has large light yield, low energy threshold and with pulse shape discrimination characteristics for $`\gamma `$/$`\alpha `$ separation. Its characteristics of being only slightly hygroscopic implies that it can be machined easily and does not require hermetic seal (that is, passive materials) in a large detector system. In addition, large (40 tons) electromagnetic calorimeter systems  have been constructed and made operational in high energy physics experiments, making this technology affordable and realistic to scale up. Considering all the associated costs, the price of CsI(Tl) is in fact less than that for NaI(Tl). In order to produce positive and definite evidence of the WIMPs, an accurate measurement of the annual modulation (where the maximal effects are only 7%) would be necessary such that the availability of large target mass is a very desirable feature.

One of the key issues to realize a Dark Matter search experiment with CsI(Tl) crystal scintillator is the studies of the experimental signatures of nuclear recoils due to WIMP-nuclei elastic scatterings. Nuclear recoils produce large charge density (dE/dx) such that the scintillating light yield is “quenched” and the timing profile of pulse shape is different relative to the same energy deposition by minimum ionizing particles .

These signatures are the same as the nuclear recoil events produced by elastic scattering of neutrons on nuclei. A measurement was performed with the neutron facility at the 13 MV Tandem accelerator at the China Institute of Atomic Energy at Beijing. The results of the quenching factors measurement and of pulse shape discrimination studies are reported in this article. They extend the work from the other recent measurements on CsI(Tl) , as well as on other scintillating crystals such as NaI(Tl) and CaF2(Eu) .

## 2 Experimental Setup and Procedures

The experiment was performed at HI-13 tandem accelerator at the China Institute of Atomic Energy (CIAE) in Beijing. A schematic diagram is shown in Figure 1. A pulsed deuteron beam at 5.6 MeV interacted with a deuterium gas target in a cell at a pressure of 6 bar and the dimensions of 1 cm in diameter and 3 cm in length. Neutrons at 8 MeV kinetic energy were produced at zero degree. The CsI(Tl) sample was located 2.02 meters away from the deuterium target. Neutrons at zero degree were selected by a 32 mm by 35 mm collimator of length 1.2 m. The collimator were surrounded by shielding materials like iron, paraffin, lead and polyethylene to reduce background. During this data taking, the repetition rate of the pulsed beam was set at 4 MHz and the average beam current was at 1 $`\mu `$A.

The scattering target, which also functioned as a detector, was a CsI(Tl) crystal scintillatorProducer: Unique Crystal, Beijing with diameter 3 cm and length 3 cm wrapped with teflon sheets, aluminium foil and black plastic tape. To minimize supporting structures, the detector was hung at the correct position by a piece of string. The readout was achieved by a 29 mm diameter photo-multiplier tube (PMT)CR110, Hamamatsu Photonics, China The PMT signals passed through an amplifier and shaper and were digitized by a 20 MHz, 8-bit, Flash Analog-to-Digital Convertor (FADC)  developed for the Kuo-Sheng reactor neutrino experiment . The trigger was provided from a time of flight (TOF) system to tag the elastically scattered neutrons. The $`\gamma `$-response of the CsI(Tl) detector was calibrated with an LED pulser operated at single photo-electron intensity, as well as standard 109Cd and 133Ba $`\gamma `$-sources. The calibration is 4 photo-electrons per keV of electron-equivalence energy.

Neutron tag was provided by liquid scintillator§§§Co-261, ST-451 detectors 105 mm in diameter and 50 mm in length equipped with PMTPhilips XP-2041 readout. The liquid scintillator provides pulse shape discrimination (PSD) capabilities for n-$`\gamma `$ differentiation. The PSD was achieved by commercial electronicsCANBERRA 2160A. The neutron tag signal was used to define the “START” timing for the TOF system. The pulsed deuteron “pick-off” signal provided by the accelerator was delayed and used as “STOP” for the TOF. The neutron detectors were placed at a distance of 0.5 to 2 m from the CsI(Tl) target at different angles, to provide optimized TOF resolution to differentiate events due to neutron elastic scattering off hydrogen atoms or heavier nuclei in the wrapping materials.

The FADC recorded data continuously in a circular buffer of 4k size. A START-STOP sequence from the TOF system initiated a trigger which stopped the FADC digitization after 25 $`\mu `$s. The pulse shape from the CsI(Tl) crystal was recorded for a pre-trigger and post-trigger periods of 5 $`\mu `$s and 25 $`\mu `$s, respectively. With the timing settings at data taking, the CsI(Tl) pulse started at $``$1.25 $`\mu `$s if it was due to nuclear recoils of heavy nuclei. This timing information provides a powerful means of the rejection of background due to proton recoils and pile-up events.

Data taking at each scattering angle was complemented by a measurement of the TOF background spectrum with empty target which was made of the CsI(Tl) wrapping materials only. Displayed in Figure 2 is the TOF spectra at 50 with empty and complete target. The lower and upper peaks are due to neutron elastic scattering off protons and the heavier nuclei such as carbon, aluminium, cesium, and iodine, respectively. The flat background is due to random coincidence. It can be seen that with the CsI(Tl) target in place, the upper TOF peak for neutron scattering off heavier nuclei becomes enhanced. The difference between the two TOF spectra is the nuclear recoils spectra from Cs and I. The proton-recoil events tagged by the low TOF peak provides a means to evaluate the signal yield and to confirm the distance between neutron detector and CsI target.

The reconstructed energy is defined as the pedestal-subtracted integrated area of the FADC pulse over 2.5 $`\mu `$s starting from the t0 of FADC (1.25 $`\mu `$s before TOF trigger in this data taking). This restricted time window helps to minimize the pile up effect from accidental low energy $`\gamma `$’s which are abundant along the beam line. Events are categorized into two groups by the central value between the two TOF peaks and their corresponding nuclear recoil energy spectra are compared. The energy spectra at 50 are shown in Figure 3a and 3b for events with high TOF values and low TOF values, respectively.

The energy spectrum tagged by the low TOF trigger is displayed in Figure 3b. Since the trigger is due to proton recoils from the wrapping materials, there should be no correlated signals from the CsI(Tl) target. The energy spectrum which shows a distinguished peak with long tail on the high-energy side is due to events from accidental $`\gamma `$’s arriving within the integration time-window.

The detector response for CsI nuclear recoil events was obtained by subtracting Figure 3b from Figure 3a with proper normalizations, where the factors are derived by normalizing the low TOF peaks from proton recoils between the full and empty target spectra.

## 3 Measurement Results

The resulting CsI recoil spectra obtained by the “statistical background subtraction” scheme discussed in Section 2 are shown in Figure 4 for eight different scattering angles. It can be seen that this procedure is valid since the background levels away from the peaks after subtraction are consistent with zero.

The nuclear recoil energy (T) is related to the neutron scattering angle ($`\theta `$) by a simple kinematical formula:

$$\mathrm{T}=\frac{2\mathrm{A}}{(1+\mathrm{A})^2}(1\mathrm{cos}\theta )\mathrm{E}_\mathrm{n},$$
(1)
where $`\mathrm{A}`$ is the target atomic mass and $`\mathrm{E}_\mathrm{n}`$ is the incident neutron energy. The neutron scattering angle can be considered to be the same in this case in both the laboratory and center-of-mass systems. Comparisons of the nuclear recoil spectra with those due to calibration sources allow the evaluation of the quenching factors, defined as the light yield from nuclear recoils versus that from $`\gamma `$ sources.

To enhance the signal-to-noise ratio in the large background environment, the nuclear recoil pulses are integrated for a restricted 2.5 $`\mu `$s (or 50 FADC time bins), as compared to the full time window 12.5 $`\mu `$s (or 250 FADC time bins) for 100% light collection in the $`\gamma `$-sources calibration measurements. A correction factor of 1.15 is applied to the nuclear recoils light yield to account for this partial summation. This factor was derived by averaging and comparing a large number of events due to nuclear recoils and 109Cd $`\gamma `$-source, as depicted in Figure 5. The estimated uncertainty is 10%, based on studies with different integration time window from 1.5 $`\mu `$s to 12.5 $`\mu `$s.

The deviation from the central value of nuclear recoil energy due to finite detector size is at most 2% for different scattering angles. This is checked by calculating the average recoil energy with proper weights of event angular distribution and detector acceptance. A shift of detector at 1 mm scale has negligible effect on the central value of nuclear recoil energy. The beam current was stable at the 1% level during data taking, as indicated by the stable data taking rate (typical value is about 4 Hz) as well as readings from beam monitor counters. The main systematic uncertainties originate from the analysis procedures, and amount to 16% at the lowest energy point at 20. The statistical and systematic uncertainties are then combined in quadrature to provide the total uncertainties for each data point.

The results at different scattering angles are summarized in Table 1. The measured quenching factors are depicted in Figure 6a, while the variation of electron-equivalence light yield (the “visible” energy) versus recoil energies are shown in Figure 6b. The other measurements  are overlaid for comparison. There is a clear trend of less quenching towards low energy points. This measurement achieved a lower threshold and with improved uncertainty compared to Ref. . The quenching factors are of the similar range to those for NaI(Tl) , which are typically 25% for Na and 8% for I. It can be seen that the light yield is higher for interactions with the high-A isotopes in CsI(Tl), providing big advantages to probing the spin-independent interactions. In addition, Figure 6b depicts a linear relationship between the electron-equivalence light yield and the nuclear recoil energy, as indicated by a linear fit our data set. Data from Ref.  shows more quenching and deviates from the linear regime at high energy. This is different from the highly non-linear response due to threshold effects observed in liquid scintillator  at the same range of recoil energy.

The total number of observed signal events are given by the area under the Gaussian fits in Figure 4. To derive the differential cross sections, the data taking time, the acceptance (proportional to $`r^2`$) and the efficiency correction factors have to be evaluated. The efficiency factors are due to the TOF selection criteria and background subtraction scheme. The efficiency is close to 100% for the well-separated TOF peaks, but is only 42% for the worst cases where the TOF peaks overlap. The contributions from inelastic scattering leading to the low excited levels for 133Cs and 127I are expected to be less than 0.1% relative to the elastic processes.

After all the correction factors are taken into account, the measured angular distribution for neutron elastic scatterings off Cs and I is displayed in Figure 7. The shape of the solid curve is the evaluated neutron elastic scattering cross sections on 133Cs and 127I from ENDF/B-VI library  which comes from the Optical Model calculation. The absolute normalization is from a best-fit value. There is excellent agreement between the relative recoil angular distribution from this measurement and the predictions from the Optical Model. In particular, this is a kinematical regime where the size of the nuclei is comparable to the wavelength of the incident neutrons, and the diffraction pattern is well-reproduced by the data.

The measurements demonstrate the validity of the nuclear recoil event selection procedures at all angles. In addition, the results in Figure 7 confirms the Optical Model predictions on neutron elastic scatterings from a direct measurement of nuclear recoils on heavy nuclei (instead of from the scattered final-state neutrons in transmission-type experiments) $``$ that is, they show that neutrons do elastically scatter off from the nuclei but not via some anomalous processes.

## 4 Pulse Shape Discrimination

The nuclear recoil measurements provide in addition an excellent data set for the study of PSD at small photo-electron level near detection threshold. Depicted in Figure 5 are the measured nuclear recoils at 50 and 95, as well as the $`\gamma `$-pulse due to the 22.1 keV line from 109Cd, averaged over many events and normalized to the same light collecition. It can be seen that heavily-ionizing events such as those due to $`\alpha `$-particles and nuclear recoils have faster decays than those from $`\gamma `$’s $``$ opposite to the response in liquid scintillator .

The neutron induced nuclear recoil pulse shapes are similar between 50 and 95 scattering angles, corresponding to 4.9 keV and 13.4 keV electron-equivalence light yield, respectively. The shape is also similar to that of the $`\alpha `$ source data  taken at high energy. At smaller scattering angles, the accidental low-energy $`\gamma `$-background becomes important such that the intrinsic nuclear recoil pulse shape cannot be derived. In the PSD studies, the recoil pulse shapes for these data set were taken to be the same as those measured at 50. The similarity of the pulse shape between 50 and 50 justifies this assumption.

Two independent approaches to obtain PSD were studied: (1) neural network technique  and (2) the “classical” method based on mean-time plus partial charge. For the neural network approach, we used 4000 events each from the above two data sets as training samples and additional 4000 events each as testing samples. The trained neural net has 20 input nodes (that is, 5 FADC time bins combined to form one node in Figure 5) in 5 $`\mu `$s window with 10 hidden nodes. The input node value should be greater than zero after pedestal correction. A negative value is reset to zero and the sum of the 20 input values are normalized to one. The algorithm used is back propagation. The excellent separation power from neural net study is depicted in Figure 8. The trained neural net is applied to the complete 50 data set. Good separation is obtained between the 43 keV nuclear recoils and the low energy $`\gamma `$ pulses. A cut on neural output greater than 0.9 selects almost all the CsI nuclear recoil events. The number of events as well as the mean and RMS electron-equivalence energy are consistent with those determined from the statistical background subtraction method in Section 3. The PSD separation with neural network method does not apply to the 20o data set. The neutron-beam environment gave rise to large accidental background (mostly low energy $`\gamma `$-events) which produce much stronger bias for the low energy nuclear recoil signatures. This would not be a problem in a realistic Dark Matter experiment where the background rate would be extremely low.

A large simulated sample was generated at different photon-statistics situation with these recoil/$`\gamma `$ pulse shapes, and the results are displayed in Figure 9. It can be seen that a 95% $`\gamma `$ rejection with $`>`$90% efficiency for pulses with an average 20 photo-electrons can be achieved. It is projected that the performance can be further improved in a low-background quiet environment which allows a longer sampling window. These are better than the PSD capabilities for NaI(Tl) detectors at the same light yield .

The alternative and complementary approach is one based on mean-time and partial charge. Two variables, the mean time $`\mathrm{t}`$ and the partial charge ratio R, are defined for each pulse, such that

$$\mathrm{t}=\frac{\underset{\mathrm{i}=1}{\overset{50}{}}(\mathrm{A}_\mathrm{i}\mathrm{t}_\mathrm{i})}{\underset{\mathrm{i}=1}{\overset{50}{}}\mathrm{A}_\mathrm{i}}$$
(2)
and

$$\mathrm{R}=\frac{\underset{\mathrm{i}=1}{\overset{30}{}}\mathrm{A}_\mathrm{i}}{\underset{\mathrm{i}=1}{\overset{50}{}}\mathrm{A}_\mathrm{i}},$$
(3)
where $`\mathrm{A}_\mathrm{i}`$ is the FADC amplitude at time bin $`\mathrm{t}_\mathrm{i}`$. The scattered plots for $`\mathrm{t}`$ versus TOF as well as $`\mathrm{t}`$ versus R for the data set at 95 are shown in Figures 10a&b, respectively. The electron-equivalence energy for CsI recoils is 13.4 keV at this scattering angle. There are clear separations between the Cs and I nuclear recoil events, represented by the box region, and the accidental $`\gamma `$-background for events triggered by the recoils of the wrapping materials (protons as well as other heavy nuclei). The cuts in $`\mathrm{t}`$ versus R space as defined in Figure 10b are subsequently applied to the other data set at lower recoil energies to select nuclear recoil events.

The events due to the 22 keV $`\gamma `$-rays from a 109Cd source are displayed in Figure 10c. The good overlap between the box and the events justifies that the cuts are valid to select “in-time” events. However, the separation between nuclear recoils from the in-time $`\gamma `$-background is not satisfactory. While the classical method is sufficient for event identification in this measurement where the main background are the accidental low-energy $`\gamma `$’s, it would be inadequate for a Dark Matter experiment where the $`\gamma `$-background also provide the trigger and are therefore in-time. In comparison with Figure 9, the neural network method gives 95% rejection of the in-time $`\gamma `$’s at a signal efficiency of better than 90% at $`>`$20 photo-electrons (that is, larger than 50 in this data set). The neural network method, therefore, will give superior performance in a Dark Matter experiment compared to the classical method.

The electron-equivalence light yield at different scattering angles for the nuclear recoil events selected by either of the two PSD algorithms are measured. The comparison between the two PSD algorithms with the subtraction method discussed in Section 3 is depicted in Figure 11, indicating that consistent results are obtained. The deviations of the mean values among these independent measurements are also consistent with the estimated systematic uncertainties.

Therefore, the robustness of the quenching factor results in Figures 6a&b measured by the statistical subtraction method discussed in Section 3 is further established. This method involves a straight-forward normalization and subtraction scheme, and is therefore expected to give more unbiased results compared to the two PSD algorithms.

## 5 Discussion and Summary

In this article, we presented a new measurement on the quenching factors of Cs and I in a CsI(Tl) crystal scintillator based on statistical background subtraction. Lower threshold and improved accuracies are achieved compared to previous measurements. The measured differential cross-section of neutron scattering from I and Cs is in excellent agreement with Optical Model derivations, and represents the first confirmation of the Optical Model on neutron elastic scattering cross-sections with a direct measurement of nuclear recoils from heavy nuclei. The recoil differential cross sections and the quenching factors are relevant to the studies of radiation damage in materials.

Two complementary pulse shaping discrimination techniques at the low photo-statistics regime were studied based on the nuclear recoil data. It is shown that the neural network approach gives superior rejection power to in-time $`\gamma `$-background compared to the classical method based on mean-time and partial charge. Both procedures are adequate for differentiating accidental $`\gamma `$’s. The electron-equivalence light yield derived from both PSD methods are consistent with those obtained by statistical subtraction. The internal consistencies of the light yield measurements among the three different algorithms and also with the Optical Model predictions of differential cross section enhance the reliability and robustness of the results.

By optimizing the detector geometry and using green-extended photo-cathode, a detector with several kg modular mass and a light yield of a few photo-electron per keV can be realized . The background from ambient radioactivity or intrinsic radio-purity is of course crucial to all low-background experiments. Their considerations for CsI(Tl) detectors are discussed in Refs. . Levels of better than the $`10^{12}`$ g/g level in concentration for the 238U and 232Th series have been demonstrated, assuming secular equilibrium. The potential problem of the internal 137Cs contaminations can be overcome via selection of clear ore materials and careful chemical processing and purification treatment .

Displayed in Figure 12 is the sensitivity plot for spin-independent WIMP interaction cross section per nucleon with a target mass of 1 ton of CsI(Tl). The sensitivities are based on the assumptions that the detection threshold is at 2 keV (or about 16 keV recoil energy) and still allows enough photo-electron statistics for background rejection with PSD, and that the observed background rate after applying the PSD cuts is at the same as or ten times better level than the NaI(Tl) experiments (1 and 0.1 per day per kg per keV, respectively). It can be seen that the DAMA allowed region, represented by the shaded region in Figure 12, can be thoroughly probed. Given that there is matured experience of scaling up the CsI(Tl) detector to multi-ton systems, there are potentials of further improvements in the sensitivities, and in the case of positive results, performing an accurate measurement of the annual modulation.

## 6 Acknowledgements

The authors would like to thank Drs. K.W. Ng, S.K. Kim and Y.D. Kim for fruitful discussions and helpful comments, and are grateful to the technical staff from CIAE and IHEP for operating the accelerator. This work was supported by contracts CosPa 89-N-FA01-1-4-2 from the Ministry of Education, Taiwan, NSC 89-2112-M-001-056 and NSC 90-2112-M-001-037 from the National Science Council, Taiwan, and NSF19975050 from the National Science Foundation, China. The support of H.B. Li is from contracts NSC 90-2112-M002-028 and MOE 89-N-FA01-1-0 under Prof. P.W.Y. Hwang.


format: text

text:
# Comment on “Stable Quantum Computation of Unstable Classical Chaos”

In a recent letter, Georgeot and Shepelyansky gs2  claim to have shown that there is a big advantage in using a quantum computer to simulate and study classical chaotic systems, over using a conventional (classical) computer. (Actually, the authors made this claim already in the second part of an earlier letter gs1 .)

Most of the paper gs2  is devoted to showing that the quantum simulation is much less affected by errors than the classical one. This is based on a clearly inappropriate comparison of the two cases. The authors consider the computationally simple, but chaotic, “Arnold cat map”, which is an area-preserving bijective map of a two-dimensional torus onto itself. They also consider a discretized version of this map which can be simulated exactly on a digital computer. The authors then propose to compute this discretized map on a quantum computer on many possible initial values simultaneously, utilizing “quantum parallelism”.

Strangely they then compare this quantum computation of the discretized map to a classical simulation of the full continuous map. Because they imagine the classical simulation to be done with fixed-precision arithmetic, the chaotic dynamics will quickly amplify rounding errors, so that after a few iterations of the map, the results will be totally wrong. It is then clear that the quantum simulation fares better, even when considering noise acting on the qubits, which the authors do. For this comparison, see e.g. figure 1 (left classical, right quantum) or the last paragraph of the paper.

Thus, in effect, the authors reach their conclusion in favor of the quantum simulation by demanding more from the classical than from the quantum simulation, specifically, that it should do an accurate simulation for any real initial values. Actually, the issue of rounding errors is not specific to classical or quantum computation, but is typical for digital computation. Note that the authors consider the usual kind of quantum computer consisting of qubits which is digital, namely, binary.

Apart from the discussion of errors, the authors make an unsubstantiated claim just before the last paragraph. They say that their quantum simulation lets one obtain “global quantities inaccessible by classical computation”. They propose to model the initial superposition of input values according to some (classical) probability density distribution (how, exactly, they do not say). They then want to extract information of interest about the final density, for example, by applying a quantum Fourier Transformation to the final superposition before observing the quantum computer.

Although this may allow the extraction of some information about the power spectrum of the final distribution, it is by no means clear that this could not be done just as efficiently on a classical computer. In particular, note that for the kind of quantum simulation described here, we need a reversible discrete map. This would, for example, allow one to evolve a few closely spaced final values backwards (on a classical computer), thereby permitting one to obtain information about the “fine structure” of a final density.

At any rate, it is an extraordinary claim, worthy of careful justification, to have found a new application for which a quantum computer provides exponential advantage over the known classical algorithms. This is especially true for “openly given” problems, thus without “black boxes” or limitations on communication etc. So far not much more than Shor’s algorithm(s) and, not surprisingly, the simulation of quantum systems achieved that.


text:
# A new method for Monte Carlo simulation of theories with Grassmann variables.

## Abstract

A new algorithm for simulation of theories with dynamical fermions is presented. The algorithm is based on obtaining the new configuration $`\stackrel{~}{U}`$ from the old one $`U`$ by solving the equation $`M(\stackrel{~}{U})\eta =\omega M(U)\eta `$, where $`M`$ is fermionic operator, $`\eta `$ is random Gaussian vector, and $`\omega `$ is random real number close to unity. This algorithm can be used for acceleration of current simulations in theories with fermions.

From the earliest days of Monte Carlo simulations fermionic fields have caused annoying difficulty, which stems from their being anticommuting variables. Although for the most actions in use one can eliminate fermions by an analytic integration, the resulting expressions involve determinants of large matrices, making the numerical simulations expensive.

Over the years many tricks have been developed to solve this problem. Hybrid Monte Carlo (HMC)  is now considered to be the standard simulation tool. In recent years it was challenged by the Multiboson (MB) method. The most popular version of MB, proposed by M. Lüscher , has been intensively studied by many authors (see  and references therein). The MB algorithm, proposed by A. Slavnov , is less known, and was tested in . Besides HMC and MB, let me mention some interesting ideas, based on the polymer  and Jordan-Wigner occupation number  representations, the direct evaluation of Grassmann integrals , the separation of low and high eigenmodes of Dirac operator , and the direct simulation of loop expansion by means of using the stochastic estimators .

Despite the variety of different approaches to the problem, there is a common impression that all existing algorithms remain inefficient . This is clearly seen if we compare the computational cost for the theories with dynamical fermions from the one side, and purely bosonic theories from the other side. Such a situation suggests that one should not stop the attempts to obtain relatively cheap fermion algorithm.

In this contribution I present a new computational strategy for treating dynamical fermions in Monte Carlo simulations. It has the virtues of exactness, nonlocality and finite step-size of update Unlike the other exact algorithms, the CPU-time required for the update in this algorithm grows only linearly with the volume $`V`$ of the system.

Suppose that we aim at sampling the partition function of theory with two flavors of degenerate fermion fields:

$$𝒵_{ferm}=Det\left(M^{}[U]M[U]\right)𝑑U$$
(1)
where $`M`$ is fermion operator, $`U`$ denotes the bosonic fields coupled to fermions. (For simplicity of formulae I do not include the purely bosonic contribution.) Then starting from old configuration $`U`$, one executes the following instructions:

* generate random vector $`\eta `$ with Gaussian distribution: $`P_G[\eta ]=Z_G^1e^{|\eta |^2}`$
* generate $`\omega [1ϵ;\frac{1}{1ϵ}]`$ with the probability: $`𝒫[\omega ,\eta ]min(1/\omega ^2,e^{\left(1\frac{1}{\omega ^2}\right)|\eta |^2})`$
* find the new configuration $`\stackrel{~}{U}`$ by solving the equation:

$$\stackrel{~}{M}\eta =\omega M\eta $$
(2)

where $`0ϵ<1`$ is an algorithmic parameter. Here and below I use the following short notations: $`\stackrel{~}{M}M[\stackrel{~}{U}];MM[U]`$.

In Ref.  it was proved that as far as the solution of eq.(2) is found in a reversible way 111i.e. the probability to obtain configuration $`\stackrel{~}{U}`$ starting from $`U`$ at any $`\omega `$ should be equal to probability to obtain $`U`$ starting from $`\stackrel{~}{U}`$ at $`1/\omega `$, the algorithm satisfies the detailed balance condition.

It is clear that choosing the procedure of solving the eq.(2) is crucial in this approach. Firstly, because obtaining the solution $`\stackrel{~}{U}`$ is the only computationally expensive ingredient of the algorithm. Secondly, because not every procedure of solving the eq.(2) can provide ergodicity. Finding a good procedure of solving the eq.(2) presents a serious challenge for the future investigations. The potential gain worth the efforts: e.g. finding some analytic solutions can give us a very cheap fermionic algorithm, comparable in cost with the algorithms for purely bosonic theories.

Making the first attempt to test the new algorithm, I solved the eq.(2) numerically by using the local iterative minimization of the quantity

$$R|(\stackrel{~}{M}\omega M)\eta |^2$$
(3)
for fixed $`U,\omega ,\eta `$. Starting from $`\stackrel{~}{U}=U`$, the minimization proceeds until $`R<\delta `$ is reached, where $`\delta `$ determines the accuracy of solution.

The tests were done for SU(3) QCD at $`8^4`$ lattice for the partition function (1) with

$$M[U]=1k^2D_{eo}D_{oe}$$
(4)
where $`D[U]`$ is Wilson difference operator. The even-odd preconditioning was used to reduce the computational cost. The hopping parameter was $`k=0.2`$, which gave the plaquette $`𝒫=0.0089(1)`$. The minimization of functional (3) was implemented by making the random moves in each color direction for all links lexicographically.

I checked the existence of solutions of eq.(2) at $`\omega [0.9,1.1]`$ for typical equilibrated configurations $`U`$. It was found that for considered $`\omega `$ one always finds some solution. On the average solving the eq.(2) at $`ϵ=0.05`$ with precision $`\delta =510^7`$ required 42 minimization iterations (it was checked, that improving the precision did not affect the results). Since the cost of one iteration is roughly equal to $`2N_{generator}=16`$ multiplications by matrix $`M`$, the total cost of finding the solution was approximately 670 matrix multiplication. This has to be compared with the cost for generating one trajectory with HMC. At step-size $`\tau =0.033`$ and trajectory length equal to $`1`$ (this ensures 70% acceptance), one trajectory cost approximately 4420 matrix multiplications, i.e. $`6.6`$ times more expensive than solving eq.(2).

My conclusion was that the usage of minimization of functional (3) for solving the eq.(2) is not favorable, because the algorithm became nonergodic. The problem is that probability $`𝒫[\omega ,\eta ]`$ to accept $`\omega <1`$ is strongly suppressed by large factor in the exponent, if $`\omega `$ is far from 1. Indeed, the squared norm of Gaussian vector can be estimated as $`|\eta |^2I_{dof}`$, where $`I_{dof}`$ is the dimensionality of vector space on which $`\eta `$ is defined. For even-odd preconditioned SU(3) QCD at $`8^4`$ lattice one has $`I_{dof}=24576`$! Therefore, using the new algorithm alone, one is restricted to choose between two unfavorable possibilities: 1) Using $`ϵ1/I_{dof}`$. It was observed that at such small values of $`ϵ`$ the evolution of $`U`$ fields becomes very slow, because the new configuration $`\stackrel{~}{U}`$ always lies too close to the old one; 2) Using $`ϵ1/I_{dof}`$. Then $`\omega <1`$ is almost never accepted, and the configurations are becoming smother.

In the second case the trouble appears due to the nonergodicity of algorithm at large $`ϵ`$. Nevertheless, even at large $`ϵ`$ our algorithm can be used for acceleration of other fermionic algorithms like HMC or MB, since the ergodicity is provided by them.

I tested the combination of the new algorithm and HMC. After each HMC trajectory the global move in configuration space was implemented by using our algorithm at $`ϵ=0.05`$. The efficiency of this algorithmic mixture was estimated by measuring the autocorrelation times for the plaquette. The results were: 1) Pure HMC, 1000 trajectories: $`\tau _{int}=1.4(2)`$; 2) HMC+new alg., 1000 (traj + $`\omega `$-update): $`\tau _{int}=0.7(1)`$. Note, that autocorellation time was reduced almost at no cost, because the global update of the new algorithm is much cheaper, than HMC trajectory. Of course, these tests are not very illustrative, since the autocorrelation times for theory with pure fermionic term at this kappa value are rather small. The further tests in models of practical relevance are needed.

Let me emphasize, that nonergodicity of the new algorithm was produced by the particular way of solving the eq.(2), which forces us to use large $`ϵ`$. One should think of finding some other procedure, which allows to use $`ϵ=0`$ (or $`ϵ0`$). Unfortunately, in this case the search for the solution of eq.(2) by minimizing the functional (3) does not work, because one immediately gets the trivial solution $`\stackrel{~}{U}=U`$. However, in typical case (e.g. SU(3) QCD), the solutions of eq.(2) are highly degenerate 222 In SU(3) QCD the expression (2) for even-odd preconditioned matrix provides $`12V`$ real equations for $`32V`$ variables., and the subspace of nontrivial solutions is not empty. Finding the reversible procedure of nontrivial solving the equation $`\stackrel{~}{M}\eta =M\eta `$, which samples the configuration space fast enough, is the perspective subject for future investigations.

I conclude that the efficiency of the new fermion algorithm drastically depends on the clever choice of procedure of solving the eq.(2). Finding the new efficient way of solving this equation can provide a very cheap simulation method.

Acknowledgments: I am very grateful to the Organizers of LATTICE2001 for generous help and support. I thank Ph. de Forcrand for helpful comments. This work was supported by grant for young scientists \# 01-01-06089, Russian Basic Research Fund \# 99-01-00190 and president grant for support of scientific schools # 0015-96046.


text:
# HST/NICMOS Imaging Survey of the Ophiuchus (Lynds 1688) Cluster

## 1 Introduction

Among molecular clouds within 200 pc of the Sun, the Ophiuchus cloud offers the best opportunity to study the formation of a stellar cluster at close range. At a distance of 145 pc (de Zeeuw et al., 1999), it contains a high density of low-mass young stellar objects (YSOs) deeply embedded in the cloud material, as evidenced by the extremely high gas and dust column densities measured there (Loren, Wootten & Wilking, 1990; Loren & Wootten, 1986; Wilking & Lada, 1983).

Most of the stellar content of the Ophiuchus cloud is obscured at visible wavelengths. However, early infrared observations revealed a rich embedded stellar population (Grasdalen, Strom, & Strom, 1973; Vrba et al., 1973; Fazio et al., 1976; Elias, 1978; Wilking & Lada, 1983). More recent observations using infrared arrays combined with molecular line maps have shown that these stars are concentrated in three dense molecular cores within the cloud (Greene & Young, 1992; Comeron et al., 1993; Strom, Kepner & Strom, 1995; Barsony et al., 1997). Age estimates based on infrared spectroscopic studies indicate that the stars are very young, with ages less than 106 yr (Greene & Meyer, 1995; Wilking, Greene, & Meyer, 1999; Luhman & Rieke, 1999).

Recent evidence suggests that most stars form in clusters (Carpenter, 2000), but there is relatively little theoretical guidance available to indicate how clusters develop (Myers, 2000; Adams & Myers, 2001; Meyer et al., 2000). The hierarchical structure of some groups has been said to arise from the hierarchical nature of condensations formed by turbulent processes (e.g. Klessen et al. 2000), but tests of cluster formation models require better knowledge of cluster spatial structure, and its temporal development. The youngest clusters are the best regions to study for this purpose because their stars have had relatively little time to move from their formation sites, and the spatial structure of the youngest members is essentially their spatial structure at birth.

In this contribution, we present the results of an infrared imaging survey using NICMOS on the Hubble Space Telescope (HST). The superior sensitivity and spatial resolution provided by HST were used to i) obtain a more complete census of the young stellar population in the dense star-forming regions of the cloud, ii) determine the degree of clustering and the binary and multiple star frequency, iii) identify candidate brown dwarfs, and iv) resolve the morphologies of deeply embedded extended sources.

Details of the observations and data analysis are described in § 2. Positions and photometry of all detected sources are presented in Table 2, along with cross-correlations with other designations from the literature. In § 3, we present results on the clustering and multiplicity of our sample, and a search for new brown dwarf candidates. A detailed description of extended sources is provided in § 4. A list of binary and multiple systems is presented in Table 3, and new candidate brown dwarfs are noted in Table 2. These results are summarized in § 5.

## 2 Observations and Data Analysis

### 2.1 Observations

Observations were obtained with HST/NICMOS Camera 3 during the June 1998 campaign. Camera 3 has a pixel scale of 0$`\stackrel{}{\mathrm{.}}`$20 pixel-1 and a field–of–view of $`51\mathrm{}\times 51\mathrm{}`$. Images were made through the F110W and F160W filters, with integration times of 39.95 seconds and 31.96 seconds, respectively. Thirteen target positions were chosen to correspond with both high surface density of previously known YSOs and high gas column density, and are listed in Table 1 (see also Figure 1). Each was observed in a $`3\times 3`$ spiral pattern with a 48″ dither, resulting in 13 $`2\stackrel{}{\mathrm{.}}45\times 2\stackrel{}{\mathrm{.}}45`$ mosaics. Each field was imaged twice in each filter, with a 3″ offset between image pairs to aid in identification and removal of bad pixels and cosmic rays.

### 2.2 Image Processing and Calibration

Uncalibrated images were processed using the IRAF/STSDAS package CALNICA, and the latest reference files provided by the Space Telescope Science Institute (STScI). Because Camera 3 suffers from a DC offset or bias problem (a.k.a. the “pedestal effect”), corrections were made using a program in IDL (PEDTHERM) provided by L. Bergeron of STScI. In regions dominated by bright extended emission, no pedestal corrections were made. For each mosaic, a sky frame was constructed of median filtered images, using only those fields which did not contain bright extended emission. Masks were constructed to prevent stars from contributing to both the pedestal and sky determinations. Sky-subtracted image pairs were cross-correlated to determine spatial shifts, and combined using the IRAF/STSDAS package DRIZZLE (Fruchter & Hook, 1997), which re-samples the image pairs onto a common pixel grid. Images containing no apparent sources were combined with DRIZZLE using the prescribed telescope dither.

Reduced images were flux calibrated using recently determined photometric scale factors, $`2.873\times 10^6`$ Jy $`(\mathrm{ADU}/\mathrm{sec})^1`$ at 1.1 $`\mu `$m and $`2.776\times 10^6`$ Jy $`(\mathrm{ADU}/\mathrm{sec})^1`$ at 1.6 $`\mu `$m (Rieke, 1999). Corresponding zero points were calculated on the Vega system, assuming zero magnitude flux densities at 1775 and 1083 Jy and effective wavelengths of 1.104 $`\mu `$m and 1.593 $`\mu `$m for F110W and F160W respectively (Rieke, 1999).

### 2.3 Photometry and Astrometry

Magnitudes were measured using aperture photometry routines in IRAF/DIGIPHOT, and are reported in Table 2. For point sources we used an aperture of radius 0$`\stackrel{}{\mathrm{.}}`$5 (5 pixels in our drizzled images). A larger aperture radius of 3″ was used to measure the flux from extended sources, as noted in Table 2. By adding artificial stars to our data, we estimate that our photometry is 90% complete to 21.0 mag at F160W and 21.5 mag at F110W. Within those limits, we detect 165 sources at F160W and 65 sources at F110W. A census of near-infrared surveys in the literature (Luhman & Rieke, 1999; Barsony et al., 1997; Strom, Kepner & Strom, 1995; Comeron et al., 1993; Greene & Young, 1992) indicates that approximately two thirds of the sources listed in Table 2 are previously unreported.

Enough overlap exists between our survey and previous studies that we can derive transformations between the F110W and F160W magnitudes and a standard ground-based system. Comparing our photometry with that of Barsony et al. (1997) for 21 sources, we find the following linear relations: $`\mathrm{m}(\mathrm{F110W})=(1.07\pm 0.001)\times \mathrm{J}(0.44\pm 0.012)`$, and $`\mathrm{m}(\mathrm{F160W})=(1.02\pm 0.001)\times \mathrm{H}+(0.03\pm 0.007)`$, where J and H are on the CIT system. These relations were used to convert our F110W and F160W magnitudes to J and H magnitudes, allowing their placement in a color-magnitude diagram ( see Figure 7, discussed in § 3.4).

Because no guide stars were visible in the highly extincted target fields, observations were unguided, leading to some drift in the field centers (of order a few to 20″ throughout a 3x3 map). This drift, coupled with the fact that many frames have no or very few sources, complicated the determination of accurate astrometric solutions with these data. Instead, the coordinate system of Barsony et al. (1997) was adopted, since that survey covers a large fraction of the cloud and has been extensively cross-correlated with previous studies. Barsony et al. estimated 1$`\sigma `$ uncertainties of 1$`\stackrel{}{\mathrm{.}}`$2 in their absolute positions.

## 3 Results

### 3.1 Cloud membership and background stars

Given the sensitivity of our survey, we can expect to detect some background stars through the molecular cloud, so some means of distinguishing cloud members from background sources would be helpful. Unfortunately, with photometry in only 2 bands we cannot use colors to distinguish reddened background field stars from non–nebulous pre–main sequence stars embedded within the cloud. However, we can estimate the background by considering the observed distribution of objects with respect to the column density of molecular gas, and the expected source counts in this part of the sky from a model of the Galaxy.

In Figure 2, we reproduce the extinction map made by Wilking & Lada (1983), based on observations of C18O (J=1–0) and 12CO (J=1–0). On it we plot the locations of all 165 sources detected at F160W. Our HST mosaics are coincident with high $`\mathrm{A}_\mathrm{V}`$ regions, with the exception of some fields on the eastern sides of cores A, B, and F. These fields coincide with a steep gradient in the cloud column density, and in them large numbers of stars are seen. As the extinction increases toward the center of the cloud, the number of stars decreases. Clearly, a significant fraction of detected sources on the edges of the cloud cores must be background stars.

The infrared Galactic model of Wainscoat et al. (1992) predicts the existence of approximately 140 stars arcmin-2 brighter than our detection limit in the direction of the Ophiuchus cloud, assuming $`\mathrm{A}_\mathrm{V}=0`$. However, when the model prediction is convolved with the $`\mathrm{A}_\mathrm{V}`$ map, the resulting background is substantially smaller. This was done by calculating the expected background for each pixel in the $`\mathrm{A}_\mathrm{V}`$ map, then summing over all pixels to obtain the total number of expected background stars as a function of H magnitude.

In Figure 3 we plot the observed H-band distribution of stars (using the conversion from magnitude at F160W to magnitude at H given in Section 2.3). The distribution predicted by Wainscoat’s model and reddened with our $`\mathrm{A}_\mathrm{V}`$ map (assuming $`\mathrm{A}_\mathrm{H}`$ = 0.155$`\mathrm{A}_\mathrm{V}`$, Cohen et al. (1982)) is shown as a heavy dotted line. It appears that the background is small, becoming dominant only for H$`>`$20 mag. For comparison, we also plot the expected background as seen through a uniform extinction cloud of $`\mathrm{A}_\mathrm{V}=`$ 35, 40, and 45 mag. The background modelled with our extinction map agrees closely with that for a uniform extinction of 40 mag. For this reason, we shall assume that all stars lying outside the $`\mathrm{A}_\mathrm{V}=`$40 mag. contour in Figure 2 are background. Within the 40 mag. contour, the probable number of background stars given by the Wainscoat model ranges from approximately 0.50 stars arcmin-2 (for $`\mathrm{A}_\mathrm{V}`$=40 mag) to 0.02 stars arcmin-2 (for $`\mathrm{A}_\mathrm{V}`$=80 mag). Of the 165 sources detected in our survey, 108 sources lie within the 40 mag contour, and 58 of these are new detections, too faint to have been detected in previous surveys. In the discussions of clustering, multiplicity, and brown dwarfs which follow, only these 108 sources will be considered.

### 3.2 The spatial distribution of young stars in the Ophiuchus cloud

Ground–based imaging surveys at 2 $`\mu `$m revealed that the young stars in the Ophiuchus cloud are clustered in 3–4 main groups associated with dense cores of molecular gas (Luhman & Rieke, 1999; Barsony et al., 1997; Strom, Kepner & Strom, 1995; Comeron et al., 1993; Greene & Young, 1992; Barsony et al., 1989; Rieke, Ashok, & Boyle, 1989). Our HST survey targeted these peaks in the YSO surface density distribution, allowing us to examine structure within these stellar concentrations to a greater depth and with higher resolution than before.

The highest concentration of sources is in core A, where 23 stars are detected at F160W within an area $`0.08\times 0.05`$ pc in size (1$`\stackrel{}{\mathrm{.}}`$9 x 1$`\stackrel{}{\mathrm{.}}`$2 , centered on HST position 2 in Table 1). The extinction in this area is $`\mathrm{A}_\mathrm{V}=4580`$ mag, and averages 60 mag, leading to an estimated background contribution of 0.3$`\pm `$2 stars. The YSO surface density in the area, excluding background objects, is then $`5\times 10^3`$ stars pc-2. Assuming a range of depths from half to twice the area width, we estimate stellar volume densities of $`410\times 10^4`$ stars pc-3, comparable to the density of $`2\times 10^4`$ stars pc-3 determined for the core of the Orion nebula cluster111It is worth noting that, while these space densities are similar, the ONC core occupies a volume about 10 times larger than the one we describe in Oph, and contains hundreds of stars. (Hillenbrand & Hartmann 1998). Other peaks in the stellar density are found in core B (HST positions 9 and 10) and core E (HST positions 12 and 13). Volume densities there are slightly less, ranging from $`0.65\times 10^4`$ stars pc-3.

Inverting the volume density yields a mean spacing between stars in the core A peak of $`0.020.03`$ pc, or $`40005000`$ AU, and $`560011000`$ AU in the B and E cores. This is similar to the 6000 AU “fragmentation” scale identified by Motte et al. (1998) in their analysis of 1.3 mm continuum dust clumps, which are concentrated in the same three regions as the infrared sources reported here.

In fact, the embedded stars in Ophiuchus and the dust clumps detected at 1.3 mm have similar, but not identical, spatial distributions. The clumps show the same sub-clustering shown by the stars, as seen in Figure 4, where the positions of the starless clumps from Motte et al. (1998; open squares) have been plotted along with the positions of the stars (filled circles). Interestingly, the clumps are more strongly associated with the highest extinction regions of the cloud, whereas the stars appear to cluster around the edges of these regions.

This similar spatial relation of the stars and the millimeter continuum clumps in each of the three subgroups of the Ophiuchus core suggests that we are observing a real effect with the same explanation (origin?) in each subgroup. Some of the millimeter continuum clumps could be highly extincted YSOs, but this seems unlikely because most of the clumps are extended, some have associated compact molecular line emission, and the necessary extinction, $`\mathrm{A}_\mathrm{V}>1000`$ mag, is in all other known cases accompanied by emission at far- and mid-infrared wavelengths (cf. Ladd et al. 1991). Adopting the interpretation of Motte et al. (1998) that the millimeter continuum clumps are prestellar, a more likely explanation is that in each subgroup, the region of highest extinction is a “starless core” which has not yet formed stars, but has formed numerous protostellar condensations with relative spacing (4000-5000 AU) similar to that of the surrounding young stars. In this picture, the similar spacing of prestellar and young stellar objects suggests that the young stellar objects have neither concentrated nor dispersed significantly since their formation. This evidence for spatial segregation by age may offer important clues to how stars form in clusters.

### 3.3 Multiplicity and clustering

To search for apparent binary pairs and higher-order multiple systems, we performed both a visual examination of all images and a nearest neighbor calculation for all detected sources. However in order to decrease the likelihood of mistaking background stars for companions, we restricted the sample to the 108 stars projected within the $`\mathrm{A}_\mathrm{V}`$40 mag contour of the cloud (as discussed in § 3.1).

For relatively faint stars (m(F160W)$``$18), our survey is sensitive to separations comparable to the resolution of the array (0$`\stackrel{}{\mathrm{.}}`$2 pixel-1), although due to the complex NICMOS PSF, our detection of faint companions within $`\theta <0.5`$″ of bright stars may be incomplete. To guard against bias, we took our estimated limiting magnitude for faint companions (m(F160W)$``$18) and imposed this limit on the sample, further restricting our sample to 62 stars. The maximum separation considered in our search was 10″(1450 AU), chosen to coincide with previous searches for pre-main sequence binaries. Within the separation range of 0$`\stackrel{}{\mathrm{.}}`$2 to 10″, we detected seven apparent binary pairs and six apparent triple systems. They are listed in Table 3, and shown in Figure 5. The triple systems pictured in panels i, j, and k in Figure 5, are unusual in that they appear to be non-hierarchical. The other triple systems reported here (panels l and m) are organized in the usual way, having a close pair with a widely separated third member.

The multiplicity fraction, defined as the ratio of the number of binary and multiple systems detected, to the number of single, binary and multiple systems observed, or mf=B+M/(S+B+M), is then 0.30$`\pm `$0.08. Several previous multiplicity studies have included targets in Lynds 1688, each one sensitive to a specific range of companion separations and limiting magnitudes. In the lunar occultation experiment of Simon et al. (1995), observations were sensitive to stars brighter than K$``$11 and binary separations of 0$`\stackrel{}{\mathrm{.}}`$005 to 10″. Observing 35 systems, they detected 10 binaries, two triples, and one quadruple, for a multiplicity fraction of 0.37$`\pm `$0.10.

There are important differences between Simon et al. (1995; hereafter S95) and this study which should be considered when comparing their results. First, the minimum separation detectable by S95 was 0$`\stackrel{}{\mathrm{.}}`$005, whereas ours is $``$0$`\stackrel{}{\mathrm{.}}`$2. The limiting magnitude of the samples also differ: $`\mathrm{K}_{\mathrm{lim}}12`$ and $`\mathrm{H}_{\mathrm{lim}}=18`$ in S95 and this study, respectively. For projected separations between 0$`\stackrel{}{\mathrm{.}}`$2 and 10″ only, S95 obtain mf$`0.26\pm 0.08`$. When we compare this with our mf of $`0.30\pm 0.08`$ for a deeper limiting magnitude, we conclude that the two results do not differ significantly. Within the uncertainties, the multiplicity fractions of S95 and this work are in agreement.

The relationship between the separation of binary pairs and clustering on a larger scale has been explored in a number of nearby star forming regions, including Ophiuchus (Ghez, Neugebauer, & Matthews, 1993; Reipurth & Zinnecker, 1993; Strom, Kepner & Strom, 1995; Nakajima et al., 1998). Of particular interest is the mean surface density of companions (MSDC), which Larson (1995) applied to young stars in the Taurus-Auriga star-forming region, finding clustering on two distinct scales. For large-scale clustering the MSDC was found to have a power-law slope of $``$ $`0.6`$, and for small separations a slope of $``$ $`2`$ was found. The “break point” between large and small scale slopes was determined to be $``$ $`0.04`$ pc. Simon (1997; hereafter S97) computed the MSDC for Ophiuchus, Taurus, and the Trapezium, finding similar power laws for each.

We computed the MSDC for the 108 stars projected within the $`\mathrm{A}_\mathrm{V}`$40 mag contour of Figure 2. The result is shown in Figure 6 for two samples of different depths: a shallow magnitude cutoff (H$``$14) and a deep limiting magnitude (H$``$21). As seen in Figure 6, the two samples have different break points. The shallow magnitude limit was selected to approximately match that of S97 (K$``$12), and not surprisingly, the two MSDCs have the same slopes and break points. The deeper sample shows a break point on a smaller spatial scale (by an order of magnitude). Thus, in a given star forming region, the properties of the MSDC can vary according to the depth of the sample. Bate et al. (1998) cautioned against overinterpreting the MSDC, in part for this reason. We offer this comparison as further caution.

### 3.4 New candidate brown dwarfs

With completeness limits of 21.5 and 21.0 magnitudes at F110W and F160W respectively, our survey is sensitive to very low-mass objects. For example, we should have detected stars of spectral type L4 V, through as much as 17 $`\mathrm{A}_\mathrm{V}`$, assuming $`\mathrm{M}_\mathrm{J}=13.2`$ and $`\mathrm{M}_\mathrm{H}=12.3`$ (Kirkpatrick et al., 1999). Such very low-mass objects are relatively “blue” in intrinsic color, typically having $`(\mathrm{J}\mathrm{H})<2`$ (Kirkpatrick et al. 1999), and so might be distinguishable from heavily reddened background stars in an infrared color-magnitude diagram.

Using the transformations from F110W and F160W to J and H listed in § 2.3 and the fluxes in Table 2, we converted our measured magnitudes to J and H magnitudes and plotted them in the color-magnitude diagram shown in Figure 7. Also plotted are the $`3\times 10^5`$ yr isochrone and the ZAMS from the pre-main sequence stellar evolution models of D’Antona & Mazzitelli (1997). Only the 65 sources detected in both the F110W and F160W bands are shown.

Brown dwarfs cannot be positively identified using near-IR photometry alone. However, those sources which are faint and relatively blue are more likely to be low-mass objects in the cloud than background stars seen through the cloud, especially if they are coincident with high extinction regions. We identified ten candidate brown dwarfs in our sample, using the following criteria: i) (J-H)$`<`$ 3.0, ii) J$`>`$15, and iii) position coincident with cloud extinction values of $`\mathrm{A}_\mathrm{V}>50`$ mag. Spectral types have been published for three of these objects. One object, 162622$``$242409 is almost certainly sub-stellar, having been assigned a spectral type of M6.5 by both Wilking, Greene, & Meyer (1999; hereafter WGM) and Luhman & Rieke (1999; hereafter LR), and having a luminosity of 0.002-0.003 $`\mathrm{L}_{}`$. The mass of object 162622$``$242354 is less certain. WGM classify it as M8.5 with $`\mathrm{L}_{\mathrm{bol}}0.004\mathrm{L}_{}`$, whereas LR adopt a spectral type of M6.5 and a luminosity of 0.067 $`\mathrm{L}_{}`$. The third classified object, 162659$``$243556, is a star of type M4 (LR), too hot to be a brown dwarf. The remaining seven (new) brown dwarf candidates are noted in Table 2. They are located in core A, which contains the highest density of stars in the survey (§ 3.1). Two of the candidates, 162622$``$242409 and 162622$``$242408, are members of triple systems (§ 3.3 and Table 3). One of the candidates, 16265$``$242303, is just on the red edge of our color criterion, but qualifies for selection when photometric uncertainties are taken into account.

There are a total of ten known brown dwarfs in the Ophiuchus cluster (LR, WGM) and another five sources known to occupy the transition region in the HR diagram between the stellar and substellar regimes (WGM). If all of the seven new candidates identified here are indeed substellar, the number of brown dwarfs known in the Ophiuchus cluster will have increased by 50%-70%. We can assess the implications for the stellar initial mass function (IMF) in the cluster. However, because fainter, lower-mass objects cannot be seen as deeply into the cloud as higher-mass stars, we restrict our analysis to an extinction-limited sample. In Figure 7, there are 30 sources within an extinction-limited region delimited by the $`3\times 10^5`$ yr isochrone and the $`\mathrm{A}_\mathrm{V}<20`$ mag limit. These include our seven new candidate brown dwarfs, and the two confirmed brown dwarfs. If all seven of the new candidates are brown dwarfs, the total fraction of substellar objects in the extinction-limited sample is 9/30 = 30%. This is consistent with the results of LR, who argued that the IMF includes a large fraction of sub-stellar mass objects.

## 4 Morphologies of extended sources

Several interesting extended objects were covered by our survey, including the well-known infrared sources GSS 30, YLW 15A and YLW 16A, as well as some catalogued objects whose morphologies were previously unresolved. These sources are shown in Figure 8 and described in detail below.

GSS 30 (162621$``$242306): This multiple source is perhaps the most studied object in Ophiuchus. The illuminating star GSS 30–IRS 1 (162621$``$242306) is a class I source of unknown spectral type. A K-band spectrum obtained by LR showed the source to be heavily veiled (r$`>`$2) with no photospheric absorption features, and with Br$`\gamma `$ in emission. Two other infrared sources, GSS 30–IRS 2 and GSS 30–IRS 3, are located within the projected bounds of the nebula.

The nebula illuminated by IRS 1 has itself been the subject of many investigations. Near-infrared polarimetry images (Chrysostomou et al., 1997, 1996; Weintraub et al., 1993; Tamura et al., 1991; Castelaz et al., 1985) delineate a bipolar morphology coincident with the near-IR emission on the northeast side, and extending to the southwest side, where there is less infrared emission. The polarization pattern is approximately centrosymmetric around the illuminating source, and shows a “polarization disk” orthogonal to the long axis of the nebula. Zhang et al. (1997) detected a somewhat flattened molecular core in C18O and 13CO (J=1–0), in roughly the same orientation as the polarization disk. In the high resolution image presented in Figure 8 one sees the characteristic hourglass morphology of a bipolar nebula. The emission is dominated by the northeast lobe, a fact which has led previous investigators to propose that the fainter southwest lobe is obscured by a tilted circumstellar disk or toroid (Castelaz et al., 1985; Tamura et al., 1991; Weintraub et al., 1993; Chrysostomou et al., 1997).

GSS 30–IRS 2 (162622$``$242254) is a class III source of late K/early M spectral type (Luhman & Rieke, 1999; Greene & Meyer, 1995), located 20″ to the northeast of GSS 30–IRS 1. Its relatively advanced evolutionary class indicates it is probably unrelated and slightly foreground to the nebula.

GSS 30–IRS 3 (162621$``$242251) is a class I source of unknown spectral type, 15″ to the northeast of GSS 30–IRS 1, and 11″ southwest of GSS 30–IRS 2. This source has an extended point spread function and a crescent-shaped morphology. In addition, the F160W image shows a patch of obscuration (about an arcsecond in size) adjacent to the emission peak. Thus it is likely that the emission we see is scattered light from an embedded source. GSS 30–IRS 3 is also a strong radio source (LFAM 1), with emission at 6 cm (Leous et al., 1991).

GY 30 (162625$``$242303): Fan-shaped nebulosity extends to the east of this apparently low-mass object.

(162704$``$243707): Detected only at F160W, this source has a bipolar morphology strongly suggestive of a YSO outflow cavity, similar to that seen in GY 244, but fainter. This source was reported in Brandner et al. (2000), who imaged it at 2.2 $`\mu `$m.

WL 15 (162709$``$243718): WL 15 (Elias 29) is a class I source (Wilking, Lada, & Young, 1989) with a molecular outflow (Bontemps, et al., 1996; Sekimoto et al., 1997). The near-IR spectrum of El 29 resembles that of GSS 30–IRS1; a heavily veiled, featureless spectrum with the Br $`\gamma `$ line in emission (Luhman & Rieke, 1999).

Simon et al. (1987) found that the best fit to their lunar occultation observations was a two-component model having a central 7 mas component producing most of the emission at 2.2 $`\mu `$m and a larger component 0$`\stackrel{}{\mathrm{.}}`$4 in extent, roughly centered on the smaller one. Our F160W image is consistent with this. The source is slightly elongated and has a measured FWHM of 0$`\stackrel{}{\mathrm{.}}`$33 (3.3 pixels in our dithered image).

GY 244 (162717$``$242856): This source has a class I spectral energy distribution (Wilking, Lada, & Young, 1989), and a K-band spectral type of M4 (Luhman & Rieke, 1999). Our images show a bipolar nebula associated with this source. The morphology of the nebula is characteristic of the scattered light seen from YSO outflow cavities (Whitney & Hartmann, 1992) and suggests that the SW lobe is inclined somewhat toward us. Faint nebulosity in the F160W image indicates that the stars near GY 244 (GY 246, GY 247, and GY 249) also have some circumstellar envelope material associated with them.

(162724$``$244102): Detected only at F160W, this source has a faint bipolar morphology similar to that of (162704$``$243707). Also reported in Brandner et al. (2000).

YLW 15A (162726$``$244051): This embedded source has a heavily veiled, featureless 2 $`\mu `$m spectrum (Luhman & Rieke, 1999). It is included in our list of binary systems (Table 3), on the basis of a companion located $``$6″ to the northwest.

YLW 16A (162728$``$243934): We detect two non-point sources at the position of YLW 16A, separated by 0$`\stackrel{}{\mathrm{.}}`$5 (position angle 270°). The flux ratio of the two peaks is 1.5 at 1.1 $`\mu `$m and 1.1 at 1.6 $`\mu `$m, with a large uncertainty due to the extended nature of the sources. As both of these sources are extended, it unclear whether they are actually two embedded stars, or a single embedded star seen in scattered light; possibly a star/disk/envelope system.

The appearance of the two peaks is more diffuse at 1.1 $`\mu `$m than at 1.6 $`\mu `$m, as would be expected if the light we detect is mainly scattered light. In addition, the flux ratio is greater at 1.1 $`\mu `$m than at 1.6 $`\mu `$m, suggesting differential reddening to the two peaks. Lunar occultation observations at 2 $`\mu `$m of YLW 16A (Simon et al., 1987) failed to resolve this source as a binary, instead showing it to be a single extended source of size (at K) $``$0$`\stackrel{}{\mathrm{.}}`$5, which is the angular separation of the two peaks observed in our HST images. Thus it is possible that YLW 16A is a single star+disk+envelope, rather than an embedded binary system.

GY 273 (162728$``$242721): Large sigma-shaped nebula. It is unclear from our images where the illuminating source is.

## 5 Summary

We have presented a NICMOS 3 survey in the F110W and F160W filters of the dense star-forming cores of the Ophiuchus molecular cloud. Among our results:

1. The number of sources detected at F160W is 165, of which 65 were detected at F110W. Based on a search of the literature, we estimate that approximately two-thirds of the sources detected at F160W were previously uncatalogued. Of the 108 sources located within the $`\mathrm{A}_\mathrm{V}=`$40 mag. contour of the cloud, 58 were previously uncatalogued.

2. The Ophiuchus star forming region has multiple peaks in stellar density. In core A we measure $`\mathrm{n}_{}=4\times 10^4\mathrm{stars}\mathrm{pc}^3`$ within a volume 0.05-0.08 pc on a side, and in cores B and E, $`\mathrm{n}_{}=0.65\times 10^3\mathrm{stars}\mathrm{pc}^3`$ within similar volumes. Such densities are similar to those seen on larger scales in rich clusters such as the Orion Nebula cluster.

3. Thirteen apparent multiple systems (eight binary pairs and five triples) with projected separations in the range 0$`\stackrel{}{\mathrm{.}}`$2 to 10″ (29 to 1450 AU) were detected.

4. Seven new candidate brown dwarfs were identified from their positions in a color-magnitude diagram. According to our analysis of an extinction-limited sample, sub-stellar mass objects may account for as many as 30% of the sources in the core of the Ophiuchus cluster.

5. The unprecedented combination of resolution and sensitivity provided by HST has revealed new structures in the infrared sources in Ophiuchus. Bipolar structure is clearly seen in five objects.

We are grateful to Al Schultz at STScI for valuable advice, and to L. Bergeron of STScI for expert help with NICMOS data reduction. We thank Bruce Wilking for providing data which allowed us to reproduce his $`\mathrm{A}_\mathrm{V}`$ map, and John Carpenter for sharing his encoded version of Wainscoat et al.’s Galaxy model.


text:
# DESY 01-144 October 2001 The spin structure of the Λ hyperon in quenched lattice QCDTalk given by M. Göckeler at Lat01, Berlin, Germany.

## 1 INTRODUCTION

In the quark model, the $`\mathrm{\Lambda }`$ hyperon is a $`uds`$ state, and its spin is completely carried by the $`s`$ quark. Combining, on the other hand, the experimental results for the $`g_1`$ structure function of the proton with the data on hyperon $`\beta `$ decay and assuming SU(3) flavour symmetry one finds that only $`60\%`$ of the $`\mathrm{\Lambda }`$ spin is carried by $`s`$ (and $`\overline{s}`$) quarks, while $`u`$ ($`\overline{u}`$) and $`d`$ ($`\overline{d}`$) quarks contribute $`40\%`$. The spin structure of the $`\mathrm{\Lambda }`$ is of special interest, because its polarisation can easily be measured via the self-analysing weak decay $`\mathrm{\Lambda }p\pi ^{}`$. Indeed, the $`\mathrm{\Lambda }`$ polarisation has been measured at the $`Z^0`$ pole in $`e^+e^{}`$ annihilation  where it is mainly due to the $`s`$ quarks, which when produced via $`Z^0`$ decays have an average polarisation of $`0.91`$. Furthermore, $`\mathrm{\Lambda }`$ polarisation has been studied in deep-inelastic scattering of polarised positrons on unpolarised protons in the current fragmentation region, i.e. selecting $`\mathrm{\Lambda }`$’s which most likely originate from the struck quark . While polarised $`s`$ quarks are rare in the proton and their contribution to electromagnetic processes is suppressed by a factor $`Q_s^2=1/9`$, there are plenty of polarised $`u`$ quarks contributing with $`Q_u^2=4/9`$. Thus even a small correlation between the spins of the $`u`$’s and the $`\mathrm{\Lambda }`$’s could make the $`\mathrm{\Lambda }`$’s useful as probes of the polarised $`u`$-quark distribution in the proton . Needless to say, however, that the interpretation of the experiments is complicated by fragmentation effects, $`\mathrm{\Lambda }`$’s from the decay of heavier hyperons, etc.

On the theoretical side, for a spin 1/2 baryon the fraction $`\mathrm{\Delta }q`$ of the spin carried by the quarks (and antiquarks) of flavour $`q`$ is given in terms of the forward matrix element of the axial vector current:

$$p,s|\overline{q}\gamma _\mu \gamma _5q|p,s=2\mathrm{\Delta }qs_\mu ,$$
(1)
where $`s^2=m^2`$ and the states are normalised according to $`p,s|p^{},s^{}=2E_p(2\pi )^3\delta (\stackrel{}{p}\stackrel{}{p}^{})\delta _{ss^{}}`$.

In the following we shall present preliminary results of a lattice computation of such matrix elements in the $`\mathrm{\Lambda }`$.

## 2 THE SIMULATION

We have performed quenched simulations with the Wilson gauge action at $`\beta =6.0`$ using nonperturbatively $`O(a)`$ improved fermions (clover fermions) with $`c_{\mathrm{SW}}=1.769`$. The lattice size was $`16^3\times 32`$. We worked with nine combinations of hopping parameters: $`\kappa _u=\kappa _d`$, $`\kappa _s`$ $`\{0.1324,0.1333,0.1342\}`$ corresponding to bare quark masses of $`166`$, $`112`$, $`58`$ MeV, respectively. So we can extra polate to the chiral limit in $`\kappa _u=\kappa _d`$ and inter polate in $`\kappa _s`$ to the physical value $`\kappa _s^{}`$. For the critical hopping parameter $`\kappa _c`$ we take the value 0.135201 determined from the PCAC quark mass , $`\kappa _s^{}=0.134138`$ was fixed by requiring the pseudoscalar mass $`m_{\mathrm{PS}}`$ for $`\kappa _u=\kappa _d=\kappa _c`$ and $`\kappa _s=\kappa _s^{}`$ to be equal to the $`K^+`$ mass $`=494`$ MeV (with the scale set by the Sommer parameter $`r_0=0.5`$ fm). In the extra- and interpolation we assumed a linear dependence of $`m_{\mathrm{PS}}^2`$ on the quark mass.

As an interpolating field for the $`\mathrm{\Lambda }`$ we used (employing Euclidean notation from now on)

$$\underset{x,x_4=t}{}ϵ_{ijk}s_i(x)\left(u_j^T(x)C\gamma _5d_k(x)\right)$$
(2)
with the charge conjugation matrix $`C`$ ($`i`$, $`j`$, $`k`$: colour indices). The quark fields were (Jacobi) smeared to improve the overlap with the $`\mathrm{\Lambda }`$ state.

As usual, the bare matrix elements are determined from ratios of three-point functions over two-point functions. In order to reduce the cut-off effects from $`O(a)`$ to $`O(a^2)`$ also in matrix elements, the improvement of the fermionic action has to be accompanied by the improvement of the operator under study. In our case, this is the axial vector current, and for the renormalised improved operator we can take

$$A_\mu ^\mathrm{R}=Z_A^0(1+b_Aam)(A_\mu +ac_A_\mu P)$$
(3)
with the bare axial vector current $`A_\mu (x)=\overline{q}(x)\gamma _\mu \gamma _5q(x)`$, the pseudoscalar density $`P(x)=\overline{q}(x)\gamma _5q(x)`$, and the bare quark mass $`am=1/(2\kappa )1/(2\kappa _c)`$. The improvement term $`_\mu P`$ vanishes in forward matrix elements, so we do not need the coefficient $`c_A`$. For $`Z_A^0`$ and $`b_A`$ we take the values from Ref. : $`Z_A^0=0.807(2)(8)`$, $`b_A=1.28(3)(4)`$.

## 3 RESULTS

From the $`\mathrm{\Lambda }`$ masses at our nine combinations of $`\kappa _d`$, $`\kappa _s`$ we have computed $`\mathrm{\Lambda }`$ masses at $`\kappa _d=\kappa _c`$ by linear extrapolation of $`m_\mathrm{\Lambda }^2`$ in $`1/\kappa _d`$. These 12 masses are plotted in Fig.1. Our value for $`\kappa _s^{}`$ reproduces quite accurately the ratio $`m_\mathrm{\Lambda }/m_p=1.19`$.

This fits in nicely with the observation that $`m_\mathrm{\Sigma }/m_p`$ and $`m_\mathrm{\Xi }/m_p`$ are also rather well reproduced by quenched simulations (see, e.g., ). The plot shows clearly the breaking of the SU(3) flavour symmetry. In particular, the dependence of the masses on $`\kappa _d`$ is rather pronounced.

In Fig.2 we plot our (preliminary) results for $`\mathrm{\Delta }s`$ and $`\mathrm{\Delta }d`$ in the $`\mathrm{\Lambda }`$ versus $`1/\kappa _s`$. In contrast with the case of $`m_\mathrm{\Lambda }/m_p`$, the dependence on $`\kappa _d`$ is rather weak. The values corresponding to $`\kappa _d=\kappa _c`$ have been obtained by extrapolating the bare matrix elements linearly in $`1/\kappa _d`$. Interpolating the bare matrix elements for $`\kappa _d=\kappa _c`$ linearly in $`1/\kappa _s`$ to $`1/\kappa _s^{}`$ we obtain the desired $`\mathrm{\Lambda }`$ matrix elements given in the last line of the following table:

|  | $`\mathrm{\Delta }u_\mathrm{\Lambda }=\mathrm{\Delta }d_\mathrm{\Lambda }`$ | $`\mathrm{\Delta }s_\mathrm{\Lambda }`$ |
| --- | --- | --- |
| quark model | $`0`$ | $`1`$ |
| exp. + SU(3)F | $`0.17(3)`$ | $`0.63(3)`$ |
| MC + SU(3)F | $`0.016(9)`$ | $`0.65(2)`$ |
| this work | $`0.01(4)`$ | $`0.67(3)`$ |

In the third line we give the $`\mathrm{\Lambda }`$ matrix elements as they follow from our Monte Carlo results for the proton matrix elements by the use of SU(3)F. They agree quite well with the matrix elements computed directly. Since the flavour symmetry breaking effects in the matrix elements are rather small (see Fig. 2), this is hardly surprising. The second line contains the values of the $`\mathrm{\Lambda }`$ matrix elements computed from the proton spin structure under the assumption of flavour SU(3)(see, e.g., ). All these results differ markedly from the predictions of the (naive) quark model shown in the first line.

Concerning the Monte Carlo results one should keep in mind that they were obtained at a fixed value of the lattice spacing. So a continuum extrapolation is not possible. Furthermore we used the quenched approximation and neglected quark-line disconnected contributions. Therefore it might be more reasonable to compare with estimates of the valence quark contribution, for which, e.g., Ashery and Lipkin  find $`\mathrm{\Delta }u_\mathrm{\Lambda }=\mathrm{\Delta }d_\mathrm{\Lambda }=0.07(4)`$, $`\mathrm{\Delta }s_\mathrm{\Lambda }=0.73(4)`$. Compared with the numbers from the second line of our table, $`\mathrm{\Delta }d_\mathrm{\Lambda }`$ gets closer to our Monte Carlo result, while $`\mathrm{\Delta }s_\mathrm{\Lambda }`$ becomes somewhat larger.

## 4 SUMMARY

Our Monte Carlo investigation of the spin structure of the $`\mathrm{\Lambda }`$ hyperon leads to satisfactory agreement with the phenomenological numbers, although the comparison is not straightforward due to the use of the quenched approximation. SU(3) flavour symmetry appears to be much less violated in the matrix elements of the axial vector current than in the baryon masses, in agreement with the empirical observation that the hyperon semileptonic decays can be parametrised rather well assuming SU(3) flavour symmetry. The relevance of our findings for using $`\mathrm{\Lambda }`$’s as a probe of the nucleon spin structure in deep-inelastic lepton-nucleon scattering remains to be studied. But we expect them to shed some light also on the lattice results for the nucleon spin structure.

## ACKNOWLEDGEMENTS

This work is supported by the DFG (Schwerpunkt “Elektromagnetische Sonden”) and by BMBF. The numerical calculations were performed on the Quadrics computers at DESY Zeuthen. We wish to thank the operating staff for their support.


text:
# Three-jet cross sections in hadron-hadron collisions at next-to-leading order

\]http://www.cpt.dur.ac.uk/ nagyz/nlo++.html

## Abstract

We present a new QCD event generator for hadron collider which can calculate one-, two- and three-jet cross sections at next-to-leading order accuracy. In this letter we study the transverse energy spectrum of three-jet hadronic events using the $`k_{}`$ algorithm. We show that the next-to-leading order correction significantly reduces the renormalization and factorization scale dependence of the three-jet cross section.

Perturbative QCD, Jet calculation, NLO  preprint: IPPP/01/48, DCPT/01/96, hep-ph/0110315

The latest version of the experiment at Tevatron and the future collider LHC will provide precise data so that not only inclusive measurements can be used to study the physics of hadronic final state. The studies of the event shapes and multi-jet event can be important projects.

One of the important theoretical tools in the analysis of hadronic final states is perturbative Quantum Chromodynamics (QCD). In order to make quantitative predictions in perturbative QCD, it is essential to perform the computations (at least) at the next-to-leading order (NLO) accuracy. In hadron collision the most easily calculated one- and two-jet cross sections have so far been calculated at NLO level Giele et al. (1993); Kunszt and Soper (1992). At next-to-leading level some three-jet observables were calculated by Giele and Kilgore Kilgore and Giele ; Kilgore and Giele (1997) and by Trócsányi Trócsányi (1996). In this letter we present a new NLO event generator for calculating jet observables in hadron-hadron collision. We compute the three-jet cross sections using the $`k_{}`$ algorithm Catani et al. (1993); Ellis and Soper (1993) to resolve jets in the final state. With our Monte Carlo program one can compute the NLO cross section for any other infrared safe one-, two- and three-jet observables. The presented distributions are given simply as illustration.

In the case of one-jet inclusive cross section in high $`p_T`$ region the forthcoming experimental data requires the knowledge of the higher order corrections. Recently, progress has been made in calculating the two loop $`22`$ matrix elements Anastasiou et al. (2001a, b); Glover and Tejeda-Yeomans (2001). These matrix element are need to set up a Monte Carlo program for calculating the next-to-next-to-leading order (NNLO) one- and two-jet cross section. However, the three-jet NLO calculation, comprising the one-loop $`23`$ and tree-level $`24`$ processes is a necessary first step of this project. A numerically stable and fast three-jet NLO program can provide enough hope that it might be possible to develop a numerically stable Monte Carlo program for calculating one- and two-jet cross sections at NNLO level.

In the last few years the theoretical developments make possible the next-to-leading order calculation for the three-jet quantities. There are several general methods available for the cancellation of the infrared divergences that can be used for setting up a Monte Carlo program Nagy and Trócsányi (1997); Frixione et al. (1996); Catani and Seymour (1997). In computing the NLO correction we use the dipole formalism of Catani-Seymour Catani and Seymour (1997) that we modified slightly in order to have a better control on the numerical calculation. The main idea is to cut the phase space of the dipole subtraction terms as introduced in Ref. Nagy and Trócsányi (1999), the details of how to this to apply for the case of hadron-hadron scattering will be given elsewhere.

The advantages of using the dipole method are the followings: i) no approximation is made; ii) the exact phase space factorization allows full control over the efficient generation of the phase space; iii) neither the use of color ordered subamplitudes, nor symmetrization, nor partial fractioning of the matrix elements is required; iv) Lorentz invariance is maintained, therefore, the switch between various frames can be achieved by simply transforming the momenta; v) the use of crossing functions is avoided; vi) it can be implemented in an actual program in a fully process independent way.

In this calculation we used the crossing symmetric tree- and one-loop level amplitudes. The parton subprocess $`0ggggg`$ Bern et al. (1993), $`0q\overline{q}ggg`$ Bern et al. (1995), $`0q\overline{q}Q\overline{Q}g`$ Kunszt et al. (1994) and the subprocesses related to these have been computed to one-loop and $`0gggggg`$, $`0q\overline{q}gggg`$, $`0q\overline{q}Q\overline{Q}gg`$ Gunion and Kunszt (1985, 1986b, 1986a); Kuijf (1991) and the crossed processes have been computed at tree level.

We have checked numerically that in all soft and collinear regions the difference of the real and subtraction terms contain only integrable square-root singularities. Furthermore, we have also checked that our results are independent of the parameter that controls the volume of the cut dipole phase space, which ensures that indeed the same quantity has been subtracted from the real correction as added to the virtual one.

Finally, to have a further check of the computation, Our method of implementing the dipole substraction terms allows for the construction of a process independent programming of QCD jet cross sections at the NLO accuracy. We use the same program structure, with trivial modifications, to compute one-, two- and three-jet cross sections.

In order to check the base structure of the program we compare the our inclusive one-jet NLO result to prediction of the program JETRAD Giele et al. (1993). In this comparison we compare the one-jet inclusive cross section using the $`k_{}`$ jet algorithm and $`\mathrm{MRSD}_{}^{}`$ parton distribution function Martin et al. (1993). We find good agreement between the two program. The differences are within the statistical error as Fig. 1 shows.

Historically, only the cone algorithm has been used to reconstruct the jet at hadron collider. In the three or higher jet calculation at NLO level the cone algorithm is not suitable because it has a lot of difficulties: an arbitrary procedure must be implemented to split and merge the overlapping cones, and an ad-hoc parameter $`R_{sep}`$ is required to accomodate the difference between the jet definitions at parton and detector level. To avoid this uncertainty we use the $`k_T`$ algorithm which has been developed by several groups Catani et al. (1993); Ellis and Soper (1993). Our implementation is based on the Ref. Ellis and Soper (1993). The algorithm starts with a list of the particles and the empty list of the jets.

1. For each particle (pseudo-particle) $`i`$ in the list, define

$$d_i=p_{T,i}^2.$$
(1)
For each pair $`(i,j)`$ of momenta $`(ij)`$, define

$$d_{ij}=\mathrm{min}(p_{T,i}^2,p_{T,j}^2)\frac{\mathrm{\Delta }R_{ij}^2}{D^2},$$
(2)
where $`\mathrm{\Delta }R_{ij}^2=(\eta _i\eta _j)^2+(\varphi _i\varphi _j)^2`$ is square of the angular separation which is expressed in the term of the pseudo-rapidity $`\eta _i`$ and the azimuth angle $`\varphi _i`$. $`D`$ is a free parameter. The usual choice of this parameter is $`D=1`$.
2. Find the minimum of all the $`d_i`$ and $`d_{ij}`$ and label it $`d_{\mathrm{min}}`$.
3. If $`d_{\mathrm{min}}`$ is $`d_{ij}`$, remove particles (pseudo-particles) $`i`$ and $`j`$ from the list and replace them with a new, merged pseudo-particle $`p_{(ij)}`$ given by the recombination scheme. In this paper we use the $`E`$ recombination scheme which is define the new pseudo-particle as the sum of the two particle

$$p_{(ij)}=p_i+p_j.$$
(3)
4. If $`d_{\mathrm{min}}`$ is $`d_i`$, remove particle (pseudo-particle) $`i`$ from the list of particle and add it to list of jets.
5. If any particles remain, go to step 1.

The algorithm produces a list of jets, each separated by $`\mathrm{\Delta }R_{ij}>D`$.

Once the phase space integrations are carried out, we write the NLO jet cross section in the following form:

$`\sigma _{AB}^{n\mathrm{jet}}`$ $`=`$ $`{\displaystyle \underset{a,b}{}}{\displaystyle 𝑑\eta _a𝑑\eta _bf_{a/A}(\eta _a,\mu _F^2)f_{b/B}(\eta _b,\mu _F^2)}`$ (4)
$`\times `$ $`\widehat{\sigma }_{ab}^{n\mathrm{jet}}[p_a,p_b,\alpha _s(\mu _R^2),\mu _R^2/Q_H^2,\mu _F^2/Q_H^2],`$
where $`f_{i/H}(\eta ,\mu _F^2)`$ represents the patron distribution function of the incoming hadron defined at the factorization scale $`\mu _F=x_FQ_H`$, $`\eta _{a,b}`$ is the fraction of the proton momentum carried by the scattered partons $`p_{a,b}`$, $`Q_H`$ is the hard scale that characterizes the parton scattering which could be $`E_T`$ of the jet, jet mass of the event, etc and $`\mu _R=x_RQ_H`$ is the renormalization scale.

Eq. (4) shows that using the dipole method one may either compute the full cross section at NLO accuracy including the convolution with the parton distribution functions, or simply the parton level cross section $`\widehat{\sigma }_{ab}`$ which can then be convoluted with the parton densities after the Monte Carlo integration. The latter procedure is the proper one if we are interested in measurement of the parton distribution functions (to avoid the recalculation of the Monte Carlo integrals after each step of the fitting iteration).

The three-jet cross sections presented here were calculated for TEVATRON collider in proton-antiproton collision at the center of mass energy $`\sqrt{s}=1800`$ GeV. We restrict the pseudo-rapidity range and the minimum transverse energy of the jets in laboratory frame to be

$$4<\eta _{\mathrm{jet}}<4,E_T>50\mathrm{GeV}.$$
(5)
We choose the transverse energy of the leading jets

$$Q_H=E_T^{(1)},$$
(6)
as the hard scattering scale.

In Fig. 2, we plotted differential cross section in the term of the the transverse energy of leading jet convoluted with the CTEQ5M1 parton distribution functions Lai et al. (2000) and using the two-loop formula for the strong coupling,

$`\alpha _s(\mu )={\displaystyle \frac{\alpha _s(M_Z)}{w(\mu )}}\left(1{\displaystyle \frac{\alpha _s(M_Z)}{2\pi }}{\displaystyle \frac{\beta _1}{\beta _0}}{\displaystyle \frac{\mathrm{ln}(w(\mu ))}{w(\mu )}}\right),`$ (7)
$`w(\mu )=1\beta _0{\displaystyle \frac{\alpha _s(M_Z)}{2\pi }}\mathrm{ln}\left({\displaystyle \frac{M_Z}{\mu }}\right),`$ (8)
where $`\alpha _s(M_Z)=0.118`$ and $`\beta _0=(11C_A4T_RN_f)/3`$, $`\beta _1=(17C_A^26C_FT_RN_f10C_AT_RN_f)/3`$, with $`N_f=5`$ flavors. For the leading order results we used the CTEQ5L distributions and the one-loop $`\alpha _s`$ ($`\alpha _s(M_Z)=0.127`$ and $`\beta _1=0`$ in Eq. (7)). In this figure the theoretical uncertainty of the three-jet cross section is shown. Over the wide range of the value the renormalization and factorization scale ($`0.5<x_{R,F}<2`$) this uncertainty is in the next-to-leading order result is much than in the Born level calculation.

In Fig. 3 we study the scale dependences of the three-jet cross section. The strong dependence on the renormalization scale observed at LO is significantly reduced. The factorization scale dependence is not significant at NLO and does not change much. Setting the two scales equal, $`\mu _R=\mu _F=\mu `$, we can observe a wide plato peaking around $`x_R=x_F=0.7`$.

The inset figure in Fig. 2 shows the $`K`$ factor (ratio of the three-jet cross section at NLO to that at LO accuracy), indicating the relative size of the correction. We can see the size of the NLO correction is between $`10`$ and $`25\%`$ for smaller values of the transverse energy and at the end of the spectra the size of the correction is almost zero. The error bars indicate the statistical error of the Monte Carlo calculation. Because of the strong logarithmic behavior of the cross section the Monte Carlo calculation is very sensitive to the “missed binning” (when a huge positive comes from the real term and the corresponding huge negative weight form the subtraction term are filled in different histogram bins).

In Fig. 4 we study the dependence of the differential cross section on parameter $`D`$. We plotted the ratios of the cross section

$$R_D(E_T^{(1)})=\frac{d\sigma _{p\overline{p}}^{3jet}(E_T^{(1)};D)}{dE_T^{(1)}}/\frac{d\sigma _{p\overline{p}}^{3jet}(E_T^{(1)};1)}{dE_T^{(1)}},$$
(9)
for three different values of the parameter $`D`$ ($`D=0.5,\mathrm{\hspace{0.17em}1.5},\mathrm{\hspace{0.17em}2}`$). The parameter $`D`$ controls the angular separation in the jet algorithm procedure in Eq. 2. Changing this parameter we expect more resolved jets (more three-jet events) with high transverse energy and less recombination for smaller values of $`D`$ and vice-versa. This behavior can be clearly observed from in Fig. 4.

In this letter we presented a NLO computation of the three-jet cross section defined with the $`k_{}`$ clustering algorithm in in hadron-hadron collision. Our results were obtained using a partonic Monte Carlo program that is suitable for implementing any detector cuts. We found that the NLO the correction is under $`30\%`$ in the case of differential cross section but the $`K`$ factor is sensitive to the allowed kinematic region. We demonstrated that the NLO corrections reduce the scale dependence significantly. We also presented how the differential cross section depends on the angular separation parameter $`D`$ used to define the jet. The same program can be used for computing the QCD radiative corrections to the (differential) cross section of any kind of one-, two-, or three-jet cross section or event-shape distribution in hadron-hadron collision. We compared the two-jet results obtained by our program to previous results and found agreement.

I thank Nigel Glover, Adrian Signer and Zoltán Trócsányi for the helpful discussions. This work was supported in part by the EU Fourth Framework Programme “Training and Mobility of Researchers”, Network “QCD and particle structure”, contract FMRX-CT98-0194 (DG 12 - MIHT), the EU Fifth Framework Programme ‘Improving Human Potential’, Research Training Network ‘Particle Physics Phenomenology at High Energy Colliders’, contract HPRN-CT-2000-00149 as well as by the Hungarian Scientific Research Fund grant OTKA T-025482.