text:
# 1 DY dilepton production in 𝑝𝑝 scattering

## 1 DY dilepton production in $`pp`$ scattering

Although cross sections are Lorentz invariant, the partonic interpretation of the microscopic process depends on the reference frame. As pointed out in , in the target rest frame, Drell-Yan (DY) dilepton production should be treated as bremsstrahlung, rather than parton annihilation (see also ). The space-time picture of the DY process in the target rest frame is illustrated in Fig. 1. A quark (or an antiquark) from the projectile hadron radiates a virtual photon on impact on the target.

A salient feature of the rest frame picture of DY dilepton production is that at high energies and in impact parameter space the DY cross section can be formulated in terms of the same dipole cross section as low-$`x_{Bj}`$ DIS. The cross section for radiation of a virtual photon from a quark after scattering on a proton, can be written in factorized light-cone form ,

$$\frac{d\sigma (qp\gamma ^{}X)}{d\mathrm{ln}\alpha }=\underset{T,L}{}d^2\rho |\mathrm{\Psi }_{\gamma ^{}q}^{T,L}(\alpha ,\rho )|^2\sigma _{q\overline{q}}(\alpha \rho ),$$
(1)
similar to the case of DIS. Here, $`\sigma _{q\overline{q}}`$ is the cross section  for scattering a $`q\overline{q}`$-dipole off a proton which depends on the $`q\overline{q}`$ separation $`\alpha \rho `$, where $`\rho `$ is the photon-quark transverse separation and $`\alpha `$ is the fraction of the light-cone momentum of the initial quark taken away by the photon. For shortness, we do not explicitly write out the energy dependence of $`\sigma _{q\overline{q}}`$. We use the standard notation for the kinematical variables, $`x_1x_2=x_F`$, $`\tau =M^2/s=x_1x_2`$, where $`x_F`$ is the Feynman variable, $`s`$ is the center of mass energy squared of the colliding protons and $`M`$ is the dilepton mass. In (1) $`T`$ stands for transverse and $`L`$ for longitudinal photons.

The physical interpretation of (1) is similar to the DIS case. The projectile quark is expanded in the interaction eigenstates. We keep here only the first eigenstate,

$$|q=\sqrt{Z_2}|q_{bare}+\mathrm{\Psi }_{\gamma ^{}q}^{T,L}|q\gamma ^{}+\mathrm{},$$
(2)
where $`Z_2`$ is the wavefunction renormalization constant for fermions. In order to produce a new state the interaction must distinguish between the two Fock states, i.e. they have to interact differently. Since only the quarks interact in both Fock components the difference arises from their relative displacement in the transverse plane. If $`\rho `$ is the transverse separation between the quark and the photon, the $`\gamma ^{}q`$ fluctuation has a center of gravity in the transverse plane which coincides with the impact parameter of the parent quark. The transverse separation between the photon and the center of gravity is $`(1\alpha )\rho `$ and the distance between the quark and the center of gravity is correspondingly $`\alpha \rho `$. Therefore, the argument of $`\sigma _{q\overline{q}}`$ is $`\alpha \rho `$. More discussion can be found in .

The transverse momentum distribution of DY pairs can also be expressed in terms of the dipole cross section . The differential cross section is given by the Fourier integral

$`{\displaystyle \frac{d\sigma (qp\gamma ^{}X)}{d\mathrm{ln}\alpha d^2q_T}}`$ $`=`$ $`{\displaystyle \frac{1}{(2\pi )^2}}{\displaystyle d^2\rho _1d^2\rho _2\mathrm{exp}[\mathrm{i}\stackrel{}{q}_T(\stackrel{}{\rho }_1\stackrel{}{\rho }_2)]\mathrm{\Psi }_{\gamma ^{}q}^{}(\alpha ,\stackrel{}{\rho }_1)\mathrm{\Psi }_{\gamma ^{}q}(\alpha ,\stackrel{}{\rho }_2)}`$ (3)
$`\times `$ $`{\displaystyle \frac{1}{2}}\left\{\sigma _{q\overline{q}}(\alpha \rho _1)+\sigma _{q\overline{q}}(\alpha \rho _2)\sigma _{q\overline{q}}(\alpha (\stackrel{}{\rho }_1\stackrel{}{\rho }_2))\right\}.`$
after integrating this expression over the transverse momentum $`q_T`$ of the photon, one obviously recovers (1).

The LC wavefunctions can be calculated in perturbation theory and are well known . The dipole cross section on the other hand is largely unknown. Only at small distances $`\rho `$ it can be expressed in terms of the gluon density. However, several successful parameterizations exist in the literature, describing the entire function $`\sigma _{q\overline{q}}(x,\rho )`$, without explicitly taking into account the QCD evolution of the gluon density. We use the parameterization by Golec-Biernat and Wüsthoff  for our calculations, Fig. 2. This parameterization vanishes $`\rho ^2`$ at small distances, as implied by color transparency  and levels off exponentially at large separations.

In Fig. 2, we compare to E772 data  on low $`x_2`$ DY dilepton production . Most of the data are quite well described without any $`K`$-factor, which does not appear in this approach since higher order corrections are supposed to be parameterized in $`\sigma _{q\overline{q}}(\rho )`$. Moreover, the calculation in the dipole approach agrees with the next-to-leading order (NLO) parton model calculation at low $`x_2`$. Note that the dipole approach is valid only at low $`x_2`$ . At large $`x_2`$, this approach is not applicable and differs strongly from the parton model calculation. The disagreement between the data and both of the calculations in some points is probably due to systematic errors in the measured cross section. Preliminary E866 data  agree well with the NLO parton model calculation.

## 2 Proton-nucleus ($`pA`$) scattering

The main advantage of the dipole approach is its easy generalization to nuclear targets. Furthermore, it also includes some higher twist effects that are important in multiple scattering, and it provides insight into the physical mechanisms underlying nuclear effects, which are not easily accessible in the parton model .

Shadowing in DY is an interference phenomenon due to multiple scattering of the projectile quark inside the nucleus. In the target rest frame, where DY dilepton production is bremsstrahlung of massive photons, shadowing is the Landau-Pomeranchuk-Migdal (LPM) effect. These interferences occur (Fig. 1), because photons radiated at different longitudinal coordinates $`z_1`$ and $`z_2`$ are not independent of each other. Thus, the amplitudes have to be added coherently. Destructive interferences can occur only if the longitudinal distance $`z_2z_1`$ is smaller than the so called coherence length $`l_c`$, which is the time needed to distinguish between a quark and a quark with a $`\gamma ^{}`$ nearby. It is given by the uncertainty relation,

$$l_c=\frac{1}{\mathrm{\Delta }P^{}}=\frac{1}{m_Nx_2}\frac{(1\alpha )M^2}{q_T^2+(1\alpha )M^2+\alpha ^2m_q^2}.$$
(4)
Here, $`\mathrm{\Delta }P^{}`$ is the light-cone energy denominator for the transition $`qq\gamma ^{}`$ and $`q_T`$ is the relative transverse momentum of the $`\gamma ^{}q`$ Fock state. For $`z_1z_2>l_c`$, the radiations are independent of each other.

An immediate consequence of this is that $`l_c`$ has to be larger than the mean distance between two scattering centers in the nucleus ($`2`$ fm in the nuclear rest frame). Otherwise, the projectile quark could not scatter twice within the coherence length and no shadowing would be observed.

We develop a Green function technique , which allows one to resum all multiple scattering terms, similar to Glauber theory, and in addition treats the coherence length exactly. The formalism is equivalent to the one proposed in  for the LPM effect in QED. Our general expression for the nuclear DY cross section reads

$`{\displaystyle \frac{d\sigma (qA\gamma ^{}X)}{d\mathrm{ln}\alpha }}`$ $`=`$ $`A{\displaystyle \frac{d\sigma (qp\gamma ^{}X)}{d\mathrm{ln}\alpha }}{\displaystyle \frac{1}{2}}\mathrm{Re}{\displaystyle d^2b_{\mathrm{}}^{\mathrm{}}𝑑z_1_{z_1}^{\mathrm{}}𝑑z_2d^2\rho _1d^2\rho _2}`$ (5)
$`\times `$ $`\left[\mathrm{\Psi }_{\gamma ^{}q}(\alpha ,\rho _2)\right]^{}\rho _A(b,z_2)\sigma _{q\overline{q}}\left(\alpha \rho _2\right)G(\stackrel{}{\rho }_2,z_2|\stackrel{}{\rho }_1,z_1)`$
$`\times `$ $`\rho _A(b,z_1)\sigma _{q\overline{q}}\left(\alpha \rho _1\right)\mathrm{\Psi }_{\gamma ^{}q}(\alpha ,\rho _1).`$
The first term is just $`A`$ times the single scattering cross section, where $`A`$ is the nuclear mass number. The second term is the shadowing correction. The impact parameter is $`b`$ and the nuclear density is $`\rho _A`$. The Green function $`G`$ describes, how the bremsstrahlung-amplitude at $`z_1`$ interferes with the amplitude at $`z_2`$.

To make the meaning of Eq. (5) more clear, let us first consider a limiting case for $`G`$. In the simplest case, the coherence length, Eq. 4, is infinitely long and only the double scattering term is taken into account. Then $`G(\stackrel{}{\rho }_2,z_2|\stackrel{}{\rho }_1,z_1)=\delta ^{(2)}(\stackrel{}{\rho }_1\stackrel{}{\rho }_2)`$ and one of the $`\rho `$ integrations can be performed. The $`\delta `$-function means that at very high energy (infinite coherence length) the transverse size of the $`\gamma ^{}q`$ Fock-state does not vary during propagation through the nucleus, it is frozen due to Lorentz time dilatation. Furthermore, partonic configurations with fixed transverse separations in impact parameter space were identified a long time ago  in QCD as interaction eigenstates. This is the reason, why we work in coordinate space. Namely, in coordinate space, all multiple scattering terms can be resummed and in the limit of infinite $`l_c`$ one obtains

$$G^{frozen}(\stackrel{}{\rho }_2,z_2|\stackrel{}{\rho }_1,z_1)=\delta ^{(2)}(\stackrel{}{\rho }_1\stackrel{}{\rho }_2)\mathrm{exp}\left(\frac{\sigma _{q\overline{q}}(\alpha \rho _1)}{2}_{z_1}^{z_2}𝑑z\rho _A(b,z)\right).$$
(6)
The frozen approximation is identical to eikonalization of the dipole cross section in Eq. (1). Thus, the impact parameter representation allows a very simple generalization from a proton to a nuclear target, provided the coherence length is infinitely long.

