text:
# Benchmarking MILC code with OpenMP and MPI

## 1 OpenMP and MPI

Open MultiProcessing (OpenMP or OMP ) and Message Passing Interface (MPI) are two strategies for using multiple processors for a single problem. The key difference between them is that in MPI, different nodes have their own memory and they communicate with each other when needed; but with OpenMP, the memory is shared between threads.

Here is an example. Suppose we have a two dimensional lattice with 4 sites in each direction, and we are using four nodes or threads, as shown below.
1 2 3 4   5 6 7 8   9 10 11 12   13 14 15 16

Suppose node/thread 1 corresponds to sites 1, 2, 5 and 6. In MPI, node 1 has information only about sites on that node, namely 1, 2, 5 and 6. If it needs information about other sites, for example about sites 3 or 7 which are nearest neighbors of sites 2 and 6 respectively, it has to use communication routines.

Contrast this with OpenMP, where all threads have access to data for all sites, but thread 1 does computations only for sites 1, 2, 5 and 6.

Table 1 summarizes the differences between the two strategies. A trend towards shared memory parallel machines or clusters of Symmetric Multiprocessing (SMP) nodes rather than the older paradigm of Massively Parallel Processing (MPP) machines makes a study of OpenMP parallelism timely. OpenMP was designed to exploit certain characteristics of shared-memory architectures. The ability to directly access memory throughout the system (with minimum latency and no explicit address mapping) combined with very fast shared memory locks, makes shared-memory architectures best suited for supporting OpenMP. The advantage of OpenMP is that it is easier to program. Unlike MPI, one does not have to worry about passing messages between nodes. In this paper, we study how OpenMP performs relative to MPI, and whether combining the two strategies gives better performance.

## 2 OpenMP DETAILS

In this section, we give some details about how a C code that works for a single processor is changed to work on multiple threads. The number of threads is determined by an environment variable, OMP\_NUM\_THREADS. The code is executed serially, on a single thread, until a parallel construct is encountered, which is executed on multiple threads and then serial execution is resumed. To define a parallel construct, lines beginning with #pragma omp are added to the code. Such pragmas are ignored by the usual C compiler, so the code may also be run as an ordinary serial code. They are, however, interpreted by an OpenMP compiler to identify parallel regions. There are several constructs that can be made to execute in parallel; here is an example for a “for” construct.
#pragma omp parallel for
for(i=0;i$`<`$N;i++){
my\_job(i);
}

This will run the function my\_job in parallel on different threads. Note that though the memory is shared, each thread must have a private copy of some variables, like i in the above example. Loop variables are made private by default but other such variables have to be declared “private”. Some variables may need to be summed over all the sites. This is accomplished with a reduction statement. The syntax is as follows:
j=0;
#pragma omp parallel for reduction(+:j)
for(i=0;i$`<`$N;i++){
j+=my\_function(i);
}

This is equivalent to
$`\text{j}=_{\text{i=0}}^{\text{N-1}}\text{my\_function(i)}`$.
Note that the sum is performed over all threads though each thread works only on part of the total number of iterations.

Identifying private and reduction variables is necessary for getting correct results.

## 3 MILC CODE AND CHANGES

The MILC  code is a set of publicly available codes developed by the MIMD Lattice Computation (MILC) collaboration for doing QCD simulations. This code has been run on a variety of parallel computers, using MPI, for many physics projects. The files are organized in different directories — the libraries directory contains low level routines like matrix multiplication, the generic directory contains oft-needed but somewhat higher level routines, including the updating and inversion routines. Then there are various application directories. For this project, we only concentrated on the conjugate gradient inverter, file d\_congrad5.c in generic\_ks directory in version 6 of MILC code. The code uses a macro FORALLSITES defined as
#define FORALLSITES(i,s)$`\backslash `$
for(i=0,s=lattice;i$`<`$sites\_on\_node;$`\backslash `$
i++,s++)
where lattice is an array of sites, site is a structure containing variables defined at each lattice point and sites\_on\_node is the number of lattice points on a given node. We needed to redefine the macro FORALLSITES because the OpenMP compiler we used could not deal with two variables (i and s) in a parallel for statement. Here is the macro redefinition.
#define FORALLSITES(i,s)$`\backslash `$
for(i=0;i$`<`$sites\_on\_node;i++){$`\backslash `$
s=&(lattice\[i\]);

We used another macro END\_LOOP, which is just defined to be a closing brace } to match the opening brace in the above macro.

## 4 COMPILING AND RUNNING

We used the KAP/Pro toolset  for this project. It includes the following:

* guidec: OpenMP compiler for C.
* guideview: OpenMP parallel performance visualization tool. It gives details of program execution, in particular, time spent in serial and parallel execution, imbalance in different regions of the code, etc.
* assurec: Compiler to be used with debugger which works by comparing single thread and multiple thread executions.
* assureview: OpenMP programming correctness tool for viewing details of errors or conflicts which occur if different threads try to read/write the same variables at the same time.

To add OpenMP parallelism to the MILC code, the following steps were required. First, we had to redefine the macro as explained above. Then we added the parallel for pragmas, specifying private and reduction variables. For example, in the FORALLSITES loop, s was made private. We changed cc to guidec in our makefiles, and we had to modify those compiler options that guidec did not recognize. Adding --backend before a compiler option forces guidec to use cc compiler options. Then, we ran assurec and assureview to locate and remove conflicts. Finally, we ran the executable on different number of threads and verified that the output agreed with the MPI output.

Even after one has a working OpenMP code, there are some issues to consider when comparing its performance with that of MPI. Some performance problems are OpenMP issues, while others are not. If single thread OpenMP performance does not match that with a single node under MPI, that may indicate a culprit other than OpenMP. For example, thread safe compilation requires the -mt switch on Sun. If using this switch on the original serial code decreases performance substantially, then the performance issue lies with the Sun compiler and its runtime libraries, not with OpenMP. If the code uses many malloc/free pairs then thread-safe memory allocation is likely the culprit. Again, this is not an OpenMP issue, but an issue with the quality of the vendor’s thread-safe compiler/runtime implementation. We verified that the -mt option on the serial version on Sun did not affect the performance significantly. Except for the case where we combine OpenMP and MPI, no malloc/free statements are used in the region of the code where the performance is evaluated. Thus, to the best of our knowledge, this is a fair comparison between OpenMP and MPI performance.

## 5 RESULTS FOR SUN E10000 AND BLUE HORIZON

We first ran the modified codes on a Sun E10000 at Indiana University. The details of architecture for this computer can be found online . Benchmarks were done for various lattice sizes and numbers of threads. As the number of threads increased, the lattice dimensions were increased to keep the volume per thread constant at $`L^4`$. The number of threads $`N`$ was increased from 1 to 16 by factors of 2. For example, for a given $`L`$, the lattice size for 2 threads is $`L^32L`$ and for 16 threads it is $`(2L)^4`$. $`L`$ was increased from 4 to 14 in steps of 2. Reported in Fig. 1 is the performance on the Kogut-Susskind quark conjugate gradient routine in Megaflop/s per CPU for both OMP and MPI. For smaller number of threads/nodes, the OMP rates are quite comparable to MPI. They lag behind for larger number of threads. There are two factors involved here: the overhead for setting up threads and the use of the cache. For small lattice sizes, since there is only a small number of computations to be performed, the former degrades the performance, but if significant portion of the problem can fit in the cache, the execution is speeded up. On the other hand, for larger lattices the thread initialization overhead is a much smaller fraction of the total computation time, but the problem size is too big to fit into the cache. We see that OMP has a “sweet spot” at size 6, much as the MPI performance peaks at size 8. Since we keep the load per thread constant, for the same lattice size the performance monotonically decreases in most cases as we increase the number of threads.

Next we benchmarked the code on Blue Horizon . This IBM SP machine at the San Diego Supercomputer Center has 8-way SMP nodes but with the current switch can support only 4 MPI processes per node. Figure 2 contains the preliminary results from Blue Horizon. These results are qualitatively similar to the E10000 results.

## 6 COMBINING OMP AND MPI

A hybrid approach combining OpenMP parallelism within MPI processes may offer better performance than either individual approach. We tried different combinations of threads and MPI processes on Blue Horizon. The hybrid approach fared better at times. Figure 3 shows the results for a total of eight processors. It can be seen again that the MPI performance peaks at size 8 (the left-most bars in Fig. 3) and the OMP at size 6 (the right-most bars). The combination of 2 threads and 4 nodes works best for smaller sizes. The processors on Blue Horizon were upgraded after these runs. We should repeat these calculations and extend the study to a larger number of CPUs.

## 7 CONCLUSION

On both computers studied, OpenMP performance was very similar to MPI performance for a small number of threads, but it deteriorated much faster as the number of threads increased, for smaller lattice sizes. Thus, OpenMP may be a viable option for someone writing a code to be used with a modest number of processors on SMP machines. The MILC Collaboration, however, already has a working MPI code that scales well on many machines. For almost all the combinations of problem sizes and number of CPUs studied in this paper, MPI is at least as good as OpenMP, if not better. The only case where we get a considerable improvement over MPI is when we combine OpenMP and MPI on Blue Horizon for $`L=4`$ and 6. Not only does the hybrid approach give the best performance on a single SMP node, it should allow us to run multi-node jobs using all eight processors on each node rather than the limit of four with the current switch. We have added OpenMP parallelism to the MILC code only for the conjugate gradient inverter for this test project. It will require considerably more effort to modify the whole code to run on multiple OpenMP threads.

It is our pleasure to thank Bill Magro and Henry Gabbs at KAI for many useful discussions. We gratefully acknowledge the help provided by the staff at Research and Technical Services at IU, especially David Hart, Mary Papakhian and Stephanie Burks. This work was supported by the DOE under grant DE-F002-91ER 40661. We thank the San Diego Supercomputer Center and NPACI for use of Blue Horizon.


format: text

text:
# A Measurement of the Decay Asymmetry Parameters in 𝚵_𝒄^𝟎→𝚵⁻𝝅⁺

## Abstract

Using the CLEO II detector at the Cornell Electron Storage Ring we have measured the $`\mathrm{\Xi }_c^0`$ decay asymmetry parameter in the decay $`\mathrm{\Xi }_c^0\mathrm{\Xi }^{}\pi ^+`$. We find $`\alpha _{\mathrm{\Xi }_c^0}\alpha _\mathrm{\Xi }=0.26\pm 0.18\text{(stat)}_{0.04}^{+0.05}\text{(syst)}`$, using the world average value of $`\alpha _\mathrm{\Xi }=0.456\pm 0.014`$ we obtain $`\alpha _{\mathrm{\Xi }_c^0}=0.56\pm 0.39\text{(stat)}_{0.09}^{+0.10}\text{(syst)}`$. The physically allowed range of a decay asymmetry parameter is $`1<\alpha <+1`$. Our result prefers a negative value: $`\alpha _{\mathrm{\Xi }_c^0}`$ is $`<0.1`$ at the 90% CL. The central value occupies the middle of the theoretically expected range but is not yet precise enough to choose between models.

preprint:   CLNS 98/1586 CLEO 98-15

