text:
# QCD in Nuclear CollisionsTalk given at “Quark Matter 2002”, Nantes, July 17-24. This work is supported in part by the Department of Energy, Grant DE-FG02-94ER40819

## 1 Introduction

This lecture is devoted to a general semiquantitative description of the early stages of a heavy ion collision. This is a rapidly evolving subject and it is not yet clear how good a description can be given in terms of QCD.

In the earliest stages of the collision gluon occupation numbers are estimated to be moderately large suggesting a classical gauge theory description. This fits in nicely with the saturation picture (color glass condensate) of a high-energy heavy ion wavefunction. The parameters necessary for such a description seem compatible with an analogous description of small-x and moderate $`Q^2`$ HERA data.

Various discussions of RHIC data in terms of gluon saturation, which dicussions have much in common, are reviewed. The general picture appears to be coherent although perhaps not yet compelling. The reader should be warned that many of the numbers quoted in the text are theoretical estimates. There are large uncertainties present, neverthelelss, I believe the estimates to be reasonable.

## 2 General features of the early stages of high energy heavy ion collisions

### 2.1 Gluon number densities at the initial time

Consider a central collision at RHIC, say at $`\sqrt{s}=130GeV.`$ The measured transverse energy is $`\frac{dE_T}{dy}500GeV`$ at central rapidities. Since some of the initial transverse energy goes into longitudinal flow one certainly has

$$\frac{dE_T^{in}}{dy}>500GeV$$

for the initial transverse energy distribution, which we assume is predominately in gluons. Thus one can write

$$\frac{dE_T^{in}}{dy}=\left[\frac{dN_g^{in}}{dyd^3b}\mathrm{\Delta }b_z\right]A_gK500GeV$$
(1)
with $`b_{},`$ and $`b_z`$ the coordinates of a point in the overlaping nuclei and where $`\mathrm{\Delta }b_z`$ is the longitudinal width of the volume occupied by gluons at their time of production, $`A_g`$ is the area, transverse to the colliding beams, occupied by the gluons and $`K`$ is a typical gluon transverse momentum. Taking $`\mathrm{\Delta }b_z=1/K`$ and $`A_g150fm^2`$ gives

$$\frac{dN_g^{in}}{dyd^3b}17/fm^3$$
(2)
and with $`<K>2/3`$ (see later)

$$\frac{dN_g^{in}}{dyd^2b}5/fm^2$$
(3)
where (2) and (3), as well as $`<K>2/3GeV`$ represent average values, averaged over all transverse coordinates. At the center of the nucleus, $`b_{}=0,`$ the more appropriate value of $`K`$ is about $`3/2`$ the average value or $`KQ_s(b_{}=0)1GeV,`$ with $`Q_s`$ the gluon saturation momentum, so that

$$\frac{dN_g^{in}}{dyd^3b}|_{b_{}=0}25/fm^3,\frac{dN_g^{in}}{dyd^2b_{}}|_{b_{}=0}7.5/fm^2.$$
(4)
Equations (2-4) give a reasonable guess as to the initial gluon densities averaged over the production region (Equations (2,3)) and at the center of the production region (Equation (4)).

### 2.2 Occupation numbers

The initial gluon number densities discussed just above can be converted to gluon occupation numbers (phase space densities) by

$$f_g^{in}=\frac{(2\pi )^3}{2(N_c^21)}\frac{dN_g^{in}}{d^3pd^3b}\frac{(2\pi )^3}{2(N_c^21)}\frac{dN_g^{in}}{dyd^2b_{}d^2p_{}}$$
(5)
where we use $`\mathrm{\Delta }y=\mathrm{\Delta }p_z/p_z\mathrm{\Delta }b_z\mathrm{\Delta }p_z`$ and where the 2 in the denominator refers to the number of spin states available to gluons. We further make the approximation

$$\frac{dN_g^{in}}{dyd^2b_{}d^2p_{}}\frac{1}{\pi Q_s^2}\frac{dN_g^{in}}{dyd^2b_{}}.$$
(6)
Using (4) and (6) in (5) gives

$$f_g^{in}1.5$$
(7)
at $`b_{}=0.`$ Thus, it would seem that gluon occupation numbers are moderately large at initial times thus fitting well the picture of strong gluon field dominance at early times in heavy ion collisions.

### 2.3 The nuclear wavefunction

It is convenient to view a head-on (central) ion-ion collision in a frame where one of the ions is at rest and the other carries a large longitudinal momentum. We view the scattering in a light-cone gauge which furnishes the smoothiest connection between gluon quanta present in the wavefunction and the gluons produced at the time of the collision. Roughly, produced gluons at mid-rapidity are already in the wavefunction of the fast ion and are merely “freed” by the collision. The picture is illustrated in Fig.1. The initial wavefunction has a large occupation number for gluons, $`f_g1/\alpha ,`$ which results in a high density of produced gluons at the time of production. The fact that $`f_g`$ is large in the high-energy ion’s wavefunction is one of the features that led McLerran and collaborators to call it a color glass condensate.

### 2.4 The McLerran-Venugopalan model

The McLerran-Venugopalan model is a simple but maybe not unrealistic model of the small-x gluonic part of the wavefunction of a heavy ion at RHIC energies. The Feynman graphs contributing to the wavefunction are illustrated in Fig.2. The gluons in the wavefunction correpond to a quasiclassical picture because only tree graphs are involved even when one squares the wavefunction to get the gluon number density. The gluon number density takes a very simple form

$$\frac{dN_g}{dyd^2bd^2k_{}}\frac{dxG_A}{d^2bd^2k_{}}=\frac{N_c^21}{4\pi ^4\alpha N_c}\frac{d^2x_{}}{x_{}^2}e^{ik_{}x_{}}(1e^{x_{}^2Q_s^2/4})$$
(8)
where the saturation momentun, $`Q_s,`$ is given by

$$Q_s^2(b,x)=\frac{4\pi ^2\alpha N_c}{N_c^21}2\sqrt{R^2b^2}\rho xG_{nucleon}$$
(9)
as a function of the impact parameter in the nucleus and of the Bjorken variable, $`x.`$ Calling

$$\frac{dxG_A}{d^2b}=2\sqrt{R^2b^2}\rho xG_{nucleon}$$
(10)
one can write the saturation momentum as

$$Q_s^2=\frac{4\pi \alpha N_c}{N_c^21}\frac{dxG_A}{d^2b}.$$
(11)
It is straightforward to get the gluon occupation number in the wavefunction from (8), using a formula analogous to (5), and one finds

$$f_g=\frac{1}{\alpha N_c}\frac{d^2x_{}}{\pi x_{}^2}e^{ik_{}x_{}}(1e^{x_{}^2Q_s^2/4})\frac{1}{\alpha N_c}_1^{\mathrm{}}\frac{dt}{t}e^{tk_{}^2/Q_s^2}$$
(12)
where the approximate evaluation, given by the t-integral, is valid when $`k_{}^2/Q_s^2`$ is not too large. This form of the integral neglects a weak $`x_{}^2`$ singularity which is present in $`Q_s^2`$ via the scale of $`xG`$ in (9).

From (12) it is apparent that $`f_g1/\alpha `$ when $`k_{}^2/Q_s`$ is not large. Occupation numbers of size $`1/\alpha `$ are characteristic of classical field theory solutions, here using the valence quarks of the nucleons of the nucleus as the sources of color charge. The McLerran-Venugopalan model gives the gluons in the wavefunction of the nucleus as the Weizsäcker-Williams field generated by the valence quarks of the nucleus. At RHIC energies this may be a not bad approximation. At LHC energies quantun evolution in the longitudinal momentum (BFKL evolution) will also be important thus further increasing the gluon numbers by raising the saturation momentum, $`Q_s`$. For gluons having $`k_{}/Q_s`$ on the order of one, which gluons will dominate the initially produced energy in an ion-ion collision, the occupation number $`f_g1/\alpha .`$ The fact that $`f_g`$ never grows larger than a fixed constant times $`[\alpha N_c]^1`$ when $`k_{}/Q_s1`$ is the essence of the idea of gluon saturation and corresponds to the dominance of the initially produced energy density by fields $`F_{\mu \nu }1/g.`$

## 3 Relationship with HERA phenomenology

### 3.1 The Golec-Biernat Wüsthoff model

Many of the properties of the small-x wavefunction of a high-energy heavy ion are also present in high-energy protons and have been extensively studied at HERA. There is some evidence for saturation at HERA, although most of the evidence is indirect and makes use of models of the Golec-Biernat Wüsthoff type.

In the Golec-Biernat Wüsthoff model one views virtual photon-proton scattering in terms a quark-antiquark dipole, coming from the $`\gamma ^{},`$ scattering on the proton. QCD evolution is put into the proton so that the scattering is that of a bare dipole, having no evolution, colliding with a highly evolved proton wavefunction. The situation is illustrated in Fig.3 while the formulas are

$$F_2=\frac{Q^2}{4\pi ^2\alpha _{em}}d^2x_{}_0^1𝑑z\underset{F}{}e_F^2|\psi _T(x_{},z,Q)|^2\sigma _0(1e^{x_{}^2\overline{Q}_s^2/4})$$
(13)
and

$$F_2^D=\frac{Q^2}{4\pi ^2\alpha _{em}}d^2x_{}_0^1𝑑z\underset{F}{}e_F^2|\psi _T(x_{},z,Q)|^2\frac{1}{2}\sigma _0(1e^{x_{}^2\overline{Q}_s^2/4})^2+q\overline{q}g\mathrm{part}$$
(14)
for $`F_2`$ and the diffractive part of $`F_2,`$ respectively.

