text:
# 1 Introduction

## 1 Introduction

Many searches for massive new particles predicted by extensions to the Standard Model (SM) assume that these particles decay promptly at the primary interaction vertex. Such searches are not sensitive to long-lived heavy particles which do not decay within the detectors. There exist, however, a number of models which predict such long-lived particles. For example, in the Constrained Minimal Supersymmetric Standard Model (CMSSM), for certain choices of the parameters, sleptons or charginos could be long-lived . R-parity violating supersymmetric (SUSY) models  also allow for long-lived heavy particles and a fourth-generation heavy lepton could be stable . In gauge-mediated supersymmetry, if the SUSY-breaking energy scale is sufficiently high , sleptons could be long-lived. Some models beyond the SM also predict the existence of particles with fractional electric charge. Previous searches for long-lived massive charged particles have been performed by the LEP collaborations with data taken at the Z0 resonance , as well as with data taken at higher centre-of-mass energies, up to 209 GeV .

This paper describes an update to a search for long-lived particles, referred to here as X±, with $`m_\mathrm{X}>m_\mathrm{Z}/2`$, and charge $`|Q/e|`$ = 1 or 2/3, pair-produced in the reaction $`\mathrm{e}^+\mathrm{e}^{}\mathrm{X}^+\mathrm{X}^{}`$. This search has been described in detail in . In this paper we also search for particles with fractional charges 4/3 and 5/3. All fractionally-charged particles are assumed to be colourless and non-strongly-interacting. To make the search for these particles as model independent as possible, only minimal calorimetric information has been used. Due to their large mass these particles would have anomalously high or low ionization energy loss $`\mathrm{d}E/\mathrm{d}x`$ in the tracking chambers. This search is therefore primarily based on the precise $`\mathrm{d}E/\mathrm{d}x`$ measurement provided by the OPAL jet chamber. The data were collected by the OPAL detector during 1995-2000, at centre-of-mass energies from 130 GeV to 209 GeV corresponding to a total integrated luminosity of 693.1 pb-1 as reported in Table 1.

## 2 The OPAL Detector

A complete description of the OPAL detector can be found in Ref. . Here only a brief overview is given. The central detector comprised a system of tracking chambers, providing track reconstruction over 96% of the full solid angle111The OPAL right-handed coordinate system is defined such that the $`z`$-axis is in the direction of the electron beam, the $`x`$-axis points towards the centre of the LEP ring, and $`\theta `$ and $`\varphi `$ are the polar and azimuthal angles, defined relative to the $`+z`$\- and $`+x`$-axes, respectively. The radial coordinate is denoted by $`r`$. inside a 0.435 T uniform magnetic field parallel to the beam axis. It consisted of a two-layer silicon microstrip vertex detector, a high-precision vertex drift chamber (CV) with axial and stereo wires, a large-volume jet chamber and a set of $`z`$-chambers measuring the track coordinates along the beam direction.

The jet chamber (CJ) is the most important detector for this analysis. It was divided into 24 azimuthal sectors, each equipped with 159 sense wires. Up to 159 position and $`\mathrm{d}E/\mathrm{d}x`$ measurements per track were thus possible, with a precision of $`\sigma _{r\varphi }135\mu `$m and $`\sigma _z`$ 6 cm. When a track was matched with $`z`$-chamber hits and hits on the stereo wires of the vertex chamber (CV), the uncertainty on its $`z`$ coordinate was $``$ 1 mm. The tracking detectors, located inside the magnet coil, provided a track momentum measurement with a resolution of $`\sigma _p/p\sqrt{(0.02)^2+(0.0015p_t)^2}`$ for tracks with the full number of hits ($`p_t`$, in GeV, is the momentum transverse to the beam direction) and a resolution on the ionization energy loss measurement of approximately $`2.8\%`$ for $`\mu ^+\mu ^{}`$ events with a large number of usable hits for $`\mathrm{d}E/\mathrm{d}x`$ measurement .

A lead-glass electromagnetic calorimeter (ECAL) located outside the magnet coil covered the full azimuthal range with good hermeticity in the polar angle range of $`|\mathrm{cos}\theta |<0.984`$. The magnet return yoke was instrumented for hadron calorimetry covering the region $`|\mathrm{cos}\theta |<0.99`$ and was surrounded by four layers of muon chambers. Electromagnetic calorimeters close to the beam axis completed the geometrical acceptance down to 24 mrad on both sides of the interaction point.

The ionization energy loss $`\mathrm{d}E/\mathrm{d}x`$ produced by a charged particle is a function of the electric charge $`Q`$ and of $`\beta \gamma =p/m`$, where $`p`$ is the momentum and $`m`$ the mass of the particle . Figure 1 shows the distribution of $`\mathrm{d}E/\mathrm{d}x`$ as a function of the apparent momentum, $`p/Q`$. Standard particles of charge $`\pm `$1 (e, $`\mu `$, $`\pi `$, p, K) with high momentum ($`p>0.1\sqrt{s}`$ ) have $`\mathrm{d}E/\mathrm{d}x`$ between 9 and 11 keV/cm. Massive particles with charge $`\pm `$1 are expected to yield $`\mathrm{d}E/\mathrm{d}x>11`$ keV/cm for high-mass values, $`m_\mathrm{X}>0.36\sqrt{s}`$, or $`\mathrm{d}E/\mathrm{d}x<9`$ keV/cm for low-mass values, $`m_\mathrm{X}<0.27\sqrt{s}`$. The $`\mathrm{d}E/\mathrm{d}x`$ measurement therefore provides a good tool for particle identification in these high- and low-mass regions. Massive particles with charge $`\pm `$2/3 would have $`\mathrm{d}E/\mathrm{d}x>11`$ keV/cm for high mass values, $`m_\mathrm{X}>0.45\sqrt{s}`$ or $`\mathrm{d}E/\mathrm{d}x<9`$ keV/cm for low-mass values, $`m_\mathrm{X}<0.35\sqrt{s}`$. The expected $`\mathrm{d}E/\mathrm{d}x`$ for massive particles of charge $`\pm `$4/3 and $`\pm `$5/3 is greater than 11 keV/cm for all mass ranges. The search for massive particles with charge $`\pm `$1/3 was not possible because the typical $`\mathrm{d}E/\mathrm{d}x`$ deposit of these particles would be too close to the instrumental noise level.

## 3 Monte Carlo Simulation

Several different Monte Carlo programs were used to generate the signal process $`\mathrm{e}^+\mathrm{e}^{}\mathrm{X}^+\mathrm{X}^{}`$. Signal events of the type $`\mathrm{e}^+\mathrm{e}^{}\stackrel{~}{\mathrm{}}^+\stackrel{~}{\mathrm{}}^{}`$ ($`\stackrel{~}{\mathrm{}}^\pm `$ being a charged scalar lepton) were generated at several centre-of-mass energies up to 206 GeV using SUSYGEN . The generated charged scalar leptons are not allowed to decay, therefore simulating a signal from heavy charged stable scalar particles. Similarly, events of the type $`\mathrm{e}^+\mathrm{e}^{}\mathrm{L}^+\mathrm{L}^{}`$ and $`\mathrm{e}^+\mathrm{e}^{}\mathrm{Q}\overline{\mathrm{Q}}`$, where L± are stable heavy spin-1/2 leptons, and Q, $`\overline{\mathrm{Q}}`$ are colourless stable heavy spin-1/2 particles with charge 2/3, 4/3, and 5/3, were generated at the same energies, using the EXOTIC  generator. All signal samples were generated with 1000 events per mass point with mass $`m_\mathrm{X}`$ ranging from 45 GeV to the kinematic limit for the centre-of-mass energy considered. The mass points were generated every 5 GeV with a finer binning of 1 GeV in the mass regions where we expect lower selection efficiencies. For the purpose of detector simulation and particle interactions, all particles were treated as heavy muons. The simulation of $`\mathrm{d}E/\mathrm{d}x`$ in the central jet chamber of OPAL accounted for the charge, mass and momentum of the particle, as described above.

The background was estimated using simulations of all Standard Model processes (two-fermion, four-fermion and two-photon processes) for all centre-of-mass energies from 130 to 206 GeV. Small differences in the centre-of-mass energies of data and background Monte Carlo samples have a negligible effect on the analysis.

The contribution to the background from two-fermion final states was estimated using BHWIDE  for the $`\mathrm{e}^+\mathrm{e}^{}`$ final states and KORALZ  and KK2f  for the $`\mu ^+\mu ^{}`$ and $`\tau ^+\tau ^{}`$ states. Hadronic two-fermion events, $`\mathrm{q}\overline{\mathrm{q}}`$, were simulated using PYTHIA . For the two-photon background, the PYTHIA , PHOJET  and HERWIG  Monte Carlo generators were used for $`\mathrm{e}^+\mathrm{e}^{}\mathrm{q}\overline{\mathrm{q}}`$ final states and the Vermaseren  and the BDK  generators for all $`\mathrm{e}^+\mathrm{e}^{}\mathrm{}^+\mathrm{}^{}`$ final states. Four-fermion final states were simulated with grc4f , which takes into account interference between all diagrams. All generated signal and background events were processed through the full simulation of the OPAL detector ; the same event analysis chain was applied to the simulated events and to the data.

## 4 Data Analysis

Pair-produced stable or long-lived massive charged particles would manifest themselves in events with two back-to-back tracks. Assuming they would not interact strongly, these particles would not produce hadronic showers. Since they are massive, they would not produce electromagnetic showers either. For these reasons, the considered events would be very similar to $`\mu ^+\mu ^{}`$ events, the only difference being the higher mass of the particles, which yields a different $`\mathrm{d}E/\mathrm{d}x`$ for the same momentum.

A preselection similar to the one described in  is used for the analyses. Several criteria have been loosened in order to increase the sensitivity to high-mass particles. The criteria are listed below:

* Events are rejected if the total multiplicity of tracks in the central detector and clusters in the ECAL is greater than 26. Cosmic ray events are rejected as in . Bhabha events are identified by requiring two energetic and collinear clusters in the electromagnetic calorimeter, these events are then rejected.
* Events are required to contain at least two tracks in the central detector, each satisfying basic quality criteria222The distance between the beam axis and the track at the point of closest approach (PCA) must be less than 1 cm; the $`z`$-coordinate of the PCA must be less than 40 cm; the innermost hit of the track measured by the jet chamber must be closer than 75 cm to the beam axis. and having a momentum $`p>0.07\sqrt{s}`$, a momentum transverse to the beam axis $`p_t>0.025\sqrt{s}`$, a polar angle satisfying $`|\mathrm{cos}\theta |<0.97`$ and at least 20 CJ hits usable for $`\mathrm{d}E/\mathrm{d}x`$ measurement. The two selected tracks are required to have opposite electric charge.
* To reduce background from two-photon events, the total visible energy333The visible energy, the visible mass and the total transverse momentum of the event are calculated using tracks and calorimeter clusters, correcting for double counting as described in . of the event is required to be $`E_{\mathrm{vis}}>0.14\sqrt{s}`$ and the acoplanarity angle444The acoplanarity angle, $`\varphi _{\mathrm{acop}}`$, is defined as 180 minus the angle between the two tracks in the $`r\varphi `$ plane. between the two tracks is required to be $`\varphi _{\mathrm{acop}}<20^{}`$.
* Events containing an isolated ECAL cluster with an energy greater than 5 GeV are rejected to reduce background from events with initial state radiation. Isolation is defined as an angular separation of more than $`15^{}`$ from the closest track.
* It is required that $`\frac{E_1}{p_1}+\frac{E_2}{p_2}<0.2`$, where $`p_{1,2}`$ are the momenta of the two selected tracks and $`E_{1,2}`$ denote the energies of the ECAL clusters associated to the tracks. This further reduces the contribution from Bhabha scattering events. No other tracks with $`p>0.5`$ GeV and no unassociated clusters with $`E>3`$ GeV should be found in a cone of $`10^{}`$ half-opening angle around each of the two selected tracks. Criterion P5 is not used in the fractionally-charged particles analysis since their interaction properties with the calorimeters are unknown. This reduces the dependence on calorimeter response around the candidate tracks.

