format: text

text:
# Properties of hadronic final states in diffractive deep inelastic 𝑒𝑝 scattering at HERA

## 1 Introduction

A class of deep inelastic scattering (DIS) events has been observed at HERA that have the characteristics of diffractive interactions. These events have a large rapidity gap  between the recoil-proton system and the produced hadronic system, and a small momentum transfer to the proton . The events can be pictured in terms of the $`t`$-channel exchange of an object that carries the quantum numbers of the vacuum, called the Pomeron ($`\mathrm{IP}`$), see Fig. 1(a). However, the nature of the Pomeron in DIS is at present far from clear. Measurements by the H1  and ZEUS  collaborations have shown that, in QCD-inspired models of the diffractive process, the Pomeron can be described as an object whose partonic composition is dominated by gluons. Alternatively, the diffractive process can be described by the dissociation of the virtual photon into a $`q\overline{q}`$ or $`q\overline{q}g`$ final state that interacts with the proton by the exchange of a gluon ladder .

The study of the hadronic final state in $`e^+e^{}`$ annihilation  has been a powerful tool in gaining information about the underlying partonic state. Similarly, the study of the partonic content of the hadronic final state in diffraction is a natural way to explore the dynamics of diffraction.

In this paper, a study is reported of the hadronic system, $`X`$, produced in the DIS process $`epeXp`$, where the diffractively scattered proton stays intact. The proton was detected and its three-momentum measured in the leading proton spectrometer (LPS) . Diffractive events are defined, for the purpose of this paper, as those events which contain a proton with more than $`97\%`$ of the initial proton beam energy. Previous results on hadronic final states in diffractive events at HERA have been obtained with the requirement of a large rapidity gap between the observed hadronic system and the scattered proton . The results obtained with rapidity-gap events were either defined in a reduced phase space by imposing a cut in rapidity , or Monte Carlo simulated events were used to extrapolate the characteristics of diffractive events over the areas of phase space removed by the rapidity-gap cut . By using the scattered proton to tag diffractive events, there is no need to rely on Monte Carlo generators to model correctly the part of the final state removed by the rapidity cuts, and the full angular coverage of the central detector can be used.

The properties of the hadronic system, $`X`$, were studied in terms of global event-shape variables such as thrust and sphericity in the center-of-mass (CMS) frame of $`X`$. This is analogous to the studies of global event-shape variables in $`e^+e^{}`$ annihilation  as a function of the CMS energy and to the analysis that led to the interpretation of three-jet events in terms of gluon bremsstrahlung . In addition to global event-shape variables, the properties of the diffractive events were described in terms of inclusive distributions such as energy flow.

## 2 Experimental setup

The data were recorded in 1997 with the ZEUS detector and correspond to an integrated luminosity of $`13.8\pm 0.3`$ pb-1. A detailed description of the ZEUS detector can be found elsewhere . A brief outline of the main detector components most relevant for this analysis is given below. Charged particles are tracked by the central tracking detector (CTD) , which operates in a magnetic field of 1.43 T provided by a thin superconducting coil. The CTD consists of 72 cylindrical drift chamber layers, organized in 9 superlayers covering the polar angle111The ZEUS coordinate system is a right-handed Cartesian system, with the $`Z`$ axis pointing in the proton beam direction, referred to as the “forward direction”, and the $`X`$ axis pointing left towards the center of HERA. The coordinate origin is at the nominal interaction point. The pseudorapidity is defined as $`\eta =\mathrm{ln}\left(\mathrm{tan}\frac{\theta }{2}\right)`$, where the polar angle, $`\theta `$, is measured with respect to the proton beam direction. region $`15^{}<\theta <164^{}`$. The transverse-momentum resolution for full-length tracks is $`\sigma (p_T)/p_T=0.0058p_T0.00650.0014/p_T`$, with $`p_T`$ in GeV. The high-resolution uranium-scintillator calorimeter (CAL)  consists of three parts: the forward (FCAL), the barrel (BCAL) and the rear (RCAL) calorimeters. Each part is subdivided transversely into towers and longitudinally into one electromagnetic section (EMC) and either one (in RCAL) or two (in BCAL and FCAL) hadronic sections (HAC). The smallest subdivision of the calorimeter is called a cell. The CAL energy resolutions, as measured under test beam conditions, are $`\sigma (E)/E=0.18/\sqrt{E}`$ for electrons and $`\sigma (E)/E=0.35/\sqrt{E}`$ for hadrons ($`E`$ in GeV).

The LPS  detects charged particles scattered at small angles and carrying a substantial fraction of the incoming proton momentum; these particles remain in the beampipe and their trajectory is measured by a system of position-sensitive silicon micro-strip detectors very close to the proton beam. The track deflection induced by the magnets in the proton beamline is used for the momentum analysis of the scattered proton. The LPS consists of six detector stations, S1 to S6, placed along the beamline in the direction of the outgoing protons, at $`Z=23.8`$ m, 40.3 m, 44.5 m, 63.0 m, 81.2 m and 90.0 m from the interaction point, respectively. In this analysis, only the stations S4 to S6 were used. These stations consist of two halves, each equipped with an assembly of six parallel planes of silicon micro-strip detectors which can be inserted to a position near the proton beam. Each detector plane has an elliptical cutout which follows the profile of the 10$`\sigma `$ envelope of the beam, where $`\sigma `$ is the standard deviation of the spatial distribution of the beam in the transverse plane. The acceptance of stations S4 to S6 of the LPS for protons close to the beam energy is a few per cent and approximately uniform for $`0.075<|t|<0.35`$ GeV2. The LPS longitudinal momentum resolution is $`\sigma (p_Z)/p_Z=5\times 10^6p_Z`$ ($`p_Z`$ in GeV). The transverse momentum resolutions, dominated by the proton beam emittance, are $`\sigma _{p_X}=35`$ MeV and $`\sigma _{p_Y}=90`$ MeV.

## 3 Kinematic variables and event reconstruction

### 3.1 Kinematic variables

The event kinematics of DIS processes can be described by the negative squared four-momentum transfer at the lepton vertex, $`Q^2=q^2=(kk^{^{}})^2`$ ($`k`$ and $`k^{^{}}`$ denote the four-momenta of the initial- and final-state positrons, respectively), and the Bjorken scaling variable, $`x=Q^2/(2Pq)`$, where $`P`$ is the four-momentum of the proton. The fraction of the energy transferred to the proton in its rest frame, $`y`$, is related to these two variables by $`y=Pq/(Pk)Q^2/xs`$, where $`\sqrt{s}`$ is the positron-proton CMS energy. The CMS energy of the $`\gamma ^{}p`$ system, $`W`$, is given by $`W^2=(q+P)^2=Q^2(1x)/x+M_p^2`$, where $`M_p`$ denotes the proton mass. Unless otherwise noted, the double angle method  is used to reconstruct the kinematic variables from the measured quantities.

Two more kinematic variables are needed to describe a diffractive DIS event, $`e(k)+p(P)e^{}(k^{})+X+p^{}(P^{})`$, where the scattered proton, $`p^{}`$, has four-momentum $`P^{}=(E_p^{},P_X^{},P_Y^{},P_Z^{})`$. The squared four-momentum transfer at the proton vertex is given by

$`t=(PP^{})^2P_T^2,`$
where $`P_T^2=P_X^2+P_Y^2`$, and the fraction of the beam momentum retained by the final proton is

$`x_L={\displaystyle \frac{|𝐏^{}|}{|𝐏|}}{\displaystyle \frac{E_p^{}}{E_p}},`$
where $`E_p=820`$ GeV denotes the incident proton beam energy. Both $`t`$ and $`x_L`$ are measured with the LPS. Other useful variables in diffractive DIS are:

$`x_{\mathrm{IP}}={\displaystyle \frac{(PP^{})q}{Pq}}={\displaystyle \frac{M_X^2+Q^2t}{W^2+Q^2M_p^2}}1x_L`$
and

$`\beta ={\displaystyle \frac{Q^2}{2(PP^{})q}}={\displaystyle \frac{x}{x_{\mathrm{IP}}}}{\displaystyle \frac{Q^2}{M_X^2+Q^2}},`$
where $`M_X`$ is the invariant mass of the final-state, $`X`$, determined as described below. In the Pomeron-exchange picture, $`x_{\mathrm{IP}}`$ is the fraction of the proton’s momentum carried by the Pomeron. For a Pomeron with partonic constituents, $`\beta `$ is then the fraction of the Pomeron’s momentum carried by the parton that absorbs the virtual photon.

### 3.2 Reconstruction of $`M_X`$

Two methods were used in this analysis to determine the invariant mass of the system $`X`$. The first method measured the mass directly from $`X`$. The individual objects belonging to $`X`$ were reconstructed by combining charged tracks measured in the CTD and neutral energy clusters measured in the calorimeter into energy flow objects (EFOs ). In this procedure, the tracks and clusters associated with the scattered positron were not used. The four-vectors of the EFOs were calculated assuming the pion mass. The invariant mass can then be reconstructed using these four-momenta as

$$M_{X,EFO}^2=(\underset{i}{}E_i)^2(\underset{i}{}P_{Xi})^2(\underset{i}{}P_{Yi})^2(\underset{i}{}P_{Zi})^2.$$
(1)
The second method inferred $`M_X`$ from the measurement of the final-state proton and the scattered positron, which enter into the reconstruction of the kinematic variables $`x`$ and $`y`$ via the double angle method. $`M_X^2`$ is then given by