At Fermilab fixed-target energies ($`\sqrt{s}=38.8`$ GeV for E772), this last condition is not fulfilled and one has to take a finite $`l_c`$ into account. The problem is however, that $`l_c`$, Eq. (4), depends on the relative transverse momentum $`q_T`$ of the $`\gamma ^{}q`$-fluctuation which is the conjugate variable to the size $`\rho `$ of this Fock-state and therefore completely undefined in $`\rho `$-representation. The quantum mechanically correct way to treat the $`q_T^2`$ in Eq. (4) is to represent it by a two-dimensional Laplacian $`\mathrm{\Delta }_T`$ in $`\rho `$-space. The Green function which contains the correct, finite coherence length and resums all multiple scattering terms fulfills a two-dimensional Schrödinger equation with an imaginary potential,

$`\left[\mathrm{i}{\displaystyle \frac{}{z_2}}+{\displaystyle \frac{\mathrm{\Delta }_T\left(\rho _2\right)\eta ^2}{2E_q\alpha \left(1\alpha \right)}}+{\displaystyle \frac{\mathrm{i}}{2}}\rho _A(b,z_2)\sigma _{q\overline{q}}\left(\alpha \rho _2\right)\right]G(\stackrel{}{\rho }_2,z_2|\stackrel{}{\rho }_1,z_1)`$
$`=\mathrm{i}\delta \left(z_2z_1\right)\delta ^{\left(2\right)}\left(\stackrel{}{\rho }_2\stackrel{}{\rho }_1\right),`$ (7)
where $`\eta ^2=(1\alpha )M^2+\alpha ^2m_q^2`$. For details of the derivation, we refer to .

The imaginary potential accounts for all higher order scattering terms. The Laplacian implies that the Green function is no longer proportional to a $`\delta `$-function. This means the size of the $`\gamma ^{}q`$ fluctuation is no longer constant during propagation through the nucleus. One can say that an eigenstate of size $`\rho _1`$ evolves to an eigenstate of size $`\rho _2\rho _1`$, so transitions between eigenstates occur.

Calculations with Eqs. (5) and (7) are compared to E772 data  in Fig. 3. Note that the coherence length $`l_c`$ at E772 energy becomes smaller than the nuclear radius. Shadowing vanishes as $`x_2`$ approaches $`0.1`$, because the coherence length becomes smaller than the mean internucleon separation. It is therefore important to have a correct description of a finite $`l_c`$ in this energy range. The eikonal (frozen) approximation, Eq. (6), does not reproduce the vanishing shadowing toward $`x_20.1`$. The curves in Fig. 3 are somewhat different from the ones in , because we used a different parameterization of the dipole cross section. Note that for heavy nuclei, energy loss  leads to an additional suppression of the DY cross section.

Nuclear effects on the $`q_T`$-differential cross section calculated at RHIC energy are shown in Fig. 4. See  for details of the calculation. The differential cross section is suppressed at small transverse momentum $`q_T`$ of the dilepton, where large values of $`\rho `$ dominate. This suppression vanishes at intermediate $`q_T2`$ GeV. The Cronin enhancement that one could expect in this intermediate $`q_T`$ region  is suppressed due to gluon shadowing .

A nuclear target provides a larger momentum transfer than a proton target and harder fluctuations are freed, which leads to nuclear broadening. Note, that not the entire suppression at low $`q_T`$ is due to shadowing. Some of the dileptons missing at low $`q_T`$ reappear at intermediate transverse momentum. At very large transverse momentum nuclear effects vanish.

## 3 Summary

We express the DY cross section in terms of the cross section for scattering a $`q\overline{q}`$ dipole off a proton. This is the same dipole cross section that appears in DIS. At low $`x_2`$ and for proton-proton scattering, calculations in the dipole approach agree with calculations in the NLO parton model. Some E772 data points are not well described by either of the approaches, which is probably due to a systematic error in the measured cross section.

At very high energy, the dipole approach is easily extended to nuclear targets by eikonalization. At lower fixed target energies (E772) the eikonal approximation is no longer valid, because the size of a Fock state varies during propagation through the nucleus. Therefore, transitions between interaction eigenstates (i.e. partonic configurations with fixed transverse separations) occur.

We develop a Green function technique, which takes variations of the transverse size into account and resums all multiple scattering terms as well. For light nuclei, calculations with the Green function technique are in good agreement with DY shadowing data from E772. For heavier nuclei, also energy loss becomes important.

We have also calculated nuclear effects in the transverse momentum distribution of DY pairs at RHIC energy. The DY cross section is suppressed at low transverse momentum. The expected Cronin enhancement at intermediate $`q_T2`$ GeV is reduced because of gluon shadowing. Nuclear effects vanish at very large $`q_T`$.

Acknowledgments: The author wishes to thank the Gesellschaft für Schwerionenforschung (GSI), Darmstadt, Germany for financial support during this summer school. This work was supported in part by the U.S. Department of Energy at Los Alamos National Laboratory under Contract No. W-7405-ENG-38.

I am indebted to Jörg Hüfner, Mikkel Johnson, Boris Kopeliovich, Jen-Chieh Peng and Alexander Tarasov for valuable discussion.


format: text

text:
# Effect of boundaries on the force distributions in granular media

## Numerical method

Our 2D packings consist of frictionless particles under gravity; the particles interact through normal Hertzian forces, where $`fd^{3/2}`$ and $`d`$ denotes the overlap distance hertz . Unless noted otherwise, the material constants and gravity are chosen such that a particle deforms $`0.1\%`$ under its own weight, and the particle radii are drawn from a flat distribution between $`0.4<r<0.6`$. Masses are proportional to the radii cubed. The container has a width of $`24`$, employs periodic boundary conditions in the horizontal direction and has a bottom consisting of a fixed hard support. The data shown in this paper were obtained from $`1100`$ realizations containing $`1180`$ particles each. We construct our stationary packings by letting the particles relax from a gas-like state by introducing a dissipative force that acts whenever the overlap distance $`d`$ is nonzero.

## Distribution of weights

In Fig. 1c we show the weight distribution $`𝒫(w)`$ for the bottom particles (dots) which differs qualitatively from the bulk distributions (open circles). Moreover, we observe that the transition is remarkably sharp: in the slice $`2<h<3`$, the weight distribution is already bulk-like (full curve). In the inset of Fig. 1c we plot $`w^2`$ which quantifies the width of $`𝒫(w)`$ as a function of height $`h`$. The sharp transition of $`𝒫(w)`$ near the bottom is clearly visible. Additionally, in the bulk, $`w^2`$ slowly increases with height, due to the deformations of the particles; this effect disappears for harder particles prep . Finally, near the top layer $`𝒫(w)`$ becomes sharply peaked and $`w^2`$ decreases.

To understand the change of $`𝒫(w)`$ near the bottom, consider the typical packing of Fig. 2a. The support aligns the bottom row of particles and thus strongly affects the directions of the interparticle forces. The forces between neighboring bottom particles are almost purely horizontal, so these hardly contribute to either the force on the particle from above \[the “weight” in (1)\] or to the force needed to support it. This approximation becomes better the smaller the polydispersity is. Thus, the average value of $`n_c`$, the number of particles that press on the particle from above and hence contribute to its weight, is on average lower at the bottom than in the bulk. Intuitively it is clear that the probability of finding a small value of $`w`$ increases with smaller $`n_c`$ for non-tensile forces. This statement can be made precise by considering Eq. (1) for fixed $`n_c`$ and analyzing $`𝒫_{n_c}(w)`$, the weight probabilities restricted to particles of given $`n_c`$. As long as the joint probability distribution of the interparticle forces remains finite for small forces, it follows from a phase-space argument that

$$𝒫_{n_c}(w)w^{n_c1}\mathrm{for}w0,$$
(2)
for all $`n_c1`$. The particles which do not feel a force from above, $`n_c=0`$, give a $`\delta `$–like contribution at $`W=1`$; for deep layers this occurs for $`w1`$.

To check the validity of this idea, we have determined $`𝒫_{n_c}(w)`$ both in the bulk and near the bottom by determining $`n_c`$ for each particle and decomposing the weight distribution $`𝒫(w)`$ into the $`𝒫_{n_c}(w)`$’s,

$$𝒫(w)=\underset{n_c}{}\rho _{n_c}𝒫_{n_c}(w).$$
(3)
The crucial difference between bottom and bulk are the fractions $`\rho _{n_c}`$ of particles that feel $`n_c`$ other particles pressing on them from above. This can be seen in Figs. 2b-c, where we show the decomposition of $`𝒫(w)`$ according to Eq. (3). This picture is also confirmed by Figs. 2d-e, which show that the individual distribution functions $`𝒫_1(w)`$ and $`𝒫_2(w)`$ do not differ significantly between bulk and bottom (to compare these, we have to normalize them not to the total average weight in each layer, but to the average weight of the particles with the same $`n_c`$). Note also that Figs. 2b-e show that Eq. (2) is valid, except for $`𝒫_1(w)`$ at the bottom for small values of $`w`$; this deviation comes from neglecting the intralayer, “almost horizontal” forces, whose small vertical components eventually make $`𝒫_1(w)0`$ for very small weights.

## Interparticle forces

The distributions of the interparticle forces are much more robust than the weight distributions. The results for the distribution $`P(f)`$ of $`|\stackrel{}{f}|`$ are shown in Fig. 3a. The only difference between bulk and bottom distributions is that the small peak around $`f=0.7`$ for bulk forces becomes a plateau for $`P(f)`$ for the forces close to the bottom (see Fig. 2a). It is intriguing to note that this change from a plateau to a peak is reminiscent of what is proposed as an identification of the jamming transition jamnote .

We also have measured the distribution of angles $`\phi _{ij}`$, which define the orientation of the $`\stackrel{}{f}_{ij}`$, and find that these angles are uniformly distributed and independent of the absolute value of $`\stackrel{}{f}`$ in the bulk, see Fig. 3c. Thus, in the bulk our packing is isotropic. Near the boundary, however, this isotropy is broken strongly: in agreement with our scenario for the influence of the packing geometry, the angles of the forces between bottom particles and those in the layer above are concentrated around $`\pi /3`$ and $`2\pi /3`$, as Figs. 2a and 3d show. Near the bottom, therefore, the interparticle forces naturally divide up into almost horizontal intralayer forces and “layer-to-layer” forces.