S. Chan,1 G. Eigen,1 E. Lipeles,1 J. S. Miller,1 M. Schmidtler,1 A. Shapiro,1 W. M. Sun,1 J. Urheim,1 A. J. Weinstein,1 F. Würthwein,1 D. E. Jaffe,2 G. Masek,2 H. P. Paar,2 E. M. Potter,2 S. Prell,2 V. Sharma,2 D. M. Asner,3 A. Eppich,3 J. Gronberg,3 T. S. Hill,3 D. J. Lange,3 R. J. Morrison,3 H. N. Nelson,3 T. K. Nelson,3 D. Roberts,3 B. H. Behrens,4 W. T. Ford,4 A. Gritsan,4 H. Krieg,4 J. Roy,4 J. G. Smith,4 J. P. Alexander,5 R. Baker,5 C. Bebek,5 B. E. Berger,5 K. Berkelman,5 V. Boisvert,5 D. G. Cassel,5 D. S. Crowcroft,5 M. Dickson,5 S. von Dombrowski,5 P. S. Drell,5 K. M. Ecklund,5 R. Ehrlich,5 A. D. Foland,5 P. Gaidarev,5 L. Gibbons,5 B. Gittelman,5 S. W. Gray,5 D. L. Hartill,5 B. K. Heltsley,5 P. I. Hopman,5 D. L. Kreinick,5 T. Lee,5 Y. Liu,5 N. B. Mistry,5 C. R. Ng,5 E. Nordberg,5 M. Ogg,5,***Permanent address: University of Texas, Austin TX 78712. J. R. Patterson,5 D. Peterson,5 D. Riley,5 A. Soffer,5 B. Valant-Spaight,5 A. Warburton,5 C. Ward,5 M. Athanas,6 P. Avery,6 C. D. Jones,6 M. Lohner,6 C. Prescott,6 A. I. Rubiera,6 J. Yelton,6 J. Zheng,6 G. Brandenburg,7 R. A. Briere,7Permanent address: Carnegie Mellon University, Pittsburgh, Pennsylvania 15213. A. Ershov,7 Y. S. Gao,7 D. Y.-J. Kim,7 R. Wilson,7 T. E. Browder,8 Y. Li,8 J. L. Rodriguez,8 H. Yamamoto,8 T. Bergfeld,9 B. I. Eisenstein,9 J. Ernst,9 G. E. Gladding,9 G. D. Gollin,9 R. M. Hans,9 E. Johnson,9 I. Karliner,9 M. A. Marsh,9 M. Palmer,9 M. Selen,9 J. J. Thaler,9 K. W. Edwards,10 A. Bellerive,11 R. Janicek,11 P. M. Patel,11 A. J. Sadoff,12 R. Ammar,13 P. Baringer,13 A. Bean,13 D. Besson,13 D. Coppage,13 R. Davis,13 S. Kotov,13 I. Kravchenko,13 N. Kwak,13 L. Zhou,13 S. Anderson,14 Y. Kubota,14 S. J. Lee,14 R. Mahapatra,14 J. J. O’Neill,14 R. Poling,14 T. Riehle,14 A. Smith,14 M. S. Alam,15 S. B. Athar,15 Z. Ling,15 A. H. Mahmood,15 S. Timm,15 F. Wappler,15 A. Anastassov,16 J. E. Duboscq,16 K. K. Gan,16 C. Gwon,16 T. Hart,16 K. Honscheid,16 H. Kagan,16 R. Kass,16 J. Lee,16 J. Lorenc,16 H. Schwarthoff,16 A. Wolf,16 M. M. Zoeller,16 S. J. Richichi,17 H. Severini,17 P. Skubic,17 A. Undrus,17 M. Bishai,18 S. Chen,18 J. Fast,18 J. W. Hinson,18 N. Menon,18 D. H. Miller,18 E. I. Shibata,18 I. P. J. Shipsey,18 S. Glenn,19 Y. Kwon,19,Permanent address: Yonsei University, Seoul 120-749, Korea. A.L. Lyon,19 S. Roberts,19 E. H. Thorndike,19 C. P. Jessop,20 K. Lingel,20 H. Marsiske,20 M. L. Perl,20 V. Savinov,20 D. Ugolini,20 X. Zhou,20 T. E. Coan,21 V. Fadeyev,21 I. Korolkov,21 Y. Maravin,21 I. Narsky,21 R. Stroynowski,21 J. Ye,21 T. Wlodek,21 M. Artuso,22 E. Dambasuren,22 S. Kopp,22 G. C. Moneti,22 R. Mountain,22 S. Schuh,22 T. Skwarnicki,22 S. Stone,22 A. Titov,22 G. Viehhauser,22 J.C. Wang,22 S. E. Csorna,23 K. W. McLean,23 S. Marka,23 Z. Xu,23 R. Godang,24 K. Kinoshita,24,§§§Permanent address: University of Cincinnati, Cincinnati OH 45221 I. C. Lai,24 P. Pomianowski,24 S. Schrenk,24 G. Bonvicini,25 D. Cinabro,25 R. Greene,25 L. P. Perera,25 and G. J. Zhou25

1California Institute of Technology, Pasadena, California 91125

2University of California, San Diego, La Jolla, California 92093

3University of California, Santa Barbara, California 93106

4University of Colorado, Boulder, Colorado 80309-0390

5Cornell University, Ithaca, New York 14853

6University of Florida, Gainesville, Florida 32611

7Harvard University, Cambridge, Massachusetts 02138

8University of Hawaii at Manoa, Honolulu, Hawaii 96822

9University of Illinois, Urbana-Champaign, Illinois 61801

10Carleton University, Ottawa, Ontario, Canada K1S 5B6

and the Institute of Particle Physics, Canada

11McGill University, Montréal, Québec, Canada H3A 2T8

and the Institute of Particle Physics, Canada

12Ithaca College, Ithaca, New York 14850

13University of Kansas, Lawrence, Kansas 66045

14University of Minnesota, Minneapolis, Minnesota 55455

15State University of New York at Albany, Albany, New York 12222

16Ohio State University, Columbus, Ohio 43210

17University of Oklahoma, Norman, Oklahoma 73019

18Purdue University, West Lafayette, Indiana 47907

19University of Rochester, Rochester, New York 14627

20Stanford Linear Accelerator Center, Stanford University, Stanford, California 94309

21Southern Methodist University, Dallas, Texas 75275

22Syracuse University, Syracuse, New York 13244

23Vanderbilt University, Nashville, Tennessee 37235

24Virginia Polytechnic Institute and State University, Blacksburg, Virginia 24061

25Wayne State University, Detroit, Michigan 48202

The weak interactions underlying hadronic charm quark weak decay are straightforward to describe theoretically, but complications arise because the quarks are bound inside hadrons by the strong force. These interactions, which are described by the theory of quantum chromodynamics (QCD), are very difficult to predict using perturbative methods because the strong coupling is large at the typical energies of charm decays. Compared to charm mesons, charm baryons offer new information for two reasons: first, non-factorizable contributions to the decay amplitude are important; W-exchange diagrams can contribute without the helicity suppression that decreases their contribution to pseudoscalar meson decays, and internal W emission is expected to be significant. The relative importance of these effects is believed to be responsible for both the observed lifetime hierarchy of the charm baryons and the relatively shorter lifetime of charm baryons compared to charm mesons. Second, parity violation in charm baryon decays is readily observable because the decay of the daughter hyperon also violates parity. In consequence a variety of models  are used to predict both the decay rate and degree of parity violation in charm baryon decays. To constrain the models it is important to provide as much experimental information as possible.

Parity violation occurs in hyperon and charm baryon $`\frac{1}{2}^+\frac{1}{2}^+0`$ decays due to the existence of two orbital angular momentum decay amplitudes of opposite parity. The experimental observable is an asymmetry in the angular decay distribution due to interference between the two amplitudes. In the decays $`\mathrm{\Lambda }_c^+\mathrm{\Lambda }\pi ^+`$ and $`\mathrm{\Xi }^{}\mathrm{\Lambda }\pi ^{}`$ both followed by $`\mathrm{\Lambda }p\pi ^{}`$, the $`\mathrm{\Lambda }`$ is produced with a polarization equal to

$$𝐏_𝚲=\frac{(\alpha _B+\widehat{𝚲}𝐏_𝐁)\widehat{𝚲}\beta _B(\widehat{𝚲}\times 𝐏_𝐁)\gamma _B\widehat{𝚲}\times (\widehat{𝚲}\times 𝐏_𝐁)}{(1+\alpha _B\widehat{𝚲}𝐏_𝐁)}$$
(1)
where $`𝐏_𝐁`$ is the parent baryon polarization, $`\alpha _B,\beta _B`$ and $`\gamma _B`$ are the parent baryon asymmetry parameters and $`\widehat{𝚲}`$ is a unit vector along the $`\mathrm{\Lambda }`$ momentum in the parent baryon frame . If the parent baryon polarization is unobserved, or if the parent baryon is not polarized, Equation 1 reduces to $`𝐏_𝚲=\alpha _B\widehat{𝚲}`$. The angular distribution of a proton from the decay of a $`\mathrm{\Lambda }`$ is therefore

$$\frac{dN}{d\mathrm{cos}\theta _\mathrm{\Lambda }}1+\alpha _B\alpha _\mathrm{\Lambda }\mathrm{cos}\theta _\mathrm{\Lambda }$$
(2)
where $`\theta _\mathrm{\Lambda }`$ is the angle between the proton momentum vector in the $`\mathrm{\Lambda }`$ rest frame and $`\widehat{𝚲}`$, and $`\alpha _\mathrm{\Lambda }`$ is the $`\mathrm{\Lambda }`$ decay asymmetry, which measures the degree of parity violation in $`\mathrm{\Lambda }`$ decay. Experimentally, a fit to equation 2 determines a product of decay asymmetries. For example, in $`\mathrm{\Xi }^{}\mathrm{\Lambda }\pi ^{}`$, $`\alpha _B=\alpha _\mathrm{\Xi }^{}`$; since $`\alpha _\mathrm{\Lambda }`$ is well measured ($`\alpha _\mathrm{\Lambda }=0.642\pm 0.013`$  ), $`\alpha _\mathrm{\Xi }^{}`$ can be extracted.

Equation 2 has also been used to extract the decay asymmetry of the $`\mathrm{\Lambda }_c`$ in several decay modes. For example in $`\mathrm{\Lambda }_c^+\mathrm{\Lambda }\pi ^+`$, $`\alpha _B=\alpha _{\mathrm{\Lambda }_c^+}`$. In the massless fermion limit the chirality of the weak interaction corresponds to $`\alpha =1`$, however, the interplay between the strong and weak interactions in the weak decay alters the value of $`\alpha `$. The decay asymmetry in $`\mathrm{\Lambda }_c^+\mathrm{\Lambda }\pi ^+`$ is relatively well measured, the world average value is $`\alpha _{\mathrm{\Lambda }_c^+}=0.98\pm 0.19`$ , which is consistent with the naive V-A expectation. This result had been predicted by Heavy Quark Effective Theory (HQET) in conjunction with the factorization hypothesis . Most pole and quark models , also predict large negative values of $`\alpha _{\mathrm{\Lambda }_c^+}`$.

In this letter we measure the $`\mathrm{\Xi }_c^0`$ decay asymmetry parameter in the decay $`\mathrm{\Xi }_c^0\mathrm{\Xi }^{}\pi ^+`$ for the first time. For $`\mathrm{\Xi }_c^0\mathrm{\Xi }^{}\pi ^+`$ theoretical predictions cover a larger range: $`1<\alpha _{\mathrm{\Xi }_c^0}<0.38`$ . The decay of $`\mathrm{\Xi }_c^0`$ is a three step process where parity violation occurs at each decay stage: $`\mathrm{\Xi }_c^0\mathrm{\Xi }^{}\pi ^+`$, $`\mathrm{\Xi }^{}\mathrm{\Lambda }\pi ^{}`$, $`\mathrm{\Lambda }p\pi ^{}`$. The differential rate is therefore given by :

$`{\displaystyle \frac{dN}{d\mathrm{cos}\theta _\mathrm{\Xi }d\mathrm{cos}\theta _\mathrm{\Lambda }}}1`$ $`+`$ $`\alpha _{\mathrm{\Xi }_c^0}\alpha _\mathrm{\Xi }\mathrm{cos}\theta _\mathrm{\Xi }+\alpha _\mathrm{\Lambda }\alpha _\mathrm{\Xi }\mathrm{cos}\theta _\mathrm{\Lambda }`$ (3)
$`+`$ $`\alpha _\mathrm{\Lambda }\alpha _{\mathrm{\Xi }_c^0}\mathrm{cos}\theta _\mathrm{\Xi }\mathrm{cos}\theta _\mathrm{\Lambda }`$ (4)
where in addition to the terms previously defined: $`\theta _\mathrm{\Xi }`$ is the angle between the $`\mathrm{\Lambda }`$ momentum vector in the $`\mathrm{\Xi }^{}`$ rest frame and the $`\mathrm{\Xi }^{}`$ momentum vector in the $`\mathrm{\Xi }_c^0`$ rest frame as shown in Fig. 1. Equation 4 reduces to the familiar form of Equation 2 when integrated over $`\mathrm{cos}\theta _\mathrm{\Lambda }`$.