$$M_{X,LPS}^2=sy(1x_Lx).$$
(2)
Since the LPS method has better resolution at higher values of $`M_X`$ and the EFO measurement is better at lower values, $`M_X`$ was evaluated as the weighted average of the values obtained with the two methods. Before combining the two $`M_X`$ values, correction factors obtained from Monte Carlo studies were applied. This resulted in a resolution on $`M_X`$ of about 25% at low $`M_X`$, improving to 15% at high $`M_X`$.

### 3.3 Boost into the $`\gamma ^{}\mathrm{IP}`$ center-of-mass system

The boost into the CMS of $`X`$ ($`\gamma ^{}\mathrm{IP}`$), Fig. 1(b), was determined from the four-momenta of the $`\gamma ^{}`$ and $`\mathrm{IP}`$. The $`\gamma ^{}`$ four-momentum was calculated from the energy and angles of the scattered positron. The best reconstruction of the $`\mathrm{IP}`$ four-momentum was obtained by combining the information from the LPS and the ZEUS main detector. The $`X`$ and $`Y`$ components of the $`\mathrm{IP}`$ four-momentum were taken to be equal to the negative of the $`X`$ and $`Y`$ momentum components of the scattered proton measured with the LPS. The $`Z`$ component was taken to be the total $`P_Z`$ of all EFOs ($`P_Z^{tot}=_iP_{Zi}`$) minus the $`Z`$ component of the $`\gamma ^{}`$ four-momentum. Finally, the energy component was calculated by requiring that the $`\gamma ^{}\mathrm{IP}`$ invariant mass be equal to the value of $`M_X`$ determined by the combination of the information from the LPS and the EFOs, as described earlier.

The resolution on the angle between the $`\gamma ^{}\mathrm{IP}`$ axis and the ZEUS $`Z`$-axis, as estimated by Monte Carlo studies, was approximately inversely proportional to $`M_X`$ and equal to 10 mrad at $`M_X=20`$ GeV. An improvement in the accuracy of the measurement by a factor of three was achieved using the LPS compared to a measurement without its use, i.e. substituting $`p_X=p_Y=0`$.

## 4 Event selection

Deep inelastic scattering events were selected in a manner similar to those of the ZEUS $`F_2`$ analysis  by requiring:

* the presence of a fully contained and isolated positron with energy greater than 10 GeV and matched to a charged track, if within the acceptance of the CTD;
* a value of $`y`$ greater than 0.05, evaluated using the Jacquet-Blondel method ;
* a value of $`\delta =E^{tot}P_Z^{tot}`$ in the range 40 $`<\delta <`$ 65 GeV, where $`E^{tot}`$ and $`𝐏^{\mathrm{𝐭𝐨𝐭}}`$ are, respectively, the total energy and momentum of the event in the main detector, including the scattered positron;
* the $`Z`$ coordinate of the event vertex within 50 cm of the nominal interaction point.

Diffractive events were then selected by requiring a well reconstructed LPS track  carrying more than $`95\%`$ of the beam energy ($`x_L>0.95`$). The LPS track was required to pass no closer than 0.04 cm to the HERA beampipe, and the quantity $`\delta ^+=E^{tot}+P_Z^{tot}+2E_px_L`$ was required to be less than $`2E_p+20`$ GeV to be consistent with longitudinal momentum conservation after taking the LPS resolution into account. Since the results presented here were found to be independent of $`t`$, no explicit cut on $`t`$ was applied.

To provide a sample for which the acceptance is large and uniform, to remove events measured with low resolution and to reject events corresponding to exclusive production of vector mesons, the following kinematic cuts were applied:

* $`70<W<250\mathrm{GeV}`$;
* $`0.0003<x_{\mathrm{IP}}<0.03`$;
* $`4<Q^2<150\mathrm{GeV}^2`$;
* $`4<M_X<35\mathrm{GeV}`$.

In addition, at least four EFOs in the system $`X`$ were required. These cuts define the kinematic range for all results presented in this paper, unless otherwise noted. The final data sample contained 2355 events.

## 5 Models of the diffractive final state

The data were compared with three Monte Carlo generators based on different theoretical models. For all generators, hadronization was simulated using the Lund string model as implemented in JETSET 7.4 . The models considered here produce either a $`q\overline{q}`$ pair or a $`q\overline{q}g`$ final state at the parton level. However, the dynamics of the production of these states is different for each model, thus yielding predictions that differ in their relative contributions of $`q\overline{q}`$ and $`q\overline{q}g`$ states, as well as in the final-state topology. Common to all three models is that the $`q\overline{q}g`$-type events dominate the final state at high masses. In such events, the gluon usually travels in the direction of the Pomeron.

### 5.1 Resolved Pomeron model

In the resolved Pomeron approach , the exchanged particle is assumed to have a partonic structure consisting of quarks and gluons. A sample of resolved Pomeron events was produced with the RAPGAP 2.08/06 generator . The hadronic final state is simulated in analogy to ordinary DIS. The Pomeron parton density functions (PDFs) used were those determined by the H1 collaboration from their measurement of $`d^3\sigma /dx_{\mathrm{IP}}dQ^2d\beta `$ .

When the virtual photon in a diffractive event interacts with a quark in the Pomeron, the lowest-order final state is a $`q\overline{q}`$. At $`𝒪(\alpha _s)`$, $`q\overline{q}g`$ final states can be produced by the radiation of an extra gluon via the QCD-Compton (QCDC) process. In addition, $`q\overline{q}g`$ final states are produced via the boson-gluon fusion (BGF) process between the virtual photon and a gluon in the Pomeron. The Pomeron remnant consists of a quark in the $`q\overline{q}`$ and QCDC events and of a gluon in BGF events. Additional parton radiation from the remnant partons is suppressed. The following RAPGAP options were used:

* the leading-order version of the H1 Fit 2 was used. In the fit, more than 80 % of the Pomeron momentum is carried by the gluons.
* the final-state parton system in $`q\overline{q}`$-type events was simulated with limited transverse momentum with respect to the $`\gamma ^{}\mathrm{IP}`$ axis (intrinsic transverse momentum, $`k_T`$), randomly distributed according to $`\mathrm{exp}(5.5k_T^2)`$ ($`k_T`$ in GeV);
* to avoid divergences in the matrix elements for the $`𝒪(\alpha _s)`$ processes for massless quarks, the transverse momentum squared of any outgoing parton was required to be larger than 3 GeV2;
* higher-order QCD radiation was simulated using initial- and final-state parton showers (MEPS) ;
* the generator was run for light flavors and charm production (produced via the BGF process) separately, and the final sample was obtained by mixing the two according to their relative cross sections as predicted by RAPGAP.

The curves corresponding to this model are labeled “RG resolved IP” in later figures. The model labeled “RG $`q\overline{q}`$ part only” in later figures was obtained by selecting events from the above sample that were not produced via the BGF process.

### 5.2 Photon dissociative model

In the rest frame of the proton, diffractive scattering can be viewed as the dissociation of the virtual photon into a $`q\overline{q}`$ pair (Fig. 1(c)) well before the interaction with the proton. For higher masses, the $`q\overline{q}g`$ final state (Fig. 1(d)) becomes important. The dissociated photon system couples to the proton by color-singlet exchange.

The simplest realization of this color singlet is the exchange of two gluons with opposite color charge (two-gluon exchange model). The RIDI2.0 program  implements a two-gluon model following the approach of Ryskin , where the diffractive dissociation is treated in the framework of the leading log approximation (LLA) of perturbative QCD. The cross section is proportional to the square of the gluon density of the proton, which was taken from the CTEQ4M  parameterization of the proton structure function with an appropriately chosen cut-off parameter for the transverse momentum of the final-state gluon. The contributions of both transversely and longitudinally polarized photons are included. A large theoretical uncertainty in the relative contributions of $`q\overline{q}`$ and $`q\overline{q}g`$ states remains. This uncertainty is reflected in the wide range allowed for the ($`\alpha _s`$-dependent) $`k`$-factors, which determine the relative cross sections of the $`q\overline{q}`$ and $`q\overline{q}g`$ contributions.

Another model based on photon dissociation is that of Golec-Biernat and Wüsthoff , in which the virtual photon splits into a $`q\overline{q}`$ or $`q\overline{q}g`$ color dipole. The interaction of this dipole with the proton can be described by an effective dipole cross section taking the dynamics of saturation into account. The dipole cross section was parameterized using the HERA measurements of the DIS total cross section between $`Q^2=0.1`$ GeV2 and $`Q^2=100`$ GeV2. The production of $`q\overline{q}`$ and $`q\overline{q}g`$ final states is calculated using these parameters. The relative contribution of $`q\overline{q}`$ and $`q\overline{q}g`$ events is fixed by QCD color factors. This model is implemented in the SATRAP Monte Carlo generator  which uses the additional approximation of strong ordering in the transverse momenta: $`k_T(g)k_T(q,\overline{q})`$. In the course of the ZEUS analysis of diffractive 3-jet events , it was found that the modeling of higher-order QCD processes was inadequately implemented in SATRAP. Most notably, no initial-state parton cascades were included, and the final-state QCD radiation from the gluon in the dominant $`q\overline{q}g`$ contribution was suppressed. A new implementation of higher-order QCD processes in SATRAP was carried out , in which the color-dipole model (CDM) was implemented in a similar fashion to that in RAPGAP. This model is referred to as SATRAP-CDM.

### 5.3 JETSET