Since the weight distribution is determined by the $`z`$-components of the forces, let us also investigate the distribution $`P(f_z)`$. According to Fig. 3b, $`P(f_z)`$ also remains non-zero for small $`f_z`$, both in the bulk and near the bottom. There is a substantial difference, however, associated with the difference in packing. In the bulk, we have seen in Fig. 3c that there is no noticeable correlation between the angles $`\phi _{ij}`$ and the force strength $`f_{ij}`$. Hence in the bulk $`P(f_z)𝑑\phi 𝑑fP(f)P(\phi )\delta (f_zfsin(\phi ))`$ with $`P(\phi )=const`$. Indeed, the distribution obtained by numerical integration of this relation with $`P(f)`$ from Fig. 3a and a uniform angle distribution yields the solid line in Fig. 3b, which closely follows $`P(f_z)`$ as measured in the bulk. Near the bottom, on the other hand, the value of $`\mathrm{sin}(\phi )`$ is concentrated around $`\frac{1}{2}\sqrt{3}0.866`$: in the approximation that the distribution of $`\mathrm{sin}(\phi )`$ is sharply peaked at this value, the shape of $`P(f_z)`$ is close to that of $`P(f)`$, which is confirmed by direct comparison of the dotted datasets of panels a and b of Fig. 3.

## Manipulating $`𝒫(w)`$

In the previous paragraphs we have shown that the weight distribution $`𝒫(w)`$ is very sensitive to the local packing geometry, while the distribution of interparticle forces is robust. This allows one to manipulate $`𝒫(w)`$ at the bottom by changing the boundary conditions. We illustrate this by simulations of weakly polydisperse particles, $`0.49<r<0.51`$, both with a flat bottom (Fig. 4a) and a curved bottom, consisting of two circle segments of radius 20 glued together (Fig. 4b). For the flat bottom, the particles form an almost perfect hexagonal packing, leading to particles with mostly $`n_c=2`$ (horizontal contacts again not included). In agreement with (2), $`𝒫(w)`$ increases linearly for small $`w`$ in this case. The weakly curved bottom locally disturbs this crystalline structure (Fig. 4b), causing a dramatic change in the fractions $`\rho _{n_c}`$, and correspondingly in $`𝒫(w)`$ (Fig. 4c). Note that the small difference between the distributions of the convex and the concave part of the bottom are reflected in the $`\rho _n`$ as well. Interestingly, $`P(f)`$ is in all these cases indistinguishable from $`P(f)`$ in the strongly polydisperse case prep .

## Perspective

The message that emerges naturally from the above analysis is clear: in experiments in which the forces on a boundary are probed, one measures effectively the weights $`w`$. These weights, however, are not the most fundamental quantities of a granular packing, as they are derived from the interparticle forces $`\stackrel{}{f}_{ij}`$. These capture the full microscopic structure, and the distribution function $`P(f)`$ of the force strength is quite insensitive to the packing, in contrast to $`𝒫(w)`$. Simple phase space considerations show that a grain with $`n_c`$ contacts with particles that press downwards on it makes a contribution to the weight distribution $`𝒫(w)`$ which scales as $`w^{n_c1}`$ as $`w0`$. Thus, the small weight distribution is dominated by particles with few such contacts, in particular by those with just one, and the change in geometry near boundaries leads to an atypical $`𝒫(w)`$. In addition, we occasionally observe a small peak at $`w1`$ due to “loose” particles ($`n_c=0`$) which do not feel a force from above; there are indications for such a peak of $`𝒫(w)`$ at $`w=0`$ in recent precise experiments (see Fig. 5 of network2 ).

In the standard $`q`$-model, weights are randomly distributed over a fixed number of neighbors one layer below; in the simplest version there are two such neighbors that receive a fraction $`q`$ and $`1q`$ of the weight qm . Due to a fixed connectivity, this model cannot be expected to capture the behavior of $`𝒫(w)`$, especially near boundaries. The product of the weight $`w`$ and $`q`$, however, could be interpreted as an interparticle force. Interestingly, in the simplest case of a uniform probability distribution of the $`q`$’s, the probability distribution $`P(qw)`$ is a pure exponential prep ; snoeijer . In simulations of the $`q`$-model with random connectivities, a variety of $`𝒫(w)`$’s can be obtained, similar to what we found here for the frictionless spheres prep .

It is possible to test our framework in experiments. For $`𝒫(w)`$ we expect that curvature effects, such as shown in Fig. 4, only play a role when they break highly ordered packings; indeed carbon-paper measurements at the sidewalls in cylinders give a very similar $`𝒫(w)`$ as near the bottom network . A simpler way to change the boundary conditions may be to include a layer of larger particles at the boundary; their value of $`n_c`$ will be higher, and we expect that $`𝒫(w)`$ for small $`w`$ will decrease for larger $`n_c`$. Furthermore, there is a strong need of direct determination of $`P(f)`$, both in the bulk and near the boundaries, since the present data concerns $`𝒫(w)`$ only network ; network2 ; exp2 . An interesting observation in the context of “jamming” jamnote  is that in our simulations the peak in $`P(f)`$ appears to vanish near the boundaries (Fig. 3a); are granular materials no longer jammed here, and is this relevant for the localization of shear bands near boundaries?

An important issue for future study is clearly the role of friction and dimensionality. Our numerical study has been done in two dimensions with frictionless spheres; however, recent studies indicate makse  that the coordination number for 3D packings with friction is similar to those of 2D frictionless packings. Qualitatively, the picture we have advanced is therefore expected to capture the realistic case of three dimensions with friction, because our phase space arguments are independent of dimension.

We finally note an interesting open issue. In the $`q`$-model, $`𝒫(w)`$ approaches its asymptotic expression as a function of depth algebraically slow snoeijer ). In our packings, the convergence appears to be much faster — is this a real discrepancy and if so is it related with the same packing issues?

We are grateful to Martin Howard and Hans van Leeuwen for numerous illuminating discussions.


text:
# Quaternionic Hamilton equations

## 1. Introduction

In the last years, a better understanding of the geometrical framework of the classical hamiltonian dynamics led to some interesting generalizations. This is for instance the case of the Nambu dynamics (, , ), where the evolution of the system in a $`n`$\- dimensional manifold is related to $`(n1)`$ hamiltonian functions, or of the Liouville dynamics, whose properties are studied in .
Very recently, in  the authors proposed a different extension of the classical hamiltonian dynamics to spaces which are endowed by three symplectic structures, having the same commutation rules as the quaternionic imaginary units. Roughly speaking, this new hyper–hamiltonian dynamics is the sum of three different hamiltonian dynamics, each one corresponding to a single symplectic structure and to an arbitrary hamiltonian function. The hyper–hamiltonian dynamics retains most of the appealing features of the classical hamiltonian one, included a variational formulation (, ).

Here the reader will not be faced with any new application of the hyper–hamiltonian dynamics, nor with any discussion about its properties. The aim of this paper is to point out that the extension given in , which is natural from the geometrical point of view, is also natural from the point of view of analysis. More precisely, if the classical hamiltonian dynamics is reinterpreted in the light of the complex analysis, and then extended by replacing the complex analysis with the quaternionic one, what one gets is again the same hyper–hamiltonian dynamics introduced in .

It should be stressed that the former step does not consist in writing the Hamilton equations by using complex coordinates, but rather in having the linear notion underlying complex holomorphy (i.e. complex linearity) coming into play. This is done in section 3.
When the complex numbers are replaced by the quaternions, it is well known that the complex holomorphy correspondingly changes in the so called quaternionic regularity (, ). This is done in section 4, where the results in  are used to tune the translation process to our needs.
Finally, in section 5 the linear quaternionic regularity replaces the complex linearity into the equation governing the dynamics. To conclude, it is shown that the resulting equations coincide with the one already obtained in  by a completely different approach.
Section 2 is devoted to introduce the right context where the previous program can be worked out. The real needs are metrics and complex structures, which calls for Kähler manifolds instead of just symplectic one, as it may be expected.

Notations. The symbols $``$ and $``$ denote the complex numbers and the quaternions; their standard bases are $`1,i`$ and $`1,i_1,i_2,i_3`$ respectively. Repeated indexes, like in $`i_\alpha /x^\alpha `$, allude to a summation on $`\alpha =1,2,3`$; when misunderstandings are possible, explicit summations or some extra information are added. Finally, if $`X`$ is a vector field over a differential manifold $`M`$, by $`X\text{  }\text{ }t`$ one intends its contraction with the tensor field $`t`$.

Acknowledgements. The authors would like to thank G. Gaeta for many helpful discussions on the subject.

## 2. General setting

Let consider a Riemannian manifold $`M`$ with a metric $`g`$. On $`M`$, take a $`g`$–invariant complex structure $`J`$, namely a linear map $`TMTM`$ such that:

$$J^2=I\text{and}g(JX,JY)=g(X,Y)X,Y.$$

The manifold $`M`$ is said to be a Kähler manifold when the 2–form defined by:

(2.1)
$$\omega (X,Y)=g(JX,Y)$$

is closed. Since $`\omega `$ is nondegenerate and anti–symmetric, it defines a symplectic structure on $`M`$. Given a function $`h:M`$, the Hamilton equations may be introduced:

(2.2)
$$X\text{  }\text{ }\omega =dh.$$

Of course one can define Hamilton equations on any symplectic manifold, regardless its symplectic structure originates from a metrics or it doesn’t. However, this restriction is essential for all the discussion in the next section.
In  the Hamilton equations are extended to cover the case where $`M`$ is an hyper–Kähler manifold. Instead of a single one, $`TM`$ has now three $`g`$–invariant complex structures such that:

(2.3)
$$J_1^2=J_2^2=IJ_1J_2=J_2J_1J_3=J_1J_2,$$

and all the 2–forms $`\omega _\alpha `$ associated to $`J_\alpha `$ via (2.1) are closed (see e.g. ). The hyper–hamiltonian mechanics introduced in  then depends on the choice of three distinct hamiltonian functions $`h^\alpha `$, its dynamics being defined by the following equations:

(2.4)
$$X\text{  }\text{ }\mathrm{\Omega }=\frac{1}{(2n1)!}_\alpha dh^\alpha (\omega _\alpha )^{2n1}.$$

where $`\mathrm{\Omega }`$ is the volume form on $`M`$, whose dimension is $`4n`$, and $`(\omega _\alpha )^{2n1}`$ denotes the wedge product of the differential form $`\omega _\alpha `$ by itself for $`2n1`$ times.
A decomposition formula for the vector field $`X`$ is then proved in :

(2.5)
$$X=X_1+X_2+X_3\text{where}X_\alpha \mathrm{}\omega _\alpha =dh^\alpha \alpha ,$$

suggesting the most general hyper–hamiltonian dynamics is the sum of three independent hamiltonian dynamics.

## 3. Hamilton equations and complex analysis

Consider first classical hamiltonian mechanics. By making use either of the metrics $`g`$ and of the complex structure $`J`$ on $`TM`$, one can rewrite the Hamilton equations (2.2) in the form:

$$g(X,Y)=dh(JY)Y.$$

The key point is now to give a closer look at the 1–form lying at the right hand side, namely:

$$\theta (Y)=dh(JY)$$

from the point of view of complex analysis. First of all, notice that $`M`$ is an almost complex manifold in a natural way. Indeed defining

$$(zY)_p:=tY_p+xJ_pY_p\text{where}z=t+ix$$

makes $`T_pM`$ a complex linear space, for every $`pM`$. As a consequence, it makes sense to talk about complex linear maps on it, and one easily realizes that $`\theta `$ is the only real valued 1–form for which the complex valued 1–form $`\theta +idh`$ is in fact a complex linear map. Such a completion property deserves the ad hoc notation:

$$\widehat{idh}=\theta .$$

The obvious identity $`idh=d(ih)`$ then suggests the following approach. First of all think of the pure imaginary complex function:

$$H=ih$$

as the hamiltonian function describing the system, and then rewrite the Hamilton equations in the equivalent form:

(3.1)
$$X\text{  }\text{ }g=\widehat{dH}.$$

Until now, nothing really new. However, everyone can now guess how to extend the approach from complex numbers to quaternions. Assume $`M`$ is an hyper–Kähler manifold, and just think of the hamiltonian as a function:

$$H:MH=h^\alpha i_\alpha $$

with pure imaginary values (summation on $`\alpha `$ is implicit). Then think of $`\widehat{dH}`$ as the completion of its differential with respect to the natural notion of holomorphy in the framework of quaternionic analysis: the so called regularity (, ). This will be done in the next section, where the following formula:

$$\widehat{dH}(Y)=dh^\alpha (J_\alpha Y)$$

will be proved.
Equation (3.1) makes sense also in this new framework, and selects a unique vector field $`X`$. The important point is that, maybe surprisingly, $`X`$ is exactly the same vector field defined by (2.4). This is the main result of the paper, and it will be proved in the last section.

## 4. Regularity and completion to a regular form

In 1935 Fueter proposed a definition of right regularity for a function $`\phi :`$, by adapting the Cauchy–Riemann equations to the increased number of independent imaginary units in $``$, namely ():

(4.1)
$$\frac{\phi }{t}+\frac{\phi }{x^\alpha }i_\alpha =0q=t+x^\alpha i_\alpha .$$

