text:
# Shock geometry and inverse Compton emission from the wind of a binary pulsar

## 1 Introduction

PSR B1259$``$63 is one of only three known radio pulsars which have a main sequence star binary companion. Such systems provide a unique environment for inverse Compton scattering because of the presence of an enormous density of low energy (optical) photons which serve as targets for electrons and positrons in the pulsar wind. The modest spin-down luminosity, and relatively high distance to two of these systems together imply that they are unlikely to be detectable sources of inverse Compton emission. In contrast, PSR B1259$``$63 is a galactic radio pulsar with a high spin-down luminosity ($`L_\mathrm{p}=8.3\times 10^{28}`$W) and is relatively nearby ($`1.5`$kpc). If just 0.1% of the wind luminosity is scattered into hard $`\gamma `$-ray photons the resulting flux should be detectable using current $`\gamma `$-ray telescopes.

PSR B1259$``$63 is in a highly eccentric orbit ($`e0.87`$) around SS2883, a B2e star of radius $`R_{}6R_{}`$ and luminosity $`L_{}8.8\times 10^3L_{}`$ \[Johnston et al. 1992, 1994, 1996\]. At periastron the pulsar is only $`23R_{}10^{11}\mathrm{m}`$ from its companion. Around periastron there are $`3.4\times 10^{11}\mathrm{cm}^3`$ Be-star photons at the pulsar, corresponding to an energy density of $`U_{\mathrm{rad}}6\times 10^{11}\mathrm{eV}\mathrm{cm}^3`$. This is some 11 orders of magnitude larger than the typical background target density available for inverse Compton scattering by the winds of isolated pulsars.

The Crab is presently the only pulsar for which observations place any constraints on the physical parameters of a pulsar wind. Observations of the PSR B1259$``$63 system aimed at detecting hard $`\gamma `$-ray inverse Compton emission are planned during 2000. If successful they should provide the first direct probe of the freely-expanding region of the wind of any rotation-powered pulsar, and serve to constrain wind models and parameters.

In the absence of an efficient deceleration process the relativistic wind of a pulsar will expand freely until it attains pressure balance with the surrounding medium. If the wind is still supersonic at the point of pressure balance it will be bounded by a termination shock. Kirk, Ball & Skjæraasen  investigated the $`\gamma `$-ray emission from PSR B1259$``$63 resulting from inverse Compton scattering by the shocked pulsar wind, downstream of the termination shock, where the pulsar wind electrons and positrons have been accelerated and isotropised as the radial flow of the wind is disrupted. The effects of the cooling of the shocked wind by inverse Compton scattering were included by Tavani & Arons  in their model for the synchrotron emission from the system. Ball & Kirk \[2000a\] considered the effects of inverse Compton scattering by the unshocked pulsar wind, upstream of the termination shock. They calculated the deceleration of the wind that occurs as a result of ‘inverse Compton drag’ as energy and momentum are transferred from the wind to the scattered photons. In particular it was shown that the inverse Compton losses were unlikely to contain the wind of PSR B1259$``$63 before it attains the radius at which pressure balance with the Be-star outflow is likely to occur. Ball & Kirk \[2000a\] then calculated the emission from the freely-expanding wind subject to the assumption that the wind was not terminated by pressure balance.

In this paper we recalculate the inverse Compton emission from the unshocked wind of PSR B1259$``$63 including the effects of termination of the radially-expanding pulsar wind as a result of pressure balance with the Be-star wind. The relative strengths of the pulsar and Be-star winds are unknown. If the Be-star wind dominates then the termination shock will be close to, and wrapped around, the pulsar. This will have the effect of decreasing the emission from the unshocked portion of the wind, and the effects vary over the binary orbit of the system. We hope that inverse Compton $`\gamma `$-rays will be detected from PSR B1259$``$63 in the near future, providing data that can be compared with such models.

The calculation of the geometry of the termination shock is outlined in §2. In §3 we present $`\gamma `$-ray spectra and light curves from the terminated wind, and compare them with those calculated for the unterminated wind by Ball & Kirk \[2000a\]. Our conclusions are presented in §4.

## 2 Shock geometry

The winds of the pulsar and its companion star will generally be separated by a pair of termination shocks separated by a contact discontinuity. The shock structures enclose the star with the weaker wind. Giuliani  developed a general description of axisymmetric flows that can be applied to the interaction of two such winds. Girard & Willson \[1987; hereafter GW87\] followed this description and derived the equations that describe the shape and position of the boundary separating two stellar winds subject to the following assumptions:

1. The shape of the shock is not affected by the orbital motion of the two stars. This implies azimuthal symmetry about the line joining the two stars.
2. The shock region has negligible thickness.
3. Both winds are emitted radially with no angular dependence.
4. The winds are non-relativistic.

The resulting set of three ordinary differential equations (ODEs) can be integrated directly to find the position of the boundary, which can be loosely referred to as the ‘shock’, given the assumption of negligible thickness. The winds are characterised by their mass loss rate $`\dot{M}`$ and (constant) radial velocity $`u`$ and these enter the equations that describe the shock geometry only in terms of the ratios

$`m={\displaystyle \frac{\dot{M_1}}{\dot{M_2}}}\mathrm{and}w={\displaystyle \frac{u_1}{u_2}}`$ (1)
where the subscripts refer to the different stars. GW87 argued that the results obtained depended only on the quantity $`\eta =mw`$, although $`m`$ and $`w`$ enter the ODEs separately. In an alternative formulation presented by Huang & Weigert  the quantities $`m`$ and $`w`$ only enter the ODEs as their product $`\eta `$.

Assumption (4) listed above is almost certainly not valid for a pulsar wind, which is likely to be highly relativistic. Modelling of the nebula produced by the wind of the Crab pulsar implies that the wind has a bulk Lorentz factor of $`\gamma 10^6`$ at radii well beyond the light cylinder \[Kennel & Coroniti 1984\]. Melatos, Johnston & Melrose  considered the interaction of a relativistic wind of PSR B1259$``$63 comprised primarily of electrons and positrons, with a wind from SS2883 dominated by much more massive ions. They showed that when the ions are dominant in both number and energy density, as is expected in PSR B1259$``$63/SS2883, the dependence of the geometry only on the product of $`m`$ and $`w`$ is lost. Nevertheless, for a given value of $`\eta `$ the opening angles calculated by Melatos, Johnston & Melrose  are within a factor of $`2`$ of the results of GW87 for a wide range of values of $`w/m`$.

In the following we approximate the shock position in the PSR B1259$``$63/SS2883 system using the formulation of GW87. This minimises the number of poorly known wind parameters that affect the results. The shock position is represented by its radial distance from the pulsar $`r_\mathrm{T}(\alpha )`$ where $`\alpha `$ is the angle between the lines joining the pulsar to the shock and the pulsar to the Be star.

Figure 1 shows the position of the termination shock (TS) obtained from numerical solution of the ODEs of GW87 for six different values of $`\eta `$. The quantity $`\eta `$ is the ratio of the rates at which the stars are transferring momentum to their winds. When star 1 is taken to be the pulsar and star 2 is the Be star it follows that

$$\eta =\frac{L_p/c}{\dot{M}v}$$
(2)
where $`\dot{M}`$ and $`v`$ are the mass loss rate and speed of the Be-star wind, and $`L_\mathrm{p}`$ is the spin-down luminosity of the pulsar. When $`\eta <1`$ the Be-star wind dominates the pulsar wind, and the TS wraps around the pulsar. As $`\eta `$ increases the apex of the TS moves further from the pulsar, the opening angle increases, and for $`\eta >1`$ the pulsar wind is dominant and the TS wraps around the Be star.

The standoff distance from the pulsar to the apex of the TS is

$`r_\mathrm{T}(0)={\displaystyle \frac{\sqrt{\eta }}{1+\sqrt{\eta }}}D`$ (3)
where $`D`$ is the stellar separation. At large distances from the stars the momenta of the winds are almost parallel and the shock tends asymptotically to a cone characterised by a half-opening angle $`\psi `$ and distance to the apex $`\rho _\mathrm{c}`$. A useful empirical approximation for $`\psi `$ as a function of $`\eta `$ is \[Eichler & Usov 1993\]:

$`\psi =2.1\left(1{\displaystyle \frac{\overline{\eta }^{\frac{2}{5}}}{4}}\right)\overline{\eta }^{\frac{1}{3}}`$ (4)
where $`\overline{\eta }=\mathrm{min}(\eta ,\eta ^1)`$. Figure 2 shows the geometry of the termination shock and illustrates the parameters used to describe it. We use $`\theta `$ to denote the value of $`\alpha `$ corresponding to the line of sight from Earth to the pulsar; i.e. $`\theta `$ is the angle between the line of sight and the line from the pulsar to the Be star.

The Be star companion of PSR B1259$``$63 possesses an excretion disk which the pulsar passes through near periastron \[Johnston et al. 1996, 1999; Ball et al. 1999\]. This disk will dominate the termination of the pulsar wind for 50 days or so near periastron, greatly complicating the physics of the pulsar wind system at these epochs. Its effects on unpulsed X-ray emission close to periastron have been considered by Tavani & Arons . However its presence is irrelevant for the majority of the 1237 day binary orbit and it is neglected here. The results presented below can therefore not be taken to be accurate at epochs close to periastron.

## 3 Inverse Compton scattering

Ball & Kirk \[2000a; Equation 15\] derived an expression for the radiation transfer that determines the intensity of the radiation emitted at the normalised energy $`ϵ_{\mathrm{out}}=E_{\mathrm{out}}/m_\mathrm{e}c^2`$. Furthermore, it was shown that absorption of the scattered photons due to pair production on the Be-star photons is negligible. It follows that when the termination of the wind is included, the scattered intensity can be written as a line of sight integral of the scattering due to electrons moving radially out from the pulsar,

$`I(ϵ_{\mathrm{out}})={\displaystyle _0^{r_\mathrm{T}(\theta )}}\left({\displaystyle \frac{L_\mathrm{p}}{4\pi }}\right){\displaystyle \frac{ϵ_{\mathrm{out}}}{\gamma _\mathrm{w}(0)m_\mathrm{e}c^3}}{\displaystyle \frac{\mathrm{d}N_\gamma }{\mathrm{d}ϵ_{\mathrm{out}}\mathrm{d}t}}ds,`$ (5)
where the fraction of the wind momentum carried by ions is assumed to be negligible. The quantity $`\gamma _\mathrm{w}(0)`$ is the initial Lorentz factor of the pulsar wind and $`\mathrm{d}N_\gamma /\mathrm{d}ϵ_{\mathrm{out}}\mathrm{d}t`$ is the differential rate of emission of inverse Compton scattered photons by a single electron. In the calculations of Ball & Kirk \[2000a\] the TS was assumed to be very distant from the pulsar, so that the upper limit of the integration was effectively infinity.

The termination of the pulsar wind will necessarily have the largest effect on its inverse Compton emission for small values of $`\theta `$, since regardless of the value of $`\eta `$, $`r_\mathrm{T}(\alpha )`$ is a monotonically increasing function. Furthermore, these effects will be most significant if $`\eta `$ is small, since in such cases the pulsar wind terminates close to the pulsar. The orbit of the PSR B1259$``$63 system is inclined at $`i=35^{}`$ to the plane of the sky, so an observer samples angles $`\theta `$ between $`90i=55^{}`$ and $`90+i=125^{}`$ over the binary period. The variation of $`\theta `$ over the orbit is shown in Figure 3. If $`\eta 1`$ the pulsar wind will terminate very close to the pulsar for small $`\theta `$, decreasing the observable emission from the unshocked wind most significantly at binary phases which correspond to values of $`\theta `$ that are close to the minimum, i.e. between around day $`200`$ and periastron.

### 3.1 Spectra

Figure 4 shows the spectra of inverse Compton emission for two values of $`\eta `$ for which the Be-star wind dominates the pulsar wind, which therefore terminates quite close to the pulsar for small $`\theta `$. The quantity plotted is the spectral energy distribution $`EF_E`$ where $`F_E`$ is the energy flux at Earth of inverse Compton scattered photons of energy $`E`$. The solid lines show the results when the termination of the pulsar wind is included, and the dashed lines show the spectra calculated by Ball & Kirk \[2000a\] when the shock termination is very distant from the pulsar. The spectra are shown assuming that the separation of the pulsar and the Be star is that which applies at periastron, and have been calculated for the extreme values of $`\theta `$ that can be sampled from Earth. Both of these extreme values occur very close to periastron, as can be seen from Figure 3, so the assumption of periastron separation is appropriate. Furthermore, it is assumed that the spectrum of target photons from the Be star can be approximated as a monochromatic distribution at a normalised energy $`ϵ_0=2.7k_\mathrm{B}T_{\mathrm{eff}}/(mc^2)10^5`$ where $`T_{\mathrm{eff}}`$ is the star’s effective temperature. Calculations including a more realistic Be-star spectrum have been presented by Ball & Kirk \[2000b\] and Kirk, Ball & Skjæraasen . The effects of including a better approximation to the target photon spectrum are generally small, but may warrant closer investigation if the system is detected in hard $`\gamma `$-rays.

All the spectra in Figure 4 are sharply peaked at an energy close to $`5\times 10^5\mathrm{MeV}`$, somewhat lower than the Thomson limit of $`\gamma _\mathrm{w}(0)^2m_\mathrm{e}c^2ϵ_05\times 10^6\mathrm{MeV}`$ because of the importance of Klein–Nishina effects. In the Klein–Nishina regime, i.e. when $`\gamma ϵ_0\stackrel{>}{}\mathrm{\hspace{0.17em}1}`$, a scattering particle loses a significant fraction of its energy to the scattered photon and the maximum scattered energy is then $`ϵ_{\mathrm{out}}^+\gamma `$. The upper limit of the emission spectrum is unaffected by the inclusion of the TS, since the scatterings that produce photons of this energy involve undecelerated electrons (with Lorentz factors $`\gamma _\mathrm{w}(0)`$) and occur very close to the pulsar. As the pulsar wind propagates away from the pulsar the wind electrons lose energy via the inverse Compton scattering process. Lower energy scattered photons are generated farther from the pulsar, and so tend to be suppressed by the termination of the wind.

For small $`\theta `$ and small $`\eta `$ the TS is very close to the pulsar where inverse Compton scattering is still very efficient because the target photon density is large and the scatterings are close to head on. The wind is terminated before it has been decelerated to Lorentz factors significantly below $`\gamma _\mathrm{w}(0)`$, so the scattering electrons are essentially monochromatic, as are the target photons in our approximation. This has the effect of dramatically reducing the emission at all energies below the upper cutoff, producing a scattered spectrum that is even more sharply peaked than that resulting from the unterminated wind. Such a case is illustrated in the top left panel for $`\eta =0.01`$ and $`\theta =55^{}`$, for which the scattered spectrum is essentially monochromatic at the upper cutoff.

The TS radius increases rapidly with increasing $`\alpha `$ (for a given $`\eta `$), with the rate of increase approaching infinity as $`\alpha `$ approaches $`180^{}\psi `$ where $`\psi `$ is the half-opening angle of the cone and $`\eta <1`$ is assumed. At energies near the upper cutoff the scattered spectrum from the terminated wind therefore rapidly approaches that from the unterminated wind as $`\theta `$ increases. The suppression at lower energies, which are produced by scatterings well away from the pulsar where the wind has been significantly decelerated, persists to larger $`\theta `$. This effect is dramatically illustrated by a comparison between each top panel of Figure 4 with the corresponding lower panel. For the case where $`\theta =125^{}`$ and $`\eta =0.01`$ the emission at and just below the cutoff energy is essentially unaffected by the inclusion of the TS, because the scatterings that produce photons of these energies occur closer to the pulsar than does the TS.

For larger values of $`\eta `$ the TS occurs further from the pulsar. It follows that for a given $`\theta `$, the difference between the inverse Compton emission from the terminated and the unterminated wind decreases as $`\eta `$ increases. For the case shown in the lower right panel of Figure 4, with $`\eta =0.1`$ and $`\theta =125^{}`$, the TS occurs so far from the pulsar that it has essentially no effect on the scattered emission in the energy range shown, though it does reduce the scattered flux at even lower energies.

### 3.2 Light curves

For $`\eta >1`$ the TS wraps around the Be star and the line of sight to the pulsar does not intersect the shock if $`\theta >\psi `$. Equation (4) implies that for $`\eta \stackrel{>}{}\mathrm{\hspace{0.17em}7}`$, $`\psi \stackrel{<}{}\mathrm{\hspace{0.17em}55}^{}`$ which is the minimum observable value of $`\theta `$ over the binary orbit of PSR B1259$``$63. Thus if $`\eta \stackrel{>}{}\mathrm{\hspace{0.17em}7}`$ the line of sight doesn’t intersect the TS at any phase of the binary orbit, the TS has no effect on the observable inverse Compton emission from the freely-expanding wind, and the results of Ball & Kirk \[2000a\] require no modification.