If the diffractive DIS hadronic system, $`X`$, is produced by a virtual photon splitting into a $`q\overline{q}`$ pair, its properties would be expected to be similar to those of the hadronic system produced by $`e^+e^{}`$ annihilation at a CMS energy $`\sqrt{s}=M_X`$. The $`e^+e^{}`$ final state, as simulated by the JETSET program , was used as a baseline to which both the data and the diffractive Monte Carlo generators were compared. JETSET is known to describe accurately many details of the final state in $`e^+e^{}`$ collisions.

## 6 Data correction and background

Monte Carlo simulations were used to correct the data for the resolution and acceptance of the main detector and the LPS, and to estimate the size and influence of the background. To estimate the model dependence of these corrections, two different generators were used: RAPGAP  and RIDI . RAPGAP events were generated with the H1 QCD Fit 2  for the IP structure function, as described in Section 5.1, except that CDM rather then MEPS was used for higher-order initial- and final-state QCD radiation. This model gives a poor description of the hadronic final state as well as of some kinematic variables as measured at the detector level. To obtain a sample which gives a good description of all observed distributions, the RAPGAP generator was reweighted to reproduce the data. The corresponding distributions are labeled “mod. RG” in Figs. 2 to 5. A similar reweighting procedure was carried out for RIDI. Both samples were passed through a GEANT 3.13  simulation of the ZEUS detector, subjected to the same trigger requirements as the data and processed by the same reconstruction programs.

In addition, changes in the HERA beamline parameters during the running period and the finite resolution of the motors that determine the position of the LPS, neither of which was simulated in the Monte Carlo, were taken into account by reweighting and smearing the simulated $`x_L`$ distribution.

The primary source of background in this analysis is the accidental overlap of a DIS event with an unrelated beam-halo proton measured in the LPS. This contribution was assumed to be constant over the running period. To estimate its size, a sample of such background events was studied. Background protons were identified by selecting DIS events with $`E^{tot}+P_Z^{tot}>100`$ GeV (using only CAL) and an LPS track with $`x_L>0.9`$. From energy conservation allowing for detector resolution, such events must result from an accidental overlap. The LPS information from these events was then combined with a sample of non-diffractive DIS Monte Carlo events generated using DJANGOH . The resulting sample was normalized to match the upper tail of the $`\delta ^+`$ distribution for the data shown in Fig. 2. With this method, the level of background contamination, after all selection cuts, was estimated to be 5.2%. Its contribution was statistically subtracted in all results presented below.

Figure 3 shows the distributions of $`Q^2`$, $`W`$, $`x_L`$ and $`t`$ for events which pass all selection cuts except the ones imposed on the variables shown (indicated by the arrows). Also shown in the figure are the distributions from the modified RAPGAP simulation plus background which was normalized so that the weighted sum of the modified RAPGAP and the background Monte Carlo events describes the $`Q^2`$ distribution in the selected range, see Fig. 3 (a). The kinematic distributions of the data, with the exception of the $`t`$ distribution, are well described by the modified RAPGAP Monte Carlo sample in combination with the background sample.

Figure 4 (a-c) shows some of the properties of the hadronic final state: the distribution of the number of EFOs, the EFO energy in the CMS frame of the $`\gamma ^{}\mathrm{IP}`$ system and $`\eta _{max}`$, where $`\eta _{max}`$ is the pseudorapidity of the most-forward calorimeter cluster of energy greater than 400 MeV. The data are well described by the modified RAPGAP sample in combination with the background sample. Figure 4 (d) shows the average $`M_X`$ reconstructed using the LPS (Eq. (2)) as a function of the invariant mass calculated using EFOs (Eq. (1)). The correlation seen in the data is well reproduced by the Monte Carlo simulation. The fact that the correlation is approximately linear in the selected mass range confirms that the final state is well contained in the ZEUS detector.

Figure 5 shows the distributions of $`M_X`$, $`x_{\mathrm{IP}}`$ and $`\beta `$. After the reweighting procedure, all data distributions are well described by the modified RAPGAP Monte Carlo sample. This reweighting, which mainly affected the low $`M_X`$ and high $`x_{\mathrm{IP}}`$ regions, had little effect on the $`\beta `$ distribution which, although peaking at low values, is well described by the relatively flat Pomeron PDFs used in RAPGAP. RIDI (not shown in Figs. 2 to 5) also reproduces the data satisfactorily.

## 7 Systematic uncertainties

The systematic uncertainties were obtained by studying the sensitivity of the results to the uncertainties in the understanding of the ZEUS detector response and by modifying the analysis procedures as listed below.

Uncertainties related to the understanding of the ZEUS main detector include:

* the uncertainty on the absolute energy scales in each major section of the calorimeter (FCAL, BCAL, and RCAL), which for this analysis was understood to a level of 3% for hadrons and 2% for the scattered positron;
* the uncertainty on the survey measurements of the positions of the major sections of the calorimeter with respect to each other and to the HERA beamline, which are accurate to 1–2 mm.

Uncertainties related to the LPS include:

* possible shifts in $`x_L`$, $`p_X`$, and $`p_Y`$ of $`\pm 0.003`$, $`\pm 10`$ MeV, and $`\pm 50`$ MeV, respectively, as determined from an analysis of elastic $`\rho `$ photoproduction ($`epe\rho p`$);
* the uncertainty on the beampipe position, as determined by alignment studies, was taken into account by changing the cut on the distance of closest approach of the proton track to the beampipe by $`\pm 400`$ $`\mu \mathrm{m}`$;
* the uncertainty in the Monte Carlo simulation of the proton track reconstruction was taken into account by applying tighter cuts on the quality of the Monte Carlo track.

To check for the effect of possible particle losses into the forward beampipe, and to account for the differences in modeling this forward region in the Monte Carlo generators, the energy deposited in the inner ring of the FCAL, which covers approximately the pseudorapidity range of $`2.7<\eta <3.9`$, was scaled by $`\pm 25\%`$, as suggested by MC studies.

The normalization of the background from an accidental overlap of a DIS event with an unrelated proton was estimated using the $`\delta ^+`$ distribution of Fig. 2. The uncertainty caused by the background subtraction was determined by changing the background normalization by $`\pm 50\%`$.

To estimate the model dependence of the results, the data were also corrected using RIDI. The difference between the RIDI and modified RAPGAP results was taken as an estimate of the model uncertainty. This uncertainty was assumed to be symmetric with respect to the nominal results obtained using the modified RAPGAP sample.

The largest contributions to the systematic uncertainty typically originate from the model dependence and the uncertainty in the $`x_L`$ reconstruction. All systematic uncertainties were assumed to be independent and were calculated separately for positive and negative variations with respect to the nominal value. The total positive and negative systematic uncertainties were calculated as the corresponding sums in quadrature.

## 8 Global event-shape variables

The results presented here are corrected to the hadron level in the kinematic range defined in Section 4 with at least four stable particles in the final state. A particle is considered stable if its lifetime is longer than $`3\times 10^{10}`$ s. If the lifetime is shorter than $`3\times 10^{10}`$ s, the daughter particles (with lifetime $`>3\times 10^{10}`$ s) are included. This corresponds to the final-state definition used by $`e^+e^{}`$ experiments .

### 8.1 Thrust and sphericity

The event-shape-variable thrust ($`T`$) is calculated by determining the unit vector $`\widehat{𝐧}`$ which maximizes

$`T(\widehat{𝐧})={\displaystyle \frac{_i|\widehat{𝐧}𝐩_i|}{_i|𝐩_i|}},`$
where $`𝐩_i`$ is the three-momentum of a final-state particle and the sum is over all particles belonging to the system under study. The resulting axis $`\widehat{𝐧}`$ is called the thrust axis and $`T(\widehat{𝐧})`$ is the thrust value. For collimated two-jet events, the value of $`T`$ approaches 1, while events with an isotropic shape yield values close to 0.5.

The sphericity ($`S`$) is defined as

$`S={\displaystyle \frac{3}{2}}(\lambda _2+\lambda _3),`$
where $`\lambda _2`$ and $`\lambda _3`$ are the two smallest eigenvalues (corresponding to orthogonal eigenvectors) of the sphericity tensor

$`S^{\alpha \beta }={\displaystyle \frac{_ip_i^\alpha p_i^\beta }{_i|𝐩_i|^2}}\alpha ,\beta =x,y,z.`$ (3)
For isotropic events, $`S`$ approaches 1, and for collimated two-jet events, $`S`$ is close to 0.

The sphericity and thrust distributions in the CMS of the final state in diffractive DIS are shown in Fig. 6 for two bins of $`M_X`$, and are compared to the ones observed in $`e^+e^{}`$ annihilation from the TASSO collaboration  at $`\sqrt{s}=`$$`M_X`$. The diffractive events show thrust (sphericity) distributions which are broader and shifted to lower (higher) values compared to the $`e^+e^{}`$ results, indicating that they are more isotropic. The diffractive events become more collimated as $`M_X`$ increases, a trend also observed in $`e^+e^{}`$ annihilation. If the virtual photon in diffractive DIS fluctuates only into a $`q\overline{q}`$ state, the resulting hadronic final state should develop in a manner similar to $`e^+e^{}q\overline{q}`$ at a CMS energy $`\sqrt{s}=`$$`M_X`$. Deviations from this $`e^+e^{}`$-type of behavior are expected, however, since there is a significant contribution from the $`q\overline{q}g`$ diagram in the color field of the proton (Fig. 1(d)), which is not present in $`e^+e^{}`$ collisions, where gluons can only be produced via higher-order QCD radiation.