After the preselection, the background is dominated by $`\mathrm{e}^+\mathrm{e}^{}\mu ^+\mu ^{}`$ events, with a small contribution from $`\mathrm{e}^+\mathrm{e}^{}\tau ^+\tau ^{}`$ and two-photon $`\mathrm{e}^+\mathrm{e}^{}\mu ^+\mu ^{}`$ events. The effect of the preselection cuts on the all data samples and various Monte Carlo background processes, for the search for $`|Q/e|=`$ 1 particles, is shown in Table 2.

### 4.1 Search for particles with unit charge

The search strategy has been simplified with respect to  and now relies entirely on $`\mathrm{d}E/\mathrm{d}x`$ information. A sample $`\mathrm{d}E/\mathrm{d}x`$ distribution for data and simulated Monte Carlo events is shown in Figure 2. The kinematic selection is no longer used. This new strategy has been applied to all data samples. The preselected events are retained if they satisfy the following requirements on $`\mathrm{d}E/\mathrm{d}x`$:

* Both high-momentum tracks must have either $`\mathrm{d}E/\mathrm{d}x>11`$ keV/cm or $`\mathrm{d}E/\mathrm{d}x<9`$ keV/cm.
* The probability that the $`\mathrm{d}E/\mathrm{d}x`$ measurements for either track were consistent with one of the standard particles (e, $`\mu `$, $`\pi `$, p, K) must be less than 30%. This removes background from poorly measured SM particles.

The effect of cuts A1-A2 for all data and simulated events can be seen in Table 2. No candidate event is found. The expected backgrounds at the various centre-of-mass energies are shown in Table 1. The total background, summed over all energies, is estimated to be 1.1$`\pm `$1.3 events, where the uncertainty quoted includes both statistical and systematic effects.

The detection efficiency for spin-0 particles varies between 75 and 90% for masses $`m_\mathrm{X}<0.27\sqrt{s}`$ or $`m_\mathrm{X}>0.36\sqrt{s}`$. The efficiency drops significantly in the region $`0.27\sqrt{s}<m_\mathrm{X}<0.36\sqrt{s}`$, but as data sets collected at different centre-of-mass energies ($`\sqrt{s}`$) are combined, a reasonable selection efficiency is achieved for all mass values up to close to the kinematic limit. The selection efficiencies of all Monte Carlo samples at all $`\sqrt{s}`$ are parametrised as a function of $`\beta \gamma =p/m_\mathrm{X}=\sqrt{s/4m_X^21}`$ of the particle. This parametrisation is used to calculate the efficiency for all masses at each centre-of-mass energy, using linear interpolation. For spin-1/2 particles, the efficiencies are 2-9% lower than for spin-0 particles due to the different angular distribution of the tracks. We analyse each centre-of-mass energy separately, then combine the results for the final cross-section limits.

### 4.2 Search for particles with fractional charges

To search for particles with fractional charges of 2/3, 4/3, and 5/3 the selection criteria P1 through P4, followed by A1 and A2 are used. The results after each cut for all data samples and various Monte Carlo background processes are shown in Table 3. For charge 2/3 the selection efficiency is between 75 and 90% for most of the mass range, while for charges 4/3 and 5/3 the efficiency is above 90% for the whole mass range. After this selection, one candidate survives in the data sample at $`\sqrt{s}=183`$ GeV, one at $`\sqrt{s}=189`$ GeV, and one at $`\sqrt{s}=207`$ GeV, while no candidate is left in any of the other data sets. The masses of the candidates are reconstructed for charge 2/3 using the $`\mathrm{d}E/\mathrm{d}x`$ and momentum information, while for charge 4/3 and 5/3 kinematic information only is used. The reconstructed masses, track momenta and $`\mathrm{d}E/\mathrm{d}x`$ values of the candidate events are reported in Table 4. The selected events and the expected background at each centre-of-mass energy are shown in Table 1. The total background, summed over all energies, is estimated to be 3.2$`\pm `$2.4 events, where the uncertainty quoted includes both statistical and systematic effects.

### 4.3 Systematic uncertainties

The main sources of systematic uncertainties for this analysis are discussed below and reported in Table 5, for both the signal efficiency and the background estimate:

* The errors on the Monte Carlo modeling of $`\varphi _{\mathrm{acop}}`$, $`E_{\mathrm{vis}}`$ and $`\mathrm{d}E/\mathrm{d}x`$ are estimated by comparing the distributions of these variables for data and background Monte Carlo. The relative difference between the averages of the distributions is used to increase or decrease the value of the cut on the relevant variable in order to decrease the overall signal efficiency and background estimate. The difference between the reduced efficiency and the one obtained with the nominal selection is taken as the systematic uncertainty due to the modeling of the variable under consideration. This estimate of the error is more conservative than the one obtained by smearing the $`\mathrm{d}E/\mathrm{d}x`$ values of the tracks.
* The MC statistical uncertainty, due to the limited number of signal events generated, has been computed using a binomial formula. For background processes the large relative statistical uncertainty is due to the limited number of background events selected.
* The uncertainty due to the linear interpolation of the signal detection efficiency is estimated as the difference between the interpolated values and the efficiency obtained at mass points where MC signal samples were generated, when that mass point was omitted from the interpolation procedure.
* For background processes of the type $`\mathrm{e}^+\mathrm{e}^{}\mathrm{}^+\mathrm{}^{}`$, the Vermaseren  generator has been used as the reference generator. The difference in the background expectation obtained by using BDK  instead of Vermaseren has been considered as the uncertainty on the MC generator.
* Two tracks separated by a distance smaller than 2.5 mm could be unresolved and reconstructed as a single track with a high $`\mathrm{d}E/\mathrm{d}x`$ value. Events with unresolved double tracks are potential backgrounds for this search. The systematic uncertainty introduced by the modeling of the double track resolution in the Monte Carlo samples has been estimated to be 2.9 percent.

The absolute uncertainty on the background is reported by centre-of-mass energy in Table 1. The uncertainty introduced on the integrated luminosity  is also reported in Table 5. At a given centre-of-mass energy the different systematic uncertainties are assumed to be independent, so that the total systematic uncertainty is calculated as the quadratic sum of the individual uncertainties. The Monte Carlo modeling of $`\mathrm{d}E/\mathrm{d}x`$, the Vermaseren-BDK generator, and luminosity uncertainties are assumed to be correlated between centre-of-mass energies, while the other systematic uncertainties are assumed to be independent.

## 5 Results

The numbers of candidates found in the search for particles with charge $`|Q/e|`$=1 and fractional charges are summarised for all energies in Table 1, together with the expected backgrounds. The data show no significant excess above the expected background from Standard Model processes.

Model-independent cross-section upper limits have been computed for the pair-production of massive charged long-lived particles, combining the results from all centre-of-mass energies, assuming s-channel production. The cross-section dependence on the energy is taken to be $`\beta ^3/s`$ for spin-0 particles and $`\frac{\beta }{s}(1\frac{\beta ^2}{3})`$ for spin-1/2 particles, where $`\beta =p/E\sqrt{14m_\mathrm{X}^2/s}`$. In evaluating upper limits, the candidates are counted in mass intervals centred on their central mass values and $`\pm 2\sigma `$ wide (where $`\sigma `$ is the error on the mass estimate of the candidates as reported in Table 4). In the case in which the mass could not be reconstructed the candidates were considered in the whole mass range (from 45 GeV to the kinematic limit for that event). A likelihood-ratio method  was used to determine an upper limit for the cross-section. The total systematic error is incorporated into the limits following the prescription of Ref. .

In Figure 3(a), the 95% C.L. upper limit on the cross-section at $`\sqrt{s}=206.6`$ GeV is shown for spin-0 particles of charge $`\pm `$ 1. The 95% C.L. upper limit on the pair-production cross-section varies from 0.005 to 0.028 pb in the mass range $`45<m_\mathrm{X}<101`$ GeV. The bump in the mass range of $`52<m_\mathrm{X}<70`$ GeV is due to the low efficiencies described in section 4.1. The cross-section limits are compared with the predicted cross-sections  for pair-production of right- and left-handed smuons and staus to determine mass limits. For these two slepton species, the production cross-section does not depend on the CMSSM parameters but only on the slepton mass. The 95% C.L. lower mass limits are 98.0 GeV and 98.5 GeV for the mass of right- and left-handed smuons and staus, respectively, as shown in Figure 3(a).

Figure 3(b) shows the 95% C.L. upper limit on the cross-section at $`\sqrt{s}=206.6`$ GeV for spin-1/2 particles of charge $`\pm `$1. The limit varies from 0.005 to 0.024 pb in the mass range $`45<m_\mathrm{X}<100`$ GeV. This limit must be compared with the predicted cross-sections for chargino production  and heavy charged lepton production . For the chargino limits, the CMSSM parameters have been chosen to minimise the predicted chargino cross-section at every chargino mass value (assuming a heavy sneutrino, $`m_{\stackrel{~}{\nu }}>500`$ GeV), without any restriction on the mass of the lightest neutralino. A 95% C.L. lower limit on the masses of long-lived charginos, of 102.0 GeV, is obtained for every choice of the CMSSM parameters. The 95% C.L. lower limit on the heavy charged lepton mass is also 102.0 GeV.

Figure 4(a) shows the 95% C.L. upper limit on the cross-section at $`\sqrt{s}=206.6`$ GeV for spin-1/2 particles of charge $`\pm `$2/3. The limit varies between 0.005 and 0.020 pb in the mass range $`45<m_\mathrm{X}<101`$ GeV. Figures 4(b) and 4(c) show the 95% C.L. upper limit on the cross-section at $`\sqrt{s}=206.6`$ GeV for spin-1/2 particles of charge $`\pm `$4/3 and of charge $`\pm `$5/3, respectively. For spin-0 particles with fractional charge the cross-section upper limits are slightly higher than for spin-1/2 particles. This is due to the difference in angular distributions.

All results obtained are valid for non-strongly-interacting colourless particles with a lifetime longer than 10-6 s. This lifetime restriction is obtained by considering the heaviest (and therefore slowest) particles excluded by this search, and then requiring that the decay probability of these particles at a flight distance larger than 3.0 m be greater than 95%. For lower mass values the results are also valid for shorter lifetimes.

## 6 Summary and Conclusions