For $`\eta <1`$ the TS wraps around the pulsar and is intersected by the line of sight to the pulsar if $`\theta <180^{}\psi `$. It follows from Equation (4) that if $`\eta \stackrel{<}{}\mathrm{\hspace{0.17em}0.14}`$, $`180^{}\psi \stackrel{>}{}\mathrm{\hspace{0.17em}125}^{}`$ which is the maximum observable value of $`\theta `$, and thus the line of sight to the pulsar intersects the TS at all binary phases.

Figure 5 shows the orbital variation of the integrated flux density $`F_EdE`$, expected from PSR B1259$``$63 for $`\gamma _\mathrm{w}(0)=10^6`$. The dashed curve shows the emission from the unterminated wind as calculated by Ball & Kirk \[2000a\]. The solid line shows the emission when the wind is terminated by a shock at the position shown in Figure 1 for $`\eta =0.1`$. The TS is wrapped around the pulsar and is intersected by the line of sight to the pulsar at all binary phases. However, between days $`+50`$ and $`+400`$ the intersection point is sufficiently distant from the pulsar that it is beyond those radii where the majority of the inverse Compton scattering occurs. The termination of the wind therefore has little effect in reducing the inverse Compton emission in this period, and the light curve from the terminated wind is indistinguishable from that of an unconfined wind.

The effect of the termination of the wind is asymmetric because of the asymmetry in the dependence of the line-of-sight angle $`\theta `$ about periastron. The angle $`\theta `$ decreases very rapidly towards its minimum value just prior to periastron, which implies that the radius at which the line of sight to the pulsar intersects the TS does the same. The increase from the minimum value of $`\theta `$ to the maximum then occurs very rapidly, starting just a few days before periastron. The termination of the wind therefore decreases the emission from the freely expanding portion of the pulsar wind most effectively prior to periastron, and thus has the effect of decreasing the characteristic asymmetry in the $`\gamma `$-ray light curve about periastron. For $`\eta =0.1`$ Figure 5 shows that the effect of the TS is to reduce the maximum integrated inverse Compton flux density by a factor of $`4`$, to make the time at which the maximum occurs a few days later, and to reduce the asymmetry about periastron measured by the ratio of the integrated fluxes at days $`\pm 50`$, from $`3`$ to $`1.6`$.

Figure 6 shows the light curve for a shock-terminated wind with $`\eta =0.01`$. In this case the shock is so tightly wrapped around the pulsar that its intersection with the line of sight to the pulsar always occurs at a radius where inverse Compton scattering is still effective. The emission from the terminated wind is thus substantially less than that from an unconfined wind throughout the binary orbit. The observable flux in this case is reduced by a factor of between 1.5 and 3 over roughly 80% of the orbital period. The maximum integrated hard $`\gamma `$-ray flux density is a factor of $`10`$ lower than that from an unterminated wind and occurs $`6`$ days later, and the $`\pm 50`$ day asymmetry is actually reversed to $`0.95`$.

## 4 Discussion and conclusions

Ball & Kirk \[2000a\] argued that inverse Compton emission from the freely-expanding portion of the wind of PSR B1259$``$63 should be detectable at energies somewhere between GeV and TeV, depending on the wind parameters. The results presented in Section 3 suggest that even if the wind of the Be star dominates that of the pulsar, terminating it in a shock which wraps around the pulsar, the scattered emission should still be detectable for a wide range of wind parameters. However, the TS has the effect of decreasing the peak integrated $`\gamma `$-ray flux from the freely-expanding wind, and of decreasing the asymmetry which otherwise produces a characteristic light curve which shows a higher $`\gamma `$-ray flux before periastron than after.

Even if the Be-star wind dominates to such a degree that the unshocked pulsar wind is not detectable, the models of Kirk, Ball & Skjæraasen  indicate that inverse Compton emission from the shocked pulsar wind should be above the thresholds of current detectors.

The scattered $`\gamma `$-ray signal from this system is most likely to comprise a combination of emission from the shocked and unshocked regions of the pulsar wind. The reduced asymmetry and periastron to apastron ratio of the emission from the terminated wind may make it harder to deconvolve the two contributions, but only detailed modelling of real data will provide such answers.

The best opportunity for detecting such emission in the near future is afforded by the new CANGAROO II imaging Cherenkov telescope located in Australia \[Yoshikoshi et al. 1999\]. When the present lack of operational $`\gamma `$-ray telescopes is relieved by the launch of the GLAST and INTEGRAL observatories, inverse Compton emission from this pulsar system may well be detectable at lower energies in the MeV–GeV range.

The next periastron of PSR B1259$``$63 occurs in October 2000. CANGAROO II observations of the system will not be possible at times near periastron because at that time the object will be too close to the Sun. In any event, the light curves shown in Figures 5 and 6 are not realistic between about days $`25`$ and $`+25`$ because we have not modelled the interaction with the Be-star disk which is important at such epochs. Observations of the PSR B1259$``$63 system with CANGAROO II are planned for July and December 2000, when our models of the inverse Compton emission from the system should be applicable. If such observations are successful in detecting hard $`\gamma `$-ray emission from this unique pulsar system they should provide valuable insights into the properties of the pulsar wind. In particular, detection of inverse Compton $`\gamma `$-ray emission from the system should help to constrain values of the bulk Lorentz factor of the pulsar wind and the fraction of the spin-down luminosity carried by electrons and positrons compared to ions.

## References

Ball, L., & Kirk, J. G., 2000a, Astropart. Phys., 12, 335

Ball, L., & Kirk, J. G., 2000b, ASP Conf., 202, 531

Ball, L., Melatos, A. M., Johnston, S., Skjæraasen, O., 1999, ApJ, 541, L39

Eichler, D., Usov, V., 1993, ApJ, 402, 271

Girard, T., Willson, L. A., 1987, A& A, 183, 247

Giuliani, J. L., 1982, ApJ, 256, 634

Huang, R. Q., Weigert, A., 1982, A&A, 112, 281

Johnston, S., Manchester, R. N., Lyne, A. G., Bailes, M., Kaspi, V. M., Qiao, G., D’Amico, N., 1992, ApJ, 387, L37

Johnston, S., Manchester, R. N., Lyne, A. G., Nicastro, L., Spyromilio, J., 1994, MNRAS, 268, 430

Johnston, S., Manchester, R. N., Lyne, A. G., D’Amico, N., Bailes, M., Gaensler, B. M., Nicastro, L., 1996, MNRAS, 279, 1026

Johnston, S., Manchester, R. N., McConnell, D., Campbell-Wilson, D., 1999, MNRAS, 302, 277

Kennel, C. F., Coroniti, F. V., 1984, ApJ, 283, 710

Kirk, J. G., Ball, L., Skjæraasen, O., 1999, Astropart. Phys., 10, 31

Kirk, J. G., Ball, L., Skjæraasen, O., 2000, ASP Conf., 202, 527

Melatos, A., Johnston, S., Melrose, D. B., 1995, MNRAS, 285, 381

Tavani, M., Arons, J., 1997, ApJ, 477, 439

Yoshikoshi, T., et al., 1999, Astroparticle Phys., 11, 267


format: text

text:
# Current Status of the Numerical Simulations of d=3 SU(2) Lattice Gauge Theory in the Dual Formulation.

## 1 INTRODUCTION

We are reporting on our continuing investigations of $`SU(2)`$ lattice gauge theory in the dual formulation. A brief introduction to the formalism as well as to the techniques that have been developed for numerical simulations was presented at LATTICE’99 . As was stated there, the partition function of the conventional LGT can be converted, upon using the techniques of character expansion and group integrations, into the dual form  (for another perspective on this issue see ):

$$Z_d=\underset{\{j\}}{}(2j+1)C_{j_a}(\beta )\underset{i=1}{\overset{5}{}}\left\{\begin{array}{ccc}a_i& b_i& c_i\\ d_i& e_i& f_i\end{array}\right\}$$
(1)
This partition sum is defined over the dual lattice where $`\{j_a\}`$ live over the links and $`\{j_b\}`$ over the diagonals to the plaquettes. The convention for the diagonals is that they connect the vertices of the odd sublattice. We have collectively designated $`\{j_a\},\{j_b\}`$ by $`\{j\}`$. Each cube of the dual lattice is seen to be spanned by 5 tetrahedra of which one is spanned entirely by $`\{j_b\}`$ while four are spanned by three $`\{j_a\}`$ and three $`\{j_b\}`$. These features are illustrated in fig.1 which exhibits a two-cube cell on the dual lattice. The b-links are shown by dashed lines while the a-links are shown by solid lines. Each tetrahedron carries a weight factor which is the $`SU(2)`$ 6-j symbol $`\left\{\begin{array}{ccc}a& b& c\\ d& e& f\end{array}\right\}`$. Periodic b.c for the original lattice is crucial for this construction.