Figures 7 (a) and (c) show the average values of thrust, $`T`$, and sphericity, $`S`$, measured in six bins of $`M_X`$. The value of $`T`$ increases and $`S`$ decreases with increasing $`M_X`$ with slopes similar to those found in the $`e^+e^{}`$ data . On average, the diffractive events are less collimated than $`e^+e^{}`$ events at a similar CMS energy.

Figures 7 (b) and (d) show the same data, but now compared to several Monte Carlo models222Henceforth, all comparisons of the data are made to MC models that have not been reweighted in the manner described in Sec. 6.. Also plotted in Fig. 7 (b) are the average thrust values measured by the H1 collaboration using a diffractive DIS sample tagged with rapidity gaps . The present measurements are in good agreement with the H1 results which have been corrected to the full phase space by a MC simulation. The results are also consistent with those presented in a previous ZEUS publication , which were obtained for a small range in $`\eta `$.

The Monte Carlo models of diffractive scattering reproduce the general trend of the data. RAPGAP and SATRAP-CDM give a good description of $`T`$ over almost the full $`M_X`$ range, while RIDI fails for $`M_X\mathrm{}20`$ GeV. The diffractive models produce events more isotropic than the ones generated with JETSET and measured in $`e^+e^{}`$ annihilation. This is mainly due to the inclusion of $`q\overline{q}g`$ final states. Figures 7 (b) and (d) also show the RAPGAP contribution for the $`q\overline{q}`$ diagram alone. These curves demonstrate the need for the $`q\overline{q}g`$ contribution, especially at the higher $`M_X`$ values.

The final state measured in terms of thrust and sphericity has no strong dependence on any kinematic variable other than $`M_X`$. Figure 8 shows $`T`$ as a function of $`M_X`$ for $`x_{\mathrm{IP}}<0.01`$, where Pomeron exchange dominates, and for $`x_{\mathrm{IP}}>0.01`$, where Reggeon exchange may become more important. No significant differences are observed and all models give a reasonable description of the data. The data were also split into two samples of $`Q^2`$, $`t`$ and $`x`$; the dependence on these variables was less pronounced than that for $`x_{\mathrm{IP}}`$.

### 8.2 Transverse momenta in and out of the event plane

The shape of the system $`X`$ was also studied in terms of two more variables, $`P_{T,in}^2`$ and $`P_{T,out}^2`$, which measure the transverse momentum in and out of the event plane. The event plane is defined by the eigenvectors of the sphericity tensor associated with the two largest eigenvalues, $`\lambda _1`$ and $`\lambda _2`$ (see Eq. (3)). These transverse momenta have also been studied in $`e^+e^{}`$ experiments, since they are sensitive to gluon bremsstrahlung. They are evaluated by multiplying the average momentum squared of the $`n`$ particles in the event by the two smallest sphericity eigenvalues $`\lambda _2`$ and $`\lambda _3`$, respectively:

$`P_{T,in}^2`$ $`=`$ $`{\displaystyle \frac{_ip_i^2}{n}}\lambda _2,`$ (4)
$`P_{T,out}^2`$ $`=`$ $`{\displaystyle \frac{_ip_i^2}{n}}\lambda _3,`$ (5)
where $`p_i^2`$ is the squared momentum of the $`i`$th particle in the event. By definition, $`P_{T,in}^2>P_{T,out}^2`$. The observation of a difference in the $`M_X`$ behavior of these two variables could be explained by the presence of events of planar shape, such as events with three partons in the final state.

Figures 9 (a) and (b) show the average values of $`P_{T,in}^2`$ and $`P_{T,out}^2`$ for the diffractive DIS sample as a function of $`M_X`$. On average, $`P_{T,in}^2`$ is about a factor of three larger than $`P_{T,out}^2`$. The dependence of $`P_{T,in}^2`$ on $`M_X`$ is almost linear over the range studied, while the $`M_X`$ dependence of $`P_{T,out}^2`$ becomes weaker as $`M_X`$ increases, indicating that the diffractive events become more planar. Figure 9 also shows the predictions of the diffractive DIS generators and JETSET. The RAPGAP and SATRAP-CDM samples give a good description of $`P_{T,in}^2`$ for low masses, but underestimate its value for $`M_X\stackrel{>}{}15`$ GeV. RIDI overestimates $`P_{T,in}^2`$ for most of the measured $`M_X`$ region. Transverse momentum out of the event plane can be produced by the fragmentation process. While RAPGAP and SATRAP-CDM give a reasonable description of the data, RIDI slightly overestimates the data at low $`M_X`$ values. The transverse momentum distributions in and out of the event plane predicted by JETSET are in excellent agreement with the diffractive scattering data.

In summary, the diffractive hadronic final state shows trends in $`T`$ and $`S`$ different from those observed in $`e^+e^{}`$ annihilation. These differences can be explained by the presence of $`q\overline{q}g`$ events from BGF (RAPGAP) or from photon dissociation (SATRAP-CDM, RIDI) which are absent in $`e^+e^{}`$ annihilation and JETSET, where the only source of gluons in the final state is from QCD radiation. The transverse momentum out of the event plane obtained from diffractive Monte Carlo models and JETSET reproduces the data, indicating that the hadronization process is consistently modeled and is similar in diffractive DIS and $`e^+e^{}`$ annihilation.

## 9 Thrust axis orientation

The orientation of the entire hadronic final state, $`X`$, with respect to the $`\gamma ^{}\mathrm{IP}`$ axis can be studied by measuring the angle $`\theta _{\mathrm{Thrust}}`$ (as shown in Fig. 1 (b)) between this axis and the thrust axis.

The values of $`\theta _{\mathrm{Thrust}}`$ and the transverse momentum, $`P_\mathrm{T}^{\mathrm{Thrust}}`$, relative to the thrust axis, defined by

$`P_\mathrm{T}^{\mathrm{Thrust}}={\displaystyle \frac{M_X}{2}}\mathrm{sin}\theta _{\mathrm{Thrust}},`$ (6)
are sensitive to various sources of transverse momentum. These include the intrinsic transverse momentum, $`k_T`$, of the partons in the proton, which was found to be 1.69$`\pm 0.18_{0.20}^{+0.18}`$ GeV in a recent ZEUS publication , intrinsic transverse momentum of the partons in the Pomeron, and contributions of hard processes, e.g. final-state gluon emission .

Figure 10 (a) shows the average thrust angle $`\theta _{\mathrm{Thrust}}`$ as a function of $`M_X`$. The angle is largest at small masses and falls off steeply as $`M_X`$ increases, indicating that the event becomes more aligned with the $`\gamma ^{}\mathrm{IP}`$ axis as the CMS energy increases. This general trend is described by all Monte Carlo models shown, but only RIDI produces a thrust angle as large as that observed in the data. SATRAP-CDM and RAPGAP fall about equally short of the observed angles. Given that the hadronization is well described by all Monte Carlo models (see Section 8.2), the measurement of $`\theta _{\mathrm{Thrust}}`$ implies that additional sources of transverse momentum as discussed above must be important. Indeed, the approach implemented in RIDI favors the production of partons with relatively large transverse momenta, typically of the order of 1 GeV .

Figure 10 (b) shows that the average transverse momentum, $`p_\mathrm{T}^{\mathrm{Thrust}}`$, produced in diffractive scattering is almost independent of $`M_X`$ for $`M_X>10`$ GeV, with an average value of about 2 GeV. In this region of $`M_X`$, the independence is correctly reproduced by all Monte Carlo models, but the RAPGAP prediction is too low. The need for the $`q\overline{q}g`$ contribution is illustrated by the curve showing only that part of the RAPGAP prediction corresponding to the $`q\overline{q}`$ final state.

## 10 Energy flow

Another measurement of the event topology is the distribution of energy as a function of the pseudorapidity of a hadron with respect to the $`\gamma ^{}\mathrm{IP}`$ CMS axis. This distribution, commonly referred to as the energy flow, is shown in Fig. 11 for data in three $`M_X`$ ranges. The data for low $`M_X`$ have a Gaussian-like shape, but for $`M_X>7.5`$ GeV a structure with two peaks develops. This structure becomes more pronounced as the mass increases. A slight asymmetry develops in the data, with more energy being produced in the $`\gamma ^{}`$ hemisphere. For comparison, the predictions from the RAPGAP, SATRAP-CDM and RIDI Monte Carlo generators are also shown. The SATRAP-CDM and RAPGAP generators predict too broad a rapidity distribution, and display a separation between the $`\gamma ^{}`$ and IP hemispheres in the lowest $`M_X`$ bin that is not exhibited by the data. RIDI gives a reasonable description of the data in the lowest $`M_X`$ bin but has a different shape at higher $`M_X`$. The asymmetry indicated in the data is slightly larger than that produced by the diffractive Monte Carlo models. For those values of $`|\eta |`$ at which the H1 collaboration has also published data , the energy flow is in good agreement, except at the highest $`M_X`$ where the H1 data are somewhat narrower.

Figure 12 shows the distribution of transverse energy, $`E_T=E_i\mathrm{sin}\theta _i`$, as a function of $`\eta `$. Similar discrepancies between the data and the Monte Carlo events were observed for the transverse-energy flow as was observed for the energy-flow distribution.

## 11 Seagull distribution

The small asymmetry between positive and negative $`\eta `$ observed in the energy-flow plot of Fig. 11 can be further investigated using the distribution of transverse momentum of the particles belonging to the system $`X`$.