The various factors in (13) have ready interpretations. The $`\frac{Q^2}{4\pi ^2\alpha _{em}}`$ is the factor connecting $`F_2`$ to the scattering cross section of an (off-shell) transversely polarized virtual photon on a proton. $`\psi _T`$ is the lowest order wavefunction for a photon of virtuality $`Q`$ to go into a quark-antiquark pair having transverse coordinate separation $`x_{}`$ and with longitudinal momentum fractions $`z`$ and $`1z`$ for the two members of the dipole. $`e^{x_{}^2\overline{Q}_s^2/4}=S`$ represents the impact parameter averaged S-matrix for the elastic scattering of a dipole on a proton while $`\sigma _0`$ represents the size (area) of the proton.

There are three parameters in the model

$$\sigma _0=2.3mb,x_0=3x10^4\mathrm{and}\lambda =0.3$$
(15)
where

$$\overline{Q}_s^2=(x_0/x)^\lambda ,$$
(16)
in units of $`GeV^2,`$ with $`\overline{Q}_s^2`$ having the interpretation of an impact parameter averaged quark saturation momentum of the proton where the gluon and quark saturation momenta are related by

$$\overline{Q}_s^2=\frac{C_F}{N_c}Q_s^2=4/9Q_s^2.$$
(17)
The interpretation of (14) goes much the same way. There is also a term involving the scattering of a quark-antiquark-gluon system on the proton, but this introduces no new parameters.

With the three parameters of the model given as in (15) the G-B W model does a good job in fitting the moderate $`Q^2`$ and low $`x`$ HERA data. Saturation plays an important role because the dipole-proton S-matrix

$$S=e^{x_{}^2\overline{Q}_s^2/4}$$
(18)
becomes zero when $`x_{}2/Q_s`$ thus limiting the contribution of large dipole configurations in (13) and (14) to that allowed by unitarity. This limitation is especially important in diffractive scattering which in perturbation theory has a strong divergence coming from large dipole sizes.

In the HERA range the value of $`\overline{Q}_s^2`$ from (16) is between 1 and 2 $`GeV^2`$ giving $`Q_s^2(24)GeV^2.`$ As we shall see later this value may be a little large and a more realistic model can be expected to give somewhat lower values of $`Q_s^2.`$

### 3.2 Going to impact parameter space

One can get information on the S-matrix for dipole-proton scattering as a function of impact parameter by using $`\gamma `$ proton $`J/\psi `$ proton and $`\gamma `$ proton $`\rho `$ proton quasielastic scattering data. The process is pictured in Fig.4 for the production of longitudinally polarized $`\rho `$-mesons at momentum transfer $`\mathrm{\Delta }.`$ To analyze the process one also needs a model of the $`\rho `$-meson wavefunction which is taken from earlier phenomenological analyses. In the approximation that the dipole-proton scattering amplitude is purely imaginary one can write the amplitude in terms of the square root of the differential cross section. It is straightforward to arrive at

$$S(b,x_{},x)=1\frac{1}{2\pi ^{3/2}N}d^2\mathrm{\Delta }e^{i\mathrm{\Delta }b}\sqrt{\frac{d\sigma }{dt}}$$
(19)
where

$$N=(\psi _{\rho _L},\psi _\gamma ^{})$$
(20)
is the overlap of the wavefunctions for a longitudinal virtual photon with a longitudinal $`\rho `$meson. $`S`$ depends on the impact parameter, $`b,`$ on the Bjorken $`x`$variable, and on the transverse size of the dipole pair, $`x_{},`$ which is fixed, approximately, in terms of $`Q.`$

In order to interpret (19) in terms quantities we have used in the Golec-Biernat Wüsthoff model write $`S`$ as in (18). (This form for $`S`$ is not generally true, but it is not unreasonable when gluon number densities are large.) From fitting the data at $`x=5x10^4`$ and $`Q^2=3.5GeV^2(x_{}0.4fm)`$ one gets

$$\overline{Q}_s^2(b=0)0.5GeV^2\mathrm{or}Q_s^2(b=0)1.1GeV^2.$$
(21)
Although the values given in (21) have considerable uncertainity, mainly from deciding what value of $`x_{}`$ to assoicate with a given $`Q^2,`$ I think they are reasonable. They are somewhat smaller than those that come from the Golec-Biernat Wüsthoff model. Using (11) to write

$$\overline{Q}_s^2(b,x)=\frac{2\pi ^2\alpha }{N_c}\frac{dxG}{d^2b}$$
(22)
we are led to

$$\frac{dxG}{d^2b}|_{b=0}5/fm^2\mathrm{and}xG6$$
(23)
for the values of $`x`$ and $`Q^2`$ quoted above.

### 3.3 From HERA to RHIC

Now let’s be really adventurous and scale (23) back to RHIC energies, $`x0.01,`$ by using the scaling law

$$\overline{Q}_s^2(x)=(\frac{x_0}{x})^{0.3}\overline{Q}_s^2(x_0)$$
(24)
found in the Golec-Biernat Wüsthoff analysis. This gives

$$xG|_{RHIC}2.5$$
(25)
a little larger than most perturbative QCD fits, but not a crazy number at all.

Now that we have $`xG`$ for the proton at RHIC energies we can use (9) and (10) to get the corresponding gluon density and saturation momentum. For Au at RHIC energies one finds from (9), (10) and (25)

$$Q_s^2|_{b=0}1.25GeV^2\mathrm{and}\frac{dxG_A}{d^2b}|_{b=0}5/fm^2$$
(26)
not too different from our estimate given in Sec.2.1 for the values of the typical momentum and number density of gluons produced in a heavy ion collision. We shall come back a little later to the question of the relationship between the number density of gluons in the wavefunction of an ion just before collision and the number density of initially produced gluons.

## 4 What dynamics can (should) one use?

For the wavefunction of a heavy ion, and at RHIC energies, the McLerran-Venugopalan model should be adequate. However, by LHC energies, gluon evolution will certainly become important and may even be significant at RHIC energies. The density of gluons in the wavefunction (26) and initially produced (2-4) is so large at RHIC that the idea of viewing the wavefunction and the production process in terms of hadrons and in terms of nucleon-nucleon collisions, respectively, does not make sense. The collision must be viewed in terms of quarks and gluons.

The McLerran-Venugopalan model uses classical gluon dynamics with random color sources given by the valence quarks. Because gluon occupation numbers are reasonably large the freeing (production) of gluons in an ion-ion collision should also be amenable to a description in terms of classical gluon dynamics. So far there have been two important approaches to this problem.

Krasnitz and Venugopalan have developed a numerical procedure for calculating the early stages of a heavy ion collision in terms of discretized classical Yang-Mills equations. The results of these calculations will be discussed shortly. Kovchegov has tried to calculate the gluon production cross section by looking at Feynman graphs. His procedure for gluon production in proton-nucleus and nucleus-proton collisions gives reliable results. (Nucleus-proton scattering is not manifestly the same as proton-nucleus scattering in Kovchegov’s calculation because of the use of a particular light-cone gauge.) It is still an open question whether his procedure is valid for nucleus-nucleus collisions, although the results which come out are quite reasonable.

### 4.1 Krasnitz, Nara, Venugopalan calculation

Starting with a McLerran-Venugopalan wavefunction as an initial condition Krasnitz, Nara and Venugopalan numerically solve discretize Yang-Mills equations as the nuclei in a heavy ion collision overlap. Gluons appear to be freed on a time scale $`\tau 3/Q_s.`$ Just before the collision the distribution of gluons in one of the nuclear wavefunctions is given by (8) and with gluon occupation numbers given by (12). Just after the gluons are freed the distribution of liberated gluons fits well to the empirical formula

$$\frac{dN_g^{in}}{dyd^2b_{}d^2p_{}}=\frac{N_c^21}{4\pi ^3}\frac{1}{\alpha N_c}\frac{0.11}{e^{\sqrt{p_{}^2+m^2}/T_{eff}}1}$$
(27)
with $`T_{eff}1.1Q_s`$ and $`m0.08Q_s`$ giving

$$f_g^{in}\frac{1}{\alpha N_c}\frac{0.11}{e^{p_{}/T_{eff}}1}$$
(28)
so long as $`p_{}/T_{eff}`$ is not too small.

Comparing (8) to (27) one sees that $`c`$, the fraction of gluons freed in the collision, obeys $`c1/2.`$ In the wavefunction $`<p_{}/Q_s>0.6`$ while in the initial distribution (28), $`<p_{}/Q_s>1.5`$ so that the freeing process seems to create a lot of additional transverse energy which is not present in the ion wavefunction. One worrisome result of the calculation is that the occupation number, given by (28), is quite small so that classical evolution cannot be expected to be reliable after the gluons have been freed and one might even worry about the validity of the freeing calculation being done at a classical level.

### 4.2 Classical field theory versus kinetic theory

When $`f_g1/\alpha `$ kinetic theory does not apply and classical field theory is the right approach. Thus in the earliest stages of a heavy ion collision, during the time of the freeing of the gluons and immediately after, classical field theory is the natural dynamics to use. There may be a practical problem, however, because the Krasnitz-Nara-Venugopalan calculation suggests that there is a rapid lowering of the occupation number as the gluons are being liberated. Thus, although there is a $`1/\alpha `$ factor in (28) the actual value of $`f_g`$ is not large in any practical circumstance. If it should turn out that the classical calculation is not valid numerically during the freeing of gluons the correct dynamics would then be a full quantum field theoretic calculation starting from just before the collision up to the time gluons have been freed, a very difficult calculation.

After gluons have been freed and in the regime (perhaps purely a formal regime) where $`1<f_g<1/\alpha `$ classical field theory and the Boltzmann equation are equivalent procedures for calculating the time evolution of the system.

When $`f1`$ classical field theory, which only keeps the $`f_g`$ in statistical factors $`(1+f_g)`$ is an underestimate of the interation strength of gluons while the Boltzmann equation is systematic.

The “bottom-up” picture of evolution uses the Boltzmann equation, including elastic and inelastic collision terms, and follows the gluon system from shortly after gluons have been produced up to the time when equilibration occurs.

