format: text

text:
# Isotopic Scaling in Nuclear Reactions

## Abstract

A three parameter scaling relationship between isotopic distributions for elements with Z$`8`$ has been observed that allows a simple description of the dependence of such distributions on the overall isospin of the system. This scaling law (termed iso-scaling) applies for a variety of reaction mechanisms that are dominated by phase space, including evaporation, multifragmentation and deeply inelastic scattering. The origins of this scaling behavior for the various reaction mechanisms are explained. For multifragmentation processes, the systematics is influenced by the density dependence of the asymmetry term of the equation of state.

The availability of high intensity radioactive beams facilitates the exploration of the isospin degree of freedom in nuclear reactions. Understanding the connection between the entrance channel isospin and the isotopic distribution of reaction products is important for studying the charge asymmetry term of the nuclear equation-of-state\[1-3\], obtaining information about charge equilibration\[4-6\], providing stringent tests for reaction models and optimizing the production of rare isotopes far from stability. In this letter, we demonstrate that isotopic distributions for statistical production mechanisms follow scaling laws. We also find circumstances where the values for the scaling parameters are influenced by the density dependence of the asymmetry term of the nuclear equation of state, a quantity that influences many important properties of neutron stars.

The scaling laws in question relate ratios of isotope yields measured in two different nuclear reactions, 1 and 2, $`R_{21}(N,Z)=Y_2(N,Z)/Y_1(N,Z)`$. In multifragmentation events, such ratios were shown to obey an exponential dependence on the neutron and proton number of the isotopes characterized by three parameters $`\alpha ,\beta `$ and $`C`$ :

$$R_{21}(N,Z)=C\mathrm{exp}(\alpha N+\beta Z)$$
(1)
Here we choose the convention that the isospin composition (neutron to proton ratio) of system 2 is larger than that of system 1. The systematics described by Eq. 1 occur naturally within the grand-canonical ensemble \[7-9\]. As shown in Ref. , the parameters $`\alpha `$ and $`\beta `$ in that limit are the differences between the neutron and proton chemical potentials for the two reactions (i.e., $`\alpha =\mathrm{\Delta }\mu _n/T`$ and $`\beta =\mathrm{\Delta }\mu _p/T`$ ) and $`C`$ is an overall normalization constant.

The accuracy of the iso-scaling described by Eq. 1 can be compactly displayed if one plots the scaled isotopic ratio,

$$S(N)=R_{21}(N,Z)\mathrm{exp}(\beta Z)$$
(2)
as a function of $`N`$. For all elements, $`S(N)`$ must lie along a straight line on a semi-log plot when Eq. 1 accurately describes the experimental data. The data points marked as ”multifragmentation” in Figure 1 show values of $`S(N)`$ extracted from isotope yields with $`1Z8`$ measured for multifragmentation events in central $`{}_{}{}^{124}Sn+^{124}Sn`$ and $`{}_{}{}^{112}Sn+^{112}Sn`$ collisions at E/A= 50 MeV . Selection of central events ensures that the average excitation energies and temperatures in the participant source should be nearly identical . The observed iso-scaling is a necessary condition for the applicability of equilibrium models; such models have described other aspects of these collisions quite well .

Rather surprisingly, iso-scaling is also observed for strongly damped binary collisions ($`{}_{}{}^{16}O`$ induced reactions on two targets $`{}_{}{}^{232}Th`$ and $`{}_{}{}^{197}Au`$)  and evaporative compound nuclear decay ($`{}_{}{}^{4}He+^{116}Sn`$ and $`{}_{}{}^{4}He+^{124}Sn`$ collisions), for which Grand-Canonical Ensemble approaches would appear to have little relevance. Our studies suggest that iso-scaling is obeyed where 1.) both reactions 1 and 2 are accurately described by a specific statistical fragment emission mechanism and 2.) both systems are at nearly the same temperature. Indeed for deeply inelastic reactions emitted forward of the grazing angles, ($`{}_{}{}^{22}Ne+^{232}Th`$ and $`{}_{}{}^{22}Ne+^{97}Zr`$ at $`\theta =12^{}`$ and $`E/A=7.9MeV`$ ) and for reactions with different temperatures ($`{}_{}{}^{124}Sn+^{124}Sn`$  and $`{}_{}{}^{4}He+^{124}Sn`$ ), iso-scaling is violated. Conditions 1 and 2 are met by the three reactions shown in Figure 1. Why iso-scaling is specifically observed in these cases and what aspects of statistical physics such scaling probes are examined below.

We first examine the strongly damped collisions, where iso-scaling is reasonably well respected at low incident energies ($`E/A<10MeV`$) and at relatively backward angles , i.e., when equilibrium is established between the orbiting projectile and target. In such cases, the isotopic yields follow the ”$`Q_{gg}`$-systematics”, and can be approximated by

$$Y(N,Z)\mathrm{exp}((M_P+M_TM_P^{}M_T^{})/T)$$
(3)
where $`M_P`$ and $`M_T`$ are the initial projectile and target masses, and $`M_P^{}`$ and $`M_T^{}`$ are the final masses of the projectile- and target-like fragment. Here, $`T`$ can be interpreted as the temperature. Eq. 3 reproduces the systematics shown in Fig. 1. To show why this is so, we have expanded the nuclear binding energy contributions to the masses in Taylor series in $`N`$ and $`Z`$. Expressing explicitly only the leading order terms that depend on $`N`$ and $`Z`$, we obtain a relatively accurate leading order approximation to Eq. 3 :

$$R_{21}\mathrm{exp}[(\mathrm{\Delta }s_nN\mathrm{\Delta }s_pZ)/T].$$
(4)
where $`\mathrm{\Delta }s_n`$ and $`\mathrm{\Delta }s_p`$ are the differences of the neutron and proton separation energies for the two compound systems. Thus, the difference in the average separation energies in Eq. 4 plays a corresponding role to the difference in chemical potentials in the grand canonical approach, an intriguing result when one considers that $`\mu s`$ in the low temperature limit . The straightforward dependence of Eq. 4 on temperature suggests that it may provide information relevant to the temperatures achieved in strongly damped collisions. The data of Figure 1 imply a temperature of 2.7 MeV, not inconsistent with values derived from alternative analyses .

Next we consider the yields from processes involving the formation of a composite system and its subsequent decay via the evaporation of different isotopes. Corresponding scaled isotopic ratios for fragments detected at backward laboratory angles ($`\theta =160^o`$) in $`{}_{}{}^{4}He+^{116}Sn`$ and $`{}_{}{}^{4}He+^{124}Sn`$ collisions at $`E/A=50MeV`$  are shown in Figure 1, next to the label ”evaporation”.

To explore the factors that govern the relevant evaporation rates, we utilize the formalism of Friedman and Lynch  which invokes statistical decay rates derived from detailed balance . When the yields are dominated by emission within a particular window of source-mass or source-temperature, the relative yields of a fragment with neutron number $`N`$ and proton number $`Z`$ are directly related to the instantaneous rates

$`dn(N,Z)/dt`$ $``$ $`T^2\mathrm{exp}(V_c/T+Nf_n^{}/T+Zf_p^{}/T`$ (6)
$`\{BE(N_i,Z_i)BE(N_iN,Z_iZ)BE(N,Z)\}/T)`$
where $`V_c`$ gives the Coulomb barrier, the terms $`f_n^{}`$ ($`f_p^{}`$) represent the excitation contribution to the free energy per neutron (proton), BE is the binding energy and $`N_i`$ and $`Z_i`$ identify the neutron and proton numbers of the parent nucleus.

Applying Eq. 6 to the calculations of $`R_{21}`$ for two systems at the same temperature, we find that the binding energies of the emitted fragments cancel and the systematics shown in Fig. 1 can be reproduced. To understand why, we again expand the binding energy of the residue with neutron number $`N_iN`$ and proton number $`Z_iZ`$ to leading order in a Taylor series to obtain:

$$R_{21}(N,Z)exp[\{(\mathrm{\Delta }s_n+\mathrm{\Delta }f_n^{})N+(\mathrm{\Delta }s_p+\mathrm{\Delta }f_p^{}+e\mathrm{\Delta }\mathrm{\Phi }(Z_iZ))Z\}/T]$$
(7)
where $`\mathrm{\Phi }(Z)`$ is the electrostatic potential at the surface of a nucleus with neutron and proton number $`N`$ and $`Z.`$ Aside from the second order term from the electrostatic potential, which is small for the decay of large nuclei, all factors in the exponent are proportional to either $`N`$ or $`Z`$, consistent with Eq. (1). The corresponding scaling parameters $`\alpha `$ and $`\beta ,`$ are functions of the separation energies, the Coulomb potential and small contributions from the free excitation energies. Using the functions of $`\mathrm{\Delta }s_n`$ and $`\mathrm{\Delta }f_n^{}`$ from Friedman and Lynch , one finds that a fixed temperature of about 3.7 MeV is required in Eq. 7, to obtain the experimental value of $`\alpha =0.6.`$ Running a full evaporation chain, using the procedure of Ref. , provides an average fragment emission temperature of about 3.3 MeV. These temperature values are comparable to those extracted by other techniques .

The Expanding Evaporating Source (EES) model  provides an alternative description of multifragmentation. Within the context of that model, additional insights can be obtained. The EES model utilizes a formula for the particle emission rates which is formally identical to that of Eq. 6 but can differ significantly in its predictions because the residue may expand to sub-saturation density . In this circumstance, the term enclosed in brackets ”{}”containing three binding energies in Eq. 6 may vanish or become negative, enhancing the emission rate of fragments with $`3Z20.`$ Detailed examination reveals that $`\mathrm{\Delta }f_n^{}`$ in Eq. 7 is usually much smaller than $`\mathrm{\Delta }s_n,`$ and the volume, surface, and Coulomb contributions to $`\mathrm{\Delta }s_n`$ largely cancel, leaving the asymmetry energy term, $`Sym(\rho )(NZ)^2/A,`$ alone as the dominant contribution to $`\alpha `$. For simplicity, we assume a power law dependence for $`Sym(\rho ),`$ i.e. $`Sym(\rho )=C_{sym}(\rho /\rho _0)^\gamma `$ where $`\gamma `$ is a variable and $`C_{sym}=23.4`$ MeV is the conventional liquid drop model constant .