By putting the units $`i_\alpha `$’s on the left, one obtain the so called left regular functions, whose theory, clearly, is completely equivalent to the one of right regular functions.
To a large extent, the quaternionic analysis based on this notion mimics the richness of complex analysis (see (), though the non commutativity of $``$ introduces some important differences: for instance, the product or the composition of regular functions need not to be regular.
Notice that equation (4.1) is in fact a condition on the differential of $`d\phi `$. Indeed, $`\phi `$ is right regular if and only if $`d\phi `$ belongs to a special class of 1–forms $`T`$, namely the $`\theta `$’s which satisfy the condition:

(4.2)
$$\theta (1)+_\alpha \theta (i_\alpha )i_\alpha =0.$$

As pointed out in the introduction, the interest here is essentially in a more general type of 1–forms, that is:

$$\theta :TM.$$

where $`M`$ is an hyper–Kähler manifold. What does it mean regularity in this framework, is a question which may be answered using the results in ().
In that paper, the notion of regularity is extended to functions acting between general quaternionic spaces, and the extension is shown to be strictly dependent on if $``$ acts on the involved spaces either on the left or on the right. For instance, in  is proved that the above considered notions of regularity make sense only when $``$ acts differently on the two spaces. In particular, the needs for the right regularity are a domain which is a left space, and a codomain which is a right one.
Now, the way $``$ acts on the domain $`TM`$ of $`\theta `$ is in fact implicit in our definition of hyper–Kähler manifold. Indeed, due to the commutation rules (2.3), the position:

$$(qY)_p=tY_p+x^\alpha (J_\alpha )_pY_p\text{where}q=t+x^\alpha i_\alpha $$

makes every $`T_pM`$ a left quaternionic space. Hence, as pointed out above, the codomain $``$ has to be considered as a right space over itself, and the only applicable notion of regularity is the right one.
Following , the correct definition of regularity for $`\theta `$ is then:

(4.3)
$$\theta (Y)+_\alpha \theta (i_\alpha Y)i_\alpha =0Y.$$

That (4.3) is really an extension of (4.2) depends on the quaternionic anti–linearity of its left–hand member (). Because of that indeed, when $`M=`$ the equality in (4.3) holds for any $`YT`$ if and only if it holds for $`Y=1`$.
Before going further, a remark is in order about the invariance of the notion of regularity. Looking at (4.3), this notion seems to depend on the choice of the imaginary units $`i_\alpha `$’s. Now, in $``$ there is an $`S^2`$ of imaginary units, and any choice of a triple of units, having the right commutation properties, is an admissible one: a natural quaternionic object must be independent on that choice. The following lemma says this is the case for the notion of regularity.

###### Lemma 4.1.

If the triple $`j_1`$, $`j_2`$ and $`j_3`$ of quaternions is chosen in such a way that:

$$j_1^2=j_2^2=1j_3=j_1j_2,$$

then, for any quaternionic valued 1–form $`\theta `$, the following identity holds:

$$_\alpha \theta (i_\alpha Y)i_\alpha =_\alpha \theta (j_\alpha Y)j_\alpha Y.$$

###### Proof.

Define a $`3\times 3`$ real matrix $`C`$ by means of:

$$j_\alpha =C_\alpha ^\beta i_\beta \alpha .$$

The commutation rules in the hypothesis say that $`C^TC=I`$. The conclusion follows from straightforward computations.∎
In the next proposition the regularity will be exploited to introduce the notion of completion to a regular form, for 1–forms $`\xi `$ which take their values in the imaginary quaternions, namely:

(4.4)
$$\xi :TM\text{with}\xi +\overline{\xi }=0.$$

###### Proposition 4.2.

If $`\xi `$ satisfies condition (4.4), then there is a unique real valued 1–form $`\widehat{\xi }T^{}M`$ such that $`\widehat{\xi }+\xi `$ is a regular 1–form. Moreover, the completion $`\xi \widehat{\xi }`$ is a real linear operation and, if $`\xi =\xi ^\alpha i_\alpha `$ with $`\xi ^\alpha T^{}M`$, then:

$$\widehat{\xi }(Y)=\xi ^\alpha (i_\alpha Y)Y.$$

The proof is a straightforward consequence of the following lemma, where regularity of the general 1–form $`\theta `$ is rewritten using its components.

###### Lemma 4.3.

Assume $`\theta =\theta ^0+\theta ^\alpha i_\alpha `$, where $`\theta ^0`$ and each $`\theta ^\alpha `$ lie in $`T^{}M`$ (namely are real valued 1–forms). Then $`\theta `$ is regular if and only if:

$$\theta ^0(Y)=\theta ^\alpha (i_\alpha Y)Y.$$

###### Proof.

Just write down the sum in (4.3), separating the contributions along $`1`$ and each $`i_\alpha `$. The four conditions obtained in this way are all equivalent: the statement expresses the vanishing along $`1`$. ∎

###### Remark 4.4.

Like the notion of regularity, the notion of completion is also independent on the choice of the imaginary units $`i_\alpha `$.

## 5. Quaternionic Hamilton equations

Coming back to dynamics, take an hamiltonian function:

$$H:M\mathrm{with}\overline{H}+H=0.$$

on the hyper–Kähler manifold $`M`$. Proposition 4.2 applies to the 1–form $`dH`$, and its completion may be used to introduce the Hamilton like equation:

(5.1)
$$X\text{  }\text{ }g=\widehat{dH}.$$

###### Remark 5.1.

In some special situations, the quaternionic Hamilton equations coincide with the classical ones. This is the case when the hamiltonian function is directed along a single pure imaginary quaternion, namely:

$$H(p)=h(p)uh\text{real valued,}u^2=1$$

Indeed, by means of Proposition 4.2 and the remark thereafter:

$$\widehat{dH}(Y)=dh(uY)Y,$$

so that equation (5.1) becomes:

$$X\text{  }\text{ }\omega _u=dh.$$

Here, of course, $`\omega _u`$ is the symplectic form associated to the complex structure $`J_uY=uY`$.

To complete the program announced in the introduction, it remains to show that both the equations (2.4) and (5.1) define the same vector field $`X`$. This is a consequence of the following result and formula 2.5.

###### Proposition 5.2.

Assume $`H=h^\alpha i_\alpha `$. Then the vector field $`X`$ satisfies equation (5.1) if and only if it decomposes as:

$$X=X_1+X_2+X_3$$

where:

$$X_\alpha \text{  }\text{ }\omega _\alpha =dh^\alpha \alpha .$$

###### Proof.

Using Proposition 4.2, condition (5.1) rewrites as:

$$g(X,Y)=dh^\alpha (i_\alpha Y)Y.$$

The decomposition suggested in the statement clearly gives rise to a a vector field $`X`$ which satisfies the previous condition. On the other hand, there is only one vector field which has this property. ∎
Summing up, two possible extensions of classical hamiltonian mechanics have been considered. One is the hyper–hamiltonian mechanics introduced in , whose dynamics is defined by (2.4), and which is natural from the point of view of symplectic geometry. The other one is the quaternionic hamiltonian mechanics introduced in this note, whose dynamics is defined by (5.1), and which is natural from the point of view of complex and quaternionic analysis. Though they arise from apparently unrelated approaches, the two extension coincide.
Moreover, in the quaternionic approach the exact meaning of the sum in (2.5) also becomes clear. The intrinsic object is the hamiltonian function $`H`$, whereas the $`h^\alpha `$ are just its components with respect to a given basis in the space of imaginary quaternions. Changing the base, the components of $`H`$ change, and so do the vector fields $`X_\alpha `$: however, its sum is intrinsic inasmuch as it does not depend on any choice.


text:
# Contact tracing and epidemics control in social networks

\[

## Abstract

A generalization of the standard susceptible-infectious-removed (SIR) stochastic model for epidemics in sparse random networks is introduced which incorporates contact tracing in addition to random screening. We propose a deterministic mean-field description which yields quantitative agreement with stochastic simulations on random graphs. We also analyze the role of contact tracing in epidemics control in small-world networks and show that its effectiveness grows as the rewiring probability is reduced.

\]

Properties of complex networks recently attracted much attention in physical community . Although perhaps it was prompted by the advent of the Internet and World-Wide Web, the importance of this subject goes far beyond computer networks. Indeed, daily commute, power and goods traffic, wired and wireless communication, disease spreading occur within certain physical or social networks. The theory of disease spreading, which is known as mathematical epidemiology, has a long and rich history (see, e.g. ). However, until recently the epidemiological studies have been mostly concerned with so-called mean-field description of epidemics, in which it is assumed that at any time the probability to get infected is the same for all individuals. In some other works, spreading of a disease on relatively simple lattices of individuals has been studied within so called “forest-fire” models . Only recently, the studies which elucidate the role of underlying network structure in the disease spreading began to appear in the literature .

Most of the epidemiological models are based on several simple assumptions regarding disease contracting and cure. In particular, the most common mechanism of infection is through a contact with another infectious individual, and the mechanism of recovery is either deterministic or purely stochastic with a certain typical time of recovery. In the simplest Susceptible-Infectious-Susceptible (SIS) model, a recovered individual immediately becomes susceptible again, while in a more complicated Susceptible-Infectious-Removed (SIR) model, cured individuals become immune and effectively excluded from further dynamics.

While these models give a good description of evolution of many common infectious diseases, they usually neglect the role of intelligent strategies to stop nascent epidemics. Few epidemiological models take into account prevention strategies such as, for example, mass and ring vaccination . In practice, one of the main counter-epidemics measures is the contact tracing, when individuals which have been in contact with an infected (and identified) individual, are found and thoroughly checked. It applies, among others, to the treatment of sexually transmitted diseases, tactics of law-enforcement organizations trying to uncover criminal or terrorists networks, cleaning of computer virus infection, etc. We only are aware of one theoretical paper  where a model of this kind has been studied. The model  is based on the assumption that infection is a slow branching process, while the contract tracing occurs at a much shorter time scale. This leads to a familiar SIR-type model with rescaled parameters and similar dynamics. In this Letter we consider a more realistic model in which infection and contact tracing occur concurrently, and their interplay determines the dynamics of the system.

Stochastic model. We assume that the population consists of $`N`$ hosts whose connections to each other form a fixed graph. The hosts are enumerated with index $`n=1,\mathrm{},N`$. A node $`n`$ is said to have a degree $`k(n)`$ if it is connected to $`k`$ other hosts. In case of random graphs the degree distribution is Poissonian with a certain mean degree $`K=k(n)`$.

For simplicity we assume that there in no spontaneous recovery, an infectious individual can only be disinfected externally through screening. Immediately upon disinfecting, the individual becomes traced (T) for a certain period of time during which its neighbors are checked for possible infection. After that time, the individual spontaneously becomes removed, and its neighbors are no longer traced.

Infection $`SI`$. Initially, the whole population except for one host is assumed to be susceptible to infection. The probability of host infection depends on the state of its nearest neighbors. The infection dynamics is modeled as a simple contact process: if a susceptible node $`n`$ has $`k_i(n)`$ infectious neighbors, the probability that it becomes infectious during a small $`\mathrm{\Delta }t`$ time interval is $`\alpha k_i(n)\mathrm{\Delta }t`$.

Tracing $`IT`$. The process of infection elimination consists in finding infectious hosts and then curing them. Hosts are being checked with certain probability $`\beta `$ which depends on the state of its neighbors. We postulate that if an infectious host is checked, it is immediately cured, eliminated or at least isolated so it cannot infect other hosts. We introduce two non-exclusive strategies of checking for infectious hosts: random checking and contact tracing. Random checking means choosing an arbitrary host with probability $`\beta _r\mathrm{\Delta }t`$, while contact tracing of host $`n`$ is done with probability $`\beta _tk_t(n)\mathrm{\Delta }t`$ where $`k_t(n)`$ is the number of neighbors of $`n`$ which are in the traced state $`T`$. The random checking process is equivalent to the removal process of general epidemics.

Removal $`TR`$. With certain probability $`\gamma \mathrm{\Delta }t`$, traced hosts are transformed into the removed state, in which they also cannot be infected, but they are no longer under observation, so they do not initiate contact tracing.

Stochastic simulations of the described process were performed using an event-driven scheme. It is significantly superior over synchronous and asynchronous update schemes both in terms of accuracy and computational speed. In the event-driven scheme we select the time lapsed between two consecutive events from a Poisson distribution with a combined probability of all events (infection, tracing, recovery, etc.), then choose a node (each node has its own probability to be chosen depending on its state and the states of its neighbors), and apply the transition from one state to another according to the ratio of individual transition probabilities.

In our simulations with random graph based networks, we typically built networks with average degree $`K=10`$ and $`1000`$ nodes. For every random graph we ran $`100`$ simulations starting every time from a single (but different in each run) infected host. Then we averaged the results for $`50`$ random graphs. The time bin was $`\mathrm{\Delta }t=10^6`$ for most simulations. In all simulations we varied the tracing parameters $`\beta _r`$ and $`\beta _t`$, while the infection constant was set at $`\alpha =0.1`$, and the transition rate from $`T`$ to $`R`$, $`\gamma =0.5`$. The latter parameter is important for the effectiveness of the targeting elimination, because the longer a node remains in the traced state, the more probable it is to trace its neighboring infectious nodes. However, since tracing presumably bears a significant cost, an optimal choice of the tracing parameters ($`\beta _r,\beta _t`$, and $`\gamma `$) for a given epidemics is an important issue.

In Fig. 1 we present the “prevalence” of epidemics (the fraction of infectious nodes in the whole population $`i=I/N`$) as a function of time for several values of $`\beta _r`$ and $`\beta _t`$. When $`\beta _r=\beta _t=0.0`$ we have a simple $`SI`$ process, and all the nodes eventually get infected (thick solid line in Fig. 1). Other curves show the fraction of infectious nodes as a function of time for $`\beta _r=0.02`$ and different values of $`\beta _t`$. The ratio $`\alpha /\beta _r`$ is chosen to be above the epidemics threshold . For $`\beta _t=0`$ we obtain the classical $`SIR`$ process with randomly removed infectives (dashed line in Fig. 1). The epidemic eventually saturates, and the fraction of infectious nodes decays exponentially. The lower lines display the evolution of the infection fraction for values of $`\beta _t`$ ranging from $`0`$ to $`2.5`$ with a step value of $`0.1`$. The initial (exponential) phase of the epidemics growth is nearly independent of $`\beta _t`$, because the contact tracing process is intrinsically nonlinear (it requires the presence of I-T connected pairs and therefore only begins after the first infected node is randomly screened). As expected, the tracing process significantly reduces the magnitude of the epidemics (maximal value of $`i`$), but at large times the infection decays with the same exponential rate as for $`\beta _t=0.0`$ (again, because we return to the linear regime at small $`i`$). The most interesting feature of the process at large $`\beta _t>0.35`$ is the presence of a second maximum of $`i`$ which indicates a second epidemic. Due to this second epidemic, the percentage of the infectious population at large times $`t>40`$ may actually increase with increase of $`\beta _t`$. It means that the range of $`\beta _t`$ values from $`0.4`$ to $`0.9`$ are not better to control the epidemics than values smaller than $`0.4`$.

Mean-field equations. At the first sight, it seems that the mean-field approach cannot be applied to the contact tracing, since it does not take into account the non-uniform distribution of infection in the population. Nevertheless, a more sophisticated mean-field approach which operates not only with the mean densities of states, but also with the densities of links connecting nodes with different states, can be devised. In the derivation we follow the method of Rand (see also ). Let us introduce the number of nodes $`A`$, the number of connected pairs $`[AB]`$ and triples $`[ABC]`$ of nodes, where $`A,B`$ and $`C`$ stand for any of the types $`S,I,T,R`$. For example, the number of connected pairs of infectious and traced nodes is denoted $`[IT]`$. Note that $`[AB]=[BA]`$ and each pair in $`[AA]`$ is counted twice. For large $`N`$, the ratios $`A/N,[AB]/N,[ABC]/N`$ approach deterministic limits which we label $`a,[ab],[abc]`$, respectively.

The dynamics of the model is described by the following set of rate equations

$`\dot{s}`$ $`=`$ $`\alpha [si],\dot{i}=\alpha [si]\beta _ri\beta _t[i\tau ],`$ (1)
$`\dot{[ss]}`$ $`=`$ $`2\alpha [ssi],`$ (2)
$`\dot{[si]}`$ $`=`$ $`\alpha ([ssi][isi][si])\beta _r[si]\beta _t[si\tau ],`$ (3)
$`\dot{[ii]}`$ $`=`$ $`2\alpha ([isi]+[si])2\beta _r[ii]2\beta _t[ii\tau ],`$ (4)
$`\dot{[i\tau ]}`$ $`=`$ $`\beta _r[ii]+\beta _t([ii\tau ][\tau i\tau ][i\tau ])\gamma [i\tau ]\beta _r[i\tau ].`$ (5)
Here we used the notation $`\tau =[T]/N`$ to avoid confusion between the density of traced nodes and time $`t`$. Note that we omit here the equations for $`\tau ,[\tau \tau ],[s\tau ]`$, as well as any combinations involving the removed state, since they do not affect the dynamics of the infectious population.

The meaning of these equations is rather straightforward. For example, the terms in the r.h.s. of the last equation can be explained as follows. A $`(p,q)`$ pair becomes $`[i\tau ]`$ through random screening of the infectious node $`q`$ in an $`[ii]`$ pair, or through contact tracing of node $`q`$ from node $`r`$ in a $`(p,q,r)`$ triple in the state $`[ii\tau ]`$. On the other hand, we can lose an $`[i\tau ]`$ pair by contact tracing of the $`p`$ node in a $`[\tau i\tau ]`$ triple $`(r,p,q)`$, by direct tracing $`p`$ by $`q`$, by removing of $`q`$, or by random screening of $`p`$. Other equations can be obtained from similar arguments.

This set of equation is not closed, as the equations for the pair densities contain triple densities. We need to introduce a closure rule. Similarly to Refs., we can use the approximation $`[abc]=[ab][bc]/b`$, which follows from the condition that the influence of a node on the state of its second neighbor in a triple is negligible.

Using this closure rule, we arrive at the following set of equations

$`\dot{s}`$ $`=`$ $`\alpha s\widehat{i},`$ (6)
$`\dot{i}`$ $`=`$ $`\alpha s\widehat{i}\beta _ri\beta _ti\widehat{\tau },`$ (7)
$`\dot{\widehat{i}}`$ $`=`$ $`(\alpha Ks\alpha \beta _r)\widehat{i}\beta _t\widehat{i}\widehat{\tau },`$ (8)
$`\dot{\widehat{\tau }}`$ $`=`$ $`(\beta _t\widehat{p}\beta _t\gamma \alpha {\displaystyle \frac{s}{i}}\widehat{i})\widehat{\tau }+\beta _r\widehat{p},`$ (9)
$`\dot{\widehat{p}}`$ $`=`$ $`\alpha {\displaystyle \frac{s}{i}}(2\widehat{i}+2\widehat{p})\widehat{i}(\beta _r+\beta _t\widehat{\tau })\widehat{p}.`$ (10)
where $`\widehat{i}=[is]/s`$ is the mean number of infectious neighbors per susceptible node, $`\widehat{\tau }=[i\tau ]/i`$ is the mean number of traced neighbors per infectious node, and $`\widehat{p}=[ii]/i`$ is the mean number of infectious neighbors of an infectious node. Notice that the equation for $`[ss]`$ dropped out as $`[ss]=Ks^2`$ at all times. We used the initial conditions $`s(0)=1i_0,i(0)=i_0,\widehat{i}(0)=(K1)i_0,\widehat{\tau }(0)=0,\widehat{p}(0)=0`$ which correspond to a small set of disconnected infectious nodes.

During the early stage of an epidemic the contact tracing can be neglected ($`\tau =0`$), and Eqs. (6)-(10) are reduced to a set of three equations for $`s,i,\widehat{i}`$ which coincide with the model has been studied in . Independently of $`\beta _t`$, the initial epidemics growth is characterized by the basic reproduction number $`K\alpha /(\alpha +\beta _r)`$. However, as the number of traced individuals grows, the growth rate is reduced and the epidemics is saturated. In Fig.2 the dynamics of the epidemics calculated from Eqs.(6)-(10) are shown for different values of $`\beta _t`$. As seen from the figure, the maximum number of infectious nodes is drastically reduced with increase of $`\beta _t`$. In the same figure we show the results of direct stochastic simulations for $`\beta _t=0,0.5`$. The most important question is whether the contact tracing is capable of arresting the exponential growth of the epidemics before it engulfs a finite portion of the total population.

To answer this question, we consider the limit of small epidemics in a large populations ($`i,\widehat{i}<<1`$), then we can set $`s(t)=1`$ and drop Eq.(7). We also observe that $`\widehat{i}/i=K1`$ and drop Eq.(7) as well. A simple calculation shows that the critical value of $`\beta _t`$ at which the exponential growth of infection is arrested,

$$\beta _{cr}=\frac{(\alpha (K1)+\gamma )(\alpha (K1)\beta _r)}{\beta _r}.$$
(11)
For $`\beta >\beta _{cr}`$, epidemic remains small at all times, and so there is no major outbreak of the epidemic, the maximum number of infectious nodes is independent of the population size, and depends on the initial size of infection $`i_0`$. For the parameter values of our stochastic simulations, $`\alpha =0.1,K=10,\beta _r=0.02,\gamma =0.5`$, we obtain $`\beta _{cr}=61.6`$. In Fig.3 we show the maximum number of infectious nodes during the epidemic as a function of $`\beta _t`$ at different initial epidemic sizes $`i_0`$. In agreement with the above argument, for $`\beta _t<\beta _{cr}`$ the value of $`i_{max}`$ is almost independent on $`i_0`$ until $`i_{max}2i_0`$ at a certain $`\beta _t^{}`$, after which it remains nearly constant. The value of $`\beta _t^{}`$ approaches $`\beta _{cr}`$ for $`i_00`$.

It is easy to see that the epidemic threshold depends only on the average node degree and not on the specific topology of the underlying graph. However, the subsequent development of the epidemics should be significantly affected by the network structure, in particular by the average minimum path length and the clustering coefficient. Low average minimum path means that any node in the graph can be exposed to infection in a short time after an epidemics begins. If the clustering coefficient is large, the infection propagates faster within a certain community but then it make it easier to trace the epidemics. Sparse random networks studied above represent a particular class of networks with short average minimal path and small clustering coefficient. Many social networks are characterized by a relatively large clustering coefficient while keeping the average minimal path low. We studied numerically the effect of the network structure on the contact tracing of epidemics within the small-world model . Changing the re-wiring probability $`p`$ allows us to scan the range of networks from regular ($`p=0`$) to random ($`p1`$) through the small-world range

($`0.001<p<0.1`$) which exhibits a short average minimal path and a large clustering coefficient typical for many social networks. We used the same number of nodes and edges as for the random graph simulations, and fixed the parameter values at $`\alpha =0.1,\beta _r=0.02,`$ and $`\gamma =0.5`$. Fig.4 shows the dependence of the maximum prevalence $`i_{max}`$ on $`p`$ for several different $`\beta _t`$. As we can see, $`i_{max}`$ changes mostly within the SW range ($`0.001<p<0.1`$) where the clustering coefficient and the average path undergo large variations.

In this Letter we studied the role of contact tracing as a part of the epidemics control strategy in complex networks. We demonstrated that by applying this strategy, a major outbreak can be significantly reduced or even eliminated at a small additional cost. Based on the pair correlation approach due to Rand , we developed the mean-field model of contact tracing for the case of random graphs. We also studied the influence of network topology on the contact tracing using the small-world model with variable re-wiring probability $`p`$, and found that its effectiveness grows as the reviring probability is reduced. The main change occurs within the small-world regime at $`p10^2`$.

This work was supported by the U.S. DOE grant DE-FG03-96ER14592, ARO MURI-ARO grant DAAG55-98-0269, MCyT BFI2000-157, and NATO Collaborative Linkage Grant PST.CLG.978512.


text:
# Strange Disoriented Chiral Condensate

## I Introduction

Relativistic nuclear collisions can produce matter in which chiral symmetry is restored. One possible consequence of the restoration and the subsequent re-breaking of chiral symmetry is the formation of disoriented chiral condensates (DCC) – transient regions in which the average chiral order parameter differs from its value in the surrounding vacuum DCCreview ; qmrev1 ; qmrev2 . Previous efforts to describe DCC signals have focused on pion production.

In Ref. GavinKapusta , Kapusta and I have explored the possible influence of DCC on kaon production, inspired by the explanation due to Kapusta and Wong KapustaWong  of measurements of $`\mathrm{\Omega }`$ and $`\overline{\mathrm{\Omega }}`$ baryon enhancement SPSbaryons  at $`17`$ $`A`$ GeV at the CERN SPS in terms of the production of many small DCC regions within individual collision events. If true, Ref. KapustaWong  implies that the evolution of the condensate can have a significant effect on strange particle production. The importance of strange degrees of freedom in describing chiral restoration has been long appreciated PisarskiWilczek ; Columbia ; GGP2 ; Karsch ; Lenaghan , but simulations of the three flavor linear sigma model had suggested that strange kaon fields are much less important than the pion fields Randrup . Nevertheless, the $`\mathrm{\Omega }`$ and $`\overline{\mathrm{\Omega }}`$ data demand that we explore without prejudice techniques for measuring kaon fluctuations.

A further implication of Ref. KapustaWong  is that the DCC regions must be rather small, with a size of about $`2`$ fm. Such a size is consistent with predictions based on dynamical simulations of the two flavor linear sigma model GGP1 . More importantly, the DCC search for anomalous event structure in neutral and charged pions by the WA98 Collaboration at the SPS has revealed no evidence of large DCC ”domains” qmrev2 . Note that I loosely use the word ”domain” to refer to spatial regions in which the condensate is somehow coherent. I stress that no thermodynamically stable domain structures are expected in theoretical descriptions of DCC.

In this paper I study kaon isospin fluctuations in the presence of many small DCC. In the next section I discuss how DCC may lead to kaon fluctuations. Following Ref. GavinKapusta , I compute the probability distribution that describe the DCC contribution to these fluctuations and combine the DCC fluctuations with a contribution from a random thermal background. Pion fluctuations due to many small DCC have been addressed by Amado and Lu AmadoLu  and Chow and Cohen ChowCohen . Here we focus on kaon fluctuations; pion as well as kaon results are presented in Ref. GavinKapusta . In sec. 4 I assess robust statistical observables that can be used to measure the impact of many small DCC at RHIC and LHC. In sec. 5 I present work in progress on event generator simulations to understand the magnitude and centrality dependence of the statistical observables in the absence of DCC. In particular, I obtain a dynamic isospin fluctuation observable analogous to the dynamic charge observable used to measure net charge fluctuations at RHIC Voloshin . Of the quantities considered, this observable isolates the DCC effect from other sources of fluctuations best.

## II Strange DCC

To illustrate how a strange DCC can form, let us first consider QCD with only up and down quark flavors. Equilibrium high temperature QCD respects chiral symmetry if the quarks are taken to be massless. This symmetry is broken below $`T_c150`$ MeV by the formation of a chiral condensate $`\sigma \overline{u}u+\overline{d}d`$ that is a scalar isopin singlet. However, chiral symmetry implies that $`\sigma `$ is degenerate with a pseudoscalar isospin triplet of fields with the same quantum numbers as the pions. In reality, chiral symmetry is only approximate and the 140 MeV pion mass is different from the $`800\pm 400`$ MeV mass of the leading sigma candidate pdg . Nevertheless, lattice calculations exhibit a dramatic drop of $`\sigma `$ near $`T_c`$ at finite quark masses.

A DCC can form when a heavy ion collision produces a high energy density quark-gluon system that then rapidly expands and cools through the critical temperature. Such a system can initially break chiral symmetry along one of the pion directions, but must then evolve to the $`T=0`$ vacuum by radiating pions. A single coherent DCC radiates a fraction $`f_\pi `$ of neutral pions compared to the total that satisfies the probability distribution

$$\rho _1(f_\pi )=\frac{1}{2f_\pi ^{1/2}}\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0}<f_\pi 1,$$
(1)
Anselm ; Blaizot ; Bjorken . Such isospin fluctuations constitute the primary signal for DCC formation in the pion sector. The enhancement of baryon-antibaryon pair production is a secondary effect due to the relation between baryon number and the topology of the pion condensate field KapustaSrivastava .