### 4.3 “Bottom-up” picture;Kharzeev, Levin,Nardi picture

The “bottom-up” picture of thermalization applies to a collision of heavy ions where the freed transverse energy is dominated by gluons having transverse momenta well into the perturbative scale. At RHIC energies the transverse momenta which dominate are in the region of one GeV and so one might expect the picture to be qualitatively and even semiquantitatively correct but precise calculations are always going to be hard to come by. The picture does not require that the initial state of the system have large occupation numbers but we shall describe it in the context of the initial conditions being given by gluon saturated heavy ion-wavefunctions.

Consider a head-on high energy heavy ion collision. Our description will refer to gluons in the central unit of rapidity in the center of mass system. The time over which the gluons are freed is $`\tau 1/Q_s.`$ This is also the time scale on which the higher rapidity gluons separate, in the longitudinal spatial direction, from the gluons in the central unit of rapidity with which we are concerned. We presume that a finite fraction, c, of the gluons in the initial wavefunction have been freed in the collision as discussed in sec.2.3 giving initial gluon occupation numbers on the order of $`1/\alpha .`$

As $`Q_s\tau `$ grows the occupation number decreases due to elastic scattering during the longitudinal expansion of the system. When $`f_g`$ is significantly less than $`1/\alpha `$ the Boltzmann equation, the basic dynamics in the bottom-up picture, begins to be applicable. When $`Q_s\tau `$ gets as large as $`\alpha ^{3/2}`$ the occupation number for gluons becomes less than one.