The data sample in this study was collected with the CLEO II detector  at the Cornell Electron Storage Ring (CESR). The integrated luminosity consists of 4.83 $`fb^1`$ taken at and just below the $`\mathrm{{\rm Y}}(4S)`$ resonance, corresponding to approximately 5 million $`e^+e^{}c\overline{c}`$ events.

We search for the decay $`\mathrm{\Xi }_c^0\mathrm{\Xi }^{}\pi ^+`$ in $`e^+e^{}c\overline{c}`$ events by reconstructing the decay chain $`\mathrm{\Xi }_c^0\mathrm{\Xi }^{}\pi ^+`$, $`\mathrm{\Xi }^{}\mathrm{\Lambda }\pi ^{}`$, $`\mathrm{\Lambda }p\pi ^{}`$***Unless otherwise noted, throughout this letter charge conjugation is implied.

The $`\mathrm{\Lambda }`$ is reconstructed by requiring two oppositely charged tracks to originate from a common vertex. The positive track is required to be consistent with a proton hypothesisHadronic particles are identified by requiring specific ionization energy loss measurements ($`dE/dx`$), combined with time-of-flight (TOF) information when available. The two measurements are combined into a joint probability for the particle to be a pion, a kaon or a proton. A charged track is defined to be consistent with a particle hypothesis if its probability is greater than 0.003.. The momentum of the $`\mathrm{\Lambda }`$ candidate is calculated by extrapolating the charged track momenta to the secondary vertex. The invariant mass of $`\mathrm{\Lambda }`$ candidates is required to be within three standard deviations (3$`\sigma `$ = 6.0 MeV/c2) of the known $`\mathrm{\Lambda }`$ mass. Track combinations which satisfy interpretation as $`K_s^0\pi ^+\pi ^{}`$ are rejected. Combinatoric and $`B`$ decay backgrounds are reduced by requiring the momentum of $`\mathrm{\Lambda }`$ candidates be greater than 800 MeV/c.

The $`\mathrm{\Xi }^{}`$ is reconstructed in the $`\mathrm{\Lambda }\pi ^{}`$ decay mode. $`\mathrm{\Xi }^{}`$ candidates are formed by combining each $`\mathrm{\Lambda }`$ candidate with a negatively charged track consistent with a pion hypothesis. The $`\mathrm{\Xi }^{}`$ candidate vertex is formed from the intersection of the $`\mathrm{\Lambda }`$ momentum vector and the negatively charged track. To obtain the $`\mathrm{\Xi }^{}`$ momentum, and $`\mathrm{\Lambda }\pi `$ invariant mass, the momentum of the charged track is recalculated at the new vertex. The invariant $`\mathrm{\Lambda }\pi ^{}`$ mass is fit to a double Gaussian signal shape with parameters fixed from a GEANT  based Monte Carlo (MC) simulation of the detector, and a first order Chebyshev polynomial to describe the combinatorial background. We find $`8434\pm 109`$ events consistent with $`\mathrm{\Xi }^{}\mathrm{\Lambda }\pi ^{}`$ with an invariant mass of $`1322.0\pm 0.03`$ MeV/c2 and width $`\sigma _{av}=2.5`$ MeV/c2, where $`\sigma _{av}`$ is the weighted average width of the double Gaussian. The mean and width are in agreement with the MC simulation.

To reconstruct $`\mathrm{\Xi }_c^0`$ candidates all $`\mathrm{\Xi }^{}`$ candidates are combined with positively charged tracks consistent with a pion hypothesis. Since charm fragmentation is a relatively hard process, the $`\mathrm{\Xi }_c^0`$ spectrum from $`e^+e^{}c\overline{c}`$ is also fairly stiff. We therefore use a scaled momentum cut, $`x_p>0.5`$ $`x_p=p/p_{max}=p/\sqrt{E_{beam}^2m_{\mathrm{\Xi }_c^0}^2}`$  to reduce combinatoric background. This cut also eliminates charmed baryons from decays of B mesons. Fig. 2 shows the invariant mass distribution of $`\mathrm{\Xi }^{}\pi ^+`$ combinations. The $`\mathrm{\Xi }^{}\pi ^+`$ invariant mass distribution is fit to a double Gaussian to describe the signal and a first order Chebyshev polynomial to describe the combinatorial background, where the parameters of the double Gaussian are fixed by MC simulation. The shaded histogram is wrong sign (WS) random $`\mathrm{\Xi }^{}\pi ^{}`$ combinations. There are less WS than right sign (RS) random combinations due to charge conservation.

The excess of right sign events over wrong sign events below M($`\mathrm{\Xi }^{}\pi ^+`$) $`=2.33`$ GeV/c2 is due to feedthrough from $`\mathrm{\Xi }_c^{+,0}\mathrm{\Xi }^{}n(\pi ),n>1,`$ where one or more pions are missing. The kinematic limit for feedthrough is $`M(\mathrm{\Xi }_c^0)M(\pi )=2.3304`$ GeV/c2. To simplify the fit to M($`\mathrm{\Xi }^{}\pi ^+`$) the feedthrough region is excluded. We find $`138\pm 14`$ $`\mathrm{\Xi }_c^0\mathrm{\Xi }^{}\pi ^+`$ candidates with a mass of $`2470.6\pm 1.0`$ MeV/c2 and width $`\sigma _{av}=10`$ MeV/c2. The mean and width are in agreement with MC simulation. We require that the invariant mass of $`\mathrm{\Xi }_c^0`$ candidates be within 30 MeV/c2 (3 $`\sigma _{av}`$) of the known $`\mathrm{\Xi }_c^0`$ mass.

MC simulation shows that the difference between generated and reconstructed values of $`\mathrm{cos}\theta _\mathrm{\Xi }`$ and $`\mathrm{cos}\theta _\mathrm{\Lambda }`$ has both Gaussian and symmetric non-Gaussian components. The resolution in $`\mathrm{cos}\theta _{\mathrm{\Xi },\mathrm{\Lambda }}`$, $`\sigma (\mathrm{cos}\theta _\mathrm{\Xi },\mathrm{\Lambda })`$, defined to be the average of the rms variance of the Gaussian and of the non-Gaussian component, weighted by the relative normalizations of the two components §§§For $`\mathrm{cos}\theta _\mathrm{\Xi }`$ and $`\mathrm{cos}\theta _\mathrm{\Lambda }`$, the non-Gaussian component comprises 60% and 32% of the distribution, respectively is: $`\sigma (\mathrm{cos}\theta _\mathrm{\Xi })=0.014`$ and $`\sigma (\mathrm{cos}\theta _\mathrm{\Lambda })=0.030`$.

The decay asymmetry parameter in $`\mathrm{\Xi }_c^0\mathrm{\Xi }^{}\pi ^+`$ is measured in a two dimensional unbinned maximum loglikelihood fit to the two-fold decay angular distribution of Equation 4, in a manner similar to . This technique enables a multi-dimensional likelihood fit to be performed to variables modified by experimental acceptance and resolution. The probability function of the signal, $`\mathrm{\Gamma }_s,`$ is determined by generating one high statistics MC sample of $`\mathrm{\Xi }_c^0\mathrm{\Xi }^{}\pi ^+`$, $`\mathrm{\Xi }^{}\mathrm{\Lambda }\pi ^{}`$, $`\mathrm{\Lambda }p\pi ^{}`$ with a known value of $`\alpha _{\mathrm{\Xi }_c^0}`$ and the world average values of $`\alpha _\mathrm{\Xi }`$ and $`\alpha _\mathrm{\Lambda }`$. The generated events are processed through the detector simulation, off-line analysis programs, and selection criteria. Using the generated angles, accepted MC events are weighted by the ratio of the decay distribution for a trial value of $`\alpha _{\mathrm{\Xi }_c^0}`$ to that of the generated distribution. By such weighting, a likelihood may be evaluated for each data event for trial values of $`\alpha _{\mathrm{\Xi }_c^0}`$, and a fit performed. The probability for each event is determined by sampling $`\mathrm{\Gamma }_s`$ using a search area centered on each data point. The size of the area is chosen such that the systematic effect from a finite search area is small and the required number of MC events is not prohibitively high.

Background is incorporated into the fitting technique by constructing the log-likelihood function:

$$\mathrm{ln}=\underset{i=1}{\overset{N}{}}ln(P_S\mathrm{\Gamma }_S+P_B\mathrm{\Gamma }_B)$$
(5)
where $`N`$ is the number of events in the signal region and $`P_S`$ and $`P_B`$ are the probabilities that events in this region are signal and background respectively. The probability distribution of background in the signal region, $`\mathrm{\Gamma }_B`$, is determined from ($`\mathrm{\Xi }^{}\pi ^+`$) mass sidebands above and below the signal region. The sidebands are: 2.5004 $`<`$ M($`\mathrm{\Xi }^{}\pi ^+`$)$`<`$2.6904 GeV/c2 and 2.3304 $`<`$ M($`\mathrm{\Xi }^{}\pi ^+`$)$`<`$2.4404 GeV/c2.

CP conservation requires $`\alpha _{\mathrm{\Xi }_c^0}=\alpha _{\overline{\mathrm{\Xi }_c^0}}`$. As the number of signal events is small the analysis is insensitive to the presence of CP violation, therefore CP is assumed to be conserved. Since Equation 4 depends on the products $`\alpha _{\mathrm{\Xi }_c^0}\alpha _\mathrm{\Xi }`$, $`\alpha _{\mathrm{\Xi }_c^0}\alpha _\mathrm{\Lambda }`$ and $`\alpha _\mathrm{\Xi }\alpha _\mathrm{\Lambda }`$ which have the same sign for particle and anti-particle states, particle and anti-particle distributions are combined.

The validity of the analysis procedure is determined by constructing artificial data sets consisting of MC signal events, generated with known asymmetries, and background events taken from data. The generated asymmetry was varied over the range: $`1<\alpha _{\mathrm{\Xi }_c^0}<0`$, in all cases the fit returns unbiased values of $`\alpha _{\mathrm{\Xi }_c^0}`$.

Applying the fit procedure to the data we find $`\alpha _{\mathrm{\Xi }_c^0}\alpha _\mathrm{\Xi }=0.26\pm 0.18`$, $`\alpha _\mathrm{\Xi }\alpha _\mathrm{\Lambda }=0.57\pm 0.21`$ and $`\alpha _{\mathrm{\Xi }_c^0}\alpha _\mathrm{\Lambda }=0.26\pm 0.37`$. The product $`\alpha _{\mathrm{\Xi }_c^0}\alpha _\mathrm{\Lambda }`$ is second order in decay angles, consequently the fit is less sensitive to this quantity than to $`\alpha _{\mathrm{\Xi }_c^0}\alpha _\mathrm{\Xi }`$ and $`\alpha _\mathrm{\Xi }\alpha _\mathrm{\Lambda }`$. We provide a value of $`\alpha _{\mathrm{\Xi }_c^0}\alpha _\mathrm{\Lambda }`$ for completeness only. The $`\mathrm{cos}\theta _\mathrm{\Xi }`$ and $`\mathrm{cos}\theta _\mathrm{\Lambda }`$ distributions for data and projections of the fit are shown in Fig. 3.

The value of $`\alpha _\mathrm{\Xi }\alpha _\mathrm{\Lambda }`$ found in the fit is in reasonable agreement with a CLEO measurement of $`\alpha _\mathrm{\Xi }\alpha _\mathrm{\Lambda }=0.291\pm 0.021\text{(stat)}_{0.015}^{+0.019}\text{(syst)}`$ obtained with approximately 8,000 $`\mathrm{\Xi }^{}\mathrm{\Lambda }\pi ^{}`$ events  and with the current world average value $`\alpha _\mathrm{\Xi }\alpha _\mathrm{\Lambda }=0.293\pm 0.007`$ . Constraining $`\alpha _\mathrm{\Xi }\alpha _\mathrm{\Lambda }`$ to the world average value does not reduce the statistical error on $`\alpha _{\mathrm{\Xi }_c^0}\alpha _\mathrm{\Xi }`$.