One immediate advantage of the dual approach is that variables are now discrete and for all practical purposes can be represented by short integers. Other advantages of this method have already been commented upon in .

### 1.1 Conventional Observables

Usually when one goes to the dual description, observables in the original formulation become very complicated in the dual picture and vice versa (for a discussion of the Wilson loop in the dual picture see ). It is therefore curious to see that, in the present context, a whole class of conventional observables are just as easily measured in both the formulations. In the conventional formulation, the expectation value of the normalized plaquette is given by

$$<P>=\frac{1}{N_p}\frac{\mathrm{log}Z(\beta )}{\beta }$$
(2)
where $`N_p`$ is the number of plaquettes in the lattice. In the dual picture it can be easily verified that this expectation value is given by

$$<P>=<\frac{d\mathrm{log}C_j(\beta )}{d\beta }>$$
(3)
where the average on the r.h.s is over the distribution of $`j`$-values on the a-links of the dual lattice. This formula can be generalized to the case of any function of the plaquette variable and more importantly to the two-point correlation functions of plaquettes or of functions of plaquette variables. This allows for a direct comparison between the dual and the conventional pictures.

## 2 UPDATES

In  we introduced two types of updating algorithms: local and quasi-local. It was already pointed out there that the local updates are not ergodic as they are unable to change integer spins to non-integer spins and vice versa. We illustrate the local update algorithm with the help of fig.1.

The b-link shown as a dotted line is contained in 6 tetrahedra (two of them containing only b-links); it is also contained in 6 triangles again two of which contain only b-links. For each triangle containing this link, one forms $`j_{max}^i=j_1^i+j_2^i`$ and $`j_{min}^i=|j_1^ij_2^i|`$ where $`j_{1,2}^i`$ are the remaining two spins in the $`i`$-th triangle containing the link to be updated. One then forms $`j_{max}=min_i\{j_{max}^i\}`$ and $`j_{min}=max_i\{j_{min}^i\}`$ and updates the link to a value in the range $`(j_{min},j_{max})`$. Clearly this does not change half-integers to integers and vice versa and the move is not ergodic by itself.

In addition to the problem of ergodicity, there are two additional problems with this method. Every now and then one may face a situation where $`j_{min}=j_{max}`$ in which case the move can not change the value of the link. Our code keeps track of this possibility through a counter irange and whenever this happens one skips the particular site and moves on to the next one. The other problem has to do with the fact that 6j’s occasionally vanish even when all the triangle relations hold. Very rarely it so happens that every new value of the link is such that the product of 6j’s over the tetrahedra (four in the case of local updates of a-links and six in the case of b-links) vanishes and the new configuration is to be eliminated because it does not contribute to either the partition sum or to any of the expectation values. This is kept track of by the counter ispecial.

### 2.1 Quasi-local Updates

As already discussed in , if we were to change a link from half-integer (henceforth referred to as fermionic) to integer (bosonic) or vice versa, in general the change proliferates throughout the lattice as at least two links have to be $`Z_2`$-flipped in every triangle. What one needs are flippings involving the smallest cluster of links.Two categories of such quasi-local updates were identified in . At an even site of the lattice, the update algorithm called globeven involves the interior links of a cluster of 8 tetrahedra as shown in fig.2.

Likewise, at odd sites the corresponding update is called globodd and it involves the 6 a-links and 12 b-links in the interior of the 8-cube cluster involving 32 tetrahedra as shown in fig.3.

In both these schemes of quasi-local updates, one keeps all the links on the outer surface fixed and $`Z_2`$-flips all the interior links. The changes should be so made as to satisfy all the triangle inequalities. In , the basic idea of achieving this goal was briefly mentioned. This is done using variables originally introduced by Bargman and developed by . This will be discussed in detail in . While this method based on the K(agome)-variables is the systematic approach to the quasi-local updates, it is pretty involved and time consuming. A simplified algorithm was proposed in  which consisted of simply shifting all spin values on the interior links by $`\frac{1}{2}`$ which can be seen to preserve all the triangle inequalities. But this can result in substantial changes in the Boltzman weight resulting in abrupt changes. In the current investigations we modified this by following the fixed $`Z_2`$-flip by randomly selecting a set of configurations satisfying the triangle inequalities and then performing heatbath or metropolis in the subspace of these configurations.

## 3 THE SIGN PROBLEM

It was already emphasised in  that one of the main drawbacks of the dual formulation is the lack of positivity of the Boltzman weight which arises out of the 6j’s not being of a particular sign. There we had proposed to use the trick of generating important samples with probabilities that are the modulus of the weights occurring in eqn (1) but by absorbing the sign of the configuration into the observables. It was soon pointed out to us by Wosiek  that no matter how small the local probability of a wrong sign, in the thermodynamic limit there will be an almost equal fraction of configurations with the two signs. This meant that the average value of the sign is very close to zero making the procedure questionable. One way out is to regroup the terms in eqn(1) so that for each group the weight is positive. But it is not clear how to do this. To make some progress we decided to adopt another strategy based on the fact that sign changes in the configurations come about through completely local sources. To see this, let us look at the expression for the average of a generic observable O:

$$O_p=\frac{Oϵ_{\stackrel{~}{p}}}{ϵ_{\stackrel{~}{p}}}$$
(4)
where $`\{\stackrel{~}{p}_i=\frac{|p_i|}{P}\}`$ with $`P=_i|p_i|`$ and $`p_i`$ are the weights occurring in eqn(1) and are not necessarily positive. Further, $`ϵ_i`$ is the sign of $`p_i`$. Now if it so happens that when thermal equilibrium is reached the average value of $`O`$ (in the usual sense) is insensitive to the sign of the configuration in the sense that the average over positive sign configurations is very close to the average over negative sign configurations, the numerator in eqn(4) is also proportional to $`ϵ_{\stackrel{~}{p}}`$ and we can simply ignore the sign problem.

We wish to argue that this is very likely to happen for relatively local observables and also for the correlation functions of relatively local observables. The reason for this is the fact that the source for the change of sign of the configuration is ultralocal. Combining this with the expectation that both signs become equally likely as the system size increases, one is led to surmise that the configurational average of reasonably local observables as well as the correlation functions of reasonably local observables are insensitive to the sign of the configuration. So far most of our runs have borne out this expectation. We found that whenever there was a reasonable fraction of configurations with each sign, the averages of inner and outer tetrahedra were indeed insensitive to the sign in the limit of thermalisation (to be more precise, in the limit of nearly stationary expectation values). We are investigating this with other classes of observables. We also propose to make a systematic study of the sign problem by studying the way the signs fluctuate locally as well by studying their correlations over the lattice.

## 4 METROPOLIS AND HEATBATH

For the local updates both the metropolis and heatbath methods work in fairly standard ways because the allowed range of values for the new link are very small. But for the quasi-local moves the situation is quite different. Even in the case of the quasi-local updates at even sites, we are dealing with the space of allowed spin values on six links. Truncating the maximum value of spins at $`\frac{9}{2}`$ as we did, we have $`10^6`$ values to check for triangle inequalities and subsequently picking a new element from among the list. This is a hopeless task. Even using the K-variable method is time consuming. So what we did was to apply the fixed $`Z_2`$-flip and select the first a-link randomly from its allowed possibilities, the second a-link randomly from its allowed values (which depends on the random choice made for the first a-link) and so on. For metropolis algorithm, after selecting a sextuplet of values, the new set was accepted or rejected according to the usual metropolis rule. For the heatbath method we likewise selected a certain number of sextuplets before the $`Z_2`$-flip and an equal number after the flip. The heatbath algorithm was then applied to this manageable subset. At present we do not have a more rigorous argument to justify our heatbath method. For the globodd update the problems are even more severe and we have again adopted the same pragmatic approach.

### 4.1 INITIAL CONFIGURATIONS