A search was performed for pair-production of stable and long-lived massive particles not subject to strong interactions, with charge $`|Q/e|`$ = 1 or fractional charges of $`2/3`$, $`4/3`$, and $`5/3`$. The primary tool used in this search was the precise $`\mathrm{d}E/\mathrm{d}x`$ measurement provided by the OPAL jet chamber. No evidence for the production of such particles was observed. For s-channel production, the upper limits on the cross-section vary between 0.005 and 0.026 pb in the mass range explored for particles of charge $`\pm `$1. Within the framework of the CMSSM, lower mass limits have been obtained: 98.0 GeV for right-handed and 98.5 GeV for the left-handed smuons and staus. Charged long-lived massive leptons and long-lived charginos with masses smaller than 102.0 GeV are excluded. For particles with fractional charge $`\pm `$2/3, $`\pm `$4/3 and $`\pm `$5/3, the upper limits on the production cross-section vary between 0.005 and 0.020 pb in the range $`45<m_\mathrm{X}<95`$ GeV. The above limits are valid at the 95% C.L. for particles with lifetimes longer than 10-6 s.

## Acknowledgements

We particularly wish to thank the SL Division for the efficient operation of the LEP accelerator at all energies and for their close cooperation with our experimental group. In addition to the support staff at our own institutions we are pleased to acknowledge the
Department of Energy, USA,
National Science Foundation, USA,
Particle Physics and Astronomy Research Council, UK,
Natural Sciences and Engineering Research Council, Canada,
Israel Science Foundation, administered by the Israel Academy of Science and Humanities,
Benoziyo Center for High Energy Physics,
Japanese Ministry of Education, Culture, Sports, Science and Technology (MEXT) and a grant under the MEXT International Science Research Program,
Japanese Society for the Promotion of Science (JSPS),
German Israeli Bi-national Science Foundation (GIF),
Bundesministerium für Bildung und Forschung, Germany,
National Research Council of Canada,
Hungarian Foundation for Scientific Research, OTKA T-038240, and T-042864,
The NWO/NATO Fund for Scientific Reasearch, the Netherlands.


format: text

text:
# Monopole clusters at short and large distances

## I Introduction

The interest in monopoles in non-Abelian gauge theories is mostly due to the dual superconductor mechanism of the confinement, for review see reviews . The mechanism assumes condensation of magnetic monopoles in the vacuum. The idea is supported by the lattice data and the phenomenology of the lattice monopoles is quite rich. Especially, in case of $`SU(2)`$ gluodynamics which we will concentrate on in this paper.

Detailed theoretical interpretations of the data appear, however, difficult because monopoles are defined not directly in terms of the original $`SU(2)`$ fields but rather in terms of projected fields. The use of a projection is rooted in the fact that monopoles are intrinsically $`U(1)`$ objects and there are infinitely many ways to select a $`U(1)`$ subgroup for the monopole definition. In particular, monopoles of the maximal Abelian projection (MAP) are defined in the following sequence of steps (which we describe in a somewhat simplified way). First, one finds the maximal Abelian gauge maximizing the functional:

$$R_{Abel}=d^4x\left[\left(A_\mu ^1(x)\right)^2+\left(A_\mu ^2(x)\right)^2\right]$$
(1)
where $`A_\mu ^a(x)`$ are components of the gauge field. Finally, monopoles are defined as singularities of the $`A_\mu ^3(x)`$ fields which violate the Bianchi identities:

$$j_\mu ^{mon}(x)=ϵ^{\mu \nu \rho \sigma }_\nu _\rho A_\sigma ^3(x),$$
(2)
where $`j_\mu ^{mon}(x)`$ is the monopole current and all the expressions are well defined on the lattice.

Geometrically, monopoles are represented by closed trajectories (on the dual lattice) and one usually discusses properties of the monopole clusters. In particular, it is important to distinguish between the percolating cluster and finite clusters. The percolating cluster fills in the whole of the lattice and is in a single copy for each field configuration. Finite clusters are characterized, in particular, by their spectrum as function of the length $`l`$. Commonly, one introduces the corresponding densities of the percolating and finite clusters:

$$l_{perc}4\rho _{perc}a^4N_{sites},l_{fin}4\rho _{fin}a^4N_{sites},$$
(3)
where $`l_{perc},l_{fin}`$ is the total length of the corresponding clusters, $`N_{sites}`$ is the number of lattice sites and $`a`$ is the lattice spacing.

An important question is whether position of the monopoles is of any physical significance. On one hand, monopole trajectories are distinguished by singularities of the projected fields, see (2), and they are point like in terms of the projected fields. On the other hand, one can suspect that these singularities are artifacts of the projection.

Probably, a priori the latter possibility looks more reasonable since the monopole definition involves a non-local gauge fixation, see (1). However, there are accumulating lattice data which indicate that in the maximal Abelian gauge the monopoles, which by construction have size of the lattice spacing $`a`$, are physical. In particular, the density of the percolating monopoles scales:

$$\rho _{perc}=0.65(2)\sigma _{SU(2)}^{3/2},$$
(4)
where $`\sigma _{SU(2)}`$ is the string tension 111In estimates we will also use $`\mathrm{\Lambda }_{QCD}`$ although we study actually $`SU(2)`$ gluodynamics, $`\sigma _{SU(2)}\mathrm{\Lambda }_{QCD}^2`$. and we quoted the data from Ref. mueller  where references to earlier papers can also be found.

Let us emphasize that the result (4) implies that the corresponding monopole trajectories are meaningful even on the scale $`a`$. Indeed, the total length of the percolating clusters which scales in the physical units is added up from small steps of size $`a`$ and finally is independent on $`a`$. Thus, one could argue, see, e.g., hart , that there are gauge invariant objects which are detected through the projection. It is of course an intriguing hypothesis which is worth to be thoroughly checked. Unfortunately, this suggestion cannot be scrutinized theoretically since the anatomy of the lattice monopoles in terms of the non-Abelian fields is largely unknown, for review see, e.g., coleman . Rather one should apply further phenomenological tests.

In particular, one can measure non-Abelian action associated with the monopoles, see anatomy  and references therein. For monopoles in the percolating clusters it turns out that, at least at presently available lattices, the action corresponds to a monopole mass which diverges in the ultraviolet:

$$M(a)\frac{const}{a}.$$
(5)
$`M(a)`$ is the lattice analog of the ‘magnetic mass’ in the continuum, $`M(a)𝐇^2d^3r`$ where $`𝐇`$ is the magnetic field. The result (5) strongly suggests that the lattice monopoles are associated with singular non-Abelian fields and the singularity (2) is not an artifact of the projection.

The next question is how singular fields can be important at all since their contribution in the limit $`a0`$ is suppressed by an infinite action 222Hereafter we tacitly assume that the $`a`$ dependence observed at the presently available lattices will persist at smaller lattice spacings as well. The smallest $`a`$ available now is about $`a(3GeV)^1`$.. An apparent answer to this question is that the action and entropy factor are both divergent in the ultraviolet but cancel each other in such a way that, say, the total density of monopoles in the percolating cluster scales in the physical units, see (4). An example of such a fine tuning is provided by $`U(1)`$ theory, see, in particular, polyakov . In that case the fine tuning is ensured by choosing tuned values of the electric charge. In case of non-Abelian theories the hypothesis on the fine tuning vz ; maxim  is rather phenomenological and the mechanism of the fine tuning is to be understood yet.

To summarize,

(1) there is evidence that at presently available lattices the monopoles defined within the maximal Abelian projection are associated in fact with singular non-Abelian fields fine tuned to the corresponding entropy factors.

(2) This means that the structure of the non-perturbative vacuum fluctuations is quite different from the standard picture of ‘bulky’ fields of the size $`\mathrm{\Lambda }_{QCD}^1`$.

(3) Although the conclusion on the fine tuning is strongly suggested by the data like (4), (5) further checks of it are highly desirable. Indeed, the whole basis is the phenomenology of the lattice monopoles and further data can bring further insight.

Motivated by these considerations, we have investigated in more detail the scaling properties of the monopole clusters 333The preliminary numerical results were presented at the “Lattice 02” conference latt02 .. The outline of the paper is as follows. In Sect. II we present detailed studies of the geometrical elements of the percolating cluster. Namely, the cluster consists of self-crossings and segments between the crossings. We report on the measurements of distribution in the length of the segments, long-range correlations between directions of the links occupied by monopole currents, scaling properties of the segments characteristics. In Sect. III we first present data on the spectrum of the finite size clusters (subsection III.1) and then on their density (subsection III.2). In subsection III.3 we present results on the correlations of the monopole trajectories at large distances which are sensitive, at least in principle, to the glueball masses. Discussions of the observations made in this work are in Sect. IV. Finally, some details of simulations are given in the Appendix.

## II Segments of the percolating cluster

### II.1 Distribution in the segments length

The percolating cluster consists of segments (that is, trajectories between crossings) and crossings. The segments, in turn, are made from the links on the dual lattice. We will consider the length of the segment, $`l_{segm}`$ and the Euclidean distance between the end points of the segments, $`d`$.

In Fig. 1 we show the normalized distribution, $`N(l)/N`$, of the lengths of the segments; $`N(l)`$ is the number of the segments with the length $`l`$ in lattice units, $`N`$ is the total number of segments.

The oscillations of $`N(l)/N`$ are clearly seen, the number of trajectories with the even $`l`$ is systematically larger than that with the odd $`l`$. This effect is due to the large number of the closed trajectories (loops), which are segments with $`d=0`$ and $`l0`$. In Fig. 2 we show $`N(l)/N`$ only for not closed trajectories, $`d0`$.

The oscillations are substantially reduced compared with Fig. 1.

### II.2 Long-range correlation between links directions

Next, we turn to angular correlation between the links on the monopole trajectory. Let us call link $`C_0`$ belonging to the percolating cluster as the ”initial” one. Then one can measure the probability that the link $`C_l`$, connected to the link $`C_0`$ by the monopole trajectory of length $`l`$, has the same direction. Of course this probability is a decreasing function of $`l`$, but it occurs that there exists a ”long memory” of the initial direction. This fact is illustrated in Fig. 3 where the correlation of the direction of the initial link $`C_0`$ with the direction of the link $`C_l`$ is shown. The direction of $`C_l`$ can be the same as that of $`C_0`$, or opposite, or else (neither the same nor the opposite). We normalize all three probabilities in such a way that they are equal to unity for random walk if $`l0`$ (for $`l=0`$ the opposite direction is forbidden).

From Fig. 3 it is seen that even for links separated by 17 lattice steps the most probable direction is the same as the direction of the initial link. For smaller values of $`\beta `$ the correlations of the directions exist for not so large values of $`l`$. The deviation from unity of the probability to have the same direction falls off exponentially: $`P_{same}=1+A_se^{\mu _sla}`$. The probability to have the opposite direction behaves as follows: $`P_{opposit}=1A_oe^{\mu _ola}`$. These fits are shown by dashed lines in Fig. 3. The masses $`\mu _s`$ and $`\mu _o`$ coincide within the error bars and are independent of the lattice spacing, see Table 1.

### II.3 Scaling properties of the segments

It was already mentioned in the Introduction that the density of the percolating cluster scales, that is $`\rho _{perc}`$ has a well defined value in the continuum. Here we will present further data on the scaling properties of geometrical elements of the percolating cluster.