This two flavor idealization only applies if the strange quark mass $`m_s`$ can be taken to be infinite. Alternatively, if I take $`m_s=m_u=m_d=0`$, then the chiral condensate would be an up-down-strange symmetric scalar field. The more realistic case of $`m_s100`$ MeV is between these extremes, so that $`\sigma \mathrm{cos}\theta (\overline{u}u+\overline{d}d)+\mathrm{sin}\theta (\overline{s}s)`$. The mixing angle $`\theta `$ is highly uncertain since it depends on the sigma mass together with the $`\pi ,K,\eta `$ and $`\eta ^{}`$ masses and the $`\eta \eta ^{}`$ mixing angle GGP2 . A disoriented condensate can evolve by radiating $`\pi ,K,\eta `$ and $`\eta ^{}`$ mesons, with the neutral pion fraction satisfying (1). Randrup and Schäffner-Bielich find that the kaon fluctuations from a single large DCC satisfy Randrup

$$\rho _1(f_K)=\mathrm{1\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0}f_K1,$$
(2)
where $`f_K=(K^0+\overline{K}^0)/(K^++K^{}+K^0+\overline{K}^0)`$. Moreover, the condensate fluctuations can now produce strange baryon pairs KapustaWong . Linear sigma model simulations indicate that pion fluctuations dominate three-flavor DCC behavior, while the fraction of energy imparted to kaon fluctuations is very small due to the kaons’ larger mass. On the other hand, domain formation may be induced by other mechanisms such as kaon condensation at high baryon density Kaplan , bubble formation qm95  or decay of the Polyakov loop condensate Pisarski .

