’ J": at“.
"3:5: i7:

,ngu
:4x‘ .: ..
nr\kalv.2...
a) zlttsnta
. 7..

22.1.1.-.
, ‘:J:xa
‘

l ‘ V
. .: Ivar-:33} ,
arﬁxaﬁJ14..§ .
:9. .r. rrrvLit. u

7:73.)..33... a

5.3.1.. I

I... 3:.

Mann}. 1
:v,

5!)... v

.f).a.rl.u..:5r.
. ‘ u. Ila‘Liibp‘I-Ia l I. .t l
1.1.3:... 51.? .I Lt...:..i-lf.u~l
Ivthxh.|~\ls|v .lPtp.
‘ i. as... n ,x.
v»... .s {:1 . .
:35! 33:1: .
53.51 c, 2: #3.»...12... ‘
1.11:... t: r315}-.. , int x
(t AVqu¢I

f

15$? 9.1:!
{1‘53}
\i‘VA|.C
...\L

.\IL|..r|.| 7......- .

iii-3‘1: .11:
\w . r 5 .1.» . ‘
th I .

...‘.|;\ ”US$52 1
rt.» I: ’V-\..5

.1!
x: .7.
J..l!.111«11
»\ yVV.»:a.uw.\A

i. to:

|Ialhxt'.
.l/ [2.x F .c
(a . r ‘
r t .1
9812/.
y. .Iriﬁnyitui
5.29:,

 

 

. ~. -‘ ‘
0" S
{\fﬁ ;_'_::a‘\

3 1293 01563 4805

This is to certify that the

dissertation entitled

Computer Simulations
of High-Energy Heavy Ion Collisions

presented by

Gerd Kortemeyer

has been accepted towards fulﬁllment
of the requirements for

PhD degree in Phys iC S

 

 

a' pro gor-

 

Date pﬂ‘NWQo 23‘ 3?

MSU is an Afﬁrmative Action/Equal Opportunity Institution 0-12771

 

LIBRARY

Michigan State
University

 

 

 

PLACE IN RETURN sex to remaite this checkout trim} your record.
TO AVOID FINES return on or before date due.

DATE DUE DATE DUE DATE DUE

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

MSU Is An Affirmative Action/Equal Opportunity Institution
CMWWM‘WDJ

 

 

 

 

COMPUTER SIMULATIONS OF HIGH-ENERGY HEAVY ION COLLISIONS

By

Gerd Kortemeyer

A DISSERTATION

Submitted to
Michigan State University
in partial fulﬁllment of the requirements
for the degree of

DOCTOR OF PHILOSOPHY

Department of Physics and Astronomy

1997

 

 

ABSTRACT
COMPUTER SIMULATIONS OF HIGH-ENERGY HEAVY ION COLLISIONS
By

Gerd Kortemeyer

One of the still most challenging questions in nuclear physics is that of the equation
of state (EOS) of nuclear matter — how does nuclear matter change its properties
under different temperatures and pressures? Is there a phase transition to quark
gluon plasma? Heavy-ion collisions are of great importance in this study; when two
large nuclei collide, the nuclear matter is compressed, and due to two—body collisions
of the nuclear constituents the temperature rises. Unlike for macroscopic matter, for
nuclear matter this state is unfortunately not directly observable: both the length
and the time scale of nuclear collisions prohibit probing nuclear matter while it is
Still at high temperature and density. Instead, the nuclear matter already cooled,

expanded, and formed more stable conﬁgurations before detection.

However, these late reaction products still carry information that helps to recon-
struct What happened in the early stages of the collision. It is the task of simulations
to aid the reconstruction of the reaction mechanism and conditions from the late

reaction products.

In this theSiS, Computer simulation tools for nuclear reactions are used to examlne

the - . . .
COnWlllences of different microscoplc models.

 

 

 

 

To my wife Anna

iii

 

Contents

LIST OF TABLES viii
LIST OF FIGURES ix

1 Introduction
1.1 Intermediate Energy ...........................
1.1.1 Percolation Models ........................
1.1.2 Transport Models .........................
1.2 Ultrarelativistic Regime .........................
1.2.1 Experimental Efforts .......................
1.2.2 Parton—Cascades .........................

#CIOODMI—“i—‘l-I

 

2 ISOSPin dependent multi-fragmentation in 112Sn + 112Sn and 124Sn +
124Sn collisions

2.1 Introduction ................................ 6
2.2 Implementation of the Isospin Dependent Percolation Model ..... 10
2.3 Results and Conclusions ......................... 14

3 Nuclear Flow in Consistent Boltzmann Algorithm Models 20
3.1 Introduction ................................ 20
3.2 Theoretical Background ......................... 24
3-3 Results and Conclusion .......................... 30
3.4 Outlook .................................. 36

4 Causality Violations in Cascade Models of Nuclear Collisions 38
4.1 Introduction ................................ 38
42 Macroscopic Causality Violations .................. ,. - 40
4-3 Microscopic Causality Violations .................... 45
4-4 Flame-Dependencies and Time—Ordering ................ 47

iv

 

 

 

4.4.1 The Transformation Equations .................. 48
4.4.2 Finding the Boost Parameters .................. 49
4.4.3 Transformation of Momenta and Energy ............ 49
4.4.4 Transformation of the Relative Positions ............ 49

4.4.5 Ways to Calculate the Positions of the Partons in the Next
Timestep ............................. 50
4.4.6 Consequences ........................... 55
4.5 The Parton Model ............................ 59
4.6 Comparison of Methods for Dealing with Superluminous Signals . . . 61
4.6.1 Model A: Restrictions on the Signal Velocity .......... 71
4.6.2 Model B: Suppressing of Low Energy Collisions ........ 71
4.6.3 Model C: First Steps Towards a Retarded Interaction ..... 72
4.6.4 Model D: Downscaling the Cross Sections ........... 73
4.6.5 Model E: Wee-Partons ...................... 74
4.6.6 Model F: Proper time approach ................. 75
4.7 Conclusions ................................ 77
A The Parton Cascade Program 79
A.1 Flow chart ................................. 79
A2 Test runs ................................. 79
A3 Scattering ................................. 81
A4 Assignment of the particles to the boxes ................ 82
A5 Checking for collisions .......................... 82
A6 Decay of particles ............................. 83
A7 Creating and deleting particles ...................... 84
B PrOgram listing 85
B.1 STDEC .................................. 85
B.1.1 Physics parameters ........................ 85
B.1.2 Other parameters ......................... 86
B.1.3 Particle data ........................... 86
B.1.4 Other parameters ......................... 88
B.1.5 Grid information ......................... 89
B.1.6 Initialization parameters ..................... 90
B.1.7 Information block for the initialization methods ........ 90
13-1-8 Data for two particles that are scattering candidates ----- 90
B19 Already scattered ......................... 92

v

 

 

B.1.10 Statistics .............................. 92

 

B.1.11 Information for particles to be created ............. 92
B112 Permutation ............................ 93
B.1.13 Timing information ........................ 93
B.1.14 Controls .............................. 93
B.1.15 Charges .............................. 93
B.1.16 Cross sections ........................... 93
B.1.17 Parameters for the cross sections ................ 94
B.1.18 Causality information ....................... 94
B.2 The main ﬁle ............................... 94
B.2.1 Local variables .......................... 94
B22 Some initializations ........................ 94
B23 The end .............................. 98
B24 The include section ........................ 99
B.3 INIT .................................... 100
B.3.1 Subroutine init .......................... 101
B.3.2 Subroutine makenucleus ..................... 103
B.3.3 Subroutine initﬂag ........................ 104
B.3.4 Subroutine initcharge ....................... 104
B.3.5 Subroutine initcheck ....................... 105
B.3.6 Subroutine initperm ....................... 105
B.3.7 Subroutine initcross ........................ 106
8.3.8 Subroutine rancor ......................... 109
B.3.9 Subroutine restat ......................... 109
BA INITMETHl ............................... 110
B.4.1 Local variables .......................... 111
B.4.2 Initialization of the particles ................... 111
8.4.3 Creating particles ......................... 112
B.4.4 The end .............................. 114
B.5 MOVE ................................... 114
B.5.1 Local variables .......................... 114
B.5.2 Determining the new gridsize .................. 115
B.5.3 Moving the particles ....................... 115
B.5.4 Assigning the particles to the boxes ............... 117
B.5.5 The end .............................. 118
B6 SCATTER ................................. 118

Vi

 

 

 

B.6.1 Local variables .......................... 118
B.6.2 Initialization of variables ..................... 119
B.6.3 Choosing the boxes ........................ 119
B.6.4 Combination of all particles from the boxes ........... 121

B7 CAUSCATTER .............................. 122
B.7.1 Local variables .......................... 122
B.7.2 Combination of all particles from the boxes ........... 123

B8 COLLIDE ................................. 123
B.8.1 Subroutine collide ......................... 124
B.8.2 Subroutine outchannel ...................... 132
B.8.3 Subroutine totalcross ....................... 135
B.8.4 Subroutine xsecentry ....................... 136

B9 TWOTOTWO .............................. 139
B.9.1 Subroutine twototwo ....................... 139
B.9.2 Subroutine montecarlo ...................... 144
B.10 NEWMOMENTA ............................. 145
B.10.1 Local variables .......................... 145
B.10.2 Construction of a coordinate system .............. 146
B.10.3 Getting the angle (,0 ........................ 146
B.11 DECAY .................................. 148
B.11.1 Local variables .......................... 148
B.11.2 Decaying of the particles ..................... 148
B.12 CREDEL ................................. 149
B.12.1 delpar ............................... 149
B.12.2 crepar ............................... 152
B.12.3 checknumbers ........................... 154
B.13 PERMUTE ................................ 156
B.13.1 Local variables .......................... 156
B.13.2 Permutation ............................ 156
3.14 DISTRI .................................. 157
C Index 160

 

List of Tables

3.1 Impact Parameter Averaged Analysis (b S 5 fm) and Experimental Data 35

4.1 Parameters used for the comparisons .................. 64

 

viii

 

 

 

 

List of Figures

2.1

2.2

2.3

2.4

2.5

3.1

3.2
3.3

3.4

Average number of intermediate mass fragments (IMFs) versus number
of charged particles (left panels) and neutrons (right panels). The
full circles denote the experimental results (Kun96] for a 12‘1Sn+12‘1Sn
collision, the open circles the results for a 1 2Sn+1128n collision, both
at 40 MeV/ A. The solid lines represent percolation simulation results
for the heavier isotopes, the dashed lines for the lighter ones. The
top row was calculated without any stability mechanism, in the second
row, an evaporation mechanism is employed, in the third row a ﬁssion
mechanism, the fourth row was achieved with a combination of both.

BUU simulations of the 40 MeV/A Sn—collisions at different impact
parameters. The top plot was done with b = lfm, the middle one with
b = 3fm, and the bottom one with b = 5fm. The left plots show the
initial conﬁguration, the right plots the same scenario 200 fm/c later.

The nucleons in the participant (overlap) zone of the nuclei are dis-
tributed on a rectangular lattice. Bonds within the participant region
are randomly broken, for the resulting fragments, different stability
mechanisms are employed. ........................

Experimental data [Tsa93] and percolation results for a 197Au+197Au
collision. The experimental data is given for E/A=100, 250 and 400
MeV, the theoretical predictions refer to different bond—breaking prob—
abilities p0. ................................
Results of a simulation with different bond-breaking probabilities for
equal and different nucleon bindings (PM = 0.6; p08 2 0.9). ......

Overview of techniques to simulate heavy—ion collisions .........

12

17

23

Additional displacement of scattering partners due to their ﬁnite volume 25

Reduced volume due to ﬁnite particle radius. In previous implementa-
tions, the testparticles were pointlike (left), which left more free space
between the particle. The current implementation (right) takes into
account the volume taken up by the testparticles. ...........

Screening effect due to the ﬁnite volume of the testparticles. On the
left, the pointlike particles 1 and 2 could scatter with each other, par—
ticle 3 is not “in the way.” On the right, due to the ﬁnite volume, they
are prevented from scattering, particle 1 would scatter with particle 3
ﬁrst. ....................................

 

3.5 The scattering enhancement factor YE as a function of the average
cross section. Shown is YE for g = 1, 2, 390. ..............

3.6 The nucleon conﬁguration in the reaction plane at different times dur-
ing the collision, shown is the nuclear density g(z, y = 0, z) .......

3.7 Collision rate for a 250 and 400 MeV (Au,Au) collision with b : 3 fm
versus time. ................................

3.8 Average ﬁnal transverse momentum versus reduced rapidity of the pro-
tons in a 250 MeV (left panel) and 400 MeV (right panel) (Au,Au)
collision with b = 3 fm ...........................

4.1 One particle coming from the left scatters with another particle coming
from the right. The information from this event is rapidly passed on
through a chain of subsequent scattering events. ............

4.2 Distance d of scattering events versus simulation time at constant cross
sections of a 20.1, 0.2 and 0.5 fm2. ...................

4.3 Lightcone and 3—dimensional hyperplane (2—dimensional illustration) .

4.4 Minkowski Plot. Highlighted is the timestep in which the particles
are at the point of closest approach. Note the time-order of the four
events “Particle 1(2) at the beginning(end) 0f the timestep” in the
c.m.-frame. These points in time are given by the intersection of the
c.m.~t axis with the dot-dashed lines that are parallel to the c.m.—z axis.

4.5 Minkowski Plot. Highlighted is the second timestep after the particles
had reached their point of closest approach. Note the the time-order of
the four events “Particle 1(2) at the beginning(end) of the timestep”
in the c.m.-frame has changed .......................

4.6 Two partons passing their point of closest approach in the c.m.—frame.

4.7 Possible discrepancies in the time-ordering of scattering processes that
are usually prevented by the spacelike distance of the scattering partners.

4.8 Typical contributions to the gluon-gluon cross section. The four point
diagram on the far right leads to divergencies in the cross section that
cannot be handled by phenomenological cut-off masses. ........

4.9 Example of partial cross sections being used versus the c.m.—energy.

4.10 Total number of scattering events versus timestep length used for the
simulation. The dashed line interpolates the results from simulations of
(La,La) collisions, the dotted line from (Cu,Cu) collisions. The symbols
denote the actual data points generated in the respective simulations.

4.11 In the left panel the projection of the initial parton conﬁguration on the
x, z-plane is depicted. The highly Lorentz-contracted proton is moving
to the left into the equally contracted gold nucleus which is moving to
the right. The right panel shows the conﬁguration 1 fm/c later.

4.12 Results of simulations with different mechanisms to suppress superlu-
minous signal transport (Part 1). ....................

4-13 Results of simulations with different mechanisms to suppress superlu-
minous signal transport (Part 2). ....................

X

32

34

41

43
46

53

54
56

57

60
61

62

63

66

68

 

 

4.14 A ﬁrst step to the introduction of retarded interactions ......... 73

4.15 Initial conﬁguration and conﬁguration after 2 fm/c for a simulation
with wee-partons .............................. 75

A.1 BUU program ﬂow chart ......................... 80

xi

 

 

Chapter 1

Introduction

1.1 Intermediate Energy

1.1.1 Percolation Models

In intermediate energy heavy—ion collisions, the production of complex fragments still
is a puzzling aspect of nuclear physics research. Experiments have been conducted
both with proton and with heavy-ion beams, and several seemingly different phe-
nomena, such as evaporation, spallation, ﬁssion and multi—fragmentation have been

observed.

In multi-fragmentation scenarios, the hot nuclear system tends to decay into mul-
tiple intermediate mass (3 3 Z s 20) fragments (IMFs). Signiﬁcant differences in
the relationships between fragment, neutron, and charged particle multiplicities were
found between 112Sn+1128n and 124Sn+mSn collisions at 40 MeV/A [Kun96]. The re-
sults are incompatible with a universal scaling of the average number of intermediate
mass fragments versus the number of charged particles and neutrons, respectively,
and in Chapter 2, the possibility to explain this phenomenon in the framework of
percolation models is explored. It was found that the results are only reproducible in
part, which indicates contributions of nuclear phenomena beyond geometrical frag—

mentation schemes.

 

1.1.2 Transport Models

For a long time, simulations of heavy-ion collisions have been rather phenomenological
in nature, which somewhat corresponded to early experiments being conducted in
an inclusive way - phenomenological models, even if based on different assumptions,
would predict the inclusive data with oftentimes almost equal quality. As experiments
became more sophisticated and exclusive, the need for advanced simulation methods
arose. This was when microscopic models became more popular: one has to be able
to dynamically simulate the collisions without any assumption concerning thermal

equilibrium — transport models are one of those microscopic dynamical models.

Boltzmann—Uehling-Uhlenbeck (BUU) (also called Landau—Vlasov) codes are semi-
classical simulations, in which as an ansatz for the Wigner-function a product of
Delta-functions is used. Each Delta-function can be associated with a localized “test—
particle.” Inserting this ansatz into the Boltzmann equation results in equations of
motion for the testparticles that are governed by two terms: a smooth modiﬁcation of
the trajectories resulting from a mean—ﬁeld term, and a collision term that describes
two-particle collisions between testparticles. The actual solution of the Boltzmann
equation is achieved by evolving the testparticle essemble in phase—space according
to the equation of motion rather than solving a partial differential equation. Chapter

3gives a more detailed overview of this mechanism.

The semiclassical transport codes work reliably at energies beyond 100 MeV per
nucleon. Below this energy, quantum effects are too pronounced to be neglected.
The BUU codes followed earlier transport codes, so—called cascades,[Cug81] which
did I10t have a mean-ﬁeld term yet. This approximation is only appropriate at even
higher energies, above around 1 GeV per nucleon (Appendix A describes the details

Ofa Cascade code). At lower energies, many two-body collisions are blocked due to

 

 

phase-space considerations, so that the mean-ﬁeld term becomes dominant. At higher
energies, Pauli—blocking is not so much of an issue anymore, since a large phase—space
volume is available, which strengthens the inﬂuence of the collision term. At the

same time, because the particles have larger momenta, the mean-ﬁeld becomes less

inﬂuential.

Transport models have successfully described many aspects of intermediate en—
ergy heavy-ion collision dynamics. One of the observables in heavy-ion collisions is
the collective flow of the reaction products, it contains information about the com—
pressibility of nuclear matter. However, it was indicated [Ale95] that within current
implementations of BUU simulations the equation of state is not implemented in a
consistent way. In Chapter 3, a hard core equation of state for the nuclear matter is

implemented and its effect on the simulated collective ﬂow is explored.

1.2 Ultrarelativistic Regime

1.2.1 Experimental Efforts

Currently, two colliders are planned for the acceleration of heavy—ions to ultra—relativistic
energies. The Relativistic Heavy-ion Collider (RHIC), presently under construction
at Brookhaven National Laboratory in New York, is a dedicated heavy-ion collider
planned for experiments in 1999. RHIC will accelerate and collide ions from protons
to heavy nuclei, such as Au, at cm. energies up to 500 GeV for protons and 100 GeV
per nucleon for Au+Au collisions. The luminosity for Au+Au will be 2 x 1026 cm“2
3‘1- Near head-on collisions of Au+Au at RHIC are expected to produce from 500

t01500 charged particles per unit pseudorapidity at midrapidity in a single collision.

Large detector systems are being constructed to analyze the products of these

inteFactions to observe indications of the possible formation of a quark—gluon plasma

 

andaPOSSIble chiral Phase transition.

Heavy—ion physics research will also be an integral part of the program for the
Large Hadron COllider (LHC), to be constructed at CERN, the European Centre for
Nuclear Physics, in Geneva, Switzerland. For Pb nuclei, the cm. energies at the
LHC Will be 5-4 TeV per nucleon pair with luminosities of 1027 cm“2s“1. Predictions
for the charged particle densities at the LHC for near head-on collisions of Pb+Pb
range from 2000 to 8000 per unit pseudorapidity. The large uncertainty in these

numbers arises primarily from the present lack of information on the distributions of

soft gluons in nuclei.
1.2.2 Parton—Cascades

It has been indicated that RHIC and LHC run into somewhat of a “theory-vacuum.”
The interpretation of these complex collisions and the possible generation of a quark—

gluon plasma pose a major problem: what are the experimental signatures? Currently,

several possible signatures are in discussion:
0 Kinematic effects of a ﬁrst-order phase transition.

0 Electromagnetic probes e+e’, iﬁu‘, 7 measure the early phase of quark and

gluon distributions and content.
0 Hadronic probes of deconﬁnement and chiral symmetry restoration

J/i/J ——> 6+ 5 (Up, \I” suppression)
gg —> 35 (strangeness enhancement)
None of these possible signatures have been commonly acknowledged as reliable. In

an effort to aid the answering of this question, a theoretical model for the collision

processes that goes beyond a phenomenological description must be developed. One

 

 

 

 

5

POSSlble microscopic approach is to extend the semiclassical transport theory to high-
energy physiclea193, 80189.1, Sor89.2, Gei92.1, Gei94.1, Gei92.2, Gei93], which leads
to so-called “Dalton Cascade” codes. As the energies increase in these models to the
ultrarelativistic regime, Lorentz covariance and causality are not strictly respected.
The standard argument is that such effects are not important to ﬁnal results; but
they have not been seriously considered at high energies. In Chapter 4, it is pointed
out how and Why these happen, how serious of a problem they may be and ways of
reducing or eliminating the undesirable effects are suggested. The appendix gives a

commented listing of the major sections of the parton cascade code that was developed

as part of this thesis.

 

Chapter 2

lsospin dependent
multi—fragmentation in 112Sn —|—
112st; and 124Sn + 124sn collisions

2.1 Introduction

Percolation models have proven highly successful in the simulation of multi—fragmen—
tation reactions in the past [Bau84, Cam85]. Within these models, fragmentation is
described by ﬁrst distributing a set of points or sites, each representing a nucleon,
on a 3—dimensional lattice, which represents the bonds between the sites. In the
case of a simple rectangular lattice, each site is connected to six nearest neighbors,
however, it has been shown that the model is to a large degree independent of the
lattice structure [Bau84, Sta79]. In the second step, randomly some lattice bonds
are broken with a probability that in non—isospin dependent percolation models is

the only free parameter. The remaining connected clusters are identiﬁed with the

fragments of the reaction, the bond-breaking probability with the excitation energy
per nucleon [Con79].

In this chapter, the percolation model of Bauer et a1. [Bau84] is modiﬁed by

the explicit inclusion of isospin degrees of freedom, i.e., the lattice is comprised of

 

 

 

 

 

7

protons and neutrons inStead of just nucleons, in an attempt to reproduce the experi-
mental data gained in the comparison of the multi—fragmentation in 112Sn+112Sn and
124Sn+124Sn 0011-1310118 [Kun96]. We especially focus on the average number of inter-
mediate mass fragments (IMFs, 3 3 Z S 20) (NIMF) versus the number of charged
partiCleS No (Fig. 2.1, left panels), and versus the number of neutrons Nn (Fig. 2.1,
right panels) detected. The full circles denote the results for the 124Sn+1248n reac-
tion, the circles the results for the 112Sn+1128n reaction. The striking feature about
these distributions is the “Splitting” of (NIMF)(NC), and the position of the maxima in
(N,MF)(N,,). Both do not agree with common multi-fragmentation models, in which
the distributions (N,MF)(NC) should lie on top of each other, and the positions of the

maxima in (Nu/IF) (Nu) Should simply correspond to the ratio of neutron abundances

in the respective isotopes.

 

Figure 2.1: Average number of intermediate mass fragments (IMFs‘) versus number of
charged particles (left panels) and neutrons (right panels). The full circles denote the
experimental results [Kun96] for a 124Sn-i—1248n collision, the open circles the results
for a 112Sn—i—“ZSn collision, both at 40 MeV/A. The solid lines represent percolation
simulation results for the heavier isotopes, the dashed lines for the lighter ones. The
top row was calculated without any stability mechanism, in the second row, an evap-
oration mechanism is employed, in the third row a ﬁssion mechanism, the fourth row
was achieved with a combination of both.

 

Figure 2.1

:_::__

 

 

 

 

 

 

 

 

 

 

 

054 210

 

23$
Ah—ZCQV

 

 

 

10

2.2 Implementation of the Isospin Dependent Per-
colation Model

For each simulated collision event, ﬁrst the impact parameter is randomly selected.
Then With a Simple Monte Carlo integration, the number of protons and neutrons in
the overlap zone of the two nuclei is determined. We employ an approximation in
which the nucleons outside the overlap zone are neglected — we found this approxima-
tion to be appropriate by studying Boltzmann-Uehling—Uhlenbeck (BUU) simulations
[Bau86] of the collisions at different impact parameters, which clearly showed distinct
spectator regions in the ﬁnal state even for small impact parameters. Figure 2.2 shows
the outcome of the BUU simulations at different impact parameters. The simulations

were run with b = 1, 3, and 5 fm, respectively, and the spectator regions were observed

even 200 fm/c later.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

11
_ I I 1 I _
I . . . | . . . . | i i . . | . . i . | .
O 10 20 30 4O
_ I 'ﬁ ' I l I ' ' I
__ -» - ' ‘O- . , _3
. l i I i I . I
O 10 20 30 40
' I ' ' l ' ' I ' ' l
I I | I

 

 

Flgllre 2.2: BUU simulations of the 40 MeV/A Sn-collisions at different impact pa-
:ﬁmeters. The top plot was done with b = lfm, the middle one with b = 3fm, and

16 bottom one with b = 5fm. The left plots show the initial conﬁguration, the right
D 0t8 the same scenario 200 fm/c later.

 

 

 

‘Peautor

 

12

3

I:

“0%

' I

a;

$53

. <

E /

4.)
5
3'
U
33
it
04

Fragments

Spectator

 

Figure 2.3: The nucleons in the participant (overlap) zone of the nuclei are distributed
‘le a rectangular lattice. Bonds within the participant region are randomly broken,
or the resulting fragments, different stability mechanisms are employed.

 

 

13

The nucleons in the overlap zone are randomly distributed on a rectangular lattice

(see Figure 2.3).

The lattice bonds are then broken with a probability p, which we determined
intwo diﬂerent ways: in one method we set it equal to a parameter p0 which we
obtain by ﬁtting to the experimental data [Kun96], in the second method we choose 1)
according to a Gauss distribution around p0, i.e., for each simulation, p varies slightly
in order to simulate excitation energy fluctuations. A comparison between both
methods showed no signiﬁcant difference in the outcome except for better statistics in
the latter method for higher event multiplicities as large bond-breaking probabilities
were included. Since we ﬁnd the inclusion of excitation energy ﬂuctuations to be

more realistic, we settled for the latter method.

We then identify clusters of nucleons which are still connected with each other.
However, since those clusters are not necessarily a stable conﬁguration of protons
and neutrons, we experimented with several different algorithms to achieve fragment
stability (see Fig, 2.3). Methods included a re-distribution of protons and neutrons
between the fragments, further ﬁssion of the fragments, evaporation of protons and

neutrons from the fragments, and simulations with no additional stability criteria
applied. Also, the deﬁnition of stability is not obvious: the experimental lifetime
data applies to nuclei in their ground state and is not directly transferable to the
fragments of a multi-fragmentation reaction. Since we found the outcome to only be
slightly dependent on the deﬁnition used, we settled on a stability criterion where
the fragments are required to have a ground state lifetime that is long enough for
them to reach the detectors. As an unfortunate side effect of ﬁssion and evaporation
mechanisms, however, the relationship between the bond breaking probability p0 and
the excitation energy is not obvious anymore, rather the combination of the initial

bond-breaking and further mechanisms leads to an effective bond-breaking probability

 

 

 

197Au + 197Au

 

 

 

 

 

 

 

 

 

12 I I I I I I I I I I I I I
12 I I I I I I
- 0.8 -
E/A=1oo uev
10“ — 10 — — — 0.7 —
— — — - - 0.65 _.__ —
8— J 8 _ ........... 0.6 _./-'/ ’2‘: \—
5 _ /\ _ ..... '2 \_
.iiiﬁi 57-: 0.5 f,
6- 'i i; — a 6 —
i i z
tﬁl" ii 250 - V _
4— 4' ii — 4—
I I5 I
. C _. ._
C
-. 400
._ f .. _ 2
HI .. - ' . ' _.
6IIIllIllllIlllIlll o'I-I'IIIIIIIIIIIIIIII
o 20 4o 6O 80 100 o 20 4o 60 BO 100
Nc NC

Figure 2.4: Experimental data [Tsa93] and percolation results for a 197Au+197Au col—
lision. The experimental data is given for E/A:100, 250 and 400 MeV, the theoretical
predictions refer to different bond—breaking probabilities p0.

that is higher than 190 — therefore, p0 sets a lower limit for the excitation energy.

0f the order of 106 — 107 events were simulated for each set-up; for each individual
event the number of charged particles, neutrons and IMFs was recorded, where in
accordance with the experimental data we employed detector efficiencies of 0.9 for all

Charged particles and 0.65 for neutrons.

2.3 Results and Conclusions

The code was ﬁrst applied to the experimental results of Ref. [Tsa93]. Figure 2.4
shows the the average number of intermediate mass fragments (3 3 Z s 20) versus
the number of charged particles ((N,Mp)(Nc)) for 197Au+l97Au collisions. The exper-
imental results in the left panel refer to different energies per nucleon [Tsa93], the

curves in the right panel to percolation simulations at different bond—breaking proba—

 

15

bilities 100- It had been found earlier [Pha92] that percolation codes generally slightly
underpredict the number of intermediate mass fragments, which was attributed to

the possible existence of non-compact decay geometries. Overall, however, the model

was found to reproduce the data reasonably well.

In the next step, we address the question of the reproducibility of the isospin—
dependence found in Ref. [Kun96]. The top row of Fig. 2.1 shows the result of a
simulation with a bond—breaking probability distribution centered around p = 0.7
and a half-width of 0.1. In this simulation, no stability mechanism is applied. The
solid line corresponds to the simulation for the heavier isotope (experimental data
indicated by full circles), the dashed line to the lighter isotope (open circles). The
difference between the isotopes in the average number of IMFs versus the number of
charged particles could not be reproduced, the outcome basically reﬂects the trivial
autocorrelation that every IMF is a charged particle, the slope is determined by the
ratio of IMFs versus lighter fragments. The linear relationship will eventually break
down in events with high multiplicities when more and more fragments are smaller
than lMFs. The experimental data shows that in the heavy isotope relatively fewer
fragments with Z < 3 are formed, a trend that cannot be seen in the simulation. The
difference in the average number of IMFs versus number of neutrons (right panel)
can be reproduced, which, however, is not surprising: in the collision of the neutron-
richer isotopes, more neutrons are emitted. The positions of the maxima, i.e., at 25
and 29 neutrons, respectively, correspond to the ratio of neutrons in the isotopes,
62 and 74, respectively (25/29 x 0.86; 62/74 x 0.84). In the experimental data the
maxima are at 25 and 40 neutrons, the ratio of z 0.63 is in—compatible with the simple
explanation above. At high neutron numbers, statistics get rather unsatisfactory, and
we do not reproduce the high neutron multiplicities seen in the experiment for the

heavy isotope.

 

 

16

The second row 0f Fig. 2.1 shows the outcome of a simulation with an evaporation
mechanism to achieve fragment stability: protons and neutrons are broken off the
fragments until the remainder is stable. In the simulation shown, in case of N > Z,
neutrons are broken off, and vice versa. In another simulation, protons and neutrons
had been broken off randomly, which lead to slightly less IMFs. Again, differences
inthe IMF—distribution versus number of charged particles (left panel) could not be
observed, however, for both isotopes, in comparison to the simulation without stability
criterion,the ratio of IMFs to lighter particles decreased. In the IMF distribution
versus number of neutrons (right panel), naturally higher neutron multiplicities are
observed, the ratio of the maxima positions (32/40:0.8) remains compatible with the

ratio of neutrons between the isotopes, though.

The third row of Fig. 2.1 results from a simulation a ﬁssion mechanism: an unsta-
ble fragment is broken into two fragments, if such a secondary fragment is unstable, it
isagain broken up into two fragments, and so on, until only stable fragments remain.
The distribution of the secondary fragment sizes is chosen to be centered around 0.5
and falls off quadratically towards 0 at 0 and 1; a simulation with a flat distribution
yielded similar results, it had only very slightly less lMFs. Both these mechanisms ﬁt
the IMF distribution versus number of neutrons rather well, the ratio of the positions

of the maxima, 31/41 R: 0.7 is half-way between the expected 0.83 and the experi—
mental 0.6, however, the different heights of the maxima could not be reproduced.
The IMF-distribution versus charged particles again fails to show differences between

the isotopes.

The forth row ﬁnally shows the result for a combination of the two mechanisms
above, each unstable fragment undergoes evaporation or ﬁssion with equal probability.
Apparently, this mechanism produces too many light fragments, both charged parti-

cles and neutrons. The ratio of the maxima in the IMF distribution versus number

L—_____,_

 

 

 

 

 

 

 

l ' 'F I l l I l l *7 l
5 —— 0124811 ”ﬂu-'3 __ _
112 .9
A4— —--- S .O. 699955830 ,”I/
Elli .. see /’///’
3: ..0 see , / ’ ’ ,
2' .0. ’/ I
l‘ .99998
0 I; l . I l
0 10 20 30 40
NC

 

Figure 2.5: Results of a simulation with different bond-breaking probabilities for equal
and different nucleon bindings (P04 : 0.6; p06 2 0.9).

of neutrons is 33/ 4420.75.

ln deviation from the established percolation models, we also worked with dif—
ferent break—up probabilities for bonds between protons and protons, neutrons and
neutrons, and protons and neutrons, that is, a probability pee for equal pairings, and
aprobability pad for unequal pairings, 3906 2 120,1. As it turns out, in this method the
outcome depends signiﬁcantly on the distribution of protons and neutrons on the lat-
tice; in the case of a highly ordered conﬁguration where except for the excess neutrons
both types of nucleons are distributed in an alternating way (“salt crystal”), i.e. in
general every proton has six neutrons as nearest neighbors and vice versa, the ex-
cess neutrons being put in randomly as “impurities,” one notes signiﬁcant differences
between both Sn—isotopes in the distribution of intermediate mass fragments versus
charged particles — however, these are still much smaller than in the experiment, as
Fig. 2.5 illustrates this for pee = 0.9 and pod = 0.6. Even though this effect leads to
aslightly closer resemblance of the experimental data, we consider it to be an artifact
since it nearly completely vanishes with a purely random distribution of protons and

neutrons on the grid. We attribute this effect to the fact that in a pure “salt-crystal”

 

 

18

conﬁguration every bond has the break-up probability Poa; every impurity will in gen—
eral lead to the introduction of six bonds with break-up probability p05, and therefore
11353 large impact. As a result, disregarding surface effects, in the collision between
the lighter isotopes about 79% of the bonds are of type Poa and 21% of type p05, while
for the heavier isotopes the percentages are 61 and 39, respectively. In the random
conﬁguration, however, there is no such ampliﬁcation, the percentages are 58 versus

42 for the lighter isotope, and 56 versus 44 for the heavier isotope.

Also, simulations were run with a “neutron skin,” which inﬂuenced the ratio
of protons and neutrons in the overlap zone for different impact parameters: for
small impact parameters the ratio of neutrons to protons was higher than for large
impact parameters, the size of the neutrons skin were determined by a Hartree—Fock
calculation[Br096] . However, also these simulations could not improve the agreement

with experimental data.

In conclusion, the experimental results could only be reproduced in part. The

main discrepancies are:

0 The difference between the two isotopes in the average number of IMFs versus

number of charged particles could not be reproduced.

o The difference in the maximum values of the IMF distributions versus number

of neutrons could not be reproduced.

o The experimental positions of the maximum number of IMFs versus number
of neutrons is in-compatible with the trivial shifting due to higher neutron

abundance that is found in the simulation.

The discrepancies found between the data and this basically geometrical approach

indicate that effects outside of percolation theory are important. The nuclear struc—

 

 

 

19

ture 0f the fragments as well as sequential feeding might play a role. Most important
however seems the role of preequilibrium emission which may not only effect the

sorting axis but as well determines the N / Z composition of the fragmenting system.

 

 

Chapter 3

Nuclear Flow in Consistent
Boltzmann Algorithm Models

3.1 Introduction

Heavy-ion collisions are of great importance in answering the question how nuclear
matter changes its properties under different temperatures and pressures.[Cze86]
When two large nuclei collide, the nuclear matter is compressed, and due to two—
body collisions between nuclear constituents also the temperature rises to about 50 —
100 MeV (about 5 to 10 times 1011 Kelvin). Unlike for macroscopic matter, for nu-
clear matter possible phase transitions etc. are unfortunately not directly observable:
both the length (about 10 fm=10‘14 m) and the time scale (about 60 fm/c =2 x 10‘22
seconds) of nuclear collisions prohibit probing the nuclear matter while it is still at
that high temperature and density. Instead, the nuclear matter has already cooled
down, expanded, and formed more stable conﬁgurations before it hits the detectors;

every observation can only be indirect.

However, these late reaction products still carry information that helps to recon-
Struct what happened in the early stages of the collision. For example, the collective
ﬂow of the reaction products contains information about the compressibility of nu-
clear matter. With the exception of the merely theoretical case of a head—on collision

20

 

 

21

0f the two nuclei (Vanishing impact parameter), from straightforward considerations
itisexpected that in the collision nuclear matter of the two partners is pushed out-
wards perpendicular to the beam within the plane deﬁned by the impact parameter.
The harder in the sense of compressibility the nuclear matter would be, the stronger
this effect would be expected to be. Also, from this very simple consideration one
would expect the ﬂow to always increase with the beam energy. However, since the
interaction between nucleons in addition to a short range repulsion has a long range
attraction, this Simple picture does not hold true: at low energies, the long range part
of the nuclear interaction dominates, the nucleons are attracted, and ﬂow is gener—
ated by the nucleons being scattered towards each other. At energies around 70—140
MeV/nucleon depending on the mass of the reaction partners both effects cancel out
each other, the flow disappears in the experiment; simulations to—date tend to under-
predict that value.[Wes86] At high energies, the nucleons collide at short range, and

the mechanism for ﬂow indeed corresponds to the simple consideration above.

To extract quantitative information about and gain insight into the microscopic
processes within nuclear matter, however, it is necessary to have a dynamical model to
simulate the collisions themselves. Comparision between the predictions for the ﬁnal
stage gained by the Simulation on the one hand, and the actual experimental results
on the other hand are used to reﬁne the model for the early stages. The simulation
of heavy-ion collisions unfortunately is equivalent to solving a quantum—mechanical
many-body problem, which to date is not fully possible. Diﬁerent approaches have
been made to nevertheless approximate the solution. Early approaches were based
on hydrodynamics, they used a nuclear equation of state, but the mean free path
of the nucleons was assumed to be so small that the nuclei in the collision merely
resembled splashing droplets. At the other extreme, there were models that due to

their extremely long mean free path basically resembled colliding gas clouds.[Da593]

in-

22

To model the situation more realistically, testparticle based models were devel—
oped. Here, a nucleon and other nuclear constituents are represented by a number
of testparticles, their individual trajectory is followed rather than being concerned
only with the global properties of some nuclear “ﬂuid” or “gas.” Obviously, these
microscopic methods, which can be implemented in different ways, have the general

property of being computationally more intense than the macroscopic ones.

One method is that of Molecular Dynamics,[Aic91] in this method both the long-
range attractive (soft) and the short-range repulsive (hard) part of the particle interac—
tion is parametrized by potentials, the trajectories of the testparticles are continously

updated in response to the local potential.

Another semi-classical particle-based method is the Boltzmann-Uehling-Uhlen-
beck (BUU) approach. Here the soft part of the interaction is represented by mean
ﬁelds, while the hard part is given by an explicit collision term. The collision term
itself can again be represented in different ways, especially the criteria for a collision
to happen are model dependent. In many codes, this decision is based on geometrical
considerations, for example, a collision is generated at the point of closest approach
between two particles. The BUU method is being used by various groups.[Ber88,
Aic96, Bau86]

One particular implementation of the collision term in BUU codes is the Direct
Simulation Monte Carlo approach (DSMC), see for example Lang et al.,[Lan93] and
Danielewicz.[Dan95] In this approach, collisions between the testparticles are not gen—
erated through geometrical and particle—trajectory based criteria, but stochastically
in a way that the correct collision rate is reproduced. In a collision only momenta
and energy of the particles are changed, while the particles themselves stay in place

until the next advection step a the particles are assumed to be pointlike.

23

 

macroscopic Hydrodynamic models etc. I

Molecular Dynamics
Geometric scattering

criterium
Hard-sphere

   
 
   
   

 

Hard interactions
through
potentials

 
 
 
 
 
 
 
 

 
  
    
 
  

microscopic

Explicit collision

  
 

Stochastic
scattering criterium equation of state

Figure 3.1: Overview of techniques to simulate heavy—ion collisions.

It was suggested that in order to reproduce a Hard-Sphere Boltzmann Equation,
the DSMC approach should be extended by an additional advection that should
take place after any collision,[Ale95] and by a modiﬁcation the collision probability

itselfiAle95, Res77] (Enskog Theory).
Figure 3.1 summarizes these different techniques.