First, we observe that the average length of the segment of the monopole trajectory between crossings, $`l_{segm}`$, does not depend on the lattice spacing, see Fig. 4. Making fit by a constant for $`\beta >2.35`$ we obtain

$$l_{segm}=1.60(1)fm.$$

The way used to convert lattice results into physical units is explained in Appendix. The lattice spacing values $`a`$ for various $`\beta `$’s are also given in Appendix, see Table III.

In the average Euclidean distance between the crossings $`d`$ the violations of the scaling are more significant, see Fig. 5. However, the deviations from the scaling can be approximated for $`\beta >2.35`$ by linear in $`a`$ corrections. Then in the continuum limit $`a0`$, $`d`$ has a non-vanishing value ($`0.20fm`$), as it is seen from Fig. 5.

In Fig. 6 the number of crossing points per unit length of the monopole trajectory is shown by open circles. This number weakly depends on $`a`$ and seems to have the continuum limit ($`a0`$) $`0.3fm^1`$.

The situation changes if we exclude the closed loops of finite lattice length connected to the percolating cluster. The number of crossings reduces and is compatible with zero in the continuum limit if we exclude the loops of the length up to 8 (or up to larger length). This fact is illustrated in Fig. 6.

The continuum limit value $`N_{cross}/l_{perc}=0.3fm^1`$ for the full percolating cluster together with the data for the percolating cluster density (Fig. 8) corresponds to approximately 10 crossings in percolating cluster per hypercube $`1fm^4`$. To complete the picture of the structure of the percolating monopole cluster in the continuum limit the data for $`l_{segm}`$, $`d`$ and $`N_{cross}/l_{perc}`$ can be compared with the average monopole radius $`\rho _m0.05fm`$ and average inter-monopole distance $`R0.5fm`$ anatomy .

## III Finite clusters

### III.1 Length spectrum of finite clusters

In Fig. 7 we show the number of the finite clusters of the length $`l`$ (in the lattice units) vs. $`l`$. The dashed line in this figure is the fit of the data by function $`Const/l^\alpha `$.

It is important that $`1/l^3`$ behavior seems to be valid for all considered values of the lattice spacings as it can be seen from Table 2 where we show results of the fit for various $`\beta `$ values. As we discuss in Section IV.3 the dependence $`1/l^3`$ means that finite clusters develop in the four-dimensional space.

### III.2 Monopole densities

As we discussed in the introduction the magnetic monopole clusters fall into two different classes uvir ; hart : ”small”, or finite clusters which have finite size in lattice units, and ”large” clusters which percolate through the lattice. It was demonstrated hart  that the percolating cluster is responsible for the string tension. If the size of the lattice is large enough, there is a natural distinction between ”small” and ”large” clusters, since in each configuration there exists only one cluster which is much longer than all others and there is a gap in the cluster length distribution hart . If the lattice size is not large enough the percolating cluster decays into one or more large clusters plus several monopole currents which wind through the boundaries of the lattice mueller .

It occurs that the monopole density for the sum of these large clusters and the winding trajectories scales mueller . We simply call this collection of the monopole currents as percolating cluster.

In Fig. 8 we show the monopole density $`\rho _{perc}`$ as a function of the lattice spacing $`a`$. Our results for $`\rho _{perc}`$ obtained with higher statistics than in mueller  agree within error bars with values obtained in mueller  for coinciding values of $`\beta `$ on $`a`$ is rather weak and we fit it by a constant for $`\beta >2.35`$ The resulting density of the percolating monopoles in the continuum limit is

$$\rho _{perc}=7.70(8)fm^3.$$
(6)
This value is shown by the solid line in Fig. 8. It is in agreement with value obtained in mueller  within error bars.

The density of the finite clusters, on the other hand, is divergent for $`a0`$ and it can be fitted (dashed curve on Fig.8) as

$$\rho _{fin}=C_1+\frac{C_2}{a},$$
(7)
where

$$C_1=6.1(5)fm^3,C_2=1.55(4)fm^2.$$

The negative value of the constant $`C_1`$ means that the fit is not valid for large values of the lattice spacing.

We made also the local fits of $`\rho _{fin}`$ by function $`Aa^\alpha `$, where $`A`$ and $`\alpha `$ are fit parameters. The fit results for $`\alpha `$ are shown in Fig. 9. To get a point on this figure used for the fit three neighbor data points in Fig. 8. Corresponding fit intervals are depicted in Fig. 9 as X axes error bars. It can be seen from Fig. 9 that the values of $`\alpha `$ are changing from the value close to 3 down to the value approximately equal 1, i.e. at small $`a`$ this fit becomes consistent with the fit Eq.(7).

Another possible fit, suggested in Ref. hart , $`\rho _{fin}=Aa^\alpha `$ gives $`A`$ and $`\alpha `$ strongly dependent on $`a`$, see Fig. 9. However, at smallest $`a`$ we have the same $`1/a`$ behavior of the density.

### III.3 Correlator of the monopole currents

We turn now to discussion of correlation functions of the monopole currents. We will consider three definitions:

* $`G_1(xy)`$ is the probability that point $`x`$ and point $`y`$ are connected by a monopole trajectory.
* $`G_2(xy)`$ is the probability that point $`x`$ and point $`y`$ are connected by a monopole line belonging to the percolating cluster.
* $`G_3(xy)`$ is the probability that the monopole current crosses point $`x`$ and point $`y`$.

Since the small clusters have finite size in the lattice units (see Sect. III.1) the functions $`G_1(xy)`$ and $`G_2(xy)`$ coincide with each other for large $`d=|xy|`$. The asymptotic behavior for $`d=|xy|\mathrm{}`$ is similar for all correlation functions and was suggested first in IvPoPo  (see also reinhardt ):

$$G_k(xy)C_k+A_k\mathrm{exp}\left\{m_kd\right\},$$
(8)
where $`k=1,2,3`$ and $`m_k`$ are expected to be the mass of the lightest glueball. Moreover, since by definition (see (3)) the monopole density is the probability to find the monopole current on the given link, the constants $`C_k`$ are equal to the squares of the corresponding densities. It is obvious that $`C_{1,2}=\left(\rho _{perc}\right)^2`$ and $`C_3=\left(\rho _{perc}+\rho _{fin}\right)^2`$.

The results of the measurements are presented in Figs. 10, 11. In Fig. 10 we show the results of the rotational invariance check of the correlators $`G_k`$ at large enough $`d`$. Namely, the dependence of $`G_3(xy)`$ on the angle between 4-vector $`xy`$ and one (arbitrary) axis in polar coordinates is shown. For illustrative purposes, we use different normalizations of $`G_3`$ on different plots in Fig. 10. The data show nice rotational invariance. Note that the gap in the angles seen for large $`d`$ is due to the geometry of the lattice.

One of the most important properties of the correlators is their scaling. The data do indicate that the correlators do have the scaling behavior. This is illustrated in Fig. 11, where we show $`G_1(d)`$ in physical units for various $`\beta `$. Note that the dimension of $`\rho `$ is $`fm^3`$. Respectively, the dimension of the correlation functions $`G_k`$ is $`fm^6`$ and they are sensitive to possible scaling violations. Comparing the data for various volumes at $`\beta =2.4`$ we found the lattice volume independence of the correlators $`G_k(d)`$.

Finally, we determined the masses $`m_k`$ in Eq. (8). The fits made for the distance $`d>0.3`$ fm show that the masses $`m_k`$ do have scaling behavior, but we do not have the high precision in their determination:

$$m_{1,2}=1.87(16)GeV,m_3=2.4(5)GeV.$$

Due to large errors we can claim only that the results do not contradict the prediction that $`m_k`$ coincide with the mass of the lightest $`SU(2)`$ $`0^{++}`$ glueball as measured, in particular, in Refs. glueball .

## IV Percolation laws at short and large distances

In preceding sections we presented data on the phenomenology of the monopole clusters and have chosen to postpone comments on their implications till this section. It might worth emphasizing from the very beginning that not every piece of information that we got has a straightforward interpretation. Moreover, it is useful to distinguish between the properties of the clusters at relatively short distances, $`\mathrm{\Lambda }_{QCD}^1la`$ where $`l`$ is the length along the monopole trajectory, and at large distances $`l\mathrm{\Lambda }_{QCD}^1`$. As we shall see, it is the short distance behavior which is quite well understood while for larger distances one could rather speak of accumulation of the data. This might look paradoxical since traditionally the non-perturbative fluctuations are deprived of interesting short-distance behavior.

### IV.1  Monopoles as free particles at short distances

One of our main results is observation of scaling for a few observables :

* average distance between the self-crossings of the percolating cluster, as measured along the trajectory, scales (see Fig. 4);
* average Euclidean distance between the same points approximately scales, see Fig. 5;
* correlation length of direction alignment of links along the trajectory of the percolating cluster scales, see Table 1;
* correlators of the vacuum monopole currents scale, as illustrated in Fig. 11;
* we confirm scaling of the density of the percolating cluster, see Fig. 8.

Each of these checks is highly non-trivial since the most natural explanation of any scaling behavior is that the monopoles of size $`a`$ are physically meaningful. Thus, each example of the scaling behavior confirms, at least indirectly, fine tuning of the lattice monopoles.

Now, let us discuss briefly the geometry of the monopole cluster at various distances. At the very short distances, of order $`a`$, monopole trajectories can be thought of as (approximate) random walk. At these distances both the action and the entropy are of order $`1/a`$. Let us illustrate this statement by a well known equation (see, e.g., ambjorn ) for the propagating mass of a scalar particle $`m_{prop}^2`$:

$$m_{prop}^2a\left(M(a)\frac{\mathrm{ln}7}{a}\right),$$
(9)
$`M(a)`$ is the field theoretical mass of the same particle, discussed in the Introduction on the monopole example. A more general definition of $`M(a)`$ is that the action of the particle is

$$S_{cl}=M(a)l,$$

where $`l`$ is the length of a path connecting two points between which the particle propagates. Moreover, the $`\mathrm{ln}7`$ factor is due to the entropy. Eq. (9) is derived in the approximation of neglecting neighboring trajectories. Note a dramatic difference between $`M(a)`$ and $`m_{prop}`$. Namely, to ensure a finite propagating mass $`m_{prop}`$ one needs an ultraviolet divergent mass $`M(a)`$ (that is, the mass which determines the classical action).

Upon accounting for possible cancellations between the action and entropy one can introduce an ‘effective action’,

$$S_{eff}=l\mu ,(la),$$

which governs the behavior of the trajectories at distances much larger than $`a`$. The very fact that the density of the percolating cluster scales implies that at distances of order $`a`$ the divergent factors in the action and entropy cancel each other and $`\mu \mathrm{\Lambda }_{QCD}`$. This is just the fine tuning hypothesis which we heavily rely on.

Eq. (9) is useful for orientation in mass scales which came out from various fits to the data. In particular fits to the correlation of the links resulted in the attenuation factor of order $`\mathrm{exp}(\mu _sl)`$ where $`\mu _s`$ is about $`300MeV`$ and $`l`$ is measured along the trajectory, see Table 1. At first sight, $`300MeV`$ might look too low a mass. However, if we apply for estimates Eq. (9) we would obtain rather $`m_{prop}^2(300MeV)/a`$. In other words, we do have evidence that the most singular part of the action (as measured on the lattice, see the Introduction) is indeed cancelled by the entropy factor. But the value $`\mu _s300MeV`$ does not necessarily corresponds directly to a physical mass.