## III DCC Mesons from Many Small Domains

Why does the DCC’s size matter? Pion measurements in individual collision events can distinguish DCC isospin fluctuations from a thermal background only if the disoriented region is sufficiently large qmrev1 . DCC can then be the dominant source of pions at low transverse momenta, since $`p_t1/R`$ for a coherent region of size $`R`$. Experiments focusing on low $`p_t`$ can study neutral and charged pion fluctuations Bjorken , wavelet Ina  and HBT signals qmrev1 ; HBT  to extract detailed information. In contrast, for small domains ($`R<3`$ fm qmrev1 ) DCC signals are hidden by fluctuations due to ordinary incoherent production mechanisms. This holds even if many such regions are produced per event. DCC mesons from small regions may have momenta of a few hundred MeV, nearer the $`pp`$ mean value. Different regions would not add coherently to alter HBT, nor would their small spatial structures affect wavelet analysis.

Importantly, baryon pair enhancement KapustaWong  is substantial only if there are many small incoherent regions. The large winding numbers that produce baryon-antibaryon pairs require many small regions with random relative orientations of the pion field. To describe strange antibaryon enhancement, Kapusta and Wong assume roughly 100 DCC regions of size roughly $`2`$ fm KapustaWong . Topological models of baryon-antibaryon pair production successfully describe $`e+e`$ and hadronic collision data EllisKowalski . The connection of DCC to topological pair production was pointed out in Ref. KapustaSrivastava ; see also DeGrand .

To compute the distribution of kaons due to many small DCC regions, define $`f=(K^0+\overline{K}^0)/(K^++K^{}+K^0+\overline{K}^0)`$. To an excellent approximation the number of neutral kaons is equal to twice the number of short-lived neutral kaons $`K_S`$ which are more readily measurable in high energy heavy ion collisions. The fraction $`f`$ ranges from 0 to 1. The statistical distribution in $`f`$ for a single domain is $`\rho _1(f)=1`$. The distribution for $`n`$ randomly oriented, independent domains is

$$\rho _n(f)=\underset{k=1}{\overset{n}{}}df_k\rho _1(f_k)\delta \left(f\frac{1}{n}\underset{j=1}{\overset{n}{}}f_j\right).$$
(3)
In GavinKapusta  I obtain

$$\rho _n(f)=n^2\underset{0k<n(1f)}{}(1)^k\frac{[n(1f)k]^{n1}}{k!(nk)!}.$$
(4)
In the limit that $`n1`$, this distribution tends toward a Gaussian of mean $`f=1/2`$ and standard deviation $`\sigma =\{12n\}^{1/2}`$. Results for pions are presented in GavinKapusta .

In a more realistic scenario some kaons will come from the decay or realignment of DCC domains and some will come from more conventional sources. I shall refer to the latter as random or thermal, even though that may be a bit of a misnomer. What I mean by random or thermal is that the distribution of kaons from non-DCC sources is

$$\rho _0(f_0)=\frac{1}{2\pi \sigma _0^2}\mathrm{exp}\left[(f_01/2)^2/2\sigma _0^2\right].$$
(5)
For a completely random source the width $`\sigma _0`$ is related to the total number $`N_{\mathrm{random}}`$ of non-DCC kaons by

$$\sigma _0^2=\frac{1/2(11/2)}{N_{\mathrm{random}}}=\frac{1}{4N_{\mathrm{random}}}.$$
(6)
Now let us assume that a fraction $`\alpha _K`$ of all kaons come from non-DCC sources and the remaining fraction $`\beta _K=1\alpha _K`$ come from $`n1`$ independent DCC domains. Letting $`N`$ denote the total number of kaons, I have $`N_{\mathrm{random}}=\alpha _KN`$ and $`N_{\mathrm{DCC}}=\beta _KN`$. Folding together two Gaussians gives a Gaussian.

$`\rho _K(f)`$ $`=`$ $`{\displaystyle 𝑑f_0𝑑f_n\rho _0(f_0)\rho _n(f_n)\delta (f\alpha _Kf_0\beta _Kf_n)}`$ (7)
$`=`$ $`{\displaystyle \frac{1}{\sqrt{2\pi \mathrm{\Delta }_K^2}}}\mathrm{exp}\left[(f1/2)^2/2\mathrm{\Delta }_K^2\right]`$
For a thermal source plus $`n`$ DCC domains, the net width is

$$\mathrm{\Delta }_K^2=\frac{\alpha _K}{4N}+\frac{\beta _K^2}{12n}=\frac{1}{4N}+\left\{\frac{\beta _K^2}{12n}\frac{\beta _K}{4N}\right\}.$$
(8)
The expression in curly brackets at the end represents the difference between the actual width and the width the distribution would have if there was no contribution from DCC kaons. This change in the width may be positive or negative, depending on the parameters.

## IV Statistical Analysis

Detection of small incoherent DCC regions in high energy heavy ion collisions requires a statistical analysis in the $`\pi ^0\pi ^\pm `$ or the $`K_S^0K^\pm `$ channels. Neutral mesons can be detected by the decays $`\pi ^0\gamma \gamma `$ or $`K_S^0\pi ^+\pi ^{}`$. The analysis we propose in GavinKapusta  is sensitive to correlations due to isospin fluctuations. We expect these correlations to vary when DCC regions increase in abundance or size as centrality, ion-mass number $`A`$, or beam energy are changed. Correlation results combined with other signals, such as baryon enhancement KapustaWong , can be used to build a circumstantial case for DCC production.

Correlations of $`\pi ^0\pi ^\pm `$ and $`K_S^0K^\pm `$ can be determined by measuring the robust isospin covariance,

$$R_{c0}=\frac{N_cN_0N_cN_0}{N_cN_0},$$
(9)
where $`N_0`$ and $`N_c`$ are the number of neutral and charged mesons. I take $`N_0=N_{\pi ^0}`$ and $`N_c=N_{\pi ^+}+N_\pi ^{}`$ for pion fluctuations and $`N_0=2N_{K_S^0}`$ and $`N_c=N_{K^+}+N_K^{}`$ for kaon fluctuations. The ratio (9) has two features that are convenient for experimental determination. First, this observable is independent of detection efficiency as are the “robust” ratios discussed in Taylor . Robust observables are useful for DCC studies because charged and neutral particles are identified using very different techniques and, consequently, are detected with different efficiency. Observe that robust quantities are not affected by the unobserved $`K_L^0`$, since the strong-interaction eigenstates $`K^0`$ and $`\overline{K}^0`$ are a superposition $`K_L^0`$ and $`K_S^0`$ until their decay well outside the collision region. Second, since (9) is obtained from a statistical analysis, individual $`\pi ^0\gamma \gamma `$ or $`K_S^0\pi ^+\pi ^{}`$ need not be fully reconstructed in each event. This feature is crucial because it would be extraordinarily difficult – if not impossible – to reconstruct a low momentum $`\pi ^0`$ in heavy ion collisions except on a statistical basis.

Next I define robust variance

$$R_{aa}=\frac{N_a^2N_a^2N_a}{N_a^2},$$
(10)
where $`a=c`$ or 0. To see why (10) is robust, denote the probability of detecting each meson $`ϵ`$ and the probability of missing it $`1ϵ`$. For a binomial distribution the average number of measured particles is $`N_a^{\mathrm{exp}}=ϵN_a`$ while the average square is $`N_a^2^{\mathrm{exp}}=ϵ^2N_a^2+ϵ(1ϵ)N_a`$. I then find

$$R_{aa}^{\mathrm{exp}}=R_{aa},$$
(11)
independent of $`ϵ`$ Pruneau ; the proof that (9) is robust is similar. The ratios (9) and (10) are strictly robust only if the efficiency $`ϵ`$ is independent of multiplicity. Further properties and advantages of these and similar quantities are discussed in Pruneau .