This advection is supposed to push the collision partners away from each other
according to their hard-sphere radius, and the collision probability is adjusted to take
into account the excluded volume of the hard spheres and screening and shadowing ef-
fects. The modiﬁed DSMC method is called Consistent Boltzmann Algorithm (CBA).
The modiﬁcations ensure a non-vanishing second virial and change the equation of
state for the scattering process from that of an ideal gas to that of a hard-sphere gas;
their effect on the calculated value of directed nuclear collective flow in heavy ion

collisions is analysed, and it is found that the flow slightly increases.

«~AM

 

24

3.2 Theoretical Background

In the DSMC approach, the positions and momenta of the particles are evolved in
a two-step process, namely advection and collisions, corresponding to one timestep
of the simulation. During the advection step the particles are propagated according
to their momenta. During the collision step ﬁrst the particles are sorted into spatial
cells of volume V. Then out of the n particles within a given box, at random, m

combinations are chosen and scattered with the probability

o(\/§)v,e1At n(n — 1)/2

W:
NV m ’

(3.1)

where o(\/§) is the energy-dependent elementary hadron-hadron cross section, N is
the number of testparticles representing one nucleon in a full—ensemble testparticle
algorithm [Wel89], At is the timestep length, and vrel is the relative velocity of the
particle pair [Lan93]. In the limit V —> 0, At —> 0, N —> 00, the solutions of this
method have been shown to converge to the exact solution of the Boltzmann equation

(without mean ﬁeld contribution)[Bad89],

8
atfl + ’U ' 'af—rl = (atf1)coli ‘ (3-2)