For illustration, we have performed calculations for the decay of the composite systems found in $`{}_{}{}^{124}Sn`$ \+ $`{}_{}{}^{124}Sn`$ and $`{}_{}{}^{112}Sn`$ \+ $`{}_{}{}^{112}Sn`$ collisions assuming, for simplicity, initial systems of ($`Z_{tot},`$ $`A_{tot})`$ of (100, 248$`)`$ and (100, 224), respectively, initial thermal excitation energies of $`E_{ther}^{}=9.5MeV,`$ and initial collective radial expansion energies of $`E_{Coll}/A=2.5MeV.`$ The fragments are emitted from these systems as they expand from an initial density $`\rho /\rho _0=1`$ to $`\rho /\rho _0=0.1.`$ The left panel of Figure 2 shows two iso-scaling functions, $`S(N)`$, calculated with the predicted isotope yields for $`3Z6`$ for $`\gamma =1`$ and $`\gamma =0`$. In spite of a rather complex interplay of expansion and fragment emission, the EES model predictions still exhibit an approximate scaling. Some of the calculated deviations from iso-scaling may, in fact, be due to the rather schematic treatment of the Coulomb barrier penetration  and to the incomplete sequential decay information used in the model. Values for $`\alpha ,`$ shown in the right panel, are determined from the slopes of best fit of the lines to the predicted iso-scaling functions.

For large values of $`\gamma ,`$ the asymmetry term, $`C_{sym}(\rho /\rho _0)^\gamma `$ decreases more rapidly with density and becomes negligible as the residue expands. At low density where fragments are predominately emitted , the difference in isotopic yields from the two reactions 1 and 2 become smaller, resulting in flatter scaling functions and smaller values of $`\alpha `$. (In the extreme case of identical yields from the two systems, $`S(N)`$ becomes a horizontal line corresponding to $`\alpha =0`$.) The dot-dash line in the right panel of Fig. 2 joining the solid points shows the EES prediction that $`\alpha `$ decreases with increasing $`\gamma `$ values. The multifragmentation data in Fig. 1 can be fairly well reproduced by $`\gamma `$ $``$0.6.

When the emission process ends at $`\rho =0.1\rho _0`$, a low density residue (LDR) may remain. The final fate of this residue is not predicted by the EES model. A potential ambiguity may result if this residue is large and if its eventual disintegration produces a significant fraction of fragments with $`Z3`$. Calculations using the Statistical Multifragmentation Model  with the mass and energy of the final LDR predict $`\alpha `$-values that increase with $`\gamma `$ (dashed line) — opposite to the trend predicted for emission during the expansion. When $`\gamma `$ decreases, the $`N/Z`$ of the LDR from the two reactions, $`{}_{}{}^{124}Sn`$ \+ $`{}_{}{}^{124}Sn`$ and $`{}_{}{}^{112}Sn`$ \+ $`{}_{}{}^{112}Sn`$ become more similar, resulting in smaller $`\alpha `$-values – a trend also predicted from isospin dependent transport theory . This accounts for the behavior of the dashed line in the right panel of Figure 2. If high energy fragments are emitted primarily from the expanding system, and low-energy fragments come from the instantaneous disintegration of a low-density residue, then the results in Fig. 2 suggest that the predicted difference in the iso-scaling for low and high energy fragments should be observable for $`\gamma <0.8`$.

In summary, we have observed a scaling between isotopic distributions which allows a simple description of the dependence of such distributions on the overall isospin of the measured systems in terms of three parameters, $`\alpha ,\beta `$ and $`C.`$ This scaling seems to apply to a broad range of statistical fragment production mechanisms, including evaporation, strongly damped binary collision, and multifragmentation. We have shown how this systematics arises within models frequently applied to such processes. In one such model, the EES model of multifragmentation, we find that the iso-scaling parameters are sensitive to the density dependence of the asymmetry term of the EOS.

This work was supported by the National Science Foundation under Grant Nos. PHY-95-28844 and PHY-96-05140.

Figure Caption:

Figure 1: The scaled isotopic ratio, $`S(N)`$ is plotted as a function of the neutron number, N, using the best fit value of $`\beta `$ obtained from fitting isotopes with $`Z3`$. The data points plotted next to the label ”multifragmentation” in Figure 1 denote $`S(N)`$ extracted from multifragmentation events in central $`{}_{}{}^{124}Sn+^{124}Sn`$ and $`{}_{}{}^{112}Sn+^{112}Sn`$ collisions with $`\alpha =0.37,\beta =0.40`$,. The scaling behavior for evaporation process is illustrated by the reactions $`{}_{}{}^{4}He+^{116}Sn`$ and $`{}_{}{}^{4}He+^{124}Sn`$  plotted next to the label ”evaporation” with $`\alpha =0.60,\beta =0.82`$. Systematics of the strongly damped binary collisions is represented by the data of $`{}_{}{}^{16}O`$ induced reactions on two targets $`{}_{}{}^{232}Th`$ and $`{}_{}{}^{197}Au`$  plotted next to the label ”deeply inelastic” with $`\alpha =0.74,\beta =1.1`$.

Figure 2: Theoretical EES model predictions for the scaling functions $`S(N)`$ (left panel) and scaling parameter $`\alpha `$ (right panel) of fragments emitted in $`{}_{}{}^{124}Sn`$ \+ $`{}_{}{}^{124}Sn`$ and $`{}_{}{}^{112}Sn`$ \+ $`{}_{}{}^{112}Sn`$ collisions. Same convention as Figure 1 applies to the various symbols used in the left hand panel.

References:

1. Bao-An Li et al., Phys. Rev. Lett. 78, 1644 (1997).

2. I. Bombaci, et al., Phys. Rep. 242, 165 (1994).

3. J.M. Lattimer and M. Prakash, Ap. J. (in press).

4. R. Laforest et al., Phys. Rev. C59 2567 (1999) and refs. therein.

5. R. Wada et. al., Phys. Rev. Lett. 58, 1829 (1987).

6. W.U. Schroeder and U.J. Huizenga, Treatise on heavy Ion Science, Ed. D.A. Bromley (Plenum Press, 1984) and references therein.

7. H. Xu et al., Phys. Rev. Lett. 85, 716 (2000).

8. J. Randrup and S.E. Koonin, Nucl. Phys. A 356, 223 (1981).

9. S. Albergo et al., Nuovo Cimento A 89, 1 (1985).

10. G. J. Kunde et al., Phys. Lett. B416, 56 (1998).

11. M. D’Agostino et al., Phys. Lett. B 371, 175 (1996); B.A. Li et al., Phys.Lett.B303, 225 (1993); A.S. Botvina et al., Nucl. Phys. A 584, 737 (1995).

12. V.V. Volkov, Phys. Rep. 44, 93, (1978).

13. J. Brzychczyk et al., Phys. Rev. C47, 1553 (1993).

14. C.K. Gelbke et. al., Phys. Rep. 42, 311 (1978).

15. S.R. Souza, W.P. Tan, R. Donangelo, C.K. Gelbke, W.G. Lynch, M.B. Tsang, Phys. Rev. C62, 064607 (2000).

16. W. Friedman and W. Lynch, Phys. Rev. C 28, 950 (1983).

17. This approach is similar to the Weisskopf model of V. Weisskopf, Phys. Rev, 52, 295 (1937).

18. V. E. Viola, K. Kwiatkowski, W. A. Friedman, Phys. Rev. C 59, 2660 (1999)

19. W.A. Friedman, Phys. Rev. Lett. 60, (1988) 2125; and Phys. Rev. C42, 667 (1990).

20. Note: $`f^{}(T/ϵ_f)/T`$ remains roughly constant during expansion.

21. A. Bohr & B. R. Mottleson, ”Nuclear Structure, Vol II”, W.A. Benjamin Inc., (1998).

22. J. Tõke, W. Gawlikowicz, and W. U. Schröder, Phys. Rev. C 63, 024606 (2001).

23. A. Majumder and S. Das Gupta, Das Gupta, Phys. Rev. C 61, 034603 (2000), Scott Pratt, Subal Das Gupta, Phys. Rev. C 62, 044603 (2000)

24. W. Tan et al, private communications and to be published.


text:
# Excluded Volume in Protein Sidechain Packing

## Introduction

One remarkable characteristic of proteins is that sidechains comprising the hydrophobic core are as closely packed as organic crystals . In most proteins with known 3D structures, the core residues are rarely disordered and adopt one of a small number of alternative conformations. This almost unique packing is effected by a combination of steric interactions (excluded volume effects) and energetic stabilization (hydrophobic, polar, and charge interactions). The proportions by which these interactions contribute to the overall stability is unknown, but several studies suggest that steric and hydrophobic interactions are of primary importance .

Sidechain packing has been studied intensively for various reasons, notably: 1) It is thought to be a crucial piece of the protein folding puzzle ; 2) The selection of protein sequences through evolution may have been influenced by how well a sequence can be packed for a particular fold ; 3) Accurate packing algorithms are necessary for complete protein structure prediction; and 4) Existing threading and homology modeling algorithms may be significantly improved by a better understanding of how sidechains are stabilized in the core . A plethora of computational methods for modeling protein sidechains and energetics exist , and yet which ones best capture the underlying physics is unclear. This is due, in part, to the use of energy functions containing many different types of interactions.

In this paper, we examine the role of excluded volume in isolation of all other types of interactions. This allows us to address the importance of geometry in determining the native state of protein sidechains. We use an all-atom, rotamer-based model of a protein to obtain repacked conformations of its interior sidechains. By choosing the most realistic representation of sidechain and backbone geometries, we eliminate errors encountered in coarse-grained models.

Because there is no consensus on how excluded volume interactions should be modeled in the context of proteins, we consider three different models. The simplest model treats all heavy atoms as hard spheres which are not allowed to overlap. The second model further restricts the sterically allowed space by adding a second shell around the hard spheres which tolerates only a limited number of overlaps throughout the entire protein. The third model uses a continuous $`r^{12}`$ potential, which is the repulsive contribution from a Lennard-Jones potential. Since a potential which decays faster than our chosen $`r^{12}`$ potential would allow more sidechain conformations, and since a slower decaying function is generally not used to model sterics , the three we have chosen should cover most of the possibilities as far as models featuring spherically symmetric, atom-based potentials.

For each model, we obtain the distribution of rms values for sidechain conformations satisfying excluded volume constraints. This amounts to a full characterization of the ground-state ensemble of protein sidechains subject to excluded volume effects only. We find a vast number of significantly different conformations with native-like packing density that meet excluded volume constraints, implying that excluded volume alone cannot stabilize the native conformation of sidechains even when the backbone is fixed. We provide an estimate for the amount of additional stabilization that must be provided through interactions other than sterics. The results obtained using three different models are consistent with each other, suggesting that our main conclusions are model-independent. We also examine the correlation between our simulations and the distribution of rotamers in the pdb database, and discuss implications for sidechain prediction. Our main goal is to elucidate the role of excluded volume in establishing the ensemble of conformations compatible with the folded state of a protein.

## Methods

All-Atom Protein Representation. Three proteins were used in this study: photoactive yellow protein (1F9I), subtilisin (1GCI), and concanavalin (1NLS). We chose these proteins because they have very high resolution crystal structures, and are of different secondary-structure classes: subtilisin is a 0.78 Å resolution structure, and is an all $`\alpha `$-helical protein; concanavalin (0.94 Å resolution) is all $`\beta `$-sheet; and photoactive yellow protein (1.1 Å resolution) is $`\alpha `$/$`\beta `$. All heavy atoms of the proteins present in the crystal structure were represented. Each sidechain torsion was allowed to take on one of three rotamer values  plus a random noise of $`\pm 10^{}`$. Rotamer probabilities  were not used in the simulation. At all times, the atoms of the entire backbone were kept fixed in their original crystal structure positions.

Models of sterics. The excluded volume interactions of protein sidechains were modeled in three ways:

1) In the *hard-sphere model* atoms were treated as impenetrable spheres of given radii. These hard-sphere radii are necessarily smaller than the van der Waals (vdW) radii, and were defined by scaling the vdW radii by some factor $`\alpha _h`$. All atom-atom interactions at distances smaller than the sum of the hard-sphere radii (“hard clashes”) were strictly forbidden. A pair of atoms $`i`$ and $`j`$ separated by a distance $`r_{ij}`$ is said to be a hard clash if $`r_{ij}<\alpha _h\left(r_0(i)+r_0(j)\right)`$, where $`r_0(i)`$ is the vdW radius of atom $`i`$. The set of vdW radii determined in  was used. $`\alpha _h`$ was taken to be the largest value such that all hard clashes in the native protein conformation can be eliminated within the $`\pm 10^{}`$ allowed at each torsion; we found that $`\alpha _h`$$`=0.76`$ for concanavalin and photoactive yellow protein, and $`0.77`$ for subtilisin. The steric energy, $`U_h`$, of a conformation in this model is given by the number of hard clashes. Only conformations with $`U_h=0`$ are sterically allowed in this model.

2) Since results could depend on the choice of $`\alpha _h`$, we explored a *soft-sphere model* in which atoms consisted of two radii - the hard-sphere radii fixed by $`\alpha _h`$, and somewhat larger, soft-sphere radii. Soft radii were defined by scaling vdW radii by a parameter $`\alpha _s`$, where $`\alpha _h<\alpha _s1`$. A soft clash occurs when $`\alpha _h\left(r_0(i)+r_0(j)\right)<r_{ij}<\alpha _s\left(r_0(i)+r_0(j)\right)`$. All hard clashes were forbidden as before, while only a limited number of soft clashes are allowed over the entire protein. The steric energy, $`U`$, of a conformation is given by $`U=U_h+U_s`$, where $`U_h`$ is the number of hard atom-atom clashes, and $`U_s`$ is the number of soft atom-atom clashes above threshold. The threshold number of soft clashes allowed is a function of $`\alpha _s`$, as described in Fig. 1. Only conformations with $`U=0`$ are sterically allowed in this model. Since van der Waals radii are strictly greater than excluded volume radii, as demonstrated in numerous studies of gas phase organic molecules , $`\alpha _s`$ must be strictly less than 1.0, and we include $`\alpha _s`$$`=1.0`$ in our plots only as an upper limit on the radii.

3) The third model used the repulsive term of a *Lennard-Jones* potential. For each pair of atoms we fit the LJ potential, $`A/r^{12}1/r^6`$, so that the minimum coincides with the sum of the vdW radii taken from . With this model, the steric energy of a conformation is given by $`U_{LJ}=_{ij}A_{ij}/r_{ij}^{12}`$ where the sum is taken over all pairs of atoms, and $`A_{ij}=(r_0(i)+r_0(j))^6/2`$. We considered a conformation to be sterically allowed if its energy was less than or equal to the energy of the native crystal structure.

Monte Carlo Packing Algorithm. Initial random sidechain conformations were packed to a given rms interval $`[R_{min},R_{max}]`$ by a Metropolis Monte Carlo simulation . At each step of the simulation, a random position was selected and changed to another rotamer. The move was accepted with probability $`p`$ given by