We have considered the following sources of systematic error and give our estimate of the percentage error on $`\alpha _{\mathrm{\Xi }_c^0}\alpha _\mathrm{\Xi }`$ in parentheses. The statistical error in the MC sample is estimated by varying the size of the MC sample used in the fit $`(\pm 2\%)`$. The error associated with the uncertainty in the $`\mathrm{\Xi }_c^0`$ fragmentation function  Based on a measurement of the $`\mathrm{\Xi }_c^0`$ fragmentation function at CLEO, we use a Peterson fragmentation function  with $`ϵ_q=0.15`$. is estimated by varying this function $`(+1\%,14\%)`$. To determine the effect of incomplete knowledge of the background shape and the effect of statistical fluctuations in the sideband sample used to model the background in the signal region, we vary the size of the lower and upper sidebands used in the loglikelihood fit. We also repeat the fit using wrong sign $`\mathrm{\Xi }^{}\pi ^{}`$ events in both the signal mass region and sideband regions to model the background shape $`(+16\%,1\%)`$. The error associated with MC modeling of slow pions from $`\mathrm{\Xi }`$ and $`\mathrm{\Lambda }`$ decay is obtained by varying the reconstruction efficiency according to our understanding of the CLEO II detector $`(+9\%,6\%)`$. Using a large sample of $`\mathrm{\Xi }_c^0\mathrm{\Xi }^{}\pi ^+`$ MC generated with $`\alpha _{\mathrm{\Xi }_c^0}=0.5`$ and the world average values of $`\alpha _\mathrm{\Xi }`$ and $`\alpha _\mathrm{\Lambda }`$ and including a randomly generated background based on the shape of the data sideband, we measure the effect of varying the size of the area element used to determine $`\mathrm{\Gamma }_S(\mathrm{cos}\theta _\mathrm{\Xi },\mathrm{cos}\theta _\mathrm{\Lambda };\alpha _{\mathrm{\Xi }_c^0},\alpha _\mathrm{\Xi },\alpha _\mathrm{\Lambda })`$ and the background shape $`\mathrm{\Gamma }_B(\mathrm{cos}\theta _\mathrm{\Xi },\mathrm{cos}\theta _\mathrm{\Lambda })`$ in the loglikelihood fit $`(+1\%,2\%)`$. The fit method is checked by integrating over $`\mathrm{cos}\theta _\mathrm{\Lambda }`$ ($`\mathrm{cos}\theta _\mathrm{\Xi }`$) and performing a one dimensional binned fit to $`\mathrm{cos}\theta _\mathrm{\Xi }`$ ($`\mathrm{cos}\theta _\mathrm{\Lambda }`$). The results are consistent with the maximum likelihood fit. Possible background from $`\mathrm{\Sigma }_c^0(2450)\mathrm{\Lambda }_c^+\pi _{\mathrm{slow}}^{}`$, $`\mathrm{\Lambda }_c^+\mathrm{\Lambda }\pi _{\mathrm{fast}}^+`$ is determined to be negligible by MC simulation. This measurement is insensitive to production polarization, $`𝐏_{𝚵_𝐜^\mathrm{𝟎}}`$, and no systematic error has been included from this source . Adding all sources of systematic error in quadrature $`(+19\%,15\%)`$, we find $`\alpha _{\mathrm{\Xi }_c^0}\alpha _\mathrm{\Xi }=0.26\pm 0.18\text{(stat)}_{0.04}^{+0.05}\text{(syst)}`$. A similar study for the systematic uncertainty on the measurement of $`\alpha _\mathrm{\Xi }\alpha _\mathrm{\Lambda }`$ in $`\mathrm{\Xi }_c^0\mathrm{\Xi }^{}\pi ^+`$ yields $`(8\%,+4\%)`$ and consequently $`\alpha _\mathrm{\Xi }\alpha _\mathrm{\Lambda }=0.57\pm 0.21\text{(stat)}_{0.02}^{+0.05}\text{(syst)}`$.

From the measurement of $`\alpha _{\mathrm{\Xi }_c^0}\alpha _\mathrm{\Xi }`$ and the PDG evaluation: $`\alpha _\mathrm{\Xi }=0.456\pm 0.014`$ , we obtain $`\alpha _{\mathrm{\Xi }_c^0}=0.56\pm 0.39\text{(stat)}_{0.09}^{+0.10}\text{(syst)}`$. The physically allowed range of a decay asymmetry parameter is $`1<\alpha <+1`$. Our result prefers a negative value: $`\alpha _{\mathrm{\Xi }_c^0}`$ is $`<0.1`$ at the 90% CL. The central value is in the middle of the theoretically expected range but is not yet precise enough to choose between models.

We note that $`\alpha _\mathrm{\Xi }`$ has, so far, always been determined as a product of decay asymmetries using Equation 2. In principle, given sufficient statistics, the loglikelihood method in conjunction with the two-fold joint angular distribution of Equation 4 allows the direct measurement of all three asymmetry parameters: $`\alpha _{\mathrm{\Xi }_c^0}`$, $`\alpha _\mathrm{\Xi }`$, and $`\alpha _\mathrm{\Lambda }`$.

In conclusion, from a sample of $`138\pm 14`$ $`\mathrm{\Xi }_c^0\mathrm{\Xi }^{}\pi ^+`$ decays we have measured $`\alpha _{\mathrm{\Xi }_c^0}\alpha _\mathrm{\Xi }=0.26\pm 0.18\text{(stat)}_{0.04}^{+0.05}\text{(syst)}`$, from which we obtain $`\alpha _{\mathrm{\Xi }_c^0}=0.56\pm 0.39\text{(stat)}_{0.09}^{+0.10}\text{(syst)}`$. To the best of our knowledge this is the first measurement of a charm strange baryon decay asymmetry parameter.

We gratefully acknowledge the effort of the CESR staff in providing us with excellent luminosity and running conditions. This work was supported by the National Science Foundation, the U.S. Department of Energy, Research Corporation, the Natural Sciences and Engineering Research Council of Canada, the A.P. Sloan Foundation, the Swiss National Science Foundation, and the Alexander von Humboldt Stiftung.


text:
# Electromagnetic 𝐾⁺ production on the deuteron with hyperon recoil polarization

## 1 INTRODUCTION

Recent rigorous calculations of light hypernuclei  have contributed interesting insight into low-energy properties of the $`YN`$ interaction above the $`\mathrm{\Lambda }`$ threshold. However, no clear understanding of the $`YN`$ interaction has emerged around the $`\mathrm{\Sigma }`$ threshold. Electro- and photoproduction processes of $`K^+`$ on light nuclei offer a unique possibility for studying the $`YN`$ interaction in the continuum, especially near the $`\mathrm{\Sigma }`$ threshold. An inclusive $`d(e,e^{}K^+)YN`$ experiment has already been performed, and the data for $`d(\gamma ,K^+Y)N`$ and 3He$`(\gamma ,K^+Y)N`$ are being analyzed at TJNAF.

We have analyzed the inclusive $`d(\gamma ,K^+)`$ and exclusive $`d(\gamma ,K^+\mathrm{\Lambda }(\mathrm{\Sigma }))`$ processes , and report here preliminary results of the electroproduction process $`d(e,e^{}K^+)`$. This study aims to investigate the coupled $`\mathrm{\Lambda }N\mathrm{\Sigma }N`$ interaction in the final state and incorporates the modern $`YN`$ interactions of the Nijmegen group, NSC97f  and NSC89 which have been found to give a reasonable binding energy for the hypertriton. Kaon photoproduction on the deuteron is also important since it allows access to the elementary cross sections on the neutron, such as $`\gamma +nK^++\mathrm{\Sigma }^{}`$, in kinematic regions where final-state interaction effects are small.

## 2 PHOTOPRODUCTION

Numerical results of the inclusive $`d(\gamma ,K^+)`$ cross sections, using an updated production operator , are shown as a function of lab momentum $`P_K`$ in Fig. 1. For the details of our theoretical formulation, we refer the reader to Ref. . The incident photon energy is 1.3 GeV, while the outgoing kaon angle is fixed to 1 degree. The two pronounced peaks around $`P_K=945`$ and $`809`$ MeV/c are due to the quasifree scattering between photon and one of the nucleons. The results with the final state $`YN`$ interaction NSC97f are compared to the PWIA results. Sizable FSI effects are seen around both $`\mathrm{\Lambda }`$ and $`\mathrm{\Sigma }`$ thresholds, and to a lesser degree at the two quasifree peak positions.

For the same $`E_\gamma `$ and $`\theta _K`$, the exclusive $`d(\gamma ,K^+\mathrm{\Lambda })n`$ cross section and double polarization observable $`C_z`$ at $`P_K=870`$ MeV/c are shown in Figs. 2 and 3, respectively. Figures 4 and 5 depict these observables for $`d(\gamma ,K^+\mathrm{\Sigma }^{})p`$ at $`P_K=810`$ MeV/c. As indicated in Fig. 1, the former value of $`P_K`$ is close to the $`K^+\mathrm{\Sigma }N`$ threshold, while the latter one corresponds to the $`\mathrm{\Sigma }`$ quasifree peak position. While the values for $`C_z`$ in PWIA are almost 100%, the FSI results show dramatic deviations. Furthermore, the two $`YN`$ forces of NSC97f and NSC89 become clearly distinguishable for this observable. Experimentally, measuring this observable involves using circularly polarized photons along with detecting the recoil polarization of the hyperon in the final state.

## 3 ELECTROPRODUCTION

Here we present first preliminary results for the electroproduction process $`d(e,e^{}K^+)`$. Two sets of results are shown in Figs. 6 and 7. The incident electron energy $`E_e`$ and momentum transfer $`Q^2`$ in Fig. 6 is set to reproduce the conditions of the recent Hall C experiment at TJLAB . As in the case of the photoproduction $`d(\gamma ,K^+)`$, $`YN`$ FSI effects are seen near both $`\mathrm{\Lambda }`$ and $`\mathrm{\Sigma }`$ threshold. However, in Fig. 6, the FSI has effects in a wide range above the $`\mathrm{\Sigma }`$ threshold. Figure 7 shows a prominent enhancement around the $`\mathrm{\Sigma }`$ threshold which is not a simple threshold effect but is caused by a $`YN`$ $`t`$-matrix pole in the complex momentum plane .

## 4 OUTLOOK

We investigate cross sections and hyperon polarization for the $`K^+`$ photoproduction on the deuteron, and find large hyperon-nucleon FSI effects in the double polarization observable $`C_z`$. Also, in the electroproduction, for suitable $`Q^2`$ values, cross sections show a prominent enhancement around the $`\mathrm{\Sigma }`$ threshold. A systematic analysis for a wide range of kinematics for both photo- and electroproduction processes is in progress. Future studies will investigate final-state interaction effects in kaon photo- and electroproduction on the $`A=3`$ system.


text:
# Etching of random solids: hardening dynamics and self-organized fractality

## Introduction

Corrosion of solids is an everyday phenomenon with evident practical consequences and applications uhlig . Recent developments of theoretical tools for the study of disordered systems and fractals in the context of statistical mechanics, has triggered and outburst of activity in this subject.

When an etching solution is put in contact with a disordered etchable solid, the solution corrodes the weak parts of the solid surface while the hard parts remain un-corroded. During this process new regions of the solid (both hard and weak) come into contact with the etching solution. If the volume of the solution is finite and the etchant is consumed in the chemical reaction, then the etching power of the solution diminishes progressively and the corrosion rate slows down. When the solution is too weak to coorrode any part of the hardened solid surface, the dynamics spontaneously stops. Note that the etchant concentration at the arrest time is not zero, and, as shown below, its final value is strictly related to the percolation critical point. One of the most interesting aspects of this phenomenon is that the final connected solid surface has a fractal geometry up to a certain characteristic scale $`\sigma `$ given by its thickness. This is precisely the qualitative phenomenology that has been recently observed in a nice experiment on pit corrosion of aluminum thin films Balazs . As we show below this fractality falls in the universality class of percolation theory Aharony .