As initial configurations we chose from among three classes: i) cold start - here we set all the a-links and b-links to $`j=1`$. This configuration clearly satisfies all the triangle inequalities. ii) extreme cold start - here all a-links and b-links were set equal to $`j=0`$. It is worth noting that the extreme cold start configuration is a fixed point of the local update algorithms. iii) hot start - here a-links and b-links are randomly chosen in a way that all triangle inequalities are satisfied. A practical way of reaching such configurations is to first start with a cold start configuration and go on choosing links randomly in their allowed ranges. After sweeping the lattice sufficiently many times one gets to a hot start configuration.

## 5 RESULTS

We carried out our simulations on a $`8^3`$ lattice at various $`\beta `$-values. We simulated with the absolute values of the weights given by eqn(1) and monitored the fluctuations in the sign of the configurations. In a few runs we noticed that the sign changed very few times. In those cases, the average values of the inner (outer) tetrahedra for the positive sign and negative sign configurations did not approach the same values even after many iterations. But for these cases $`ϵ_{\stackrel{~}{p}}`$ is not close to zero and one can safely use eqn (4). On the other hand for a majority of runs, the sign fluctuated quite a bit and even when $`ϵ_{\stackrel{~}{p}}0.5`$, the average values for the two cases were very close to each other bearing out the surmise made earlier. This was so for all the classes of update algorithms, initial states and metropolis/heatbath. We plan to study this issue for a larger class of observables.

In our simulations using only the local updates, we found that while the results of metropolis and heatbath converged for given initial states, the results did not converge for trials with different initial states. We have interpreted this as due to the non-ergodic nature of the local updates. The results for the quasi-local updates were not so clear cut. We also ran the code with the updating algorithms being chosen randomly; at the first instance we have run the code such that all updating algorithms are equally likely. Here we find that while the runs for the cold-start and hot-start converged, the results for metropolis and heat bath are widely divergent. We have not really understood the reason for this.

## 6 CONCLUSION

Clearly a lot more work needs to be done before establishing the dual framework as a viable alternative for numerical simulations of Gauge and Spin systems. The most pressing problem is to gain a satisfactory understanding and control of the sign problem. While our surmise in this direction is well borne out by the initial runs, we are yet to show its correctness in a controlled manner. One interesting possibility is to study the model defined by the absolute value of weights in its own right and we plan to undertake that study shortly. But we still have to understand the lack of convergence between the different runs before we can take any of the conclusions seriously. To this end, we are going to investigate the efficacy of our metropolis and heat bath methods by studying them in contexts where we can get results by deterministic as opposed to stochastic methods.

Another line of attack we are contemplating is to look for simple models that are analytically tractable but which nevertheless have all the essential ingredients of the dual formulation studied here. Indeed, such models were encountered in the study of analytic variational methods for euclidean lattice gauge theories  and in the analytic evaluation of the partition function of the unit hypercube for $`d=4`$ $`SU(2)`$ lattice gauge theories . We are hoping that the insights so gained will enable us to complete the investigations reported here.

There are a large class of very important problems that belong to the same calculational scheme as what we have attempted here. The most notable are: the $`d=4`$ version of these results as obtained by Halliday and Suranyi, the class of models considered by Ooguri  and their non-topological descendants etc.Language and techniques very similar to the ones used in the dual picture are also becoming increasingly popular in treatments of quantum gravity . Replacing the 6j’s by q-deformed 6j’s would immediately relate these developments to the well known Turaev-Viro models and their generalizations. Clearly it is of great importance to make further progress in the numerical simulations of theories in the dual formulation.

## 7 ACKNOWLEDGEMENTS

NDH wishes to thank Dick Haymaker, J. Wosiek, Erhard Seiler, Peter Weisz, Srinath Cheluvaraja, A. Hart, N. Kawamoto and Maarten Golterman for many interesting discussions. We are thankful to all those who visited our posters at LATTICE2000 for the lively discussions.


text:
# 1 Parton momentum distribution, summed over quarks and gluons. Our simulation is compared with the GRV parameterization []. (a) At the small scale 𝑄₀=0.48 GeV. (b) Evolved to 𝑄²=10 GeV2, with a cut so that branchings do not produce partons below 𝑥_{𝑚𝑖𝑛}=10⁻⁴.

LU TP 00-47

hep-ph/0011282

November 2000

A Toy Model of Colour Screening
in the Proton
Johann Dischler and Torbjörn Sjöstrand
Department of Theoretical Physics,
Lund University, Lund, Sweden

Abstract
In hadronic collisions, the mini-jet cross section is formally divergent in the limit $`p_{}0`$. We argue that this divergence is tamed by some effective colour correlation length scale of the hadron. A toy model of the hadronic structure is introduced, that allows an estimate of the screening effects, and especially their energy dependence.

The perturbative parton–parton interaction cross section in a hadronic collision is divergent roughly like

$$\frac{\mathrm{d}\widehat{\sigma }}{\mathrm{d}p_{}^2}\frac{\alpha _\mathrm{s}^2(p_{}^2)}{p_{}^4}$$
(1)
for $`p_{}0`$. This is to be convoluted with parton distributions $`f_i(x,Q^2p_{}^2)`$ to give a hadronic interaction cross section. The dominant contributions come from scattering by $`t`$-channel gluon exchange. With an artificial lower cut-off scale $`p_{\mathrm{min}}`$,

$$\sigma _{\mathrm{int}}(p_{\mathrm{min}})=\underset{i,j,k,l}{}dx_1f_i(x_1,Q^2)dx_2f_j(x_2,Q^2)_{p_{\mathrm{min}}^2}dp_{}^2\frac{\mathrm{d}\widehat{\sigma }_{ijkl}(\widehat{s}=x_1x_2s)}{\mathrm{d}p_{}^2}.$$
(2)
The jet/mini-jet cross section is twice this, since each scattering gives two jets in the lowest-order approximation considered here. Studying e.g. Tevatron collider energies, $`\sqrt{s}2`$ TeV, $`\sigma _{\mathrm{int}}(p_{\mathrm{min}})`$ exceeds the total $`\mathrm{p}\overline{\mathrm{p}}`$ cross section $`\sigma _{\mathrm{tot}}`$ for $`p_{\mathrm{min}}\text{ }\stackrel{<}{}\text{ }3`$ GeV. Since 3 GeV still is well above the $`\mathrm{\Lambda }_{\mathrm{QCD}}`$ scale of $`0.2`$ GeV, there is no obvious reason why perturbation theory should have broken down, and one would seem to be in trouble. The resolution of this paradox probably comes in several steps, and at least in two.

First of all, the mini-jet cross section above is inclusive. Thus, if an event contains two parton–parton interactions , it counts twice in $`\sigma _{\mathrm{int}}(p_{\mathrm{min}})`$ but only once in $`\sigma _{\mathrm{tot}}`$. Thereby $`\sigma _{\mathrm{int}}(p_{\mathrm{min}})>\sigma _{\mathrm{tot}}`$ becomes allowed. Multiple parton–parton interactions is the concept that, based on the composite nature of hadrons, indeed several parton pairs may scatter in a typical hadron–hadron collision . Over the years, evidence for this mechanism has accumulated , such as the direct observation by CDF . The events studied experimentally, with two parton pairs at reasonably large $`p_{}`$, only form the tip of the iceberg, however. One may expect that most interactions are at lower $`p_{}`$, where they do not produce visible jets, but only contribute to the underlying event structure. As such, they are then believed to be at the origin of a number of key features, like the broad multiplicity distributions, the significant forward–backward multiplicity correlations, and the pedestal effect under jets .

While the effects of multiple interactions on event properties are smaller for lower-$`p_{}`$ scatterings, in the models we studied so far the interaction rate increases faster with decreasing $`p_{}`$ than the effect per scattering, so there is no stability in the limit $`p_{\mathrm{min}}0`$.

The second necessary aspect is likely that a regularization of the jet cross section should occur at small $`p_{}`$ from the fact that the incoming hadrons are colour singlets — unlike the coloured partons assumed in the divergent perturbative calculations — and that therefore the colour charges should screen each other in the $`p_{}0`$ limit. Thus $`p_{\mathrm{min}}`$ could take on a physics meaning, roughly, as the inverse of some colour screening length in the hadron. Of course, one would not expect a sharp cut-off of the mini-jet cross section at $`p_{\mathrm{min}}`$, but rather a smooth dampening. Nevertheless the $`p_{\mathrm{min}}`$ parameter provides a useful first approximation. Fits to data typically give $`p_{\mathrm{min}}2`$ GeV.