In inclusive DIS, $`epeX`$, the fragmentation of the hadronic system occurs between the struck quark (forming the system X) and the remaining quarks in the proton (forming the proton remnant). Measurements in the Breit frame  have shown that, whereas the particle multiplicity and momentum distributions in the hemisphere of the struck quark are roughly consistent with those measured in $`e^+e^{}q\overline{q}`$, particles are produced with smaller average transverse momentum in the proton-remnant hemisphere.

In diffractive DIS, an asymmetry in the momentum distribution between the $`\gamma ^{}`$ and IP hemispheres could be observed if the Pomeron behaves as an extended object like the proton and produces a remnant after the scattering process. This asymmetry is usually investigated by studying single-particle distributions as a function of the scaled longitudinal momentum, $`x_F`$, and the momentum transverse to the photon direction, $`p_T`$. The photon direction defines the longitudinal axis in the $`\gamma ^{}\mathrm{IP}`$ CMS as well as in the $`\gamma ^{}p`$ CMS. If $`p_T`$ and $`p_L`$ are the momentum components of a final-state hadron perpendicular and parallel, respectively, to this axis, the variable $`x_F`$ is given by:

$`x_F=p_L/p_L^{max},`$
where positive $`x_F`$ is in the direction of the $`\gamma ^{}`$ and $`p_L^{max}`$ is the maximum kinematically allowed value of $`p_L`$. In the $`\gamma ^{}\mathrm{IP}`$ CMS, $`p_L^{max}=M_X/2`$. For the $`\gamma ^{}p`$ case, $`p_L^{max}=W/2`$.

The average $`p_T^2`$ of particles as a function of $`x_F`$, commonly referred to as the “seagull plot”, is shown in Fig. 13 for $`11<M_X<17.8`$ GeV. Also plotted are $`\gamma ^{}p`$ data from the EMC  collaboration at $`W=14`$ GeV, equal to the average value of $`M_X`$ in this bin. The EMC DIS $`\mu p\mu X`$ data indicate a suppression of the average $`p_T^2`$ associated with a proton remnant which is not as apparent in the diffractive data in this mass range.

Shown in Fig. 14 (a)-(c) is the seagull plot for three different $`M_X`$ bins compared with the predictions from RAPGAP, SATRAP-CDM and RIDI. The data exhibit a growing asymmetry as $`M_X`$ increases. This asymmetry can also be seen in the ratios of the average squared transverse-momentum in the $`\gamma ^{}`$ and IP hemispheres as a function of $`|x_F|`$ (Fig. 14 (d)-(f)). The data are well reproduced by both RAPGAP and RIDI, while for $`M_X>16`$ GeV, the transverse momentum generated by SATRAP-CDM is smaller than in the data. RAPGAP and SATRAP-CDM describe the size of the asymmetry correctly in all $`M_X`$ bins, while for $`M_X>16`$ GeV RIDI slightly underestimates the transverse momentum in the IP direction, resulting in a slightly larger asymmetry than that seen in the data. RAPGAP produces this asymmetry by including a Pomeron remnant. RIDI and SATRAP-CDM, in contrast, produce the asymmetry via the $`q\overline{q}g`$ diagram.

## 12 Conclusion

A study of the hadronic system, $`X`$, in the reaction $`epeXp`$ has been reported for the kinematic range $`4<M_X<35`$ GeV, $`4<Q^2<150`$ GeV2, $`0.0003<x_{\mathrm{IP}}<0.03`$ and $`70<W<250`$ GeV. The use of the LPS allows diffractive events to be tagged without applying cuts on the system $`X`$. It also provides a powerful constraint on the diffractive kinematics, allowing, for example, an accurate determination of the $`\gamma ^{}\mathrm{IP}`$ axis in the center-of-mass frame of the system $`X`$.

The diffractive hadronic final state becomes more collimated as the invariant mass $`M_X`$ of the system increases. This trend is similar to the one observed in $`e^+e^{}`$ annihilation. However, on average the diffractive final state is more isotropic. This can be attributed to contributions not present in $`e^+e^{}`$ annihilation, such as the boson-gluon fusion process in the resolved Pomeron approach, or $`q\overline{q}g`$ production from the dissociation of the virtual photon.

The mean transverse momentum out of the event plane is similar to that found in $`e^+e^{}`$ annihilation, indicating the universality of the hadronization. Even after considering the broadening effects of hadronization, it is apparent that more $`k_T`$ than is usually associated with the resolved Pomeron Monte Carlo approach is required to accommodate the large thrust angle and narrow energy flows at low $`M_X`$.

Particle production becomes asymmetric along the $`\gamma ^{}\mathrm{IP}`$ axis as $`M_X`$ increases, resulting in more average transverse momentum in the virtual-photon hemisphere. This asymmetry is consistent with both the concept of a remnant in the resolved Pomeron model and with the production of $`q\overline{q}g`$ final states in the photon-dissociation approach. The comparison with the Monte Carlo models suggests that a dominant gluon contribution to the partonic final state is necessary.

The invariant mass, $`M_X`$, of the hadronic system is the only variable upon which the characteristics of the system was found to depend. The system is independent of the momentum transfer at either the positron vertex, $`Q^2`$, or at the proton vertex, $`t`$. Neither is there any dependence on $`x`$ or on the fractional momentum of the Pomeron, $`x_{\mathrm{IP}}`$.

Many models of diffraction are able to reproduce the measured diffractive cross sections. However, none of the models discussed here is able to describe all aspects of the data. It is clear, therefore, that measurements of the detailed characteristics of diffractive events, such as presented in this paper, will become more and more crucial in understanding the underlying physics of diffraction in deep inelastic scattering.

## Acknowledgments

This measurement was made possible by the inventiveness and the diligent efforts of the HERA machine group. The strong support and encouragement of the DESY directorate has been invaluable. The design, construction, and installation of the ZEUS detector has been made possible by the ingenuity and dedicated effort of many people from inside DESY and from the home institutes who are not listed as authors. Their contributions are acknowledged with great appreciation. We acknowledge support by the following: the Natural Sciences and Engineering Research Council of Canada (NSERC); the German Federal Ministry for Education and Science, Research and Technology (BMBF), under contract numbers 057BN19P, 057FR19P, 057HH19P, 057HH29P, 057SI75I; the MINERVA Gesellschaft für Forschung GmbH, the Israel Science Foundation, the U.S.-Israel Binational Science Foundation, the Israel Ministry of Science and the Benozyio Center for High Energy Physics; the German-Israeli Foundation, the Israel Science Foundation, and the Israel Ministry of Science; the Italian National Institute for Nuclear Physics (INFN); the Japanese Ministry of Education, Science and Culture (the Monbusho) and its grants for Scientific Research; the Korean Ministry of Education and Korea Science and Engineering Foundation; the Netherlands Foundation for Research on Matter (FOM); the Polish State Committee for Scientific Research, grant no. 2P03B04616, 620/E-77/SPUB-M/DESY/P-03/DZ 247/2000 and 112/E-356/SPUB-M/DESY/P-03/DZ 3001/2000, and by the German Federal Ministry for Education and Science, Research and Technology (BMBF); the Fund for Fundamental Research of Russian Ministry for Science and Education and by the German Federal Ministry for Education and Science, Research and Technology (BMBF); the Spanish Ministry of Education and Science through funds provided by CICYT; the Particle Physics and Astronomy Research Council, UK; the US Department of Energy; the US National Science Foundation.


text:
# Fast Tree Search for Enumeration of a Lattice Model of Protein Folding

## I INTRODUCTION

The protein folding problem  has long attracted the attention of scientists from various disciplines. The relationship between the amino-acid sequence and the three-dimensional structure of a protein is not only an extremely important and practical problem in biology, but also a fundamental problem in science. Despite a tremendous amount of effort and progress over many decades, the problem remains essentially unsolved. At least part of the difficulty arises from the intrinsic complexity of the protein-folding problem. Since the seminal work of Anfinsen  about 40 years ago, it has been demonstrated that the native state of a small, single domain protein is the global minimum of the free energy. However, the minimum-free-energy conformation of a polypeptide chain is “hiding” in a large space of $`z^N`$ conformations, where $`N`$ is the length of the chain and $`z`$ is the effective coordination number. Even if we count only the compact conformations for which $`z2`$ (see below) and take $`N=100`$ for typical small proteins, the number of conformations is huge, $`z^N2^{100}10^{30}`$. On top of this huge conformational space lies the heterogeneity of amino acids. The 20 natural amino acids differ in size, hydrophobicity, and other physical and chemical properties. This heterogeneity is coupled with two intrinsic features of polymers: chain connectivity and the excluded volume effect. The free energy is a sensitive and complicated function in this huge conformational space with complex constraints.

In the last decade or so, there has been increasing interest in studying simple lattice models of protein folding. In these models, polypeptide chains are represented by self-avoiding walks on a regular lattice (e.g., Fig. 1), greatly simplifying the conformational space. Very often the sequence space and hence the heterogeneity is also simplified by using only two types of amino acids: hydrophobic (H) and polar (P). These so-called “HP lattice models”  nonetheless capture some essential features of the protein folding problem. Simple lattice models have been applied to a wide range of problems including collapse and folding transitions , the influence of packing on secondary-structure formation , and differences in the designability of structures . The advantage of HP lattice models is that they are simple enough to be amenable to thorough theoretical study. These studies can provide fruitful insights to feed back to or test against realistic models and experiments.