Once this cancellation is confirmed by the data, one can make predictions about short clusters with the length $`\mathrm{\Lambda }_{QCD}^1la`$ maxim . Namely, at such length one can neglect the mass factor and the spectrum of the closed loops should correspond to free particles. Then one can derive (see maxim  and references therein):

$$N(l)\frac{const}{l^3},R(l)\sqrt{l},$$
(10)
where $`R(l)`$ is the radius of the cluster of the length $`l`$. Let us emphasize that (10) is a consequence of free field theory in Euclidean space-time (actually, in $`d=4`$ Eq. (10) survives adding Coulomb like interaction as well maxim ). Eq. (10) is derived most naturally in the so called polymer representation, see, e.g., polyakov1 ; maxim .

Measurements of $`N(l)`$ and $`R(l)`$ were reported in hart  and are in perfect agreement with (10). Our analysis confirms (10) on larger statistics and for smaller values of $`a`$. In particular, the data on the $`l`$ distribution are summarized in the Table 2. We did not include the data on the $`R(l)`$ which we have but they do confirm (10) for all $`a`$ tested.

### IV.2 Percolating cluster

Existence of a single infinite cluster is typical for all percolating systems in the supercritical phase, see, e.g., grimmelt . There is a general theorem on the uniqueness of the percolating cluster at $`p>p_c`$ where $`p`$ is the probability to have a link ‘open’ (in our case, to belong to a monopole trajectory) and $`p_c`$ is the point of the phase transition. The percolating cluster is characterized by the probability $`\theta (p)`$ of a given link to belong to the cluster (in our case this probability is equal to the cluster density, $`\rho _{perc}`$). Generically, $`\theta (p)(pp_c)^\beta ,\beta >0`$. In our case,

$$\theta (p)(a\mathrm{\Lambda }_{QCD})^31.$$
(11)
Smallness $`\theta (p)`$ implies closeness to the point of phase transition to the percolation. In the limit $`a0`$ we hit exactly the point of phase transition. Expression (11) can be considered as formulation of the hypothesis on fine tuning of the monopoles in terms of the percolation theory.

Another general property of the percolating cluster is that its fractal dimension coincides with the dimension of the space:

$$D_{fr}^{perc}=d.$$
(12)
Which is trivially satisfied in our case, since the length of the percolating cluster is proportional to the whole volume, $`V_4`$.

However, to the best of our knowledge, the percolation theory is not powerful enough to predict or explain the data which we have on such characteristics as $`l_{segm},d`$.

It is worth emphasizing that geometrical characteristics of the percolating cluster at length of order $`l_{segm}`$ differ quite drastically from the corresponding characteristics of the finite clusters. First of all, the monopole trajectories are no longer random walks. Indeed, for a random walk we would have had:

$$l\frac{1}{2}\frac{d}{m_{prop}a},$$
(13)
where $`l`$ is the distance between two points measured along the trajectory and $`d`$ is Euclidean distance between the same points.

Our results on $`l_{segm}`$ and $`d`$ are given in Fig. 4 and 5, respectively. The data certainly rule out (13) for any finite $`m_{prop}`$. Moreover if we assume that it is only the leading power of $`1/a`$ which is cancelled in Eq. (9) then $`m_{prop}\sqrt{\mathrm{\Lambda }_{QCD}/a}`$. Our data are not consistent with (13) for such mass either.

Also, existence of correlation between directions of the links, (see Sect. II.2), could not be reconciled with the random walk.

### IV.3 Long-range correlation of finite clusters

Turn now to discussion of the data on the density of the short clusters, see Sect. III.2. Expectations for $`\rho _{fin}`$ are easy to formulate. Indeed, the spectrum $`N(l)l^3`$ implies that the average length of the finite clusters is saturated in the ultraviolet, i.e. at length of order $`a`$. Then the density of the finite clusters should not depend on $`\mathrm{\Lambda }_{QCD}`$ at all and on the dimensional grounds $`\rho _{fin}a^3`$.

Instead, we find that the most singular piece in the density $`\rho _{fin}\mathrm{\Lambda }_{QCD}^2a^1`$. The only interpretation of this striking observation is that even short, or finite clusters are not independent at large scales. Again, on dimensional grounds alone it is clear that the finite clusters are in fact associated with two-dimensional surfaces whose area scales.

At first sight, this conclusion looks too bizarre. However, in fact it is known to be true from independent measurements on the P-vortices. Namely, it was observed first in Ref. giedt  that monopoles are associated predominantly with P-vortices, whose total area scales. This conclusion was recently reinforced by measurements for smaller lattice spacings kovalenko .

Thus, we are invited to think about percolation of monopoles as a two-step process: monopoles percolate on the surface of dimension $`d=2`$ and the surface itself percolates in $`d=4`$ space-time. As far as we know, this type of percolation has not been studied theoretically at all. And we can, therefore, suggest only very preliminary considerations.

First the probability $`\theta (p)`$, see Eq. (11) is to be thought of as a product of two factors:

$$\theta (p)(a\mathrm{\Lambda }_{QCD})^3(a\mathrm{\Lambda }_{QCD})\frac{V_2}{V_4},$$
(14)
where the first factor is to be interpreted as a power of $`(p^{(2)}p_c^{(2)})`$ for percolation of monopoles on a surface while the second factor accounts for suppression of the ‘phase space’ available on a $`d=2`$ space, $`V_2`$, as compared to the whole lattice volume $`V_4`$. Note that for probabilities defined on the surface we still have according to (14)

$$(p^{(2)}p_c^{(2)})(a\mathrm{\Lambda }_{QCD})1,a0.$$

Another major change, is that finite clusters now occupy a finite fraction of the links on the surface. Furthermore, one could also speculate that for finite monopole clusters of length $`l\mathrm{\Lambda }_{QCD}^1`$ we would still have

$$Rl^{1/d}\sqrt{l},d=2$$

As the last remark, let us notice that the interplay between the monopole trajectories and P-vortices is more complicated than simply saying ‘monopoles belong to P-vortices’. Indeed, for short clusters the length spectrum is sensitive to the dimension of the space polyakov1 ; maxim :

$$N(l)\frac{1}{l^{d/2+1}},$$
(15)
and the spectrum $`l^3`$ corresponds to $`d=4`$. In this sense, at short distances the monopole trajectories are ‘primary’ objects. Thus, one can expect that for larger $`l`$ where the infrared behavior of the clusters sets in the spectrum is changed into $`1/l^2`$.

## V Conclusions

There is accumulated evidence that monopoles defined within the Maximal Abelian projection appear physical even if studied with resolution of the lattice spacing $`a`$. In particular, we have seen that there are amusingly simple scaling laws for various observables which depend only on the product $`(a\mathrm{\Lambda }_{QCD})`$ and look perfectly gauge-invariant. Thus, non-perturbative fluctuations in the lattice $`SU(2)`$ appear to have non-trivial structure in the ultraviolet as well as in the infrared.

What is specific for the field, is that phenomenology seems to be far ahead of theory. In particular, we have argued that there are indications that percolation of the monopoles is to be considered as combination of percolation of the trajectories over a surface $`d=2`$ and of percolation of the surface over the $`d=4`$ space. To the best of our knowledge, no theoretical framework is available to describe such a percolation.

## Acknowledgements

The authors are thankful to M. Müller-Preussker for contribution at the initial stage of this work. We are obliged to F.V. Gubarev, M.N. Chernodub, A.V. Kovalenko and S.N. Syritsyn for numerous fruitful discussions. V. G. B., P. Yu. B. and M. I. P. are partially supported by grants RFBR 02-02-17308, RFBR 01-02-17456, DFG-RFBR 436 RUS 113/739/0, INTAS-00-00111 and CRDF award RPI-2364-MO-02.

## Appendix

The list of configurations used in this work is given in Table 3. To fix the Maximal Abelian gauge we use the Simulating Annealing algorithm SA  and study 10 randomly generated gauge copies for each configuration to avoid the Gribov copy problem.

Most of the quantities in the paper are given in physical units (e.g. in fm). In order to express the lattice data in physical units we use the data of Refs. physunit  for the $`SU(2)`$ string tension and assume that $`\sqrt{\sigma }=440MeV`$. The corresponding values of the lattice spacing are given in Table 3.


text:
# Synchronous vs. asynchronous dynamics of diffusion-controlled reactions

## I Introduction

Recently, the modelling of diffusion-controlled reactions on lattices has attracted renewed interest due to synergies with some research areas like e.g. heterogeneous catalysis alban ; zhdan , trapping problems rice , spin models racz , game theory fell , population dynamics dick  or biological problems faum . An issue common to these systems is the important role played by the coexistence of different intrinsic time scales, the lattice characteristics (size and dimensionality) and many-body effects. The interplay of these ingredients may strongly affect the efficiency with which statistical processes such as front propagation on catalytic substrates, the spread of an infection, kink propagation in magnetic systems, or exciton trapping in photosynthetic cells take place. In order to obtain analytical insight for such systems, it is often necessary to make simplifying assumptions which nevertheless preserve their main generic features.

In this spirit, a prototypical diffusion-reaction system consists of two interacting walkers performing nearest neighbor jumps on a lattice fish . The walkers are assumed to react with each other whenever they meet at the same lattice site or attempt to exchange positions. Each of such “encounters” results in instantaneous reaction, the process therefore being diffusion-controlled. Regardless of the particular outcome of the reaction, a measure for its efficiency will be clearly given by the mean encounter time of both walkers.

In a pioneering work, Montroll investigated a simplified version of the problem in which one of the walkers is stationary, thereby playing the role of a fixed trap mont1 . In a recent article nickoz1 , Montroll’s results have been extended with the help of a Markovian method to account for the possibility of simultaneous displacement of the walkers. In the present work, the problem is further extended to the case in which the motion of the walkers consists of a random sequence of two different events: at each tick of the clock, a synchronous event takes place with probability $`p`$, i.e., both walkers hop simultaneously to randomly chosen nearest neighbor sites. Alternatively, an asynchronous event takes place with probability $`1p`$, i.e. one of the walkers performs a nearest neighbor jump while the other remains immobile. Thus, the parameter $`p`$ interpolates between the one walker plus trap case studied by Montroll ($`p=0`$) and the case of two simultaneously moving walkers ($`p=1`$). A central question we shall investigate here concerns the influence of the degree of synchronicity $`p`$ of the walkers and the size $`N`$ of the lattice in which they are embedded upon the reaction’s efficiency.

Our model may be of interest in several contexts. One can e.g. easily accomodate it to allow for events in which both walkers remain immobile in the original reference frame. Thus, three qualitative different joint events become possible for the two-walker system, which can e.g. be interpreted as resulting from a combination of two internal states for each of the walkers, namely a diffusing and an immobile state. Such two-state random walks weiss2 ; roer ; vanden  are frequently used to model chromatographic gand  or electrophoretic separation processes, in which the propagation of charged particles in an external field may be occasionally stopped by entanglement with the molecules of the substrate. Such systems, along with hopping models for transport and recombination of carriers in solids nool ; fister ; schmid , might provide additional motivation for considering a generalized version of the random walk.