One key issue is the energy dependence of $`p_{\mathrm{min}}`$; this may be relevant e.g. for comparisons of jet rates at different Tevatron energies, and even more for any extrapolation to LHC energies. The question actually is more pressing now than at the time of the study in , since nowadays parton distributions are known to be rising more steeply at small $`x`$ than the flat $`xf(x)`$ behaviour normally assumed for small $`Q^2`$ before HERA. This translates into a more dramatic energy dependence of the multiple-interactions rate for a fixed $`p_{\mathrm{min}}`$, and thence to a charged multiplicity rising faster than data . Based on such considerations, the $`p_{\mathrm{min}}`$ in the Pythia program  was made explicitly energy-dependent some time ago:

$$p_{\mathrm{min}}=(1.9\mathrm{GeV})\left(\frac{s}{1\mathrm{TeV}^2}\right)^ϵ,$$
(3)
with $`ϵ=0.08`$. This value is picked to agree with the $`s^ϵ`$ behaviour assumed in the parameterization of the total cross section in hadron–hadron collisions  that, via reggeon phenomenology, should relate to the behaviour of parton distributions at small $`x`$ and $`Q^2`$. A study of the energy dependence of $`\mathrm{d}N_{\mathrm{charged}}/\mathrm{d}\eta |_{\eta =0}`$ also confirms that an $`ϵ`$ in this range gives sensible agreement with data .

In the following, we will develop an extremely simple model for the proton structure, where the amount of screening may be studied . Our basic approach is to assign an exclusive parton configuration to each proton. This configuration is — for a given evolution scale $`Q`$ — defined in terms of the longitudinal momentum fraction $`x`$ values, the transverse coordinate positions and colour charges of ‘all’ the partons present in a proton. The traditional inclusive parton-distribution picture is obtained as the average over many such proton ‘snapshots’. A simple plane wave represents the exchange of a gluon between two colliding hadrons. Longitudinal coordinates of the partons are not specified, so the exchanged gluon is only characterized by its transverse momentum. The exclusive proton picture then allows a comparison of a coherent versus an incoherent sum of colour charges coupling to the plane wave. It is this ratio, coherent over incoherent, that will provide our measure of screening between partons in the proton. In particular, the long wavelength limit then explicitly corresponds to a vanishing of interactions between two colourless hadrons. As the collision energy is increased, more partons at smaller $`x`$ values become accessible and the screening effects increase in importance.

We certainly know this ansatz to be wrong in many of its details, but still hope that it will catch enough of the spirit to provide some insight. It is clear that a full-fledged description — e.g. with lattice QCD calculations — is way beyond the current capability, so toy models is all that could be offered today. The approach we present here should be viewed as only one possible aspect or formulation of the required dampening effect. The vanishing net colour charge of hadrons has been used in models since long , as has the spatial size of the colour singlet wave function . Our model also has some similarities with the dense-packing view of the proton , although without the possibility of parton recombinations. Among alternatives, one could mention the possibility of a nontrivial structure of the QCD vacuum , the solution of a Dyson-Schwinger equation for the gluon propagator , the introduction of an effective gluon mass from lattice QCD results , and the possibility of several hard scatterings within a single parton chain , e.g. formulated in terms of non-integrated structure functions and the Linked Dipole Model description of non-$`p_{}`$-ordered chains .

In our model, the $`Q^2`$ evolution of the parton distributions is given by the conventional DGLAP equations . As usual, the starting configuration requires nonperturbative input. A valence-like ansatz at a small $`Q_0`$ scale  limits the number of required free parameters. We have therefore used the GRV distributions as a reference , where the three valence quarks together with two ‘valence gluons’ almost completely define the proton at the small scale $`Q_0=0.48`$ GeV. In our simulation, these five partons are first selected according to a single-particle distribution

$$f(x)x^\alpha (1x)^\beta ,$$
(4)
with tuned parameters $`\alpha =0.4`$ and $`\beta =1.2`$. Thereafter all five $`x_i`$ are normalized so that their sum corresponds to unity:

$$(x_i)_{\mathrm{norm}}=\frac{x_i}{_{j=1}^5x_j}.$$
(5)
This means that the original distribution is pushed towards the middle, and it is especially difficult to tune the low-$`x`$ tail, at least with the simple $`f(x)`$ ansatz above. A reduction of the $`Q_0`$ scale from 0.48 to 0.44 GeV helps improve the small-$`x`$ behaviour, Fig. 1a.

From there on, the evolution follows the standard leading-order DGLAP formalism, expressed in terms of Monte Carlo-generated branchings of partons, $`\mathrm{q}\mathrm{qg}`$, $`\mathrm{g}\mathrm{gg}`$ and $`\mathrm{g}\mathrm{q}\overline{\mathrm{q}}`$. Unfettered, such an approach gives an infinity of gluons in the $`x0`$ limit. For the exclusive approach we have in mind, this is not so physical. In the study of scatterings above some $`p_{\mathrm{min}}`$ scale, only partons above some related $`x_{\mathrm{min}}`$ scale can at all interact. Partons below $`x_{\mathrm{min}}`$ then have a dubious existence. To simplify the picture, we assume they are not resolved, i.e. that only branchings producing both daughters above $`x_{\mathrm{min}}`$ are allowed. This reduces the amount of evolution, but not as dramatically as might be guessed at first glance. In Fig. 1b we see that, in a typical case, the parton density with an $`x_{\mathrm{min}}`$ constraint is at most 20% above the no-constraints normal evolution in the region $`x>x_{\mathrm{min}}`$, while the former obviously vanishes for $`x<x_{\mathrm{min}}`$. Disagreements are more visible — but also less interesting for us — in individual parton species, especially for sea quarks, which are not part of our simple ansatz at $`Q_0`$ but are in the GRV one.

The kinematics of a hard scattering requires that $`x_1x_2=\widehat{s}/s4p_{}^2/s`$. For collisions at central rapidities this suggests $`x_{\mathrm{min}}2p_{}/E_{\mathrm{CM}}`$. Away from the center, $`x_{\mathrm{min}}`$ would be smaller on one side, down to $`4p_{}^2/E_{\mathrm{CM}}^2`$, and correspondingly larger on the other, which to some extent should cancel out. Therefore we will have as standard scenario that $`x_{\mathrm{min}}`$ should increase proportionately to the $`p_{}`$ considered, but will include as alternative a fixed $`x_{\mathrm{min}}`$. Furthermore, we will make the association that $`E_{\mathrm{CM}}1/x_{\mathrm{min}}`$, for $`x_{\mathrm{min}}`$ evaluated at some fixed reference $`p_{}`$ like 1 GeV, but allow as alternative that the relation could be more like $`E_{\mathrm{CM}}1/\sqrt{x_{\mathrm{min}}}`$.

Next we consider the assignment of transverse position coordinates. The original five partons are here selected according to a Gaussian shape, with a radius of 0.7 fm in each of the two transverse dimensions. This gives a minimal correlation between the position coordinates of these five partons, and thus probably errs on the side of simplicity. Attempts with somewhat more correlated forms gave fairly similar results, however, so the Gaussian ansatz is not so critical.

The spatial extent of a fluctuation like $`\mathrm{q}\mathrm{qg}\mathrm{q}`$, $`\mathrm{g}\mathrm{gg}\mathrm{g}`$ and $`\mathrm{g}\mathrm{q}\overline{\mathrm{q}}\mathrm{g}`$ is of order $`1/Q`$ according to the uncertainty relation. Therefore, when a DGLAP branching occurs at a scale $`Q`$, the two daughters are assumed to have time to fluctuate a distance of this order away from their production coordinates before they branch in their turn, or are probed by our plane gluon wave. The daughter partons are assumed to move out in opposite directions from the production vertex, with a uniform azimuthal distribution. The two distances are picked independently of each other, uniformly between 0 and $`1/Q`$, alternatively between 0 and $`2/Q`$. We will see that this choice is relevant, whereas results are essentially equivalent e.g. for exponentially dampened distributions with the same mean. For a parton which takes the momentum fraction $`z`$ in a branching, an additional factor $`\sqrt{1z}`$ is included, giving an approximate $`p_{}/Q`$ dampening factor for the separation in transverse rather than longitudinal direction. Here the relation $`p_{}=\sqrt{1z}Q`$ follows in light-cone kinematics assuming the mother and recoiling partons to be massless. It is of some relevance that a very soft gluon emission does not affect the position of the hard parton unduly, since else multiple soft gluons could lead to a too rapid increase of the proton radius with $`Q`$. For the final result here, however, the $`\sqrt{1z}`$ factor turns out not to be crucial.