When $`\alpha ^{3/2}<Q_s\tau <\alpha ^{5/2}`$ the dominant process is soft gluon emission from the hard gluons. During this period the density of soft gluons (having $`p_{}\sqrt{\alpha }Q_s)`$ is catching up with the density of hard gluons (having $`p_{}Q_s`$), although both densities are decreasing in time due to longitudinal expansion.

When $`Q_s\tau \alpha ^{5/2}`$ the density of soft gluons reaches that of hard gluons. At later times two important changes happen. (i) The soft gluons stimulate further emission of soft gluons so that now the number density of soft gluons grow with time in spite of the longitudinal expansion. (ii)The soft gluons equilibrate amongst themselves with the hard gluons serving as a source of energy causing the temperature of the soft gluons to grow with time as $`T\alpha ^3Q_s\tau .`$ The rate of flow of energy from hard gluons continues to grow until, at $`Q_s\tau \alpha ^{13/5},`$ the hard gluons disappear, having lost all their energy, and one is left with a completely equilibrated system. At later times the temperature decreases $`\tau ^{1/3}`$ as the equilibrated system continues to expand. The various stages in the equilibration process are shwon in Fig.5. Finally, it is important to note that it is probably overly optimistic to believe that these separate regimes are sharply distinct at RHIC energies.

Now let us try to simply characterize charged hadron production at RHIC. Take a central collision. Then at mid-rapidity the initial gluon density in the McLerran-Venugopalan model is

$$\frac{dN_g^{in}}{dyd^2b}=c2\sqrt{R^2b^2}\rho xG$$
(29)
where $`c`$ is the fraction of gluons liberated. Eq.(29) is easily obtained by integrating (8) over $`k_{},`$ using (9), and inserting the gluon freeing factor $`c.`$ Now at late times the experiments measure the rapidity distribution of charged hadrons. From (29) we get

$$<\frac{2}{N_{part}}\frac{dN_{ch}}{dy}>=\frac{2}{3}RcxG(x,<Q_s^2>)3.8$$
(30)
where $`N_{part}`$ is the number of participants, the factor $`2/3`$ accounts for the fraction of charged particles and R is an inelasticity factor as the partonic system evolves to equilibrium. (R also accounts for any, presumably small, change in particle numbers as quarks and gluons go to hadrons.) The experimental value of 3.8 corresponds to $`\sqrt{s}=200GeV.`$

Kharzeev, Levin and Nardi take a very simple picture where $`R=1`$ so that $`\frac{2}{3}cxG3.8.`$ The KLN analysis was the first detailed analysis of particle production at RHIC using the saturation (color glass condensate) picture and so it has great historical importance. I believe it is probably difficult to get $`cxG`$ as large as 5.7.

The Krasnitz, Nara, Venugopalan picture is also a saturation picture where, however, one of the primary achievements was to calculate $`c`$ and find that $`c1/2`$ at RHIC energies. Moderately large values of $`R`$ and $`xG`$ are necessary to agree with the data.

The “bottom-up” picture does not give a precise value for $`R,`$ and of course there is no calculation of $`c`$ or of $`xG.`$ However, given values of $`R,c`$ and $`xG`$ one does get values (estimates) for the time at which equilibration occurs and the temperature at that time. For $`R3,x1`$ and $`xG2`$ one gets $`T_{equil}230MeV`$ and $`\tau _{equil}3.6fm,`$ values which seem reasonable.

The KLN, KNV and “bottom-up” pictures are all based on the same ideas, and differ only in the focus of the calculation and in the various details.

## 5 Surprises, puzzles, connections

### 5.1 Gluon liberation, proton-nucleus vs nucleus-nucleus

Consider gluon emission in quark-nucleus scattering where the coherence time of the gluon is larger than the nuclear radius. The process is best viewed in the rest system of the nucleus and is shown in Fig.6 in lowest order. The graph in part (a) of the figure corresponds to a gluon in the wavefunction of the quark as the quark reaches the nucleus. Graph (b) corresponds to gluon emission after the quark passes the nucleus. The calculation is straightforward to perform. The square of graph (a) liberates all gluons having $`k_{}^2/Q_s^21`$ while the square of graph (b) gives the “shadow” of graph (a) producing gluons with a $`dk_{}^2/k_{}^2`$ spectrum for $`k_{}^2/Q_s^21.`$ The interference terms do not contribute when $`k_{}^2/Q_s^21`$ and they cancel the direct terms for $`k_{}^2/Q_s^21.`$ The calculation can be generalize to an incident proton on a nucleus with a similar result. Thus one gets

$$\frac{dN_g}{dy}=2xG_{proton}(x,Q_s^2)$$
(31)
corresponding to a gluon liberation factor of 2. The fact that the gluon liberation factor is 2 rather than 1 corresponds to the “shadow” term and has a quantum mechanical interpretation analogous to the result of $`\sigma =2\pi R^2`$ for the scattering of a particle off a hard sphere in non-relativisstic quantum mechanics.

Now, what is surprising is that the Krasnitz, Nara, Venugopalan calculation gives a much smaller number, $`c1/2.`$ The lack of a “shadow” term in heavy ion collisions is perhaps natural. Then we are left with trying to understand why $`c`$ is near to $`1/2`$ rather than 1. In the case of proton-A scattering most of the freed gluons coming from the (a) term in Fig.6 have a momentum, $`k_{},`$ in the proton’s wavefunction obeying $`k_{}/Q_s1,`$ parametrically for large $`Q_s.`$ For A-A collisions the freed gluons in the wavefunction have $`k_{}/Q_s`$ on the order of one so perhaps it is not so surprising that some of the gluons having $`k_{}`$ on the order of $`Q_s`$ are not freed. This indeed, might be expected if it were not for the fact that the freeing process, as described in Sec.4.1 creates so much transverse momentum.

### 5.2 $`v_2(p_{})`$

One of the strongest indications of an early equilibration of the quark gluon plasma in ion-ion collisions at RHIC is elliiptic flow. Detailed calculations using the equations of ideal hydrodynamics suggest that thermalization must occur by $`\tau 0.6fm`$ in order to get a sufficiently large elliptic flow. On the other hand it appears extremely difficult to get such early equilibration from perturbative QCD; the interaction strength does not seem strong enough to cause equibibration at such early times. One can see the difficulty very simply. To match RHIC data the thermalized quark-gluon plasma at $`\tau 0.6`$ has a temperature about $`330MeV`$ corresponding to a typical quark and gluon momentum of about $`1GeV.`$ Perturbative calculations should be reliable to test whether gluons having $`2GeV`$ can be in equilibrium. The mean free path for a $`2GeV`$ gluon to have a collision with momentum transfer $`2GeV`$ is on the order of 10fm at $`T=330MeV.`$ Thus it is extremely difficult to see how 2 GeV gluons can be thermalized at $`\tau 0.6fm.`$ On the other hand gluons of $`500MeV`$ are still in the process of being produced at $`\tau 0.6fm.`$ Indeed, Krasnitz and Venugopalan find that the dominant energy density only becomes produced when $`\tau 3/Q_s.`$ Thus, while it is not impossible that an effective thermalization of gluons in the regime $`0.51GeV`$ could occur at $`\tau 0.6fm`$ it would be puzzling from the QCD dynamics point of view.

### 5.3 Saturation and $`p_{}`$broadening

Experiments on $`\mu `$pair production in proton-A reactions observed a broadening in the transverse momentum, $`p_{},`$ spectrum of the $`\mu `$pair. On the theoretical side one expects the change in the width of the $`p_{}^2`$spectrum to be

$$\mathrm{\Delta }p_{}^2=\frac{2\pi ^2\alpha }{N_c}L\rho xG$$
(32)
where L is the length of nuclear material the annihilating antiquark passes through, $`\rho `$ is the nuclear density, and $`xG`$ is the nucleon gluon distribution. With $`xG1`$ (32) is in reasonable agreement with data. Now recalling (9) and (17) we see that (9) and (32) are really the same formula. Thus we can get another formula for $`Q_s^2`$ as

$$Q_s^2(b)=\frac{N_c}{C_F}2\sqrt{R^2b^2}\frac{d\mathrm{\Delta }p_{}^2}{dz}.$$
(33)
Now take the experimental number

$$\frac{d\mathrm{\Delta }p_{}^2}{dz}0.02GeV^2/fm$$
(34)
for the broadening of the $`\mu `$pair spectrum at fixed target energies. This gives

$$Q_s^2|_{b=0}0.6GeV^2.$$

This number is about a factor of 2 smaller than what one might expect at RHIC, but at RHIC energies $`\frac{d\mathrm{\Delta }p_{}^2}{dz}`$ may well be significantly larger than at fixed target energies. A comparison of the $`\mu `$pair spectrum of proton-proton and proton-nucleus reactions could give a good direct measurement of $`\overline{Q}_s^2`$ while a comparison of the $`J/\psi `$ production transverse momentum spectrum could give a good direct measurement of $`Q_s^2.`$


format: text

text:
# Competition among molecular fragmentation channels described with Siegert channel pseudostates

## Abstract

To describe multiple interacting fragmentation continua, we develop a method in which the vibrational channel functions obey outgoing wave Siegert boundary conditions. This paper demonstrates the utility of the Siegert approach, which uses channel energy eigenvalues that possess a negative imaginary part. The electron scattering energy in each such channel is rotated upward, giving it an equal and opposite imaginary part. This permits a natural inclusion of vibrational continua without requiring them to appear as explicit channels in the scattering matrix. Calculations illustrate the application of this theory to photoionization, photodissociation, and dissociative recombination.

In this Letter, we propose a general method for describing coupling between electronic and dissociative continua, based on a Siegert pseudostate basis representation of the vibrational degree of freedom. The underlying rationale for this idea traces back to the original recognition by Kapur and Peierls Kapur1938  that the narrow resonances of a scattering spectrum can be described in terms of a complex energy eigenstate, with the imaginary part of the energy defining a resonance width parameter. This proposal was further developed by Siegert Siegert1939  in two important ways. First, Siegert demonstrated that a set of eigenstates of any Hamiltonian could be chosen to satisfy pure incoming or outgoing wave boundary conditions in the asymptotic limit. Second, and of importance to our work, Siegert’s derivation allowed for overlapping resonances of arbitrary width, giving a smooth background term in the cross-section. Siegert eigenstates formally correspond to S-matrix poles in the complex plane. Sharply resonant features associated with bound states can be identified with poles lying on the real axis, while broad background scattering can be described by closely spaced eigenstates with finite imaginary parts that serve as a discretized approximation to the true continuum.

Traditionally, the use of Siegert states has been complicated by the nonlinearity of the associated eigenvalue problem. Because the wavenumber appears linearly in the boundary condition but quadratically (as the energy) in the eigenvalue, the eigenproblem is quadratic, and in the past could only be solved iteratively. Tolstikhin *et al.* Tolstikhin1998 ; Tolstikhin1997  recently demonstrated how this difficulty may be circumvented using finite range *Siegert pseudostates*.

Initially, the true asymptotic boundary condition is replaced by a finite range approximation,

$$\left(\frac{d}{dR}ik\right)\varphi (R)|_{R=R_0}=0,$$
(1)
where $`R_0`$ is a value beyond which the value of the potential is negligible. We seek a solution expanded in terms of some primitive basis set

$$\varphi (R)=\underset{j=1}{\overset{N}{}}c_jy_j(R),0RR_0.$$
(2)
Here $`N`$ is the dimension of our basis, and we have selected a non-orthogonal B-spline basis for the $`y_j(R)`$. Inserting this into the Schrödinger equation premultiplying by $`y_j^{}`$, and employing the boundary value (1), we find a matrix equation for the coefficients $`c_j`$

$$\begin{array}{cc}& \underset{j=1}{\overset{N}{}}(\frac{1}{2}_0^{R_0}\frac{dy_j^{}}{dR}\frac{dy_j}{dR}dR\frac{ik}{2}y_j^{}(R_0)y_j(R_0)\hfill \\ & +_0^{R_0}y_j^{}(R)\mu [V(r)E]y_j(R)dR)c_j=0.\hfill \end{array}$$
(3)
Note that we have used a Green’s theorem identity before substituting in the boundary condition, and that the Hamiltonian has been multiplied through by the reduced mass $`\mu `$. Written more concisely in matrix notation, we have a system of the form

$$(\stackrel{~}{𝐇}ik𝐋k^2𝐎)\stackrel{}{c}=0,$$
(4)
where $`L_{j,j^{}}`$ is the surface matrix $`y_j(R_0)y_j^{}(R_0)`$, $`\stackrel{~}{H}_{j,j^{}}`$ is the matrix $`2\mu H_{j,j^{}}+y_j(R_0)\frac{d}{dr}y_j^{}(R_0)`$, and $`𝐎`$ is the overlap matrix for the spline basis set.

This equation is manifestly nonlinear, but the method of Tolstikhin *et al.* allows it to be “linearized” by recasting it as a new eigensystem in a basis of doubled dimension. (Friedman1956 ; Huestis1975  discuss related techniques for solving differential equations where the eigenvalue appears in a boundary condition.) We define $`d_i=ikc_i`$, yielding a trivial second equation $`ik𝐎\stackrel{}{c}=𝐎\stackrel{}{d}`$. Substituting this into the original eigenequation now gives a linear equation in the doubled basis space

$$\left(\begin{array}{cc}\stackrel{~}{𝐇}& 0\\ 0& 𝐎\end{array}\right)\left(\begin{array}{c}\stackrel{}{c}\\ \stackrel{}{d}\end{array}\right)=ik\left(\begin{array}{cc}𝐋& 𝐎\\ 𝐎& 0\end{array}\right)\left(\begin{array}{c}\stackrel{}{c}\\ \stackrel{}{d}\end{array}\right).$$
(5)
This is an equation for the eigenvalue $`\lambda =ik`$, giving $`2N`$ solutions lying either on the (Re $`\lambda `$)-axis or in conjugate pairs in the right half of the complex plane.

In their work, Tolstikhin *et al.* used completeness properties of the Siegert state set to construct a Green’s function of the Hamiltonian, a scattered solution, and the associated scattering matrix, for a variety of single channel model problems. For systems with multiple channels, we instead appeal to the well-understood machinery of multichannel quantum defect theory (MQDT).

For resonance series corresponding to high electronically excited intermediates (Rydberg states) of diatomic molecules, the most natural description of the system is one with quantum defect parameters defined in terms of a fixed internuclear distance $`R`$ and a well-defined projection of the orbital angular momentum $`\mathrm{\Lambda }`$ onto the axis of symmetry. This is because the electron spends most of its time far from the nuclear core, and when it does penetrate into the core, it gains enough speed from falling through the Coulomb potential that the nuclei are essentially frozen on the time scale of its motion. The quantum defect functions $`\mu _\mathrm{\Lambda }(R)`$ in this representation, the so-called “body-frame”, may either be calculated from highly accurate *ab initio* techniques, or extracted from a semi-empirical fitting of experimental data Jungen1997 . In order to connect them with the true asymptotic ionization channels defined in terms of Siegert pseudostates of the residual core, $`j=\{v^+,N^+\}`$, a *frame transformation* must be performed Fano1970 ; Greene1985 , where $`N^+`$ is the ionic rotational momentum, and $`v^+`$ is the vibrational quantum number of the pseudostates. In our procedure, we directly evaluate the S-matrix by the frame transformation integral

$$\begin{array}{cc}\hfill S_{j,j^{}}=& \underset{\mathrm{\Lambda }}{}N^+|\mathrm{\Lambda }_0^{R_0}\varphi _j(R)e^{2i\pi \mu _\mathrm{\Lambda }(R)}\varphi _j^{}(R)𝑑R\mathrm{\Lambda }|N^+^{}\hfill \\ & +i\underset{\mathrm{\Lambda }}{}N^+|\mathrm{\Lambda }\frac{\varphi _j(R_0)e^{2i\pi \mu _\mathrm{\Lambda }(R_0)}\varphi _j^{}(R_0)}{k_j+k_j^{}}\mathrm{\Lambda }|N^+^{}.\hfill \end{array}$$
(6)
The surface term in (6) also arises in the orthonormality relation Tolstikhin1998 . A similar transformation converts the body-frame transition dipole elements $`D_\mathrm{\Lambda }(R)`$ into reduced dipole matrix elements in the same S-matrix representation,

$$\begin{array}{cc}\hfill D_j^S=& (2J+1)\underset{\mathrm{\Lambda }}{}\mathrm{\Lambda }|J_0^{(J)}\mathrm{\Lambda }|N^+\hfill \\ & \times _0^{R_0}\varphi _0(R)D_\mathrm{\Lambda }(R)e^{i\pi \mu _\mathrm{\Lambda }(R)}\varphi _j(R)dR.\hfill \end{array}$$
(7)
Here $`\varphi _0(R)`$ is the initial vibrational wavefunction, and $`J_0`$ and $`J`$ are the total angular momenta of the initial and final states of the system, respectively. Note that the Siegert pseudostates are *never* conjugated in these expressions, even when they formally belong to the dual (“bra”) space. In particular, this means that the quantity labeled as $`\stackrel{}{D}^S`$ below is calculated by conjugating only $`e^{i\pi \mu _\mathrm{\Lambda }(R)}`$ in the definition above, and not the dipole matrix elements directly.

At this stage of the calculation no information about the long-range behavior of the channels has yet been included, and since the body-frame quantum defects are nearly energy independent, the resulting S-matrix is typically a smooth and fairly weak function of energy. The method of *channel elimination* Seaton1983 ; Aymar1996  systematically eliminates flux in all electronic channels below the energy threshold for electron escape (the “closed-channel subspace”) to form a “physical” S-matrix $`𝐒^{phys}`$, by taking the proper linear combination of short-range solutions that ensures exponential decay at infinity. For a long-range Coulomb potential, this procedure gives

$$𝐒^{phys}=𝐒_{oo}𝐒_{oc}(𝐒_{cc}e^{2i\beta })^1𝐒_{co}.$$
(8)
Here, $`\beta `$ is a diagonal matrix of the usual Coulomb long-range phase parameter $`\pi (\nu _j)`$ where $`\nu _j`$ is the (possibly complex) effective quantum number in the $`j`$th channel, $`𝐒`$ is the scattering matrix, and the subscripts indicate partitions of the matrices into closed and open subspaces Aymar1996 .

For a Siegert state basis, this physical scattering matrix is in general *not* unitary, but rather subunitary, reflecting the loss of flux at the boundary $`R_0`$ via coupling to the Siegert pseudo-continuum states. It may be used to calculate partial cross-sections by means of conventional formulae, but with the departure from unitarity, $`1_j|𝐒_{j,j^{}}^{phys}|^2`$, identified as the probability $`|𝐒_{d,j^{}}^{phys}|^2`$ for scattering into the dissociative continuum. This method also provides all quantities necessary to find the partial photoionization cross-section into any open channel, $`\sigma _j`$; see Eq. 2.59 of Aymar1996  for further details. The contributions from all open channels may then be summed to give the total cross-section for photoionization.

Alternatively, the total photoabsorption cross-section may be found directly from a *“preconvolution” formula* first derived by Robicheaux to handle the energy smoothing of densely spaced resonances Robicheaux1993 ; Granger2000 ,

$$\begin{array}{cc}\hfill \sigma _{total}(E)=& \frac{4\pi ^2\alpha \omega }{3(2J_0+1)}\text{Re}\stackrel{}{D}^S\hfill \\ & \left[\mathrm{𝟏}𝐒e^{2i\beta }\right]^1\left[\mathrm{𝟏}+𝐒e^{2i\beta }\right]\stackrel{}{D}^S\hfill \end{array}$$
(9)
where Re signifies taking the real part of everything that follows, and the $``$ here conjugates only the operator, not the entire matrix element. The diagonal matrix written as $`e^{2i\beta }`$ has a nontrivial definition in terms of the quantum defect parameters, it may be approximated quite well by taking $`\beta _j=i\mathrm{}`$ for “closed” channels with $`(E)<\text{Re}E_j`$, and $`\beta _i=\pi \nu _j`$ for “open” channels with $`(E)>\text{Re}E_j`$. Here $`E`$ is the total energy of the system, $`E_j`$ is the threshold energy for channel $`j`$, and $`\nu _j=1/\sqrt{2(E_jE)}`$ on the branch where Im $`\nu >0`$. The utility of this expression lies in recognizing that the value of the cross-section at a complex energy in the above formula is equivalent to the cross-section at a real energy, smoothed over a channel-dependent width $`\mathrm{\Gamma }_j=2\text{Im}ϵ_j`$. Within the Siegert state formulation, the electron energy $`ϵ_j=EE_j`$ will naturally take on a complex value in any channel where the channel eigenenergy $`E_j`$ is itself complex, while $`E`$ remains real.

Given $`𝐒`$ and $`ϵ_j`$, either of the two cross-section formulae above may be evaluated, with appropriate allowances for the possibility of complex energy eigenvalues. Note that the first procedure simply gives a sum over the flux into specific ionization channels, while the second gives a single value for the total photoabsorption cross-section. This means that the latter will contain information about the solution wavefunction along the $`R=R_0`$ boundary not contained in any of the open ionization channels. In general, the value of $`\sigma _{total}`$ will be equal to *or greater than* the sum over the individual $`\sigma _j`$, and any difference may be attributed to the effect of coupling to high-lying Siegert states in the continuum. Thus, the difference between these two formulae at any energy provides the dissociative cross-section.

In order to test the validity of this hypothesis, we began by defining a set of Siegert pseudostates for the $`H_2^+`$ internuclear potential. The eigensolutions fall into three classes, as shown in Fig. 1. Those lying along the positive (Im $`k`$)-axis are associated with negative eigenenergies on the physical sheet of the $`E`$-plane, the bound states of the potential. These are the channel thresholds to which the Rydberg autoionization series of the ionization spectrum converge, and so we include all of their states. The solutions along the negative (Im $`k`$)-axis lie on the unphysical energy sheet, and we reject them as antibound states arising from the doubling of the dimension space. The remainder of the solutions fall above and below the (Re $`E`$)-axis, corresponding to conjugate solution pairs of the eigenvalue parameter $`\lambda =ik`$. We select only those with negative Im $`E_j`$, ensuring that they obey outgoing wave boundary conditions. (This amounts to selecting only those states that contribute to the retarded, rather than the advanced, Green’s function. See Tanabe2001  for related discussion.) For MQDT matrix elements it is also acceptable to reject states lying very high in the continuum, since their Franck-Condon overlap with the bound states is negligible.

Tolstikhin *et al.* discuss the unusual completeness relation obeyed by the full set of Siegert pseudostates, with an additional factor of $`2`$. Our restricted subset of Siegert pseudostates does not, of course, obey that doubled completeness relation. We have confirmed through numerical tests, however, that our restricted subset behaves like a complete set, to at least $`10^{12}`$ accuracy, for representing either $`L^2`$ functions confined within the boundary or functions with purely outgoing wave character at the boundary. For an impressive demonstration (in a somewhat different context) of the convergence properties of a similarly truncated Siegert basis also used to describe smooth continuum physics, see Seideman1991 .

In the region of the ungerade H2 spectrum between 127200 and 127800 cm-1 there are several strongly predissociated resonances, members of the $`np\pi ,v^+=8`$ and $`np\pi ,v^+=5`$ series. In each case, our calculated spectrum correctly reproduces them in the total absorption cross-section, but shows them as weak or absent in the ionization. Comparisons of our results with other theoretical and experimental values GlassMaujean1987 ; Jungen1997  for the relative yields of selected resonances appear in Table I. Note particularly that we are able to correctly describe the strong rotational dependence of the $`4p\pi ,v^+=5`$ branching ratio, a nontrivial consequence of subtle channel interactions.

As a test of the method in an entirely different energy regime we considered the problem of dissociative photoionization, a three-body breakup channel accessible only at much higher energies. Experimental measures of the ratio between pure ionization and dissociative ionization have been performed since the 1970s by a number of researchers Browning1973 ; Backx1976 ; Chung1993 , along with at least one early theoretical calculation Ford1975 . Since our ionization spectrum is a sum over individual channels, we can easily distinguish between contributions from channels above and below the dissociative threshold. Our results, plotted against those of past experiment and theory, are presented in Fig. 2.

Finally, we have performed a model calculation demonstrating the utility of our method for treating dissociative recombination, particularly in systems where indirect channels (those involving scattering into intermediate autodissociating Rydberg states) play an important role. Fig. 3 shows the dissociative recombination spectrum of a simplified H2 model potential (neglecting rotation and with R-independent quantum defects), compared with the familiar approximation of O’Malley for smooth background scattering by direct processes OMalley1966 . Our spectrum accurately reproduces this background, and also describes complex interference effects from the series of resonances converging to each Rydberg threshold.

Some aspects of the Siegert MQDT method remain poorly understood, and would benefit from greater clarification. For example, the utility of a subset of the Siegert basis for MQDT depends on ability of that subset to represent all energetically accessible regions of the continuum. While this requirement appears from our numerical tests to be reasonably easy to satisfy, we have not yet rigorously derived it from the relevant completeness relations. Also, it is not presently clear how to extend the energy-smoothed formula to include non-Coulombic long-range electronic potentials.

Other avenues of investigation could provide insight concerning the applicability of our method to more complex systems. Polyatomic molecules, for example, might be handled either by reduction to hyperspherical coordinates Kokoouline2001 , or by a multidimensional generalization of the Sierget state boundary conditions on an arbitrary hypersurface. Since our method yields a full solution to the Schrödinger equation along the boundary at R0, it should also be possible to project onto the continuum functions of different dissociative channels, and explicitly resolve partial dissociation cross-sections Yaris1979 . Even in its current form, however, we believe the Siegert MQDT method offers a simple description of the flux escaping into dissociative channels, by working within a channel basis that obeys a physically motivated boundary condition.

This work is supported in part by a grant from the National Science Foundation. We thank B. Esry for assistance in the early stages. Discussion with R. Santra, M. Baertschy, M. S. Child, C. W. McCurdy, J. M. Hutson, T. N. Rescigno, and B. I. Schneider has also been helpful.

TABLE I. Photoionization and photodissociation yields for select *ungerade* resonances in H2 for which the relative yields have been experimentally observed GlassMaujean1987 .

| State | Source | Energy | % Ion. | % Diss. |
| --- | --- | --- | --- | --- |
| $`3p\pi ,v^+=8,R(0)`$ | Observed | 127248.2 | 10(5) | 95(5) |
|  | Theory$`^{\text{Jungen1997}}`$ | 127246.9 | 1 | 99 |
|  | Present | 127242.2 | 1 | 99 |
| $`5p\sigma ,v^+=4,R(0)`$ | Observed | 127599.4 | 90(10) | 10(10) |
|  | Theory$`^{\text{Jungen1997}}`$ | 127602.2 | 88 | 12 |
|  | Present | 127606.8 | 76 | 24 |
| $`4p\pi ,v^+=5,R(0)`$ | Observed | 127667.6 | 82(5) | 18(5) |
|  | Theory$`^{\text{Jungen1997}}`$ | 127665.4 | 93 | 7 |
|  | Present | 127666.6 | 97 | 3 |
| $`4p\pi ,v^+=5,R(1)`$ | Observed | 127599.4 | 30(10) | 70(10) |
|  | Theory$`^{\text{Jungen1997}}`$ | 127758.4 | 17 | 83 |
|  | Present | 127759.5 | 29 | 71 |


text:
# Three-site quantum lattice with thermal bath

## 1 The physical model

The Hosltein hamiltonian for a lattice connected to independent site oscillators can be distinguished into three terms

$`h`$ $`=`$ $`h_a+h_b+h_{ab}`$ (1)
corresponding to the *primary* system $`a`$, representing those states allowed to the quantum particle, the *secondary* system $`b`$, representing the vibrating lattice masses ($`m`$) and the interaction energy $`ab`$ \[Holstein-59a\].

In a three-site lattice, assuming the nearest neighbours approximation hold, the quantum states are described by the hamiltonian

$$h_a=V\left(|12|+|21|+|23|+|32|\right)$$
(2)
where $`|n`$ are the eigenstates of the isolated sites and $`V`$ is the overlap integral between neighbour sites.

Vibrations are described by harmonic oscillators of equal mass $`m`$ and frequency $`\omega `$, with respective position $`r_n`$ and conjugated variable $`p_n`$, (where $`n`$ indicates the location on the lattice):

$$h_b=\frac{1}{2m}\left(p_1^2+m^2\omega ^2r_1^2\right)+\frac{1}{2m}\left(p_2^2+m^2\omega ^2r_2^2\right)+\frac{1}{2m}\left(p_3^2+m^2\omega ^2r_3^2\right)$$
(3)
Interaction terms are assumed to be linear with respect to $`r_n`$

$`h_{ab}`$ $`=`$ $`ϵ_1|11\left|+ϵ_2\right|22\left|+ϵ_3\right|33|`$ (4)
$`ϵ_n=\chi r_n`$ , with $`\chi `$ real positive, being the coupling energy at site $`n`$. Such physical model has been studued by Hennig \[hennig-91\] and Kenkre & Andersen \[kenkre-dimer-86\], but without thermal environment.

### 1.1 Lattice with free boundary conditions

When loose ends boundary conditions apply to the primary system, we are lead to the following equation for the density matrix:

$`\dot{\rho }_{11}`$ $`=`$ $`i\omega _o\left(\rho _{21}\rho _{12}\right)`$
$`\dot{\rho }_{12}`$ $`=`$ $`i\omega _{12}\rho _{12}+i\omega _o\left(\rho _{22}\rho _{11}\rho _{13}\right)`$
$`\dot{\rho }_{13}`$ $`=`$ $`i\omega _{13}\rho _{13}+i\omega _o\left(\rho _{23}\rho _{12}\right)`$
$`\dot{\rho }_{21}`$ $`=`$ $`i\omega _{12}\rho _{21}+i\omega _o\left(\rho _{11}+\rho _{31}\rho _{22}\right)`$
$`\dot{\rho }_{22}`$ $`=`$ $`i\omega _o\left(\rho _{12}+\rho _{32}\rho _{21}\rho _{23}\right)`$ (5)
$`\dot{\rho }_{23}`$ $`=`$ $`i\omega _{23}\rho _{23}+i\omega _o\left(\rho _{13}+\rho _{33}\rho _{22}\right)`$
$`\dot{\rho }_{31}`$ $`=`$ $`i\omega _{13}\rho _{31}+i\omega _o\left(\rho _{21}\rho _{32}\right)`$
$`\dot{\rho }_{32}`$ $`=`$ $`i\omega _{23}\rho _{32}+i\omega _o\left(\rho _{22}\rho _{31}\rho _{33}\right)`$
$`\dot{\rho }_{33}`$ $`=`$ $`i\omega _o\left(\rho _{23}\rho _{32}\right)`$
where $`\omega _o=V/h`$ is the frequency for the oscillation of the quasi-particle between two neighbour sites and $`\omega _{kn}=\omega _{nk}=\left(ϵ_nϵ_k\right)/\mathrm{}=\chi \left(r_nr_k\right)/\mathrm{}`$.

Using operators

$`u_1`$ $`=`$ $`\rho _{11}\rho _{22}`$
$`u_2`$ $`=`$ $`\rho _{22}\rho _{33}`$
$`v_1`$ $`=`$ $`i\left(\rho _{12}\rho _{21}\right)`$
$`v_2`$ $`=`$ $`i\left(\rho _{32}\rho _{23}\right)`$ (6)
$`v_3`$ $`=`$ $`i\left(\rho _{31}\rho _{13}\right)`$
$`w_1`$ $`=`$ $`\rho _{12}+\rho _{21}`$
$`w_2`$ $`=`$ $`\rho _{23}+\rho _{32}`$
$`w_3`$ $`=`$ $`\rho _{13}+\rho _{31}`$
a new set of equations, that will be used to perform numerical simulations, can be obtained for the primary system.

As a fluctuation-dissipation process has been finally attached to each oscillator, the complete dynamical model can be written:

$`\dot{u}_1`$ $`=`$ $`\omega _o\left(2v_1+v_2\right)`$
$`\dot{u}_2`$ $`=`$ $`\omega _o\left(v_1+2v_2\right)`$
$`\dot{v}_1`$ $`=`$ $`\omega _{12}w_1+\omega _o\left(2u_1+w_3\right)`$
$`\dot{v}_2`$ $`=`$ $`\omega _{23}w_2\omega _o\left(2u_2w_3\right)`$
$`\dot{v}_3`$ $`=`$ $`\omega _{13}w_3\omega _o\left(w_1w_2\right)`$
$`\dot{w}_1`$ $`=`$ $`\omega _{12}v_1+\omega _ov_3`$
$`\dot{w}_2`$ $`=`$ $`\omega _{23}v_2\omega _ov_3`$
$`\dot{w}_3`$ $`=`$ $`\omega _{13}v_3\omega _o\left(v_1+v_2\right)`$ (7)
$`\dot{r}_1`$ $`=`$ $`p_1/m`$
$`\dot{r}_2`$ $`=`$ $`p_2/m`$
$`\dot{r}_3`$ $`=`$ $`p_3/m`$
$`\dot{p}_1`$ $`=`$ $`m\omega ^2r_1{\displaystyle \frac{\chi }{3}}\left(c+u_2+2u_1\right)\gamma _1p_1+f_1\left(t\right)`$
$`\dot{p}_2`$ $`=`$ $`m\omega ^2r_2{\displaystyle \frac{\chi }{3}}\left(c+u_2u_1\right)\gamma _2p_2+f_2\left(t\right)`$
$`\dot{p}_3`$ $`=`$ $`m\omega ^2r_3{\displaystyle \frac{\chi }{3}}\left(cu_12u_2\right)\gamma _3p_3+f_3\left(t\right)`$
where $`\gamma _n`$ are the damping coefficients and $`f_n\left(t\right)`$ models gaussian noise, with $`\delta `$-shaped time correlation, satisfying the following relations:

$`f_n(t)`$ $`=`$ $`0`$ (8)
$`f_n(t)f_n(t^{})`$ $`=`$ $`2\gamma \theta \delta _{\left(tt^{}\right)}`$ (9)
having defined thermal energy

$$\theta =k_BT$$
(10)
($`k_B`$being the Boltzmann constant); hence

$$f_n^2(t)^{1/2}=\left(2\gamma \theta \right)^{1/2}$$
(11)
It can be readily verified that

$$c=\rho _{11}+\rho _{22}+\rho _{33}$$
(12)
is a constant of motion. A further constant, that will be used to monitor numerical simulations, is given by

$$K=\frac{4}{3}\left(u_1^2+u_2^2+u_1u_2\right)+v_1^2+v_2^2+v_3^2+w_1^2+w_2^2+w_3^2$$
(13)
The dynamical problem is described by a system of stochastic real valued, first order, differential equations (7) that must be numerically solved.

Notice that, restricting the number of sites to $`n=1,2`$ we are left with the spin-boson equation with noise:

$`\dot{u}_1`$ $`=`$ $`2\omega _ov_1`$
$`\dot{v}_1`$ $`=`$ $`\omega _{12}w_1+2\omega _ou_1`$
$`\dot{w}_1`$ $`=`$ $`\omega _{12}v_1`$
$`\dot{r}_1`$ $`=`$ $`p_1/m`$ (14)
$`\dot{r}_2`$ $`=`$ $`p_2/m`$
$`\dot{p}_1`$ $`=`$ $`m\omega ^2r_1{\displaystyle \frac{\chi }{3}}\left(c+u_2+2u_1\right)+f_1\left(t\right)`$
$`\dot{p}_2`$ $`=`$ $`m\omega ^2r_2{\displaystyle \frac{\chi }{3}}\left(c+u_2u_1\right)+f_2\left(t\right)`$
that has been studied by many researchers; in particular, the author’s work has been inspired the papers\[kenkre-dimer-86, grigo-dimer-91\], referenced to in the bibliography.

## 2 Numerical integration method

Stochastic differential equations of first order in time for a vector variabile $`𝐱(t)`$ take the from,

$$\dot{𝐱}=𝐅+\mathrm{𝐠𝐟}$$
(15)
where $`𝐅=\left(F_i\right)_{1,n}=\left(F_i(𝐱,t)\right)_{1,n}`$ is the vector of deterministic fields, $`𝐟=\left(f_i\right)_{1,n}=\left(f_i(t)\right)_{1,n}`$ the gaussian stochastic force, satisfying the conditions:

$`f(t)`$ $`=`$ $`0`$
$`f(t)f(t^{})`$ $`=`$ $`\delta (tt^{})`$ (16)
while th parameters set $`𝐠=\left(g_i\right)_{1,n}=\left(g_i(𝐱,t)\right)_{1,n}`$ represents the interaction of the system with the thrmal bath.

Equations given in (15) are *entirely coupled,* since, in general, every components $`F_i`$ and $`g_i`$ are functions of the whole variable $`𝐱`$,

The formal solution of (15), assuming the *implicit* time dependence $`F=𝐅(𝐱(t))`$ e $`g=𝐠(𝐱(t))`$, is given by \[spde\]

$`x_i(t)x_i(0)={\displaystyle _0^t}𝑑t^{}F_i+{\displaystyle _0^t}𝑑t^{}g_if_i`$ $`i=1,\mathrm{}n`$ (17)
Given a time interval $`[0,h]`$ that can be considered infinitely small for $`F_i`$ e $`g_i`$ functions, their value at the arbitrary instant $`t^{}[0,h]`$, can be approximated by the first $`\kappa `$ terms of a Taylor expansion, centered at $`t=0`$:

$`F_i\left(t\right)`$ $`=`$ $`F_i^o+\delta ^\kappa F_i`$ (18)
$`g_i\left(t\right)`$ $`=`$ $`g_i^o+\delta ^\kappa g_i`$ (19)
The index $`\kappa `$ in (18,19) indicates the truncation order of the series:

$$\delta F_i=\underset{j}{}\left[F_i\right]_j\delta x_j(t^{})+\frac{1}{2}\underset{jk}{}\left[F_i\right]_{jk}\delta x_j(t^{})\delta x_k(t^{})+\frac{1}{3!}\underset{jk}{}\left[F_i\right]_{jkl}\delta x_j(t^{})\delta x_k(t^{})\delta x_l+\mathrm{}$$
(20)
$$\delta g_i=\underset{j}{}\left[g_i\right]_j\delta x_j(t^{})+\frac{1}{2}\underset{jk}{}\left[g_i\right]_{jk}\delta x_j(t^{})\delta x_k(t^{})+\frac{1}{3!}\underset{jk}{}\left[g_i\right]_{jkl}\delta x_j(t^{})\delta x_k(t^{})\delta x_l+\mathrm{}$$
(21)
in which $`F_i`$ and $`g_i`$, time derivatives evaluated at $`t=0`$, have been represented by the following notation:

$$\left[F_i\right]_j\frac{F_i}{x_j}(0)\left[F_i\right]_{jk}\frac{^2F_i}{x_jx_k}(0)\left[F_i\right]_{jkl}\frac{^3F_i}{x_jx_kx_l}(0)$$
(22)
$$\left[g_i\right]_j\frac{g_i}{x_j}(0)\left[g_i\right]_{jk}\frac{^2g_i}{x_jx_k}(0)\left[g_i\right]_{jkl}\frac{^3g_i}{x_jx_kx_l}(0)$$
(23)
Truncation at $`\kappa =0`$ the general solution obtained after substitution of (18,19) into (17):

$$x_i(h)x_i(0)=_0^h𝑑t^{}\left(F_i^o+\delta ^\kappa F_i\right)+_0^h𝑑t^{}\left(g_i^o+\delta ^\kappa g_i\right)f_i$$
(24)
yields

$$x_i(h)x_i(0)=F_i^oh+g_i^o_0^h𝑑t^{}f_i$$
(25)
where $`F_i=\left[F_i\right]=F_i^o`$ and $`g_i=\left[g_i\right]=g_i^o`$ have been defined.

Since the gaussian integral

$$Z_1(h)=_0^h𝑑t^{}f(t^{})$$
(26)
is of $`h^{1/2}`$ order, the second term in (25) is the lowest order approximation of the trajectory (24):

$$\delta x_i^{(1/2)}=g_i^oZ_1$$
(27)
having defined $`\delta x_i=x_i(h)x_i(0)`$.

Substituting (27) into (20,21), hence into (17), and retaining only contribution of the order $`h`$ in (24), we obtain the first order correction

$`\delta x_i^{(1)}`$ $`=`$ $`F_i^oh+{\displaystyle _0^h}𝑑t^{}{\displaystyle \underset{j}{}}\left[g_i\right]_j\delta x_j^{(1/2)}f_i(t^{})`$ (28)
$`=`$ $`F_i^oh+{\displaystyle \underset{j}{}}\left[g_i\right]_jg_j^o{\displaystyle _0^h}𝑑t^{}Z_1(t^{})f_i(t^{})`$
$`=`$ $`F_i^oh+{\displaystyle \frac{1}{2}}{\displaystyle \underset{j}{}}\left[g_i\right]_jg_j^oZ_1^2`$
The corresponding displacement on the trajectory is thus given by

$$\delta x_i=\delta x_i^{1/2}+\delta x_i^1$$
(29)
Again substituting the first order correction (28) into the series (20,21), hence into (17), then retaining only terms of the order $`h^{3/2}`$ in (24), we get:

$`\delta x_i^{(3/2)}`$ $`=`$ $`{\displaystyle \underset{j}{}}\left[F_i\right]_jg_j^oZ_1h+{\displaystyle _0^h}𝑑t^{}{\displaystyle \underset{j}{}}\left[g_i\right]_j\delta x_j^{(1)}f_i(t^{})`$ (30)
$`=`$ $`{\displaystyle \underset{j}{}}\left[F_i\right]_jg_j^oZ_1h+{\displaystyle _0^h}𝑑t^{}{\displaystyle \underset{j}{}}\left[g_i\right]_j\left(F_i^oh+{\displaystyle \frac{1}{2}}{\displaystyle \underset{k}{}}\left[g_i\right]_kg_k^oZ_1^2(t^{})f_i(t^{})\right)`$
$`=`$ $`{\displaystyle \underset{j}{}}\left[F_i\right]_jg_j^oZ_1h+{\displaystyle \underset{j}{}}\left[g_i\right]_j\left(F_i^oh+{\displaystyle \frac{1}{2}}{\displaystyle \underset{k}{}}\left[g_i\right]_kg_k^{}^o{\displaystyle _0^h}𝑑tZ_1^2(t^{})f_i(t^{})\right)`$
$`=`$ $`{\displaystyle \underset{j}{}}\left[F_i\right]_jg_j^oZ_1h+{\displaystyle \underset{j}{}}\left[g_i\right]_j\left(F_i^oh+{\displaystyle \frac{1}{3!}}{\displaystyle \underset{k}{}}\left[g_i\right]_kg_k^{}^oZ_1^3\right)`$
In the simplified case in which $`g`$ does not depend on $`𝐱`$, the numerical algorithm to reach the order $`h^2`$ in approximating the exact solution of (15) can then be written

$$x_i(h)x_i(0)=g_i^oZ_1+F_i^oh+\underset{j}{}\left[F_i\right]_jg_j^oZ_2+\frac{1}{2}\underset{j}{}\left[F_i\right]_j\left[F_j\right]h^2i=1,\mathrm{}n$$
(31)
or, in vectorial notation:

$$𝐱(h)𝐱(0)=𝐠^o𝐙_1+𝐅^o𝐡+\mathrm{𝐅𝐠𝐙}_2+\frac{1}{2}𝐅𝐅h^2$$
(32)
### 2.1 The algorithm for the three-site lattice

Equations (7) belong to the class (15); therefore the approximation method that has been recalled in the previous chapter can be implemented.

To better expose our procedure, it is convenient switching to the vectorial notation; the system variable $`𝐱`$, taking the form:

$$𝐱=\left[\begin{array}{c}u_1\\ u_2\\ v_1\\ v_2\\ v_3\\ w_1\\ w_2\\ w_3\\ r_1\\ r_2\\ r_3\\ p_1\\ p_2\\ p_3\end{array}\right]$$
(33)
evolves in time according to the equation

$$\dot{𝐱}=𝐅+𝐆$$
(34)
where $`𝐅`$ is the vector of the deterministic forces

$`𝐅`$ $`=`$ $`\left[\begin{array}{c}\omega _o\left(2v_1+v_2\right)\\ \omega _o\left(v_1+2v_2\right)\\ \omega _{12}w_1+\omega _o\left(2u_1+w_3\right)\\ \omega _{23}w_2\omega _o\left(2u_2w_3\right)\\ \omega _{13}w_3\omega _o\left(w_1w_2\right)\\ \omega _{12}v_1+\omega _ov_3\\ \omega _{23}v_2\omega _ov_3\\ \omega _{13}v_3\omega _o\left(v_1+v_2\right)\\ p_1/m\\ p_2/m\\ p_3/m\\ m\omega ^2r_1\frac{\chi }{3}\left(c+u_2+2u_1\right)\gamma _1p_1\\ m\omega ^2r_2\frac{\chi }{3}\left(c+u_2u_1\right)\gamma _2p_2\\ m\omega ^2r_3\frac{\chi }{3}\left(cu_12u_2\right)\gamma _3p_3\end{array}\right]`$ (49)
while $`𝐆`$ contains the average values of the stochastic forces:

$$𝐆=\left[\begin{array}{c}0\\ 0\\ 0\\ 0\\ 0\\ 0\\ 0\\ 0\\ 0\\ 0\\ 0\\ \left(2\gamma _1\theta _1\right)^{1/2}\\ \left(2\gamma _2\theta _2\right)^{1/2}\\ \left(2\gamma _3\theta _3\right)^{1/2}\end{array}\right]$$
(50)
The approximated solution up to the second order in time can be computed according to the formula (32):

$$𝐱(h)𝐱(0)+𝐠𝐙_1+𝐅𝐡+𝐅𝐠𝐙_2+\frac{1}{2}𝐅𝐅𝐡^2$$
(51)
where $`Z_1`$ and $`Z_2`$ are gaussian integrals and the force gradient is given by:

$$𝐅=\left[\begin{array}{cccccccccccccc}0& 0& 2\omega _o& \omega _o& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& \omega _o& 2\omega _o& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 2\omega _o& 0& 0& 0& 0& \omega _{12}& 0& \omega _o& \chi _1w_1& \chi _2w_1& 0& 0& 0& 0\\ 0& 2\omega _o& 0& 0& 0& 0& \omega _{23}& \omega _o& 0& \chi _2w_2& \chi _3w_2& 0& 0& 0\\ 0& 0& 0& 0& 0& \omega _o& \omega _o& \omega _{13}& \chi _1w_3& 0& \chi _3w_3& 0& 0& 0\\ 0& 0& \omega _{12}& 0& \omega _o& 0& 0& 0& \chi _1v_1& \chi _2v_1& 0& 0& 0& 0\\ 0& 0& 0& \omega _{23}& \omega _o& 0& 0& 0& 0& \chi _1v_2& \chi _3v_2& 0& 0& 0\\ 0& 0& \omega _o& \omega _o& \omega _{13}& 0& 0& 0& \chi _1v_3& 0& \chi _3v_3& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 1/m_1& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 1/m_2& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 1/m_2\\ 2\chi _1/3& \chi _1/3& 0& 0& 0& 0& 0& 0& m_1\omega _1^2& 0& 0& \gamma _1& 0& 0\\ \chi _2/3& \chi _2/3& 0& 0& 0& 0& 0& 0& 0& m_2\omega _2^2& 0& 0& \gamma _2& 0\\ \chi _3/3& 2\chi _3/3& 0& 0& 0& 0& 0& 0& 0& 0& m_3\omega _3^2& 0& 0& \gamma _3\end{array}\right]$$
(52)
The implementation of the algorithm according to the rule (51) is displayed in the last section . Ada95 language has been the language of choiche, since it is maintained that a rather easy and dependable parallel code, that could be required to perform long lasting simulations can be worked out in that environment, .

Marsaglia & Tsang \[marsaglia\] fortran routine for random numbers generation has been used to compute $`Z_1`$ and $`Z_2`$ integrals.

## 3 Results and Conclusions

First run of the program allowed to plot variables : $`u_1=\rho _{11}\rho _{22}`$ and $`u_2=\rho _{22}\rho _{33}`$, using the set of physical parameters:

* inter-site interaction $`V=0.09`$
* free oscillators frequencies $`\omega ^2=2x10^3`$
* vibration-quantum state coupling $`\chi =4x10^2`$ ;
* $`\gamma =0.2`$ that imply $`\omega ^2/2\gamma ^2=0.025`$ and $`\chi ^2/2m\omega ^2=0.1`$

Time integration step has been set to $`10^3`$ and the reliability of the computation has been monitored through the constant (13), whose changes have kept below $`10^7`$.

Picture of the site occupation differences $`u_1=\rho _{11}\rho _{22}`$ (lower) and $`u_2=\rho _{22}\rho _{33}`$ (upper) evolving in time (arbitrary units) shows the localized motion of the quantum state.

The possibilty to model the dynamics of a three-site lattice coupled to stochastic oscillators allows one to study how a quantum state propagates in a real lattice, when vibrational motion can be described classically, then to investigate the duration of out of equilibrium states and the conditions under which energy transfer without dissipation may occur.

## 4 Code listing

Here included are the main tokens of code from the iterator procedure, in which the numerical integration at time step dh is performed:

* procedure iterator (x : in out coordinates ; xh : out coordinates ; time\_step : in double ; -- time integration step inner\_cycle : in integer ; -- number of steps between next -- print outs

outer\_cycle : in integer ; -- total number of trajectory -- points to print out

v,am1,am2,am3,gam1,gam2,gam3,om1,om2,om3: in double; chi1,chi2,chi3,t1,t2,t3,e1b,e2b,e3b : in double;

output\_file : in string )

is

Fortran random numbers generation routine rnor is imported here:

function rnor(dum : in Fortran.Fortran\_Integer) return Fortran.double\_precision;

pragma Import (fortran, rnor, "rnor\_");

Definition of numerical parameters:

dh : double := time\_step ; -- time integratin step

dh2 : double := time\_step * time\_step; -- squared time integration step

dh25 := dh2 * 0.5;

Definition of physical parameters:

v2: double := v * v ; --

ak0 := c1\*c1 + c2\*c2 + c3\*c3;

qom1 := om1\*om1;

qom2 := om2\*om2;

qom3 := om3\*om3;

qgam1 := gam1 * gam1;

qgam2 := gam2 * gam2;

qgam3 := gam3 * gam3;

norm := (1.0/3.0)\*( x(1)\*x(1) + x(2)\*x(2) + x(1)\*x(2) );

norm := norm + x(3)\*x(3) + x(4)\*x(4) + x(5)\*x(5);

norm := norm + x(6)\*x(6) + x(7)\*x(7) + x(8)\*x(8);

norm := (3.0/4.0) * norm;

Here the iteration starts:

for i in 1..outer\_cycle loop

for j in 1..inner\_cycle loop

Definition of some useful expressionss:

time := time + time\_step;

w11 := double(rnor(dum));

w12 := double(rnor(dum));

w21 := double(rnor(dum));

w22 := double(rnor(dum));

w31 := double(rnor(dum));

w32 := double(rnor(dum));

z11 := a11\*w11 ;

z12 := a21\*z11 + a22\*w12;

z21 := a11\*w21 ;

z22 := a21\*z21 + a22\*w22;

z31 := a11\*w31 ;

z32 := a21\*z31 + a22\*w32;

om12 := de12 - chi2\*x(10) + chi1\*x(9);

om23 := de23 - chi3\*x(11) + chi2\*x(10);

om13 := de13 - chi3\*x(11) + chi1\*x(9);

qom12 := om12 * om12;

qom13 := om13 * om13;

qom23 := om23 * om23;

akey21 := - chi2\*x(13)/am2 + chi1\*x(12)/am1 ;

akey23 := - chi2\*x(13)/am2 + chi3\*x(14)/am3 ;

akey31 := - chi3\*x(14)/am3 + chi1\*x(12)/am1 ;

ajay1 := chi1\*(ak0+x(1)/2.0+x(2))/3.0;

ajay2 := chi2\*(ak0+x(1)/2.0-x(2)/2.0)/3.0;

ajay3 := chi3\*(ak0-x(1)-x(2)/2.0)/3.0;

omx3 := om23 + om13;

omx1 := om12 + om13;

Variables incrementation:

xh(1) := x(1) + 2.0\*v\*( 2.0\*x(7) + x(6) ) * dh +

v2\*(x(2)-2.0\*x(1)+3.0\*x(5)) * dh2 +

v\*(2.0\*om23\*x(4) - om12\*x(3)) * dh2;

xh(2) := x(2) - 2.0\*v\*( x(7) + 2.0\*x(6) ) * dh +

v2\*(x(1)-2.0\*x(2)-3.0\*x(5)) * dh2 +

v\*(2.0\*om12\*x(3) - om23\*x(4)) * dh2;

xh(3) := x(3) + ( om12\*x(6) + v\*x(8) ) * dh -

qom12\*x(3) * dh25 + v2\*(x(4)-x(3)) * dh25 +

v\*(om12\*x(2) + omx1\*x(5)) * dh25 +

akey21\*x(6) * dh25 ;

xh(4) := x(4) - ( om23\*x(7) + v\*x(8) ) * dh -

qom23\*x(4) * dh25 + v2\*(x(3)-x(4)) * dh25 +

v\*(om23\*x(1) - omx3\*x(5)) * dh25 +

akey23\*x(7) * dh25 ;

xh(5) := x(5) - ( om13\*x(8) + v\*x(6) + v\*x(7) ) * dh -

qom13\*x(5) * dh25 + v2\*(x(1)-x(2) - 2.0\*x(5)) * dh25 +

v\*(omx1\*x(3) - omx3\*x(4)) * dh25 -

akey31\*x(8) * dh25 ;

xh(6) := x(6) - ( om12\*x(3) - v\*x(2) - v\*x(5) ) * dh -

qom12\*x(6) * dh25 - v2\*(5.0\*x(6)+3.0\*x(7)) * dh25 -

v\*omx1\*x(8) * dh25 -

akey21\*x(3) * dh25 ;

xh(7) := x(7) + ( om23\*x(4) - v\*x(1) + v\*x(5) ) * dh -

qom23\*x(7) * dh25 - v2\*(3.0\*x(6)+5.0\*x(7)) * dh25 -

v\*omx3\*x(8) * dh25 -

akey23\*x(4) * dh25 ;

xh(8) := x(8) + ( om13\*x(5) - v\*x(3) + v\*x(4) ) * dh -

qom13\*x(8) * dh25 - 2.0\*v2\*x(8) * dh25 -

v\*(omx1\*x(6) + omx3\*x(7)) * dh25 +

akey31\*x(5) * dh25 ;

xh(9) := x(9) + (x(12)/am1) * dh -

qom1\*x(9) * dh25 -

dh25 * (ajay1 + gam1\*x(12))/am1 +

(D1/am1) * z12 ;

xh(10) := x(10) + (x(13)/am2) * dh -

qom2\*x(10) * dh25 -

dh25 * (ajay2 + gam2\*x(13))/am2 +

(D2/am2) * z22 ;

xh(11) := x(11) + (x(14)/am3) * dh -

qom3\*x(11) * dh25 -

dh25 * (ajay3 + gam3\*x(14))/am3 +

(D3/am3) * z32 ;

xh(12) := x(12) + D1\*z11 -

( am1\*qom1\*x(9) + ajay1 + gam1\*x(12) ) * dh +

qom1\*(am1\*gam1\*x(9) - x(12)) * dh25 +

(chi1\*v\*x(6) + gam1\*ajay1 + qgam1\*x(12)) * dh25 -

gam1\*D1\*z12;

xh(13) := x(13) + D2\*z21 -

( am2\*qom2\*x(10) + ajay2 + gam2\*x(13) )\*dh +

qom2\*(am2\*gam2\*x(10) - x(13)) * dh25 +

(-chi2\*v\*(x(7)+x(6)) + gam2\*ajay2 + qgam2\*x(13)) * dh25 -

gam2\*D2\*z22;

xh(14) := x(14) + D3\*z31 -

( am3\*qom3\*x(11) + ajay3 + gam3\*x(14) )\*dh +

qom3\*(am3\*gam3\*x(11) - x(14)) * dh25 +

(chi3\*v\*x(7) + gam3\*ajay3 + qgam3\*x(14)) * dh25 -

gam3\*D3\*z32;

for ii in 1..config\_space loop

x(ii) := xh(ii);

end loop;

end loop;

end loop;