A different approach suggests to regard the fluctuations in the diffusivity of the walkers as a result of random fluctuations of lattice sites switching between a conducting and a stopping state, whereby the translational invariance of the lattice is preserved on average. An alternative formulation of the problem in terms of fluctuating dichotomous barriers between sites also seems possible beni . Such models for dynamic random media are relevant for the description of several physiochemical processes like e.g. ligand diffusion in proteins elber  or proton migration in water tucker .

On the other hand, if one goes back to the original picture of two dissimilar walkers $`A`$ and $`B`$, the model may be regarded as a schematized starting point for describing the dynamics of exciton absorption in biological light-harvesting systems whit . In photosynthetic cells, a photon is absorbed by the pigment molecules (e.g. chorophyll) in the cell and may give rise to an excited energy state. The exciton hops by resonant energy transfer through a network or lattice of 200-500 pigment molecules (antenna system) and can be eventually trapped at a reaction center vang , which is usually considered to be immobile within the time scale for trapping (a few hundred picoseconds). The exciton energy is then used to trigger a series of redox processes in the chain of chemical reactions leading to the production of sugars and carbohydrates. Thus, one of the time-limiting steps for the production of oxygen is precisely the absorption of the exciton by the trapping center. If we allow for a certain mobility of the reaction center (thereby generalizing Montroll’s approach for exciton trapping), we can identify the latter with a slowly hopping walker, say the $`A`$ walker, while the propagating exciton would play the role of the $`B`$ walker. As long as the hopping rate of the $`A`$ walker remains small, the situation may be identified with the almost purely asynchronous case (small $`p`$); this is normally the case in in vivo light-harvesting systems. However, a modification of the physical properties of the antenna system so as to make exciton propagation slower might, at least in principle, lead to interesting antiresonance phenomena as observed in our model.

The paper is organized as follows. In section 2, we define the two-walker system in detail and show how it can be recast into an equivalent one-walker system with absorbing sites. In section 3, we report analytic work on the $`1d`$ case and identify the principal differences between the purely asynchronous ($`p=0`$), the purely synchronous case ($`p=1`$) and the mixed case ($`0<p<1`$). These results are also supported by Monte Carlo simulations. In Section 4, similar results are found for the $`2d`$ and the $`3d`$ cases by means of simulations In section 5, the results are compared with the predictions of the continuum approximation valid for large lattices. Finally, section 6 summarizes the main conclusions and discusses possible extensions.

## II Formulation of the problem: two-walker vs. one-walker picture

The starting point is a 1d periodic lattice with $`N`$ sites and discrete time dynamics (cf Fig. 1a). We place two walkers $`A`$ and $`B`$ on two distinct lattice sites and let them evolve at each time step as follows:

1. with probability $`p`$, both walkers hop simultaneously to randomly chosen nearest neighbor sites (synchronous event).
2. with probability $`1p`$, one of the walkers (no matter which one) remains at rest while the other performs a jump to a nearest neighbor site (asynchronous event).

The walkers are assumed to be unbiased, i.e., their jumps are symmetric. We additionally assume that their jump directions are completely uncorrelated. An encounter takes place when both walkers ‘land’ on the same site or attempt to exchange positions. Each encounter triggers instantaneously an irreversible reactive interaction, say the annihilation reaction $`A+B0`$. The encounter time can thus be regarded as the characteristic reaction time governing the diffusion-controlled two-particle annihilation.

In principle, the mean encounter time can be computed for a single initial configuration of the walkers; in many practical situations, though, one has little knowledge about the initial conditions. We shall therefore give preference to a definition of the mean encounter time which contains an additional coarse graining over all possible initial configurations. In the sequel, we shall denote this quantity by $`n`$. The smaller $`n`$, the more efficient the reaction will be.

Note that, due to the translational invariance of the lattice, the physical distinguishability of the walkers is irrelevant for the computation of $`n`$: in other words, it does not matter which of the walkers $`A`$ or $`B`$ hops more often, since $`n`$ depends only on the relative motion of both walkers, the latter being fully characterized by $`p`$. Therefore, we shall assume for simplicity and without loss of generality that $`A`$ and $`B`$ are physically identical walkers and drop the labels $`A`$ and $`B`$, as has been done in Fig. 1a.

For $`p=0`$, only reactions through simultaneous site occupancy are possible. As far as $`n`$ is concerned, this case is equivalent to the one treated by Montroll mont1 . However, if $`p>0`$, reaction through position exchange becomes possible. This technical obstacle makes the analytical treatment of the problem more difficult.

A first step to circumvent this difficulty is to take advantage of the translational invariance of the lattice and reformulate the problem in terms of a single walker. To do so, we must switch over to a new comoving reference frame in which one of the walkers remains stationary, thereby playing the role of a fixed perfect trap $`T`$ (Fig. 1b). Clearly, an asynchronous event in the original two-walker system is equivalent to nearest neighbor hopping in the single-walker system, while a synchronous event results either in a walker’s jump by two lattice sites (if both walkers hop in opposite directions in the original reference frame) or its remaining at rest (if they hop in the same direction). In this picture, reaction takes place any time the walker reaches or overreaches the trap $`T`$.

We can now take advantage of the simple lattice geometry and the fact that the trap is perfect to unfold the $`N`$-site lattice into an equivalent one with $`N+1`$ sites and two perfect traps $`T`$ sitting at each end site, as shown in Fig. 2 for $`N=7`$. This transformation does not of course affect the characteristics of the walk; a walker jump in anticlockwise direction will be equivalent to a jump to the left by the same number of sites in the transformed lattice. Next, we replace each trap $`T`$ by two (fictitious) reactive sites $`r`$ (Fig. 3). The random walk will instantaneously terminate when the walker lands on any of these $`r`$-sites. We shall therefore occasionally call the $`r`$-sites “absorbing” in what follows.

Let us assign the coordinates $`1`$ to $`N1`$ to each of the non-absorbing sites, as shown in Fig. 3. Additionally, we term the $`r`$-sites at the left end $`1`$ and $`0`$ and those at the right end $`N`$ and $`N+1`$. An attempt to overcome, say the left trap, triggered by a synchronous event will be equivalent to a jump from site $`1`$ to site $`1`$. In contrast, the walker’s landing on the trap can be realized either by an asynchronous event (jump from site $`1`$ to site $`0`$) or by a synchronous event (jump from site $`2`$ to site $`0`$). If the dynamics is purely asynchronous, only one $`r`$-site at each end will be needed, since jumps by two sites are not possible in this case. Thus, events involving position exchange can be dealt with by going over to a one-walker picture and introducing fictitious absorbing sites. The mean duration of the walk averaged over the initial positions $`1,\mathrm{},N1`$ of the walker plays the role of $`n`$ in the original two-walker system.

A standard approach for the calculation of $`n`$ is to formulate the problem for the restricted walk between the absorbing sites in terms of a conditional first passage problem for an unrestricted walk on an infinite $`1d`$ lattice weiss . The starting point is the Markovian master equation

$$P_{n+1}(j)=\frac{p}{4}\left[P_n(j2)+2P_n(j)+P_n(j+2)\right]+\frac{1p}{2}\left[P_n(j1)+P_n(j+1)\right],$$
(1)
where $`P_n(j)`$ is the probability to find the walker at a given site $`j`$ in the infinite lattice after $`n`$ time steps. The first term on the right hand side of eq. (1) is the contribution due to the synchronous events, by which the walker either remains at rest or it moves (here) two lattice sites either to the left or to the right. The second term describes jumps by one lattice site as a consequence of asynchronous events. The mean time to absorption $`n`$ can now be viewed as a first-passage time and computed via a generating function approach weiss . Specifically, $`n`$ can be expressed in terms of the generating function for the probability of the walker arriving for the first time at a given $`r`$-site after a given number of steps without having been previously in neither of the other three $`r`$-sites. Unfortunately, the resulting expressions for $`n`$ are not very transparent and their analytical dependence on $`N`$ and $`p`$ is not easy to determine.

Another possibility is to make use of the single step probabilities appearing as coefficients in the right hand side of eq. (1) to compute the transition probabilities between the states of the underlying absorbing Markov chain. In this approach, $`n`$ can be related to the row sums over the elements of the fundamental matrix for the transition to the absorbing state or, alternatively, to its smallest eigenvalue kem . The disadvantage of this method is that it requires the inversion of increasingly large matrices as $`N`$ becomes large.

Finally, the method we shall further develop here exploits the analogy of our random walk problem with the classical ruin problem studied by Feller fell . Even though this approach has the disadvantage of being difficult to generalize to higher dimensions, it provides an elegant solution for the 1d problem.

## III Connection with ruin problem and analytical solution

We first recall briefly the classical gambler’s ruin problem. Consider a single walker (in our setting, this would correspond to the limit p=0), whose position $`z`$ is viewed as the capital of a gambler playing against an adversary whose capital is $`Nz`$. At each time step, a trial is made, as a result of which the gambler wins or loses one euro. Thus, the gambler’s winning corresponds to a nearest neighbor jump of the walker to the right, while losing the trial corresponds to a jump to the left. The game goes on until the gambler’s capital is reduced to zero or increases to $`N`$ (absorption of the walker at sites $`0`$ or $`N`$). One is interested in the mean duration of the game $`n_z`$ when the gambler starts with a given capital $`z`$. This quantity can be shown to be finite as long as $`N<\mathrm{}`$ and obeys the following difference equation fell :

$$n_z=\frac{1}{2}n_{z+1}+\frac{1}{2}n_{z1}+1,0<z<N$$
(2)
with the boundary conditions

$$n_0=0\text{and}n_N=0.$$
(3)
Eq. (2) states that the walker has no memory of where it has been at earlier time steps; if it is initially at site $`z`$, it will either jump to site $`z+1`$ or to site $`z1`$ with probability $`1/2`$. Once at any of these sites, the walker will continue his walk without remembering its previous position $`z`$. It is as though the walker started a new walk from $`z+1`$ or $`z1`$ with equal probability, except that the expected value of the mean time to absorption must be increased by one unit walsh . The boundary conditions (3) reflect the fact that, if the walker is initially placed at a $`r`$-site, it is immediately absorbed. In the original two-walker system, this is equivalent to placing both walkers at the same site.

Eqs. (2)-(3) can be solved by standard methods, e.g. by writing the general solution as the sum of the general solution of the corresponding homogeneous equation plus a particular solution jord . One obtains in this way

$$n_z=z(Nz),$$
(4)
and

$$n=\frac{1}{N}\underset{z=1}{\overset{N1}{}}n_z=\frac{N(N+1)}{6}.$$
(5)
The opposite limit of the above ($`p=1`$), corresponding in our setting to a purely synchronous motion of two walkers, is somewhat less standard. In this case, only jumps by zero or two lattice sites may occur. Depending on its initial position, the walker may land on any of the four $`r`$-sites depicted in Fig. 3. Therefore, two additional boundary conditions for the absorbing sites $`1`$ and $`N+1`$ are needed.

In the gambler’s jargon, a jump by two sites means doubling the stake of each trial, i.e., the player wins or loses two euros with probability 1/4. Besides, a tie occurs with probability 1/2 (no jump). The game terminates when either the player or his adversary reaches or overreaches a capital of $`N`$ euros. Now, the difference equation for the duration of the game is a fourth-order one:

$$n_z=\frac{1}{4}n_{z+2}+\frac{1}{2}n_z+\frac{1}{4}n_{z2}+1,0<z<N.$$
(6)
The solution of this equation requires four boundary conditions, namely

$$n_1=0,n_0=0,n_N=0\text{and}n_{N+1}=0.$$
(7)
A particular solution of (6) is given by $`\frac{1}{2}z^2`$. The characteristic equation of the homogeneous equation has two double roots $`\lambda _{1,2}=\pm 1`$. Therefore, the full inhomogeneous solution reads

$$n_z=\frac{1}{2}z^2+c_1+c_2z+(c_3+c_4z)(1)^z,$$
(8)
Substituting into eq. (6) one obtains a system of equations for $`c_1`$ to $`c_4`$ whose solution yields

$`n_z`$ $`=\left(z(N1z)+N\right)/2`$  for $`z`$ odd, $`N`$ odd,
$`=z(N+1z)/2`$ for $`z`$ even, $`N`$ odd, (9)
$`=\left(z(Nz)+(N+1)\right)/2`$ for $`z`$ odd, $`N`$ even,
$`=z(Nz)/2`$ for $`z`$ even, $`N`$ even.
The spatially averaged time to absorption $`n`$ is again easily computed. One obtains

$`n`$ $`=N(N+1)(N+2)/12(N1)`$ for $`N`$ even,
$`=(N+1)(N+3)/12`$ for $`N`$ odd. (10)
For notational convenience, let us rename the value of $`n`$ obtained from eqs. (III) as $`n^{(1)}`$, and the corresponding result of Montroll for the one-walker problem (eq. (5)) as $`n^{(0)}`$. A comparison between $`n^{(0)}`$ and $`n^{(1)}`$ is shown in Fig. 4. As expected, the synchronous case becomes more efficient as soon as the lattice gets sufficiently large ($`N5`$). In the limit of a very large lattice, it is asymptotically twice as efficient as the asynchronous case \[see Fig. 5)\].

It is also worth comparing the $`z`$-distribution of the encounter time for the case $`p=1`$ (eq. (III)), with the earlier known result for $`p=0`$ (eq. (4)). For odd $`N`$, the spatial profile in $`z`$ displays some qualitative similarities in both cases (cf. Figs. 6a and 6b). The highest value of $`n_z`$ is attained at those sites $`z`$ that maximize the distance $`d_z=\text{min}(z,Nz)`$ from the closest absorbing site. For $`p=0`$, the encounter time always increases with increasing $`d_z`$, and it becomes maximum when $`z=(N\pm 1)/2`$. In contrast, the behavior is no longer strictly monotonic for $`p=1`$, since $`n_z`$ either increases or remains constant as $`d_z`$ becomes larger, giving rise to a “staircase” profile for high enough values of $`N`$ (Fig. 6b).

For even lattices, there are more marked differences between the cases $`p=0`$ and $`p=1`$. For $`p=0`$, the discrete spatial profile resembles an inverted parabola like in the odd lattice case, the maximum now being located at $`z=N/2`$. However, if $`p=1`$, one observes a series of alternating valleys and peaks in the distribution (cf Figs. 7a and 7b). There are two subcases here: if $`N`$ is divisible by $`4`$, the highest values of $`n_z`$ are attained for $`z=(N\pm 2)/2`$ (Fig. 7a). Otherwise, the maximum value corresponds again to $`z=N/2`$ (Fig. 7b).

Note that for $`p=0`$ and even $`N`$ the time to absorption when the walker is started at sites $`1`$ or $`N1`$ is smaller than in the purely synchronous case. An intuitive argument points to the fact that, in the former case, the walker is absorbed one out of two times after the first time step, while absorption takes place only one out of four times if $`p=1`$. However, this argument should be taken with care, since it fails for odd $`N`$. Besides, as we shall see later on, the minimum of $`n_1=n_{N1}`$ corresponds to a process with $`p>0`$.

We finally turn to the general case of a mixed walk. Assume that a given event is synchronous with probability $`0<p<1`$ and asynchronous with probability $`1p`$. The difference equation for $`n_z`$ now reads:

$$n_z=\frac{p}{4}n_{z+2}+\frac{1p}{2}n_{z+1}+\frac{p}{2}n_z+\frac{1p}{2}n_{z1}+\frac{p}{4}n_{z2}+1,0<z<N.$$
(11)
The boundary conditions are again given by (7). A particular solution of eq. (11) is given by $`z^2/(1+p)`$. The roots of the underlying characteristic equation are

$$\lambda _{1,2}=1,\lambda _{3,4}=\frac{1\pm \sqrt{1p^2}}{p}.$$
(12)
The full solution reads

$$n_z=z^2/(1+p)+c_1+c_2z+c_3\lambda _3^z+c_4\lambda _3^z,$$
(13)
where the constants $`c_1c_4`$ are now given by the linear system

$`{\displaystyle \frac{1}{1+p}}+c_1c_2+\lambda _3^1c_3+\lambda _4^1c_4=0,`$ (14a)
$`c_1+c_3+c_4=0,`$ (14b)
$`{\displaystyle \frac{N^2}{1+p}}+c_1+Nc_2+\lambda _3^Nc_3+\lambda _4^Nc_4=0,`$ (14c)
$`{\displaystyle \frac{(N+1)^2}{1+p}}+c_1+(N+1)c_2+\lambda _3^{N+1}c_3+\lambda _4^{N+1}c_4=0.`$ (14d)
We prefer to omit the rather lengthy analytic form of the coefficients $`c_1`$-$`c_4`$. From the above equations, $`n_z`$ and $`n`$ can be explicitly computed. For $`n`$ one has an expression of the form

$`n`$ $`=`$ $`{\displaystyle \frac{N(2N1)}{6(1+p)}}+c_1(N,p)+{\displaystyle \frac{Nc_2(N,p)}{2}}`$ (15)
$`+{\displaystyle \frac{1}{N1}}\left(c_3(N,p)[\lambda _3(p)]^{N1}+c_4(N,p)[\lambda _4(p)]^{N1}\right).`$
Table 1 displays the analytical expressions of $`n`$ as a function of $`p`$ for increasing values of $`N`$. These are rational functions of $`p`$ whose complexity increases with $`N`$. For $`N4`$ and small $`p`$, one can obtain a first-order approximation to the exact solution by Taylor-expanding eq. (15):

$$n=n^{(0)}\left[\mathrm{\hspace{0.17em}1}+\frac{N3}{N}p+𝒪\left(p^2\right)\right].$$
(16)
Figs. 8a and 8b contain several numerical plots of $`n`$, as a function of $`p`$ for different values of $`N`$. For $`N=2`$, the purely asynchronous case is more effective than the purely synchronous case <sup>1</sup><sup>1</sup>1Obviously, the purely asynchronous case has a maximum efficiency in this case, since the walker is always trapped after the first step.. As synchronicity is turned on, a monotonic increase of $`n`$ as a function of $`p`$ is observed. For $`N=3`$, $`n`$ does not depend on $`p`$. At each individual time step, the probability of reaction is always $`1/2`$, regardless of whether both walkers hop or only one of them. In a sense, the cases $`N=2`$ and $`N=3`$ are non-generic, since they only involve a single symmetry-distinct initial condition. For $`N4`$, a new parity effect appears: if $`N`$ is even, the most effective process is observed for an intermediate value of $`p`$ associated to a minimum of the function $`n`$ in the physically acceptable $`p`$-interval $`[0,1]`$; in contrast, for odd lattices, $`n`$ is a monotonically decreasing function of $`p`$, i.e., the most effective process is always the purely synchronous one. For $`N=4`$, any intermediate value of $`p`$ makes hopping more effective than in both limiting cases, and a minimum of $`n`$ is obtained for $`p_{min}=2/3`$. For higher, even values of $`N`$, the minimum is rapidly shifted to the right ($`p_{min}0.86`$ for $`N=6`$) and the $`p`$-interval for which processes are more efficient than the purely synchronous case shrinks dramatically. For large $`N`$, $`p_{min}`$ gets arbitrarily close to $`1`$ (see Table 2).

A series of Monte Carlo simulations for the periodic two-walker system has been carried out to confirm our analytic results based on the one-walker description. For two different lattice sizes, namely $`N=7`$ and $`N=8`$, we have performed a series of statistical runs to compute $`n`$, each run thereby comprising a whole set of symmetric-distinct nonreactive configurations. Due to the large variability $`n^2n^2`$ characteristic of first-passage problems, a relatively high number of runs was needed ($`10^6`$) to obtain an accurate value for $`n`$. The $`n`$-values obtained from simulations (with an accuracy of three significant digits) are listed in Tables 3 and 4 for the cases $`N=7`$ and $`N=8`$, respectively. The agreement with the analytical predictions is good, the maximum observed deviation is off the theoretical value by about $`1\%`$ only.

In order to obtain additional insight in the even-odd transition mechanism, we have studied qualitatively $`n_z`$ as a function of $`p`$ and $`N`$ for each single initial position $`z`$ of the walker. Let us characterize each site $`z`$ by its distance $`d_z`$ to the closest $`r`$-site. For not too large values of $`p`$, the behavior of the encounter time is roughly the same for all $`z`$ values regardless of the parity of $`N`$, i.e., a decrease of the encounter time is observed (Figs. 9a and 9b). However, the qualitative $`p`$-dependence of $`n_z`$ in the large $`p`$ limit becomes different for initial positions with even or odd values of $`d_z`$: for even values of $`N`$ and initial positions with odd $`d_z`$, $`n_z`$ begins to increase sharply, as a result of which a minimum of the curve is observed (cf Fig. 9a). Even though the contribution to the global efficiency $`n`$ arising from the $`(N/2)1`$ sites with even $`d_z`$ decreases strongly in this regime, this effect is overcome by the increase of the contribution yielded by the $`N/2`$ sites with odd $`d_z`$, thus giving rise to a net increase in $`n`$. In contrast, the sensitivity to the initial condition is less systematic and less important for odd values of $`N`$ (Fig. 9b): again, a monotonic decrease is observed for even $`d_z`$. Even though minima are still observed for $`d_z=1`$, for all other odd values of $`d_z`$, they are either absent <sup>2</sup><sup>2</sup>2This is e.g. the case for $`N=7`$ and $`d_z=3`$ (not shown here). or very flat (see curve for $`d_z=3`$ in Fig. 9b).

## IV Monte Carlo results in two and three dimensions

To complement the analytic $`1d`$ results for the dependence of $`n`$ on the value of the parameter $`p`$, we have also investigated how $`n`$ depends on $`p`$ for dimensions $`2d`$ and $`3d`$ using Monte Carlo simulations. Figures 10 and 11 show the dependence of $`n`$ with respect to $`p`$. As shown above for $`1d`$, there is distinctively different behavior based on whether $`N`$ is even or odd. In higher dimensions this even-odd effect becomes much more pronounced. It is clear that when $`N`$ is an even number, the maximal efficiency is attained at some intermediate value between 0 and 1. Also, one notices that for the $`6\times 6`$ square lattice and for all $`3d`$ even lattices which we studied $`n^{(0)}<n^{(1)}`$. This is a surprising result which we do not see in the $`1d`$ case for lattices with $`N>4`$. This suggests the existence of a crossover effect in the large size limit when switching from $`2d`$ to $`3d`$ lattices, implying that in the former case the two simultaneously moving walkers are more effective than two asynchronous walkers, while in $`3d`$ the opposite holds. As in the $`1d`$ case, the value of $`p_{min}`$ in $`2d`$ and $`3d`$ tends toward 1 in the limit of large lattice size, but it is interesting to note that in the largest $`3d`$ lattice which we studied ($`N=1000`$), the difference between $`p=0.999`$ and $`p=1`$ is $`600,`$ which is about $`30\%`$ of $`n^{(1)}`$ in that case. This shows that for even lattices, a minute amount of asynchronicity allows for a much greater efficiency. Similar arguments to those given above for $`1d`$ exist for in $`2d`$ and $`3d`$ when attempting to determine why such even-odd behavior arises, and further analysis of this striking behavior is given in Ref. bentz .