Finally we come to the colour space picture. Here we have picked a simple planar representation, with the three primary colours placed in a triangle around the origin. Thus the full phase information is lost; i.e. it is possible to ensure that the proton state is colour neutral but not that it is in a singlet. This certainly is a major simplification of the real world, and one which it is difficult to estimate the impact of. A gluon emission corresponds to a quark changing colour, e.g. $`\mathrm{q}(r)\mathrm{q}(b)+\mathrm{g}(r\overline{b})`$. The emitted gluon thus carries one colour and a different anticolour, i.e. the two colour diagonal gluon states are not populated; another simplification. Further, the ratio of gluon to quark charge is $`\sqrt{3}`$ versus $`\sqrt{N_C/C_F}=3/2`$ in QCD. Also in a gluon branching to two, a new colour is picked, e.g. $`\mathrm{g}(r\overline{b})\mathrm{g}(r\overline{g})+\mathrm{g}(g\overline{b})`$, in order to avoid diagonal gluons. The gluon branching to quarks is unambiguous, $`\mathrm{g}(r\overline{b})\mathrm{q}(r)+\overline{\mathrm{q}}(\overline{b})`$. For colour assignments in the starting configuration, the three valence quarks are first picked to be $`r+g+b`$, and then the two gluons are handled as if emitted from two of the quarks, picked at random.

An effort is thus made to retain overall vanishing colour for all partons in a proton that are seen by an exchanged gluon. In case of multiple interactions, it could be argued that only the ‘first’ exchange would be between singlet hadrons, while ‘subsequent’ ones would be between already coloured hadrons. Clearly, the imposition of such a time ordering can be questioned but, even if taken at face value, colour screening could still provide a significant dampening of the naive perturbative answer.

To summarize, an explicit parton configuration of the proton is constructed for any value of the evolution parameter $`Q`$ above some low $`Q_0`$ starting scale. Each parton $`k`$ is characterized by its momentum fraction $`x_k`$, its transverse position coordinate $`𝐫_k`$ and its two-dimensional colour charge $`q_k`$. If such a proton is probed by a gluon plane wave function of transverse momentum $`𝐩_{}`$, we may define a ratio

$$A=\frac{|_kq_ke^{i𝐫_k𝐩_{}}|^2}{_k|q_k|^2}$$
(6)
between a coherent and an incoherent sum of the colour charges. Note that we assume the exchanged gluon to be a superposition of all possible colour-anticolour pairs, so that we do not need to single out a special direction in colour space.

The screening effect is most easily observed if one fixes the $`Q`$ of the probed proton and then varies the $`p_{}`$ of the probing gluon, Fig. 2a. The suppression factor $`A`$ then has to vanish quadratically as $`p_{}0`$, but in the interesting $`p_{}`$ range the dominant aspect is instead the slower rise to the asymptotic limit $`A=1`$. The choice of a fixed $`Q`$ scale is representative for the way parton distributions de facto are used with $`Q_0`$ as scale for $`Q<Q_0`$. For popular distributions like the current CTEQ and MRST ones , where $`Q_0^21`$ GeV2, the inclusion of a suppression factor then significantly modifies the picture for $`Q<Q_0`$. Actually $`Q_0`$ has tended to come down with time; some years ago $`Q_0^24`$ GeV2 was the norm. Also remember that the cross section suppression is $`A^2`$, i.e. with one factor of $`A`$ for each hadron beam.

A conventional choice, but not a unique one, is to pick the $`Q`$ scale of parton distributions to agree with the $`p_{}`$ scale of the hard interaction. Then the approach $`A1`$ is much slower at large $`p_{}`$, since the branchings that occurred at scales $`Q`$ not so much smaller than $`p_{}`$ are not fully resolved, Fig. 2b. We remind that the proton seen by a gluon probe is a much more ‘busy’ place than the one seen by a photon in Deeply Inelastic Scattering processes, where the gluon content is not probed directly. Thus the destructive interference we introduce here does not invalidate the conventional structure function measurements and partonic interpretations. In hadronic physics, higher order perturbative corrections could partly contain, and partly mask, the physics of an $`A`$ somewhat below unity.

More interesting is the behaviour at small $`p_{}`$. We here note a spike in $`A`$ at the scale $`p_{}=Q_0=0.44`$ GeV, Fig. 2b. At this scale, the proton consists of the original five partons, and these are normally well separated in space. Above it, the partons emitted in the evolution process tend to have a more clustered spatial distribution and therefore screen more. Below it, there is no further evolution of the proton but only a straightforward destructive interference. Thus the spike is an artifact of how we discontinuously in $`Q`$ go from a fixed proton picture to a very rapid evolution; remember that the evolution rate behaves like $`\alpha _\mathrm{s}(Q^2)\mathrm{d}Q^2/Q^2`$ and thus is as largest just at the lower cut-off. In order to provide a smoother and physically maybe more sensible picture, we therefore do not probe the proton at scales below 1 GeV, i.e. put $`Q=\mathrm{max}(p_{},1\mathrm{GeV})`$. The original $`Q_0`$ scale and evolution below 1 GeV still lives on in terms of a sensible $`x`$/$`𝐫`$/colour configuration at 1 GeV.

Fig. 2b is based on a fixed $`x_{\mathrm{min}}`$ scale. However, an $`x_{\mathrm{min}}`$ scaling $`p_{}`$ behaves very similarly, Fig. 2c, if the two are matched for $`p_{}=1`$ GeV. If instead $`Q`$ is fixed, the shape is more changed.

We now turn to the prime objective of this paper, namely to study the energy dependence of $`p_{\mathrm{min}}`$. Actually, the sharp cut-off $`\theta (p_{}p_{\mathrm{min}})`$ is not very physical. A more likely behaviour is a dampening factor something like $`p_{}^2/(p_{}^2+p_0^2)`$, where $`p_0`$ takes over the role of free parameter. With one such factor per incoming proton, the $`1/p_{}^4`$ singularity of the perturbative parton-parton cross section is regularized.

Whereas the proposed dampening form does not completely reproduce the curves shown above, qualitatively there is agreement. Instead of doing a fit to the curve shape, $`p_0`$ is extracted as the $`p_{}`$ value for which $`A=1/2`$. This would also have been a sensible definition of a $`p_{\mathrm{min}}`$, backed up by experience with the impact of the two regularization procedures on fits to data, where $`p_0p_{\mathrm{min}}`$ is obtained .

The resulting $`p_0`$ values are plotted as a function of $`1/x_{\mathrm{min}}`$ in Fig. 2d. Three scenarios are compared:

* The proton probed at a fixed $`Q=2`$ GeV scale, travel distance of a parton is distributed between 0 and $`1/Q`$ of the branching where it is produced.
* Also proton probed at a fixed $`Q=2`$ GeV scale, but travel distance twice as large.
* The proton is probed at a running $`Q=p_{}`$, travel distance between 0 and $`2/Q`$ as in 2.

In all these three cases, representative for several more studied by us, a clear energy-dependence of $`p_0`$ is observed. If one tries to make a fit to an exponential form $`p_01/x_{\mathrm{min}}^\delta `$, typically one obtains a $`\delta `$ in the range 0.05 to 0.08. However, two comments are in order. One is that the curves tend to rise less steeply at higher energies than implied by the fit form. The other is that, if one attempts the association $`1/x_{\mathrm{min}}E_{\mathrm{CM}}`$, then eq. (3) would have predicted $`\delta =2ϵ=0.16`$, i.e. a much steeper slope. The alternative association $`1/x_{\mathrm{min}}E_{\mathrm{CM}}^2`$ would give nicer agreement, but we iterate that it may be the less plausible ansatz.

In summary, we have presented a very simpleminded model of the longitudinal momentum, transverse coordinates and colour structure of the proton, in order to be able to estimate colour screening effects at small $`p_{}`$. Many of the details turn out to be less critical, but there are two key assumptions. The first is that the proton, by virtue of being a colour singlet, should couple to a soft gluon with reduced strength. The second is that, at higher energies, the proton can effectively be resolved into more partons that, in principle, can interact to give scatterings at a given $`p_{}`$ scale. We show that these two assumptions lead to the conclusion that the perturbative ansatz for the cross section is only valid above some $`p_{\mathrm{min}}`$ scale that increases with energy. If this increase is fitted to a form $`p_{\mathrm{min}}E_{\mathrm{CM}}^\kappa `$, any value $`\kappa `$ in the range 0.05 to 0.16 could be obtained, which certainly overlaps with an experimentally acceptable behaviour. However, one lesson is that the ansatz in eq. (3) may tend to overestimate the rate of increase at large energies, and thereby e.g. lead to an underestimate of multiplicities at LHC energies. Conversely, the study of minimum bias events at LHC may help improve the understanding of very low-$`p_{}`$ physics processes.