One approach for calculating thermodynamic and other properties of lattice models is to enumerate all possible sequences and conformations . Since native globular protein structures and presumably most of the low-energy states are compact, enumeration studies are usually done on compact conformations only. The number of compact conformations $`C_N`$ scales with the chain length $`N`$ as $`C_Nz^N`$, where $`z`$ is an effective coordination number. Numerical estimation gives $`z1.47`$ for 2D square lattices  and $`z1.86`$ for 3D cubic lattices , in good agreement with mean-field calculations of $`Z/e`$  and $`(Z1)/e`$ , respectively, where $`Z`$ is the coordination number of the lattice and $`e=2.718\mathrm{}`$ is the base of the natural logarithm. (For real peptide chains the number of “distinct” states an amino acid can take, as estimated very roughly from Ramachandran plots of dihedral-angle frequencies, is about 5 or 6, which gives $`z2`$.) For HP models in which there are only two types of amino acids, the number of sequences is $`2^N`$. If in the enumeration study the energies of every sequence folded into every compact conformation are evaluated, the total number of energy calculations is $`2^N\times C_N`$. The largest system previously evaluated in this way is an $`N=27`$ ($`3\times 3\times 3`$) cubic lattice model  where $`2^N\times C_N=2^{27}\times 10334610^{13}`$. Several interesting results were found in the enumeration of the 27-mer, in particular the idea of designability and its relation to thermodynamic stability . However, $`N=27`$ is still small compared with typical protein sizes. It would be very desirable to enumerate larger systems if at all possible. In this paper, we report results of a complete enumeration of an $`N=36`$ ($`4\times 3\times 3`$) cubic-lattice model. The number of compact conformations is $`84,731,192`$ , so, naively, the total number of energy calculations is $`2^{36}\times 84,731,1926\times 10^{18}`$. The task was made possible by using a binary model, a fast tree-search algorithm yielding a speed-up factor of 1600, and a 53-processor 200MHz Pentium Pro cluster.

## II THE MODEL

The protein folding model we use in the enumeration study is the solvation model discussed in Ref. : Denote a sequence of amino acids by $`\{\sigma _i\}`$. We take only two types of amino acids: hydrophobic H ($`\sigma =1`$) and polar P ($`\sigma =0`$). A “structure” is the set of all reflections and rotations of a given compact conformation. The energy of a sequence folded into a structure is taken to be the sum of the contributions from each amino acid upon burial away from water:

$$E=\underset{i=1}{\overset{N}{}}\sigma _is_i,$$
(1)
where $`s_i`$ is a structure-dependent number characterizing the degree of burial of the $`i`$-th amino acid in the chain. Larger $`s_i`$ corresponds to smaller surface area accessible to solvent. For a $`4\times 3\times 3`$ structure on a cubic lattice there are 4 different kinds of sites: center, face, edge, and corner (see Fig. 1). So in principle there could be 4 different values of $`s_i`$. To simplify the calculation, we take only two values for $`s_i`$: we define a string $`\{s_i\}`$ for each structure with $`s_i=1`$ if the $`i`$-th site is a “core” (center or face) and $`s_i=0`$ if it is a “surface” (edge or corner). Thus each compact structure of $`4\times 3\times 3`$ is mapped into a string of 1’s and 0’s, $`\{s_i\}`$, with 12 1’s (“cores”) and 24 0’s (“surfaces”). The surface-to-core ratio is 2, close to the values for small natural proteins. With this simplification, the energy in Eq. (1) is just minus the dot product of two binary strings. For a given sequence $`\{\sigma _i\}`$, a ground-state structure is one that minimizes Eq. (1). A sequence may have more than one ground state structure, but we will be primarily interested in sequences with unique ground states. Out of the $`84,731,192`$ compact structures, the number of distinct structure strings is $`14,062,236`$, among which $`2,662,050`$ (corresponding to $`1,331,025`$ lattice conformations) each represent exactly one structure. Each of the remaining $`11,400,186`$ 36-bit strings represents multiple structures. We also analyzed a $`3\times 3\times 3`$ model with 7 “cores” (1 center and 6 faces) and 20 “surfaces” (12 edges and 8 corners). In this case, there are $`103,346`$ structures and $`6,291`$ distinct structure strings, among which only $`120`$ (corresponding to 60 lattice conformations) represent exactly one structure apiece.

## III Tree Searching Algorithm

Because the protein chains in our model are considered to be directed, there are generally two structures for each geometrical conformation, related by reversal of the direction of the chain (see Fig. 1). A small subset of structures are their own reversals. A structure string can be its own reversal even if its associated structure is not reversal symmetric. We found that among the $`14,062,236`$ structure strings in the $`4\times 3\times 3`$ model, there are $`2,850`$ which are their own reversals. The remaining $`14,059,386`$ strings form $`7,029,693`$ pairs, with the two members of each pair being the reversals of each other. To reduce memory use, we keep only one member of each such pair, with an extra bit tagged on the string to indicate that it actually represents two strings: itself and its reversal. There are thus $`7,029,693+2,850=7,032,543`$ distinct strings which we keep in the calculation.

Each of the $`7,032,543`$ distinct 36-bit structure strings $`\{s_i\}`$ has exactly 12 1’s and 24 0’s. Our goal is to find a way, given a 36-bit sequence string $`\{\sigma _i\}`$, to find if there is a unique entry in the table which maximizes the dot product of the structure string with $`\{\sigma _i\}`$, or, equivalently, which minimizes the energy of the sequence $`\{\sigma _i\}`$ according to Eq. (1).

To do this rapidly, we organized the strings in the table into a binary tree (Fig. 2). First we describe how the tree is organized, and then later how the tree was actually constructed.

Each node in the tree represents a subset of the 7 million strings in the table. For each node, the following information is maintained:

* Known-ones: A 36-bit string which has a 1 at the $`i`$’th position if and only if all the table entries corresponding to this node have 1’s at the $`i`$’th position.
* Undecided: A 36-bit string which has a 1 at the $`i`$’th position if and only if there is a table entry in this node which has a 1 at the i’th position and there is another table entry in this node which has a 0 at the $`i`$’th position.
* Missing-ones: Each string in a node will have 1’s at some undecided positions. For each string, Missing-ones is the sum of these 1’s. By construction each string has exactly 12 1’s, so Missing-ones is equal to 12 minus the sum of Known-ones. That is, Missing-ones is a single integer no greater than 12 for each node.

If the node is not a leaf, then it also contains a position number $`i`$ and two child nodes. These children partition the entries in the parent according to the value of the indicated position $`i`$: one child has all the parent entries where $`i=1`$ and the other child has all of the entries where $`i=0`$. Each leaf node at the end of the tree contains a small list of structure strings – in practice we found 16 strings per leaf to work best.

Given a 36-bit sequence string $`\{\sigma _i\}`$ and a node of the tree, what bounds can we place on the dot product of $`\{\sigma _i\}`$ and all structure strings represented by the node? Clearly, for all strings in the node the total dot product is at least as big as the dot product of $`\{\sigma _i\}`$ with Known-ones. On the other hand, the total dot product is at most $`\{\sigma _i\}`$ dotted with Known-ones plus the dot product of $`\{\sigma _i\}`$ with Undecided. Another upper bound for the total dot product is $`\{\sigma _i\}`$ dotted with Known-ones plus the integer Missing-ones.

Given such a tree, where the root corresponds to all of the 7,032,543 entries, and a 36-bit sequence string $`\{\sigma _i\}`$, here is how we search the table:

Compute an upper bound for the total dot product using the smaller of the two upper bounds described above. Call this the “goal” $`G`$. If there are any entries which achieve this goal, we are done. If not, repeat with the reversed version of the sequence. (This is necessary because our table contains only one member of each reversed pairs of structure strings. Whenever a ground state structure string is found for a sequence, the reversed sequence necessarily has the reversed structure string as a ground state.) Again, it we achieve the goal $`G`$, we are done. If not, decrease the goal $`G`$ by 1 and try again. Repeat until the goal is satisfied.

Given a goal $`G`$, we search the tree as follows, starting at the root node:

If the bound on the node indicates that goal $`G`$ is unachievable, return failure.

If the node is a leaf node, check each entry. If one is found that satisfies the goal, return success. If not, then return failure.

If the node is not a leaf node, try each of the children. We first try the child that matches $`\{\sigma _i\}`$. That is, if the children split on the value at the $`j`$’th position, then we refer to bit $`j`$ of $`\{\sigma _i\}`$. If it is a 1, then we first do the child having all 1’s at position $`j`$, and second do the child having all 0’s at position $`j`$. Similarly, if bit $`j`$ of $`\{\sigma _i\}`$ is a 0, then we first do the child having all 0’s at position $`j`$.

The essential advantage of the tree structure is that, typically, we do not have to check many structure strings for each sequence because nodes high up in the tree get eliminated by the upper bound. An additional advantage accrues because we are only interested in sequences with unique ground states. Therefore, as soon as two strings are found that satisfy the goal $`G`$, the search can be stopped for that sequence. Our protocol of following the branches that match the sequence $`\{\sigma _i\}`$ is intended to quickly identify strings which satisfy $`G`$.

We now discuss how the tree is actually built. In order to take advantage of the natural clustering of structure strings, we choose to split each node at the position that makes each of the two child nodes as tightly clustered as possible. We measure this clustering for each child as the sum of the “entropies” for each bit, with the tightest clustering corresponding to the minimum entropy. Specifically, for each child we evaluate the total entropy $`S`$ of its set of structure strings as