## The model

Recently, a simple dynamical model of etching, capturing the aforementioned phenomenology, has been proposed and studied model ; GBS .

The model is reviewed in Fig. 1 and described as follows:
(1) The solid is represented by a lattice of sites exhibiting random “resistances to corrosion” $`r_i[0,1]`$ uniformly distributed. It has a width $`L`$ and a given fixed depth $`Y`$ (in the case infinite). The etching solution has a finite constant volume $`V`$ ($`VL`$) and contains an initial number $`N_{et}(0)`$ of etchant molecules. Therefore the initial etchant concentration is given by $`C(0)=N_{et}(0)/V`$. At any time-step $`t`$ the “etching power” $`p(t)`$ of the solution is proportional to the etchant concentration $`C(t)`$. Hereafter we take the proportionality constant equal to one without loss of generality, and we choose $`p_c<p(0)1`$, where $`p_c`$ is the percolation threshold of the lattice.
(2) The solution is initially in contact with the solid through the bottom boundary $`y=0`$. At each time-step $`t`$, all surface sites with $`r_i<p(t)`$ are dissolved, and a particle of etchant is consumed for each corroded site. Hence, the concentration of the solution decreases with $`t`$. Moreover, depending on the connectivity criterion chosen for the lattice, $`m(t)`$ new solid sites, previously in the bulk, come into contact with the solution. Note that at the next time-step, only these last sites can be corroded by the solution, as the other surface sites have already resisted to etching and $`p(t)`$ decreases with $`t`$. Calling $`n(t)`$ the number of dissolved solid sites at time $`t`$, one can write

$$p(t+1)=p(t)\frac{n(t)}{V}.$$
(1)
The corrosion dynamics spontaneously stops at the time-step $`t_f`$, when all the surface sites have a resistance $`r>p(t_f)`$.

The main characteristics of the model are:
(i) The final value $`p_f=p(t_f)`$ is slightly smaller than the percolation threshold $`p_c`$. The difference $`|p_fp_c|0`$ as $`(L/V)^{\alpha _p}`$ with $`\alpha _p0.45`$, when the limits $`L+\mathrm{}`$ and $`V/L+\mathrm{}`$ are taken in the right way GBS . This fact can be explained in terms of percolation theory. Indeed, this theory implies that, for $`p(t)p_c`$, the probability to have a connected path of solid sites each one with $`r>p(t)`$, then stopping corrosion, is equal to one in the right thermodynamic limit.

(ii) The global temporal behavior of $`p(t)`$ can be divided (see Fig. 2-a) into two regimes: a smooth regime when $`p(t)>p_c`$, and a critical regime when $`p(t)p_c`$. The duration of the smooth regime is measured approximately by $`t_c`$, defined by $`p(t_c)=p_c`$. The quantity $`t_c`$ is found to scale with the ratio $`L/V`$ in the following way $`t_cV/L`$. The duration of the critical regime is given by $`(t_ft_c)`$, and scales as $`(t_ft_c)(V/L)^{\alpha _{t_f}}`$ with $`\alpha _{t_f}0.55`$ (apart from a further linear dependence on $`\mathrm{log}L`$ of the scaling coefficient due to the extremal nature of $`t_f`$ frac ).
(iii) The global surface of the solid in contact with the solution displays a peculiar roughening dynamics (three snapshots at almost initial, intermediate, and final time are reported in Fig. 2-b). First of all, during the process, the solution can surround and detach finite solid islands. Consequently, the global solid surface is composed by the sum of the perimeter of the “infinite” solid, here called corrosion front, and the surfaces of the finite islands. In the first smooth regime, the corrosion is well directed and the corrosion front becomes progressively rougher and rougher, while finite detached islands are quite small. In the second critical regime, the correlations revealed by the hardening process become important: the dynamics becomes locally isotropic generating a fractal front. Detached islands in this last period are quite big with a fractal perimeter too.

(iv) The final corrosion front is fractal with a fractal dimension $`D_f=1.753\pm 0.05`$ (see Fig. 3-a) up to a characteristic scale $`\sigma `$ which is the final thickness of the front. Beyond this length the corrosion front acquires a one-dimensional shape. Also $`\sigma `$ is found to scale with the ratio $`L/V`$ in the following way: $`\sigma (L/V)^{\alpha _\sigma }`$ with $`\alpha _\sigma 0.571/D_f`$. As shown below, all these properties can be explained in the framework of percolation theory and in particular of Gradient Percolation GP .
(v) Another important feature captured by this simple model is the progressive hardening of the solid surface, due to the corrosion only of sites with a resistance weaker than the etching power of the solution, which, in its turn, decreases with time. This hardening effect is described quantitatively by the behavior of the normalized resistance histogram $`\varphi (r,t)`$ of the global surface sites. Obviously, $`\varphi (r,0)=1`$ with $`r[0,1]`$, while $`\varphi (r,t_f)`$ is given in Fig. 3-b. As one can see, the final global surface is much more resistant to corrosion than the initial one.

## Analytical result and percolation

In order to have a mathematical insight in the previous numerical results, and a more clear explanation in terms of percolation, it is useful to study the temporal evolution of the above defined resistance histogram. Let be $`h(r,t)dr`$ the number of surface sites with resistance in the interval $`[r,r+dr]`$ at time $`t`$ ($`\varphi (r,t)`$ is the normalized version of $`h(r,t)`$). Obviously, the number of corroded sites at time $`t`$ is $`n(t)=_0^{p(t)}𝑑rh(r,t)`$. Let us define $`\omega (t)=m(t)/n(t)`$. Then, using Eq. (1), we can write the two coupled equations:

$`h(r,t+1)=h(r,t)h(r,t)\theta [p(t)r]+\omega (t){\displaystyle _0^{p(t)}}𝑑r^{}h(r^{},t)`$
$`p(t+1)=p(t){\displaystyle \frac{_0^{p(t)}𝑑r^{}h(r^{},t)}{V}},`$ (2)
where $`\theta (x)`$ is the usual step Heavyside function. In the first relation of Eq. (2), the second term on the right side represents the number of just corroded sites with resistance between $`r`$ and $`r+dr`$, while the third one gives the contribution to the new histogram coming from the $`m(t)`$ new surface sites. All the geometrical features of the system related to the behavior of $`p(t)`$ and $`h(r,t)`$ are described by the coefficient $`\omega (t)`$. In order to close this system of equations, we introduce two different approximations of $`\omega (t)`$ for the smooth and the critical regimes. One can show GBS  that, since in the smooth regime $`n(t),m(t)1`$ and fluctuations are small, the right approximation is $`\omega (t)=1/p(t)`$. Using this relation in Eq. (2), one gets

$$p(t)=p(0)\mathrm{exp}(t/\tau ),$$
(3)
with $`\tau =1/\mathrm{ln}(1L/V)V/L`$. The good accuracy of this approximation in the smooth regime is shown by Fig. 2-a. Eq. (3) gives $`t_cV/L`$ in good agreement with simulations. In the critical regime, the only simple approximations is a mean-field-like relation $`\omega (t)=1/p_c`$ GBS , forgetting fluctuations which are important in this regime. Putting this relation inside Eq. (2) one obtains, $`(p_cp_f)\sqrt{L/V}`$ and $`(t_ft_c)\sqrt{V/L}`$ which are not very far from the simulation behaviors. Finally, the histogram $`h(r,t)`$ obtained by Eq. (2) within these approximations is in excellent agreement with simulations (see Fig. 3-b).

The link between this dynamical etching model and percolation criticality can be shown by analyzing more deeply Eq. (3) valid in the smooth regime. From Eq. (3), one has $`t_cV/L\mathrm{ln}[p(0)/p_c]`$. This implies that the condition for smooth regime is $`t<V/L`$. Then in this temporal range one can linearize Eq. (3):

$$p(t)p(0)(1tL/V).$$

Since in this regime the solution advances mainly layer by layer into the solid, one can say that the layer etched at time $`t`$ is at a depth $`y=t`$. We can map then our dynamical model in a static model defined by the following rule: a general solid site belonging to the layer $`y`$ is replaced by the solution with a $`y`$-dependent probability $`p(y)p(0)(1yL/V)`$. This static model, if extended to any value of $`p(y)>0`$, is known under the name of Gradient Percolation (GP) GP . In our case the gradient $`p(y)`$ of replacement (etching) probability satisfies $`pL/V`$. Why this mapping can be extended also out of the smooth regime is more clear from field theory arguments Munoz  . The GP model is characterized by the following features: the boundary (the corrosion front in our case) between the connected etched region and the connected non-etched one has fractal features in the universality class of critical percolation. In particular it has a fractal dimension $`D_f=7/4`$ coinciding with the hull fractal dimension of infinite percolation cluster at criticality Duplantier . The fractality is not extended to the whole corrosion front, but it is restricted to regions of size $`\sigma p^{1/D_f}`$, which, considering $`pL/V`$, is in good agreement with our numerical results. Furthermore, the range of $`p`$’s on this boundary is given by $`\mathrm{\Delta }p\sigma p(L/V)^{11/D_f}`$, which is well confirmed by simulation considering the quantity $`p_cp_f`$. Note that from previous scaling relations, one can write $`\sigma \mathrm{\Delta }p^\nu `$, where $`\nu =1/(D_f1)=4/3`$ is the usual correlation length exponent in percolation. This shows that in the dynamical etching model $`\sigma `$ can be interpreted as the correlation length of the corrosion front: that is, the fact that $`p_f`$ is not exactly equal to $`p_c`$ introduce a finite correlation length in the problem.

This work has been supported by the European Community TMR Network “Fractal Structures and Self-Organization” ERBFMRXCT980183.


text:
# Simulations of deep pencil-beam redshift surveys

## 1 Introduction

The redshift distribution of galaxies in the pencil-beam survey of Broadhurst et al. (1990, hereafter BEKS) displayed a striking periodicity on a scale of 128$`h^1`$Mpc. This result has attracted a good deal of interest over the subsequent decade, and the significance and nature of periodicity or regularity in the distribution of galaxies has remained the subject of a stimulating debate in both observational and theoretical cosmology. Although a number of studies have been devoted to the BEKS pencil-beam survey and other similar surveys, several fundamental questions remain unanswered.

From the theoretical viewpoint, it is important to decide whether such apparently periodic galaxy distributions can occur with reasonable probability in a Cold Dark Matter universe, or require physics beyond the standard paradigm. Performing large simulations can directly address this question. The first simulation specifically designed for pencil-beam comparisons was that of Park and Gott (1991, hereafter PG). Their rod-shaped CDM simulation allowed them to create twelve quasi-independent mock pencil-beam surveys similar in length to that of BEKS. One of their samples appeared ‘more periodic’ than the BEKS data according to the particular statistical test they used for comparison. Other authors (Kurki-Suonio et al.; Pierre 1990; Coles 1990; van de Weygaert 1991; SubbaRao & Szalay 1992) have used purely geometrical models such as cubic lattices and Voronoi foams to explore the implications of apparent regularities similar to those found by BEKS. In particular, SubbaRao and Szalay (1992) presented a sequence of Monte Carlo simulations of surveys of Voronoi foams, showing that such a model can successfully reproduce the data as judged by a variety of statistical measures, for example, the heights, positions and signal-to-noise ratios of the highest peaks in the power spectra. Kaiser & Peacock (1991) argued that the highest such peaks in the BEKS data are not sufficiently significant to be unexpected in a CDM universe, but did not support this conclusion with detailed simulations. Dekel et al. (1992) introduced other statistics, more similar to those of PG, and again concluded that the apparent periodicity seen in the real data is not particularly unlikely in any of the toy models they used for comparison. Their models include Gaussian models with an extreme initial power spectrum with power only on scales $`100h^1`$Mpc. They found regular ‘galaxy’ distributions a few per cent of the time and concluded that the BEKS data do not rule out all Gaussian models. However, these theoretical studies did not give any clear answer to the question posed above: is the BEKS regularity compatible with the standard CDM paradigm? We attempt to answer this below using versions of all the statistical tools developed in earlier papers.