To study DCC fluctuations I define the dynamic isospin observable

$$\nu _{\mathrm{dyn}}^{c0}=R_{cc}+R_{00}2R_{c0}.$$
(12)
Analogous observables have been employed to study net charge fluctuations in particle physics Whitmore ; Boggild  and were considered in a heavy ion context in Voloshin  and Mrowczynski . This quantity can be written in terms of

$$\nu ^{c0}=\left(\frac{N_0}{N_0}\frac{N_c}{N_c}\right)^2.$$
(13)
To isolate the dynamical isospin fluctuations from other sources of fluctuations, one obtains (12) by subtracting from (13) the uncorrelated Poisson limit $`\nu _{\mathrm{stat}}^{c0}=N_0^1+N_c^1`$. Indeed, we show in (19) below that the quantity (12) depends primarily on the fluctuations of the neutral fraction $`f`$, while the individual ratios (9) and (10) have additional contributions.

I illustrate the effect of DCC on the dynamic isospin fluctuations by writing $`N_0=fN`$ and $`N_c=(1f)N`$. Small fluctuations on $`f`$ or $`N`$ results in the changes

$$\frac{\mathrm{\Delta }N_0}{N_0}=\frac{\mathrm{\Delta }N}{N}+\frac{\mathrm{\Delta }f}{f},\frac{\mathrm{\Delta }N_c}{N_c}=\frac{\mathrm{\Delta }N}{N}\frac{\mathrm{\Delta }f}{1f}.$$
(14)
I obtain the average

$$\frac{\mathrm{\Delta }N_0^2}{N_0^2}=v+\frac{2c}{Nf}+\frac{\mathrm{\Delta }^2}{f^2}.$$
(15)
Here the contribution of the variance of the total number of mesons is $`v\mathrm{\Delta }N^2/N^2`$ and the charge-total covariance is $`c\mathrm{\Delta }N\mathrm{\Delta }f`$. DCC formation primarily effects the charge fluctuation contribution, $`\mathrm{\Delta }^2(\mathrm{\Delta }f)^2`$, from (15) or (17). Similarly,

$$\frac{\mathrm{\Delta }N_c^2}{N_c^2}=v\frac{2c}{N(1f)}+\frac{\mathrm{\Delta }^2}{(1f)^2},$$
(16)
and

$$R_{c0}=v+\left(\frac{1}{f}\frac{1}{1f}\right)\frac{c}{N}\frac{\mathrm{\Delta }^2}{(1f)^2}$$
(17)
where $`R_{c0}`$ is given by (9). Using (21) I get

$$\nu _{\mathrm{dyn}}^{c0}=\frac{1}{f(1f)}\left(\frac{\mathrm{\Delta }^2}{f(1f)}\frac{1}{N}\right).$$
(18)
This observable isolates the isospin fluctuations, whereas the individual $`R_{ab}`$ depend on the fluctuations in total meson number, $`v`$ and $`c`$ as well.

I estimate the effect of DCC on the dynamical fluctuations (18) using (6) and (8). I take $`N=N_K`$ for kaons and $`N=N_\pi `$ for pions; these are the total number of mesons of the indicated kind. For kaons

$$\nu _{\mathrm{dyn}}^{c0}(\mathrm{K}\mathrm{DCC})=4\beta _K\left(\frac{\beta _K}{3n}\frac{1}{N_K}\right).$$
(19)
These quantities can be positive or negative depending on the magnitude of $`\beta `$ compared to the number of domains per kaon. In fact the dynamical fluctuation may even be positive for one kind of meson and negative for the other.

## V Fluctuations in the Absence of DCC

Let us now discuss work in progress in which Abdel-Aziz and I use event generators to simulate conventional sources of kaon fluctuations hamlet . In the absence of DCC, $`\alpha =1`$ and $`\beta =0`$ so that (19) implies $`\nu _{\mathrm{dyn}}^{c0}0`$. On the other hand, incomplete equilibration may result in dynamical correlations in nuclear collisions not described by (5). Little is known from $`pp`$ experiments about kaon fluctuations. Event generators such as HIJING and URQMD models both yield negative values of $`\nu _{\mathrm{dyn}}^{c0}`$ in pp collisions. HIJING simulations of central Au+Au at 200 $`A`$ GeV in the rapidity range $`0.5<y<0.5`$ yield $`\nu _{\mathrm{dyn}}^{c0}0.002`$ for 47 $`K^+`$ and 44 $`K_S^0`$ on average HIJING .

The onset of DCC formation can substantially change the value of $`\nu _{\mathrm{dyn}}^{c0}`$. To search for this onset it is useful to study fluctuations as a function of collision centrality. I use HIJING to estimate the influence of conventional collision geometry and dynamics on the centrality dependence. However, I find that one can understand the HIJING results quite simply using the wounded nucleon model. In a multiple collision models such as the wounded nucleon model, one describes a nucleus-nucleus collision as a superposition of M independent nucleon-nucleon sub-collisions. The robust variance and covariance satisfy

$$R_{ab}=\frac{r_{ab}}{M}+\frac{M^2M^2}{M^2},$$
(20)
where the $`r_{ab}`$ are coefficients describing the fluctuations in each sub-collision and $`a,b`$ can equal either $`c`$ or $`0`$. The dynamical isospin observable,

$$\nu _{\mathrm{dyn}}^{c0}=R_{cc}+R_{00}2R_{c0}=\frac{\nu _0}{M},$$
(21)
is independent of the contribution from the fluctuations of $`M`$. In the wounded nucleon model, $`M`$ is given at each impact parameter by the number of participant nucleons.

In figs. 1 and 2, I show the kaon covariance (9) and dynamical isospin observable (12) as functions of the number of participants $`M`$ computed from 50,000 HIJING events for Au+Au at 200 A$`GeV`$ in the rapidity range $`0.5<y<0.5`$ compared to wounded nucleon model calculations (20). The number of participants at impact parameter $`b`$ for a symmetric Au+Au collision is computed using $`M(b)=2𝑑sT(s)\{1e^{\sigma _{NN}T(bs)}\}`$, where $`T(b)=\rho (z,b)𝑑z`$ is the familiar nuclear thickness function and $`\rho `$ is the three-parameter Fermi nuclear density for Au. In fig. 1, the solid curve is determined using (20) with a coefficient $`r_{c0}`$ computed from 50,000 S+S collisions, while the dashed curve is obtained by varying $`r_{c0}`$ to fit Au+Au. The HIJING results scale roughly as $`R_{c0}M^1`$ as expected, but the difference from the wounded nucleon model are rather large.

Figure 2 shows the dynamical isospin observable as a function of $`M`$ from HIJING. The solid curve is obtained from (21) with $`\nu _0`$ determined from HIJING S+S. I find that the wounded nucleon model is in excellent agreement with HIJING, suggesting that the correlations that increase $`R_{ab}`$ compared to the wounded nucleon model are similar for all kaon charge states. The agreement of HIJING and the wounded nucleon model in fig. 2 is likely due to the following. First, baryon stopping is unimportant at the high RHIC energy. Different numbers of protons and neutrons would alter the isospin balance. Second, high $`p_T`$ aside, the way HIJING describes soft interactions is quite similar to the wounded nucleon model, since it doesn’t incorporate final state cascading. To understand the effect of cascading, we are currently studying URQMD collisions hamlet .

To illustrate the possible scenario for the onset of DCC effects, I assume that DCC kaons add to the kaons from multiple sub-collisions according to

$$\nu _{\mathrm{dyn}}^{c0}=\beta ^2\nu _{dcc}+(1\beta ^2)\nu _{wnm},$$
(22)
where $`\beta `$ is the fraction of DCC kaons, $`\nu _{dcc}`$ is the DCC contribution given by (10) and $`\nu _{wnm}`$ is given by (21). I assume that DCC production above an impact parameter $`b_0`$ exhibits a threshold behavior, $`\beta =\beta _0[1(b/b_0)^2]`$, where $`b_0`$ and $`\beta _0`$ are ad hoc constants. In fig. 4, I show estimates assuming that 10 domains contribute kaons in the range $`0.5<y<0.5`$ for $`b_06`$ fm, varying the dcc fraction $`\beta _0`$ between 10% and 20%. Note that the DCC contribution to $`\nu `$ is positive for these values.

## VI Discussion and Conclusion

Reference KapustaWong  argued that the anomalous abundance and transverse momentum distributions of $`\mathrm{\Omega }`$ and $`\overline{\mathrm{\Omega }}`$ baryons in central collisions between Pb nuclei at 17 $`A`$ GeV at the CERN SPS is evidence that they are produced as topological defects arising from the formation of many domains of disoriented chiral condensates (DCC) with an average domain size of about 2 fm. Motivated by this interpretation, Kapusta and I studied the effect of DCC on the distribution of the fractions of neutral kaons and pions in GavinKapusta .

The DCC pioneers Anselm ; Blaizot ; Bjorken ; DCCreview  had hoped that a large percentage of pions might be emitted from just a few big domains, on the order of 5 to 8 fm (kaons were not considered). Such large domains have been ruled out at SPS qmrev2 , but remain possible at RHIC. More conservatively, as the number of domains grow and their average size diminishes, the impression left on the fluctuations in the neutral fraction becomes more subtle and less unique. For many small domains, statistical measurements of both neutral kaons (pions) and charged kaons (pions) are needed to observe the rather small isospin fluctuations. We have identified robust observables for that purpose. In particular, we have shown that the dynamical isospin observable (12) is sensitive to DCC but not to thermal or multiple-collision sources of fluctuations, as discussed in the text near eqs. (19) and (21) respectively. While HIJING simulations support this conclusion, more work remains. For example, one can use URQMD to study the effect of final state scattering on this observable. Isospin fluctuations can appear as changes in the magnitude of the dynamical isospin observable as centrality is varied. We emphasize that similar consequence may follow from any mechanism that produces many small domains that decay to pions and kaons, such as the Polyakov Loop Condensate Pisarski .

I am grateful to Joe Kapusta for a very enjoyable collaboration and thank M. Abdel-Aziz, R. Bellwied, C. Pruneau, S. Voloshin and, especially, H. Stöcker for many useful discussions. This work was supported in part by U.S. DOE Grant number DE-FG02-92ER40713.