## V Comparison with continuum approximation

It is instructive to compare the above results in $`1d`$ with the continuum approximation valid for large $`N`$. To do so, consider the one-walker system with absorbing sites and a fixed lattice length $`L=N\mathrm{\Delta }x`$, where $`\mathrm{\Delta }x`$ is the intersite distance (lattice constant). According to this definition, $`L`$ is the distance between the inmost $`r`$-sites. The walker’s distance to site $`0`$ is $`x=z\mathrm{\Delta }x`$. We perform the continuum limit by letting $`\mathrm{\Delta }x`$ and $`\mathrm{\Delta }t`$ simultaneously go to zero under the additional requirement that the diffusive combination $`(\mathrm{\Delta }x)^2/\mathrm{\Delta }t`$ tend to a finite constant. Since $`L`$ is fixed, this implies that one lets the number of sites $`N`$ go to infinity while the spatial and temporal resolutions $`\mathrm{\Delta }x`$ and $`\mathrm{\Delta }t`$ are scaled as $`1/N`$ and $`1/N^2`$, respectively. Let us next replace $`n_z`$ by a function $`n_x`$ varying smoothly in the space interval $`[0,L]`$. The mean elapsed time $`t_x`$ to absorption will then simply be the mean number of steps $`n_x`$ times the time unit $`\mathrm{\Delta }t`$. From eq. (11), we have:

$$\frac{p}{4}\left[t_{x+2\mathrm{\Delta }x}2t_x+t_{x2\mathrm{\Delta }x}\right]+\frac{1p}{2}\left[t_{x+\mathrm{\Delta }x}2t_x+t_{x\mathrm{\Delta }x}\right]+\mathrm{\Delta }t=0$$
(17)
with the boundary conditions <sup>3</sup><sup>3</sup>3In the case $`p=0`$, one only has two boundary conditions, namely $`t_0=t_L=0`$.

$$t_{\mathrm{\Delta }x}=t_x=t_L=t_{L+\mathrm{\Delta }x}=0.$$
(18)
We now divide eq. (17) by $`\mathrm{\Delta }t`$ and expand the expressions in the brackets in $`\mathrm{\Delta }x`$. Taking the diffusive limit in the resulting equation yields

$$D\frac{d^2t_x}{dx^2}=1,$$
(19)
where the relative diffusion coefficient $`D`$ is given by

$$D=(1+p)\underset{\mathrm{\Delta }x,\mathrm{\Delta }t0}{lim}\frac{(\mathrm{\Delta }x)^2}{\mathrm{\Delta }t}=(1+p)D^{(0)}.$$
(20)
In the rightmost equation, $`D^{(0)}`$ is the value of the diffusion coefficient for $`p=0`$. For $`p>0`$, $`D^{(0)}`$ is increased by the prefactor $`1+p`$, i.e., the variance of the single-step probabilities of the random walk. In this limit, the four boundary conditions (18) coalesce into two distinct ones, namely $`t_0=0`$ and $`t_L=0`$. The solution of (19) which fulfils these boundary conditions is

$$t_x=\frac{x(Lx)}{2D}.$$
(21)
The spatially averaged reaction time $`t`$ is obtained by integrating over $`x`$:

$$t=\frac{1}{L}_0^Lt_x𝑑x=\frac{L^2}{12D}.$$
(22)
As expected, $`t`$ is proportional to the squared lattice length and inversely proportional to the relative diffusion coefficient $`D`$. For the special cases $`p=0`$ and $`p=1`$, this result is recovered by directly taking the diffusive limit in the discrete solutions (5) and (III). The relation (22) leads to the asymptotic law

$$\frac{t^{(0)}}{t}=\underset{N\mathrm{}}{lim}\frac{n^{(0)}}{n}=1+p,$$
(23)
where $`t^{(0)}=L^2/(12D^{(0)})`$ and $`n`$ is given by eq. (15). Eq. (23) shows that in the continuum limit, the efficiency of the reaction increases linearly with $`p`$. The continuum approximation applies for sufficiently large $`N`$, namely when the typical displacement of the walker $`\mathrm{\Delta }l\sqrt{1+p}\mathrm{\Delta }x`$ at each time step is small compared to the lattice length $`N\mathrm{\Delta }x`$. If this condition is not fulfilled, the approximation gets significantly worse at large values of $`p`$ (cf Fig. 12).

A generalization of (19) for higher order moments $`t^j_x`$ can be obtained by writing down the difference equations for the discrete quantities $`n^j_z`$ and taking the diffusive limit thereof. One then gets the coupled set of equations

$$D\frac{d^2t^{j+1}_x}{dx^2}=(j+1)t^j_x,$$
(24)
Eqs. (24), not to be further dealt with here, are well known from the theory of first-passage problems, where they are usually obtained from the adjoint Fokker-Planck equation for the underlying diffusion process weiss ; cox . Again, deviations from the dynamics dictated by (24) are expected for small lattices.

## VI Conclusions and outlook

We have seen that the $`1d`$ problem of computing the mean reaction time between two diffusing co-reactants can be reduced to a trapping problem for a single walker. The latter can be viewed as a generalized ruin problem, the duration of the game plays thereby the role of the mean time to absorption.

In the diffusive limit, equivalent to the limit $`N\mathrm{}`$ if the lattice length $`L`$ is held fixed, the reaction efficiency increases linearly with $`p`$, but important deviations are observed for not too large values of $`N`$. Beyond the crossover size $`N=4`$, a new parity effect is observed. For odd values of $`N`$, the reaction time still increases monotonically (but no longer linearly) with $`p`$, while for even values of $`N`$, the efficiency is optimized for an intermediate value $`p_{min}<1`$. In higher dimensions, this parity effect is even more pronounced, i.e. for even lattices there is a drastic increase in efficiency when a tiny amount of synchronicity is introduced. In contrast to the $`1d`$ case, the effect is enhanced with increasing $`N`$.

Let us briefly comment on the $`1d`$ results from the perspective of the ruin problem. Assume that one of the gamblers is successively given an starting capital of $`1,2,\mathrm{},N1`$ euros at each round, while his adversary gets $`N1,N2,\mathrm{},1`$ euros. Let us further suppose that the gamblers can choose between two kinds of trials: the stake for the first one is one euro and there is no tie. In the second trial, either nobody wins or one of the gamblers wins two euros. According to our results, for an odd capital $`N`$, the gamblers minimize the average playing time if they always make a two-euro bet. However, if $`N`$ is even, they should make a small amount of one-euro bets in order to finish the round as soon as possible.

Our work can be extended in many different ways. Perhaps the most straightforward one is the detailed characterization of the whole distribution $`P(n)`$ in terms of $`p`$ and $`N`$. A generalization of results such as equation (23) to higher dimensional integral and fractal lattices is also of interest, since it may further clarify the role of dimensionality and the lattice coordination number. According to our results, in one-dimension two synchronously moving walkers are asymptotically twice as efficient as when they hop one after the other. In a 2d square planar lattice, Kozak et al. have shown that the purely synchronous case is $`\sqrt{2}`$ times more efficient than the purely asynchronous one in the large $`N`$ limit for odd lattices nickoz1 . The question is whether or not the relative efficiency of both processes in lattices with fractal dimension $`1<d_f<2`$ lies between $`2`$ and $`\sqrt{2}`$. Preliminary calculations on a Sierpinski gasket (with fractal dimension $`d_f=1.585`$) seem to indicate an asymptotic relative efficiency higher than $`2`$ in this case, despite the fact that the lattice has (up to the three vertex sites) the same coordination number as a 2d square planar lattice bentz . The reason for this may be the important role played by the specific form of the lattice boundaries, even in the limit of a large lattice. This may motivate the study of boundary conditions other than periodic ones for the two-walker system. However, the analytical treatment of this case is considerably harder, at least in the framework of the method of difference equations, since the lattice is no longer translationally invariant.

As a further extension of our work, one can also consider more complex reactive schemes nickoz2  involving more than two walkers to study the combined effect of synchronicity and many-body effects.

###### Acknowledgements.

We are indebted to K. Karamanos, M. Plapp, A. Provata and F. Vikas for fruitful discussions. The authors gratefully acknowledge financial support from the Université Libre de Bruxelles. This work is also supported by a NATO Cooperative linkage grant PST.CLG.977780 and by the European Space Agency under contract number 90043. One of us (JLB) acknowledges support from the Arthur P. Hellwig Award.

### Figure Captions

Figure 1: a) two-walker system on a seven-site periodic lattice (walkers represented by black circles). b) Equivalent one-walker system with trap. For convenience, the walker labels $`A`$ and $`B`$ in Fig. 1a have been left out. The arrows in Fig. 1a indicate that both walkers perform a synchronous step. In Fig. 1b this corresponds to a two-site jump of the walker.

Figure 2: Lattice transformation for the one-walker system displayed in Fig. 1b.

Figure 3: Replacement of the traps $`T`$ at each end site of the transformed lattice by two absorbing sites $`r.`$

Figure 4: Mean encounter time as a function of the lattice size for the purely synchronous case $`n=n^{(1)}`$ (circles) and the purely asynchronous case $`n=n^{(0)}`$ (crosses).

Figure 5: Ratio $`n^{(0)}/n^{(1)}`$ as a function of $`N`$. Note that the value of the ratio is the same for two consecutive odd and even values of $`N`$. The inset displays the behaviour for small $`N`$.

Figure 6: $`z`$-distribution of the encounter time in the cases $`p=0`$ (dashed lines) and $`p=1`$ (continuous lines) for a) $`N=3,5,7`$ and b) $`N=9,11`$. For $`N=3`$ both cases display the same flat distribution.

Figure 7: $`z`$-distribution of the encounter time for $`p=0`$ (dashed lines) and $`p=1`$ (continuous lines) for a) $`N=6,10`$ and b) $`N=4,8.`$

Figure 8: Mean encounter time as a function of $`p`$ for a) $`N=2,\mathrm{},5`$ and b) $`N=6,\mathrm{},9.`$

Figure 9: a) $`p`$-dependence of $`n_z`$ for all possible values of $`d_z`$ on a lattice with a) 10 sites b) $`9`$ sites.

Figure 10: Mean encounter time as a function of $`p`$ for two walkers on a $`2d`$ square planar lattice with periodic boundary conditions.

Figure 11: Mean encounter time as a function of $`p`$ for two walkers on a $`3d`$ cubic lattice with periodic boundary conditions.

Figure 12: Ratio $`n^{(0)}/n`$ as a function of $`p`$.