There have been several interesting observational developments after BEKS. Willmer et al. (1994) found that, if the original BEKS deep survey at the North Galactic Pole had been carried out 1 degree or more to the west, many of the peaks would have been missed. On the other hand, Koo et al. (1993) added new data from a wider survey to the original BEKS data and found the highest peak in the power spectrum to be further enhanced. They also analysed another set of deep pencil-beam surveys and found a peak of weaker significance on the same scale, 128 $`h^1`$Mpc. This raises another question: is 128 $`h^1`$Mpc a preferred length scale for the galaxy distribution? Further support for such a preferred scale has been presented by Tully et al. (1992), Ettori et al. (1997) and Einasto et al. (1997). Thus one can wonder whether a single scale could be indicated with such apparent consistency within the CDM paradigm.

With the important exception of the work of PG there has been surprisingly little comparison of the BEKS data with direct simulations of standard CDM cosmogonies. Even before the BEKS discovery, White et al. (1987) had shown that pencil-beams drilled through periodic replications of their CDM simulations frequently showed a kind of ‘picket fence’ regularity in their redshift distribution. Frenk (1991) confirmed this result and concluded that regular patterns similar to that seen in the BEKS data are easy to find in their simulations. However, it is clearly dangerous to make use of periodic replications of a simulation when assessing the significance of apparent periodicities in the redshift distribution. It is preferable to simulate a volume large enough to encompass the whole survey. Furthermore, since many independent artificial surveys are needed to establish that the real data are highly unlikely in the cosmogony simulated, the simulated volume must be fully three-dimensional (unlike that of PG) to allow the creation of many quasi-independent lines-of-sight. A final consideration is that the BEKS data reach to redshifts beyond 0.3, so that evolution of clustering along the survey may not be negligible.

In this paper we investigate the distribution of ‘galaxies’ along the past light-cones of hypothetical observers. Particle positions and velocities on these light-cones were generated as output from the Hubble Volume Simulations (Evrard et al. 2000). These very large CDM $`N`$-body simulations were recently performed by the Virgo consortium and each used $`10^9`$ particles to follow the evolution of the matter distribution within cubic regions of an $`\mathrm{\Omega }=1`$ $`\tau `$CDM ($`\mathrm{\Omega }=0.3`$ $`\mathrm{\Lambda }`$CDM) universe of side 2000 $`h^1`$Mpc (3000 $`h^1`$Mpc). Such large volumes allow many independent light-cones to be generated out to $`z1`$. The light-cone output automatically accounts for clustering evolution with redshift. The principal uncertainty lies in how to create a ‘galaxy’ distribution from the simulated mass distribution. We employ Lagrangian bias schemes similar to those of White et al. (1987) and Cole et al. (1998). Individual particles are tagged as galaxies with a probability which depends only on the smoothed initial overdensity field in their neighbourhood. The parameters of these schemes are adjusted so that the present-day correlations of the simulated galaxies match observation. Many quasi-independent mock pencil-beam surveys can then be created adopting the geometry and the galaxy selection probability with distance of the BEKS surveys.

Our discussions focus primarily on the significance of the BEKS data in comparison with our CDM samples. We begin by following the methods used originally by BEKS, namely, redshift counts, pair separation distributions, and the one-dimensional power spectrum. Redshift counts are translated into a distribution in physical distance assuming the same cosmological parameters as BEKS. For the one-dimensional power spectra, the height of the highest peak is the most important statistic. Szalay et al. (1991) show that the statistical significance of the highest peak of the BEKS data is at $`10^4`$ level, based on the formal probability for the peak height. This calculation was disputed by Kaiser& Peacock (1991) because of the difficulty in estimating the appropriate noise level. We calculate relative peak-to-noise ratios of the highest peaks in the power spectra in identical ways for real and simulated data and so can compare the two without needing to resolve this issue. We also apply two additional statistical tests for regularity, the $`\mathrm{\Delta }`$ test of PG and a ‘supercluster’ statistic designed by Dekel et al. (1992).

Our paper is organised as follows. In section 2 we present details of the $`N`$-body simulations from which our pencil-beam samples are drawn. In section 3 we explain our bias scheme. In section 4 we describe our mock pencil-beam surveys which mimic as closely as possible the actual observations of BEKS. Our main results for power spectra are given in section 5. Results are given in section 6 for the $`\mathrm{\Delta }`$ test, and in section 7 for the ‘supercluster’ statistic. We present our conclusions in section 8.

## 2 $`N`$-body simulation

The simulation data we use are the so-called “light-cone outputs” produced from the Hubble Volume simulations (details are in Evrard et al. 2001). The basic simulation parameters are tabulated in Table 1, where $`L_{\mathrm{box}}`$ is the box size in $`h^1`$Mpc, $`\mathrm{\Gamma }`$ stands for the shape parameter of the initial power spectrum and $`m_\mathrm{p}`$ is the mass per particle; other notations are standard.

The light-cone outputs are created in the following way. We define an observer at a specific point in the simulation box at the final time. The position and velocity of each particle is recorded whenever it crosses the past light-cone of this observer, and these phase-space coordinates are accumulated in a single data file. The evolution of clustering with lookback time (distance from the observer) is automatically included in such data. As we require mock pencil-beam surveys which reach $`z`$0.5 (spanning $``$1000$`h^1`$Mpc in physical scale), such light-cone output is both realistic and desirable. We use stored data from two different light-cone outputs for each cosmology. Each covers one octant of a sphere, and they emanate in opposite directions from the same point. Figure 1 illustrates this geometry. In each case we use data out to a comoving distance of 1500$`h^1`$Mpc, corresponding to redshift 0.77 in the $`\tau `$CDM model and $`0.58`$ in the $`\mathrm{\Lambda }`$CDM model. For the $`\tau `$CDM case the total length covered is larger than the side of the simulation box, but this has a negligible effect on the mock BEKS surveys we construct.

## 3 Galaxy selection

To create realistic mock surveys we have to select particles as galaxies with the same distribution in depth as the real data and with an appropriately ‘biased’ distribution relative to the dark matter. We do this in two stages. First we identify a biased subset of the particles chosen according to the value of the smoothed linear mass overdensity at their position at high redshift. The parameters defining this identification are chosen so that the two-point correlation function of the identified ‘galaxies’ at $`z=0`$ matches the observed correlation function of low redshift galaxies. For the $`\tau `$CDM model, we are able to achieve this while retaining about two-thirds of the simulation particles as ‘galaxies’. The resulting comoving ‘galaxy’ number density is 0.08 $`h^3`$Mpc-3. For the $`\mathrm{\Lambda }`$CDM model we get a number density of ‘galaxies’ in the range 0.02 to 0.033 $`h^3`$Mpc-3 depending on the bias scheme. This lower number density is due to the low number density of the dark-matter particles in this model. The second stage is to mimic the effect of the apparent magnitude limits of the real galaxy surveys by including ‘galaxy’ particles into the final mock catalogues with a probability which depends on distance from the observer. Since this stage is independent of the first, we are effectively assuming that the clustering of galaxies is independent of their luminosity. Our radial selection function is based on those directly estimated for the BEKS surveys.

### 3.1 Lagrangian Bias

Cole et al. (1998) developed and tested a set of bias schemes to extract ‘galaxies’ from $`N`$-body simulations. The procedure we use for the first stage of our galaxy selection is similar to their Model 1, but has a different functional form for the probability function. Since we need a bias factor greater than unity for the $`\tau `$CDM model and less than unity for the $`\mathrm{\Lambda }`$CDM model (a result of the differing mass correlations in the two cases) ‘galaxies’ need to avoid regions of low initial density in the $`\tau `$CDM simulation and to avoid regions of high initial density in the $`\mathrm{\Lambda }`$CDM simulation. We begin by smoothing the density field at an early time with a Gaussian, $`\mathrm{exp}(r^2/2r_o^2)`$ with $`r_o=3h^1`$Mpc and assigning the corresponding overdensity $`\delta `$ to each dark-matter particle. Then a normalised overdensity $`\nu =\delta /\sigma _s`$ is computed, where $`\sigma _s`$ is the root mean square value of particle $`\delta `$-values. Finally, we define a probability function $`P(\nu )`$ which determines whether a particle is tagged as a ‘galaxy’. We random-sample dark-matter particles for tagging as galaxies based on this probability. Once tagged as a galaxy in this way, particles remain tagged throughout the simulation, and so become potentially visible in our mock surveys whenever their world-line crosses a light-cone.

For all the bias models described below the ‘galaxy-galaxy’ correlation function was calculated in real space within a cubic box of side 200 $`h^1`$Mpc with the observer at one corner. These correlation functions are shown in Figure 2.

$`\tau `$CDM model bias t1: For the $`\tau `$CDM model, we chose a simple power law form $`P(\nu )(\nu \nu _c)^{0.2}`$ for the probability function. We impose a threshold at $`\nu =\nu _c=0.55`$ below which the probability is set to be zero. This suppresses the formation of ‘galaxies’ in voids. These parameters were determined by matching the present-day two-point correlation function of the ‘galaxies’ to the observational result for the APM survey (Baugh 1996, see Figure 2) on length scales from 0.2 $`h^1`$Mpc to 20$`h^1`$Mpc. We note here that in our $`N`$-body simulations the gravitational softening length is 0.1$`h^1`$Mpc.

$`\mathrm{\Lambda }`$CDM bias model L1: For the $`\mathrm{\Lambda }`$CDM model, we must ‘anti-bias’ because the predicted mass correlations on small scales are substantially larger than observed galaxy correlations (see, for example, Jenkins et al. 1998). We set a sharp upper cut-off at $`\nu _c=1.34`$, above which $`P(\nu )`$ is zero. All particles below this threshold are equally likely to be ‘galaxies’ ($`P=`$ const). Although this may seem unphysical, more realistic modelling of galaxy formation in $`\mathrm{\Lambda }`$CDM models does indeed produce the anti-bias required for consistency with observation, albeit through a more complex interplay of statistical factors (Kauffmann et al. 1999; Benson et al. 2000). We use a simpler scheme in order to produce the desired two-point correlation function; on scales of interest here, only a small anti-bias is necessary.

$`\mathrm{\Lambda }`$CDM bias model L2 : For comparison purposes, we applied a second bias model to the $`\mathrm{\Lambda }`$CDM simulation. We set an additional lower threshold at $`\nu _{low}=0.7`$ below which we again set $`P=0`$. Thus the probability takes a non-zero (and constant) value only in the range $`\nu _{low}\nu \nu _{high}`$, where now $`\nu _{high}=0.9`$. This model fits the observed correlations of galaxies just as well as L1 but enhances the emptiness of voids

Figure 3 illustrates bias effects by comparing the distribution of dark-matter particles and of ‘galaxies’ in a thin slice through part of the simulation box at $`z=0`$. For model L1 the effect is difficult to detect visually, whereas the effect of the lower cut-off in model L2 is obvious. Similarly, for the $`\tau `$CDM model, underdense regions(voids) are clearly accentuated in the ‘galaxy’ distribution relative to the dark matter distribution. Cole et al. (1998) show similar plots to demonstrate how strong bias in high density model universes maps underdense regions in the mass distribution onto voids in the galaxy distribution. Such contrasted voids are generally seen in strongly biased models regardless of the functional form of $`P(\nu )`$.

We can quantitatively study the bias in our models by measuring the nonlinear biasing parameters introduced by Dekel & Lahav (1999) (see also Sigad et al. 2000; Somerville et al. 2000). We compute the slope $`\widehat{b}`$ and nonlinearity $`\stackrel{~}{b}`$ following the procedure described in Somerville et al. (2000). Figure 4 show the biasing relation between the ‘galaxy’ density field $`\delta _g`$ and the dark matter density field $`\delta _d`$, smoothed with a 8$`h^1`$Mpc scale top-hat filter. For each of our bias models, the mean biasing function $`b(\delta _d)`$, and its moments