$$S=N_{\mathrm{child}}\underset{i}{}(p_i\mathrm{ln}p_i+q_i\mathrm{ln}q_i),$$
(2)
where $`p_i`$ is the probability of the $`i`$th position being a 1, and $`q_i`$ is the probability of the $`i`$th position being 0, averaged over all $`N_{\mathrm{child}}`$ structure strings in the child node. We choose to split each node at the position that minimizes the combined entropy of the two children.

The only remaining decision is when to stop splitting. When a node contains sufficiently few entries, it is fastest to just examine each entry in the node. Different stopping sizes were tried, and 16 seemed to be optimal to minimize search time per sequence. Another consideration was memory. The memory used by the entries themselves is 7,032,543 * 5 = 35,162,715 bytes = 33.5 megabytes, but each node takes additional space. With the given stopping criterion, there were 1,384,679 nodes in the tree, and the total space required was 70.5 megabytes.

Constructing a tree in this manner took a few hours on a 200 MHz Pentium Pro. Once the tree was constructed, carrying out each search took, on average, about 800 microseconds per sequence, but this time varied widely from sequence to sequence.

For comparison, we also implemented a naive search algorithm in which the energy of a sequence is computed for each structure string. The tree-search algorithm ran approximately 1600 times faster than our best variant of the naive approach.

## IV COMPUTING ALL THE GROUND STATES

The tree-search algorithm allows us to quickly compute ground states for each protein sequence in turn, and to record those sequences with unique ground states, together with the corresponding structure string and energy value. When all this is done, we would also like additional statistics, such as the designability of each structure, i.e., for how many sequences it is the unique ground state. This and other statistics can be computed, after the fact, if the unique ground-state solutions are stored.

Because our protein chains are directed, i.e. the two ends are not considered identical, both a structure and its oppositely directed partner are allowed. (Sometimes these are the same structure.) As a result, if a particular sequence $`\{\sigma _i\}`$ has a particular structure as a unique ground state, then the reversed sequence must have the reversed structure as its unique ground state. For the $`4\times 3\times 3`$ problem, there are therefore $`(2^{35})+(2^{17})=34,359,869,440`$ possible sequences that need to be considered, counting all 36-bit binary strings but rejecting reversed strings. There are 7,032,543 distinct structure strings. Since we are interested in unique ground states, each structure string is tagged with an additional bit to indicate whether it represents exactly one, or more than one, geometrical structure.

The overall computation is trivially parallelizable because calculating the ground state for each sequence can be done independently. In order to manage the calculation of ground states for all sequences in parallel, it proved useful to divide the space of sequences into “bundles”, and use these as the unit of parallelism. Instead of organizing these bundles by fixing some high-order bits of the 36-bit sequence $`\{\sigma _i\}`$ and varying the rest, we fixed an equal number of low- and high-order bits, varying the bits in the middle. This way entire bundles could be eliminated as reversals of the sequences in another bundle. We performed our computations with 14 bits fixed and 22 bits varying, which produced a total of $`(2^{13})+(2^6)=8,256`$ bundles.

These 8256 bundles were executed on the NECI Large Array Multiple Processors (LAMP) system, which is a collection of 28 computers, all containing 200MHz Intel Pentium Pro microprocessors. Three machines were uniprocessors, the rest had two processors each. Every machine ran the Linux operating system and had at least 128 MB of memory. Each computation typically required about 10 MB of memory for itself, plus about 70 MB for the (read-only) tree of protein structure information, which was shared in memory on the multiprocessors. The dual processors handled two independent bundles concurrently, and sharing the tree was necessary to keep the total memory “footprint” of the jobs small enough to fit together without conflict.

The distribution of bundles to “workers” was handled by a single “master” machine running shell and AWK scripts to poll the others, start new bundles, collect results, and detect any crashes that might occur. The scripts were written to be restartable with minimal lost effort in the event of a failure affecting their own operation. As each bundle was completed, a compactly coded binary output file was produced. At the end, all 8,256 output files were merged into a single 450 MB result file. Auxiliary programs were written to extract human-readable data from these binary files.

The complete computation ran for about 198 hours, with an average of 39 processors running at any one time, giving a total of 7,805 CPU hours. Bundle execution times varied from 23 seconds to almost 19 hours, with a mean of 56.7 minutes.

## V Results

Using the tree-searching algorithm and the LAMP system, we were able to completely enumerate the $`4\times 3\times 3`$ HP lattice protein model. We found that $`114,572,949`$ sequences have unique ground states, which is about $`0.3\%`$ of all sequences. In comparison, $`0.09\%`$ of the sequences in the $`3\times 3\times 3`$ model have unique ground states.

We associate to each structure a quantity called designability . The designability of a structure is the number of sequences having that structure as the unique ground state. By this definition, if two or more structures share the same string representation then they have zero designability since those structures can never be the unique ground state of any sequence. Thus, only the $`1,331,025`$ structures (and their reverse paths) each of which has its own, unshared string representation can have nonzero designability. For these structures, the average designability is $`114,572,949/1,331,025=86`$. However, the designability of these structures has a very broad range: from $`1`$ to $`4,466`$. (The minimum designability is 1 because the sequence with the same bit-string as the structure is guaranteed to have that structure as a unique ground state.) In Fig. 3, we plot the number of structures with a given designability versus the designability. One sees a long tail in the high designability region, consistent with previous results on the $`3\times 3\times 3`$ model  and on various two-dimensional models .

Also consistent with these previous works, there are noticeable geometrical differences between the highly designable $`4\times 3\times 3`$ structures and the less designable ones. In Fig. 4, we have plotted the average structure string $`s_i`$ for highly designable structures, and for all structures. Since $`s_i=1`$ for a core site, and $`s_i=0`$ for a surface site, the ensemble average $`s_i`$ gives the probability that the $`i`$th monomer on the chain occupies a core site. It is seen in Fig. 4 that for highly designable structures the first few monomers on the chain tend to occupy core sites, while there is no such tendency for the average compact structure. The average of $`s_i`$ for each chain is exactly 1/3 because every $`4\times 3\times 3`$ structure has exactly 12 core sites ($`s_i=1`$) and 24 surface sites ($`s_i=0`$). Therefore the tendency of the ends of highly designable structures to occupy core sites must be balanced by a tendency for the rest of the structure to occupy surface sites. This is also seen in Fig. 4 – the central third of the chain for highly designable structures has an increased probability to occupy surface sites, on average.

In Fig. 5 we have plotted the two-point correlation function of structure strings, $`C(i,j)=s_is_js_is_j`$, averaged over highly designable structures and over all structures. There is a clear correlation of site types, with a range of roughly one monomer in either direction along the chain. That is, if the $`i`$th monomer of a chain occupies a core site, there is an enhanced probability for monomers $`i\pm 1`$ to occupy core sites. This simply represents a general geometrical property of compact, self-avoiding structures . The correlation length is slightly shorter for highly designable structures (Fig. 5(b)), implying more frequent transitions between surface and core sites. In Fig. 6, we plot the number of transitions $`t`$ between surface and core sites versus designability $`N_S`$. For clarity, only structures with selected $`N_S`$’s are included in the plot. We see a weak positive correlation between $`t`$ and $`N_S`$, with a large variance for a given $`N_S`$. A much stronger positive correlation between surface-core transitions and designability was found in a two-dimensional $`6\times 6`$ lattice model , possibly reflecting the larger more compact core in the two-dimensional model.

As an example, the topmost designable $`4\times 3\times 3`$ structure is plotted in Fig. 7. The geometry of this structure is consistent with the results shown in Figs. 4 and 5. The 12 core sites are equally divided between the two ends of the chain, with the center part of the chain, $`i=926`$, consisting entirely of surface sites. Moreover, while the core sites tend to cluster, the longest run of core sites is only three, $`i=3335`$, consistent with the average correlation length shown in Fig. 5.

The complete enumeration of sequences and structures allows us to identify all sequences with a given structure as their unique ground state. We can therefore analyze the statistical properties of sequences that design a particular structure. For example, in Fig. 8, we have plotted the probability that a hydrophobic monomer occupies position $`i`$, averaged over all 4466 sequences that design the structure in Fig. 7. Fig. 8 therefore represents the complete mutation pattern of the topmost designable structure. The 12 core sites are easily identified since the probability of a hydrophobic monomer at these sites is nearly one. That is, nearly all of the 4466 sequences that design this structure have 1’s at these 12 positions. Similarly, the first three surface sites at the beginning of the chain and the last four surface sites at the end of the chain are always occupied by polar monomers. Interestingly, the monomers in the central parts of the chains, $`i=1026`$, have a roughly 1/3 chance of being hydrophobic even though all of these are surface sites. The mutation pattern therefore has the nontrivial feature that the monomers at some sites are critical for the stability of the ground state while the monomers at other sites are freely mutable.

## VI DISCUSSION

The existence of highly designable structures emerged from a study of $`3\times 3\times 3`$ and $`6\times 6`$ HP lattice protein models with interaction energies that included both solvation and segregation components. A later study verified the existence of highly designable structures for HP lattice models in two dimensions including only solvation energies . Moreover, in two dimensions, there was little change in the qualitative behavior of designability with increasing structure size, suggesting that highly designable structures persist up to realistic protein chain lengths. Within the $`3\times 3\times 3`$ solvation model, however, only 120 ($`0.1\%`$) out of the 103,346 total structures have nonzero designability. (These 120 (60 lattice conformations) stand out as highly designable structures with the largest average gaps in the solvation plus segregation model of Li et al. .) It has remained unclear whether the solvation model in three dimensions can produce a significant fraction of highly designable structures at realistic protein chain lengths. The current study, extending the solvation model up to $`4\times 3\times 3`$ structures, offers strong evidence that the existence of highly designable structures is a general feature of solvation models in three dimensions.