text:
# PSCz vs. 1.2 Jy velocity fields: a Spherical Harmonics Comparison

## 1 Introduction

Nusser & Davis (NS) show that in linear gravitational instability (GI) theory the peculiar velocity field in redshift space is irrotational and thus can be expressed in terms of a potential: $`\stackrel{}{v}=\mathrm{\Phi }(\stackrel{}{s})`$. The angular dependencies of the potential velocity field and the galaxy overdensity field \[both measured in redshift space and expanded in spherical harmonics, $`\mathrm{\Phi }_{lm}(s)`$ and $`\widehat{\delta }_{lm}(s)`$, respectively\] are related by a modified Poisson equation:

$$\frac{1}{s^2}\frac{d}{ds}\left(s^2\frac{d\mathrm{\Phi }_{lm}}{ds}\right)\frac{1}{1+\beta }\frac{l(l+1)\mathrm{\Phi }_{lm}}{s^2}=\frac{\beta }{1+\beta }\left(\widehat{\delta }_{lm}\frac{1}{s}\frac{d\mathrm{ln}\varphi }{d\mathrm{ln}s}\frac{d\mathrm{\Phi }_{lm}}{ds}\right),$$
(1)
where $`\varphi (s)`$ is the selection function and $`\beta \mathrm{\Omega }^{0.6}/b`$, where $`\mathrm{\Omega }`$ is the matter density and $`b`$ the bias of the galaxy distribution. To solve this differential equation we first compute the density field on an angular grid using cells of equal solid angle and 52 bins in redshift out to $`s=\mathrm{18\hspace{0.17em}000}`$ kms-1.

$$1+\widehat{\delta }_j(\stackrel{}{s}_n)=\frac{1}{(2\pi )^{3/2}\sigma _{1.2n}^3}\underset{i}{\overset{N_j}{}}\frac{1}{\varphi (s_i)}\mathrm{exp}\left[\frac{(\stackrel{}{s}_n\stackrel{}{s}_i)^2}{2\sigma _{1.2n}^2}\right]$$
(2)
where the sum is over all the galaxies within the catalogue $`j`$, $`N_j`$. The Gaussian smoothing width for the cell $`n`$ at redshift $`s_n`$, $`\sigma _{1.2n}`$, is given by $`\sigma _{1.2n}=\text{max}\{100,[\overline{n}_{1.2}\varphi _{1.2}(s_n)]^{1/3}\}`$ km s$`^1`$, where $`\overline{n}_{1.2}`$ and $`\varphi _{1.2}`$ are the 1.2 Jy mean number density and selection function, respectively.

## 2 Datasets

The IRAS PSCz catalogue has been recently completed and contains 15 500 IRAS PSC galaxies with a 60 $`\mu `$m flux larger than 0.6 Jy. The average depth of this survey is $``$ 100 $`h^1`$Mpc . In our analysis we will restrict to the PSCz sub-sample of 11 206 galaxies within 20 000 km s-1 from the Local Group. The unobserved region is modeled by the angular mask of Saunders et al. , which leaves unmasked 80$`\%`$ of the sky. This region is filled-in with the cloning procedure described in Branchini et al. (B99).

The 1.2-Jy catalogue (Fisher et al. ) contains 5 321 IRAS PSC galaxies with a 60 $`\mu m`$ flux limit of 1.2-Jy within 20 000 km s-1 of the Local Group. This catalogue has a slightly larger sky coverage of $``$ 87.6$`\%`$ and a smaller median distance of $``$ 84 $`h^1\mathrm{Mpc}`$. This catalogue is supplemented with “synthetic” objects in the ZoA and other excluded regions following the same technique.

## 3 Results, Discussion and Conclusions

We compare the line of sight peculiar velocities relative to the Local Group, $`u(\stackrel{}{s})[\stackrel{}{v}(\stackrel{}{v})\stackrel{}{v}_{LG}]\widehat{s}`$, obtained from the IRAS PSCz and 1.2-Jy redshift surveys gravity fields. The comparison is performed in redshift-space. To perform the decomposition in spherical harmonics coefficients we apply $`\beta =0.6`$.

In Fig. 1 we display the monopole of the velocity field, $`u_{00}`$. In the left panel, the estimate of the velocity monopole in the 1.2-Jy (dashed line) is systematically larger than the PSCz one (continuous line) in the range $`\mathrm{1\hspace{0.17em}000}<s<\mathrm{8\hspace{0.17em}000}`$ km s-1.

In Fig. 2 we show the three velocity dipole components of the two surveys along with the total amplitude (bottom left panel). The various dipole components exhibit good agreement, specially $`u_{1,1}`$ (top right panel). Davis, Nusser and Willick (DNW) find systematic discrepancies between the Mark III and IRAS 1.2-Jy–predicted flow fields, particularly a difference in the $`m=1`$ component of the dipole. Estimating the velocity field from the IRAS PSCz catalogue allows us to understand the influence of sparse-sampling on the 1.2-Jy inferred fields. Since both dipoles are so similar, we might conclude that the discrepancies with the Mark III velocities do not disappear when the shot noise is reduced by a factor $`\sqrt{3}`$, like in the PSCz velocity field.

Willick et al. applied their VELMOD machinery to compare the observed velocities of the Mark III galaxies with those predicted by the IRAS 1.2-Jy model within $`s=\mathrm{3\hspace{0.17em}000}`$ $`\mathrm{km}\mathrm{s}^1`$. They found that the residuals are well modeled by a quadrupole with an amplitude of $`3\%`$ that of the Hubble flow.

In Fig. 3, the quadrupole components inferred from both IRAS catalogues agree quite well, except that the $`m=0,2`$ magnitudes are somewhat larger in 1.2-Jy, therefore a quadrupole mismatch between Mark III and the PSCz velocities is also to be expected. Note that when all the multipoles are considered (i.e. when performing a full v-v comparison) the PSCz and 1.2-Jy gravity fields look fully consistent (B99).

In all plots, shaded and hatched regions represent 1-$`\sigma `$ uncertainties. These incorporate sources of errors in the velocity field: i) the shot-noise error due to sparse-sampling of the underlying density and velocity fields; ii) systematic and random errors caused by a less than perfect method of reconstruction. The uncertainties are evaluated summing in quadrature the errors due to shot-noise and filling-in procedure. We compute the shot-noise error by generating 100 bootstrap realizations of the observed distributions of IRAS galaxies and computing the velocity fields from these realizations. To quantify random and systematic errors caused by the filling-in procedure we use a suite of 20 mock-catalogues that mimic the main properties of the PSCz and 1.2-Jy redshift surveys. A complete description of the prescription followed to quantify the velocity uncertainties can be found in Teodoro, Branchini and Frenk .

Where does the discrepancy between the monopoles of the PSCz and 1.2-Jy surveys come from? As shown in the plots, the difference is larger than that expected from the shot–noise (included in the error budget). Tadros et al. have suggested that the PSCz catalogue may be incomplete for fluxes $`0.7`$ Jy. If true, then we would expect that the velocity monopole for the PSCz with a flux cut at 0.7 Jy (PSCz0.7) would be in good agreement with the 1.2-Jy survey. The dotted line ($`u_{00,\text{PSC}\text{z}_{0.7}}`$) in the top panel of Fig. 1 shows shows that objects fainter than 0.7 Jy are only partially responsible for the monopole mismatch that therefore cannot be ascribed to incompleteness at low fluxes however that this is not case. It is only in cutting the PSCz catalogue at a flux level of 1.2 Jy that, as expected, the discrepancy disappears. This is clearly seen in the right panel of Fig. 1 in which the dashed line indicates PSCz1.2 velocity monopole.

We conclude that a better sampled IRAS catalogue cannot resolve alone the mismatches between the Mark III peculiar velocity field and the IRAS predicted gravity field.

## Acknowledgments

LT has been supported by the grants PRAXIS XXI/BPD/16354/98 and PRAXIS/C/FIS/13196/1998. Thanks do Enzo Branchini and Carlos Frenk for advice on this work as well as comments on the text.