$$\widehat{b}=\frac{}{\sigma ^2},\stackrel{~}{b}^2=\frac{}{\sigma ^2},$$
(1)
are given in Figure 4. In the above expression we have used $`\sigma ^2=<\delta _d^2>`$ for the standard deviation. Strong biasing in t1 and anti-biasing (for $`1+\delta _d>0`$) in L1 and L2 are clearly seen, and reflected in the values for the effective slope; $`\widehat{b}=1.44`$ for t1, $`0.84`$ for L1 and $`0.90`$ for L2.

## 4 Survey Strategy

### 4.1 Geometry

We construct artificial surveys with a geometry very similar to that of the data analysed by BEKS. This consists of four surveys – a deep and a shallow survey near each Galactic Pole. The northern deep survey lies within a cone of 40-arcmin diameter about the pole and is made up of a set of roughly circular patches each 5-arcmin in diameter. About 10 small patches were surveyed but not all were completed by the time of writing so that the exact number of patches used in BEKS is unclear. For our artificial surveys we choose 9 circular patches within the 40-arcmin diameter cone, each of diameter 5 arcminutes. We place these irregularly and ensure no overlaps between them. For model L2 the number density of ‘galaxies’ is too small to match the observations, so we increased the diameter of our patches to 7-arcmin. Although this widening results in a slight increase in the effective survey volume, the small patches still lie well within the larger cone of 40-arcmin diameter. The volume increase is compensated in the radial selection we decribe below, so that the resulting ‘galaxy’ distribution is consistent with the desired distribution given in BEKS. For the deep-narrow pencil-beams, the transverse length scale is much smaller (the cone diameter is $`4h^1`$Mpc at z=0.2 where the radial selection function takes its maximum value) than the $`100h^1`$Mpc scale we address, so the increase in the patch width does not affect our results. In all cases the redshift counts in all patches were binned together to create a single deep survey. The northern shallow survey has a simpler geometry. A square area of about 14 square degrees is selected near the Galactic Pole, but with its centre offset by 7 degrees. The magnitude limit of the shallow survey is about 5 magnitudes brighter than that of the deep survey.

Towards the South Galactic Pole both surveys are centred very close to the pole itself. The deep survey is confined within a cone of 20-arcmin diameter,while the shallow survey covers an area of 14 square degrees and has a magnitude limit about 4 magnitudes brighter.

When making an artificial survey we choose a random direction in the simulation as the Galactic polar axis and then define all areas on the artificial sky with reference to this direction. The light-cone outputs from our Hubble Volume simulations cover enough ‘sky’ to allow us to make well over 1000 near-independent artificial surveys.

### 4.2 Radial selection

‘Galaxies’ projected in our survey regions are assigned weights for selection depending on their distances. We use the estimated galaxy distributions given in Figure 1 of BEKS to define the relevant selection functions for each survey. The data are read off in redshift bins of width $`\mathrm{\Delta }z=0.005`$ for the deep surveys and $`\mathrm{\Delta }z=0.001`$ for the shallow surveys. We then derive a smoothed model galaxy distribution $`dN/dz`$ for each survey and compute the corresponding comoving number densities from the number counts and the volume elements given by the survey geometry and the assumed cosmology. As explained above, we had to increase the size of the patches in the northern deep survey in case L2 in order to get the correct mean counts. This is easily accounted for by appropriate renormalisation. These radial selection functions are used as sampling probabilities to determine whether a particular ‘galaxy’ is included in a catalogue or not. We normalise our probabilities by matching the mean number of ‘galaxies’ in each survey to the number of galaxies in BEKS data. This matching is done for the 4 surveys independently. For consistency, the normalisation coefficients obtained are then kept constant when constructing all realisations for a particular model.

### 4.3 Peculiar velocities

The peculiar velocities of ‘galaxies’ must be taken into account to create realistic mock redshift surveys. We simply assign our ‘galaxies’ the peculiar velocities of their corresponding dark-matter particles. Thus, while the spatial distribution of ‘galaxies’ is biased, there is no additional bias associated with their peculiar velocities. On small scales peculiar velocities lead to ‘finger-of-God’ effects which suppress power in the apparent spatial distribution at high wavenumber. In our mock catalogues the root mean square values of the ‘galaxy’ line-of-sight peculiar velocities are 342 km/sec in t1, 358 km/sec for L1 and L2. The redshift bin width shown in BEKS is $`\mathrm{\Delta }z`$=0.005 for the deep surveys, which translates $``$1500 km/sec in recession velocity. Therefore, the assigned peculiar velocities of ‘galaxies’ do not smear out the true width of clumps in one-dimensional distributions, while they reflect properly the underlying velocity field. At the intermediate and small wavenumbers corresponding to the linear and quasi-linear regime, they increase the apparent power (e.g. Kaiser and Peacock 1991). These line-of-sight distortions reflect the enhanced contrast produced by infall onto superclusters. It is thus important to include the peculiar velocities when comparing simulations to the structures seen in the BEKS data.

## 5 Galaxy distribution and power spectrum analysis

The geometry of our light-cone datasets allows the axis of our artificial surveys to lie anywhere within one octant of the ‘sky’ (see Figure 1). We construct an ensemble of mock surveys with axes distributed uniformly across this octant in such a way that the areas covered by the corresponding deep surveys do not overlap. We end up with 1936 quasi-independent deep surveys for our $`\tau `$CDM model. For the $`\mathrm{\Lambda }`$CDM case an additional pair of light-cone outputs were stored, allowing us to construct 2625 disjoint deep surveys. To these deep pencil beams we add shallow surveys, whose volumes then have slight overlaps with those of neighbouring surveys. In practice, however, rather few ‘galaxies’ appear in more than one of our mock catalogues.

A series of plots of the redshift distribution and derived statistics are given for selected ‘mock BEKS surveys’ in Figure 6, which consists of 6 sets of 3 figures. These can be compared with Figure 5, which is actually for the real BEKS data, which we reproduce here for comparison with our simulation results. We read these data from Figure 2 in Szalay et al. (1991) where they are given as a histogram of bin width 10 $`h^1`$Mpc; when necessary for our analysis, we assume that the galaxies in each bin are uniformly distributed across the bin. The particular mock surveys in the following 6 plots were chosen to illustrate a variety of points made in the following sections.

### 5.1 One-dimensional distribution

In each set of plots in Figure 6, the top panel shows the distance histogram of ‘galaxies’ in the combined deep and shallow surveys. The total number of galaxies in these combined surveys is given in this panel. We have assumed an Einstein-de Sitter universe for both of our models when converting redshift to physical distance, although the actual value of $`\mathrm{\Omega }`$ is 0.3 in the $`\mathrm{\Lambda }`$CDM case. This apparent inconsistency is needed to allow a direct comparison with the analysis in BEKS where $`\mathrm{\Omega }=1`$ was also assumed. Szalay et al. (1991) noted that using low values of $`\mathrm{\Omega }`$ to convert redshift to distance reduces the significance of the apparent periodicity in the BEKS data. Throughout this paper we assume $`\mathrm{\Omega }=1`$ for this conversion.

In Figure 6, panel (a) shows one of the best catalogues in our t1 ensemble in that it gives the impression that ‘galaxies’ are distributed periodically and, in addition, the 1-D power spectrum shows a sharp and high peak. Panel (b) shows the same features but with a smaller characteristic spacing. In each of these plots we mark the best periodic representation of the data in the same way as BEKS. We determine the characteristic spacing from the position of the highest peak in the power spectrum, and we adjust the phase to match the positions of as many big clumps as possible. The characteristic spacing is indicated by the vertical dashed lines in the top panels in Figure 6. Panel (c) shows a good example whose power spectrum has a very high peak while the actual distance distribution does not show a periodic feature (discussed in section 5.3). Panels (d) and (e) show the best examples from our model L1 and L2, respectively, which show a good visual impression that ‘galaxies’ are spaced regularly. Finally panel (f) shows an example from model L2, which has a large characteristic length scale of $`200h^1`$Mpc.

### 5.2 Pairwise separation histograms

From the apparent distance distributions of the ‘galaxies’, it is easy to produce histograms of pairwise distance differences which can be used to search for characteristic scales in the structure within our pencil-beam surveys. Such pair counts are shown in the middle plot of each panel in Figure 6. These counts typically display a series of peaks and valleys which are particularly prominent in panels (a), (e) and (f), and, as noted by BEKS themselves, in the original BEKS data. Notice that these peaks appear regularly spaced as indicated by the dashed lines in these panels. The contrast between peaks and valleys can be used as a measure of the strength of the regularity. For the BEKS data the height difference between the first peak and the first valley is about a factor of 3, while the corresponding numbers are 2.4, 2.2 and 3.4 in panels (a), (e), and (f), respectively. Many of our artificial samples show a more complex behaviour, however. In panel (b) there is a deep valley at 150$`h^1`$Mpc followed by a high peak at 200$`h^1`$Mpc; the contrast is a factor of 5.3 despite this uneven spacing. A robust and intuitive definition of contrast is difficult to find. An alternation of small-scale peaks and valleys can coexist with apparently significant larger scale variations as is clearly seen in panel (c). If we focus specifically on the strongest peaks and valleys, their ratio, and indeed even their identification can depend on the specific binning chosen for the histograms. Because of these ambiguities we will not use these distributions further for quantitative analysis in this paper.

### 5.3 Power spectra

In order to compare our results directly with BEKS we calculate one-dimensional power spectra for our samples using the method described in Szalay et al. (1991). Each galaxy is represented by a Dirac delta-function at the distance inferred from its redshift (including its peculiar velocity). The power in each Fourier component is then

$$f_k=\frac{1}{N_g}\underset{n}{}\mathrm{exp}(2\pi ikr_n),$$
(2)
$$P_k=|f_k|^2$$
(3)
where $`N_g`$ is the total number of galaxies in the sample. The power spectra calculated in this way for our various samples are shown in the bottom plots of each panel in Figure 6. Our units are such that the wavelength corresponding to wavenumber $`k`$ is 1000/$`k`$ $`h^1`$Mpc. In panels (a) and (b) visual impression of the separation of clumps is consistent with the wavelength inferred from the power spectra. The highest peak is at $`k`$=8.0 in panel (a) and $`k=`$16.0 in panel (b), giving wavelengths of 125$`h^1`$Mpc and 62.5$`h^1`$Mpc, respectively. As we shall see, there is no unique length scale inferred from the power spectra.

If a pencil-beam penetrates a rich cluster, an interesting feature can arise. For example, in panel (c) there is a single large cluster at 600 $`h^1`$Mpc. Together with a few other clumps of moderate size, it produces a very high peak in the power spectrum without the distance distribution as a whole giving a visual impression of regularity (c.f. the top panel of panel (c)). Many of our samples in both the $`\tau `$CDM and $`\mathrm{\Lambda }`$CDM ensembles show high peaks in the power spectra with no apparent periodicity. Thus a very high peak in the power spectrum, particularly at low wave-number, is a poor indicator of the kind of regularity which is so striking in the original BEKS data. Interestingly, as Bahcall (1991) pointed out, if one of the BEKS survey beams passed near the centre of a rich cluster, the galaxy count in the corresponding distance bin would have been much larger than the maximum of 22 seen in the actual BEKS data (see also Willmer et al. 1994). (For comparison, the maximum bin count in the histogram of panel (a) is 23.) We note that the Poisson sampling noise in each power spectrum can be estimated as $`1/N_g`$. As a result, a big clump raises the statistical significance of ‘structure’ both by enhancing the strength of peaks and by lowering the estimated noise.

In order to compare samples with different total numbers of ‘galaxies’, we calculate the signal-to-noise ratio of the highest peak in the power spectrum following the procedure of Szalay et al. (1991). We define the peak-to-noise ratio of a sample as, $`X`$=(peak height)/(noise level) where the noise is estimated from the sum in quadrature of the Poisson sampling noise and the clustering noise,