Here, f1(r, 'v; t) is the one particle distribution function , which, according to Boltz—

mann, is deﬁned in such a way that
f1(r,v;t) d'r do

is the number of particles that at a time t are within a volume element dr around 7'
and have a velocity do around '0. The inﬁnitisimal volume elements have to be taken
with a grain of salt, since taken literally, within a mathematically inﬁnitisimal volume
element, f1 could either be ’1’ or ’0’ v the pointlike particle could either be inside the

inﬁnitisimal volume, or not. Rather, the volume element has to be a physical element

‘Qx

 

 

Figure 3.2: Additional displacement of scattering partners due to their ﬁnite volume

that is microscopically inﬁnite — small enough so that physical properties do not vary

appreciably over its extension, but big enough to contain a large number of particles.

This approach does not take into account the ﬁnite size of the nucleons; the
testparticles are point-like, and if it was not for the contribution of the mean ﬁeld, they
would be following an ideal gas equation of state. It has therefore been suggested by

Alexander et a1. [Ale95] to include an extra displacement d of the collisions partners,

1 vi" _ v?” 0(\/§) (33)

:2|v’T—v,| 7r ’

 

v, = '01 — 112 being the velocity difference before, and v’, : v’1 — v’2 being the velocity
difference after the collision. Particle 1 is displaced by d and particle 2 by —d. This
additional advection pushes the nucleons apart according to their hard-sphere radius.

Figure 3.2 illustrates this mechanism.

Within the Boltzmann equation, this displacement corresponds to change of the

collision term

(atf1,2)c011=/dv1/dvg/dviw(v2ivlivl2avi) [fi,2fi,1 — f1,2f1,1] , (3.4)

 

 

 

26

W again being the scattering probability, and the f ’s being the respective one-particle

distribution functions before and after the scattering for the two particles, i.e.,
f1,1 = f1(r,v1;t), ff), = f1(r,'v’1;t) for Particle 1, (3.5)
and f1,2 = f1(r,'v2;t), fig 2 f1(r,v’2;t) for Particle 2. (3.6)
With the additional displacement, it is
f1 (r,v2; t) f1 (7301) —> f1(r + d,v2;t)f1(r — d, '01; t) . (3.7)
The same replacement has to be done with d —> —d for the inverse term f{)2f{,1.

It is not obvious right away how this displacement scales with N. However, as
in the mean free path 1/((o/N)(Ng)) of a testparticle, 9 being the nuclear density,
there is no N dependence, the average number of collisions that a certain testpar-
ticle is involved in is independent of N. Therefore, in order to achieve the same
total displacement during the course of the simulation, the individual displacement
per collision should not depend on N either. The testparticles are therefore pushed
apart according to the nucleonic radius, and not according to the effective testparticle

radius.

The ﬁnite radius of the particles also makes it impossible for one particle to be
within the “hard core” of the others, and thereby from the available volume V a
fraction

% - giro3 = Vg- gm? (3.8)
is occupied, where a is the average “hard core radius” of the nucleon, and n/N is
the number of nucleons in the respective box. If one only looks at the centers of the
nucleons, they cannot be closer together than 2a ~ in that respect, another view of
the scenario would be to assign each “point nucleon” a “sphere of inﬂuence” with

radius 2a. We obtain a by randomly picking it particle combinations for a respective

box and calculating their cross sections o,(\/s—,).

27

The Pauli principle is approximately taken into account in the following way:
since collisions between particles from the same nucleus that had never scattered
before are assumed to be Pauli-forbidden, in the determination of the average cross
section only particle combinations with partners that are not from the same nucleus,
unless at least one of them had scattered before, are taken into account. According
to the Golden Rule, the cross section is proportional to the available phase space,
however, since the particle combinations chosen for the determination of the average
cross section were chosen without any further restrictions, a momentum space volume
of 4/37r(p3 +pp)3, where p3 is the momentum of the beam per nucleon, and pp is the
Fermi momentum, was implicitly assumed. In reality, the two nuclei occupy part of
the momentum space, namely, 2-4/37rp32. Thus, we adopt the following approximation

for the effective radius a, i.e.,

. = iii (ggwa) . (1 _ 2(———~)—) . (3.9)

For small Ebb, this effective radius is about 0.84 fm, in the range between 100—400

 

MeV it is about 0.47 fm. This is slightly larger than what has recently been suggested
in Ref. [Dan96], there, the effective radius derived from delays in elementary processes

is about 0.6-0.8 fm and 0.15-0.3 fm, respectively.

Due to the reduced volume

V = (1— g- grail) V (3.10)

alone, a modiﬁed higher scattering probability

W: ZW (3.11)
V

has to be used — the particle have less free space between them and therefore scatter

more often, see Figure 3.3.

28

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Figure 3.3: Reduced volume due to ﬁnite particle radius. In previous implementa-
tions, the testparticles were pointlike (left), which left more free space between the
particle. The current implementation (right) takes into account the volume taken up
by the testparticles.

However, the scattering probability is lowered again by another effect: the particles
are screening each other. A particle might not be available for scattering with another
particle because there might be a third particle in between, an effect that is illustrated

in Figure 3.4.

Through this effect, a portion AS of the total 2—dimensional surface projection
S = 7r(2a)2 = 47ra2 of the “spheres of inﬂuence” of the particles is not available for
scattering, which lowers the collision frequency by a factor (1 — A—S/S), where A—S is
the average screened surface area. With :1: being the distance of the center of spheres
1 and 3 in Fig. 3.4. Allowing the distance to vary, within a range (as: + dzr), the
average number of “No. 3 - spheres” is 47rgzczd1r, each of them screening a surface
ring area 27ra(a — :1: / 2) of sphere 1. Therefore, the average screened area is

_ 2a 11 2 5
AS =/ dx 47rgm2 - 27m (a — g) 2 ¥ . (3.12)

Recalling the change in collision frequency, this effect leads to a reduction of the

29

3 5 . 3
J
1 .
2 1 2
Figure 3.4: Screening effect due to the ﬁnite volume of the testparticles. On the left,
the pointlike particles 1 and 2 could scatter with each other, particle 3 is not “in the

way.” On the right, due to the ﬁnite volume, they are prevented from scattering,
particle 1 would scatter with particle 3 ﬁrst.

scattering probability by a factor of

<1 — g) = (1 — 9- EMS) . (3.13)

Again, the product g - a3 is independent of the number of testparticles per nucleon

N. Including this factor, the modiﬁed scattering probability is

 

W’ = YEW , (3.14)
where
1— g- ﬂm3 1—11bEg/8 2
E _ 12 __ E _ 3
y _ 1_ Q. §Wa3 __ 1* 2bEg , b — 37m , (3.15)

b5 being the second virial coefﬁcient as yielded by the Enskog Theory of the dense

hard-spheres ﬂuid [Res77]; see ﬁgure 3.5.

The implementation of the modiﬁcations makes the second virial of the hard part
of the interaction non-vanishing; bE is positive and therefore leads to an increase in
preasure. This is partly compensated by the negative virial b5 that is due to the soft

(mean ﬁeld) part of the interaction; see for example [Dan96]. The equation of state

 

 

 

 

 

 

_ I r T j I I ‘1' IV 1' 1" 4
4 L a
J
3 — —<
m

>4 i
2 C— _
1 ; —
> I l I I I I l I I I I l I I I I I I l I I I I I

0 20 40 60 BO 100

Average Cross Section (mb)

Figure 3.5: The scattering enhancement factor YE as a function of the average cross
section. Shown is YE for g = 1, 2, 3m.

deviates from that of an ideal gas by both contributions, i.e.,
P = air (1 + 9(bE + as) + . . .) . (3.16)

One should note at this point that the Enskog Theory is non-relativistic; both the
advection vector (1 and the excluded volume, therefore also YE, are calculated in a

frame-dependent way.
3.3 Results and Conclusion

Our numerical calculation is based on the MSU BUU-code by Bauer et a1. [Bau86]
which was modiﬁed from a geometrical formulation of the collision term to a stochastic
formulation according to Ref. [Lan93]. Only N N collisions were taken into account,
which for the energies considered turned out to be a justiﬁed approximation. A full-
ensemble and a parallel—ensemble implementation proved to have similar results; re-

sults for different systems were compared with both Refs. [Bau86] and Ref. [Dan95];

they were found in satisfactory agreement. An interesting side-result at this stage

 

31

however was that changing the algorithm from geometrical to stochastic scattering
increased the frame—of-reference dependence of the result: when running the simu—
lation Within the lab-frame we found an asymmetry in the ﬂow which we attribute
to the fact that in this frame the projectile is Lorentz-contracted and the target is
not. Therefore within a spatial box inside of the overlap—zone of the two nuclei there
are many more testparticles originating from the projectile than from the target, re-
sulting in an asymmetry of the respective scattering rates. We are currently trying
to overcome this problem by y-dependent modiﬁcations of the scattering probabil-
ities, however, so far with only little success. The ﬂow—asymmetry vanishes when
the calculation is performed in the c.m.-frame, which is what we did in this work.
The geometry-based code did not appear to be sensitive to this asymmetry, however,
Lorentz invariance certainly still is an issue, see Chapter 4. Finally the full—ensemble
version of the stochastic code was modiﬁed according to Ref. [Ale95], and the scaling

of the collision rate and the ﬂow with the number of testparticles N was checked.

We simulated an (Au,Au)—collision at projectile energies of 250 and 400 MeV
per nucleon, and b : 3 fm over a total time of 70 fm/c. These energies are large
enough so that the repulsive (hard) part of the interaction dominates, but, as already
pointed out, low enough so that inelastic N N scattering can still be neglected. The
timestep length was 0.1 fm/c, the volume V was approximately 1.8 fm3, the number
of testparticles per nucleon N was 250, m in equation (3.1) was chosen to be n2,
and we used a soft momentum—dependent equation of state. Figure 3.6 shows the

evolution of the conﬁguration in the reaction plane.

32

 

m
o

—25o MeV

 

 

400 Mév

 

 

 

 

 

 

 

 

 

 

 

 

 

Figure 3.6: The nucleon conﬁguration in the reaction plane at different times during
the collision, shown is the nuclear density 9(m, y = 0, z).

 

 

 

33

 

 

    

 

 

 

 

 

 

   

60:"“l""l""|"": ."r'TﬁﬁHF'H'IH'i
so;— 250 MeV —; 60-. 400 MeV _
E 40:— —. I
\ : : -
8 30} —j 40.— ‘
43 : I
S 20:— i -
U : 20f —
10.— i I
0'. .. ....i.... 1. 0.11 . A...
O 20 4O 60 O 20 4O 60
t (fm/c) t (fm/c)

Figure 3.7: Collision rate for a 250 and 400 MeV (Au,Au) collision with b = 3 fm
versus time.

The result for the unmodiﬁed algorithm is shown on the left, the result for the
modiﬁed algorithm on the right, respectively. As expected, the nuclei disintegrate
slightly more violently due to the additional advection, resulting in a lower nuclear
density. A calculation with the additional advection alone revealed that therefore
also the collision rate decreases. However, this is partly compensated by the modiﬁed
scattering probability Eq. (3.15). On the other hand, a calculation with the modiﬁed
scattering probability alone shows a strongly enhanced collision rate, as expected.
The left panel of ﬁgure 3.7 illustrates this for the 250 MeV collision, the right panel
shows the collision rates for the 400 MeV collision. The solid curve refers to the
unmodiﬁed simulation, the dashed curve to the modiﬁed one. For the 250 MeV
collision two dotted lines were added, the upper line refers to a calculation where only
the scattering enhancement was taken into account, the lower line to a calculation

that only incorporated the additional advection.

Figure 3.8 shows the average ﬁnal transverse momentum versus the reduced rapid-

 

 

 

 

34
9 ,,,,, _ 11111
§ r(Au,Au), b=3 fm, 250 MeV ,E ' '(Au,Au), b=3 fm, 400 MeV m, ,nr‘
v . . . .... ,_
e 100— “i“ilooe .. —
:1 t zo ' l’ /
+3 M
e ‘ u ‘ i
g _ i
2 o 0
vi ' / e
g ... .
.. . _
E—l ' O ‘ " II ‘
go —1001§:I! i—IOO" mg —
H — ’Emﬂx _
d.) » l _
> f
< 1....11..l l ..|.. .[....l.. | l

 

 

 

 

 

 

—l.O —0.5 0.0 0.5 1.0 —l.0 —0.5 0.0 0.5 1.0
Reduced Rapidity (COM) Reduced Rapidity (COM)

Figure 3.8: Average ﬁnal transverse momentum versus reduced rapidity of the protons
in a 250 MeV (left panel) and 400 MeV (right panel) (Au,Au) collision with b = 3

fm.

ity as an indicator for nuclear ﬂow. The reduced rapidity is the rapidity ch devided
by the rapidity of the beam, which for the 250 MeV collision is approximately 0.36,
and for the 400 MeV collision approximately 0.45. From the slope of a linear ﬁt
around the origin, one can determine the nuclear ﬂow. The circles and the solid ﬁt

refer to the unmodiﬁed, the stars and the dashed ﬁt to the modiﬁed algorithm.

For the 250 MeV collision, at b = 3 fm the flow is z 147 MeV/(c-Unit of Red.
Rap.) for the unmodiﬁed, and z 170 MeV/(c-Unit of Red. Rap.) for the modiﬁed
algorithm. A calculation that only took into account the scattering enhancement
resulted in z 166 MeV/(c-Unit of Red. Rap), a calculation that only incorporated
the additional advection in z 159 MeV/(c-Unit of Red. Rap). For the 400 MeV
collision, at b = 3 fm the ﬂow for the unmodiﬁed algorithm is z 185 MeV/(c-Unit of

Red. Rap); for the modiﬁed algorithm it is z 195 MeV/(c-Unit of Red. Rap).

One concludes that in this speciﬁc simulation the ﬂow increases by approximately

16% for the 250 MeV collision, and 5% for the 400 MeV collision, with the introduction

of the new algorithm. Introducing only the additional advection, as already pointed

 

35

Table 3.1: Impact Parameter Averaged Analysis (b S 5 ﬁn) and Experimental Data

 

Nuclear Flow
Energy Unmodiﬁed Modiﬁed Plastic Ball[Gus88] EOS[Par95]

 

 

(MeV) (MeV/(c-Unit of Reduced Rapidity))
250 132 148 130 119
400 166 185 169 151

 

out, the collision rate slightly decreases, however, the nuclear ﬂow increases by about
8% for the 250 MeV collision. With the introduction of the modiﬁed scattering
probability alone, the ﬂow increases by about 13%. The fact that the contributions
from both modiﬁcations do not “add up” indicates that a perturbative approach to

their representation would not be justiﬁed.

An impact parameter averaged analysis for 250 MeV collisions with b g 5 fm
resulted in approximately 132 MeV/(c-Unit of Reduced Rapidity) nuclear ﬂow for the
unmodiﬁed, and 148 MeV/(c-Unit of Reduced Rapidity) for the modiﬁed algorithm
(12% increase). The Plastic Ball data indicate approximately 130 MeV/(c-Unit of
Red. Rap.) [Gus88] nuclear ﬂow, the EOS data only 119 MeV/(c-Unit of Red. Rap.)
[Par95]. For 400 MeV collisions the same calculations resulted in in approximately 166
MeV/(c-Unit of Red. Rap.) nuclear ﬂow for the unmodiﬁed, and 185 MeV/(c-Unit
of Red. Rap.) for the modiﬁed algorithm (11% increase; Plastic Ball: % 169 MeV/
(c-Unit of Red. Rap.) [Gus88]; E08: 3 151 MeV/(c-Unit of Red. Rap.) [Par95]).

Table 3.1 summarizes these results.

Overall we found the effect of the additional advection and the modiﬁed scatter-

 

 

36

ing probability to be signiﬁcant, but not crucial. Their implementation moves the
outcome of the simulations away from the experimental results. This indicates th<
need for an in-medium reduction of the N N cross section. This type of reduction was
ﬁrst found to be needed in studies of the disappearance of ﬂow [Wes86] and later alsr
in theoretical studies based on thermodynamic T-matrix theory at ﬁnite temperaturr
[Alm95]. These results were obtained by algorithm with closest approach techniques
If one wishes to address the question of the nuclear equation of state with a DSMC
algorithm, however, the corrections discussed in the present paper should be taker

into account.

3.4 Outlook

In this chapter, the virial corrections to heavy ion collisions in analogy with the theory

of gases have been studied numerically.

Recently, Morawetz [Mor97] introduced a more complete set of effective shifts that
represent mean values of various non—localities of the scattering integral. These shift:
enter the scattering integral in a form known from the theory of gases, however, th<
set of shifts is larger than the one intuitively expected. “Particle diameters” an(
other non-localities of the scattering integral are given in form of derivatives of tilt
phase shift in binary collisions. It can be shown that the shifts reduce to a “particlt
diameter” d for the classical limit of hard spheres collisions. In the case of nuclea

matter, the three distinguished space shifts maintain.

The derived nonlocal shifts are possible to be interpreted within a classical picturv

of scattering:

0 As already indicated in this chapter, one can imagine that a hard sphere shif

Ad has to be included in order to simulate hard sphere behaviour of the collidin;

 

37

particles, because the two particles cannot approach each other nearer than the

hard core distance Ad.

Next one has to account for the delay of asymptotic scattering with respect

to the real trajectory of two colliding particles. Danielewicz and Pratt have
recently discussed these time delays for equilibrium processes [Dan96]. This
leads us to the picture of two particle sticking together and moving as a molecule

with their center of mass velocity for a time A).

The two particles are forming a molecule which can rotate during their corre-

lated motion. This gives rise to an additional rotational shift A).

An extension of the above approach to incorporate the additional shifts is desirable.

Chapter 4

Causality Violations in Cascade
Models of Nuclear Collisions

4.1 Introduction

Experiments like the upcoming nucleus—nucleus collisions at Brookhaven’s Relativistic
Heavy Ion Collider (RHIC) will without doubt open new possibilities to study the
properties of nuclear matter at extreme pressures and temperatures. One major
focus of the research is the question of whether or not there is a phase transition
between hadronic and quark-gluon matter. The major open question that comes
with this undertaking is: If there is a quark—gluon plasma, what will it look like in
the detectors? What are its experimental signatures? The latter question can only be
answered reliably if a theoretical model for the collision processes that goes beyond
a phenomenological description is in place. One possible microscopic approach is to
extend the semiclassical transport theory (as for example described in Chapter 3) to

high-energy physics[Ka193, Sor89.1, Sor89.2, Ge192.1, Gei94.1, Gei92.2, Gei93].

Simulations of ultrarelativistic heavy ion collisions inevitably involve both soft and
hard processes. Below some energy scale, soft or nonperturbative phenomena neces-
sitate phenomenological description. Uncertainty is duly noted but for the foreseeable
future a rigorous theory for soft physics is beyond understanding. However, through-

38

39

out much of the collision hard processes dominate which are described by elementar
interactions between quarks, antiquarks and gluons (partons) as essentially semiclassi
cal particles. After parton initialization according to some reasonably chosen nucleo
structure functions, spacetime propagation is accomplished by discretizing time int
units Atstep and updating phase space densities according to relativistic transpor

equations including a crucial collision term.

Scattering processes in this theoretical framework are assumed to be non-retarde
which, on a microscopic level, leads to information transport with velocities that ca
approach ﬁ/Atstep, where o is the parton-parton cross section. That is, it can in
crease without reasonable bound. As energies for the individual interactions decreas
the corresponding cross sections increase. The timestep Atstep is for reasons of con
vergence chosen to be less than the parton mean free paths divided by their velocities
of the order a few thousandths of a fm/c. Two problems occur immediately: On
macroscopic level a series of subsequent causality violating interactions can lead t
shock-waves propagating faster than the speed of light. This is clearly unphysica
On a microscopic level the time-ordering of the incoming and outgoing partons of
scattering process becomes frame-dependent and Lorentz covariance is lost — this als
is unphysical. These two problems are especially serious in ultrarelativistic part0
cascades since collision rates are high, and the mean free paths approach the inter
particle distance. Problems arise from describing quantum-dynamical processes in
semi-classical picture, and from the demand of Lorentz—invariance in an equal-time
character simulation. Where is an acceptable compromise? And what does acceptabl

mean?

This chapter is organized in the following way: In Sect. 4.2 the origin of super
luminous information transport on a macroscopic scale is described in some deta

and some ﬁrst steps to eliminate it are being suggested. Sect. 4.3 moves towar

40

microscopic physics and provides mathematical details for origins of unphysical ef-
fects. Whenever cross sections are ﬁnite and action-at-a—distance inﬂuences particle
trajectories these problems inevitably arise. Sect. 4.4 deals with the resulting frame-
of—reference dependence of the simulation. In Sect. 4.5 the used version of a parton
cascade implementation is described. Then in Sect. 4.6 several different schemes
which reduce or eliminate superluminous transport are compared. Methods include
simply blocking collisions or truncations, scaling the cross sections downward while
increasing the number of particles, and suppressing low-energy collisions. We also
consider so-called wee-partons as a different way to deﬁne the initial conditions of
the simulation, and point out how much this scheme affects causality violating mech-
anisms. Finally, Sect.4.7 concludes by brieﬂy summarizing, and by discussing the

outlook for future studies.

4.2 Macroscopic Causality Violations

Superluminous macrosc0pic information transport occurs mainly in the transverse
(perpendicular to the beam) direction. The signal can travel over the diameter of
the cross section a in a single timestep and then continue this propagation from
timestep to timestep. Transverse signal velocities can therefore reach ﬁ/Atstep over
several timesteps, one is reminded of a chain of falling dominoes. The situation is
depicted in Figure 4.1 for such a transport. As part of the model the two initially
scattered particles cannot scatter again in this timestep. However, one timestep later
that scattered particle from the right scatters with another particle coming from the
right. As a result information from the ﬁrst scattering, such as particle momenta and
type, has travelled to the second one. This is a general problem of all transport
codes, and is worsened because the gluon-gluon cross section becomes rather large for

low energies (instead of vanishing). Along with the large gluonic cross sections are

41

 

 

 

V Signal
' . _ N CAtstep
\ V .
c { '--~
7r ' .

 

Figure 4.1: One particle coming from the left scatters with another particle comir
from the right. The information from this event is rapidly passed on through a cha
of subsequent scattering events.

relatively large gluon densities which result in very high probabilities for scatterir
in subsequent timesteps. Therefore information transport could be supported OVI

relatively large distances without damping.

In many existing cascade codes ﬁrst steps for dealing with the problem of supe
luminous signals have been taken: Most of the codes, including the one presentei
allow only one interaction per particle per timestep. This restriction is conseque1
and clearly justiﬁed since a timestep is the shortest scale in the model. This restrictic

prevents signals from avalanching over huge volumes within only one timestep.

A second restriction implemented in many codes is the “closest approach” crit
rion. For a scattering process it demands in addition to “spatial distance within tot
cross section” also “the two particles have reached their point of closest approai
assuming their current trajectories”. Usually one would look at |r“:r,,| with 33,, beii
the four—vector distance between the two particles involved, and demand that th
quantity is minimal. As the two particles’ positions are taken in the same timesti
this quantity reduces to the spatial distance squared. Even though this looks li

a Lorentz-invariant criterion, it is not, since the conceptual necessity of taking bo

42

particles’ position at the beginning of a timestep (i.e. with vanishing time-separation
in the lab-frame) is a non-Lorentz invariant restriction, see section 4.4. Nevertheless,
the “closest approach” restriction prevents causality violating signal transport in the
longitudinal (parallel to the beam) direction, but unfortunately has no effect on the

transverse direction.

To demonstrate those causality violating shock-waves, the parton cascade simu-
lation of a 100 GeV/nucleon (p,Au)-collision was run with different constant cross
sections and S—wave scattering, the distance of of scattering events from the beam
axis versus the simulation time is shown in Figure 4.2. The solid line indicates the
distance that could be reached by a signal that travels with the speed of light. One
can clearly see the outwards travelling shock-wave and how it gets damped out with
time. For the smaller cross sections the causality violation is rather small and only
present in the initial stages of the interaction. One should note that for realistic
energy-dependent cross sections the initial scattering events take place with a lower

cross section than the later ones because the c.m.-energy is much higher.

In the lower right panel the maximum distance dmax from the beam axis at which
a scattering event occured during the simulation is plotted versus the cross section 0.
The horizontal line indicates a distance that during the simulation time could only

be reached with the speed of light.

 

in1|I..rerruulvvvulnuvulvviv

4 0:0.1 fmg

’

 

W.l .‘.1...-.l..'.

0’ . .
0.0 0.5 1.0 1.5 2.0 2.5
time (fm/c)

.1
i

1

.l.

 

n

3.0

 

lllil!VVV[VVVVIVIIYIIVIIIVII

0.0 0.5 1.0 1.5 2.0 2.5 3.0

time (fm/c)

  

43

 

 

 

IllvllililiilllIlYT‘l'liiiKill

4 0:0.2 fmz —

 

1.
time (fm/C)

vavuluuuulviriliuuv'vvv’rlr‘l

/_

 

T23 fm/c_
11111
0.0 0.1 0.2 0.3 0.4 0.5
Cross Section (fmg)

|..;

 

Figure 4.2: Distance d of scattering events versus simulation time at constant cros.

sections of o 20.1, 0.2 and 0.5 fmz.

 

.

44

While during the initial stages of the collisions the outmost scattering events
occur at distances larger than those allowed by causality arguments, the information
transport soon appears to be damped. It turns out that this phenomenon is not — as
one would expect — mainly due to a dropping of the collision rate, but rather to a form
of random-walk: only a fraction of the individual signal propagations lead outwards,
others have the opposite effect, which results in an effective damping. In fact it
turns out that the expectation value of the distance of scattering events from the
beam axis is roughly proportional to the square—root of the simulation time, which
is a characteristic feature of random-walk mechanisms, compare the top panels of
Figure 4.2. This mechanism is a valid description until the cross sections are so large
that the parton distribution basically appears solid, in that limit the information
travels outwards proportional to the simulation time, see the panel for a = 0.5 fmz.
By the nature of random-walks, the information expands fastest in the initial phase
of the collision. With realistic energy—dependent cross sections however, those initial
scattering events happen at lower cross sections than the ones in the later stages due
to the higher c.m.-energy, which partly suppresses this initial outburst. Overall for
020.2 fm2 the causality violations are still rather moderate, for 0.5 fm2 the effect

becomes dominant.

As already pointed out, Figure 4.2 also shows the maximum distance of scattering
events from the beam axis as they occurred for different cross sections. As it turns
out, this value depends linearly on the cross section for an extended range, beyond
a z 0.4 fm2 the outmost scattering events occur at a distance that can only be

reached with superluminous signal velocities.

45

4.3 Microscopic Causality Violations

In addition to macroscopic causality violations such as superluminous shock-waves,
both the non-retarded interactions and the model inherent propagation of the whole
particle conﬁguration from one timestep to the next lead to causality violations on

the level of elementary scattering processes.

A Lorentz-invariant simulation of parton scattering would require a truly 4—dimen—
sional conﬁguration space, in which an interaction can be established between any
space-time point of one parton’s trajectory and any other space—time point on another
parton’s trajectory. In the framework of classical ﬁelds, the scattering of two particles
would be represented by their continuous change of trajectory in the retarded ﬁeld of

the respective other particle.

The actual realization of scattering events in the simulation only agrees in a certain
limit with this model, namely if the two partons are coming from inﬁnity with inﬁnite
rapidity and opposite directions: with increasing rapidity, in the lab-frame the ﬁelds
are distorted to pancakes with an orientation perpendicular to the trajectory, reducing
the dominant part of the interaction between the particles to a shorter and shorter
region along their trajectories. In the limit of inﬁnite rapidity, the interaction between
two partons travelling with a non-vanishing impact parameter in opposite directions
would be reduced to a momentary interaction at the point of their closest approach,
making them change their trajectories at equal-times. This interaction would take
place at a spacelike distance equal to the impact parameter , why does this not violate
causality? The reason is, that the parton’s ﬁeld—pancake e even though it is travelling
along with its present position — is build up from contributions out of the parton’s
past, the respective other parton does not scatter from the ﬁeld generated at the ﬁrst

parton’s present position, but from a ﬁeld component that was generated along the

46

 

Figure 4.3: Lightcone and 3—dimensional hyperplane (2-dimensional illustration)

ﬁrst parton’s trajectory at a lightlike distance.

This limit agrees both with the claim of Refs. [Sor89.1, Sor89.2] that the interac-
tion distance should be spacelike, because otherwise an interaction would inﬂuence
the absolute past of one of the partons, and with the actual realization of scattering
processes in the simulation: in the framework of a 3+1 dimensional transport simula—
tion the distance of the testparticles is by the model itself determined to be spacelike:
The information available at the beginning of a timestep is no more than the points
of the trajectories on a certain any, z-hyperplane, as well as the corresponding mo—
mentum 4—vectors of the particles, also only at that time-coordinate. The lightcone
of any one particle extends both before and after that hyperplane (see Fig. 4.3), but
has no extension into that hyperplane itself, causing all other particles present in the

simulation to have a spacelike distance from it.

But this limit soon loses its applicability: Staying in the model of interactions
between the partons due to classical ﬁelds, as already pointed out, the ﬁeld trav-

elling with the partons is in fact build up from contributions out of their history,

which for partons coming from inﬁnity results in Lorentz-contracted ﬁeld-pancakes.

 

 

47

But in a cascade simulation it is highly unlikely that seen from a parton A, parton
B already was on the same trajectory a lightlike distance ago. In fact, as the par-
tons are nearly travelling with the speed of light themselves, spacetime points with
lightlike distances along their trajectories can happen to be very long ago, which in
connection with the high interaction rates makes the idealized picture of the ﬁeld
pancakes inapplicable. Changes in trajectories can make more than one point along
the trajectory of parton B have a lightlike distance to parton A, the acceleration
connected to the changes in trajectory of parton B generates additional ﬁelds. Even
worse, through the mechanism of parton generation and absorption, parton B might
not even have existed at lightlike distances from parton A, still within the framework

of the transport simulation, scattering would be possible — and causality violating.

Also, as in the idealized picture the c.m.-frame is moving parallel to the partons’
trajectories and perpendicular to their distance at the point of closest approach, the
equal-time character is true in both the lab- and the c.m.-frame. As soon as the two
partons are not travelling into opposite directions, their spacial distance also leads to

a time-separation, making the scattering time different in the lab- and c.m.-frame.

4.4 Frame-Dependencies and Time-Ordering

An obvious consequence of the model’s causality violations is its frame-of-reference
dependence: The simulation will for example certainly lead to different results when
run in the rest-frame of the target. This was already shown for the example of a
12C'+12(3' internuclear cascade calculation in Ref. [Kod84], the authors demonstrated
that both the total number of collisions and the individual time-ordering of given
collisions were frame-dependent. As obvious as that seems from the simple problems

pointed out in sections 4.2 and 4.3, for a parton cascade it is not yet the whole story:

48

The simulation cannot even be run in the rest-frame of any one nucleus, because there
is no possibility to create appropriate initial conditions. The parton distribution in
nucleonic matter is frame-dependent, the higher the energy of a nucleon the higher
the number of virtual partons. To put it in other words, usually components of the
nucleonic wavefunction are considered a parton if they carry a momentum fraction
beyond a threshold given by a minimum art-value, the latter one being a parameter
of the model. However, by boosting the nucleonic wavefunction into an accelerated
frame of reference, a procedure described by the model—parameter Q2 of the parton
distribution function f (:3, Q2), more and more components of it cross that threshold

m-value, leading to higher and higher parton numbers.

Therefore, when asking questions about how the collision of two partons looks in
their rest-frame as compared to the lab-frame, another valid question is: Do those
two partons even exist in the other frame? Or are in that frame other partons around

that the collision partners would be much more likely to scatter with instead?

However, as this question cannot be addressed within the framework of a transport
simulation, it is worthwhile to study the effects of a transformation of two given

partons into their c.m.-frame.

4.4.1 The Transformation Equations

The transformation of a four vector ((10, a) with a boosting vector ,6 and '7 2: V1]?

 

is given by the following set of equations:

 

(10m = 7(001ab — Balab) (4-1)
acm = alab + 75—21(16alab)13 _ VIBGOlab (42)
aOlab : 7(a0cm + ,Ba'cm) (4'3)
alab = acm + %(ﬁacm)ﬁ + 7ﬁ00cm (44)

 

 

49

4.4.2 Finding the Boost Parameters

The c.m.-frame is deﬁned by
pem,1 : ~pcm,2 ’ (4'5)
pm],l and pm,2 being the momenta of the partons in the c.m.—frame. This leads to

._ P1ab,1 + P1ab,2

.-~ . 4.6
Elab,1 + Elab,2 ( )

4.4.3 Transformation of Momenta and Energy

Using equation (4.1) with ,3 from equation (4.6) gives the energies of the partons in

the c.m.-frame:

E1cm,1 : 7(Elab,1_ﬂplab,1) (4'7)

Ecm,2 = ”7(E1ab,2 “ 511mm) - (4-8)

The relative momentum can be calculated by means of equation (4.2) with both plab,1

and plam, however, parton 1 is chosen:

”y — 1
pcm :2 plab,1 + 7(ﬁplab,l)ﬂ _ VIBEIabJ (4-9)

The total energy \/E in the c.m.-frame is given by means of equation (4.3) and the
four vector of the total momentum in the c.m.—frame (Em), + Ecmg, pom1 + pcm,2) 2:
(Ecm, 0):

_ Elab,l + Elab,2

\/E :2 Ecm — (4.10)

4.4.4 Transformation of the Relative Positions

The spatial distance between the two partons is deﬁned as the vector from parton 2

to parton 1. The time separation is to be assumed between zero and the timestep

L:

 

 

50
length:
Arlab I= 7'1ab,1 — ’f‘iab,2 (4-11)
0 S Atlab S Atstep (412)

If however one assumes the particles’ positions to be taken at a certain moment inside
of the timestep interval, for example at the beginning, the time separation Atlab is

zero .

The distance vector in the c.m.-frame is then given by

 

AT‘cm = Arlab + ’7/6—2 1(ﬁA7‘1abM'3 - (4-13)

In addition to that, the partons will now have a time separation in the c.m.-frame:

Atcm : _716Arlab (414)

As soon as this time-separation is larger than the transformed timestep length,
Lorentz-invariance is clearly violated: The partons’ positions are not taken within

the same timestep inside the c.m.-frame anymore.

4.4.5 Ways to Calculate the Positions of the Partons in the
Next Timestep

Propagation Within the c.m.-Frame

A ﬁrst idea is to propagate the partons in their c.m.-frame. The distance vector then
changes to

Ar’cm : Arm, + (15$ + Pc_m> Atstepycm . (4.15)

cm,1 Ecm ,2

In order to have the time Atstepycm correspond to a four vector with the time compo-
nent Atstep in the lab frame, one has to choose

At 6
Atsteppm : _,:l_t_p , (4.16)

 

51

and therefore

Arém 2 Arm + ($7 + %> Atstep . (4.17)
cm, cm,

This strategy however completely neglects the spatial components of the timestep

that occur during the reverse transformation.
Propagation Within the Lab—Frame

A more reﬁned strategy is the propagation of the partons in the lab-frame and subse—
quent transformation of the new positions to the c.m.—frame. The propagation vectors

for the partons take the form
P i
(Atstep: ﬁgAtstep> -
With Aao(a) being the spatial Lorentz-transformation into the c.m.-frame for a

vector (a0, a), one gets

I plab,1
Arcm = A0+Atstep—0—Atstep (rlabd + E

 

 

P
Atstep _ rlab,2 _ Elab,2 Atstep) (418)
ab,1 lab,2

 

 

_ Arm, + A0 <Plab,1 Atstep _ Piab,2 Atstep)
Elab,1 ElabZ

Y

— A"hem ‘1' A0 (‘ZﬂAtstep> _ A0 (MAtstep> -
lab,1 17,2

 

Ela
Because of
A0(plab,i) : AElab,i(plab,i) + VﬁElabﬁ : pcm,i + VIBElabJ (419)
the result is
P P
Ar'm = Arc + < m + m > At, 6 . 4.20
C m 7(Ecm,1 + ﬁpcm) 7(Ecm2 ‘- Ichm) t P ( )

The following Minkowski plots again show the above in one dimension. The dashed

lines denote the world lines of the two partons, in one dimension their closest approach

 

 

 

 

52

can only be an intersection. The two solid lines denote the coordinate system in the

partons’ c.m.—frame.

 

 

 

 

 

 

 

 

53
zlab [fm]
.10 rP .03 .00
o o O c_>
'4 r '3 a
O O
| | I I z
00 03 CO 00 ,_a.
0° 5
x I \ I :5:
\\ \ 2
to \ \ U)
C? \ W
l—¥ _ 1—1-
0|
w E
O
:s <-+-
O
s. N _ a
O
. E
.. 2
75 o
r.535 B m
...C.’ m (D
99 E; _ m
i ,1, ﬁ
g“ >
\oo “C
&o "d
c. _ r1
0| 0
w 5%
:3”
H
'<.> \ ~-
3— \ \T D
I "* \
N o \ |—‘
n E, \\ 9
in H' '
t;_ E \ _ B
0| o \ (D
N H \ :5
I \ 3
O
:3

Figure 4.4: Minkowski Plot. Highlighted is the timestep in which the particles are at
the point of closest approach. Note the time-order of the four events “Particle 1(2)
at the beginning(end) of the timestep” in the c.m.-frame. These points in time are
given by the intersection of the c.m.-t axis with the dot-dashed lines that are parallel
to the c.m.-z axis.

54

 

 

 

 

 

 

 

 

2 lab [fm]
.N 51‘ .03 .00
C? Q 0 C?
H H I—x H
o o o o
| I I I g
to co on no ,_..
O D
O
\
\ \ l W
O
\ \
\ 2
_N \ \ U)
C? \ w
H __ \ \ _ l-"
l O \ \
l \ "U
to \ ._.
\ \ O
11> \ \ ri-
C.’ \ \
... _ \\ N — 9s
0I \ °
to \\ B Q
\. \\ '2 O
f"' g \'\ \ \ a U)
n—a - \ \ \ (D
a s - w
h 3,, ﬁ
...,
a >
\oo "C
&Q "d
... ’1
CI 0
w g
_.. tr
.0 ....
H
O D
(L, I—x
,_. U
19 H'
H B
C; (D
N b
9.
O
D

Figure 4.5: Minkowski Plot. Highlighted is the second timestep after the particles
had reached their point of closest approach. Note the the time-order of the four events
“Particle 1(2) at the beginning(end) of the timestep” in the c.m.—frame has changed.

 

 

I,
I 55

‘, The time step length in these examples is 2 - 10‘3 fm. The dot—dashed lines run
1 through simultaneous time-space points in the respective frames. In ﬁgure 4.4 the
time step with the point of closest approach is selected, one can clearly see the time
separation of the two partons both at the beginning and at the end of the time
step — this separation is rather small compared to the transformed time step length.
However, two time steps (see Fig. 4.5) later those two time separations overlap, in the
c.m.-frame the faster parton has reached it’s position at the end of the lab frame’s
time step before the slower parton has reached it’s position at the beginning of the lab
frame’s time step. The smaller the lab frame’s time step the more likely this happens,
which is counterintuitive to the usual fact that a simulation gets better with smaller
timestep lengthes. In one dimension this reversed time-ordering is impossible in the
time step of closest approach. In two or three dimensions this conﬁguration can occur
c.a.

as soon as Z c at the point of the closest approach distance [Arc'a' in other

Cm—

words, in two or three dimensions the reversed time-ordering is a problem even at the

critical timestep of closest approach.
4.4.6 Consequences

As it had turned out, by boosting the partons with ,8 into their c.m.—system, there

they have a time—separation of
Atom : “VIB'A’rlab ' (4-21)

The partons’ positions, being taken at equal time in the lab—frame, will, unless ﬂ
and Ahab are perpendicular, not be at equal time in their c.m.—frame and vice versa.
The construction of a “Lorentz-invariant CMS distance” as in Refs. [Sor89.1, Sor89.2]
seems doubtful. This time-separation can very well be larger than the transformed

timestep length, leading to a situation where the partons’ positions are not taken

 

 

 

 

 

56

pom. pc rn
Aram. Artm.
_pc.m. _pc.m.

Figure 4.6: Two partons passing their point of closest approach in the c.m.-frame.

within the same timestep within the c.m.-frame anymore. As a result, the closest
approach criterion can only provide a means to get an averaged point in time for an

actual scattering event.

For this criterion, described in section 4.2, one can choose to deﬁne it either in
the particles’ c.m.—frame or in the lab-frame (a third approach will be described in
subsection 4.6.6). Because of the non—vanishing time—separation of the partons in
their c.m.—frame, the two possibilities are not equivalent. Both methods had been
compared, but nevertheless no signiﬁcant change in the collision rate was found, even
though the individual collisions happening in both simulations are different. For
further simulations the c.m.-frame was chosen, since intuitively this frame seems to
be more signiﬁcant for the individual scattering events. In this frame, the scalar
products of the momentum and the distance vector both at the beginning and at
the end of the timestep are to be calculated. If this scalar product changes sign,
the position of closest approach is reached within this timestep. In Figure 4.6 Arm,
denotes the distance vector at the beginning of the timestep, Ar’cm at the end. The

criterion is (pcm - Arcm)(pcm ' Artm) S 0-

The non-vanishing c.m.-time-separation unfortunately also makes the calculation

of Arém ambiguous.

57

2!

A

 

 

 

 

—t—>

Lab—Frame c . m. -Fram9

Figure 4.7: Possible discrepancies in the time-ordering of scattering processes that
are usually prevented by the spacelike distance of the scattering partners.

In the model used, the c.m.-positions at the end of the timestep are calculated
by propagating the partons within the lab-frame and boosting the result into the

c.m.-frame. The result, expressed in c.m.-quantities, is

Ar’

Cm

: Arcm (4.22)

 

 

pcm pcm
+ + Ats e -
(fl/(Ecmﬁ + g ' pcm) 7(Ecm,2 _ :6 ' pom)) t P

Due to the spacelike distance of the partons however, the deﬁnition of “beginning”

and “end of a timestep” becomes frame-dependent.

Consider a process where the partons A and B scatter and produce a parton pair
C and D, see Figure 4.7. In the lab-frame this happens within one timestep, at
the beginning of the lab-frame’s timestep (denoted by thin vertical lines in the left
panel) partons A and B are present, at the end partons C and D. In the transformed
timestep within the c.m.-frame however (denoted by thin vertical lines in the right
panel), this set-up can easily be distorted to a situation where in the beginning of
that timestep partons B and C are present, and at the end partons C and D. Seen

from yet another frame, in the beginning A and D could exist.

One consequence is that in the simulation re-scattering processes have to be ex-

 

58

plicitly forbidden: In spite of the closest-approach criterion re-scattering would still
be possible if two partons are scattered towards each other — as the distance between
partons involved in a scattering event can very well be greater than Atstepc, it is
possible for those two scattered partons to again reach a point of closest approach in

a later timestep. In the case of a retarded interaction, this would not be possible.

More relevant then the time-scale given by the timestep, which is a model-depen-
dent parameter, is the time—seale given for example by the time it takes for the two
nuclei to cross through each other, which is about 0.2 fm/c. This time-interval corre—
sponds to an impact parameter of 0.2 fm beyond which the time-ordering of scattering
events is frame-dependent. As the interactions in question happen between particles
travelling into the same direction, those would also be low-energy interactions with
cross sections very well in the range of such impact parameters. The effects of this
frame-dependent time-ordering of the incoming and outgoing particles of a scattering

event have to be subject of extensive examination.

An interesting approach is that of Refs. [Lan93, Dan95]: In this method, con-
ﬁguration space is divided into boxes. Within the boxes the scattering partners are
randomly chosen according to a probability that is a function of the cross section —
convergence of this method towards the solution of the Boltzmann equation for an
inﬁnitely small box size and timestep length, as well as an inﬁnite number of test-
particles, is shown in Ref. [Bad89]. The code of Rev. [Dan95] was successfully made
more efficient by only considering a fraction of the possible pair combinations within
each box, and compensating the otherwise resulting loss in collision rate by an en-
hanced probability for the collision of the chosen pairs. The authors of Ref. [Lan93]
claim that this stochastic method of formulating the collision term is covariant since
it is dealing with transition rates instead of geometrical interpretations, therefore no

problems connected with the time-ordering of processes would occur. In fact, since in

59

the model the time—order of processes is chosen randomly anyway, the model has no
“right” time-order that could be distorted by relativistic effects. However, this new
approach does not overcome the problem of superluminous shock-waves. Within a
given cell of longitudinal size 67‘“ and perpendicular size 671 any two particles have a
chance of colliding. Since the subsequent advection step may carry some of these into
neighboring cells, the maximum transverse velocity for information transport will be
(Sn/Atstep. Since in general 67% is several times larger than the cross section a, the
resulting maximum possible causality violations in transverse direction are even larger
than in the method discussed above. In addition, the stochastic method also allows
for superluminous information propagation in longitudinal direction, with velocity

6T||/Atstep-

4.5 The Parton Model

Classical simulations, which have been successfully applied to heavy ion collisions
at intermediate energies [Ber88], are now [Ka193, Sor89.1, Sor89.2, Gei92.1, Gei94.1,
Gei92.2, Gei93] being extended to high energies. The main step in the extension of
this microscopic model is using a parton based picture of the nuclei rather than a
hadronic picture; consequently the interactions between the testparticles are to be

described in the framework of QCD, leading to so-called parton cascades.

The code works in 3+1 dimensions using fully relativistic kinematics for the par-
tons, where the quarks are consequently treated as massive particles. Both quarks
and gluons can be off—shell. The initial conditions are determined by standard parton
distribution functions f (2:, Q2) [Bot93], where the value for Q2 and the minimum :1:
are parameters of the model. Technically stability of the incoming nuclei is guaran-

teed both by a coherent motion of all partons in longitudinal direction, and by the

60

Figure 4.8: Typical contributions to the gluon-gluon cross section. The four point
diagram on the far right leads to divergencies in the cross section that cannot be
handled by phenomenological cut—off masses.

restriction that particles from the same nucleus cannot scatter with each other before

at least one of them has scattered with a particle from the other nucleus.

In this preliminary version of the code only QCD processes with two partons
in the incoming and two partons in the outgoing state are implemented [Mur93].
Phenomenological screening or cut-off masses have been added into the propagators
to avoid divergent total cross sections. However, the gluon—gluon scattering cross
section includes a four point diagram which does not contain a propagator, Fig. 4.8

shows typical contributions to the gluon-gluon cross section.

This means that the divergence in this cross section must be handled differently
from the other cross sections. One usually regularizes the gluon—gluon cross section
by setting it to a constant value below a certain cut-off energy, in this case a a cutoff
of s z 0.25 GeV2 has been chosen corresponding to a S 0.45 fm2. This cutoff seems
appropriate because it agrees both with a reasonable c.m.—energy below which Per-
turbative QCD loses validity, and with the limiting cross section for superluminous
shock-waves found in section 4.2. Other methods of cutting off these cross sections

while maintaining the correct physics are being investigated. Also, no medium mod-

 

61

100 I I I I l I I I I I I I I I I I I I I I

k "rig—>99

119—>719
uu—>uu

 

 

118—>115“
uﬁ—Iuﬂ
uﬁ—Igg
ggeuﬁ

It'll—>55

 

 

 

lo—B I I I I l I I I I l I I I 1 l I I I I l 1 1 I I
O 2 4 6 8 10

\/s

 

Figure 4.9: Example of partial cross sections being used versus the c.m.—energy.

iﬁcations to the elementary cross sections are taken into account. Figure 4.9 shows
the cross sections being used. The analytic expressions for cross sections have been
checked in the massless quark limit against results from the literature [Com77], some

of them also in the massive case [Eic84].

4.6 Comparison of Methods for Dealing with Su-
perluminous Signals

This section contains a comparison of different methods for dealing with the problem
of superluminous signal transport and is the outcome of simulating central collisions

of a 100 GeV proton on a 100 GeV per nucleon gold nucleus.

62

Effect of Timestep Length

 

 

 

 

25000 | l Illllll I I IIIIIII | I IIIIIII I I IIIIIII I I | Illll
I +____ :
20000— D ‘W--—»1\+ _
:3 : :
\

If: ' \ ‘
DD .. \\ _
E 15000— \ —I
0 _ \ a
3 - \ .
(2‘3 _ \+ a
u—u F ID\ ‘
0 10000— x —
t ' \ ‘
.0 — E _
E! _ o o o _
2 — x x x _
5000—— —
I x I
0"] I llIlIII I IIIIIIII I I IIIIIII I I IIIIIIl I llllllf‘

10—5 10‘4 10'3 10'2 10—1

Timestep Length [fm/c]

Figure 4.10: Total number of scattering events versus timestep length used for the
simulation. The dashed line interpolates the results from simulations of (La,La)
collisions, the dotted line from (Cu,Cu) collisions. The symbols denote the actual
data points generated in the respective simulations.

63

 

 

x (fm)

 

 

 

 

I.
.|....|.I..|....|....|..

—-3| —2 —l 0 l 2 I '25 —3J
2 (fm)

 

 

Figure 4.11: In the left panel the projection of the initial parton conﬁguration on
the x, z-plane is depicted. The highly Lorentz-contracted proton is moving to the left
into the equally contracted gold nucleus which is moving to the right. The right panel
shows the conﬁguration 1 fm/c later.

This set—up seems suitable since it allows observation of the information transport
from the incoming (comparably small) particle in nuclear matter. Table 4.1 summa-
rizes the parameters used for the simulations. The timestep length of 0.0002 fm/c was
determined by running the simulation with different timestep lengths and observing
the total number of collisions. Figure 4.10 shows the outcome of those simulations,
0.0002 fm/c can be identiﬁed as the longest timestep length possible without leaving
the plateau of nearly constant total collision numbers. Choosing the longest possi—
ble timestep length for which the simulation still converges minimizes the effect of

causality violations due to “instantaneous” reactions over a ﬁnite distance.

With the values chosen for Q2 and the minimum m—value approximately 9000
partons per testrun are generated. For the following calculations, the longitudinal
extent of the nuclei in the model is determined by the Lorentz—contraction only,

Figure 4.11 shows the initial conﬁguration in the lab-frame and 1 fm/c later.

The effects of “wee-partons” are examined in the last subsection. The cross sec-

64

Table 4.1: Parameters used for the comparisons

 

 

Parameter Chosen Value
Total time of simulation 3 fm/c
Timestep length 0.0002 fm/c
Number of timesteps 15000
Proton energy 100 GeV
Energy per nucleon Au 100 GeV
Impact parameter 0 fm
Minimum m-value for f (3:, Q2) 0.005

Q for f(:I:, Q2) 25 GeV
Bag—radius for nucleons 0.9 fm

as 0.2

Cut—off mass for gluon propagators 1.0 GeV
Cut—off mass for quark propagators 0.2 GeV
Number of parallel test runs 5

tions were energy-dependent according to section 4.5 and Figure 4.9.

The simulation was ﬁrst run for S—wave scattering in the expectation that this
would be the worst case. Compared to the realistic forward peaked angular distribu—
tions, S—wave scattering favors transverse signal propagation. However, it turned out,
that the angular distribution hardly makes any difference as far as causality violating
effects are concerned, the only noticeable change was in the rapidity distributions,
where it turned out that the gap in rapidity between unscattered and scattered par-
tons was wider for S-wave scattering: the partons lose more units of rapidity in
the initial collisions. The small effect of the angular distributions is understandable
from the fact that even isotropic distributions in the c.m.-frame are strongly forward

peaked in the lab-frame.

The top row of Figures 412,413 summarizes the results for a simulation with

65

a realistic angular distribution. The top left histogram of Figure 4.12 shows the
signal velocity distribution regarding directly subsequent collisions, see Figure 4.1. In
other words, with T5, being the position of the ith collision of the nth particle, the
distribution of

U501: 7') :2 ITnJ _ I"71,z'—lI/(tn,z' " tn,i—1) (423)

for all n and i > 1 is shown. The solid histogram shows the radial component of those
velocities, the dashed histogram includes the longitudinal component — the strong
peak at 0 results from processes for which the signal transport was dominated by the
partons’ motion with the respective nuclei. On the order of 50 isolated scattering

events with velocities from 500 up to 13600 have been suppressed in the plot.

66

Figure 4.12: Results of simulations with different mechanisms to suppress superlumi-
nous signal transport (Part 1).

67

VA/C

0.4 1.0

vS/c

0.4 1.0

Figure 4.12

1

O.

0.1

Wommwmbnm

Kong .9

EOQE 0

30%; U

gong H

 

4.0 10.0 40.0

 

 

 

 

 

 

 

 

 

 

 

0.4 1.0

 

 

4.0 10.0 40.0

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

4.0 10.0 40.0

0.4 1.0

4.0 10.0 40.0

vA/c

0.1

vS/c

0.1

68

Figure 4.13: Results of simulations with different mechanisms to suppress superlumi-
nous signal transport (Part 2).

Figure 4.13

<v>/c

<v>/c

<v>/c

<v>/c

<v>/c

<v>/c

 

Nurnber of Events
0 51015 20

69

 

I T' *“T'

10.0
5.0 -

 

 

 

 

 

 

 

 

Number of Events

 

 

Radial Distance Radial Distance Radial Distance Radial Distance
r N
U‘ C

\>;:§§§iéll\
| l .

Radial Distance

Radial Distance

tnne(hn/®
oo 05 10 15 20 25 ao
2.5 _ I I I I I _

 

 

 

...
01
|

I

 

 

.1"
UI
I
I

 

1.0 ﬁx I I I I
25 '

 

I

1.0

W
2.5 — i

1.0
2.5

2.0

 

 

 

1.0
0.0 0.5 1.0 1.5 2.0 2.5 3.0

tune(ﬁn/®

CI IGPOIAI O IGPOW H IQPOIAI V IQPOIAI 99H913J9H

ﬂ [9POIAI

70

The right panel shows a more general impression of the signal velocities. Instead
of just taking into account the signal velocity occurring within one single scattering
event, in this histogram the signal velocity from the very ﬁrst to the very last scat-
tering event that a parton was involved in is calculated. With the notation from
above and N (n) being the number of collisions that the parton n was involved in, the

distribution of
11,407.) :: Irn,N(n) _ rn,1I/(tn,N(n) _ tn,1) (4.24)

for all n is shown. Again, the solid histogram shows the radial component of the signal

velocities, the dashed histogram also takes into account the longitudinal components.

Because of damping effects the average velocities are much smaller than the Veloc-
ities from one event to the next, the left panel of Figure 4.13 shows how subsequent
scattering events are leading to this overall damping of peak velocities occurring in
single scattering events: again starting from the very ﬁrst scattering event of a parton,

the signal velocity to the ith following scattering event of the partons

’U(’I’L, Z) = Irmi - THAI/(117,”, — tn,1) (4.25)

is calculated, and the average (v)(z') over all particles n that had an 2th collision
is shown; again both the total and only the radial component of the velocities are

plotted.

Finally, the right panel shows how information is travelling outwards. The outer
boundary of its spreading is after a very fast expansion in the ﬁrst few tenth of a
fm/c travelling with approximately 0.170. The maximum distance of a scattering
event from the beam axis is 2.5 fm. A more detailed analysis of the overlay of the
average c.m.-energy and the distance of scattering events from the beam axis shows
that the expansion of a shock-wave begins rapidly immediately after the c.m.-energies

get smaller, but then is damped out — only the initial stages of the collision seem

71

to be problematic. However, one should note that in the later stages the individual
collisions happen with very high signal velocities, and that it is rather by the random-
walk character of the signal propagation over longer distances that the information

does not travel with the same high velocity, compare Sec. 4.2.

4.6.1 Model A: Restrictions on the Signal Velocity

An obvious way to avoid superluminous signal velocities is to systematically suppress
all scattering events that would lead to signal velocities faster than the speed of light,

the second row of Figures 412,413 summarizes the results of this simulation.

It is not surprising that the velocity restriction leads to a huge reduction in the
collision rate: While in the ﬁrst stages of the collision the collision rate drops to
fourtynine percent of the original rate, it drops in the later stages to only about eight
percent of the original rate, there is no outwards travelling information at all — it
would be unwise to apply such an unjustiﬁed prescription given the huge effect, even
more so, as from the extrapolation of intermediate energy results one would expect a

rather strong outgoing effect of the proton impact.

4.6.2 Model B: Suppressing of Low Energy Collisions

Once a parton has scattered with a parton from the other nucleus it can scatter with
fellow partons from the same nucleus at much lower c.m.-energies. Soft processes
between partons from the same nucleus are problematic since the elastic gluon-gluon
cross sections at low energies are large; systematically suppressing these reactions by
placing a c.m.-energy cut-off on the scattering events should also suppress high signal
velocities. For this model calculation the cross sections are set to zero below 0.4 GeV,
see the third row of Figures 412,4.13. The outcome of this simulation proves that

indeed the shock-wave travels outwards through those low energy collisions, only the

72

initial outburst remains.

In fact this model provides a way to address the problem of the energy-dependent
parton resolution described in the beginning of Sec. 4.4: Disregarding interactions
below a certain c.m.-energy cutoff could also be viewed as disregarding the interaction
partners. By smearing out the step-function of the cutoff to a smoother function that
determines the probability of scattering events in dependence of the c.m.-energy, one
could try to simulate the behavior of the parton distribution function in dependence
of the resolution parameter Q2. However, other physics would have to take the place
of the QCD interactions, in this lower-energy regime nucleonic interactions take over

and hadronization patterns must be applied.

4.6.3 Model C: First Steps Towards a Retarded Interaction

As already pointed out, the main reason for the superluminous signal velocities is
the instantenous character of the elementary interactions. To consequently overcome
this problem no theoretical framework with time delays has been developed so far for
parton systems, although some simple hadronic examples have been worked out. To
estimate the effects of retardation it nevertheless seems appropriate to introduce a
simple - but unphysical — way of time delay in the interactions: After an interaction
took place the momentum transfer to both partons involved is delayed by lArlabl / (2c),
where Ahab is the distance vector in the lab frame. In the meantime the partons
are considered not being able to interact, see Figure 414. A technical problem in
the calculations however is that the time delay has to be implemented for integer
numbers of timestep lengths and therefore does not necessarily fully correspond to its

supposed length.

It turns out that the propagation of the outgoing shock-wave becomes more ho-

mogeneous, the initial outburst is suppressed while the distances of the outmost

73

 

Figure 4.14: A ﬁrst step to the introduction of retarded interactions.

scattering events are comparable to the simulations with instantaneous interactions,

see the forth row of Figures 412,413.

4.6.4 Model D: Downscaling the Cross Sections

The maximum signal velocity is proportional to the maximum cross section. Down-
scaling the cross sections to 1 / nth of its original value and compensating this by the
initialization of n times as many testpartons leads to a reduction of the maximum
signal velocity by 1 /\/17. Due to limitations given by the increasing computation time,
for this model calculation a factor of only 71. = 5 has been chosen which corresponds

to a reduction of the maximum signal velocity by a factor of a: 0.45.

It turns out that this technique radically changes the characteristics of the sim-
ulation, the outgoing shock-wave is strongly suppressed, see the ﬁfth row of Fig—
ures 412,413. This outcome is surprising, since Ref. [Wel89] found that for internu—
clear cascade codes the outcome of this “full-ensemble” method (n times as many test
particles) is comparable the usual “parallel-ensemble” method (77. parallel test runs).

In the case of parton cascades apparently this does not hold true anymore due to the

reduction of superluminous signals.

74

Before favoring a method however that leads to such strong changes in the outcome

of the simulation, further analysis is mandatory.

4.6.5 Model E: Wee-Partons

When initializing the parton conﬁguration the Lorentz-contraction leads to very thin
pancakes of nuclear matter, the contracted thickness of a nucleon is around 8 - 10’3
fm, that of the gold nucleus 6 - 10’2 fm. Thereby the z-coordinate of the partons
is ﬁxed rather precisely. However, the momentum of the low—as-components of the
nucleonic wave-function is also determined within the order of a few hundred MeV.
This leads to a violation of the uncertainty principle, which lead to the suggestion
that in any frame of reference the nucleonic wave-function should be smeared out to
a pancake thickness of at least approximately 1 fm, where the low-x-components are
situated further outside and only the valence-quarks are actually within the highly
contracted pancake [Gei92.1, Gei94.1, Gei92.2, Gei93, Gei942, Mc182]. Applying this
kind of initial conﬁguration naturally leads to a smaller parton density; in this case
approximately by a factor of 20. Figure 4.15 shows the parton conﬁguration both at

the beginning of the simulation and at 2 fm/c.

One’s hope can be that the increased width and smaller parton density of the nu-
clei compared to the fully Lorentz-contracted model helps to eliminate some causality
violating effects, see the bottom row of Figures 412,413: The initial outburst of a
shock-wave at impact is eliminated because the nuclei enter each other rather gradu-
ally, an outgoing shock-wave within the respective nuclei is stronger damped because
of a decrease in collision rate, and ﬁnally time-ordering problems on the scale of the
time the nuclei are passing each other are rariﬁed because that time-scale becomes

much longer.

However, there are now other inconsistencies: Although this concept given by the

75

 

 

 

 

 

 

 

 

 

 

- I I I I I . ' I ' I I I I ' -

_ ﬁ -4

5 — —1 5 p— --I
A h : .fi; ::2 i
E ';E I? i 315;}

«... ° _ T o _ :' 3:5 —
v I .. 5.31:2"
' air-in?"-
X ‘it' ‘ - £17113]:

-5 I— —1 -5 — —.

. , ‘ _ ‘ ' .1

l l l I l l l l I l J L LI I J l l I l l l l I l l I l l l l l I I I l l I l 1 #1 L i 1 J 1 I 1 l L l I l l l I
-7.5 -5.0 —2.5 0.0 2.5 5.0 7.5 —7.5 —5.o —2.5 0.0 2.5 5.0 7.5

2 (f m) z ( f In)

Figure 4.15: Initial conﬁguration and conﬁguration after 2 fm/c for a simulation with
wee-partons.

uncertainty relation is truly valid in any frame of reference, within this model there
is no way of implementing it in a Lorentz—invariant way. In reality the components of
the wave-functions transform, in different frames different components of the wave-
function are considered a parton. In the model used only the parton coordinates and
momenta transform, while the partons themselves, once generated, exist in any frame.
While for the initial collisions the c.m.-frames of individual parton collisions nearly
coincides with the lab-frame, in subsequent collisions the c.m.-frames are closer to the
rest-frames of the respective nuclei —- in its rest-frame however, a nucleon, smeared

out to 1 fm in the lab—frame, has a longitudinal radius of 100 fm.

4.6.6 Model F: Proper time approach

As already pointed out, a fully covariant description of a particle collisions would
require a fully 4—dimensional conﬁguration space. This is not possible without giving
up the equal-time character of the simulation. Ref. [Kod84] however proposes a
causality preserving scheme that while retaining the unique global time character of

the simulation minimizes the frame-dependence of the choice of collision partners for

76

the particles, but not the frame-dependent time-ordering of those collisions. Each

particle i is considered to have its own clock showing its proper time 7,,

d7,- = dtI/l — [3.20), 7.0:) = /0tdt’ 1 — 530) . (4.26)

Since the particles do not change their momenta between collisions, the above integral
can be reduced to a sum of products of the type ﬁmAta- With 736(3) being the proper
i-time of the collision of particle 2' with particle j, and 7,0 being the proper i-time of

the most recent collision of particle i, let
674]) :2 Tic(j) — 7'10 (4.27)

be the proper time distance between those two events, a Lorentz-invariant quantity.
The collision instant 736(3) is deﬁned individually within the rest-frame of particle i
through the closest approach to particle j — in the other methods, the closest ap-
proach is deﬁned either within one frame, usually the lab—frame or the c.m.-frame
of the nuclei, or within the respective c.m.-frame of a particle pair. In the previous
discussions, the latter mechanism had been chosen. To co—relate the individual clos-
est approach tests, for collisions only the particle pairs (2', j) are considered for which
both

673(3) =Min {6T,-(l) >0, l=1,...,N;l;éi} (4.28)

and

67,71): Min {6T,-(l) > 0, I = 1,. . .,N;l :ﬁj} , (4.29)

N being the total number of particles in the simulation. Through the restriction
67 > 0 only collisions in the absolute future of each particle are considered. What the
two above conditions mean is that for particle i the very next possible collision (in
the sense of proper time) is with particle j and vice versa. This algorithm only allows

collisions that have no risk of not happening in another frame. Suppose particle 1

 

 

77

has a minimum proper time distance to particle 2, but particle 2 has its minimum

proper time distance to particle 3, then no collision will take place at all.

The search for the very next collision partner of every particle in the simula-
tion is an un-avoidable N 2-problem. Since each of these tests involves a change of
coordinate system, this mechanism is computationally very intense. Therefore, it
was only possible to simulate the initial stages of one single (p,Au)-collision. It was
observed that the number of collisions in this phase dropped dramatically, in this
testrun to about 10% of the number in the other simulations; due to computational
limitations however one is not able to claim signiﬁcant statistics on this percentage.
This outcome is compatible with the conjecture of the authors of Ref. [Kod84] that
the mechanism might underestimate the number of collisions [Kod94]. In smaller
simulations of (p,p)-collisions situations were observed where particle 1 had particle
2 as its closest collision partner, particle 2 had particle 3, and particle 3 again had
particle 1. “Ring”-conﬁgurations like this, possibly spanning over even more parti-
cles, might lead to the underestimation of the collision rate. This effect might have
been enhanced by the high particle density and the extreme Lorentz-contraction of
the nucleons as compared to lower energy internuclear cascades the mechanism was
developed for. Further investigation of this method is required, but due to computa-
tional limitations was not possible within this work. It is however not expected that

the above mechanism can overcome the macroscopic problem of shock-waves.

4.7 Conclusions

The influence of superluminous signal velocities on the signal propagation in parton
cascade codes was found to be smaller than expected. This is mainly due to two

effects: In the initial stages of the interaction the energy-dependent cross sections

 

 

 

 

78

tend to be small. Without this beneﬁcial effect the signal propagation has a threshold
in the region between 0.4 and 0.5 fm2 from where on the velocity of the outgoing
shock-wave reaches the speed of light. In the later stages of the interaction shock-
waves get damped out not because of a lack of interactions but because of signal
propagation that resembles a random walk: not every collision actually leads to an
outward propagation, in fact, the equally probable inward propagation leads to a

virtual damping of the shock-wave.

APPENDIX

Appendix A

The Parton Cascade Program

A. 1 Flow chart

After initialization of the two pancakes by placing the test particles into the con—
ﬁguration space and setting their initial momenta, the program mainly repeats the
steps

0 Moving the particles according to their momenta
0 Checking for collisions between the particles.

0 Decaying of particles.

This is done for a given amount of timesteps.
The program structure is shown in the flow chart Figure A.1.

A.2 Test runs

The program is able to simulate the behaviour of more than one particle system at
a time. Each of the particle systems is called a ’test run’. Since scattering is done
by Monte Carlo methods each of the conﬁgurations which are initially similar will
behave differently on the long run. The test runs however are set up to start with
slightly different initial conditions, too.

The test runs can be used to ﬁnd averaged properties of the systems. They
can also be used to calculate mean ﬁelds for the particles, but this is not currently
implemented. As long as no mean ﬁeld is taken into account, the test runs are
independent from each other.

79

 

 

80

  

  
     

No formal
reasons against
scattering.

 

 

Read parameters

 

|

 

Initialization

 

I

 

Outputs

 

I

 

 

Check constants

 

 

 

 

 

 

Permutations

 

 

 

Move particles
according to
momenta and
put them into
boxes

 

 

 
 
   

 

 

Detemiine ﬁrst box I

 

l

Determine second box

 

 

Already checked?

Check all
narticle comb.

All neighbor boxes?
no es

no All boxes checked?

yes

 

 

Time for outputs?

no yes

  

 

I Outputs

 

 

/

structure

 

 

I
/ [ Check particle data
I

 

 

I Check constants

 

no

   

Multistage

decision

They will collide

 

Collide particles

 

I

 

Set ﬂags

 

 

 

_____ _____._______-..__

 
  

     
   
 

Creating
and deleting
particles

BUU for high energy
collisions
01/31/94

 

More timesteps?

> Ready

Figure A.1: BUU program ﬂow chart

81

A.3 Scattering

As already explained running the simulation consists mainly of two steps: Moving
the particles according to their momenta as one step, scattering the particles as
another step. An additional step can be the decaying of partons. In principle the
scattering step involves the examination whether or not scattering takes place for
every particle pair in the test runs. This would result in an N 2-time dependence of
this step. However since the spatial range of the interactions taken into consideration
is limited, most of the time particle pairs that will never have a chance of being
scattered would be examined.

To overcome this waste of time the coordinate space is divided into boxes con-
structed in a way such that only particle combinations with members of the same
or neighboring boxes have to be taken into account. The boxes establish a grid in
conﬁguration space.

This implies that initialy in the ac and y direction the edges of the boxes have

0 . . . . .
to be longer than 2‘/-BEE With amax being the max1mum cross section taken into
7r

account. In the z direction the edges have to be longer than 2tstep - c with tstep being
the timestep length. During program execution the 2 extension has to be adjusted as
particles are scattered from their original z or —z direction respectively. With gmax
being the maximum angle any one particle has from the initial direction, the grid

 

. 0 .
extension in z direction is adjusted to 2tstep - c + 2‘/ max Sln(9max).
II

The total extension of the grid has to be chosen such that during the whole
simulation no particle leaves the grid, the latter has to be checked during the particles
are moved. The assignment of the particles to the boxes is also done during the move-
step.

In the scattering step all box combinations between neighboring boxes have to be
selected, and for each box combination all particle combinations have to be examined.
This is done in the following way:

To prevent the occurence of preferential directions always ﬁrst one box is chosen
at random, then in a random order one of the neighboring boxes or itself is chosen as
a second box. It is made sure that each box combination is only taken into account
once. Then all particle combinations with one particle from the ﬁrst box and the
other one from the second box are checked.

Before the physics of a possible scattering is checked, a few formal knock out
criteria are applied:

0 Neither of the particles has already scattered in this timestep. This is done
to avoid possible violations of causality: One timestep length is considered the
smallest scale on which any one thing can happen to a particle.

0 Not both of the particles have had the respective other one as their last scatter—
ing partner. This is done to avoid multiple scattering, a phenomenon already
taken into account on the level of elementary processes.

 

 

82

0 Not both of the particles haven’t scattered with a particle from the respective
other nucleus of origin before. This is done to prevent scattering between parti-
cles from the same nucleus before any collision with the other nucleus had taken
place — since no mean ﬁeld guarantees the stability of an unperturbed nucleus,
stability is provided in this artiﬁcial way by the lack of any interaction. In prin-
ciple the same effect can be achieved by having a cutoff energy that suppresses
interactions with a too small c.m.-energy.

A.4 Assignment of the particles to the boxes

Since during the scattering process always ﬁrst the boxes are chosen, it is necessary
that, given a box number, all particles in that box can be accessed rapidly. However,
since the particle data is stored independently of any box membership, assignment
to the boxes is rather done by a pointer structure being established each time new
when the particles are moved.

This movement of the particles is done in a random order; each box possesses one
pointer that points to the ﬁrst particle being found in the box area after moving, or
to nowhere if the box area turned out to be empty.

Each particle itself possesses a pointer that points to the next particle being found
in the same box, or to nowhere if it was the last particle being found there.

As a result all particles in a certain box can be accessed by following this chain of
pointers, starting with the boxes’ own pointer and continuing through the particles’
’next-particle’-pointer, until any of those points to nowhere.

A.5 Checking for collisions

After two particles have been admitted for scattering by formal reasons (see A.3 on
page 81) the physics criteria for scattering or not are checked, in particular:

0 Are the two particles close enough for the particular c.m.-energy dependent
total cross section of the reaction?

0 Have they reached their point of closest approach in the c.m.-system?

While the ﬁrst condition can be checked very quickly, performing the other two tests
takes much more computation time since it involves a transformation into the c.m.-
system of the two particles.

To avoid this as long as possible, a multistage decision is applied: Weaker cri-
teria which take less computation time are checked ﬁrst, the decision sequence is
immediately aborted as soon as one of its requirements is not met.

 

83

A.6 Decay of particles

Off-shell particles that don’t belong to a certain nucleon anymore and haven’t collided
with another particle in the actual timestep are assumed to be able to decay. Due to

the uncertainty relation
AE - At ~ 1

their lifetime tlife is to be chosen according to
Am ' tlife N 1. .

Am is the off-shellness of the particle and calculated according to

Am = II/PuP" — mon-shelII -

momsheu is the on-shell mass of the particle, 19,, its four momentum.

In order to ﬁnd the probability p for them to decay in the actual timestep the
livetime tlive is assumed to be the expectation value of possible lifetimes:

00
Z: pn ntstep

- n=0

tlive = 00
Z 1)..
n=0

Here pn is the probability for the particle to survive It timesteps with a total lifetime
of ntstep. This must be equal to the probabilty of being not “killed” in n timesteps:

 

pn = (1 _ 19)"
Therefore
in _t 2(1 -p)”n
1e "‘ t *—
8 ep 2(1 _ p)"
and with
1 (9 a:
n = — => " = —— < 1
:1: l—zr’ 8:1: 237” (1—a:)2 I33I
ﬁnally
1 _
tlive = tstep—p
and
_ Amtstep
p 1+Amt,,e,, '

In the decay step for all particles in question, a random number between 0 and 1 is
determined and compared to p. If it is smaller than p the particle decays into two
other particles.

 

84

A.7 Creating and deleting particles

Since during execution oftentimes particles are created or deleted, an efficient method
had to be found for this.

As all particle information is stored in static ﬁxed-size arrays, the simplest way
would be to always have those arrays ﬁlled up to a variable index. This maximum
index would be stored elsewhere. Whenever a new particle is created, it would be
placed on top of the others and the maximum index would be increased, whenever a
particle is deleted, the stack above it would have to fall down and ﬁll the gap, and
the maximum index would have to be decreased by one. However the “falling down”
involves a lot of copying, also, a lot of pointer information would have to be updated.

The other extreme would be not to care about maximum indices and instead
only to provide a flag for each array entry whether or not it is used. When a new
particle is to be created one would use some free entry in the array, write in the
information and ﬂag it ’used’. Deleting a particle would simply be done by marking
the corresponding entry ’unused’. No pointer structures would be updated, one would
always run through the whole array or pointer chains and just skip ’deleted’ particles.
However, if the arrays are only ﬁlled rarely, this would be a waste of time.

So a combination of both methods is chosen:

0 The array entries can be marked ’unused’.

o The maximum used index number is stored elsewhere, all loops only have to
run up to this index.

0 The pointer chain for the box contents is updated.

When a particle is created, beginning from the bottom of the array the ﬁrst unused
position is searched for and the information written to it. The entry will be marked
’used’. If it is above the maximum index so far (which can only be by ’1’) the
maximum index is updated. The new particle is not incorporated into any pointer
structures since as it already had interacted it would not take part in any other
interaction in this timestep anyway.

When a particle is deleted its array position is marked ’unused’. The particle is
removed from the box pointer chains, the maximum array index is updated.

Appendix B

Program listing

This is not a complete listing of the program, but only a selection of the most relevant
subroutines and ﬁles.

B.1 STDEC

Gerd Kortemeyer
Last revision of ﬁle: 08/31/94

This is the include ﬁle with all declarations of the common blocks used for communication between
the subroutines.

It also sets

c Type check on
implicit none
c

which forces the programmer to declare any local variables explicitly. Some physics parameters are
also set.

The standard declarations are to be included in any subroutine that is supposed to manipulate
global data. This is done simply by adding

include ’stdec.for’

Parts of this ﬁle should not be copied directly into subroutines since this might result in a lot of
confusion if at some stage the deﬁnition of common-blocks is changed in the standard declarations
ﬁle. The compiler will not be able to ﬁgure out that the deﬁnition of some common block is different
in different subroutines and map the variables wrong.

B.1.1 Physics parameters

The include ﬁle deﬁnes the following constants:

85

 

 

86

c Mass of nucleon [GeV/c**2]
double precision manuc
parameter (manuc=0.939)

c hbar*c [GeV*fm]
double precision hbc
parameter (hbc=0.197327)

c r0 [fnﬂ
double precision r0
parameter (r0=1.2)

B.1.2 Other parameters

Parameters for the dimensions of arrays:

0000

Maximum number of quarks
integer mquark
parameter (mquark=15000)

c Maximum number of test runs

integer mtest
parameter (mtest=5)

c Grid-dimensions
integer mgridxy,mgridz
parameter (mgridxy=20)
parameter (mgridz =250)

c Dimensions for cross section array
integer maxpoints
parameter (maxpoints=3000)
double precision sqstep
parameter (sqstep=1000./(maxpoints-1.)**4)

c Maximum particle type index
integer maxtype
parameter (maxtype=6)

c Maximum number of channels
integer maxchannel
parameter (maxchannel=6)

c Very small number

double precision eps

parameter (eps=.0000000001)

c The number PI

double precision pi
parameter (pi=3.1416)

B.1.3 Particle data

In the common block particles all information that goes with each of the single particles is collected.

The particles are identiﬁed by two numbers: The number of the test run they are in, ranging from
1 to mtest, and the particle number in that run, ranging from 1 to mquark.

 

 

87

As the arrays are dimensioned staticly the two variables ntr and nq(mtest) are included. ntr (int)
holds the number of test runs really active, it is set in subroutine parmin from the information in
the parameter ﬁle.

nq(mtest) (int) gives the maximum particle-index per test run. It is ﬁrst set in subroutine init,
but can vary when new particles are produced or other particles are absorbed. This is done in the
procedures crepar (see B.12.2, p. 152) and delpar (see B.12.1, p. 149).

So if one wants to construct a loop over all particles a typical structure would be
do 10 test=1,ntr
do 20 particle=1,nq(test)
This doesn’t necessarily mean that all array entries between 1 and nq(test) are actually used.
The particle information includes
a tq(mtest ,mquark) (int), the parton type. This is 0 for gluons, 1 for up-quarks, 2 for down—
quarks, 3 for strange-quarks, 4 for charm-quarks, 5 for bottom-quarks and 6 for top—quarks.
The accompaining anti—quarks have got the same but negative number, so for example -1

stands for anti—up. This parton code corresponds to the codes deﬁned in the CTEQ—routines
for parton-distributions.

rq (mtest ,mquark , 3) (double), the position of the particles. The last array index corresponds
to the z, y and 2 components of the position vector respectively. The unit is fm. The time
component of the position is deﬁned by the time step.

pq(mtest,mquark,3) (double), the momentum of the particles. Again the last array index
corresponds to the spatial components of this vector. The unit of the spatial components is
GeV/c. The energy component is given separately by the array eq(mtest ,mquark) (double),
its unit is GeV.

nnq (mtest ,mquark) (int), the number of the nucleon the particle was assigned to during the
initialization. The total number of nucleons is given by nntar and nptar for the number of
neutrons and protons in the target, and by nnpro and nppro for the number of neutrons and
protons in the projectile. The numbers nntar etc. are set in parmin by the user, nnq is set
in init.

A particle that was scattered by subroutine collide will have nnq=0 afterwards, a particle
created by subroutine crepar will have nnqz-l.

fq(mtest,mquark) (int) is a flag wether the particle belongs to the projectile (fqz—l) or
to the target (fqzl), they are ﬁrst set in subroutine init. This ﬂag is used in subroutine
scatter to prevent particles from the same nucleus to scatter with each other before at least
one of them has scattered with a particle from the other nucleus. After their first scattering
with any particle their ﬂag is set to fq=0. A particle with fq=0 can then scatter with any
other particle.

last (mtest ,mquark) (int) is a pointer to the particle in the same test run the actual particle
has last scattered with, so values for last naturally range from 1 to mtest.
But last has a second usage: last=-1 means that the particle has been deleted in the last
interaction. So last also flags empty array positions.
However particles that haven’t scattered before carry last (i , j)=j as preset by the subrou-
tines init or crepar.
So to write a loop over all existing particles one should code

do 10 test=1,ntr

do 20 particle=1,nq(test)
if (last(test,particle).ne.-1) then

 

88

The values for last are set in subroutine collide and used in subroutine scatter to prevent
particles from being engaged in multiple scattering with each other. Two subsequent scatter-
ing events between the same particles are forbidden, however it is possible to scatter particles
1 and 2, then 2 scatter with 3, ﬁnally 1 again with 2. This results in a kind of 3-body-force.

nextgrid (mtest ,mquark) (int) points to another particle of the same test run that is in the
same space box as itself. If it is the last particle in that chain the value for nextgrid is -1.
The procedure is decribed in the section for procedure move, B.5 on page 114.

lastpoint (mtest ,mquark , 3) (double) The most recent point in the lab frame that the parton
scattered at. The scattering point is assumed in the middle between both partons involved.

1asttime(mtest ,mquark) (int) The most recent timestep-number that the parton scattered
in. Both lastpoint and lasttime are needed to calculate the signal velocity from one
scattering event to the next.

fdue (mtest,mquark) (int), feq(mtest,mquark) (double), qu(mtest,mquark,3) (double)
This is data for experiments with retarded interactions, f due is the amount of timesteps that
the parton is considered to be involved in the interaction, feq and qu are the energy and
the momentum that the parton is supposed to get.

B.1.4 Other parameters

In the common block others control parameters for the program run are collected.

In particular there are

nt ime (int) The actual number of the timestep being executed at the moment.

timstp (double) The length of one timestep in fm/c. This value is read from the parameters
ﬁle.

proxO (double) The a: coordinate in fm of the projectiles starting position, deﬁned in the
parameters ﬁle.

pron (double) The 2 coordinate in fm of the projectiles starting position, deﬁned in the
parameters ﬁle.

iseed (double) The starting number for the random number generator, deﬁned in the pa-
rameters ﬁle.

nmts (int) The total number of timesteps to be performed, deﬁned in the parameters ﬁle.

whenout (int) The number of timesteps to be performed between two ’big’ outputs, for ex-
ample plots of the particle position in coordinate Space, output of statistics. . .. Read from
the parameters ﬁle.

elabpnpro (double) The kinetic energy in GeV per nucleon of the projectile, deﬁned in the
parameters ﬁle.

elabpntar (double) The kinetic energy in GeV per nucleon of the target, deﬁned in the
parameters ﬁle.

maxsig (double) The maximum total cross section in fm2 as being found during the initial-
ization of the cross section tables xsec. . . , see subroutine initcross (B.3.7, p. 106).

gridspxy (double) The a; and y length in fm of the box edges. The value is determined
according to the maximum cross section maxsig.

89

o gridspz (double) The 2 length in fm of the box edges. The value is re-determined after every
timestep according to the maximum angle a parton is headed away from it’s forward direction
in connection with maxsig and the timestep length timstp. This is done in subroutine move
(see B.5 on page 114):

gridspz=2 . * (t imstp-I-sqrt (maxs ig/ pi) *sin (mtheta) )
o gridz (double) is the bin-width for conﬁguration plots.

The items elabtar and elabpro became obsolete and were removed with the introduction of a new
initialization method for the particles. The item lowcut became obsolete with a new treatment of
low energy collisions.

B.1.5 Grid information

In the common gridinfo information for the spatial grid is stored.

In the integer array
gridroot (mtest,-mgridxy:mgridxy,-mgridxy:mgridxy,-mgridz:mgridz)

the pointers to the ﬁrst particles found in the boxes are stored. A ’-1’ stands for ’no particle’,
meaning ’the box is empty’.

The ﬁrst index represents the test run number, the following indizes the box coordinate. In the
a: and y direction the boxes are numbered from -mgridxy to mgridxy with box number -mgridxy
being the leftmost or lowest boxes respectively. The 2 coordinate is described by the last index.

The integer arrays
permx ( -mgridxy :mgridxy) , permy (-mgridxy :mgridxy) , permz (-mgr idz : mgr idz)

are permutations of the numbers from -mgridxy to mgridxy or from -mgridz to mgr idz respectively
and are used to have access to the boxes in a random way.

Instead of accessing gridroot (test , i , j ,k) during scattering rather
gridroot(test,permx(i),permy(j),permz(k))

is accessed. The random order of accessing the boxes is necessary because with each box also the
neighboring boxes will be accessed and a ﬁxed order for that would result in preferential directions.

The integer arrays
deltx(3),delty(3),deltz(3)

have the same function for accessing the neighboring boxes in a random way, they have a permutation
of the numbers from —1 to 1.

90

B.1.6 Initialization parameters

In the common block initparms information concerning the different methods of initialization are
stored as they are being read from the parameter ﬁle.

B.1.

initmeth (int) In this variable the number of the initialization method to be used is stored.
It is read in parmin and used in makenucleus (see B.3.2 on page 103).

initxmin (double) The minimum a: value to be used in the parton distribution function
f(:r, Q2), see B.14 on page 157.

initqsq (double) The Q2 value to be used in the parton distribution function f(:1:, Q2), see
B.14 on page 157.

initrnuc (double) The nucleon (’bag’) radius assumed.

7 Information block for the initialization methods

In the common block create all necessary information for the different initialization methods is
stored It has to be ﬁlled before makenucleus is called. In particular:

B.1.

crelabpn (double) The kinetic energy per nucleon of that nucleus in the lab frame.
crexO (double) The initial 3: position of the nucleus’ center in the lab frame.

crezo (double) The initial 2 position of the nucleus’ center in the lab frame. The y position
is chosen to be zero always.

crenp (int) The number of protons in the nucleus.
crenn (int) The number of neutrons in the nucleus.

cremodir (double) The direction for the momentum: ’1’ means ’travelling in positive 2 di-
rection’, ’-1’ means ’travelling in negative 2 direction’.

creflag (double) The flag the particles from this nucleus should have. This corresponds to
the array fq in the particle data block.

8 Data for two particles that are scattering candidates

In the common block twopart gradually information about two particles that are deﬁned by sub-
routine scatter (B.6, p. 118) for being candidates for scattering is collected.

The data is mainly set and used in subroutine collide (B.8, p. 123), but also in twototwo (B.9,
p. 139).

In particular the entries are

tppl (int) The number of the ﬁrst particle being chosen.

tpp2 (int) The number of the second particle being chosen.

tprn (int) The number of the test run the particles are from.

tpsq (double) The total energy squared in GeV2 0f the two particles in their c.m.-system.
tpr1(3) (double) The distance vector of the two particles in the lab frame in fm.

tpeql , tpeq2 (double) The energies of the two particles in GeV in the c.m.—system.

 

91

tpmol (3) (double) The vector of momentum in the c.m.-system for the ﬁrst particle in fm/e
before scattering.

tpbeta (double) The transformation vector 3 to enter the c.m.-system.

tpeql (double) The total energy of the two particles in the lab frame.
tpbq (double) The value of [32.

_ 1

o tpgamma (double) The value 7 — M.
o tpsigma (double) The total cross section of the two particles in fm2 as being set by subroutine
totalcross (B.8.3, p. totalcross). However, for technical reasons after xsecentry (B.8.4,
p. 136) has been called, the variable includes the partial cross section for the speciﬁed process.

tpgambet (double) The value of 3331.

tptd (double) The time difference between the two particles when seen in the c.m.-frame.

tprc(3) (double) The spatial distance vector in the c.m.-frame.

tpcos (double) The cosine of the scattering angle 0 in the c.m.-frame.

tpsc1(3) (double) The momentum vector of the ﬁrst particle in the c.m.-system after scat-
tering. It is produced from tpmo1 by subroutine newmomenta.

tpsceql, tpsceq2 (double) The energies of the two particles after scattering.

tpchin (int) The channel number of the two incoming particles. The channels are coded in
the following way:

Channel number | particle content

99

‘19

qq

qt?

qq’

q’ tj’ — outgoing channel only

 

0301;:me

The prime denotes a quark different from the unprimed quark, the bar denotes the respective
antiparticle.

tpchout (int) The channel number of the two outgoing particles.

tpcrindex (int) The energy index for the arrays xsec. . ., see B.1.16 on page 93.

tpdq (double) The squared distance in fm2 between the two particles in the lab frame.

tpdqc (double) The squared distance in fm2 between the two particles in the c.m.-frame.

tinl , tin2 (int) The parton types of the two incoming partons as being set in totalcross
(B.8.3, p. 135).

toutl , tout2 (int) The parton types of the two outgoing partons as being set in outchannel
(B.8.2, p. 132).

The item tpeql became obsolete.

 

92

B. 1.9 Already scattered

In the common block already there is only one array alr (mtest ,mquark) (int) that ﬂags for each
particle if it had already been scattered in this timestep or not. ’1’ denotes ’scattered’, ’0’ denotes
’not scattered’.

It is reset in subroutine scatter, set in twototwo and used in both scatter and decay.

B.1.

10 Statistics

The common block stat includes some statistics information that is being printed out to the protocol
ﬁle during ’output’.

minsigma, maxsigma (double) The minimum and maximum values in fm2 for the cross
sections occurred between the outputs.

minsq, maxsq (double) The minimum and maximum values for the c.m.-energy squared in
GeV2 occurred between the outputs.

mincos, maxcos (double) The minimum and maximum cosine of the scattering angle in the
c.m.—frame as determined by the montecarlos.

nstaO (mtest) , . . . , nsta4 (mtest) (int) Counters for the multistage check: They are incre-
mented each time one particle combination passes the ﬁrst, second,. . . stage of the test. This
can be used to supervise the efﬁciency of the test. For each testrun there exists a separate
counter.

timefault (int) This counter is incremented each time the time-distance between two par-
ticles that occurres during transformation to the c.m.-system is bigger than the transformed
timestep length.

anglestat (maxchannel,6) (int) With this array the c.m. angular distribution of scattering
events in the respective output channels is examined. The second index belongs to bins with
cos 6 from -1 to —%, —% to —%, . . . The contents of this array are written to the protocol ﬁle
at every ’big’ output.

labang1e(6) (int) With this array the lab angular distribution of scattering events is exam-
ined. The index belongs to bins with c030 from —1 to —%, ~37 to —%, ...The contents of
this array are written to the protocol ﬁle at every ’big’ output.

chscat (maxchannel,maxchanne1 (int) In this array the scattering events from the channel

number in the ﬁrst index to the channel number in the second index are counted.

It also includes information about the number of particles created or deleted:

B.1.

crenumber, delnumber (int) The number of particles created / deleted since the last output.

lastnumber (int) The number of particles at the previous output, used for checks.

11 Information for particles to be created

The common block newpart has to be used when a new particle is to be created by subroutine
crepar, see B.12.2 on page 152.

In particular the following information has to be provided: The position nequ(3), the momentum
newpq(3), the energy neweq, the parton type newtq, and the testrun newtest that the particle
should be created in.

The item newpscat became obsolete with the removal of the array pointscat.

93

B.1.12 Permutation

The common block perm only has one entry, the array perm(mquark) (int). It is initialized to hold
all numbers from 1 to mquark, but in a permuted way. This array can be used to access the particles
in a random way by ﬁrst building a subpermutation with the appropriate number of particles and
then accessing the particles always through the subpermutation, see B.5.3 on page 116.

B.1.13 Timing information

The common block time has the variable to (real) that can be used to ﬁnd the total CPU-time so
far by calling secnds (t0).

B.1.14 Controls

The common block control has integer entries determining wether or not certain outputs are writ-
ten. ’0’ means ’no’, ’1’ means ’yes’. The values are read from the parameters ﬁle. In particular these
are

0 ploto (int) Wether or not the conﬁguration and rapidity plots are done.

a scattering (int) Wether or not information about scattering events is written out, see sub-
routine twototwo on page 139.

Both the items causalpl and velocitypl became obsolete with the introduction of scattering.

B.1.15 Charges

In the common block charges there is one array charge that has the charge in units of the electron
charge for every parton type.

B.1.16 Cross sections

In these arrays the partial cross sections for different processes are stored. The entries are of type
double precision:

xsecgg22(maxpoints),
xsecqq22(maxtype,maxpoints),
xsecqg22(maxtype,maxpoints),
xsecqu22(maxtype,maxpoints),
xsecqqp22(maxtype,maxtype,maxpoints),
xsecqugg(maxtype,maxpoints),
xsecggqu(maxtype,maxpoints),
xsecquqqup(maxtype,maxtype,maxpoints)

The last index always refers to the energy mesh the cross sections are tabulated on. The preceding
indices refer to the parton types involved. The channel combination is coded in the name of the
arrays. For a description of these cross section arrays see B.3.7 on page 106.

For the incoming channel 99 the open output channels are 99 or qq. The partial cross sections
for two incoming gluons are added up in the order gg —> 99 (xsecgg22), gg —> mi, 99 —> dd, ...,

gg —> tt- (xsecggqu).

94

For the incoming channels qg and qq only the elastic channel is open.

For the incoming channel qrj the partial cross sections are added up in the following order: qrj —) qr]
(xsecqu22), qtj —> gg (xsecqugg), qrj —) q’q’ (xsecquqqup). The cross section array entries for
(16 —> q’ (7’ are added up in the order ch —) m], qrj —> dd and set to zero when q = q’.

As one can see, the partial cross sections always add up to the total cross section in the last entry
of that sequence, the difference between the entries and therefore the relative probabilty for the
channel are the partial cross sections themselves. For the initialization of the xsec. . . arrays see
B.3.7 on page 106.

B.1.17 Parameters for the cross sections

The following parameters are read from the parameters ﬁle and used in the cross section and mon-
tecarlo routines:

0 cm (double) is the cutoff for the gluon propagator in the t and u channels.

0 cm (double) is the cutoff for the gluon propagator in the 3 channels

0 crk (double) is the cutoff for the quark propagators in any channels

0 cralphas is the value for as.

B.1.18 Causality information

The common block causal with the items steps, stymax, stcount, veymax, vc(0 :mvelocity) and
vgrid became obsolete.

B.2 The main ﬁle

Gerd Kortemeyer
Last revision of the ﬁle: 04/27/94

This ﬁle consists the main program and includes the other major parts.

B.2. 1 Local variables

c Date/Time
character dat*9,tim*9

c Output-Counter, do variables
integer outp, i, j

c For initialization of random number generator
double precision dummy

B.2.2 Some initializations

First, date and time are read out and the timer for the CPU-seconds is reset

 

95

c Get date/time
call date(dat)
call time(tim)

c Reset timer
t0=secnds(0.0)

Next, the protocol ﬁle is opened:

c
c File for all control outputs
c

open(1,file=’buu_protocol.dat’,status=’unknown’)
c

write(1,*) ’ *** BUU for high energy collisions ***’
write(1,*) ’ ’

write(1,*) ’National Superconducting Cyclotron Laboratory’
write(1,*) ’G.K. 1993, 94’

write(1,*) ’Date: ’,dat,’ Time: ’,tim

write(1,*)

& i )

 

The parameters are read in from the parameter ﬁle:

c
c Initialization of commons
c

call parmin()

Now, the random number generator is initialized with the iseed given in the parameter ﬁle. The
random number generator is designed to start a new sequence of numbers when a negative value is
given as the argument:

c
c Reset random numbers
c
iseed=~abs(iseed)
dummy=ran1(iseed)

The particle data and some other data is initialized by subroutine init.

c
c Initialization of particle data
c

call init()

The arrays for the causality and velocity plots are reset:

steps=int(nmts/50.)
stcount=0
stymax=0

c Reset causal information

100

96

do 100 i=0,mgridxy
nc(i)=0
continue

vgrid=mvelocity*timstp/sqrt(maxsig/pi)

veymax=0

c Reset velocity information

110

do 110 i=0,mvelocity
vc(i)=0
continue

The graphics output ﬁles belonging to the output options given in the parameters ﬁle are opened:

C

c File for graphics output

C

if (ploto.ne.0) then
open(3,fi1e=’buu_plots.plt’,status=’unknown’)
write(3,*) ’SET STORAGE X Y Z SIZE=1000OOO’
endif

if (causalpl.ne.0) then
open(4,file=’buu_causal.plt’,status=’unknown’)
write(4,*) ’NEW FRAME’
write(4,*)
& ’TITLE 1 9.7 SIZE 2 ”HIGH ENERGY BUU ’,dat,’ ’,tim,’CAUSALITY”’
write(4,*) ’CASE ” LLL LLLLL ’,
&’ LLLLLLLL”’
write(4,*) ’SET FONT DUPLEX’
write(4,*) ’SET LABELS SIZE=2’
write(4,*) ’READ MESH’
write(4,*) ’FOR Y=’,(i*gridspxy,i=0,mgridxy)
endif

if (velocitypl.ne.0) then
open(5,file=’buu_velocity.plt’,status=’unknown’)
write(5,*) ’NEW FRAME’
write(5,*)
& ’TITLE 1 9.7 SIZE 2 ”HIGH ENERGY BUU ’,dat,’ ’,tim,’VELOCITY”’
write(5,*) ’CASE ” LLL LLLLL ’,
&’ LLLLLLL”’
write(5,*) ’SET FONT DUPLEX’
write(5,*) ’SET LABELS SIZE=2’
write(5,*) ’READ MESH’
write(5,*) ’FOR Y=’,(i/vgrid,i=0,mvelocity)
endif

The procedure plotout is called to plot a picture of the initial particle distribution in conﬁguration
space and the initial rapidity distribution:

if (ploto.ne.0) call plotout()

 

 

97

The procedure constant calculates some constants of motion and writes them to the protocol ﬁle:

call constant()

Now the big loop over all timesteps starts. The counter outp counts the timesteps since the last
’big’ output, the parameter whenout determines after how many steps output is done.

c
c Loop for all timesteps
c

outp=1
c

write(*,*) ’Starting: ’,dat,’, ’,tim
c

do 10 ntime=1,nmts
c

This is the kernal of the loop: Calling move, scatter and decay. However, since 1 —) 2 parton
processes are not yet implemented, the corresponding call statement is commented out.

call permute()

call move()

call scatter()
* call decay()

Now it is checked if it is time for ’big’ output:

c Graphics and other output?
if (outp.eq.whenout) then

 

 

c
call date(dat)
call time(tim)
c
write(*,*) ’Timestep’,ntime,’ of’,nmts,’: ’,dat,’, ’,tim
c
write(1,*)
&’ l
write(1,*)
&1 )
write(1,*) ’Timestep:’,ntime
write(1,*) ’Date: ’,dat,’ Time: ’,tim,’ CPU:’,secnds(t0)
write(1,*) ’-’
if (ploto.ne.0) call plotout()
write(1,*) ’-’
write(1,*) ’Multistage collision check statistics:’
do 20 i=1,ntr
write(1,*) ’Multistage: ’,
& nsta0(i),nsta1(i),nsta2(i),nsta3(i),nsta4(i)
20 continue

write(1,*) ’-’
write(1,*) ’Channel combinations:’
do 30 i=1,maxchannel

 

98

write(1,*) (chscat(i,j),j=1,maxchannel)
30 continue
write(1,*) ’Angular distributions, cos(theta) in c.m.framez’
do 40 i=1,maxchannel
write(1,*) ’ --------------------------------------- ’
write(1,*) ’Outchannel’,i,’:’
write(1,*) ’-’
write(1,*) (anglestat(i,j),j=1,6)
40 continue
write(1,*) ’Angular distributions, cos(theta) in lab frame:’
write(1,*) (labangle(j),j=1,6)
write(1,*) ’Timefaults ’,timefault
write(1,*) ’Min/Max cm sq.: ’,minsq,maxsq
write(1,*) ’Min/Max Sigma : ’,minsigma,maxsigma
write(1,*) ’Min/Max Cosine: ’,mincos,maxcos
write(1,*) ’Gridspacing z : ’,gridspz

The procedure checknumbers prints out statistics about particle numbers and checks them. It is
also able to track particles in the grid chains, but this feature is not used here. The procedure
restat resets the statistics and is part of the init major ﬁle.

call checknumbers(0,0,0)
call restat()
call constant()

outp=0
endif
outp=outp+1
c
10 continue

B.2.3 The end

Now a few ﬁnal plot commands are produced:

if (causalpl.ne.0) then

write(4,*) ’SET LIMITS Y FROM 0 TO’,(stymax+1)*gridspxy

write(4,*) ’CONTOUR LABEL=OFF’

write(4,*) ’CONTOUR 1 LABEL=OFF DOTDASH’

write(4,*) ’PLOT AXES’

write(4,*) ’TITLE X CENTER ”TIME IN FM/C”’

write(4,*) ’TITLE Y CENTER ”SQRT(X*X+Y*Y) IN FM”’
endif

if (velocitypl.ne.0) then

write(5,*) ’SET LIMITS Y FROM 0 TO’,(veymax+1)/vgrid
write(5,*) ’CONTOUR LABEL=0FF’

write(5,*) ’CONTOUR 1 LABEL=0FF DOTDASH’

write(5,*) ’PLOT AXES’

write(5,*) ’TITLE X CENTER ”TIME IN FM/C”’
write(5,*) ’TITLE Y CENTER ”VELOCITY IN UNITS OF C”’

endif

99

c
call date(dat)
call time(tim)
write(1,*) secnds(t0),’: *** Program successfully terminated ***’
write(1,*) ’Date: ’,dat,’ Time: ’,tim

c

stop ’That”s all folks!’
end

B.2.4 The include section

The major parts of the program are included here. The following example is taken from a test
version of the program:

 

Includes

 

0***

Procedure to read parameters
include ’parmin.for’
Procedure to init particle data
include ’init.for’
include ’testinit.for’
Procedure to plot data
include ’plotout.for’
Procedure to move all particles according to their momenta
include ’move.for’

0

0*

0

c Procedure to check for scattering
c with causality violation ...
* include ’scatter.for’
c and without
include ’causcatter.for’
c Procedure for particle decay

include ’decay.for’
Function to calculate on-shell masses
include ’mshell.for’
Procedure that decides wether or not two particles scatter - and scatters
include ’collide.for’
include ’testcollide.for’
Procedure that completes two to two scatterings
include ’twototwo.for’
* include ’testtwototwo.for’
Procedure to randomly exchange entries in permutation array
include ’permute.for’
Procedure to set new momenta in c.m. frame according to scattering angle
include ’newmomenta.for’
include ’testnewmom.for’
Procedure to calculate random numbers
include ’random.for’
Parton distribution function
include ’distri.for’
Procedure to check constants
include ’constant.for’
c Procedure to remove and create particles

0

0

0

0

0

0

0

0

 

100

include ’credel.for’
This is where new initialization methods have to be linked:

c Initialization procedures
include ’initmeth1.for’
include ’initmeth2.for’
include ’initmeth3.for’
c Gauss Legendre Weights and Abscissae
include ’gauleg.for’
Procedures that calculate total cross sections depending on c.m. energy
Test versions
include ’testsigqq.for’
include ’testsigqqp.for’
Procedures doing the Monte Carlo, Testversions
include ’testmonqq.for’
include ’testmonqqp.for’
Procedures for total cross sections and Monte Carlos by Joelle
include ’MONTEGG.FOR’
include ’MONTEGGQQB.FOR’
include ’testmonteggqu.for’
include ’MONTEQG.FOR’
* include ’testmontqu.for’
include ’fmontqu.for’
include ’MONTEQQ.FOR’
include ’MONTEQQB.FOR’
include ’MONTEQQBGG.FOR’
include ’MONTEQQBP.FOR’
include ’MONTEQQP.FOR’
include ’TSIGGG22.FOR’
* include ’tsiggg2205.for’
include ’TSIGGGQQB.FOR’
include ’TSIGQG22.FOR’
include ’TSIGQQ22.FOR’
include ’TSIGQQBQ2.FOR’
include ’TSIGQQBGG.FOR’
include ’TSIGQQBQPQBP.FOR’
include ’TSIGQQP22.FOR’

0**0**00

*1-

B.3 INIT

(lerd I(ortenieyer
Last revision of ﬁle: 07/21/94

In this ﬁle there are most of the routines for the initialization of the particle data except for the
routine that creates the particles. This exception was made to make it possible to have more than
one way of initialization.

Also the permutation arrays and other control data are initialized here.

From outside this ﬁle only the subroutine init is called.

 

 

101

B.3.1 Subroutine init

In the subroutine init ﬁrst the number of particles for all the testruns nq(test) is set to zero:

c
c No particles so far
c
do 10 test=1,ntr
nq(test)=0
10 continue

Then the particles for the target are initialized by the subroutine makenucleus (see B.3.2 on page
103). All the data necessary for the initialization of this nucleus is ﬁrst to be stored in the common
block create.

This is
o The kinetic energy per nucleon of that nucleus in the lab frame.
0 The initial :1: position of the nucleus’ center in the lab frame.

0 The initial 2 position of the nucleus’ center in the lab frame. The y position is chosen to be
zero always.

0 The number of protons in the nucleus.
0 The number of neutrons in the nucleus.

0 The direction for the momentum: ’1’ means ’travelling in positive 2 direction’, ’-1’ means
’travelling in negative 2 direction.

0 The ﬂag the particles from this nucleus should have. This corresponds to the array fq in the
particle data block.

c

c Initialize target

c
write(1,*) ’ —————————————————— TARGET ____________________ >

c
crelabpn=e1abpntar
crex0=0.
crezO=O.
crenp=nptar
crenn=nntar
cremodir=1.
creflag=1

c

call makenucleus ()

The same procedure for the projectile:

c
c Initialize projectile
c
write(1,*) ’ ---------------- PROJECTILE —————————————————— ’

 

 

102

crelabpn=e1abpnpro
crex0=prox0
crezO=pron
crenp=nppro
crenn=nnpro
cremodir=-1.
creflag=-1

call makenucleus()
Now the following particle data should have been initialized:

rq, pq, eq, tq, fq and nnq

Also nq should give the correct number of particles per test run.

The next thing is to initialize the other ﬂags that are independent of the initialization mechanism:

c Set some flags
call initflag()

These are
last,qqznuipointscat
The charge-array is initialized in the subroutine initcharge:

c Initialize chargearray
call initcharge()

As small checks charge and mass of the two nuclei are calculated from the particle data and compared
to values calculated from the parameter ﬁle values:

c Small checks
call initcheck()

The permutation arrays perm, permx, permy, permz, deltx, delty and deltz are initialized:

c Initialize permutations
call initperm()

The cross section tables xsec. . . are initialized:

c Initialize cross section table
call initcross()

The statistics with maxsigma, minsigma, timefault, . . . are reset:

c Reset statistics
call restat()

103

The initial gridspacings are set:

c
gridspxy=2.*sqrt(maxsig/pi)

c
write(1,*) ’Gridspacing xy [fm]:’,gridspxy
write(1,*) ’(Gridextension xy +/-’,gridspxy*mgridxy,’ fm)’

c

The gridspace in z-direction however should be at least 20 times the initial spacing between the
nuclei:

gridspz=2.*timstp
if (gridspz*mgridz.le.20.*abs(pron)) gridspz=20.*pron/mgridz

write(1,*) ’Gridspacing z [fm]:’,gridspz
write(1,*) ’(Gridextension z +/-’,gridspz*mgridz,’ fm)’
For plots the z-spacing gridz will be used:

write(1,*) ’*** This gridextension will be set for all plots ***’
gridz=gridspz

Finally the timestep number nt ime is set to zero:

c
c Time step zero
c
ntime=0
c
write(1,*) secnds(t0),’: Init finished’
write(1,*)
& 2 __________________________________________________________ 2
c
return
end

B.3.2 Subroutine makenucleus

This subroutine is calling the initialization method chosen by the user from the parameter ﬁle. It
ﬁrst prints the initialization values to the protocol ﬁle.

write(1,*) ’Nucleus initialization, method ’,initmeth,’ :’

write(1,*) ’ Kin. energy per nucleon, lab. [GeV]:’,crelabpn
write(1,*) ’ Position x in lab frame [fm] :’,crex0
write(1,*) ’ Position 2 in lab frame [fm] :’,crezO
write(1,*) ’ Number of protons :’,crenp
write(1,*) ’ Number of neutrons :’,crenn
write(1,*) ’ Momentum direction (forward/backw.):’,cremodir

write(1,*) ’ Parton flag :’,creflag

104

Now the call of the chosen initialization method is performed. This list has to be expanded as soon
as a new initialization method is to be implemented:

if (initmeth.eq.1) call initmeth1()

B.3.3 Subroutine initflag

In this subroutine the last scattering partner is set to ’itself’, and the last scattering point is set to
’nowhere’:

write(1,*) secnds(t0),’: Setting flags and charges’
c For all test runs
do 10 test=1,ntr

c
c For all existing particles in test run

do 20 i=1,nq(test)
c Last scattering with itself
1ast(test,i)=i
c Not scattered before
pointscat(test,i)=-1
20 continue

Also, for unused entries of the particle arrays the entry last is set to ’-1’. This is done for creating
and deleting particles, see B.12 on page 149.

c
c Set remaining particle flags to ’not existing’
c These entries are flaged ’last=-1’

c
do 30 i=nq(test)+1,mquark
last(test,i)=-1
30 continue
c
c Number of particles
c
lastnumber(test)=nq(test)
c
10 continue
c
return
end

B.3.4 Subroutine initcharge

In this subroutine the array charge is initialized:

charge(0 )= 0.
charge(1 )= 2./3.

 

105

charge(-1)=-2./3.
charge(2 )=-1./3.
charge(-2)= 1./3.
charge(3 )= 2./3.
charge(-3)=-2./3.
charge(4 )=-1./3.
charge(-4)= 1./3.
charge(5 )= 2./3.
charge(-5)=-2./3.
charge(6 )=—1./3.
charge(—6)= 1./3.

B.3.5 Subroutine initcheck

This subroutine calculates the total charge and the total mass of all particles and compares it to the
values being calculated from the data in the parameter ﬁle:

write(1,*) secnds(t0),’: Checking charge and mass’
do 10 test=1,ntr
qchar=0.
mass=0.
do 200 i=1,nq(test)
qchar=qchar+charge(tq(test,i))
mass=mass+sqrt(eq(test,i)**2

& —pq(test,i,1)**2-pq(test,i,2)**2-pq(test,i,3)**2)
200 continue
c
write(1,*) ’Check charge :’,qchar,’ (’,nptar+nppro,’)’
write(1,*) ’Check mass :’,mass,’ (’,
& (nptar+nppro+nntar+nnpro)*manuc,’)’
10 continue

B.3.6 Subroutineinitperm

This subroutine initializes the permutation arrays by setting each entry to its own index number.

do 400 i=1,mquark
perm(i)=i
400 continue
do 410 i=-mgridxy,mgridxy
permx(i)=i
permy(i)=i
410 continue
do 420 i=—mgridz,mgridz
permz(i)=i
420 continue
do 430 i=1,3
deltx(i)=2-i
delty(i)=2-i
deltz(i)=2—i
430 continue

 

106

B.3.7 Subroutine initcross
In this subroutine the cross section arrays xsec and the maximum cross section maxsig are initialized.

The cross sections are calculated on a mesh with meshpoints eps+sqstep* (i-1) **4, where i is the
index of the mesh point. All together the mesh has maxpoints points.

integer i,j,k

c
write(1,*) secnds(t0),’: Initializing cross sections’
c
write(1,*) ’Max. energy [GeV**2]:’,
& sqstep*(maxpoints-1.)**4
c
c Maximum total cross section to be found:

maxsig=0.

The cross section arrays are ﬁrst set to the corresponding partial cross sections. The variables tin1,
tin2 and tout1 are used by the cross section routines to determine the on-shell masses of the partons
involved.

The big loop runs over all mesh point indices:

c
c First calculate cross sections
c

do 20 i=1,maxpoints

The energy corresponding to the index is calculated using the above formula:

tpsq=eps+(i-1.)**4*sqstep
call tsiggg22()
xsecgg22(i)=tpsigma

The next loop runs over all possible incoming quark types:

do 20 tin1=1,maxtype
call tsigqq22()
xsecqq22(tin1,i)=tpsigma
if (tpsigma.gt.maxsig) maxsig=tpsigma
call tsigqg22()
xsecqg22(tin1,i)=tpsigma
if (tpsigma.gt.maxsig) maxsig=tpsigma
call tsigqu22()
xsecqu22(tin1,i)=tpsigma
call tsigqugg()
xsecqugg(tin1,i)=tpsigma

For tsigggqu the varible tout1 has to be set according to the types of the outgoing partons:

tout1=tin1

call tsigggqu()
xsecggqu(tout1,i)=tpsigma

107

In the following cross sections two different quark types are involved. For tsigqqp22 those two
different quarks are in the incoming channel, so also tin2 has to be set, for tsigquqqup one
quark type has to be set for the incoming channel and one for the outgoing channel.

20

do 20 j=1,maxtype

if (tin1.ne.j) then
tin2=j
call tsigqqp22()
xsecqqp22(tin1,tin2,i)=tpsigma
tout1=j
call tsigquqqup()
xsecquqqup(tin1,tout1,i)=tpsigma

else
xsecqqp22(tin1,j,i)=0.
xsecquqqup(tin1,j,i)=0.

endif

continue

The partial cross sections at the minimum and maximum energy of the mesh are written to the
protocol ﬁle:

C

c Control output

C

30

4O

50

60

7O

write(1,*) ’-’
write(1,*) ’Parton label (6, 5, 4, 3, 2, 1, O, -1, ...... , -6)’
write(1,*) ’ for (t, b, c, s, d, u, g, u_bar, ..., t bar)’
write(1,*) ’-’
write(1,*) ’Elastic channels at min/max energies [fm**2]:’
write(1,*) ’ ------------------- -- ---’
write(1,*) ’gg->gg’,xsecgg22(1),xsecgg22(maxpoints)
write(1,*) ’qq->qq’
do 30 i=1,maxtype

write(1,*) i,xsecqq22(i,1),xsecqq22(i,maxpoints)
continue
write(1,*) ’qg->qg’
do 40 i=1,maxtype

write(1,*) i,xsecqg22(i,1),xsecqg22(i,maxpoints)
continue
write(1,*) ’qu->qu’
do 50 i=1,maxtype

write(1,*) i,xsecqu22(i,1),xsecqu22(i,maxpoints)
continue
write(1,*) ’qqp->qqp’
do 60 i=1,maxtype

do 60 j=1,maxtype

write(1,*) i,j,xsecqqp22(i,j,1),xsecqqp22(i,j,maxpoints)

continue
write(1,*) ’Inelastic channels at min/max energies [fm**2]:’
write(1,*) ’ ----------------------------------------------- ’
write(1,*) ’qu->gg’
do 70 i=1,maxtype

write(1,*) i,xsecqugg(i,1),xsecqugg(i,maxpoints)
continue

 

108

write(1,*) ’gg->qu’
do 80 i=1,maxtype
write(1,*) i,xsecggqu(i,1),xsecggqu(i,maxpoints)
80 continue
write(1,*) ’qu->qqup’
do 90 i=1,maxtype
do 90 j=1,maxtype
write(1,*) i,j,xsecquqqup(i,j,1),
& xsecquqqup(i,j,maxpoints)
90 continue

Now for a given combination of incoming partons the partial cross sections for all combinations of
outgoing partons are added up in a certain sequence. Finally the entry for a certain partial cross
section is the cross section in question plus the sum of all preceeding cross sections in that sequence.

For the incoming channel 99 the open output channels are 99 or ch. The partial cross sections
for two incoming gluons are added up in the order 99 -—> 99 (xsecgg22), 99 ——> ml, 99 —> dd, ...,

gg-—)tf(xsecggqu).

For example in xsecggqu(3,en_ergy) one will ﬁnd the partial cross section for gg —> 35 plus the
partial cross sections for 99 —> dd, 99 —) mi and 99 —> 99.

For the incoming channels qg and qq only the elastic channel is open.

For the incoming channel qc'j the partial cross sections are added up in the following order: qrj —) qrj
(xsecqu22), qr‘j —) gg (xsecqugg), qrj —-> q’ (j’ (xsecq_qbqqup). The cross section array entries for
qr‘j —> q’ q" are added up in the order qr] —) ml, ch —> dd and set to zero when q = q’.

As a ﬁnal example in xsecquqqup(2 ,4,e_nergy) one will ﬁnd the partial cross section_for dd ——> bf)
plus the following_partial cross sections: dd —> as and dd —> mi from xsecquqqup, dd —> 99 from
xsecqugg and dd —) dd from xsecqu22.

As one can see, the partial cross sections always add up to the total cross section in the last entry
of that sequence.

c

c Now pile them up

c

c Elastic channels first, then inelastic
c

do 100 i=1,maxpoints
xsecggqu(1,i)=xsecggqu(1,i)+xsecgg22(i)
do 110 j=2,maxtype
xsecggqu(j,i)=xsecggqu(j,i)+xsecggqu(j-1,i)
if (xsecggqu(j,i).gt.maxsig) maxsig=xsecggqu(j,i)
110 continue
do 120 j=1,maxtype
xsecqugg(j,i)=xsecqugg(j,i)+xsecqu22(j,i)
xsecquqqup(j,1,i)=xsecquqqup(j,1,i)+xsecqugg(j,i)
do 120 k=2,maxtype
xsecquqqup(j,k,i)=xsecquqqup(j,k,i)+

& xsecquqqup(j,k-1,i)
if (xsecquqqup(j,k,i).gt.maxsig)
& maxsig=xsecquqqup(j,k,i)
1 20 cont inue

1 00 continue

 

109

The total cross sections, again at minimum and maximum mesh energy, are given out to the protocol
ﬁkx

c
c Control output
c
write(1,*) ’Total cross section at min/max energies [fm**2]:’
write(1,*) ’ ------------------------------------------------ ’
write(1,*) ’gg->anything’,
& xsecggqu(maxtype,1),xsecggqu(maxtype,maxpoints)
write(1,*) ’qu—>anything’
do 200 i=1,maxtype
write(1,*) i,xsecquqqup(i,maxtype,1),
& xsecquqqup(i,maxtype,maxpoints)
c
200 continue
c
write(1,*) ’Maximum total cross section found [fm**2]:’,maxsig
c
return
end

B.3.8 Subroutine rancor

This subroutine randomly chooses as, y and 2 such that
rr+f+¢2gl,

0323.31.231.

subroutine rancor(x,y,z)
double precision x,y,z

10 x=1.-2.*ran1(iseed)
y=1.-2.*ran1(iseed)
z=1.-2.*ran1(iseed)
if (x*x+y*y+z*z.gt.1.) goto 10

return
end

B.3.9 Subroutine restat

This subroutine resets the statistics values being written to the protocol ﬁle whenever there is a ’big’
output.

c

c Reset min/max

c
maxsigma=0.
minsigma=maxsig

110

maxsq=0.
minsq=1000000000.
mincos=1.
maxcos=-1.
c
c Reset multistage and credel statistics
c
do 300 i=1,ntr
nsta0(i)=0
nsta1(i)=0
nsta2(i)=0
nsta3(i)=0
nsta4(i)=0
crenumber(i)=0
delnumber(i)=0
300 continue

The scattering channel statistics chscat are reset. In this array scattering events in the different
channel combinations will be counted:

do 310 i=1,maxchanne1
do 310 j=1,maxchannel
chscat(i,j)=0
310 continue

The angular distribution bins are emptied:

do 410 i=1,maxchannel
do 410 j=1,6
anglestat(i,j)=0
410 continue

do 420 j=1,6
labangle(j)=0
420 continue
The impact parameter statistics are reset:
do 500 i=0,20
impact(i)=0
500 continue

Finally, the number of timefaults is set to zero:

timefau1t=0

B.4 INITMETHl

Gerd Kortemeyer
Last revision of the ﬁle: 01/20/94

111

This subroutine is an example for what requirements an initialization routine for nuclei must meet.
It represents a very simple version that only initializes the valenz quarks and can be used as a primer
when a more sophisticated initialization routine is to be written.

In general, an initialization routine must set the following variables:

rq, Pq, eq, tq, fq and nnq

Also nq must give the correct number of particles per test run.

B.4.l Local variables

In addition to stdec . for the following variables are deﬁned:

c Radius of nucleus in xy and 2
double precision radxy, radz
c Radius of nucleon in 2
double precision radnucz
c Momentum per nucleon
double precision monuc
c Total energy of nucleon in lab frame
double precision etlabn
c Momentum per quark
double precision moqua
c Energy per quark
double precision eqqua
c Beta and gamma for transformation from lab to rest
double precision beta, gamma
c Test run
integer test
c Nucleon counter
integer ncount
c do variables
integer i,j
c Coordinates
double precision x,y,z,xc,yc,zc

B.4.2 Initialization of the particles

First, it is important that the initialization routine clearly writes to the protocol ﬁle which kind of
initialization is chosen:

C
write(1,*) '>>> ___________________________________ :
write(1,*) ’>>> Initializing Valenz-Quarks only ...’
write(1,*) ’>>> ................................... ,
C

Now a few kinematical calculations are performed:

112

c =Ekin+m“2

etlabn=crelabpn+manuc

write(1,*) ’Total energy per nucleon, lab frame [GeV]:’,etlabn
c E‘2-p“2=m“2
c p“2=E“2-m“2

monuc=cremodir*sqrt(etlabn*etlabn-manuc*manuc)

write(1,*) ’Total momentum per nucleon, lab [GeV/c]:’,monuc
c beta=p/E

beta=monuc/et1abn

write(1,*) ’Beta for transformation from rest frame :’,beta
c gamma=1lsqrt(1-beta‘2)

gamma=1./sqrt(1.-beta*beta)
c r=r0*A“(1/3)

radxy=r0*(crenp+crenn)**(1./3.)

write(1,*) ’Radius of nucleus xy [fm]:’,radxy
c rz=r/gamma

radz=radxy/gamma

write(1,*) ’Radius of nucleus 2 (contracted) [fm]:’,radz
c
write(1,*) ’Radius of nucleon xy [fm]:’,initrnuc
radnucz=initrnuc/gamma
write(1,*) ’Radius of nucleon z (contracted) [fm]:’,radnucz
write(1,*) ’---’
c
eqqua=et1abn/3.
write(1,*) ’Energy per quark [GeV]:’,eqqua
c
moqua=monuc/3.
write(1,*) ’Momentum per quark [GeV/c]:’,moqua

B.4.3 Creating particles

Now for all testruns the particles are created:

c
c For all testruns
c
do 10 test=1,ntr
write(1,*) ’ ----- Testrun : ’,test
c Reset nucleon counter
ncount=0

From here on the initialization methods might differ. First, all protons are created, their constituents
are just being thrown into one and the same arrays — only the nucleon number nnq still identiﬁes
which nucleon they belong to, while the flag fq identiﬁes which of the two nuclei, target or projectile,
their nucleon belongs to:

c
c Do the protons

do 20 i=1,crenp
ncount=ncount+1

113

c Where to place nucleon?
call rancor(x,y,z)
xc=crex0+x*(radxy-initrnuc)
yc= y*(radxy—initrnuc)
zc=crezO+z*(radz -radnucz)
c xc, yo and zc are the center coordinates of the nucleon
c Now put 3 valenz quarks into nucleon
do 30 j=1,3
c One particle more
nq(test)=nq(test)+1
c Position
call rancor(x,y,z)
rq(test,nq(test),1)=xc+x*initrnuc
rq(test,nq(test),2)=yc+y*initrnuc
rq(test,nq(test),3)=zc+z*radnucz
c Momentum
pq(test,nq(test),1)=0.
pq(test,nq(test),2)=0.
pq(test,nq(test),3)=moqua
c Energy
eq(test,nq(test))=eqqua
c Flag
fq(test,nq(test))=cref1ag
c Nucleon number
nnq(test,nq(test))=ncount
30 continue
0 Now make them up and down
tq(test,nq(test)-2)=1
tq(test,nq(test)-1)=2
tq(test,nq(test) )=1
20 cont inue

The creation of the neutrons of course is similar except for the valenz quark types:

c:
0 Do the neutrons

do 40 i=1,crenn
ncount=ncount+1
C Where to place nucleon?
call rancor(x,y,z)
xc=crex0+x*(radxy-initrnuc)
yc= y*(radxy-initrnuc)
zc=crezO+z*(radz -radnucz)
C 3x<=, yc and zc are the center coordinates of the nucleon
C NOW put 3 valenz quarks into nucleon
do 50 j=1,3
C One particle more
nq(test)=nq(test)+1
C Position
call rancor(x,y,z)
rq(test,nq(test),1)=xc+x*initrnuc
rq(test,nq(test),2)=yc+y*initrnuc
rq(test,nq(test),3)=zc+z*radnucz

114

c Momentum
pq<test.nq(test).1)=o.

qutest.nq(test),2)=o.
qutest,nq(test),3)=moqua

c Energy
eq(test,nq(test))=eqqua

c Flag
fq(test,nq(test))=creflag

c Nucleon number
nnq(test,nq(test))=ncount

50 continue

c Now make them up and down
tq(test,nq(test)-2)=2
tq(test,nq(test)-1)=1
tq(test,nq(test) )=2

40 continue

c
10 continue

B.4.4 The end

Finally as a small check for the last testrun the number of created nucleons is compared to the

expected value:

c
write(1,*) ncount,’(’,crenp+crenn,’) nucleons.’
c
return
end

B.5 MOVE

Gerd Kortemeyer

Last revision of the ﬁle: 04/04/93
This subroutine moves all particles according to their momenta and also initializes the pointer
structure for the assignment of the particles to the boxes of the spatial grid.

B .5.1 Local variables

C Do variables
integer i,j,k,m,nx,ny,nz
C Permutation subset with numbers only from 1 to nq(test_run_number)
integer subperm(mquark)
C Array with current endparticles of gridboxes
integer endpart
.& (~mgridxyzmgridxy,-mgridxy:mgridxy,-mgridz:mgridz)
C Velocity 3 vector

real v(3)
c Counter for particles having left boxes

115

integer outb
c Maximum angle, angle
double precision mtheta,theta

B.5.2 Determining the new gridsize

In order to adjust the gridsize during runtime the greatest angle between the direction any particle
is moving in and the z-direction has to be determined:

c
c Determine new gridspacing
c

mtheta=0.
c

do 100 i=1,ntr
do 100 j=1,nq(i)
if (last(i,j).ne.-1) then
c Get angle

theta=abs(
& atan(sqrt(pq(i,j,1)*pq(i,j,1)+pq(i,j,2)*pq(i,j,2))/
& pq(i.j,3)))
if (theta.gt.mtheta) mtheta=theta
endif

100 continue

With theta: dmax being the maximum angle any particle has from the initial direction, the grid

 

extensron 1n 2 direction 18 adjusted to 2tstep ~ c + 2 max sm(6’max):
V r

gridspz=2.*(timstp+sqrt(maxsig/pi)*sin(mtheta))

B.5.3 Moving the particles

First the variable outb is set to zero here. It will count the particles that have left the grid area.

c
outb=0
c:

Over all there is a loop over all test runs:

c: Loop for all active test-runs:
do 10 i=1,ntr

For the respective test run the gridroots are set to ’empty box’.

C Reset gridroots and endpart
do 20 j=-mgridxy,mgridxy
do 20 k=-mgridxy,mgridxy
do 20 m=-mgridz,mgridz
gridroot(i,j,k,m)=-1

20 continue

116

In the integer array subperm the permuted numbers from 1 to nq(i) are initialized in the entries
1 to nq(i). i is the number of the test run, nq(i) is the number of particles in this run. This is
a ’subpermutation’ of the permutation perm that has the numbers 1 to mquark in the entries 1 to
mquark. subperm is designed for the random access of the nq(i) particles in the actual test run.

c Initialize permutation subset
j=1
do 30 k=1,mquark
if (perm(k).le.nq(i)) then
subperm(j)=perm(k)
j=j+1
endif
30 continue

If the subpermutation was initialized correctly, the value of j is nq(i)+1:

c Check subperm
if (j.ne.nq(i)+1) then
write(1,*) ’Sub-Permutation wrong, j=’,j

write(1,*) ’ Test-Run:’,i
write(1,*) ’ Number of particles:’,nq(i)
stop

endif

Now there is a loop over all particles in the test run, but the actual index j for the particle that
is to be moved is taken from the subpermutation. This is done in order to avoid that the chain of
pointers for the particles in the boxes has the same order in subsequent timesteps:

c Loop for all particles in that run:
do 10 m=1,nq(i)
c Particle numbers chosen from permutation array
j=subperm(m)

It then has to be checked if the particle is actually existing:

c Particle present?
if (last(i,j).ne.-1) then

See A.7 on page 84 and B.1.3 on page 87 for details.

The moving of the particle is done using

 

m + At) 3 m) + m) - At = m) + E

c Loop for all components of the 3-vector
do 40 k=1,3
c

rq(i,j,k)=rq(i,j,k)+pq(i,j,k)*timstp/eq(i,j)

40 continue

117

B.5.4 Assigning the particles to the boxes

The indices of the box the particle is in now are calculated as the coordinate divided by the respective
gridspacing.

c In which box is the particle now?

c

nx=nint(rq(i,j,1)/gridspxy)
c

ny=nint(rq(i,j,2)/gridspxy)
c

nz=nint(rq(i,j,3)/gridspz)

Now it is to be checked if the particles left the grid area. If so, they are put into the outmost box
and the counter outb is incremented.

c Did particles leave the box?

if (nx.gt.mgridxy) then
nx=mgridxy
outb=outb+1

endif

if (nx.1t.-mgridxy) then
nx=-mgridxy
outb=outb+1

endif

if (ny.gt.mgridxy) then
ny=mgridxy
outb=outb+1

endif

if (ny.lt.-mgridxy) then
ny=-mgridxy
outb=outb+1

endif

if (nz.gt.mgridz) then
nz=mgridz
outb=outb+1

endif

if (nz.lt.-mgridz) then
nz=-mgridz
outb=outb+1

endif

Now the particle certainly has a valid box index it is time to incorporate it into the pointer structure.

The pointer array endpart is used to always point to the last particle being found in the respective
box during initialization of the pointer structure, it is not needed once the initialization is ﬁnished.

If the actual particle is the ﬁrst particle found in its box both gridroot and endpart will be
initialized pointing to that particle.

If the box wasn’t empty the previously found particles pointer is set to the actual particle.
c Is it the first particle in that box?

if (gridroot(i,nx,ny,nz).lt.0) then
gridroot(i,nx,ny,nz)=j

118

endpart(nx,ny,nZ)=l
else
nextgrid(i,endpart(nx,ny,nz))=j
endif

Finally, the actual particles pointer is set to ’nowhere’ since it is not sure if there will be another
particle found. If there is one more particle in the box, the actual particles pointer will be set to it
at the time the next particle is found by the above statement nextgrid(i ,endpart . . ..

Also the endpart pointer is set to actual particle number.
c Set nextgrid to nowhere, endpart to j

nextgrid(i,j)=-1
endpart(nx,ny,nz)=j

B.5.5 The end

Now the if is closed and the label for the two loops over the test run and the particles is set:

endif
10 continue

Also a warning is given out if one or more particles left the grid:

c Write warnings
if (outb.ne.0) write(1,*)

& ’!!!’,outb,’ grid boundary violations in timestep’,ntime
c
Bye bye. . .
return
end

B.6 SCATTER

Gerd Kortemeyer
Last revision of ﬁle: 09/19/94

This subroutine is responsible for ﬁnding all particle combinations that are able to scatter. Can-
didates are delievered to subroutine collide where the physics tests for scattering are run and in
case the scattering is executed.

B.6. 1 Local variables

Besides stdec . for the following local variable deﬁnitions are made:

119

c do variables
integer i,j,k,1,m,n

c box numbers of center box
integer ni,nj,nk

c box numbers of second box
integer mi,mj,mk

c bases for linear index
integer base,basesq

c linear index of first and second box
integer linfirst,linsecond

c Data about first particle
integer 1astp1,flagp1

B.6.2 Initialization of variables

First the subroutine deﬁnes the variables

base=2*mgridxy+1
basesq=base*base

that will be used later to calculate an identiﬁcation number for the boxes. The compiler will treat
both variables like constants.

B.6.3 Choosing the boxes

For choosing the boxes ﬁrst a big loop over all test runs is executed. The loop variable is chosen to
be tprn so automatically the right test run number can be found in the common twopart.

c
c For all test runs active
c

do 10 tprn=1,ntr

The array alr, that is used to ﬂag particles that already have scattered in this timestep, is set to
zero, meaning ’not scattered’

c No scattering so far in this run
do 20 j=1,nq(tprn)
alr(tprn,j)=0
20 continue

Next, all possible center boxes are looked at. The mapping from the do-variables i , j ,k to the box
indizes ni ,nj ,nk is done by means of the permutation arrays permx,permy,permz:

c Choose center box, not including outmost boxes
do 30 i=-mgridxy+1,mgridxy-1
do 30 j=-mgridxy+1,mgridxy-1
do 30 k=-mgridz+1,mgridz-1
ni=permx(i)
nJ' =permy (j >
nk=permz(k)

120

The outmost boxes are excluded from the loop. The mapping through permx, permy and permz
generates no problems here since it is taken care in subroutine permute that the indices for the
outmost boxes are not permuted.

Now gridroot (tprn,ni,nj ,nk) points to the ﬁrst particle being found in that box. However, if
the box is empty the gridroot value is ’-1’. In this case it makes no sense to use this box as the
center box and the program jumps to label 40 that immediately preceeds the loop label 30: The
100p chooses the next center box.

c Center box empty?
if (gridroot(tprn,ni,nj,nk).lt.0) goto 40

Without further effort every box combination would be checked twice: One time when it is the
center box, the other time when it is the neighbor box of another center box. If the boxes had one
index instead of three the problem would be easy to solve: One could simply only use combinations
where the index of the ﬁrst box is smaller than or equal to the index of the second box.

Nevertheless a single linear index can be produced by the following mapping:

c Calculating linear index of first box
linfirst=(nk+mgridz)*basesq
& +(nj+mgridxy)*base
& + ni+mgridxy
2-mgridxy+1 is chosen as base, basesq is basez. The adding of mgridz and mgridxy is done to
make for example nj +mgridxy range from 0 to 2-mgridxy.

Now a loop for all neighboring boxes is set up. All together there are 27 box combinations to be
checked. In order to do this in a random sequence the mapping via deltax ,deltay ,deltaz is used,
the indizes for the second box are mi ,mj ,mk afterwards:

c Second box
do 50 l=1,3
do 50 m=1,3
do 50 n=1,3
mi=ni+deltx(l)
mj=nj+delty(m)
mk=nk+deltz(n)

If the second box is empty the program jumps to label 60 that immediately preceeds the loop label
50 for the second box.

c Second box empty?
if (gridroot(tprn,mi,mj,mk).1t.0) goto 60

For this box also a linear index is calculated, afterwards the indizes are compared. If the ﬁrst index
is bigger than the second index it is jumped to label 60 for the next second box:

c Calculating linear index of second box
1insecond=(mk+mgridz)*basesq
& +(mj+mgridxy)*base
& + mi+mgridxy
c Now we have two non-empty boxes with indices ni,nj,nk and mi,mj,mk.
c As you don’t want to check box combinations twice, check only the
c upper diagonal.
if (linfirst.gt.linsecond) goto 60

121

B.6.4 Combination of all particles from the boxes

As already explained the particles from the boxes are accessed by following a chain of pointers
starting with the gridroot. For the actual particle number directly the variables tppl and tpp2
from the common block twopart are used. The label 80 is used for a handmade loop that runs over
all tpp1, the label 70 is used for a loop over all tpp2.

The labels 90 for the ﬁrst particle and 100 for the second particle are used to skip the respective
particle and go on with the next one.

c Combine all particles from one box with the other boxes particles.
tpp1=gridroot(tprn,ni,nj,nk)
80 continue

Now the formal reasons against scattering explained in A.3 on page 81 have to be checked. As
already explained in the array alr for particles not scattered in this timestep one can ﬁnd ’0’, as
soon as they have scattered ’1’.

c Particle 1 has already scattered - next one!
if (alr(tprn,tpp1).ne.0) goto 90

lastp1=last(tprn,tpp1)
flagp1=fq(tprn,tpp1)

c

c Start loop for second particle
tpp2=gridroot(tprn,mi,mj,mk)

c Particle loop, tpp2

70 continue

c Particle 2 already scattered? Next particle!

if (alr(tprn,tpp2).ne.0) goto 100

In the array fq for all particles that never have scattered one can ﬁnd ’-1’ if they originally belonged
to the projectile, and ’1’ if they originally belonged to the target. Scattered particles have fq=0.

c Same nucleus?
if (flagp1*fq(tprn,tpp2).eq.1) goto 100

In the array last the number of the last scattering partner of the particles can be found.

c Particle 2 scattered with 1 recently?

if (tpp2.eq.lastp1) goto 100
c Particle 1 scattered with 2 recently?

if (tpp1.eq.last(tprn,tpp2)) goto 100
c

If all this requirements are met the subroutine collide is called.

call collide()
if (alr(tprn,tpp1).ne.0) goto 90

 

122

After the particle combination is checked, the next particles are found by the latter particles’
nextgrid pointers. It is important that the procedures for creating and deleting particles which
might be called within collide do not change the nextgrid-pointers of the chosen particles, see
B.12 on page 149. If the nextgrid-pointers have the value ’-1’ no more particles are in the respective
box:

 

 

100 tpp2=nextgrid(tprn,tpp2)
if (tpp2.ge.0) goto 70
c
90 tpp1=nextgrid(tprn,tpp1)
if (tpp1.ge.0) goto 80
c
c Second box
60 continue
50 continue
C >
c Center box
40 continue
30 continue
c Test-Runs
10 continue
Goodbye:
return
end

B.7 CAUSCATTER

Gerd Kortemeyer
Last revision of ﬁle: 05/04/94

This subroutine is responsible for ﬁnding all particle combinations that are able to scatter. Can-
didates are delievered to subroutine collide where the physics tests for scattering are run and in
case the scattering is executed.

It is a blueprint of the subroutine scatter, see B.6 on page 118 with the exception, that collisions
that would lead to signal velocities faster than the speed of light are forbidden.

This ﬁle can be linked instead of scatter from the main ﬁle.

The differences are:

B.7 .1 Local variables

c Scattering midpoint, velocity
double precision midp(3) ,velo

123

B.7 .2 Combination of all particles from the boxes

As the last check before the collision routine is called, the signal speed is calculated for both partons
involved. If one of the signals would be too fast, the collision is forbidden.

The collision is assumed to take place in the middle between the partons:

c Okay, midpoint calculation
midp(1)=(rq(tprn,tpp1,1)+rq(tprn,tpp2,1))/2.
midp(2)=(rq(tprn,tpp1,2)+rq(tprn,tpp2,2))/2.
midp(3)=(rq(tprn,tpp1,3)+rq(tprn,tpp2,3))/2.

For the actual calculation the parton data entries lastpoint and 1asttime are used. In them the
midpoint coordinates of the most recent collision the parton was involved in are stored.

c Faster than the speed of light ... ???
c
c < ----
if (nnq(tprn,tpp1).eq.0) then
velo=sqrt( (midp(1)-lastpoint(tprn,tpp1,1))**2

& +(midp(2)-lastpoint(tprn,tpp1,2))**2
& +(midp(3)-lastpoint(tprn,tpp1,3))**2 ) /
& ( (ntime - lasttime(tprn,tpp1))*timstp )
if (velo.gt.1.) goto 100
endif
c
if (nnq(tprn,tpp2).eq.0) then
velo=sqrt( (midp(i)-lastpoint(tprn,tpp2,1))**2
& +(midp(2)-lastpoint(tprn,tpp2,2))**2
& +(midp(3)-lastpoint(tprn,tpp2,3))**2 ) /
& ( (ntime - lasttime(tprn,tpp2))*timstp )
if (velo.gt.1.) goto 100
endif
c -——- >
c
c _______________________________________________________________________

c Collide them!

B.8 COLLIDE

Gerd Kortemeyer
Last revision of ﬁle: 08/17/94

This subroutine decides wether or not two particles, that have been determined by the subroutine
scatter, will collide. It does this by performing a multistage check, the subroutine is left immedi-
ately if one of the requirements is not met. If the two particles scatter the output channel is chosen
and a corresponding other routine called.

124

B.8.l Subroutine collide

Local variables

The ﬁle stdec .for is included.

c Scalar product rp in cm-frame
double precision rpcm1,rpcm2
c Distance vector cm-frame one timestep later
double precision dvcm(3)
c Transformed time step
double precision cmtimstp
c Factors in transformation formula
double precision scalbetax,a

Checks in the lab frame

To check the multistage statistics everytime the routine is entered nstaO is increased:

c
c Statistics: entered routine
c

 

nstaO(tprn)=nsta0(tprn)+1

The ﬁrst checks take place in the lab (or ’computer’) frame. The distance vector tprl is calculated
as the difference of the two particles coordinate vectors. tpdq is the square of the difference vector.

c
c Calculate distance vector lab and distance lab
C ——————————————————————————————————————————————
c
tprl(1)=rq(tprn,tpp1,1)-rq(tprn,tpp2,1)
tprl(2)=rq(tprn,tpp1,2)-rq(tprn,tpp2,2)
tprl(3)=rq(tprn,tpp1,3)-rq(tprn,tpp2,3)
c

tpdq=tprl(1)*tpr1(1)+tprl(2)*tpr1(2)+tprl(3)*tprl(3)

Statistics for minimum and maximum distance in the lab frame are being set:

if (tpdq.le.mindq) mindq=tpdq
if (tpdq . ge .maxdq) maxdq=tpdq

It is of great importance that the spatial distance in the c.m.-frame is always bigger than distance
in lab-frame! This is because the time component of the distance is zero in the lab frame.

The spatial distance in the c.m.-frame is also always bigger than or equal to the impact parameter,
because the impact parameter vector is only a projection of the distance vector in the c.m.-frame.
But in case of a possible scattering event the point of closest approach has to be reached as another
condition, so if scattering is possible, the c.m.-distance equals the impact parameter and ﬁnally the
distance in the lab frame will be smaller than the impact parameter.

125

If however the distance in the lab frame is bigger than the maximum impact parameter given by
the maximum cross section maxsig, that is determined during the initialization of the cross section
table, no scattering is possible:

 

c

c

c lst check !!
c --------- !!
c Distance in c.m.-frame is always bigger than distance in lab-frame! !!
c Distance in c.m.-frame always bigger or equal impact parameter. !!
c Scattering only when closest approach in c.m.-syste ==> scatter !!
c only when c.m.-distance equals impact parameter. !!
c !!

if (tpdq.gt.maxsig/pi) return !!
If the ﬁrst test was passed, the statistics for this test are incremented:

c Increase statistics: Passed first test in this test run !!
nstal(tprn)=nsta1(tprn)+1 !!

 

Finding the boost velocity

The transformation between c.m.- and lab frame are given by the following equations:

Ecm = 7(Elab - gﬁab) (13:1)
_. _. 7 — 1 —9_. a —o

Pcm = Plab + 77(ﬁp1abm - 7 E 03-?)

Elab = 7(Ecm + 35cm) (B3)
-0 -o 7 _ 1 #4 —o —o

plab = pcm + _,52 (,chmlﬁ ‘l‘ ’YﬁEcm (BA)

The only question is what 6' for the intended transformation is. The c.m.-frame is deﬁned by
51cm = ‘52cm a
51cm and 52cm being the momenta of the two particles in the c.m.-frame.

This leads to

3 = ﬁilab + 17sz
Ellab + E2lab

 

(B.5)

In tpeql the sum of the lab energies is stored, in tpbeta the vector 6, tpbq is set to 62.

c

c Calculate beta-vector for transformation from lab-frame to cm.
tpeql = eq(tprn,tpp1) +eq(tprn,tpp2)
tpbeta(1)=(pq(tprn,tpp1,1)+pq(tprn,tpp2,1))/tpeql
tpbeta(2)=(pq(tprn,tpp1,2)+pq(tprn,tpp2,2))/tpeql
tpbeta(3)=(pq(tprn,tpp1,3)+pq(tprn,tpp2,3))/tpeql

c Calculate beta squared
tpbq=tpbeta(1)**2+tpbeta(2)**2+tpbeta(3)**2

126

For the total energy squared in the c.m.-frame tpsq equation (BB) is used. There the vector ﬁcm is
the total momentum of the two particles in the c.m.-frame and therefore zero.

Division by '7 results in

E2
52 :2 E2", = jab = Efaba — 52) . (B.6)

 

2

c Calculate total energy c.m.
tpsq=tpeql*tpeql* (1 . -tpbq)

A few statistics are set for the total energy squared:

c Statistics

c
if (tpsq.lt.minsq) minsq=tpsq
if (tpsq.gt.maxsq) maxsq=tpsq

Further tests in the lab frame

After the distance in the lab frame had been compared to the maximum assumed cross section now
the comparision to the real cross section for the process is performed. The index tpcrindex is used
for looking up the cross section from the cross section arrays, see subroutine initcross in B.3.7 on
page 106.

 

sqstep is set in stdec, so is maxpoints as the maximum lookup-table index.

The rounding is done by truncation. This is important as otherwise it could happen that if the
energy-limit below which the collision is kinematically forbitten is just between the two indizes in
question, and rounding up leads to a forbidden collision.

c
c Calculate index for lookup table
c
tpcrindex=int(sqrt(sqrt((tpsq-eps)/sqstep)))+1
if (tpcrindex.gt.maxpoints) tpcrindex=maxpoints
c

The variables tinl and tin2 are set to the particle types of the incoming particles:

tin1=tq(tprn,tpp1)
tin2=tq(tprn,tpp2)

The total cross section tpsigma is looked up from the cross section tables, also the incoming channel
tpchinisseu

c Get total cross section
call totalcross()

The statistics for minsigma and maxsigma are set:

c
c Statistics
c
if (tpsigma.1t.minsigma) minsigma=tpsigma
if (tpsigma.gt.maxsigma) maxsigma=tpsigma

 

127

Now the distance is compared to the real cross section:

 

c 2nd check !!
c --------- !!
if (tpdq.gt.tpsigma/pi) return !!
c Increase statistics: Passed second test in this test run !!
nsta2(tprn)=nsta2(tprn)+1 !!

 

The transformation of the momenta to the c.m.-system

t amma is set to y = —1—. For the transformation the values of 37—31 and Eff“ are calculated
P8 VG:E§ g ab

in tpgambet and scalbetax respectively.

In case the two paricles already are nearly in their c.m.-system (lﬁl << 1), the non-relativistic limit
of the transformation is used:

c Now for the transformation to c.m. system!

c Calculate gamma
tpgamma=1./sqrt(1.-tpbq)

if (tpbq.gt.eps) then
tpgambet=(tpgamma-1.)/tpbq
c Transform momentum 1
scalbetax=tpbeta(1)*pq(tprn,tpp1,1)+
& tpbeta(2)*pq(tprn,tpp1,2)+
& tpbeta(3)*pq(tprn,tpp1,3)

The variable a is only an appreviation and can be used in equation (B.2) when 6 is written in front
of the two last terms:

a=tpgambet*scalbetax-tpgamma*eq(tprn,tppl)

Now for equations (B.1). . .

tpeql =tpgamma*(eq(tprn,tpp1)-scalbetax)
- . and (B.2):
tpmo1(1)=pq(tprn,tpp1,1)+a*tpbeta(1)

tpmol(2)=pq(tprn,tpp1,2)+a*tpbeta(2)
tpmol(3)=pq(tprn,tpp1,3)+a*tpbeta(3)

128

Now the same transformation could be done for particle 2, but the result of the momentum trans-
formation is clear: -—f)'2cm.

c Transform momentum 2
scalbetax=tpbeta(1)*pq(tprn,tpp2,1)+
& tpbeta(2)*pq(tprn,tpp2,2)+
& tpbeta(3)*pq(tprn,tpp2,3)

tpeq2 =tpgamma*(eq(tprn,tpp2)-scalbetax)

The momentum of course is not transformed.

The else branch of the condition does a non-relativistic transformation:

else
c
c Nonrelativistic transformation as limit for low beta
c
scalbetax=0.
tpeql =eq(tprn,tpp1)
tpm01(1)=(pq(tprn,tpp1,1)-pq(tprn,tpp2,1))/2.
tpmol(2)=(pq(tprn,tpp1,2)-pq(tprn,tpp2,2))/2.
tpmo1(3)=(pq(tprn,tpp1,3)-pq(tprn,tpp2,3))/2.
tpeq2 =eq(tprn,tpp2)
endif

 

C

A small check is done to the transformation by examining if the sum of the two particle energies
really equals the total energy previously calculated. If not, the program is aborted with a report
written to the protocol ﬁle:

c Small check
if (abs((tpeq1+tpeq2)**2-tpsq)/tpsq.gt.0.01) then
write(1,*) ’ERROR in Lorentztransformationz’
write(1,*) (tpeq1+tpeq2)**2,tpsq
write(1,*) ’1 : ’,
& eq(tprn,tpp1),pq(tprn,tpp1,1),Pq(tprn,tpp1,2),pq(tprn,tpp1,3)
write(1,*) ’Ml: ’,
sqrt(eq(tprn,tpp1)**2-pq(tprn,tpp1,1)**2
& -pq(tprn,tpp1,2)**2-pq(tprn,tpp1,3)**2)
write(1,*) ’2: ’,
& eq(tprn,tpp2),pq(tprn,tpp2,1),pq(tprn,tpp2,2),pq(tprn,tpp2,3)
write(1,*) ’M2: ’,
sqrt(eq(tprn,tpp2)**2-pq(tprn,tpp2,1)**2
& -pq(tprn,tpp2,2)**2-pq(tprn,tpp2,3)**2)
write(1,*) ’tpeq1/2: ’,tpeq1,tpeq2
write(1,*) ’tpmol :’
write(1,*) tpm01(1),tpmol(2),tpmo1(3)
write(1,*) ’M1tz’,

8?

8?

& sqrt(tpeq1**2-tpmol(1)**2-tpmol(2)**2-tpmol(3)**2)
write(1,*) ’M2t:’,
& sqrt(tpeq2**2-tpmo1(1)**2-tpmol(2)**2-tpm01(3)**2)
stop ’ERROR in Lorentztransformation’
endif

129

The transformation of the relative position to the c.m.-system

The distance vector of the two particles is also transformed with equations (B.1) and (B.2) used for
the position rather than the momentum. The time difference between the two particles in the lab
frame is zero. The time distance in the c.m.-frame is called tptd, the spatial distance tprc.

Again it has to be checked if a non-relativistic transformation should be used:

c
c Transform relative position
if (tpbq.gt.eps) then
scalbetax=tpbeta(1)*tprl(1)+
& tpbeta(2)*tprl(2)+
& tpbeta(3)*tprl(3)
a=tpgambet*scalbetax

tptd =tpgamma*scalbetax

tprc(1)=tprl(1)+a*tpbeta(1)
tprc(2)=tprl(2)+a*tpbeta(2)
tprc(3)=tpr1(3)+a*tpbeta(3)

In case no transformation is necessary the values are set directly:

 

else
scalbetax=0.
tptd=0.
tprc(1)=tprl(1)
tprc(2)=tprl(2)
tprc(3)=tprl(3)
endif

The transformation to the c.m.-frame is somewhat strange if afterwards the time distance between
the two particles is bigger than the transformed timestep length, called cmtimstp. cmtimstp is
calculated now by means of (8.3). . .

c
c Transform timestep-length
cmtimstp=timstp/tpgamma

. . . and compared to the time distance. If the time distance is bigger, the statistics counter timefault
is incremented:

c Time difference in c.m. frame > timestep length?
if (abs(tptd).gt.cmtimstp) timefault=timefault+1

The spatial distance squared is calculated and named tpdqc:

c Distance in c.m. frame
tpdqc=tprc(1)**2+tprc(2)**2+tprc(3)**2

130

Checks in the c.m.-frame

The condition ’c.m.-distance greater than cross section radius’ throws out more events than the
condition ’impact parameter greater than cross section radius’! Therefore it is more effective and
still correct because in the case when there really is scattering both conditions are equivalent — as
already pointed out because of the second condition ’closest approach’ the impact parameter equals
the distance.

 

3rd check 11
......... 1 1

000

The condition ’distance greater than cross section’ throws out !!
more events than the condition ’impact parameter greater than !!
cross section’! Therefore it’s more effective. !!
In the case when there is really scattering both conditions are !!

equivalent since ’impact parameter=distance’ for closest approach. !!
I I

000000

if (tpdqc.gt.tpsigma/pi) return !!
c Increase statistics: Passed third test in this test run !!
nsta3(tprn)=nsta3(tprn)+1 !!

 

C

 

Now the ’closest approach’ test is performed.

It is done by calculating the scalar product of the distance vector with the momentum vector both
now and one timestep later. At the point of closest approach itself the scalar product is zero.
However the changing of the sign of the scalar product from now to the next timestep is taken as
an indication that the point of closest approach will be reached in this timestep.

The ﬁrst scalar product can be calculated directly:

Closest approach test

 

Calculate scalar product of c . m . -distance and c . m. -momentum
(proportional to angle between them! )
rpcm1=tprc ( 1) *tpmol (1) +tprc (2) *tpmol (2) +tprc (3) *tpmo1 (3)

00000

This is followed by the propagation of the partons in the lab-frame and subsequent transformation
of the new positions to the c.m.-frame. The propagation vectors for the partons take the form

P1 b,’
(Atstep, 'E—a'l‘Atstep) -

lab,1

With Aao(a) being the spatial Lorentz-transformation into the c.m.-frame for a vector ((10, a), one
gets

plab 2
I _ :
Arcm — AO+Atstep_O—Atstep (712,1); + —E ab 1 Atstep - rlab,2 - Elab 2

)

 

Atstep) (B.7)

 

 

 

P1 b,1 P1 b,2
a , a )

131

 

 

Because of
Ao(P1ab,i) = Aliyah, (Plab,i) + 731531511); = Pcm,i + 7:6Elab,i (B-8)
the result is
Argm = Arcm + ( Pcm + Pcm ) Atstep . (B.9)
7(Ecm,1 + .chm) 7(Ecm,2 — ﬁpcm)

c Calculate distance vector one timestep later
c
dvcm(1)=tprc(1)+
& (tpmol(1)/eq(tprn,tpp1)+tpmol(1)/eq(tprn,tpp2))*timstp
dvcm(2)=tprc(2)+
& (tpmol(2)/eq(tprn,tpp1)+tpmol(2)/eq(tprn,tpp2))*timstp
dvcm(3)=tprc(3)+
& (tpmol(3)/eq(tprn,tpp1)+tpm01(3)/eq(tprn,tpp2))*timstp
c Calculate same scalarproduct as before one timestep later
rpcm2=dvcm(1)*tpmol(1)+dvcm(2)*tpmol(2)+dvcm(3)*tpmol(3)

 

The closest approach test is now simply done by checking the sign of the product of the two scalar
products:

c
c
c 4th check !!
c --------- !!

if (rpcm1*rpcm2.gt.0.) return !!
c Increase statistics: Passed forth test in this run !!

nsta4(tprn)=nsta4(tprn)+1 !!

 

 

If the particles have passed all this test, they will scatter:

Scattering the particles

c
c
c Okay, scatter
c

 

 

But ﬁrst it is made sure that the two particles really exist:

c Last check: ’last’ check
if ((last(tprn,tppl).eq.-1).or.(last(tprn,tpp2).eq.-1)) then
write(1,*) ’Trying to collide non-existing particles ...’
write(1,*) ’tprn:’,tprn
write(1,*) ’tpp1:’,tpp1,’ last:’,last(tprn,tpp1),’ nextz’,
& nextgrid(tprn,tpp1)

132

write(1,*) ’tpp2:’,tpp2,’ last:’,last(tprn,tpp2),’ nextz’,
& nextgrid(tprn,tpp2)
write(1,*) ’nq(tprn):’,nq(tprn)
call checknumbers(tprn,tppl,tpp2)
stop ’Trying to collide non-existing particles ...’
endif

 

Now the output channel is chosen by subroutine outchannel, see B.8.2 on page 132. The channel
number is stored in tpchout afterwards:

c
c Determine output channel tpcrout
call outchannel()

The subroutine xsecentry (see B.8.4, page 136) is called to check the cross section entry in the
arrays xsec. . .:

 

c __________________________________________________________
c Check cross section table !
call xsecentry() !

C ---- — — ———————————————————————————————————————

 

 

 

Since at the moment only 2 —> 2 processes are implemented, the subroutine twototwo is directly
called to complete this process. However, once 2 —> 1 processes are implemented, depending on
tpchout another routine could be called, too.

c

call twototwo()
c

return

end

B.8.2 Subroutine outchannel

The subroutine outchannel is responsible for choosing the output channel based on the relative
sizes of the partial cross sections:

c
C ____________________________________________________________________
c Determines outgoing particles
c
subroutine outchannel()
c

include ’stdec.for’
c Random number

double precision rd
c Do variable

integer i

The following reminds the channel numbers as descripted in B.1.8 on page 91:

133

Channels:

- gs

- qs

- qq
qu

- qqp
qqup

00000000
OlU'erCOMH
I

The following gives an overview of the cross section arrays:

c
c>> common /cross/ xsecgg22(maxpoints),

c>> & xsecqq22(maxtype,maxpoints),

c>> & xsecqg22(maxtype,maxpoints),

c>> & xsecqu22(maxtype,maxpoints),

c>> & xsecqqp22(maxtype,maxtype,maxpoints),
c>> & xsecqugg(maxtype,maxpoints),

c>> & xsecggqu(maxtype,maxpoints),

c>> & xsecquqqup(maxtype,maxtype,maxpoints)

 

For a description of this cross section arrays see B.3.7 on page 106.

A random number between 0 and the total cross section tpsigma is determined:

c
c Random number between 0 and total cross section
c

rd=ran1(iseed)*tpsigma
c

The variable tpchout will be set by this routine according to the output channel number determined.
It is preset to —1 to check later if a channel had been found:

tpchout=~1

The open output channels of course depend on the input channel, stored in tpchin and determined
by the subroutine totalcross, see B.8.3 on page 135:

if (tpchin.eq.1) then

Incoming channel 1 is 99, so the open output channels are gg or qr]. The partial cross sections for
two incoming gluons are added up in the order gg —+ 99, 99 —> ml, 99 —> dd, . . . , 99 —-> tf:

if (rd.le.xsecgg22(tpcrindex)) then
tpchout=1
tout1=tin1
tout2=tin2
else
do 10 i=1,maxtype
if (rd.le.xsecggqu(i,tpcrindex)) then
tpchout=4

134

tout1=i
tout2=~i
goto 30
endif
10 continue
endif
endif

Label 30 is just before ﬁnal checks and exiting the routine. The variables tinl and tin2 are set
to the particle types of the incoming partons, tout1 and tout2 are to be set by this routine to the
parton types of the outgoing partons.

Incoming channel 2 is qg, only the elastic channel is open:

if (tpchin.eq.2) then
tpchout=2
tout1=tin1
tout2=tin2

endif

Incoming channel 3 is qq, again only the elastic channel is open:

if (tpchin.eq.3) then
tpchout=3
tout1=tin1
tout2=tin2

endif

Incoming channel 4 is qq, the partial cross sections are added up in the following order: qq —> qq,
ch —> gg_, qq —-) q’ q’ . The cross section array entries for qq —> q’ q’ are added up in the order qq —> ml,
qq' —> dd and set to zero when q = q’.

if (tpchin.eq.4) then
if (rd.le.xsecqu22(abs(tin1),tpcrindex)) then
tpchout=4
tout1=tin1
tout2=tin2
else
if (rd.le.xsecqugg(abs(tin1),tpcrindex)) then
tpchout=1
tout1=0
tout2=0
else
do 20 i=1,maxtype
if (rd.1e.xsecquqqup(abs(tin1),i,tpcrindex)) then
tpchout=6
tout1=i
tout2=~i
goto 30
endif

 

135

20 continue
endif
endif
endif

For the incoming channel 5, qq’, again only the elastic outchannel is open:

if (tpchin.eq.5) then
tpchout=5
tout1=tin1
tout2=tin2

endif

Label 30 is the exit-label:

c
30 continue
c

As a small check it is made sure that actually one outchannel had been found:

if (tpchout.1t.0.) then
write(1,*) ’No output channel found’
write(1,*) tin1,tin2,tpchin
stop ’No output channel found’

endif

return
end

B.8.3 Subroutine totalcross

This subroutine determines the total cross section from the cross section arrays. It is also responsible
for determining the incoming channel tpchin. As the partial cross sections have been added up by
the subroutine initcross, see B.3.7 on page 106, the total cross section is found in the entry for
the last partial cross section in this sequence:

c ....................................................................
c Determines total cross section
c
subroutine totalcross()
c
c Channels:
6 1-ss
6 2-qs
c 3-qq
c 4 - qu
C 5 ' qqp
c 6-qqup
c
c

include ’stdec.for’

 

136

The total cross section will be written to tpsigma, it is set to ——1 ﬁrst for check purposes:
tpsigma=-1.
The variables tin1 and tin2 are set to the types of the incoming partons.

c
if (tin1.eq.tin2) then
if (tin1.eq.0) then
tpsigma=xsecggqu(maxtype,tpcrindex)
tpchin=1
else
tpsigma=xsecqq22(abs(tin1),tpcrindex)
tpchin=3
endif
else
if (tin1.eq.-tin2) then
tpsigma=xsecquqqup(abs(tinl),maxtype,tpcrindex)
tpchin=4
else
if ((tin1.eq.0).or.(tin2.eq.0)) then
tpsigma=xsecqg22(abs(tin1+tin2),tpcrindex)
tpchin=2
else
tpsigma=xsecqqp22(abs(tin1),abs(tin2),tpcrindex)
tpchin=5
endif
endif
endif

As small checks the total cross section is compared to the maximum cross section maxsig being de-
termined by the subroutine initcross, also, it is checked if the subroutine actually could determine

the cross section:

if (tpsigma.gt.maxsig) then
write(1,*) ’maxsig not maximal’
write(1,*) tin1,tin2,tpcrindex
stop ’maxsig not maximal’

endif

if (tpsigma.1t.0.) then
write(1,*) ’tpsigma not calculated’
write(1,*) tin1,tin2
stop ’tpsigma not calculated’
endif

return
end

B.8.4 Subroutine xsecentry

The subroutine xsecentry checks the validity of the cross section tables by calculating the partial
cross section for the given input and output channels and comparing them to the table entries. This

137

subroutine again reﬂects the structure of the cross section arrays, in order to get a partial cross
section for a process one has to subtract the previous entry in the cross section sequence from the
given one.

For a description of this cross section arrays see B.3.7 on page 106.

 

c
c ____________________________________________________________________
c Calculate cross section for given channels and compare to table
c
subroutine xsecentry()
c
include ’stdec.for’
c
c Channels:
6 1-ss
6 2-qs
c 3-qq
c 4 - qu
c 5-qqp
c 6-qqup
c
c TSIGGG22.FOR;1 TSIGGGQQB.FOR;1 TSIGQG22.FOR;1 TSIGQQ22.FOR;1
c TSIGQQB22.FOR;1 TSIGQQBGG.FOR;1 TSIGQQBQPQBP.FOR;1 TSIGQQP22.FOR;1
c
c>> common /cross/ xsecgg22(maxpoints),
c>> & xsecqq22(maxtype,maxpoints),
c>> & xsecqg22(maxtype,maxpoints),
c>> & xsecqu22(maxtype,maxpoints),
c>> & xsecqqp22(maxtype,maxtype,maxpoints),
c>> & xsecqugg(maxtype,maxpoints),
c>> & xsecggqu(maxtype,maxpoints),
c>> & xsecquqqup(maxtype,maxtype,maxpoints)
c
double precision tabsig
c

if (tpchin.eq.tpchout) then
if (tpchin.eq.1) then
call tsiggg22()
tabsig=xsecgg22(tpcrindex)
endif
if (tpchin.eq.2) then
call tsigqg22()
tabsig=xsecqg22(abs(tin1+tin2),tpcrindex)
endif
if (tpchin.eq.3) then
call tsigqq22()
tabsig=xsecqq22(abs(tinl),tpcrindex)
endif
if (tpchin.eq.4) then
call tsigqu22()
tabsig=xsecqu22(abs(tinl),tpcrindex)
endif
if (tpchin.eq.5) then
call tsigqqp22()

138

tabsig=xsecqqp22(abs(tinl),abs(tin2),tpcrindex)

endif

else
if

(tpchin.eq.1) then

call tsigggqu()

if (abs(tout1).gt.1) then
tabsig=xsecggqu(abs(tout1) ,tpcrindex)

-xsecggqu(abs(tout1)-1,tpcrindex)

else
tabsig=xsecggqu(1,tpcrindex)-xsecgg22(tpcrindex)

endif

endif

if

(tpchin.eq.4) then
if (tpchout.eq.1) then
call tsigqugg()
tabsig=xsecqugg(abs(tin1),tpcrindex)
-xsecqu22(abs(tin1),tpcrindex)
endif
if (tpchout.eq.6) then
call tsigquqqup()
if (abs(tout1).gt.1) then
tabsig=xsecquqqup(abs(tin1),abs(tout1) ,tpcrindex)
-xsecquqqup(abs(tin1),abs(tout1)-1,tpcrindex)
else
tabsig=xsecquqqup(abs(tinl),1,tpcrindex)
-xsecqugg(abs(tin1),tpcrindex)
endif
endif

endif

endif

if (tpsigma.ne.0.) then

if

(abs(tpsigma-tabsig)/tpsigma.gt..1) then
write(1,*)

& ’Table entry for cross section off by more than 10 percent’

write(1,*) ’tpsigma:’,tpsigma,’ tabsig:’,tabsig
write(1,*) ’tpsq:’,tpsq,’ tpcrindex ...:’,
eps+(tpcrindex-1.)**4*sqstep
write(1,*) ’tpchin/tpchout:’,tpchin,tpchout
write(1,*) ’tin1,2’,tin1,tin2,’ tout1,2’,tout1,tout2

stop

& ’Table entry for cross section off by more than 10 percent’
endif

else
if

(tabsig.gt.1.e-4) then

write(1,*) ’Table entry for zero cross section .gt. 1.e-4’

write(1,*) ’tpsigma:’,tpsigma,’ tabsig:’,tabsig

write(1,*) ’tpsq:’,tpsq,’ tpcrindex ...:’,
eps+(tpcrindex-1.)**4*sqstep

write(1,*) ’tpchin/tpchout:’,tpchin,tpchout

write(1,*) ’tin1,2’,tin1,tin2,’ tout1,2’,tout1,tout2

stop ’Table entry for zero cross section .gt. 1.e-4’

endif

endif

139

return
end

B.9 TWOTOTWO

Gerd Kortemeyer
Last revision of ﬁle: 09/ 02 / 94

B.9.1 Subroutine twototwo

This subroutine completes the scattering of two particles into an output channel with also two
particles. It is called by collide.

c Causal index, angle index, D0 variable
integer ic,ia,i
c Factors in transformation formula
double precision scalbetax,a
c Middle point
double precision midp(3)
c Velocities, radius
double precision velol,velo2,radius
c Old momenta and energies, new momentum squared, scalar product
double precision msl,ms2,me1,me2,m1(3),m2(3),nmsl,nms2,scp
c Angle
double precision labang

Now the cosine of the scattering angle 0 tpcos is calculated by the appropriate routines. Setting
tpcos to ’1’ is only done for safety reasons if the monte carlos are commented out for test purposes.

c Get tpcos=cos(theta) for scattering

c ___________________________________
c Default: No scattering
tpcos=1.
c
c Do the Monte Carlo
c
call montecarlo()
c

A small check if the cosine is really between -1 and 1 is performed:

c Check cosine
if ((tpcos.ge.1).or.(tpcos.le.-1)) then
write(1,*) ’Error in Monte Carlo, rien ne va plus.’
write(1,*) ’tpcos:’,tpcos
write(1,*) ’tpchin,tpchout:’,tpchin,tpchout
stop ’Error in Monte Carlo, rien ne va plus.’

140

The angular distribution statistics are set:

ia=int(tpcos*3.)+4
if (tpcos.1t.0) ia=ia-1
anglestat(tpchout,ia)=ang1estat(tpchout,ia)+1

Now the parton types can be set according to tout1 ,2:

 

c
c Set types for output
C-u— ——————
tq(tprn,tpp1)=tout1
tq(tprn,tpp2)=tout2
c

The scattering itself is then done by subroutine newmomenta:

c Scatter according to tpcos

call newmomenta()

The subroutine newmomenta places the new momentum for the ﬁrst particle into the vector tpscl

Some statistics for the scattering angle are set:

c

c Statistics

c
if (tpcos.1t.mincos) mincos=tpcos
if (tpcos.gt.maxcos) maxcos=tpcos

The backtransformations have to be done simply with the vector —[5’ instead of 6. However, ﬁrst a
few values are stored for statistics and checks:

c
c Backtransformation

c Remember values for angle statistics and checks
c
mel=eq(tprn,tpp1)
me2=eq(tprn,tpp2)
msl=0.
ms2=0.
do 100 i=1,3
m1(i)=pq(tprn,tpp1,i)
msl=msl+m1(i)*m1(i)
m2(i)=pq(tprn,tpp2,i)

141

ms2=ms2+m2(i)*m2(i)
100 continue
c
scalbetax=-tpbeta(1)*tpsc1(1)
& -tpbeta(2)*tpsc1(2)
& -tpbeta(3)*tpsc1(3)

c Momentum 1
a=-tpgambet*scalbetax+tpgamma*tpsceq1

eq(tprn,tpp1) =tpgamma*(tpsceql-scalbetax)
pq(tprn,tpp1,1)=tpsc1(1)+a*tpbeta(1)
pq(tprn,tpp1,2)=tpsc1(2)+a*tpbeta(2)
pq(tprn,tpp1,3)=tpsc1(3)+a*tpbeta(3)

c Momentum 2
a=tpgambet*scalbetax+tpgamma*tpsceq2

eq(tprn,tpp2) =tpgamma*(tpsceq2+scalbetax)
pq(tprn,tpp2,1)=-tpsc1(1)+a*tpbeta(1)
pq(tprn,tpp2,2)=-tpsc1(2)+a*tpbeta(2)
pq(tprn,tpp2,3)=-tpsc1(3)+a*tpbeta(3)

The new lab momentum squared is calculated for various purposes:

c Calculate new momenta squared
nmsi=0.
nms2=0.
do 300 i=1,3
nmsl=nmsl+pq(tprn,tpp1,i)*pq(tprn,tpp1,i)
nms2=nms2+pq(tprn,tpp2,i)*pq(tprn,tpp2,i)
300 continue

The transformation is checked now by again calculating the total energy:

c Check transformation
tpeql = eq(tprn,tpp1) +eq(tprn,tpp2)
tpbeta(1)=(pq(tprn,tpp1,1)+pq(tprn,tpp2,1))/tpeql
tpbeta(2)=(pq(tprn,tpp1,2)+pq(tprn,tpp2,2))/tpeql
tpbeta(3)=(pq(tprn,tpp1,3)+pq(tprn,tpp2,3))/tpeql
tpbq=tpbeta(1)**2+tpbeta(2)**2+tpbeta(3)**2

if (abs(tpsq-tpeql*tpeq1*(1.-tpbq))/tpsq.gt.0.001) then
write(1,*) ’ERROR in Backtransformation or Scatteringz’
write(1,*) tpsq,tpeql*tpeq1*(1.-tpbq)
write(1,*) eq(tprn,tpp1),pq(tprn,tpp1,3)
write(1,*) eq(tprn,tpp2),pq(tprn,tpp2,3)
stop ’ERROR in Backtransformation or Scattering’

endif

C

Also, it is checked if the parton had gotten to large an ’imaginary mass’ during the transformation.
’Imaginary masses’ occur due to numeric precision problems if the parton had a very little mass
before the transformation which during the transformation calculations turned into a very little
imaginary mass:

142

if (eq(tprn,tpp1)*eq(tprn,tpp1)-nmsl.lt.-.05) then
write(1,*) ’Particle has non-negligible imaginary mass.’
write(1,*) ’Old/new energy:’,me1,eq(tprn,tpp1)
write(1,*) ’Old/new momentumz’
write(1,*) m1(1),m1(2),m1(3)
write(1,*) pq(tprn,tpp1,1),pq(tprn,tpp1,2),pq(tprn,tpp1,3)
write(1,*) ’Old real mass :’,sqrt(me1*me1-msl)
write(1,*) ’Imaginary massz’,

& sqrt(-eq(tprn,tpp1)*eq(tprn,tpp1)+nmsl)

write(1,*) ’tpsq :’,tpsq
write(1,*) ’tpcos:’,tpcos
stop ’Particle has non-negligible imaginary mass.’

endif

if (eq(tprn,tpp2)*eq(tprn,tpp2)-nms2.1t.-.05) then

write(1,*) ’Particle has non-negligible imaginary mass.’
write(1,*) ’Old/new energy:’,me2,eq(tprn,tpp2)
write(1,*) ’Old/new momentumz’
write(1,*) m2(1),m2(2),m2(3)
write(1,*) pq(tprn,tpp2,1),pq(tprn,tpp2,2),pq(tprn,tpp2,3)
write(1,*) ’Old real mass :’,sqrt(me2*me2-ms2)
write(1,*) ’Imaginary massz’,

& sqrt(-eq(tprn,tpp2)*eq(tprn,tpp2)+nms2)
write(1,*) ’tpsq :’,tpsq
write(1,*) ’tpcos:’,tpcos
stop ’Particle has non-negligible imaginary mass.’

endif

 

 

After having made sure that the correction is negligible the parton is artiﬁcially given an eps-mass:

c Artificially correct imaginary masses. Above check makes sure that the
c correction is not too big. Particle given eps-mass.
c
if (eq(tprn,tpp1)*eq(tprn,tpp1).1t.nmsl)
& eq(tprn,tpp1)=sqrt(nmsl)+eps
if (eq(tprn,tpp2)*eq(tprn,tpp2).lt.nms2)
& eq(tprn,tpp2)=sqrt(nms2)+eps

The lab angle statistics are set:

c
c Lab angle statistics
c

scp=0.

do 200 i=1,3

 

scp=scp+pq(tprn,tpp1,i)*m1(i)
200 continue
labang=scp/sqrt(nmsl*msi)
ia=int(labang*3.)+4
if (labang.1t.0) ia=ia-1
labangle(ia)=1abangle(ia)+1

143

scp=0.
do 210 i=1,3
scp=scp+pq(tprn,tpp2,i)*m2(i)

210 continue

labang=scp/sqrt(nms2*ms2)

ia=int(labang*3.)+4

if (labang.lt.0) ia=ia-1

1abangle(ia)=labang1e(ia)+1

Now that the particles are successfully scattered the ﬂags have to be set accordingly. First the total
number of scattering events is increased by one, the channel statistics are updated:

c
c Particles successfully scattered, set flags
6 -_ _
c Total number
totscat=totscat+1
c Channel table
chscat(tpchin,tpchout)=chscat(tpchin,tpchout)+1

 

The variables fq are set to ’not belonging to the nucleus of origin anymore’, the ﬂags alr are set to
’scattered’ and in the particle data last is set to the respectively other scattering partner.

c Don’t belong to certain nucleus anymore
fq(tprn,tpp1)=0
fq(tprn,tpp2)=0

c Particles scattered in this run
alr(tprn,tpp1)=1
alr(tprn,tpp2)=1

c Set last scattering partners
1ast(tprn,tpp1)=tpp1
last(tprn,tpp2)=tpp2

The middle point of the scattering is calculated in the lab frame:

c Calculate middle point of scattering
midp(1)=(rq(tprn,tpp1,1)+rq(tprn,tpp2,1))/2.
midp(2)=(rq(tprn,tpp1,2)+rq(tprn,tpp2,2))/2.
midp(3)=(rq(tprn,tpp1,3)+rq(tprn,tpp2,3))/2.

Then information about the distance of the scattering event from the beam axis and signal velocities
is calculated. The signal velocity can only be calculated if there had been a previous scattering event:

c Causal information

radius=sqrt(midp(1)**2+midp(2)**2)
c Velocity information

if (nnq(tprn,tpp1).eq.0) then
velol=sqrt( (midp(1)-lastpoint(tprn,tpp1,1))**2

& +(midp(2)-lastpoint(tprn,tpp1,2))**2

& +(midp(3)-1astpoint(tprn,tpp1,3))**2 ) /

& ( (ntime - 1asttime(tprn,tpp1))*timstp )

 

 

144

else
ve101=-1
endif

if (nnq(tprn,tpp2).eq.0) then
velo2=sqrt( (midp(l)-lastpoint(tprn,tpp2,1))**2

& +(midp(2)-lastpoint(tprn,tpp2,2))**2
& +(midp(3)-1astpoint(tprn,tpp2,3))**2 ) /
& ( (ntime - lasttime(tprn,tpp2))*timstp )
else
velo2=-1
endif

After the ﬁrst scattering event the partons are not considered part of a certain nucleon anymore.
With nnq one can check that, ’1’ means ’belonging to a certain nucleon’, ’-1’ means just created in
a scattering process, ’0’ means ’already scattered itself’.

c Don’t belong to certain nucleon anymore
nnq(tprn,tpp1)=0
nnq(tprn,tpp2)=0

Finally, the information about the last scattering point is stored to the particle data:

c Set lastpoint
lasttime(tprn,tpp1)=ntime
lasttime(tprn,tpp2)=ntime

lastpoint(tprn,tppl,1)=midp(1)
lastpoint(tprn,tpp2,1)=midp(1)
lastpoint(tprn,tppl,2)=midp(2)
lastpoint(tprn,tpp2,2)=midp(2)
lastpoint(tprn,tpp1,3)=midp(3)
lastpoint(tprn,tpp2,3)=midp(3)

The inﬁmation about the scattering event is written to the scattering data ﬁle:

c
if (scattering.ne.0) write(4,*)
& ntime*timstp,midp(3),radius,velol,velo2,sqrt(tpsq)
c
return
end
c

B.9.2 Subroutine montecarlo

This procedure calls the monte carlo routine belonging to tpchin and tpchout:

c Call right monte carlo

 

145

subroutine montecarlo()

include ’stdec.for’

C

c Channels:

6 1-ss

6 2-ss

0 3-qq

c 4 - qu

c 5-qqp

c 6-qqup

C

c TSIGGG22.FOR;1 TSIGGGQQB.FOR;1 TSIGQG22.FOR;1 TSIGQQ22.FOR;1
c TSIGQQB22.FOR;1 TSIGQQBGG.FOR;1 TSIGQQBQPQBP.FOR;1 TSIGQQP22.FOR;1
c MONTEGG.FOR;1 MONTEGGQQB.FOR;1 MONTEGG.FOR;1 MONTEQQ.FOR;14
c MONTEQQB.FOR;3 MONTEQQBGG.FOR;1 MONTEQQBP.FOR;3 MONTEQQP.FOR;11
C

if (tpchin.eq.tpchout) then
if (tpchin.eq.1) call montegg()
if (tpchin.eq.2) call montqu()
if (tpchin.eq.3) call monteqq()
if (tpchin.eq.4) call montequ()
if (tpchin.eq.5) call monteqqp()
else
if (tpchin.eq.1) call monteggqu()
if (tpchin.eq.4) then
if (tpchout.eq.1) call montequgg()
if (tpchout.eq.6) call montequp()
endif
endif

 

return
end

B.10 NEWMOMENTA

 

Gerd Kortemeyer
Last revision of the ﬁle: 10/ 19/ 93

The subroutine newmomenta is responsible to rotate the momentum vectors in the c.m.-frame ac-
cording to the given scattering angle c0319 stored in tpcos.

B.10.1 Local variables

double precision pxr(3),pxra,uperp(3)
double precision pxpxr(3),pxpxra,vperp(3),eperp(3)
double precision tpsin,tpmoa,scp,phi,siphi,cophi

146

B.10.2 Construction of a coordinate system

First the cross product of the distance vector and the momentum vector is calculated. This is turned
into a unit vector perpendicular to the momentum vector:

0

Calculate cross product of tprc and tpmol

pxr(1)=tprc(2)*tpmol(3)-tprc(3)*tpmol(2)
pxr(2)=tprc(3)*tpm01(1)-tprc(1)*tpm01(3)
pxr(3)=tprc(1)*tpmol(2)-tprc(2)*tpmol(1)

c
c Calculate length of this vector ...
c

pxra=sqrt(pxr(1)*pxr(1)+pxr(2)*pxr(2)+pxr(3)*pxr(3))
c
c ... and produce unit-vector
c

uperp(1)=pxr(1)/pxra
uperp(2)=pxr(2)/pxra
uperp(3)=pxr(3)/pxra

 

Now another unitvector is calculated which is perpendicular to both the ﬁrst unit vector and the
momentum vector. It is in the plane of the momentum vector and the distance vector, and points
into the direction of the distance vector.

c
c Calculate cross product of tpmol and pxr.

c This vector is in the (tprc,tpmol) plane and perpendicular to tpmol,
c pointing away from tprc

c
pxpxr(1)=tpmol(2)*pxr(3)-tpm01(3)*pxr(2)
pxpxr(2)=tpm01(3)*pxr(1)-tpmol(1)*pxr(3)
pxpxr(3)=tpmol(1)*pxr(2)-tpm01(2)*pxr(1)

c

c Calculate length of this vector ...

c

pxpxra=sqrt(pxpxr(l)*pxpxr(1)+pxpxr(2)*pxpxr(2)+pxpxr(3)*pxpxr(3))
c ... and produce unit-vector

vperp(1)=pxpxr(1)/pxpxra

vperp(2)=pxpxr(2)/pxpxra

vperp(3)=pxpxr(3)/pxpxra

Now with uperp and vperp a coordinate system for the rotation is provided.

B.10.3 Getting the angle (,0

The angle (,0 is chosen randomly between 0 and 27r. siphi and cophi are used as appreviations for
sincp and coscp respectively:

 

147

c
c Get phi randomly between 0 and 2pi
c

phi=2.*pi*ran1(iseed)

siphi=sin(phi)
cophi=cos(phi)

Now with uperp and vperp as coordinate system a vector eperp which is perpendicular to tpmol,

but inclosing the angle <p with vperp, is constructed:

c

c Produce unit—vector with angle phi

c
eperp(1)=siphi*uperp(1)+cophi*vperp(1)
eperp(2)=siphi*uperp(2)+cophi*vperp(2)
eperp(3)=siphi*uperp(3)+cophi*vperp(3)

Calculate the sine of the scattering angle 6:

c
c Calculate sin(theta)
c

tpsin=sqrt(1-tpcos*tpcos)

Calculate the length of the momentum vector:

c
c Calculate length of momentum vector
c

 

tpmoa=sqrt(tpmoI(1)*tpmol(1)+tpm01(2)*tpmol(2)+tpmol(3)*tpmol(3))

Now everything is prepared for the rotation of the momentum vector. The rotated vector is called

tpsc1:

C

c Rotate momentum theta degrees in eperp-direction

C

tpsc1(1)=tpmol(1)*tpcos+tpmoa*eperp(1)*tpsin
tpscl(2)=tpmol(2)*tpcos+tpmoa*eperp(2)*tpsin
tpscl(3)=tpmol(3)*tpcos+tpmoa*eperp(3)*tpsin

Finally it is checked if the angle between tpmol and tpscl is really 19 degrees by calculating the

scalar product of them:

c
c Check rotation
c

 

scp=tpsc1(1)*tpmol(1)+tpsc1(2)*tpmo1(2)+tpsc1(3)*tpmol(3)
if (abs(scp/(tpmoa*tpmoa)-tpcos).gt.0.01) then
write(1,*) ’ERROR in rotation of momenta’
write(1,*) ’MO: ’,tpm01(1),tpm01(2),tpmol(3)

148

write(1,*) ’SC: ’,tpsc1(1),tpsc1(2),tpsc1(3)
write(1,*) scp/(tpmoa*tpmoa),tpcos
stop

endif

For this elastic scattering the energies didn’t change, so the energies of the scattered particles
tpeqscl and tpeqsc2 are set to the old values:

c
c Energies didn’t change
c
tpsceq1=tpeq1
tpsceq2=tpeq2
c
return
end

B.11 DECAY

Gerd Kortemeyer
Last revision of the ﬁle: 01 / 31 / 94

This subroutine simulates the decay of particles. Besides move and scatter it is one of the main
routines of the program. See A.6 on page 83 for a description.

subroutine decay()

B.11.1 Local variables

c DO variables
integer test,i

c Difference from on-shell mass times timestep length
double precision ddx

c Function for on-shell masses
double precision mshell

B.11.2 Decaying of the particles

The subroutine checks all testruns and particles:

c All testruns, all particles
do 10 test=1,ntr
do 10 i=1,nq(test)

Particles in question are — of course — existing, they have not collided in the scatter step and don’t
belong to a certain nucleon:

c Not deleted, not already collided and not hadronized
if ((1ast(test,i).ne.-1).and.(alr(test,i).ne.1).and.
& (nnq(test,i).le.0) ) then

 

 

149

The variable ddx represents the value of Amtstep:

c Calculate difference from on-shell mass times timestep length
ddx=abs(sqrt( eq(test,i) *eq(test,i)
-pq(test,i,1)*pq(test,i,1)
-pq(test,i,2)*pq(test,i,2)
-pq(test,i,3)*pq(test,i,3) )
- mshell(tq(test,i)) ) * timstp

ﬁ'ﬁﬁﬁ'

The function mshell (tq( . . .)) provides the on—shell masses of the different parton types.

Now a random number is calculated and compared to the probability for decay, see A.6 on page 83.
If it is smaller, the subroutine delpar (see B.12.1 on page 149) is called to delete the particle:

c Delete particle?
if (ran1(iseed).lt.(ddx/(1.+ddx))) call delpar(test,i)

Now two other partons should be created, however, the theory for this is not yet worked out.

endif
10 continue
c
return
end

B.12 CREDEL

Gerd Kortemeyer
Last revision of the ﬁle: 02/25/94

The subroutines in this ﬁle are responsible for creating and deleting particles. Also, a routine to
check the validity of the particles data structure is included.

B.12.1 delpar

This subroutine is to be used when a particle should be deleted.
subroutine delpar(test,n)

The two parameters test and 11, both integer, specify the testrun and the particle number respec-
tively.

The local variables are

c DO variables
integer i,j,k

c Highest existing particle index
integer highexpar

c Found particle in nextgrid chain?
integer found

150

First it is checked if the particle to be deleted exists. This is done by means of the last-ﬂag, see
BJSonpamSI

if (last(test,n).eq.-1) then
write(1,*) ’Trying to delete non-existing particle.’
write(1,*) ’test:’,test,’ n:’,n
call checknumbers(test,n,0)
stop ’Trying to delete non-existing particle.’
endif

The procedure checknumbers is described in B.12.3 on page 154.
In the next step the particle is deleted:

last(test,n)=-1

Instead of just marking the particle ’removed’ and then during scatter always to check if the chosen
particles exist it is faster to completely remove the particle from the nextgrid pointer structure.
This is simply done by having the previous particle in the chain point to the subsequent one. The
number of times the particle was found in the chain is counted in found. The result can only be ’1’
if it was inside of the chain, or ’0’ if it was a gridroot.

In addition to that in this loop the highest occupied array index is determined:

c Skip particle by removing it from the nextgrid chain.
c Find highest array index.
c

found=0

highexpar=-1

do 10 i=1,nq(test)
if (last(test,i).ne.-1) then
if (nextgrid(test,i).eq.n) then
nextgrid(test,i)=nextgrid(test,n)
found=found+1
endif
highexpar=i
endif
10 continue

Now if the loop was executed without ﬁnding the particle it must have been a gridroot. If so, the
gridroot is set to the subsequent particle:

c If the loop was executed without finding the particle ...

c
if (found.eq.0) then

c

c . maybe it was a gridroot?

c

do 30 i=—mgridxy,mgridxy
do 30 j=-mgridxy,mgridxy
do 30 k=-mgridz,mgridz

 

151

if (gridroot(test,i,j,k).eq.n) then
gridroot(test,i,j,k)=nextgrid(test,n)

found=1
goto 20
endif

30 continue

If the particle wasn’t a gridroot either there is a serious problem. One might think that it still could
be possible that the respective array position was deleted and than on the same position another
particle could have been created in the same timestep — this particle would not be member of any
particle chain, see B.12.2 on page 152. But a freshly created particle can not be deleted in the same
timestep because for this timestep it by deﬁnition has already interacted, so:

c Still searching? Well, there is a serious problem!

c
write(1,*) ’Particle to be deleted not found in nextgrid chain’

write(1,*) ’Testrun’,test,’, particle’,n
write(1,*) ’nq(test):’,nq(test)
stop ’Particle to be deleted not found in nextgrid chain’

The next endif belongs to the found.eq.0 condition.

 

endif

The next label is the one it is jumped to when the particle was found to be a gridroot. The program
is also running over this label when the particle was found as a nextgrid:

20 continue

The next run-time check is if the particle was found more than one time in the chains:

if (found.ne.1) then
write(1,*) ’Corrupted nextgrid chain.’
write(1,*) ’Particle’,n,’ test’,test,’ found’,found,’ times’
stop ’Corrupted nextgrid chain.’

endif

If no existing particle was found below nq there might be a problem with it:

if (highexpar.lt.0) then
write(1,*) ’nq(test) found corrupted while deleting particle’
write(1,*) ’Testrun’,test,’, nq(test)’,nq(test)
stop ’nq(test) found corrupted while deleting particle’

endif

As everything seems to be going pretty well, the highest particle index for this testrun is updated
according to highexpar:

nq(test)=highexpar

Finally, the number of deleted particles is updated:

152

c Statistics
c
delnumber(test)=delnumber(test)+1

Notice that the nextgrid pointer of the deleted particle is not changed! This is important because
otherwise things would get screwed up in subroutine scatter. The nextgrid of the deleted particle
still points to a valid member of the chain. The only thing that happened to the deleted particle:
It is not pointed to anymore.

return
end

B.12.2 crepar

This subroutine creates a particle with the information from the common newpart.
subroutine crepar()
Only one local variable is deﬁned:

c DO variable
integer i

The procedure ﬁrst searches for a free space in the particle arrays:

c Find empty entry

c
do 10 i=1,mquark
if (last(newtest,i).eq.-1) goto 20
10 continue

When the loop was fully executed, no free space was found:

c ... Array is full
write(1,*) ’Could not create particle, stack capacity exceeded.’

write(1,*) ’Testrun’,newtest,’, nq(newtest)’,nq(newtest)
stop ’Could not create particle, stack capacity exceeded.’

The label 20 is jumped to as soon as a free space is found:
20 continue

The free space could be a hole in the positions 1 to nq(test) or the ﬁrst entry above nq(test). In
the latter case nq has to be updated. However, the new value can at most differ by one from the
old value:

153

if (i.gt.nq(newtest)) then
if (i-nq(newtest).ne.1) then
write(1,*)

& ’nq(newtest) found corrupted while creating particle’
write(1,*) ’Testrun’,newtest,’, nq(newtest)’,nq(newtest)
stop ’nq(newtest) found corrupted while creating particle’

endif
nq(newtest)=i
endif

The variable crenumber is updated:
crenumber(newtest)=crenumber(newtest)+1

The position now is marked ’used’ by setting the last scattering partner to ’itself’ instead of ’-1’:
last(newtest,i)=i

Now the remaining data from the common block newpart is written to the free position. The nnq
entry is set to ’-1’ meaning ’not belonging to a certain nucleon’.

c Position
rq(newtest,i,1)=nequ(1)
rq(newtest,i,2)=nequ(2)
rq(newtest,i,3)=nequ(3)

c Momentum
pq(newtest,i,1)=newpq(1)
pq(newtest,i,2)=newpq(2)
pq(newtest,i,3)=newpq(3)

c Energy
eq(newtest,i)=neweq

c Type
tq(newtest,i)=newtq

c Nucleon number
nnq(newtest,i)=-1

c Nucleus number for ’already scattered’
fq(newtest,i)=0

The updating of nextgrid will be done in the next move step.

Notice that thenextgrid pointer of the created particle is not changed! This is important because
it could very well happen that the created particle is on an array-position that was occupied by
a particle that has been deleted in the same collision. If the nextgrid pointer would have been
modiﬁed things would get screwed up in subroutine scatter. The nextgrid of the created particle
still points to a valid member of the deleted particles chain. But it is not pointed to. The new
particle must not be scattered in this timestep anyway, so. . .

return
end

154

B.12.3 checknumbers

This subroutine checks the validity of the particle data structures. It can also write information
about two particles chosen by the parameters to the protocol ﬁle - this option can be used for debug
purposes and had been very useful in the past.

subroutine checknumbers(test,n1,n2)

The parameters are the testrun number and the particle numbers of the two particles to be checked.
If the parameters are set to ’0’ no particles are checked:

c Parameters
integer test,n1,n2
c Do variables
integer i,j,k,1,m
c Count variables
integer pcount, tcount,typecount(-6:6)

The procedure itself is coded pretty straightforward:

do 10 i=1,ntr

c
c Count particles in data arrays
c
pcount=0
c
do 20 j=1,mquark
if (last(i,j).gt.0) pcount=pcount+1
20 continue
c
c Count particles in pointer chain
c
tcount=0
c
do 25 j=-6,6
typecount(j)=0
25 continue
c
do 30 j=-mgridxy,mgridxy
do 30 k=-mgridxy,mgridxy
do 30 l=-mgridz,mgridz
c
m=gridroot(i,j,k,l)
if (m.eq.-1) goto 40
c
if ((m.eq.n1).and.(i.eq.test)) then
write(1,*) ’ ------------------------------- ’
write(1,*) ’Particle ’,n1,’ is gridroot!’
endif
if ((m.eq.n2).and.(i.eq.test)) then
write(1,*) ’ ------------------------------- ’
write(1,*) ’Particle ’,n2,’ is gridroot!’
endif

155

50 tcount=tcount+1
typecount(tq(i,m))=typecount(tq(i,m))+1

if ((m.eq.nl)
write(1,*)
write(1,*)
write(1,*)

endif

if ((m.eq.n2)
write(1,*)
write(1,*)
write(1,*)

endif

.and.(i.eq.test)) then

3... ___________ ’

 

’Particle:’,n1
’Gridbox :’,j,k,l

.and.(i.eq.test)) then

i _______________________________ )
’Particle:’,n2
’Gridbox :’,j,k,l

m=nextgrid(i,m)
if (m.ne.-1) goto 50

c
40 continue
30 continue
c
c Print out information
c
write(1,*)

& ’Number of particles cre/del: ’,crenumber(i),delnumber(i)

write(1,*)

& ’Total number of particles : ’,pcount

write(1,*)
& ’Grid box chains

’,tcount

write(1,*) ’Gluons:’,typecount(0)

write(1,*) ’u, ubar:’,typecount(1),typecount(-1)
write(1,*) ’d, dbar:’,typecount(2),typecount(-2)
write(1,*) ’s, sbar:’,typecount(3),typecount(-3)
write(1,*) ’c, cbar:’,typecount(4),typecount(-4)
write(1,*) ’b, bbar:’,typecount(5),typecount(-5)
write(1,*) ’t, tbar:’,typecount(6),typecount(-6)

c
c Right number of particles?
c

if (lastnumber(i)+crenumber(i)-delnumber(i).ne.pcount) then
write(1,*) ’Total number of particles incorrect.’

write(1,*) ’pcount:

’,pcount,’, lastnumber’,lastnumber(i)

stop ’Total number of particles incorrect’

endif
c
lastnumber(i)=pcount
c
10 continue
c

return
end

 

 

 

156

B.13 PERMUTE

(SendlKorunneyer
Last revision of the ﬁle: 12/ 7 / 93

This subroutine permutes the permutation arrays.

B.13.1 Local variables

Besides the variables from stdec.for the following local variables are deﬁned:

c Main loop
integer k

c Entries to be exchanged
integer i,j

c Buffer
integer exch

B.13.2 Permutation

The permutation of the respective arrays is done simply by exchanging size of array times two entries
of the arrays.

First the array perm is permuted:

do 10 k=1,mquark
c Find entry numbers
i=1+int(ran1(iseed)*(mquark-1))
j=1+int(ran1(iseed)*(mquark-1))
if ((i.gt.mquark).or.(i.lt.1).or.(j.gt.mquark).or.(j.lt.1))
& then
write(1,*) ’Error in random numbers, i,j=’,i,j
write(1,*) ’mquark:’,mquark
stop
endif
c Exchange them
exch=perm(i)
perm(i)=perm(j)
perm(j)=exch
c
10 continue

Now permx and permy are permuted. One has to be careful not to permute the position of the
outmost boxes since they can not be used as center boxes in subroutine scatter and are taken
out of the 100p there by simply starting the 100p variable at for example -mgridxy+1. This simple
means of skipping the outmost boxes would not work if the index for mgridxy would also have been
permuted:

do 20 k=-mgridxy,mgridxy
c Find entry numbers, don’t permutate outmost boxes
i=int((1.-2.*ran1(iseed))*(mgridxy-2))
j=int((1.-2.*ranl(iseed))*(mgridxy-2))

 

157

c Exchange
exch=permx(i)
permx(i)=permx(j)
permx(j)=exch
c Find entry numbers, don’t permutate outmost boxes
i=int((1.-2.*ranl(iseed))*(mgridxy-2))
j=int((1.-2.*ran1(iseed))*(mgridxy-2))
c Exchange
exch=permy(i)
permy(i)=permy(j)
permy(j)=exch
20 continue

Since permz has another size than permx and permy it is permuted separately:

do 30 k=-mgridz,mgridz
c Find entry numbers, don’t permutate outmost boxes
i=int((1.-2.*ran1(iseed))*(mgridz-2))
j=int((1.-2.*ran1(iseed))*(mgridz-2))
c Exchange
exch=permz(i)
permz(i)=permz(j)
permz (j )=exch
30 continue

Finally deltax, deltay and deltaz are permuted:

i=1+int(3.*ran1(iseed))
j=1+int(3.*ran1(iseed))
exch=deltx(i)
deltx(i)=deltx(j)
deltx(j)=exch

i=1+int(3.*ran1(iseed))
j=1+int(3.*ran1(iseed))
exch=delty(i)
delty(i)=delty(j)
delty(j)=exch

i=1+int(3.*ran1(iseed))
j=1+int(3.*ran1(iseed))
exch=deltz(i)
deltz(i)=de1tz(j)
deltz(j)=exch

return
end

B.14 DISTRI

CTEQ Collaboration

 

 

 

158

Last revision of the ﬁle: 10/7/93
These subroutines calculate the parton distribution function ctq2pd as they will be used by init.

C

0

00000 0000 O

0

By:

To

FUNCTION th2Pd (Iset, Iparton, x, Q, Irt)
Version 2 CTEQ distribution function in a parametrized form.
J. Botts, H.L. Lai, J.G. Morfin, J.F. Owens, J. Qiu, W.K. Tung & H.Weerts

avoid the proliferation of parton distribution functions, we recommend that

these distributions should replace Version 1 CTEQ distributions for all general
usage.

The

(No data tables are needed.)
returned function value is the PROBABILITY density for a given FLAVOR.

A companion function (next module), which this one depends on,
th2df (Iset, Iparton, X, Q, Irt)
gives the VALENCE and SEA MOMENTUM FRACTION distributions.

Also included are two functions Prth2 & Wlamd2 which return relevant
QCD parameters associated with these parton distributions sets.(See below)

A parallel (independent) program th2Pds (not included in this file)
in Subroutine form is also available.

It returns ALL the parton flavors at once in an array form.

See details in that separate file if you are interested.

Since this is an initial distribution, and there may be updates, it is
useful for the authors to maintain a record of the distribution list.
Please do not freely distribute this program package; instead, refer any
interested colleagues to direct their request for a copy to:
Bottthades.ifh.de or Lachteq11.pa.msu.edu

If you have detailed questions concerning these distributions, direct inquires
to Botts, Lai (see above) or Wu-Ki Tung (Tunngsupa.pa.msu.edu).

This function returns the CTEQ parton distributions f‘Iset_Iprtn/proton

where Iset (= 1, 2, ..., 5) is the set label;

Name convention for CTEQ distributions: CTEQnSx where
n : version number (currently n = 1)
S : factorization scheme label: = [M D L] for [MS—bar DIS LO]

resp.
x : special characteristics, if any

(e.g. S(F) for singular (flat) gluon, L for "LEP lambda value")

Iprtn is the parton label (6, 5, 4, 3, 2, 1, O, -1, ...... , -6)

for (t, b, c, s, d, u, g, u_bar, ..., t bar)

 

 

 

159

C X, Q are the usual x, Q; Irt is a return error code (not implemented yet).

C ——> Iset = 1, 2, 3, 4, 5 correspond to the following CTEQ global fits:
C cteq2M, cteq2MS, cteq2MF, cteq2ML, cteq2L respectively.

C --> QCD parameters for parton distribution set Iset can be obtained inside

C the user’s program by:

C Dum = Prctq2

C > (Iset, Iord, Ischeme, MxFlv,

C > Alam4, A1am5, A1am6, Amas4, Amas5, Amas6,

C > Xmin, Qini, Qmax, ExpNor)

C where all but the first argument are output parameters.

C They should be self-explanatory —- see details under ENTRY Prctq2.

C Since the QCD Lambda value for the various sets are needed more often than

C the other parameters in most applications, a special function

0

Wlamd2 (Iset, Iorder, Neff) is provided
C which returns the lambda value for Neff = 4,5,6 effective flavors as well as
C the order these values pertain to.

C The range of (x, Q) used in this round of global analysis is, approxi—

C mately, 0.01 < x < 0.75 ; and 4 GeV‘2 < 0‘2 < 400 GeV“2 for fixed target
C experiments and 0.0001 < x < 0.01 from preliminary data of HERA.

C The range of (x, Q) used in the reparametrization of the QCD evolved

C parton distributions is 10E—5 < x < 1 ; 1.6 GeV < Q < 1 TeV. The

C functional form of this parametrization is:

C A0 * x‘Al * (1-x)‘A2 * (1 + A3 * x“A4) * [log(1+1/x)]‘A5

C with the A’coefficients being smooth functions of O. For heavy quarks,

C a threshold factor is applied to A0 which simulates the proper Q-dependence
C of the QCD evolution in that region.

C Since this function is positive definite and smooth, it provides sensible

C extrapolations of the parton distributions if they are called beyond

C the original range in an application. There is no artificial boundaries

C or sharp cutoff’s.

 

 

 

 

 

Appendix C

Index

alr, 92, 119, 121, 143, 148
already, 92
anglestat, 92, 110, 140

base, 120

basesq, 120

Boltzmann equation, 20, 24
Boltzmann-Uehling-Uhlenbeck (BUU), 2
Bonds, 6

BUU, 10, 22, 30

buu scatteringdat, 144

c.m.—system, 82

Cascade codes, 2, 38
causal, 94

Causality information, 143
Causality violations, 38
causalpl, 93, 96

causcatter, 122

CBA, 23

Channels, 91, 132, 135, 137
charge, 93, 104
checknumbers, 98, 154
chscat, 92, 110, 143
Closest approach criterion, 41
Clusters, 13

cmtimstp, 129

collide, 121, 123

collisions, 82

Consistent Boltzmann Algorithm, 23
constant, 97
control, 93
CPU-time, 93
cralphas, 94
create, 90, 101
credel, 149
creﬂag, 90, 101
crelabpn, 90, 101
cremodir, 90, 101
crenn, 90, 101
crenp, 90, 101
crenumber, 92, 153, 154
crepar, 152
crex0, 90, 101
crezO, 90, 101
crk, 94

crm, 94

cm, 94

ctq2pd, 158
Cutoff, 94

decay, 97, 148
delnumber, 92, 151, 154
delpar, 149

deltax, 120, 157

deltay, 120, 157

deltaz, 120, 157

deltx, 102, 105

delty, 102, 105

deltz, 102, 105

Direct Simulation Monte Carlo, 22
Displacement, 25

160

 

   

 

distri, 157
DSMC, 22

elabpnpro, 88
elabpntar, 88
elabpro, 89

elabtar, 89

endpart, 117

Enskog theory, 23
E08, 20

eps, 86, 106, 126, 127, 142
eq, 87, 102, 111, 131
Equation of state, 20
Evaporation, 1, 13
Excitation energy, 13

fdue, 88

feq, 88

Fission, 1, 13

Flow, 20

qu, 88

fq, 87, 90, 101, 102, 111, 112, 143
Fragmentation, 6

Frame dependencies, 47

Full-ensemble testparticle algorithm, 24, 30

Gluon—gluon cross section, 40
graphics output, 93

grid, 81

grid boundaries, 81, 103
gridroot, 115, 121, 150, 154
gridspxy, 88, 103

gridspz, 89, 103, 115

gridz, 89, 103

161

Hard processes, 38

Hard-sphere Boltzmann equation, 23
Hard-sphere radius, 25

hbc,85

IMF, 1, 7

impact, 110

init, 95, 100
initcharge, 104
initcheck, 105
initcross, 106, 136
initﬂag, 104
initmeth, 90, 100, 103
initmethl, 104, 110
initparms, 90
initperm, 105

initqsq, 90

initrnuc, 90

initxmin, 90
Intermediate mass fragment, 1, 7
iseed, 88, 95

Isospin, 6

labangle, 92, 142
Landau—Vlasov, 2
last, 87, 102, 104, 121, 148, 153
lastnumber, 92
lastpoint, 88, 123
1asttime, 123

Lattice, 6

LHC, 4

Lightcone, 46

Lorentz invariance, 45
lowcut, 89

makenucleus, 90, 101, 103
manuc, 85
maxchannel, 86

 

 

 

 

 

maxcos, 92, 109, 140
maxdq, 124

maximum cross section, 81
maxpoints, 86, 106, 126
maxsig, 88, 106, 125, 136
maxsigma, 92, 102, 109, 126
maxsq, 92, 109, 126
maxtype, 86

mean ﬁeld, 79, 82

mincos, 92, 109, 140
mindq, 124

minsigma, 92, 102, 109, 126
minsq, 92, 109, 126
Molecular Dynamics, 22
montecarlo, 144

move, 97, 114

move step, 81

mquark, 86

mtest, 86

multiple scattering, 81

nc, 95

Neutron skin, 18

neweq, 92

newfq, 153

newmomenta, 140, 145
newpart, 92, 152

newpq, 92, 153

newpscat, 92

nequ, 92, 153

newtest, 92

newtq, 92, 153

nextgrid, 88, 122, 150, 153, 154
nmts, 88

nnpro, 87

nnq, 87, 102, 111, 112, 144, 148, 153
nntar, 87

nppro, 87

nptar, 87

nq, 87, 101, 102, 111, 116, 151, 152
nsta, 92, 109

nstaO, 124

nstal, 125

ntime, 88, 103

162

ntr, 87
Nuclear ﬂow, 20

One particle distribution function, 24
others, 88

outb, 115

outchannel, 132

Parallel-ensemble testparticle algorithm, 30
parameter ﬁle, 90, 95, 103
parmin, 87

Partial cross section, 40, 91, 93, 106, 135,
137

particles, 86

parton distribution, 39, 157

Pauli principle, 27

Percolation model, 6

perm, 93, 102, 105, 116, 156
permute, 156

permx, 102, 105, 119, 156

permy, 102, 105, 119, 156

permz, 102, 105, 119, 157

pi, 86

ploto, 93, 96

plotout, 96

pointer, 82

pointscat, 92, 102, 104

pq, 87, 102, 111

preferential directions, 81
projectile, 101

Proper time, 75

protocol ﬁle, 92, 95, 103, 128, 154
prox0, 88

pron, 88

qq, 102
Quark-gluon plasma, 4, 38

 

 

 

r0, 85

rancor, 109

Random walk, 44
Reduced volume, 27
restat, 98, 109
Retarded interaction, 72
RHIC, 4, 38

rq, 87, 102, 111

scalbetax, 127

scatter, 97, 118, 123, 156
Scattering, 81

scattering, 93

Scattering data ﬁle, 144
Screening, 28

Second virial coefﬁcient, 29
Signal velocity, 143

Sites, 6

Soft processes, 38
Spacelike distance, 46
Spectator region, 10
sqstep, 86, 106, 126

stat, 92

stcount, 94, 95

stdec, 85

steps, 94, 95

stymax, 94, 95

subperm, 116
Superluminous signals, 40

t0, 93, 94

target, 101

test run, 79

timefault, 92, 102, 109, 129
Timestep length, 63
timstp, 88, 131

tinl, 91, 106, 126, 134, 136
tin2, 91, 106, 126, 134, 136

163

Total cross section, 91, 94, 108, 135
totalcross, 91, 133, 135
totscat, 143

tout1, 91, 106, 134, 140
tout2, 91, 134, 140
tpbeta, 91, 125

tpbq, 91, 125

tpchin, 91, 126, 133, 135, 144
tpchout, 91, 132, 133, 144
tpcos, 91, 139, 145
tpcrindex, 91, 126

tpdq, 91, 124

tpdqc, 91, 129

tpeq1,90

tpeq2, 90

tpeql, 91, 125

tpeqscl, 148

tpeqch, 148

tpgambet, 91, 127
tpgamma, 91, 127

tpmol, 91, 131

tpmo2, 131

tppl, 90

tpp2,90

tprc, 91, 129, 131

tprl, 90, 124

tprn, 90, 119

tpscl, 91, 140, 147
tpsceql, 91

tpsceq2, 91

tpsigma, 91, 126, 133, 136
tpsq, 90, 126

tptd, 91, 129

tq, 87, 102, 111, 140
Transformation equations, 48
Transport models, 3
twopart, 90, 119
twototwo, 132, 139

uperp, 146

 

 

164

vc, 94, 95
velocitypl, 93, 96
veymax, 94, 95
vgrid, 94, 95

Virial coefﬁcient, 29
vperp, 146

Wee partons, 74
whenout, 88, 97

xsec..., 102, 106, 132

xsecentry, 91, 136

xsecgg22, 93, 108, 133, 137
xsecggqu, 93, 108, 133, 136, 137
xsecqg22, 93, 133, 136, 137
xsecqq22, 93, 133, 136, 137
xsecqu22, 93, 94, 108, 133, 134, 137
xsecqugg, 93, 94, 108, 133, 134, 137
xsecquqqup, 93, 94, 108, 133, 134, 136,
137

xsecqqp22, 93, 133, 136, 137

 

Bibliography

[Aic91]

[Aic96]

[Ale95]

[Alm95]
[Bad89]
[Bau84]

[Bau86]

[Ber88]
[Bot93]

[Bro96]
[Cam85]

J. Aichelin, Phys. Rep. 202, 233 (1991); H. Feldmeier, Nucl. Phys. A515,
147 (1990); A. Ono et al., Prog. of Theoret. Phys. 87, 1185 (1992); H.
Sorge, H. St6cker, and W. Greiner, Nucl. Phys. A498, 567C (1989); H.
Sorge, H. Stﬁcker, and W. Greiner, Ann. Phys. 192, 266 (1989).

J. Aichelin, Phys. Rev. C33, 537 (1986); H. Stécker, and W. Greiner, Phys.
Rep. 137. 277 (1986); U. Mosel, Annu. Rev. Nucl. Part. Sci. 41, 29 (1991);
W. Bauer et al., Annu. Rev. Nucl. Part. Sci. 42, 77 (1992); P. Schuck et
al., Prog. Part. Nucl. Phys. 22, 181 (1989); Y. Pang, T. Schlagel, and S.
H. Kahana, Nucl. Phys. A544, 435C (1992); D. E. Kahana, D. Keane, Y.
Pang, T. Schlagel, and S. Wang, Phys. Rev. Lett. 74, 4404 (1995).

F. J. Alexander, A. L. Garcia, and B. J. Alder, Phys. Rev. Lett. 74, 5212
(1995)

T. Alm et al., Nucl. Phys. A587, 815 (1995)
H. Badovsky, Eur. J. Mech., B/Fluids 8, 41 (1989)

W. Bauer, D. R. Dean, U. Mosel, and U. Post, in Procedz'ngs of the 7th
High Energy Heavy Ion Study (Report GSI-85-10, 1984), 701; W. Bauer,
D. R. Dean, U. Mosel, an U. Post, Phys. Lett 150B, 53 (1985); W. Bauer,
U. Post, D. R. Dean, and U. Mosel, Nucl. Phys. A452, 699 (1986); W.
Bauer, Phys. Rev. C38, 1297 (1988); W. Bauer and A. Botvina, Phys.
Rev. C52, R1760 (1996).

W. Bauer, G. F. Bertsch, W. Cassing, and U. Mosel, Phys. Rev. C34,
2127 (1986); W. Bauer, Nucl. Phys. A471, 604 (1987); W. Bauer, G. F.
Bertsch, and H. Schulz, Phys. Rev. Lett. 69, 1888 (1992).

G. F. Bertsch and S. Das Gupta, Phys. Rep. 160, No. 4, 189 (1988)

Version 2 CTEQ distribution function in a parametrized form, J. Botts,
H. L. Lai, J. G. Morﬁn, J. F. Owens, J. Qiu, W.-K. Tung and H. Weerts,
CTEQ Collaboration.

Hartree-Fock Code DENS, B. A. Brown, priv. communication.

X. Campi and J. Debois, in Proceedings of the 23rd Bormeo Conference
(Ricerca Scientiﬁca et Educatione Permanente, Milano, 1985), 497; T. S.
Biro, J. Knoll, and J. Richert, Nucl. Phys. A459, 692 (1986); X. Campi,

165

 

 

 

[Com77]
[Con79]
[Cug81]

[Cze86]

[Dan91]
[Dan95]
[Dan96]
[Da593]

[Eic84]
[Fie89]

[Gei92.1]
[Gei92.2]
[Gei93]
[Gei94.1]
[Gei94.2]
[Gus88]

[K3393]
[Kod84]

[Kod94]

166

J. Phys. A19, L917 (1986); J. Nemeth, M. Barranco, J. Debois, and C.
Ngf), Z. Phys. A325, 347 (1986); J. Debois, Nucl. Phys. A466, 724 (1987);
C. Cerruti, J. Debois, R. Boisgard, C. Ng6, J. Natowitz, and J. Nemeth,
Nucl. Phys. A476, 74 (1988); S. Das Gupta, C. Gale, and K. Haglin, Phys.
Lett. 302B, 372 (1993).

B. L. Combridge, J. Kripfganz and J. Ranft, Phys. Lett. 70B, No. 2, 234
(1977)

A. Coniglio, H. E. Stanley, and W. Klein, Phys. Rev. Lett. 42, 518 (1979);
D. W. Hermann and D. Stauffer, Z. Phys. B44, 339 (1981).

J. Cugnon, T. Mizutani, J. Vandermeulen, Nucl. Phys. A352, 505 (1981)

L. P. Czernai, and J. I. Kapusta, Phys. Rep. 131, 223 (1986); S. Das Gupta,
and G. D. Westfall, Physics Today 46(5), 34 (1993); R. Stock, Phys. Rep.
135, 259 (1986); H. Gutbrod et al, Rep. Prog. Phys. 52, 1267 (1989).

P. Danielewicz and G. F. Bertsch, Nucl. Phys. A533, 712, appendix (1991)
P. Danielewicz, Phys. Rev. C51, 716 (1995)
P. Danielewicz, and S. Pratt, Phys. Rev. CS3, 249 (1996)

S. Das Gupta and G. D. Westfall, G. D., 1993, Physics Today 46(5), 34
(1993)

E. Eichten, I. Hinchliffe, K. Lane and C. Quigg, Rev. Mod. Phys., Vol. 56,
No. 4, 579 (1984)

R. D. Field, Frontiers in Physics: Applications of Perturbative QCD,
Addison—Wesley 1989

K. Geiger and B. Miiller, Nucl. Phys. A544, 467C (1992)
K. Geiger and B. Miiller, Nucl. Phys. B369, 600 (1992)
K. Geiger, Phys. Rev. Lett. 71, No. 19, 3075 (1993)

K. Geiger, Phys. Rev. C49, 3234 (1994)

K. Geiger, private communication

H. A. Gustafson et (11. (Plastic Ball Collaboration), Mod. Phys. Lett. A 3,
1323 (1988)

U. Kalmbach, T. Vetter, T. S. Bird and U. Mosel, Nucl. Phys. A563, 584
(1993)

T. Kodama, S. B. Duarte, K. C. Chung, R. Donangelo and R. A. M. S.
Nazareth, Phys. Rev. C29, 2146 (1984)

T. Kodama, private communication

 

[Kor93]

u<un93
[Lan92]

[Lan93]
[Mcl82]

[Mor97]
[Mur93]

[Par95]
[Pha92]

[Res77]

[80189. 1]
[Sor89.2]
[Sta79]

[Tsa93]
[Wel89]

[VWssq

167
G. Kortemeyer, J. Murray, S. Pratt, K. Haglin and W. Bauer,
NSCL/ Cyclotron Laboratory Annual Report 1993, 63, 65
G. J. Kunde et al., Phys. Rev. Lett 77 (1996), 2897.

A. Lang, H. Badovsky, W. Cassing, U. Mosel, H.—G. Reusch, and K. Weber,
J. Comp. Phys. 106, 391 (1993)

A. Lang, H. Badovsky, W. Cassing, U. Mosel, Hans-Georg Reusch and
Klaus Weber, J. Comp. Phys. 106, 391 (1993)

L. McLerran, Quark Matter and Heavy Ion Collisions; Bielefeld Workshop
1982, World Scientiﬁc, 63

K. Morawetz et al., Universitiit Rostock, in preparation.

J. Murray, G. Kortemeyer, S. Pratt, K. Haglin and W. Bauer,
NSCL/ Cyclotron Laboratory Annual Report 1993, 61

M. D. Partlan et al. (EOS Collaboration), Phys. Rev. Lett. 75, 2100 (1995)

L. Phair et al., Phys. Lett. 285B, 10 (1992); L. Phair, W. Bauer and C.
K. Gelbke, Phys. Lett. 314B, 271 (1993).

P. M. Résibois, and M. De Leener, Classical Kinetic Theory of Fluids, p.
156, John Wiley & Sons, New York, 1977

H. Sorge, H. Stécker and W. Greiner, Nucl. Phys. A498, 5670 (1989)
H. Sorge, H. Stocker and W. Greiner, Ann. Phys. 192, 266 (1989)

D. Stauffer, Phys. Rep. 540, 1 (1979); J. W. Essam, Rep. Prog. Phys. 43,
883(1980)

M. B. Tsang et al., Phys. Rev. Lett 71, 1502 (1993)

G. Welke, R. Malﬂied, C. Grégoire, M. Prakash and E. Suraud, Phys. Rev.
C40, 2611 (1989)

G. D. Westfall et al., Phys. Rev. Lett. 71, 1986 (1993); D. Kakow, G.
Welke, and W. Bauer, Phys. Rev. C48, 1982 (1993).

 

 

 

 

lllllllllllllllllllllllllllllllllllllllllllllll