To understand the ubiquitous appearance of highly designable structures, it is helpful to review the geometrical interpretation of designability for the solvation model. To this end, the energy in Eq. (1) is rewritten as

$$E=\frac{1}{2}\underset{i=1}{\overset{N}{}}\left[|\sigma _is_i||\sigma _i||s_i|\right].$$
(3)
The last term is constant for a given sequence, and the second-to-last term is constant for all compact structures. Therefore, the ground state for a given sequence is determined by the first term alone, which is one-half the Hamming distance between the sequence string and the structure string. Simply put, the structure nearest to a sequence is its ground state. The designability of a structure is thus equal to the number of sequences that lie closer to it than to any other structure.

This geometrical interpretation suggests that highly designable structures are those with few nearby competing structures. To test this, we have plotted in Fig. 9 the number of neighboring structures, as a function of Hamming distance between structure strings, for the topmost designable structure and for structures of intermediate ($`N_S=100`$) and low ($`N_S=1`$) designability. It is seen that high designability implies a reduced number of neighbors, and this correlation persists out to distance 10 (structures whose strings differ by interchange of five surface and core sites). For comparison, we have also plotted in Fig. 9 the expected $`n(d)`$ if the $`14,062,236`$ structure strings were uniformly distributed on the hyperplane given by the constraints $`s_i=12`$ and $`(1)^i\times s_i=0`$. The second constraint comes from the fact that the $`4\times 3\times 3`$ cubic lattice is a bipartite lattice, so that for any structure there are six core sites with $`i`$ even and six core sites with $`i`$ odd. The total number of points in the hyperplane is $`(C_{18}^6)^2=344,622,096`$, where $`C_n^m=n!/m!(nm)!`$. The number of points in the plane at a Hamming distance $`d=2\times l`$ from a given point is simply: $`N(d)=_{k=0}^lC_6^kC_{12}^kC_6^{lk}C_{12}^{lk}`$, with $`C_n^m=0`$ if $`m>n`$. The expected $`n(d)`$ is then $`\rho \times N(d)`$, where $`\rho =14,062,236/344,622,0960.04`$ is the average density of real structure strings in the hyperplane. We see that the structure strings are more clustered than the uniform distribution – the lower the designability, the more clustered they are.

A scarcity of near neighbors corresponds to a narrow width of the distribution of neighbor distances, since every structure has the same total number of neighbors. We have therefore plotted in Fig. 10 the width of the distribution of neighbor distances over the entire range of designability. The width falls smoothly with increasing designability, indicating a general correlation between high designability of a structure and a scarcity of nearby competing structures.

The existence of highly designable structures can therefore be viewed as a geometrical property of the space of structure strings. For the two dimensional solvation model, it has been stressed  that the most highly designable structures fall in regions of low density in the space of structure strings. In this sense, highly designable structures have “atypical” patterns of surface and core sites. The current work extends this conclusion to a three dimensional solvation model with chain lengths ($`N=36`$) approaching those of real proteins.

In summary, we have employed a fast tree-search algorithm to find the ground states of all sequences for a $`4\times 3\times 3`$ lattice HP model of proteins. The results confirm the existence of highly designable structures in a three-dimensional solvation model with a surface-to-core ratio of 2-to-1, close to the values for small natural proteins. Highly designable structures are found to differ geometrically from other structures. Interestingly, structures are found to have nontrivial mutation patterns with some sites strictly conserved and others mutable. The fast tree-search algorithm is particularly well suited to lattice HP solvation models, and we hope our detailed description of the method will prove useful.


text:
# Autocorrelation functions in 3D Fully Frustrated Systems

## Abstract

We present a numerical study of autocorrelation functions of a 3D Fully Frustrated Ising model (FFIM) simulated by spin-flip Monte Carlo dynamics finding simple exponential decay for all the temperature above the critical temperature $`T_c`$ for the autocorrelation of squared magnetization and stretched exponential decay for the energy autocorrelation below a temperature $`T^{}`$ with $`T_c<T^{}T_p`$ where $`T_p`$ is the Kasteleyn-Fortuin and Coniglio-Klein percolation temperature. The results are compared to those on 2D FFIM to make light on the relevant mechanism in the onset of stretched exponential relaxation functions.

Autocorrelation functions of the form

$$f(t)=f_0\mathrm{exp}(t/\tau )^\beta $$
(1)
called stretched exponential or Kohlrausc-Williams-Watts (KWW) function are observed in glassy systems and spin glass (SG) well above the respective transition temperature and the origin of this kind of behaviour is still matter of debate.

For example in SG Randeria  *et al.*  suggested that the onset of the non exponential relaxation function coincides with the critical temperature of the ferromagnetic Ising model and base their conjecture on the presence of non frustrated ferromagnetic-type clusters of interactions.

Furthermore Huse and Fisher  on the basis of droplet model and more recently Stauffer  on the basis of percolation arguments have predicted a KWW autocorrelation function also for ferromagnetic Ising model (FIM) in 2D for $`TT_c`$ where $`T_c`$ is the critical temperature. These arguments show that no such behaviour is expected in more than two dimensions and Stauffer has corroborated his prediction with simulations in 2 and 3D.

It is useful to note at this point that in FIM the critical temperature $`T_c`$ coincides with the percolation temperature $`T_p`$ of the Kasteleyn-Fortuin and Coniglio-Klein (KF-CK) clusters.

These clusters can be obtained introducing bond with probability $`p=1e^{2J/(k_BT)}`$ between nearest-neighbor pair of spins satisfying the interaction ($`J`$ is the interaction strength, $`k_B`$ the Boltzmann constant and $`T`$ the temperature) and represent correlated spins regions.

It would be interesting to investigate whether such prediction also holds for frustrated systems where the percolation temperature and the critical temperature do not coincide.

In particular in this work we considered the fully frustrated Ising model (FFIM). Here we resume the numerical results on 2D  and present those in 3D.

The FFIM is a model with no disorder (unlike the SG and like the FIM) and with the interactions assigned in a regular way to give rise to *frustration*, *i.e.* a competition between them (like the SG and unlike the FIM).

The FFIM Hamiltonian is given by

$$H=J\underset{ij}{}(ϵ_{ij}S_iS_j1),$$
(2)
where $`ϵ_{ij}=\pm 1`$ are quenched regular assigned variables that give rice to frustration and $`S_i=\pm 1`$ are the Ising spins. In this model the KF-CK clusters are such that they do not include a  *frustrated loop*, that is a closed path of bonds which contains an odd number of interactions with $`ϵ_{ij}=1`$. These clusters overlap regions of not frustrated interactions and their percolation temperature is greater then the critical temperature. In Ref. it is studied a general Hamiltonian, with Ising spins and $`s`$-state Potts spins, that has the FFIM as a particular case for $`s=1`$. It is possible to see that for all $`s1`$ $`T_p`$ of KF-CK clusters corresponds to a thermodynamic transition of a $`s`$-state Potts model and that this thermal transition disappears for $`s=1`$, *i.e.* for FFMI. This suggests that the KF-CK clusters may have an important role also for FFIM.

We simulate this model using the standard spin flip  dynamics.

In the 2D FFIM where $`T_p>T_c=0`$ recent results  shows that stretched exponential autocorrelation functions are observed for $`TT_p`$. Therefore the same qualitatively behaviour of 2D FIM is observed but in the FFIM is clear that the onset of stretched exponentials coincides, within numerical precision, with $`T_pT_c`$, while in FIM is $`T_p=T_c`$.

To put more light on the problem we have done simulations on 3D FFIM.

To this goal we have calculated the autocorrelation function of squared magnetization $`M^2`$ and of energy $`E`$, where for a generic observable $`A`$ the relaxation function is defined as

$$f_A(t)=\frac{\delta A(t)\delta A(0)}{(\delta A)^2},$$
(3)
where $`\delta A(t)=A(t)A`$, $`t`$ is the time and $``$ is the thermal average.

We fit our data with a general function

$$f(t)=f_0\frac{\mathrm{exp}(t/\tau )^\beta }{t^x}$$
(4)
where $`f_0`$, $`\tau `$, $`\beta `$ and $`x`$ are parameters.

In Fig.1 we summarize the results for systems with linear size $`L=20`$ and for temperatures below the estimated  $`T_p(L)`$ of KF-CK clusters in cubic FFIM.

For $`f_{M^2}`$ the estimated parameters are $`\beta 1`$ and $`x0`$ for all simulated temperatures below $`T_p(L)`$. For $`f_E`$ there is a temperature $`T^{}`$, between $`T_p(L)`$ and $`T_c1.355J/k_B`$ (the estimated critical temperature ), below which is $`\beta <1`$ and $`x`$ weakly changes. Simulations on $`L=30÷50`$ confirm these data.

Therefore the 3D FFIM results show a simple exponential decay in $`M^2`$ correlation functions but a stretched exponential decay in $`E`$ correlation functions for $`TT^{}`$ with $`T_c<T^{}T_p`$. In conclusion the results on $`f_E`$ in 3D FFIM, even if less clear than those on 2D FFIM, confirm the existence of an onset $`T^{}`$ of stretched exponential for the relaxation function well above the critical temperature $`T_c`$ and almost equal (or slightly below) to the KF-CK percolation temperature $`T_p`$ that, in a generalized model with also $`s`$-state Potts variables with $`s1`$, corresponds to a thermodynamic transition of the Potts variables.