$$f_0=\frac{1}{N_g}+\frac{\xi _0}{M},$$
(4)
where, as before, $`N_g`$ is the total number of galaxies in the sample, $`\xi _0`$ is the small-scale two point correlation function averaged over a cell of depth 30$`h^1`$Mpc, and $`M`$ is the number of cells along the survey axis. We use the value $`\xi _0/M=1/80`$ as in Szalay et al. (1991). Although Szalay et al. derived this formula from a simple model with cylindrical geometry, they showed that this estimator agrees well with another internal estimator calculated from the cumulative distribution of power. To facilitate direct comparison with the earlier work we also use equation (3) to compute signal-to-noise ratios for the highest peaks in our samples. These S/N ratios are given in each of the power spectrum plots in Figure 6. Figure 7 shows both the differential and the cumulative distribution of peak-to-noise ratio in our mock surveys. The principal difference between the $`\tau `$CDM and $`\mathrm{\Lambda }`$CDM ensembles lies in the position of the peak in the differential count. For both L1 and L2 the peak is at smaller $`X`$ than in t1. This difference can be traced to the value of $`\mathrm{\Omega }`$ we assume for analysis. Adopting $`\mathrm{\Omega }=1`$ for converting redshift to physical distance causes the value of $`M`$, the number of cells of width 30$`h^1`$Mpc along the survey axis, to be underestimated for $`\mathrm{\Lambda }`$CDM. Using the noise estimator (3) with this value of $`M`$ then overestimates the noise levels for L1 and L2 (see SubbaRao and Szalay (1992) for discussion of a similar point).

We plot in Figure 8 the wavenumber distribution of the peaks whose S/N ratios are higher than that of the original BEKS data ($`X`$=11.8). We find 14 samples satisfy this condition in t1, 7 in L1 and 13 in L2. The distributions of the peaks with $`X>8.0`$ (the ‘tails’ of the number count in Figure 7) are also shown in Figure 8. By checking the distance distributions we have found that highly significant peaks at $`k5`$ are almost always due to one or two strong clumps, as noted above. Very few catalogues give a high peak on scales similar to the BEKS data. It is interesting that the frequency of such catalogues is significantly higher in t1 than in L1 and L2. The difference is primarily due to the number density of rich clusters over the redshift range surveyed. Deep pencil-beams in our L1 and L2 models have more chance than in t1 to hit a rich cluster at redshift $`0.30.5`$. Then high peaks in the power spectra tend to appear on small wavenumbers in L1 and L2 for the reason explained above.

Overall, we conclude that although roughly half a per cent of our mock surveys give a power spectrum peak stronger than that of the BEKS data, very few of these actually correspond to redshift distributions with similar regularity and similar spacing of the spikes. We now study this further by considering two additional tests for regularity which have been used on the BEKS sample.

## 6 PG $`\mathrm{\Delta }`$-test

In comparing their own simulation to the BEKS data, Park & Gott (1991) made use of a test specifically designed to probe the apparent “phase-coherence” of the series of redshift spikes. For each ‘galaxy’ they calculated the distance to the nearest tooth of a perfectly regular comb-like template. They then ratioed the mean of this distance to the separation of the teeth, and minimised the result over the period and phase of the template. Let us call the resulting statistic $`\mathrm{\Delta }`$. Then a distribution in which each galaxy is at some node of a regular grid will give $`\mathrm{\Delta }=0`$, and a uniform distribution in distance would give $`\mathrm{\Delta }=0.25`$ in the large-sample limit. In our application of this test we restrict the range of possible periods to 50 – 500 $`h^1`$ Mpc. For the BEKS data we obtain $`\mathrm{\Delta }`$=0.165 for a best period of 130$`h^1`$ Mpc. Our value of $`\mathrm{\Delta }`$ differs from that given by PG because they applied the test only to the deep surveys while we use the combined deep and shallow BEKS data. Among the 1936 samples in our t1 ensemble, 209 have lower values of $`\mathrm{\Delta }`$ than the BEKS data; for the L1 and L2 ensembles the corresponding numbers are 134/2625 and 127/2625 respectively. According to this test, therefore, the observed sample appears only marginally more regular than expected in our CDM cosmologies.

Within each of our ensembles the significance of the regularity in the BEKS data appears somewhat higher than was estimated by PG. Their simulation ensemble was made up of only 12 mock catalogues of which one had lower $`\mathrm{\Delta }`$ than the BEKS deep data. The median $`\mathrm{\Delta }`$ for these twelve was 0.1695, while we find medians of 0.176, 0.180 and 0.189 for the combined deep and shallow data in ensembles t1, L1 and L2 respectively. The difference with PG is probably small enough to be attributed to the small number of realisations in their ensemble. Figure 9 shows the distribution of periods for catalogues in each of our own ensembles with lower $`\mathrm{\Delta }`$s than the BEKS data. It is interesting that relatively small periods are favoured and that there is no preference for values in the range $`120130h^1`$Mpc. Ettori et al. (1997) used a related test, the comb-template test (Duari et al. 1992), to analyse four pencil-beam surveys near the South Galactic Pole. They found a best period near 130 $`h^1`$Mpc in two of these four directions, in apparent agreement with the BEKS result.

In summary, the difference in regularity between the BEKS sample and our CDM mock-catalogues is less significant when measured by the $`\mathrm{\Delta }`$-test than when measured using the power spectrum test of the last section. Nevertheless, for periods near $`125h^1`$Mpc there are few CDM samples which are more regular than the BEKS data. In addition, our result appears insensitive to the choice of biasing; we find essentially no difference between L1 and L2 in Figure 9. This is puzzling since Figure 3 shows clear differences in the emptiness of the voids in the two cases. Apparently the value of $`\mathrm{\Delta }`$ is more sensitive to departures from regularity in the spacing of the walls than it is to the density contrast of the voids.

## 7 Supercluster statistics

Dekel et al. (1992) proposed an alternative technique for assessing apparent periodicity in data samples like that of BEKS. In this section we use the term ‘supercluster’ to refer to clumps in the one-dimensional redshift histograms derived from such pencil-beam surveys, even though these do not correspond precisely to the superclusters (or walls or filaments) seen in fully three-dimensional surveys. The method of Dekel et al. is based on the redshift distribution of supercluster centres rather than on that of individual galaxies. The first step is to correct the galaxy redshift histogram for the survey selection function. We do this by weighting each galaxy by the inverse of the selection function for the particular survey of which it is a part (North or South, shallow or deep). This reverses the procedure by which we created our mock catalogues from the simulations. To avoid overly large sampling noise where the selection function is small, we restrict our redshift histograms to $`z0.31`$ for the SGP survey and $`z0.5`$ for the NGP survey (see Dekel et al. 1992). We smooth these histograms with a Gaussian of variance $`l^2`$ and identify supercluster centres as local maxima of the result. (Note that, following Dekel et al., no threshold is applied.) We have tried smoothing lengths $`l`$ between 20 and 40 $`h^1`$Mpc, but find our results to be insensitive to the exact value within this range. In what follows we set $`l`$ =25 $`h^1`$Mpc.

Given a distribution of the supercluster centres, the characteristic period is determined in the following way. As a first estimate we take the mean separation $`L_m`$ between neighbouring superclusters. Next we apply the PG $`\mathrm{\Delta }`$-test for periods $`p[0.5L_m,2.0L_m]`$. The value of $`p`$ in this range which minimises $`\mathrm{\Delta }`$ is taken as the characteristic period of the distribution. For this period we calculate the Rayleigh statistic $`R`$ as follows (Dekel et al. (1992) and Feller(1971)). The positions of the supercluster centres are mapped onto a circle of circumference $`p`$. Consider the $`n`$ unit vectors $`𝐮_i`$ which point from the centre of the circle towards each of the $`n`$ superclusters. Denote their vector average by $`𝐮`$, the modulus of $`𝐮`$ by $`V`$ and define $`R=1V`$. For an exactly periodic distribution the unit vectors would all be identical so that $`V=1`$ and $`R=0`$. For a distribution with no long-range phase coherence the directions of the unit vectors would be random and so in the large sample limit $`V0`$ and $`R1`$. Small values of $`R`$ are thus expected for near-periodic distributions.

For the BEKS data, we find $`R=0.33`$ for a period of 130 $`h^1`$Mpc obtained as described above. Lower values of $`R`$ are found for 46 samples in t1, for 66 samples in L1, and for 56 samples in L2. Thus, according to this test the supercluster distribution in the BEKS data is more periodic than the CDM models at the 2.4 per cent significance level for t1, the 2.5 per cent level for L1 and the 2.1 per cent level for L2. The period distribution of the samples with $`RR`$(BEKS) is shown in Figure 10. Many of these low-$`R`$ samples have periods in the range \[100$`h^1`$Mpc, 140$`h^1`$Mpc\]. Thus in CDM model universes it is common for supercluster spikes to have a typical separation similar to that seen in the BEKS data and in a few per cent of cases the spikes are just as regularly spaced as in the real data.

## 8 Discussion and conclusion

By creating a number of mock pencil-beam surveys we have compared the apparent periodicity in two CDM model universes with that observed in the data of Broadhurst et al. (1990). The power spectrum analysis alone shows that the BEKS data are significantly more periodic than the models at about the half per cent level, while the PG $`\mathrm{\Delta }`$-test shows less significance, about 10 per cent for t1 and 5 per cent for L1 and L2. The supercluster statistic gives a two per cent probability of finding a structure as regular as the BEKS data. Restricting to a length scale $``$100-150$`h^1`$Mpc, however, the number of samples which show the kind of periodicity seen in the BEKS data is extremely small for each of these statistics. Overall no sample is more regular than the BEKS data for all three statistics for a single period. The two popular CDM models we have studied here are apparently unsuccessful in reproducing the observed periodicity. From this result together with the fact that the statistical results appeared to be insensitive to the choice of the bias model, we conclude that CDM models conflict with the BEKS observation. Either the models need additional physics, or the data are a fluke or are somehow biased.

Various possible physical explanations have been proposed, such as coherent peculiar velocities (Hill, Steinhardt and Turner 1991) oscillations in the Hubble parameter (Morikawa 1991) or baryonic features in the power spectrum (Eisenstein et al. 1998) but all of them seem to require either additional mechanisms with fine tunings beyond the standard theory or cosmological parameters significantly different from currently favoured values. Intriguingly, Dekel et al. (1992) demonstrated that built-in power on a large ($`100h^1`$Mpc) length scale in the initial density fluctuation could indeed reproduce periodic features on a given scale, at least by some of the tests we have considered. If such excess power on large scales (hence still in the linear regime) exists, it will be detectable in the power spectra of the future 2dF and Sloan surveys.

Having found at least a few examples that are nearly as periodic as the BEKS data, we cannot rule out the possibility that the BEKS data (or the Galactic Pole direction) are a fluke. On the other hand, one should be aware of the complexity of the original observations –an incomplete compilation of a narrow and deep, and of a wide and shallow survey at each of the Galactic Poles. It is not clear whether such a combination constitutes a fair sample. No evidence for periodic structure on $``$130$`h^1`$Mpc has been found so far in two other deep redshift surveys, the ESO-Sculptor Survey (Bellanger and de Lapparent 1995) and the Caltech Faint Galaxy Redshift Survey (Cohen 1999). Follow-up observations to BEKS by Koo et al. (1993) did not show a strong regularity in two other directions, although around the Galactic Pole the regularity was found to be further strengthened. Our results give the a priori probability for such apparent periodicity in CDM models. Several more deep surveys might suffice to judge whether the discrepancy with BEKS reflects a major inconsistency. The planned VIRMOS Deep Survey (Le Fèvre et al. 1998, see also Guzzo 1999) will survey the range $`0.3z1`$ and will provide, together with the large volume 2dF and Sloan surveys, much larger and more complete samples in the near future.

## Acknowledgments

We are grateful to Carlton Baugh and Shaun Cole for fruitful discussions. We thank the referee, Avishai Dekel, for giving us constructive comments on the manuscript. The simulations were carried out on the Cray-T3E at the Rechenzentrum, Garching (RZG).