$$p=\{\begin{array}{cc}1\hfill & \delta H0\hfill \\ \mathrm{exp}(\delta H/T)\hfill & \mathrm{otherwise}\hfill \end{array}$$

with

$$\delta H=\delta U+\{\begin{array}{cc}C\left(R_{min}R\right)\hfill & R<R_{min}\hfill \\ 0\hfill & R_{min}RR_{max}\hfill \\ C\left(RR_{max}\right)\hfill & R>R_{max}\hfill \end{array}$$

where $`U`$ is the steric energy, $`R`$ is the rms, and $`C=10`$. All reported rms values were computed over the subset of sidechains having more than 40% occluded surface area as defined in reference 15, and at least 1 rotatable $`\chi `$ angle. All heavy atoms of a sidechain, including the $`C_\beta `$ atom, entered into the rms calculation as in reference 3. Prolines were kept fixed in their crystal structure positions throughout this work. The 40% occluded positions for each protein are given below. Additionally, by visual inspection we found positions whose $`C_\alpha `$-$`C_\beta `$ bonds were pointing toward the core, and/or whose sidechains were surrounded by other atoms on all sides. These positions are indicated in bold type.
Subtilisin:
(1) 8I,11V,17H,21L,22T,26V,27K,28V,30V
(10) 31L,32D,33T,35I,36S,38H,40D,41L,49F,50V
(20) 58D,62H,64T,65H,66V,69T,70I,73L,79V,80L
(30) 82V,87E,88L,89Y,91V,92K,93V,94L,104S,105I
(40) 109L,111W,117M,119V,121N,122L,123S,130S,133L,137V
(50) 139S,141T,145V,146L,147V,148V,151S,159I,160S,161Y
(60) 165Y,169M,171V,174T,175D,183F,184S,185Q,186Y,190L
(70) 191D,192I,193V,197V,199V,201S,202T,203Y,214T,215S
(80) 216M,218T,220H,221V,227L,228V,229K,232N,237N,240I
(90) 241R,243H,244L,245K,256L,257Y,261L,262V,263N,265E
(100)268T
Concanavalin:
(1) 4I,5V,7V,8E,9L,10D,11T,14N,17I
(10) 19D,24H,25I,27I,28D,29I,31S,32V,40W,52I
(20) 54Y,55N,61L,62S,65V,67Y,75V,79V,80D,81L
(30) 85L,89V,90R,91V,93L,94S,96S,97T,102E,103T
(40) 104N,105T,106I,108S,109W,110S,111F,113S,115L,117S
(50) 128F,130F,133F,134S,140L,143Q,145D,148T,153N,154L
(60) 156L,157T,164S,169S,170V,172R,174L,175F,179V,180H
(70) 191F,195F,196T,197F,199I,201S,208D,210I,211A,213F
(80) 214I,215S,216N,219S,225S,229L,230L,232L,233F
Photoactive Yellow Protein:
(1) 6F,11I,12E,15L,18M,22Q,23L,26L,28F
(10) 31I,32Q,33L,34D,38N,39I,41Q,42F,43N,46E
(20) 49I,50T,57V,61N,62F,63F,66V,70T,72S,75F
(30) 76Y,79F,83V,88L,90T,92F,96F,105V,106K,107V
(40) 108H,109M,110K,111K,117S,118Y,119W,120V,121F,122V
Umbrella sampling. Umbrella sampling  was used to determine the number of states $`\mathrm{\Omega }(R)`$ for each of the 3 models. First, different sterically allowed conformations were collected for rms intervals of size 0.05 Å over the entire rms range ($`04.0`$ Å). The rms values of conformations differing only by random noise were averaged. Next, the probability to observe a conformation with rms $`R`$, $`p_i(R)`$, was obtained for the $`i`$-th rms interval, $`[R_i,R_{i+1}]`$. Since the distribution of conformations $`\mathrm{\Omega }(R)`$ must be a continuous function of $`R`$, $`\mathrm{\Omega }(R)`$ can be determined by appropriately scaling the probabilities $`p_i(R)`$ to ensure continuity at interval boundaries. Specifically, the probabilities $`p_i(R)`$ are sequentially scaled by $`p_{i1}(R_i)/p_i(R_i)`$, starting from $`i=2`$. $`\mathrm{\Omega }(R)`$ is then obtained by multiplying all rescaled $`p_i`$’s by $`\mathrm{\Omega }(R_1)/p_1(R_1)`$.

The umbrella sampling scheme can only work provided that the conformational space within each bin is explored evenly during simulation. In order to move the system into a particular rms bin, we had to use $`C=10`$ and work at very low temperature ($`T=0.005`$). These conditions were necessary due to the large amount of entropy available just outside any given bin, both in rms and in steric energy (the number of clashed conformations is exponentially larger than the number of well-packed ones). The downside is that the system may end up in a part of conformational space within a given bin that is disconnected from other parts of the bin, i.e. the system behaves non-ergodically under the conditions in which it must be sampled. The rms histograms within such isolated clusters are not representative of the entire space within a bin, and can therefore lead to large errors when the full rms histogram is assembled. To overcome this difficulty, we performed 200 short runs, starting from random conformations, for each bin and computed the size of the cluster explored in each run. The cluster size is defined as the average torsional distance between all pairs of uncorrelated conformations encountered in a run. We found that most bins had many clusters of various sizes. Since the number of states within a cluster should scale exponentially with its size, the largest cluster encountered should dominate the statistics within a given bin. We therefore sampled each bin by starting a run from a conformation obtained in a short run that had entered the large cluster. For every rms bin, we collected 10,000 states in the largest cluster, and used only these states when assembling the full rms histogram.

Several tricks proved handy when trying to find the large cluster. The plot of the largest cluster’s size over rms should be relatively smooth. If we found a cluster of size 20 at rms 0.90, a cluster of size 10 at rms 0.95, and a cluster of size 24 at rms 1.00, we could be certain that the 0.95 bin would be badly sampled. We could thus focus our efforts on bins with anomalously low cluster sizes. Some bins required up to 500 short runs before the large cluster was found. We observed that there is only one cluster at low rms values. If we started a run from the lowest rms bin, and allowed it to escape into a higher rms bin, the run would always end up in the largest cluster for that bin. This worked up to rms 1.2 because entropy was steeply increasing with rms, and one cluster was dominating the statistics. For rms values between 1.2 and 2.0, several clusters of similar sizes emerged, and thus sampling the largest one was no longer adequate. In this rms range, the distribution of states was relatively flat over rms, and thus the entire range could be fully sampled as a single large bin. For rms higher than two, the method of multiple short runs worked well. The distribution of clusters as rms is varied emerged as an interesting problem in itself, but was beyond the scope of this paper.

CHARMM minimization of structures Using CHARMM, we minimized and added hydrogen atoms to a set of 50 randomly chosen repacked conformations of each of the three proteins used in this study. All such minimizations, including those performed for other randomly selected packed structures, did not change the rotamer states of the sidechains, but did eliminate close contacts (as defined by CHARMM) resulting from the addition of hydrogens. A set of 1000 repacked conformations are available for public download at http://www-shakh.harvard.edu/$``$shimada/index.html.

## Results

Figure 2 shows the distribution of *packed* (sterically allowed) conformations for three proteins, obtained over a range of rms values using umbrella sampling. For each protein, the rms values are computed over the subset of partially occluded residues (see Methods), and thus only the occluded positions contribute to the total number of states. All other positions are free to repack, but are not counted in our sampling. With $`\alpha _s`$$`=`$ $`\alpha _h`$ (the hard-sphere model), the number of packed conformations of subtilisin is approximately $`10^{32}`$. As $`\alpha _s`$ is increased using the soft-sphere model, the number of packed conformations decreases. This trend reverses when the soft radii reach the van der Waals size (see Discussion), at $`\alpha _s`$$`=1.0`$, and the number of repackings is approximately $`10^{19}`$. We emphasize that the soft-sphere model tolerates fewer clashes than the hard-sphere model, because it allows no hard-sphere clashes (at $`\alpha _h`$) and only a limited number of soft-sphere clashes (at $`\alpha _s`$). Using the steric LJ model, the number of repackings is $`10^{25}`$. Since the number of repackings scales exponentially in the number of free torsions, we give the number of rotamer states per torsion in Table 2 for each protein, calculated from the LJ model. These numbers are obtained by taking the $`n`$-th root of the total number of states, where $`n`$ is either the number of free torsions or residues. Normalizing by the number of torsions gives smaller variation over proteins than when normalizing by number of residues, due to compositional differences between proteins (variation in number of torsions per residue among proteins).

As a control, we applied the same sampling method to obtain the distribution of states for all possible sidechain conformations, packed and not packed, as shown in Figure 3. The total number of random sidechain conformations for a protein with $`n`$ free torsions is $`3^n`$. For all three proteins, the area under the random curve obtained from umbrella sampling is $`n\mathrm{log}(3)\pm 2`$, demonstrating that the sampling is reliable, and that the total error is approximately 2 orders of magnitude. For each protein we included the distribution of states using the hard-sphere model (dashed line) in order to show the dramatic reduction in available conformations once excluded volume constraints are applied.

In Figure 4 we show the sampling of packed conformations for subtilisin without umbrella constraints. The same data is shown separately over rms and over torsional distance from the native conformation. We see that the distributions over rms for the different steric models often yield two humps. When plotted over torsional distance, all models give perfect normal distributions. We found no correlation between rms and torsional distance in the range of rms that was sampled without umbrella constraints. One reason for this lack of correlation is that deviations in the $`\chi _1`$ angle yield much larger deviations in rms than the other $`\chi `$ angles. Thus, unless rms is very low, it yields almost no information about the torsional states of the protein.

The packing density of a conformation is usually defined as the volume of the van der Waals envelope of a protein divided by its Voronoi volume. The density of a random sample of 1000 repacked conformations of each protein was computed and found to vary by only $`\pm 1`$% of the native density. Given this minimal variation in packing density, we conclude that the requirement of having native-like density does not significantly reduce the number of packed sidechain conformations. In addition, we were able to build hydrogen atoms (using CHARMM ) on randomly selected, packed conformations, confirming that the conformations were not too compact.

While the number of packed conformations is large, it is conceivable that most repackings are structuraly similar to one another, and thus that the number of truly different repackings may be much smaller than observed. To determine how different these repacked conformations are, we uniquely identified each conformation by the rotamers of the partially occluded positions. The distance between two conformations was measured by counting the number of single-bond torsions which had to be rotated (by $`120^{}`$) in order to obtain one conformation from the other. For example, for the partially occluded positions of subtilisin, having a total of 174 torsions, we expect random conformations to differ by 115 ($`=174\times 2/3`$) torsion moves.

The average torsional distance between repacked conformations (Fig. 2) shows that the conformations are different from each other. The similarity between conformations depends on the particular model. In general, larger soft radii lead to more similar repackings until the soft radii reach the van der Waals size (see Discussion). For the three proteins we examined, average conformations obtained with the LJ model differ at 20%-25% of their torsions. We note that the differences between the soft-sphere models for $`\alpha _s`$$`0.90`$ and the LJ model are minimal in all three proteins.

The entropy of each sidechain was computed by sampling uncorrelated states meeting the LJ constraints for each protein (Fig. 5). The plot shows that entropy is not distributed evenly over positions. Approximately 25% of positions have nearly zero entropy, implying that they are rarely repacked. We note that only 5% of positions have exactly zero entropy, demonstrating that almost all positions can accomodate more than one rotamer. The maximum entropy possible for a sidechain with $`n`$ torsions is $`n\mathrm{ln}(3)=1.1n`$. Since we observe positions with entropy greater than 1.1, it is clear that large residues can take on several different rotamers.

To explain the uneven distribution of entropy, we tried to correlate it with various quantities. The entropy of a given position showed no correlation with either crystallographic B-factors, or with occluded surface area. We found significant but relatively low correlation ($`R=0.65`$) between the entropies we observed and entropies calculated from rotamer probabilities from the PDB database . To investigate this trend further, we plotted the rotamer probabilities observed in our simulation versus probabilities observed in the pdb database. Specifically, for each rotamer state at every partially occluded position we plotted the frequency with which it was observed in our simulation, versus the frequency with which it occurs in the pdb database at its particular $`\varphi `$ and $`\psi `$ backbone angles.

These scatter plots for the three proteins are shown in Figure 6. The plots indicate that native rotamers (red marks, corresponding to the rotamers observed in the crystal structure) are often seen with high probability in our simulations. Since we observe in simulations only a few non-native rotamers (blue marks) with probability greater than 90%, the low entropy positions in Figure 5 largely correspond to the native conformation. The pdb probabilities are obtained by averaging over many different folds and sequences, and thus sequence-specific effects are largely averaged out, leaving only the backbone $`\varphi `$ and $`\psi `$ angles as a determinant of rotamer frequency. Hence the simulation frequency of rotamers lying near the diagonal is explained well by the local backbone conformation at that position. Rotamers lying off the diagonal are being observed with anomalous frequency for their local backbone conformation. Because the backbone conformation does not explain the frequency of these positions, they must be particularly sensitive to one or more sequence-specific interactions. We observe several native rotamers with anomalously high probability (upper-left corner), indicating that sequence-specific steric interactions can completely stabilize the native rotamer even when it is rarely observed in most other proteins. We see a few non-native rotamers anomalously frequently (blue marks in the upper-left of the figure). These rotamers are preferentially selected by sequence-dependent excluded volume interactions over the native rotamers.

One can classify positions, loosely speaking, as either “well-predicted” if the correct native rotamer is seen more than 50% of the time, or “poorly-predicted” if the native rotamer is seen less than 50% of the time. Non-native rotamers are plotted as squares if they belong to “well-predicted” positions, and circles if they belong to “poorly-predicted” positions. A non-native rotamer which is seen with high probability in the PDB database, may be seen with low probability in our simulation, lying in the bottom-right corner of the figures. Most points in this region are squares, that is, non-native rotamers which are being seen too infrequently because they belong to “well-predicted” positions. There are very few circles in this region indicating that, for the most part, significant deviations from observed pdb rotamer probabilities are due to sequence-specific interactions stabilizing the native rotamer preferentially over a more popular rotamer.

Concanavalin appears to have more scatter than the other two proteins. In particular, it has more native rotamers below and near the diagonal. Concanavalin is an all-$`\beta `$ protein, while subtilisin is all-$`\alpha `$, and photoactive yellow protein is $`\alpha `$/$`\beta `$. It is likely that the $`\alpha `$-helices of subtilisin, which act to compactify the fold locally, impose more stringent constraints on the allowed rotamers at a given position. This is manifested in our observation of fewer near-diagonal native rotamers in subtilisin than in concanavalin.

To check whether simulation probabilities correlate with amino acid type, we plotted the sequence identity of the native rotamers as letters in Figure 7. We see that points lying in the lower left-hand corner are generally polar or charged amino acids. Native rotamers seen with high-probability in simulation are mostly hydrophobic. There are exceptions to both of these general trends. Since hydrophobic residues are on average more buried than charged residues, one might suggest that the separation based on polarity is actually due to solvent-exposure. As stated above, however, we found no correlation between occluded surface and the simulation probabilities. In order to verify that this was not due to some peculiarity of the definition of occluded surface, we visually classified positions as buried or partially exposed using rough criteria outlined in Methods. The purple letters in Figure 7 correspond to buried positions (as defined in Methods), while the green letters correspond to partially exposed positions. We see no significant clustering of either color in these figures, confirming that entropy in our simulations is not governed by exposure. We conclude, not surprisingly, that the importance of electrostatic interactions, rather than solvent exposure, explains the poorly-predicted native rotamers of charged amino acids.

We examined the correlation between the rotamer changes made at a pair of positions $`x`$ and $`y`$ by measuring the number of states lost, $`\mathrm{exp}\left(\mathrm{\Delta }S\right)`$, due to interactions between $`x`$ and $`y`$, where $`\mathrm{\Delta }S=S_{xy}(S_x+S_y)`$, $`S_{xy}`$ is the joint entropy, and $`S_x`$ and $`S_y`$ are the individual entropies. We found only 3 pairs of positions in subtilisin exhibiting any correlation ($`\mathrm{exp}\left(\mathrm{\Delta }S\right)`$ between 0.1 and 0.3). Similar numbers were found in the other two proteins. Repacking at a single position can therefore be accomodated by the surrounding residues in many different ways.

Given that there is a significant amount of torsional entropy in our model, it interesting to ask how well an average sidechain structure prediction based solely on excluded volume would perform. For subtilisin, we observe the rms of predictions broadly distributed between 1.4 and 1.8 angsroms. Using either the LJ model or the $`\alpha _s`$$`=0.90`$ model, we predict correct rotamers for 65% of residues, and predict $`\chi _1`$ angles correctly for 85% of positions. The best performing sidechain prediction algorithms, which incorporate several physical and knowledge-based terms, predict 93% of $`\chi _1`$ angles, so excluded volume alone performs very well in this regard.

## Discussion

Upon folding to its native backbone topology, the conformational space of a protein is cut by many orders of magnitude (Fig. 3) due to excluded volume. At the level of a single torsion, the number of rotamer states is reduced by approximately a factor two. The native backbone can still accomodate many different sidechain conformations, while maintaining its high degree of compactness. Other sources of stabilization (such as attractive van der Waals interactions, hydrogen bonding, polar interactions) are necessary to overcome the final entropic cost associated with fixing the sidechains in their native conformations.

The packing of sidechains in a native backbone has often been described as a jigsaw-puzzle: the shapes of the sidechains (the pieces) are enough to determine how they will fit in the protein interior . This argument has been used to explain the high density observed in proteins. The large number of repackings and the diversity of these conformations (Fig. 2), all with native-like density, suggest that a precise fitting of residues is not necessary for a high packing density. The rotamer frequency plots demonstrate that while there is a sequence-specific excluded volume effect that acts to single out the native rotamers, it is not sufficient to overcome the fact that the backbone generally allows several different rotamers at each position. The distribution of entropy shows that many positions admit significantly more than one rotamer. Furthermore, the correlation calculations show that a conformational change at one position does not necessitate a particular change at spatially proximate positions. We conclude that while sidechains must fit together to avoid steric clashes and maintain native-like density, these constraints alone do not give rise to precise and unique fitting of residues. Other energetic interactions must bias sidechain conformations toward the native state, leading to the observation of directional packing in protein structures .

The three models employed in our study yield consistent results. All three vastly reduce the total number of available states. The hard-shere model is consistently more forgiving than the other models. The soft-sphere models with $`\alpha _s`$ $`0.90`$ all roughly agree to $`\pm 2`$ orders of magnitude. The LJ model is only slightly more tolerant than the soft-sphere models, suggesting that a two-shell model can adequately model sterics without introducing continuous distance dependence. The soft-sphere model has two limiting behaviors. For $`\alpha _s`$= $`\alpha _h`$, the soft-sphere model is identical to the hard-sphere model. For $`\alpha _s`$$`1`$, all atom pairs in the *native* structure are considered to clash, and thus the soft-sphere model will tolerate all possible clashes. Since the hard-sphere excluded volume is the only remaining constraint, the soft-sphere model is again equivalent to the hard-sphere model in this limit. The observed rise in the number of states as $`\alpha _s`$ reaches 1.0 is the necessary turn-around between these two limiting behaviors. Sterics must therefore be modelled with $`\alpha _s`$$`<1`$. When $`\alpha _s`$$`=1.0`$, the radii have reached van der Waals sizes, and are in the regime of attraction rather than excluded volume. While increasing the radii past the hard-core value gives a significant gain in stability, we see that there is a limit on how much stabilization can be achieved from excluded volume effects.

Our ability to predict 85% of $`\chi _1`$ angles correctly shows that excluded volume is of primary importance in determining the general orientation of a sidechain with respect to its backbone. Our prediction of only 65% of positions correctly shows that sterics are less important in determining the rest of the sidechain’s rotameric conformation. Of particular interest are native rotamers that are observed in the PDB database in less than 10% of proteins. Many of these rotamers are observed with significantly higher probability in our simulation, while some are observed at their expected PDB level. The excluded volume interaction is thus exhibiting some frustration between the distribution of states compatible with the local backbone conformation, and the sequence-specific neighborhood of the residue. Other interactions must provide significant stabilization to overcome this problem.

A natural extension of our work is to add an attractive interaction to the potential which would stabilize the native packing. From the distribution of packed conformations, we estimate that the native conformation must be about $`0.6kT`$ per sidechain (see Table 2) lower in energy than the average packed conformation with native backbone. This estimate can be used to guide efforts towards designing better potentials for sidechain packing. Furthermore, the techniques presented here, combined with a stepwise increase in the complexity of potentials, can determine how much stabilization each additional interaction provides. It will be interesting to see whether a potential that significantly stabilizes the native rotamers can stabilize the native fold when backbone motions are added to the model.

## Acknowledgements

We are grateful for stimulating discussions with Gabriel Berriz, Alexey Finkelstein, Leonid Mirny, and Michael Morrissey, and especially thank Kirk Doran for his assistance. Financial support from NIH (grant 52126) and NSF (graduate fellowship, E.K.) is acknowledged.

## Figures

Figure 1
The number of soft clashes in near-native conformations as a function of $`\alpha _s`$ for each protein studied. The solid line denotes the average number observed in an ensemble of near-native conformations generated by adding $`\pm 10^{}`$ random noise to the native sidechains. The dotted line is determined by taking three standard deviations below the average, and corresponds to the soft clash threshold used in this study.

Figure 2
The figures in the left column show the density of packed conformations as a function of rms for various values of $`\alpha _s`$(indicated by the colored, solid curves). The dashed curve is the corresponding distribution using the repulsive Lennard-Jones model. Figures in the right column show the average number of bond rotations needed to bring one conformation to another as a function of rms for the various models. The average is taken over all pairs of 1000 randomly selected conformations found in a given rms window.

Figure 3
The distribution of random sidechain conformations for each of the three proteins (solid lines). Umbrella sampling over rms without any excluded volume constraints was used to construct these figures. The dashed lines correspond to the hard-sphere model (identical to the black curves from Figure 2) and are provided for comparison.

Figure 4
Uncorrelated packed conformations for subtilisin sampled using the various models with no rms constraints are histogrammed separately over rms from the native crystal structure, and over torsional distance from the native rotamers. Each histogram is obtained over a total of 240,000 states. The states shown in this plot correpond to the tops of the rms histograms from Figure 2, but were obtained without umbrella sampling. We therefore do not see any low rms states in these plots, because there are exponentially fewer states as rms is lowered. The native rotamer configuration of subtilisin was obtained by minimizing rms in simulation. The rms of this nearly-native rotameric state is 0.5 Å from the crystal structure. The reason for the non-zero rms is that deviations of more than 10 degrees from perfect rotamers are observed in the crystal structure, and our simulations allows for no more than 10 degrees of tolerance around each rotameric value for a given sidechain $`\chi `$ angle.

Figure 5
The entropy for the partially occluded positions of each protein. 50,000 uncorrelated packed conformations were sampled for each protein using the repulsive LJ model without umbrella sampling. The probability of observing each rotamer at each position was calculated, and the entropy of each position was computed using $`S=p_i\mathrm{log}p_i`$, where $`p_i`$ is the frequency of occurence of the $`i`$-th rotamer at the given position. The histograms of these entropy values are also provided for each protein.

Figure 6
Observed simulation frequency of a rotamer vs. the PDB database frequency of that rotamer, given its amino acid identity and the $`\varphi `$-$`\psi `$ angles of its backbone. Database probabilities were taken from the rotamer library described in . Each point in the plots corresponds to a single rotamer. Red points are native rotamers (there is one red point for each sequence position); blue points are non-native rotamers (there are multiple blue points for each position). We called positions “well-predicted” if the native rotamer was observed more than 50% of the time, and “poorly-predicted” otherwise. We further classified the non-native rotamers (blue points) as squares if they belonged to “well-predicted” positions or circles if they belong to “poorly-predicted” positions. For example, if a native rotamer at a given position is observed 85% of the time, all the other, non-native rotamers for that position will be squares. Conversely, if at a different position the native rotamer is seen only 15% of the time, all the non-native rotamers for that position will be circles. Simulation frequencies are calculated over 50,000 packed conformations obtained using the LJ model without umbrella sampling.

Figure 7
Simulation vs. PDB frequency for each *native* rotamer. The native rotamers for each protein are identified by their amino acid type (letter) and are colored purple if they are in a buried position, and green if they are in a partially exposed position. The buried positions are given in the methods section, indicated in bold type.

## Tables

Table 1
The atomic groups, VdW radii ($`r`$), and the hard core distances ($`\alpha _h`$$`r`$) used in our study. The atom groups and VdW radii were obtained from Table 2 of Tsai *et al*. (1999). The atomic group A$`n`$H$`m`$ refers to an atom of element A with a chemical valence of $`n`$ and $`m`$ hydrogen atoms bonded to it. A methyl carbon ($``$CH<sub>3</sub>), for example, falls under the C4H3 atomic group.

## Tables

Table 2
Statistics of repackings for photoactive yellow protein (1F9I), concanavalin (1NLS), and subtilisin (1GCI) using the repulsive Lennard-Jones model. The second column lists the number of torsions that contribute to the overall entropy. These torsions belong to residues having more than 40% occluded surface area. The number of occluded residues having one or more free torsions is given in the third column, with the total length of the protein in parentheses. The number of states and entropy (S) per $`\chi `$ and per residue are calculated over the occluded residues only.


text:
# Radiative four–fermion processes at LEP2

## 1 Introduction

One of the main goals of electroweak physics at LEP2 is the study of the properties of the $`W`$ bosons . The center of mass (c.m.) energy above the threshold of $`W`$-pair production offers the possibility to extract information about the mass and couplings of the $`W`$ boson from the analysis of four-fermion ($`4f`$) final states. As well known, because of the precision of the experimental measurements, radiative corrections to $`e^+e^{}4f`$ processes are needed in order to provide suitable theoretical predictions . From this point of view, events with four fermions plus a visible photon, i.e $`e^+e^{}4f+\gamma `$, are a building block of the full $`O(\alpha )`$ electroweak corrections to $`4f`$ processes, providing the hard bremsstrahlung contribution. Furthermore, radiative $`4f`$ processes are also an interesting physics subject by themselves, since the luminosity achieved at LEP makes them directly accessible to the experimental investigation, as recently discussed in refs. , where first results on the measurement of the $`W^+W^{}\gamma `$ cross section have been reported.

A very peculiar feature of the processes under consideration is that they give the opportunity of directly testing the non-abelian structure of the gauge boson self-interactions. Actually, as other processes studied at LEP, such as $`e^+e^{}W^+W^{}4f`$, $`e^+e^{}We\nu _e`$ (single-$`W`$ production) and $`e^+e^{}\nu \overline{\nu }\gamma `$, $`e^+e^{}4f+\gamma `$ reactions are sensitive to trilinear gauge couplings (TGC). More importantly, they can be used to test quartic gauge couplings (QGC), since they are, together with $`\nu \overline{\nu }\gamma \gamma `$ final states, the only accessible LEP2 processes that contain quartic gauge boson vertices with at least one photon at tree-level . Quartic vertices involving only massive gauge bosons give rise to six-fermion final states and are outside the sensitivity of LEP, being accessible only at the energies of a future $`e^+e^{}`$ Linear Collider (LC). Both charged current (CC) radiative $`4f`$ processes, mediated by two $`W`$-bosons, and neutral current (NC) ones, mediated by two $`Z`$-bosons, are in principle suitable to examine the effect of possibly anomalous gauge couplings (AGC). In this work particular attention is paid to CC processes, because of the larger cross section of $`WW\gamma `$ events with respect to $`ZZ\gamma `$ final states.

The first tree-level calculations of $`e^+e^{}4f+\gamma `$ processes were performed in refs. . In these papers all electroweak contributions as well as fermion mass effects were accounted for by exploiting different approaches to the automatic calculation of the exact tree-level matrix element . Since then, some event generators for the simulation of $`4f+\gamma `$ events in $`e^+e^{}`$ collisions have been developed: RacoonWW , a generator based on the calculation of all $`4f+\gamma `$ final states in the massless approximation ; CompHEP  and grc4f , which are general-purpose packages relying upon the automatic calculation of tree-level amplitudes (including fermion masses) and phase-space; Helac/Phegas , a program implementing a recursive algorithm for the calculation of the scattering amplitudes. The interested reader is referred to ref.  for more details. In addition to the above computational tools, also a calculation of the massive matrix element of $`e^+e^{}4f+\gamma `$ processes have recently appeared in the literature , accompanied by a detailed phenomenological analysis of fermion mass effects in $`4f`$ and $`4f+\gamma `$ final states at LEP2 energies . In all the theoretical studies devoted insofar to $`4f+\gamma `$ production, the effects of quartic anomalous gauge couplings (QAGC), which are a window on the mechanism of spontaneous symmetry breaking  and are presently of special experimental interest, have not been taken into account. Actually, recent phenomenological studies on the subject of QAGC at high-energy $`e^+e^{}`$ colliders have been performed by considering three-vector boson $`WW\gamma ,ZZ\gamma ,Z\gamma \gamma `$ production and treating $`W,Z`$ particles in the on-shell approximation . Anomalous quartic couplings in $`\nu \overline{\nu }\gamma \gamma `$ production via $`WW`$ fusion have been analyzed in ref. . Experimental searches for QAGC at LEP rely upon the theoretical results of refs.  and make use of the computational tools of refs. .

In the light of the present situation and in view of future measurements at LC, a full calculation of $`e^+e^{}4f+\gamma `$ processes, including the effects of AGC and of the most important radiative corrections, is desirable. This task is accomplished in the present paper, by presenting the new event generator WRAP (W Radiative process with ALPHA  and Pavia) for the simulation of $`4f+\gamma `$ processes at $`e^+e^{}`$ colliders. The program is based on the calculation of exact matrix elements, including the effect of fermion masses, both for CC and NC processes. Charged trilinear anomalous gauge couplings (TAGC) and the genuinely QAGC, i.e. those giving no contributions to trilinear vertices, are included in the calculation, as well as the large effect of initial state radiation (ISR). A tuned comparison between the predictions of a preliminary version of WRAP and of the other two event generators RacoonWW and Helac/Phegas can be found in ref. .

The paper is organized as follows. In Sect. 2 the main features of the the calculation are described. After the description of the treatment of the multi-particle phase-space, the theoretical details concerning the calculation of the exact matrix elements, the implementation of anomalous gauge couplings and ISR are given. A sample of numerical results as obtained by means of the Monte Carlo WRAP is presented in Sect. 3, paying particular attention to the contribution of fermion masses, to the impact of ISR and to the effects of AGC at LEP2 and LC energies. Conclusions and possible perspectives are drawn in Sect. 4.

## 2 Features of the calculation

### 2.1 Phase-space integration

The kinematics of the $`25`$ particles processes has been treated generating the $`5`$-body phase-space recursively, since the process can be seen as production and subsequent decay of a pair of massive gauge bosons.

Concerning CC like processes, the configurations of interest are related to a photon emitted from the initial state (see Fig. 1), from an intermediate $`W`$-boson (see Fig. 2) and from the final state charged fermions (see Fig. 3). As far as emission from the initial state is concerned, the adopted phase-space decomposition reads as follows:

$$\begin{array}{ccc}d\mathrm{\Phi }_5=\hfill & (2\pi )^6\hfill & d\mathrm{\Phi }_3(P;p_5,Q_{V_1},Q_{V_2})d\mathrm{\Phi }_2(Q_{V_1};p_1,p_2)\hfill \\ & & d\mathrm{\Phi }_2(Q_{V_2};p_3,p_4)dQ_{V_1}^2dQ_{V_2}^2,\hfill \end{array}$$
(1)
where $`V_1`$ and $`V_2`$ indicate the $`W`$ gauge bosons, the momenta $`p_i`$ with $`i=1,\mathrm{},4`$ stand for the momenta of the final state fermions and $`p_5`$ is the photon momentum. The eleven independent variables have been chosen to be:

* photon variables $`E_\gamma `$, $`\theta _\gamma `$ and $`\varphi _\gamma `$ in the c.m. frame;
* invariant mass squared $`Q_{V_1}^2`$ and $`Q_{V_2}^2`$;
* three $`\theta `$ and $`\varphi `$ angle pairs in the rest frame of each decaying “particle”, namely in the frames given by the conditions $`𝐏𝐩_5=\mathrm{𝟎}`$, $`𝐐_{V_1}=\mathrm{𝟎}`$ and $`𝐐_{V_2}=\mathrm{𝟎}`$.

When the photon is emitted from the final state, which means for a CC process from a virtual $`W`$ or from a virtual charged fermion, the following decomposition can be conveniently considered:

$$\begin{array}{ccc}d\mathrm{\Phi }_5=\hfill & (2\pi )^6\hfill & d\mathrm{\Phi }_2(P;Q_{V_1},Q_{V_2})d\mathrm{\Phi }_3(Q_{V_1};p_5,p_1,p_2)\hfill \\ & & d\mathrm{\Phi }_2(Q_{V_2};p_3,p_4)dQ_{V_1}^2dQ_{V_2}^2.\hfill \end{array}$$
(2)
In the case of the photon emitted from an internal gauge boson the independent variables can be chosen as follows:

* photon variables $`E_\gamma `$, $`\theta _\gamma `$ and $`\varphi _\gamma `$ in the c.m. frame;
* invariant mass squared $`Q_{V_1}^2`$ and $`Q_{V_2}^2`$;
* one $`\theta `$ and $`\varphi `$ $`W`$-angle pair in the c.m. frame;
* two $`\theta `$ and $`\varphi `$ angle pair for $`𝐩_1`$ and $`𝐩_3`$ in the rest frame of the bosons $`V_1`$ and $`V_2`$, respectively.

In the case of photon emission from a final state fermion the following independent variables can been adopted:

* invariant mass squared $`Q_{V_1}^2`$ and $`Q_{V_2}^2`$;
* one $`\theta `$ and $`\varphi `$ $`W`$-angle pair in the c.m. frame;
* one $`\theta `$ and $`\varphi `$ angle pair for $`𝐩_3`$ in the rest frame of the boson $`V_2`$;
* energies of $`p_1`$ and $`p_5`$ momenta in the c.m. frame;
* azimuthal angle $`\varphi `$ of $`p_1`$ in the rest frame of $`V_1`$;
* $`\varphi _\gamma `$ in the rest frame of the radiating fermion;
* $`\mathrm{cos}\theta _{\gamma f}`$ in the c.m. frame ,

where $`\theta _{\gamma f}`$ is the relative angle between the radiating fermion and the photon.

An analogous phase-space parameterization has been implemented in WRAP for the case of NC processes, neglecting of course the channels related to the photonic emission from internal lines. The above phase-pace decompositions, iterated for each possible radiation pattern, give rise to several channels, depending on the final state considered.

In the previous formulas $`d\mathrm{\Phi }_n`$ represents the element of $`n`$-body phase-space given by:

$$d\mathrm{\Phi }_n(P;p_1,\mathrm{},p_n)=(2\pi )^4\delta ^4(P\underset{i=1}{\overset{n}{}}p_i)\underset{i=1}{\overset{n}{}}\frac{d^3p_i}{(2\pi )^32E_i}.$$
(3)
The code works taking into account all the configurations discussed above according to a standard multi-channel Monte Carlo approach .

In order to perform an efficient event generation, the peaking behaviour of the matrix element has been treated in the following way:

* the squared invariant masses of the massive gauge bosons $`V`$ are sampled according to a Breit-Wigner distribution centered around $`M_V^2`$, while the photon propagator is sampled according to the $`1/Q^2`$ distribution;
* the infrared divergence is sampled according to a $`1/E_\gamma `$ distribution;
* the collinear peak arising from the photon emission due to an external charged fermion is taken under control by sampling it with a distribution proportional to $`1/(1\beta \mathrm{cos}\theta )`$, where $`\theta `$ is the separation angle between the radiating fermion and the photon, and $`\beta =\sqrt{1m^2/E^2}`$, $`m`$ and $`E`$ being the mass and the energy of the fermion respectively.

For a realistic account of gauge boson properties, and to avoid integration singularities, it is mandatory to include the gauge boson width in the propagators. The so-called fixed-width scheme  is adopted in WRAP. Actually, as shown in ref. , the fixed-width scheme, even if it violates $`SU(2)`$ gauge invariance, is a reliable $`U(1)`$ gauge-restoring method and is able to guarantee predictions for $`e^+e^{}4f+\gamma `$ processes in good numerical agreement with a scheme preserving all the relevant Ward identities, such as the complex-mass scheme .

### 2.2 Tree-level matrix element

As already mentioned, the present work is based on the calculation of the fully massive Born matrix element of $`e^+e^{}4f+\gamma `$ processes. The exact matrix elements for CC and NC $`e^+e^{}4f+\gamma `$ processes are available in WRAP. The calculation is performed by using ALPHA , an iterative algorithm for the automatic evaluation of tree-level scattering amplitudes without using Feynman graphs (see ref.  for a review of the method and of recent phenomenological applications). For the processes under consideration, a completely numerical approach turns out to be particularly convenient not only for the very large amount of contributing Feynman diagrams, but also because the calculation can be performed in the presence of fermion masses without any additional complication. This is of special importance for $`4f+\gamma `$ final states involving muons, where the separation angle between muon and photon can be realistically set to zero, and a calculation taking into account of the finite muon mass is mandatory, to avoid collinear singularities.

### 2.3 Anomalous gauge couplings

Information about the structure of TGC and QGC can be obtained by the analysis of $`4f+\gamma `$ production processes. In particular, CC radiative $`4f`$ processes, although characterized by a lower statistics, are potentially a complementary channel to the $`4f`$ final states in order to test the effect of TGC, because of the larger amount of diagrams involving trilinear gauge interactions. More importantly, $`4f+\gamma `$ processes are interesting in order to put bounds on deviations from standard quartic gauge couplings. In the following, the theoretical details of the parameterization adopted in order to keep under under control this important phenomenological issue are described.

#### 2.3.1 Trilinear anomalous gauge couplings

It is possible to take into account the effect of charged TGC (anomalous and not) by means of the following lagrangian :

$$\begin{array}{cc}iL_{TGC}\hfill & =g_{WWV}[g_1^VV^\mu (W_{\mu \nu }^{}W^{+\nu }W_{\mu \nu }^+W^\nu )\hfill \\ & +\kappa _VW_\mu ^+W_\nu ^{}V^{\mu \nu }+\frac{\lambda _V}{m_W^2}V^{\mu \nu }W_\mu ^{+\rho }W_{\rho \nu }^{}\hfill \\ & +ig_5^V\epsilon _{\mu \nu \rho \sigma }(\left(^\rho W^\mu \right)W^{+\nu }W^\mu \left(^\rho W^{+\nu }\right))V^\sigma ]\hfill \\ & V=\gamma ,Z,\hfill \end{array}$$
(4)
which represents the most general lagrangian describing trilinear $`WWV`$ gauge interactions, with the exception of the operators violating $`C`$, $`P`$ and $`CP`$ simmetries. $`L_{TGC}`$ has been implemented in ALPHA, and the presence of anomalous couplings can be studied, as done at LEP , by using the relations

$$\mathrm{\Delta }\kappa _\gamma =\frac{c_W^2}{s_W^2}(\mathrm{\Delta }\kappa _Z\mathrm{\Delta }g_1^Z)\lambda _Z=\lambda _\gamma \lambda ,$$
(5)
where $`\mathrm{\Delta }\kappa _V=\kappa _V1`$ and $`\mathrm{\Delta }g_1^Z=g_1^Z1`$. The Standard Model (SM) Lagrangian is recovered for $`g_1^V=k_V=1`$, $`\lambda =0`$, $`g_5^V=0`$. Triple anomalous neutral gauge couplings, considered in ref.  and looked for at LEP in $`e^+e^{}Z\gamma ,ZZ`$ processes , are not presently taken into account.

#### 2.3.2 Quartic anomalous gauge couplings

Quartic gauge couplings involving at least one photon are analyzed at LEP . In particular, $`W^+W^{}\gamma \gamma `$ and $`W^+W^{}Z\gamma `$ vertices are probed in $`WW\gamma 4f+\gamma `$ and $`\nu \overline{\nu }\gamma \gamma `$ final states , while $`e^+e^{}Z\gamma \gamma `$ processes  are investigated to put bounds on the $`ZZ\gamma \gamma `$ vertex, which is a gauge interaction not predicted by SM at tree level. In the present work the operators considered in ref.  for genuine anomalous quartic couplings containing at least one photon, namely $`W^+W^{}Z\gamma `$, $`W^+W^{}\gamma \gamma `$ and $`ZZ\gamma \gamma `$ vertices, have been implemented in ALPHA, upgrading the version used in ref.  for the analysis of QAGC in six-fermion final states at the energies of future linear colliders. The implemented lagrangian include all the relevant six-dimensional operators and reads as follows

$`L_{QGC}=`$ $`W_1+W_2+Z_1+Z_2`$ (6)
$`W_0^Z+W_c^Z+W_1^Z+W_2^Z+W_3^Z.`$
In the above equation the Lorentz structure of the operators is given by:

$$\begin{array}{cccccc}W_1=\hfill & a_{w1}F_{\mu \nu }F^{\mu \nu }W_\rho ^+W^\rho \hfill & & & & \\ W_2=\hfill & a_{w2}F_{\mu \nu }F^{\mu \rho }W^{+\nu }W_\rho ^{}+\mathrm{h}.\mathrm{c}.\hfill & & & & \\ Z_1=\hfill & a_{z1}F_{\mu \nu }F^{\mu \nu }Z_\rho Z^\rho \hfill & & & & \\ Z_2=\hfill & a_{z2}F_{\mu \nu }F^{\mu \rho }Z^\nu Z_\rho \hfill & & & & \\ W_0^Z=\hfill & a_{wz0}F_{\mu \nu }Z^{\mu \nu }W_\rho ^+W^\rho \hfill & & & & \\ W_c^Z=\hfill & a_{wzc}F_{\mu \nu }Z^{\mu \rho }W^{+\nu }W_\rho ^{}+\mathrm{h}.\mathrm{c}.\hfill & & & & \\ W_1^Z=\hfill & a_{wz1}F_{\mu \nu }W^{+\mu \nu }Z^\rho W_\rho ^{}+\mathrm{h}.\mathrm{c}.\hfill & & & & \\ W_2^Z=\hfill & a_{wz2}F_{\mu \nu }W^{+\mu \rho }Z^\nu W_\rho ^{}+\mathrm{h}.\mathrm{c}.\hfill & & & & \\ W_3^Z=\hfill & a_{wz3}F_{\mu \nu }W^{+\mu \rho }Z_\rho W^\nu +\mathrm{h}.\mathrm{c}.,\hfill & & & & \end{array}$$
(7)
where the $`a_i`$ are coefficients of dimension $`M^2`$. It is worth noticing that, by imposing appropriate relations between the $`a_i`$’s, symmetry properties, such as for instance $`SU(2)_c`$ custodial symmetry or $`SU(2)\times U(1)`$ gauge invariance, can be guaranteed, as shown in ref. . In the parameterization adopted in ref.  the $`a_i`$ are real coefficients whose explicit expression can be directly read off from the corresponding operator structure of ref.  itself. In particular, the coefficients $`a_0`$ and $`a_c`$, originally introduced in ref.  and related to the $`WW\gamma \gamma `$ and $`ZZ\gamma \gamma `$ structure, can be obtained from the above $`a_i`$ coefficients by means of the following relations

$$\begin{array}{ccc}a_{w1}=\hfill & \frac{e^2}{8\mathrm{\Lambda }^2}& a_0\hfill \\ a_{z1}=\hfill & \frac{e^2}{16\mathrm{cos}^2\theta _w\mathrm{\Lambda }^2}& a_0\hfill \\ a_{w2}=\hfill & \frac{e^2}{16\mathrm{\Lambda }^2}& a_c\hfill \\ a_{z2}=\hfill & \frac{e^2}{16\mathrm{cos}^2\theta _w\mathrm{\Lambda }^2}& a_c,\hfill \end{array}$$
(8)
where $`\mathrm{\Lambda }`$ represents a scale of new physics. As far as $`WWZ\gamma `$ vertex is concerned, an additional structure has been proposed in the literature , whose expression can be derived from the above $`a_i`$ coefficients by means of the following relations

$$\begin{array}{c}a_{wzc}=i\frac{e^2}{16\mathrm{cos}\theta _w\mathrm{\Lambda }^2}a_n\hfill \\ a_{wz2}=i\frac{e^2}{16\mathrm{cos}\theta _w\mathrm{\Lambda }^2}a_n\hfill \\ a_{wz3}=i\frac{e^2}{16\mathrm{cos}\theta _w\mathrm{\Lambda }^2}a_n.\hfill \end{array}$$
(9)
On the experimental side, bounds on $`a_0`$, $`a_c`$ and $`a_n`$ couplings are quoted by LEP collaborations via the analysis of $`WW\gamma `$ and $`\nu \overline{\nu }\gamma \gamma `$ final states . It is worth noticing, in passing, that, thanks to the implementation of the lagrangian of eq. (6) in the ALPHA code, an improved version of the Monte Carlo generator NUNUGPV  is also available for the study of QAGC in $`\nu \overline{\nu }\gamma \gamma `$ events.

The result of the implementation in ALPHA has been carefully cross-checked by an independent analytical calculation of all $`V_1V_2V_3V_4`$ amplitudes, with $`V_i=\gamma ,W,Z`$. The check has been performed for all the processes obtained from $`W^+W^{}Z\gamma `$ scattering, by permutating particles.

### 2.4 Initial-state radiation

In order to match the precision of LEP measurements, the most important radiative corrections have to be considered. Among them, it is well known that undetected initial-state radiation (ISR) plays a major r$`\widehat{\mathrm{o}}`$le. It can be taken into account in the leading log approximation by using the QED Structure Function (SF) approach, in terms of collinear  or $`p_{}`$-dependent SF . Following recent work done in ref. , both prescriptions are available in WRAP, for the reason explained below. When ISR is included via collinear SF, the QED corrected cross section can be written as

$$\sigma _{QED}^{4f+1\gamma }(s)=𝑑x_1𝑑x_2D(x_1,s)D(x_2,s)𝑑\sigma _0^{4f+1\gamma }(x_1x_2s),$$
(10)
by convoluting the tree-level cross section with electron SF. However, due to the presence of an observed photon in the hard-scattering matrix element, the inclusion of ISR needs some care. Actually, since the collinear SF can be viewed as the result of an integration over the angular variables of the photon radiation, an overlapping between the detected kernel photon and pre-emission photons at large angle may occur<sup>1</sup><sup>1</sup>1The same problem is discussed in detail in ref.  for the process $`e^+e^{}\nu \overline{\nu }+n\gamma `$. We refer the reader to ref.  for more details on the strategy here adopted.. The consequence is that a double counting takes place if higher-order QED corrections are naively included by using collinear SF .

On the other hand, it is expected that the bulk of the correction is well estimated by collinear SF, since the emission of a photon from an on-shell initial state fermion is almost collinear. However, in order to provide a more appropriate treatment of photon corrections and give an estimate of the double-counting effect, the SF method can be improved by means of the use of $`p_{}`$-dependent SF, i.e. by generating angular variables for the ISR photons according to $`1/(pk)`$, which is the leading behaviour for radiation of momentum $`k`$ emitted by an on-shell fermion of momentum $`p`$. In such a scheme, the QED corrected cross section can be calculated as

$`\sigma _{QED}^{4f+1\gamma }=`$ $`{\displaystyle 𝑑x_1𝑑x_2_{\mathrm{\Omega }_c}𝑑c_\gamma ^{(1)}𝑑c_\gamma ^{(2)}}`$
$`\stackrel{~}{D}(x_1,c_\gamma ^{(1)};s)\stackrel{~}{D}(x_2,c_\gamma ^{(2)};s)d\sigma ^{4f+1\gamma },`$ (11)
where $`D(x,c_\gamma ;s)`$ is the $`p_{}`$-dependent SF . According to eq. (11), an “equivalent ” photon is generated and accepted as an ISR contribution only if it satisfies a rejection algorithm based on the following requirements:

* the energy of the “equivalent” photon is below the energy threshold for the observed photon, for arbitrary angles; or
* the “equivalent” photon is collinear to a charged particle (i.e. under the minimum separation angle required in order to be detected), for arbitrary energies.

Within the angular acceptance of the observed photon, the cross section is computed by means of the exact matrix element for the considered $`4f+\gamma `$ final state. Therefore, eq. (11) applies to the signature of four fermions plus an isolated hard photon, corrected by the effect of undetected soft and/or collinear radiation. Aiming to obtain a correct evaluation of the size of the double counting effects, a limit of the present treatment of undetected radiation is that only ISR is actually considered. This issue could be addressed in a more complete way by using, for example, a QED parton shower approach as developed in ref. , in order to describe the radiation from all external charged legs, thus including the contribution of undetected final-state radiation.

## 3 Numerical results

The aim of the present section is to give some details on the technical precision of WRAP and discuss the impact of the effects due to fermion masses, ISR and AGC on observables of experimental interest.

In order to test the reliability and the theoretical accuracy of the event generators, a detailed tuned comparison between the predictions of WRAP and other available programs have been carried out in the context of the four-fermion working group of the LEP2 Monte Carlo workshop at CERN . The comparisons, referred to integrated cross sections and differential distributions of several CC processes, showed perfect technical agreement. The comparison is here extended to NC processes, as shown in Tab. 1, between the predictions of WRAP and RacoonWW with input parameters and cuts as in ref. . As can be seen, also for NC final states the agreement is excellent. A further comparison between the predictions of WRAP and those of ref.  is reported in Tab. 2, for several cross sections of CC and NC processes, in the presence of finite fermion masses and in terms of the same input parameters and cuts as adopted in ref. . Perfect agreement is registered for all the considered $`4f+\gamma `$ final states.

The phenomelogical analysis makes use of the following input parameters:

$$\begin{array}{cccc}G_F=1.1663710^5\mathrm{GeV}^2\hfill & M_Z=91.1867\mathrm{GeV}\hfill & & \\ M_W=80.35\mathrm{GeV}\hfill & \mathrm{sin}^2\theta _w=1M_W^2/M_Z^2\hfill & & \\ \mathrm{\Gamma }_Z=2.49471\mathrm{GeV}\hfill & \mathrm{\Gamma }_W=2.04277\mathrm{GeV}\hfill & & \\ m_\mu =0.10565839\mathrm{GeV}\hfill & m_s=0.15\mathrm{GeV}\hfill & & \\ m_c=1.55\mathrm{GeV}\hfill & & & \end{array}$$
(12)
The form used for the propagator of the massive gauge bosons is, according to the fixed-width scheme, $`1/(p^2M^2+i\mathrm{\Gamma }M)`$. The processes considered are the radiative semi-leptonic final states of the kind $`e^+e^{}l^+\nu _lq\overline{q}^{^{}}\gamma `$. The cuts adopted are:

$$\begin{array}{c}|\mathrm{cos}\theta _\gamma |0.985;E_\gamma 1\mathrm{GeV};\hfill \\ |\mathrm{cos}\theta _l|0.985;E_l5\mathrm{GeV};\hfill \\ \theta _{\gamma f}5^{};\hfill \\ M_{q\overline{q}^{^{}}}10\mathrm{GeV},\hfill \end{array}$$
(13)
where $`\theta _\gamma ({}_{l}{}^{})`$ is the photon(lepton) scattering angle, $`E_\gamma ({}_{l}{}^{})`$ is the photon(lepton) energy, $`\theta _{\gamma f}`$ is the angular separation between photon and final charged fermions, and $`M_{q\overline{q}^{^{}}}`$ is the $`q\overline{q}^{^{}}`$ invariant mass.

Table 3 shows the effect of fermion masses on integrated cross sections at $`\sqrt{s}=200`$ GeV for two different photon-fermion separation angles. In the first row, the angular separation $`\theta _{\gamma f}`$ between photon and all charged final state fermions is fixed at $`5^{}`$, while in the second row $`\theta _{\gamma f}=1^{}`$. As expected, the relative difference between the massless and massive cross section increases, going from $`2\%`$ of the first row to the $`7\%`$ of the second row, because of the importance of fermion mass contributions when the photon approaches the collinear region around an on–shell charged particle. In the case of a final state containing a muon, the separation angle $`\vartheta _{\gamma f}`$ can be realistically set to zero, because of different behaviour of photons and muons in the experimental apparatus. Table 4 shows the difference between massive and massless cross section, with the minimal separation between quarks and photon fixed at $`5^{}`$ and progressively relaxing the separation cut between muon and photon. It can be seen that the massless calculation is still reliable for $`1^{}`$ of minimum separation, being the relative difference around $`2\%`$, but it becomes inadequate when the separation falls at some fraction of degree, the relative difference being of the order of $`10\%`$. Therefore, in particularly stringent experimental conditions, only a massive $`4f+\gamma `$ calculation can provide reliable predictions in the presence of muons in the final state.

Figure 4 shows the line-shape of the cross sections of the radiative semi-leptonic processes $`\mu ^+\nu _\mu \overline{u}d+\gamma `$ and $`e^+\nu _e\overline{u}d+\gamma `$, as a function of the c.m. energy in the LEP2 range. The QED corrected cross section via collinear SF for the $`\mu ^+\nu _\mu \overline{u}d+\gamma `$ final state is also plotted. The comparison shows that the contribution due to the additional $`t`$-channel diagrams present in the $`e^+\nu _e\overline{u}d+\gamma `$ final state is not particularly relevant for the adopted selection criteria, small being the differences between the cross sections of the two processes. Concerning ISR in the strictly collinear approximation, its impact on the cross section is at the level of $`1015\%`$, which is a phenomenologically relevant effect in the light of the LEP experimental accuracy. It is worth noticing that this result, obtained by means of a standard treatment of ISR as tipically adopted in the experimental analysis of QAGC in radiative events at LEP , just provides the bulk of the effect due to ISR but it is affected, as previously discussed and quantified below, by a double counting because of the presence of a radiative process as hard-scattering reaction.

The contribution of initial-state photon radiation is also shown in Figs. 5, as a function of the threshold energy $`E_\gamma ^{min}`$ of the observed photon. It can be noticed that the reduction factor due to collinear ISR is around 12-13%, almost independent of the photon detection threshold. However, as previously discussed, ISR in the collinear SF scheme introduces a double counting effect when the pre-emission “equivalent” photon enters the phase-space region of the kernel photon. In order to get an estimate of this overlapping contribution, the comparison of the corrections due to the collinear SF and $`p_{}`$-dependent SF is shown. Is is observed that the two prescriptions can differ at 5% level for $`E_\gamma ^{min}`$ close to 1-2 GeV, while the becomes smaller and smaller as $`E_\gamma ^{min}`$ increases. Numerical investigation points out that, as expected, the discrepancy between collinear and $`p_{}`$-dependent SF is larger near the soft and collinear region and at the level of some per cent, thus yielding an estimate of the size of the double-counting effect at the level of ISR. Therefore, in the presence of particularly stringent experimental constraints sensitive to the soft and collinear emission, precise predictions demand a treatment of ISR able to keep under control the transverse degrees of freedom of photon radiation.

Let us come to the discussion of the effects due to (a sample of) AGC. Both integrated cross sections (Figs. 6-9) and differential distributions (Figs. 10-13) are considered. In Figs. 6-7 the (relative) effect of the TAGC $`\lambda `$ on the $`e^+e^{}\mu ^+\nu _\mu \overline{u}d\gamma `$ cross section is examined, by plotting the relative difference between the cross section in the presence of a non-vanishing $`\lambda `$ coupling and the SM cross section ($`\lambda =0`$ ), as a function of the $`\lambda `$ value at $`\sqrt{s}=192`$ GeV. Fig. 6 shows a comparison of the effect of the $`\lambda `$ coupling on the radiative $`\mu ^+\nu _\mu \overline{u}d+\gamma `$ process and the corresponding $`4f`$ final state, which, as already remarked, differ for their content of trilinear gauge interactions. The numerical results for the $`4f`$ process have been obtained by means of the program WWGENPV . For the considered $`\lambda `$ values, the relative contribution is almost the same on the two processes (obviously the cross sections are quite different), giving a difference at $`23\%`$ level only for extreme $`\lambda `$ values. Therefore, in the presence of standard cuts on the observed photon, trilinear gauge interactions due to $`W`$ radiation in radiative $`4f`$ processes doesn’t enhance the sensitivity to TAGC with respect to a pure $`4f`$ final state. This conclusion is further corroborated by the results shown in Fig. 7, where the effect of the $`\lambda `$ coupling is studied for different photon cuts, with the aim of suppressing the mostly collinear fermion radiation by imposing more and more severe cuts on the detected photon. By comparing the relative deviations shown in Fig. 6 and Fig. 7, one can conclude that in radiative $`4f`$ processes $`W`$ radiation can be hardly disentangled from the radiation off fermions, being the observed deviations almost at the same level for all the set of cuts considered.

The effect of the QAGC $`k_0^W`$, as defined in ref. , is shown in Figs. 8-9 for the process $`\mu ^+\nu _\mu \overline{u}d\gamma `$, as a function of the parameter $`k_0^W`$ at $`\sqrt{s}=200`$ GeV. For the scale of new physics $`\mathrm{\Lambda }`$, the value $`\mathrm{\Lambda }=M_W`$ is used, as conventionally done in the literature. Absolute and relative effects are shown in Fig. 8 and Fig. 9, respectively. In terms of the coefficients $`a_i`$ of eq. (7) the $`k_0^W`$ coupling can be expressed as:

$$\begin{array}{cc}a_{w1}\hfill & =\frac{e^2g^2}{2\mathrm{\Lambda }^2}k_0^W\hfill \\ a_{wz0}\hfill & =\frac{e^2g^2}{\mathrm{\Lambda }^2}\frac{\mathrm{cos}\theta _w}{\mathrm{sin}\theta _w}k_0^W.\hfill \end{array}$$

The solid line refers to the complete $`4f+\gamma `$ calculation of WRAP with input parameters and photon cuts as used in ref.  and the additional cuts on fermions as given by eq. (13). In order to compare with the results of ref. , the dash-dotted line has been obtained by means of a calculation of the process $`e^+e^{}W^+W^{}\gamma `$ (performed independently and in agreement with the results of ref. ), by taking into account the suitable branching ratios of the $`W`$ bosons decaying into $`\mu ^+\nu _\mu `$ and $`\overline{u}d`$ pairs. The dotted line is the prediction as obtained by WRAP, with additional cuts on the invariant masses of the two fermionic pairs constrained within 75 GeV and 85 GeV, in order to enhance, as much as possible, the contribution of diagrams with two resonant $`W`$ bosons. It can be clearly noticed that, even in the presence of cuts on the invariant masses of the decay products, the complete $`4f+\gamma `$ calculation differ from the prediction of the $`WW\gamma `$ approximation, thus proving the importance of a full calculation for the extraction of meaningful limits on QAGC.

In Figs. 10-13 the most important photonic distributions are displayed using the code as event generator with the cuts of eq. (13) at a typical LEP2 energy $`\sqrt{s}=192`$ GeV. In each plot, the SM Born and the QED corrected predictions are compared with those obtained in the presence of AGC. The values used for the anomalous couplings are: $`\lambda =0.25`$ and $`k_0^W/\mathrm{\Lambda }^2=0.01`$. For the sake of comparison, all the data sample are normalized to the same luminosity. The $`\mathrm{cos}\theta _\gamma `$ distribution and the distribution of the cosine of the angle between the photon and the nearest charged particle are shown in Fig. 10 and Fig. 11, respectively. Typical peaking behaviour in the close-to-collinear regions is clearly registered. In such regions, a particularly significative impact of the QACG $`k_0^W`$ is also observed. Fig. 12 and Fig. 13 refer to the $`E_\gamma `$ and transverse photon momentum $`p_{}`$ distribution, respectively, showing the characteristic infrared peak. As already noticed in ref. , these observables turn out to be particularly sensitive to the presence of a QAGC in the region of high energy and $`p_{}`$, being the operator involved of derivative type with respect to the photon field. In all the considered distributions, ISR introduce sizeable effects if compared with the deviations due to anomalous couplings.

As far as the parameterization of QAGC in terms of $`a_0,a_c,a_n`$ parameters is concerned, numerical results are shown in Tab. 5 and Tab. 6, at $`\sqrt{s}=200`$ GeV and $`\sqrt{s}=500`$ GeV, respectively. The cross sections in the presence of non-vanishing anomalous couplings are compared with the pure SM predictions. By looking at Tab. 5 and Tab. 6, it can be noticed that the sensitivity of the $`4f`$+$`\gamma `$ processes to QAGC is much higher at the energies of a future LC than at LEP2, as a priori expected and already noticed in the literature for the $`WW\gamma `$ process .

## 4 Conclusions

The production of four fermions plus an additional detected photon in $`e^+e^{}`$ collisions is studied at LEP to test electroweak gauge boson couplings and in particular to derive bounds on QAGC. In order to provide predictions of phenomenological interest, an exact calculation of $`4f+\gamma `$ processes, including the effect of fermion masses, AGC and ISR has been performed. On the basis of the experimental accuracy, the contribution of fermion masses and ISR has been analyzed in comparison with typical deviations introduced by AGC. A new Monte-Carlo event generator (WRAP) has been developed and is available for the simulation of radiative $`4f`$ events.

The main conclusions of the present study can be summarized as follows. The effect of finite fermion masses, as analyzed in the $`\mu ^+\nu _\mu \overline{c}s\gamma `$ final state, turn out to be very sensitive to the separation angle $`\vartheta _{\gamma f}`$ between photon and charged fermions, ranging from about 2% for $`\vartheta _{\gamma f}=5^{}`$ to about 7% for $`\vartheta _{\gamma f}=1^{}`$. For the realistic situation of a vanishing separation angle $`\vartheta _{\gamma \mu }`$ a massive calculation is strictly unavoidable.

Particular care has been devoted to the inclusion of ISR, as a consequence of the presence of an observed photon in the final state. The contribution of ISR has been studied in terms of collinear and $`p_{}`$ dependent SF. Numerical results illustrate that ISR introduces corrections of the order of $`1015\%`$ on the integrated cross section. However, in order to get a reliable estimate of ISR corrections and to avoid double counting, $`p_{}`$ photon effects have to be considered. It has been shown that the double counting, affecting the QED corrected cross section via collinear SF, may reach the $`5\%`$ level in a realistic event selection and hence it has to be taken into account carefully. A more accurate evaluation of double counting effects should however consider also the photonic radiation off final state charged fermions.

Both trilinear and genuinely quartic anomalous gauge couplings have been implemented in WRAP, and their effects on total cross section as well as on photon distributions have been investigated. The impact of TAGC on the considered observables does not seem to be very sensitive to the cuts imposed on the detected photon, suggesting that $`W`$ radiation is not easily disentangled from the fermion radiation. Thus $`4f+\gamma `$ final states are not the ideal place where to look for TAGC, if compared with $`4f`$ final states, which benefit of a higher statistics. On the contrary, these radiative processes are significantly affected by QAGC. In particular it has been shown that difference in the effect of QAGC are present between the predictions of the complete calculations by means of WRAP and the ones obtained in the limit of on-shell $`W`$ bosons, which is the approximation presently used in the literature. A more complete investigation of QAGC in radiative events at $`e^+e^{}`$ colliders is currently in progress.

Acknowledgements

The authors wish to thank F. Cavallari, D.G. Charlton, A. Denner, S. Dittmaier, M. Musy, C.G. Papadopoulos, M. Roth, D. Wackeroth, and the participants in the CERN LEP2 Monte Carlo Workshop for useful discussions. A particular thank goes to the authors of RacoonWW for the very helpful collaboration in numerical comparisons and to P. Bell, D.G. Charlton and M. Thomson for useful information and help concerning QAGC.