A new numerical method to solve non-linear coupled differential equations for various field theory models

Nuclear Physics A461 (1987) North-Holland. Amsterdam

539-552

A NEW

NUMERICAL

METHOD

COUPLED

DIFFERENTIAL

NON-LINEAR FOR

VARIOUS

FIELD

THEORY

U. POST, J. KUNZ

for the solution

MODELS*

Giessen, D-6300 Giessen,

Received 17 February (Revised 1 September

We present a new method

SOLVE EQUATIONS

and U. MOSEL

Znstitut fib Theoretische Physik, Universitiit

Abstract:

TO

West Germany

1986 1986)

of the coupled

differential

equations

which have to

be solved in various field-theory models. For the solution of the eigenvalue problem a modified version of the imaginary time-step method is applied. Using this new scheme we prevent the solution from running into the negative-energy sea. For the boson fields we carry out a time integration with an additional damping term which forces the field to converge against the static solution.

Some results are given for the Walecka

model and the Friedberg-Lee

model.

1. Introduction During the last decade there has been an increasing interest in field-theory models for the description of nuclei with a nucleon substructure as well as for hadrons with a quark substructure. thus yield

field coupled coupled

These

a relativistically

models invariant

start from a covariant theory le6). They

lagrangian

all are based

density

and

on a fermion

to one or more boson fields. Thus these models lead to highly non-linear

differential

and a boundary

equations,

that consist of an eigenvalue

value problem

for the meson

problem

for the fermions

fields.

There exist various approaches to solve these equations. Koeppel and Harvey ‘) convert the eigenvalue problem for the fermions into a boundary-value problem by a scaling

transformation.

algorithm.

This method

The

it does exist. A shortcoming the scalar meson coupling Therefore,

coupled

system

has the advantage

is then

that a solution

is the fact that the physical constant

gs, are obtained

solved

by a collocation

is directly

parameters,

calculated,

only as a result of the calculation.

one has to iterate the input values until the physical

values are reached.

With this method it becomes very difficult to find solutions when several states are occupied, because one obtains a separate g, for each state. This difficulty

has been

circumvented

eigenvalue

E of the Dirac hamiltonian

well. Then

E(T) is the solution

l

Work supported

by BMFI

by Dodd that depends

of an additional

and GSI Darmstadt.

0375-9474/87/$03.50 @ Elsevier Science Publishers (North-Holland Physics Publishing Division)

B.V.

if

most important

fermion

and Lohe “) who start with an on the radial coordinate,

differential

equation

r, as

de/& = 0.

540

U. Post et al. / Numerical method

The authors fixed meson

of refs. 2-3*9) alternatingly field (by inward

and the full boundary-value consistency

is reached.

solve the full eigenvalue

and outward problem

integration

and subsequent

with a fixed fermion-source

All these various methods,

problem

for a

matching)

term until

self-

efficient as they may be for systems

with spherical symmetry, cannot easily be applied to deformed systems. Reinhard et al. 13) use a modified version of the latter scheme. They transform the Dirac eigenvalue problem into an effective Schrodinger equation which is solved iteratively. In this method, however, it may happen that the solution runs into the negative-energy states, if the positive and negative parts of the spectrum are not clearly separated. In the nonrelativistic case for deformed systems the imaginary time step method lo) and the related gradient method 14) turned out to be very powerful numerical schemes. We therefore think it is both promising and essential to extend this method to the relativistic domain. Here the negative-energy part of the Dirac hamiltonian obviously requires modifications. In sect. 2 our method for solving the eigenvalue problem is described. Sect. 3 shows how the boundary value problem for the boson fields is solved. Finally, in sects. 4 and 5 we present some results for finite nuclei obtained from Walecka’s model 179)and calculations for the Friedberg-Lee soliton bag model 3V4). 2. Solution of the eigenvalue problem The solution

of the eigenvalue

problem HJ,=EJI

(2.1)

is carried out by modifying the imaginary time step method described by Davies et al. lo), that consists of the repeated operation of the time-evolution operator e-iA’.H’fi on a trial wave function I&,, using an imaginary time step At = - 23, and a subsequent renormalization of &+, = e-iAfH’h &. This is an application of Mises’ algorithm for the evaluation of that eigenvector of a given matrix that corresponds to the eigenvalue of maximum absolute value, as we show in appendix A. One can see immediately, that we have to find a modification

of the imaginary

time-step

method,

if we want

to evaluate e.g. the lowest particle state of a Dirac hamiltonian, because there exists also a negative part of the spectrum corresponding to the antiparticle states. This problem can be circumvented easily by using A=e-

S(H-d”))*

(2.2)

This defines together with eq. (A.l) a new algorithm that leads us to that eigenvalue to a rough estimate of the of H which is closest to e(O). Therefore, E(O)corresponds eigenvalue we want to calculate exactly. Analogously it is also possible to calculate any excited state by choosing a larger value of E(‘I . Obviously, the convergence rate is determined by those eigenstates that correspond to eigenvalues in the region around the desired value.

U. Post et al. / Numerical method In

practical

quarks

applications

and nucleons

of such a method

fixed boson

H usually

the hamiltonian

which obey their own equations

to relativistic

of motion.

fields in H. The solution

depends

In this section

to the full coupled

541

mean field theories

with

also on bosonic

fields

we treat only the case of problem

can be obtained

by an alternating execution of the iteration steps for the solution of the boson field equation (see sect. 3) and the fermion wave function. This method circumvents the need for an exact solution of the eigenvalue and boundary value problems before self-consistency is reached. In order to develop a computer code one has to discretize the hamiltonian in a first step. In particular this includes a special choice of a computational grid, which may be spherical, cylindrical or Cartesian, the choice of certain differential formulas for the spatial derivatives and some boundary conditions. This discretized hamiltonian, Hd, has a finite set of eigenvalues &i, Q, . . . , E, for a total number of n grid points. Not all of these eigenvalues necessarily correspond to physical solutions. For example, there are eigenstates on a spherical grid that do not exhibit the proper asymptotic behaviour -r’ where 1 is the angular momentum of the state. Furthermore, if HA is not strictly hermitian there may exist eigenvalues with nonvanishing imaginary part. Additionally, the spectrum can contain also real eigenvalues with an oscillating behaviour of the corresponding eigenvector that are also artefacts of the discretization, their occurence can be circumvented by a proper choice of the discretization (appendix B). The numerical evaluation of A=e-

S(H-F(O))2

=jgoyqH_B(o))2j

(2.3)

requires that one has to cut off the infinite sum at some value j = N. This results in an approximation AN to A. We define the function (2.4) Then we obtain

for the eigenvalues hk

Therefore,

we have to determine IL@,,S,N(E)I < 1

hk of A,: =~(“),S,N(Ek)

*

6 and N from the convergence for~E[~l,~,]-{~(o)}.

(2.5)

condition (2.6)

The convergence rate becomes the faster the more IfE(~),G,N(~)Iis peaked around E = e(O). We now have to clarify how the convergence rate depends on the values of N and 6. For the exact exponential formula, eq. (2.2), it becomes faster if 6 > 0 increases. For a finite cut off for the sum one then has to use a large value of N in order to fulfil (2.6). This requires not only a higher amount of computing time but

542

U. Post et al. / Numerical method

also numerical rounding errors become more important. some optimized combination of S and N. If we start desired

from any arbitrary

one are suppressed

wavefunction

by the iteration

all other

scheme

x2= S(& -E(O))%

Therefore,

j&k- &(O)(to be sufficiently

iteration is roughly proportional calculation is proportional to

except

)

(xi;;)’

the

(2.7) of &Jx)

with *l (e.g.

.si - E(‘))~} .

Therefore, within one iteration step an eigenstate &kf e(O) is suppressed by a factor of approximately _~W)(,~_“(O))Z e , if we assume

components

as long as

where x&:2 can be obtained from fig. 1 by the intersection x(l) 1.41). Then the optimum value of S is max=Ji= 6(N) = (x~~~)2/maxi{(

one has to find

(2.8)

corresponding

small. The computing

to N, and, hence,

to an eigenvalue

time required

for one

the total time needed

for the

N N maxi{( ei - E(O))~} a(N) (C3’ ’

(2.9)

in order to obtain a certain level of accuracy. Looking at fig. 1 we observe that the important prefactor N/(x~~~)* has its minimum for N = 1. Therefore, we have chosen this value for our calculations and determine 6 from

s < l/(Ei

-2

- E(O))2

-1 Fig.1.L,G,Nfor N=l,...,

fori=l,n.

0

(2.10)

tl 6 and E=O, S=l.

t2

U. Post et al. / Numerical

The choice

N = 1 yields an algorithm

for the solution

of the nonrelativistic

the single spectrum

particle

Hartree-Fock

of this operator

543

method

which corresponds Hartree-Fock

hamiltonian

equations

replaced

has a lower bound

to the gradient

method

in ref. r4), but with

by (H - E(‘))~. Again,

and is even positive

the

definite.

As we have seen before the convergence rate is limited by the maximum possible choice of 6, which itself is limited by the high-lying eigenstates of the discretized hamiltonian. However, the convergence can be accelerated by introducing an energy dependence of 6 to suppress admixtures of the high-lying eigenstates of (H-&(O))2 during

the iterations.

This is done by analogy A = 1 - 6(~ -

to ref. 14) if we replace

do))2

) A=l-(T_E(~))2+~:rJ(H-& (0) 2

7

I$‘““‘) - Al t,b(“)),

(2.11)

where T=-icrV is the relativistic

kinetic

energy

operator,

and

IT = 1 - Il+P’)( I/?(“‘[ projects

onto the subspace

the rate of convergence

orthogonal similar

to I$I(“)), and ISo is a parameter

to 8. However,

the optimum

that determines

choice of E. does not

depend on the mesh size. Using T instead of the full Dirac hamiltonian is sufficient because T yields the main contribution to the energy eigenvalues for the high-lying eigenstates. Furthermore, one has to carry out the time consuming matrix inversion only once. For higher-dimensional grids the inversion of T becomes difficult due to the higher dimension of the corresponding matrix. A damping algorithm, e.g. as described in ref. 14), may be used in order to suppress the high-lying eigenstates. As for the nonrelativistic

case 14) these

modifications

decrease

the number

of

iterations needed to achieve a fixed accuracy of the solution up to about one order of magnitude. The initial wavefunction may be chosen using various methods: Because eq. (2.2) yields a fast convergence, if the level density close to the approximate value is rather small, the choice of 1+4~ is not critical in most cases. Thus, it is sufficient to use a reasonable guess for the wavefunction. E.g. one may choose some analytical ansatz for the radial part that already exhibits the correct node structure. But even that is not a necessary requirement.

544

U. Post et al. / Numerical

Another

method,

up to about

which is only practical

100 (e.g. spherical

method

for a rather small number

symmetry),

consists

fixed values of the other fields. This may be carried out using standard routines.

Then those eigenvectors

of grid points

of the diagonalization

of HA for

diagonalization

are selected that exhibit the desired nodal structure

and belong to an eigenvalue close to the desired one. Naturally, this does not lead to the final solution if self-consistency of several field equations is required. However, we thus obtain initial wavefunctions and values for a(‘) that differ not too much from the final ones and can be iterated simultaneously with the other fields. 3. Solution of the meson field equation In actual

calculations

contains a coupling equation:

for relativistic

mean field problems

term to some mesonic

A4(r)

=f(&

the Dirac hamiltonian

field 4 that in turn obeys the Klein-Gordon

~(r)+~oofo~I~I+~,

p),

(3.1)

where !P represents all other fermion or meson fields, that are assumed to be fixed at first. It is straightforward to solve this coupled problem by a method that iterates between the solutions of eqs. (3.1) and (2.1). Because we do not need the exact solutions neither of the eigenvalue nor of the boundary-value problem at any intermediate step of the iteration we have chosen another way by solving the time-dependent equation d*@(r,

at2

t)

t)-f(@,

=A@(r,

IV)-~7.

@Wr,t)

It differs from eq. (3.1) by the additional time dependence on the left-hand and especially by the friction term -Ad@/dt. A > 0 is some arbitary constant.

side

Since we treat here only static problems, the time t = n At in eq. (3.2) really denotes step n in the iteration scheme. Advancing t by Ar takes much less time than a full solution of the static problem of eq. (3.1) in each iteration and is, furthermore, not restricted to one-dimensional problems. How are the solutions

of (3.1) and (3.2) related @(r, t) = 4(r)+

with ~(r, t) assumed

to be small, we obtain

d2&(r t) d=As(r, at2

t)-

v-(4,

T)

w

E(r,

to each other?

Making

t) ,

(3.3)

a new equation

of motion

as(r t) E(r, t)--h ----+O(EZ). rlt

Within a first approximation we neglect not only the second-order assume As(r, t) to be small. Thus, we obtain

a24r,t) = a-(+, w at*

a4

the ansatz

.c(r, t)-h----

Wr, t) at ’

for E (3.4) term but also

(3.5)

U. Post et al. / Numerical

This corresponds

to a damped

harmonic

545

method

oscillator

equation

for ~(r, t) for any

fixed r if we require

v-(4, w>. 84 This condition

is obviously

fulfilled

(3.6)

*

whenever

f(4, q’) = b(w1* 4

(3.7)

with a non-zero effective mass m( ?P). Therefore, the solution of eq. (3.5) vanishes everywhere for t + co. In order to obtain fast convergence A is adjusted to meet the aperiodic case, i.e. A = 2m. For our calculations we always used a constant value, but one might use also a spatially varying A(r) > 0. Taking the term AE into account yields a coupling between the oscillators at neighbouring positions. In order to see its effect on the time evolution of E( r, t) let us assume, for example, that it has a maximum at r,, for a fixed time t = to. That means for the second order derivative A&(r,,, t,)
a*E(ct) =Ae(r, at*

The general

solution

is the superposition

t)-A

Mr, t) . at

of the solutions

(3.8) rp,(r, t) given by

cpk(r, t) = eik’r+ol*, a=

-;Ak&2-k2.

(3.9)

For A > 0 it can be seen eaisly that the real part of (Y is always negative except the case k = 0. But the solution (P,_~ does not contribute to c( r, t), if one requires @(r, t) to fulfil the same boundary conditions as c,b(r), i.e. for the electromagnetic field @(r, t) = Q/jrl. Th is means that E(r, t) + 0 for 1rl + ~0. Hence, we have also in this case .c(r, t) + 0 for t + a~ and are, therefore, able to solve also the Poisson equation for the electromagnetic field with this method.

546

U. Post et al. / Numerical

The numerical exact solution

solution

E(r, t = 0). Therefore,

accurate

method

discretization

of eq. (3.2) is a minor

Q(r, t) converges

method

problem

it is not necessary

for the time integration.

for the approximation

to spend

due to the fact that the

4(r)

to the static solution

for any initial

much

effort on an extremely

We have chosen

of the first-(second-)

deviation

a simple

2-(3-)

point

order time derivative

of

@(r, t).

4. Application

to the Walecka

model

Although as in the non-relativistic case the main strength of the numerical method developed here is in its application to two- and three-dimensional problems, we illustrate its feasability and properties here only for systems with spherical symmetry. As the first example designed to describe Z=

-l&a,

we treat the lagrangian given by Walecka ‘) that has been nuclear matter properties as well as finite systems. It reads + M)l// -;((a&)‘+

m$,42) -&,F,,

- gsrcr*~ + ~gv~YJk*~, 3

-$n&$OP (4.1)

where Fw,, =

a,~,

-&so,.

(4.2)

In the mean field approximation, 4 and w+ are treated as classical fields and I,G is a Dirac spinor describing the nucleon wavefunctions. M is the free nucleon mass and m, = m, = 783 MeV. The coupling constants gs and gv and the scalar mass m, are adjusted to reproduce binding energy and saturation density of nuclear matter and to the rms radius of 40Ca [ref. “)I. The field equations derived from (4.1) are

a,w,=o.

(4.3)

Starting from the last equations we will describe now in detail how a complex solution containing several occupied states can be constructed. As an example we have chosen the spherical nucleus 40Ca with 6 occupied states. The calculation was carried out along the following steps: (i) Using a solution for I60 as an initial guess the occupation numbers were increased by a factor of two in order to obtain meson fields which are strong enough so that also bound Id,,,, Id,,, and 2s,,, states do exist. (ii) With these fixed 4 and w. fields Ho was diagonalized for K = +l, +2, -3. (iii) The proper eigenvectors corresponding for a fully selfconsistent iteration.

to the occupied

states then are used

U. Post et al. / Numericai method

547

(iv) In the final step a smaller mesh size Ar = 0.06 fm was chosen in order to achieve a better accuracy of the final solution. If the iteration scheme defined by eq. (2.11) and a discretization excluding fermion doubling are applied, this calculation requires approximately 25 seconds of CPUtime on an IBM 3084 in order to achieve an accuracy of about 0.01% of the binding energies. In fig. 2 the results of the calculation for 40Ca are shown. The parameters for this calculation correspond to the mean field values of Walecka and Serot 9, (parameter

Fig. 2. Solution for ‘%a in the Walecka model. The scalar and the are drawn versus r[fm] in the first row. In the lower part of the figure (fun line: upper component; dashed line: lower component). The are ‘): mS = 2.79 fm-‘, m, = 3.97 fm-‘, gs = 9.573, gv = 11.67 and the

vector density and the meson fields the 6 occupied states are displayed parameters which have been used free nucleon mass M = 4.767 fm-‘.

548

U. Post et al. / Numerical method

set of “QHD

I”). The density distribution

and the rms radius where

another

interaction

(3.01 fm) are similar

set of parameters

as well as the binding but differ slightly

including

energy (7.80 MeV/A)

from the results of ref. ‘)

the p-meson

and the electromagnetic

was used.

5. Application to the Friedberg-Lee soliton bag model Friedberg et al. “) and Goldflam et al. ‘) use a model lagrangian that is based on a self-interacting scalar field coupled to quarks to describe the baryons. It is given by ~=-~(r,a,+M)rCr-l(a,a)2-

U(a)-g&a,

(5.1)

where the potential U(a)

=&T4+~ba3+;ca2+p

(5.2)

represents the fourth order self-interaction of the scalar field, (+. The potential is chosen to exhibit two minima. a, b and c are free parameters and p is determined from the constraint CJ(ov) = 0 for the vacuum value a,. The convergence condition (3.6) is fulfilled only in the region around both minima. However, the algorithm for the solution of the Klein-Gordon equation converges too, due to the second-order spatial derivatives. The field equations derived from (5.1) are (YP%+M+gsa)@=O, (aP& - m&r = +g,lJ+. We obtain

for the spherically

symmetric EU =

du

+-+-

dr

du Ev=-dr-

(5.3)

ansatz: 1-K

r

v+(M+gsc+bJ,

l+K -u-(Mfg,a)u, r

(5.4) Within this model calculations have been performed using a = 236.13, b = -11 614, c = 180 000, g, = 25 and M = 0 [ref. “)I. In fig. 3 the results for the zero-node s,/=, d s/29 G/2, p112 and &I2 states are shown. the If,,, and the 2s1,= states are P3/2, unbound. 6. Conclusion Based imaginary

on the observation time-step method,

that in non-relativistic as well as the gradient

Hat-tree-Fock theory the method with additional

U. Post et al. / Numerical

method

Fig. 3. Various bound states for the Friedberg-Lee model for a = 236.13, b = -11614, c = 180 000, gs = 25 and M = 0 “). The first column shows the vector (full) and the scalar (dashed) density drawn versus r[fm], the second one gives the radial dependence of the scalar field (g#) and the third column contains the upper (full) and lower component (dashed) of the fermion wavefunction.

U. Post et al. / Numerical method

550

refinements, physical

have been proven

system

is deformed.

case. This generalization presence

to be an appropriate We have generalized

avoids

the instability

of the negative-energy

demonstrated

by applying

tool, even when the

this method

of the original

states. The feasibility

it to a number

numerical

to the relativistic method

due to the

of the new method

of one-dimensional

(spherically

has been symmetric)

cases ranging from the cases described above to a chiral bag model lagrangian 15) and one that contains explicit coupling of scalar and vector meson fields “). We have used also this method on a cylindrical grid to prepare the spherical initial states for a collision problem. Calculations for non-spherical systems using a constraint on the quadrupole moment are in progress.

Appendix A We show here that the imaginary time-step method is a direct application of Mises’ algorithm which is used for the evaluation of the eigenvalue of maximum absolute value and the corresponding eigenvector: Let A be a square matrix of order n, with eigenvalues Al,. . . , A, and x,, some arbitrary vector with n components. We define a sequence of vectors x, by

,.

x,+1 =

&a+1

A%,,

xm+l=m

with Z,, = x0.

We assume IAll > lhil for i = 2,3,. . . , n. Because A is fully diagonalizable express x0 as an expansion of the eigenvectors y, , y2, . . . , y,: x0=

i k=l

(A-1) we can

(A-2)

ffkyk,

where A,,

=

AkYk

.

(A.3)

Then we have n

%n =c

ffkA:Yk

k=l

(A.4) Using Ihk/hll < 1 for k = 2,3, . . . , n we obtain that the sum on the right-hand of eq. (A.4) becomes negligibly small for m + 00. Therefore &n+y1 x+,Ax, + Al, if (or # 0. Of course,

the first line of eq. (AS)

side

(m+a), (A.5) is only valid up to a phase factor.

U. Post et

Coming

back to the imaginary

al. / hkmericd

time-step

551

method

method

which uses

,4 = e-6.H eq. (AS) means that we will obtain of the Hamilton

operator

(A.@

Ar = e-“& 1, where now E, is the lowest eigenvalue

H. Appendix B

It is demonstrated

that

the discretized

Dirac-hamiltonian

HA may

have

real

eigenvalues with corresponding eigenvectors that exhibit an oscillating behaviour. These eigenvectors, therefore, do not represent a physically meaningful solution and their existence is due to the chosen discretization. As an example we treat the case of a fermion state in a spherical potential well V(r). For the Dirac spinor we use the ansatz 11)

(B-1) and obtain

03.2)

EC= -g-$ii-Mir+V(r)G. The discretization

might be carried _ Et&=

+

&vi=

-

vi+l

out in the following

-;i-l

K

--Ci+MUi+ ri

2Ar iii+1 - Ii-, 2Ar

It can be shown solution

immediately

for the same value

K

--iii-M~i+ ri

V(r,)&. of (B.3) we obtain

another

of E by

exists another discrete It is given by

v”;= -(-)‘& symmetry

rj; = ( - )$ ) which is connected

V(ri)Ci,

that for any solution

;; = ( - )itji ) There vanishes.

form:

of eq. (B.3) if the bare

v”; = ( - )$ )

to (B.4) by the continuous “I Ui=Vi,

-I Vi=-Uiy

_

03.4)

K’=-K.

fermion

K’=K,

mass

(B-9

transformation K’= -_K

For the solution of the eigenvalue problem there are different possibilities rid of this type of oscillating, non-physical solutions. One may either smooth

03.6)

to get u” and

552

U. Post et al. / Numerical

u” or

start from another

set of differential

method

equations,

e.g. those for u = C/I-, II= C/r.

In the latter case an additional

term occurs from the derivative

breaks the energetic

The symmetry

given by eq. (B.4) is already

Hence,

scheme

degeneracy.

if we restrict

ourselves

to a fixed

one solution

and one has to choose the proper

K.

the iteration

(l/r)a(ru)/ar converges

that broken, to only

value of E in eq. (2.2). Furthermore,

the degeneracy is broken for most practical applications by a finite fermion mass or a coupling to any other field, and the oscillating solutions can be suppressed also by some smoothing procedure. The transformation defined by eq. (B.4) is sometimes called fermion doubling. These spurious solutions can also be excluded by using two different equidistant meshes, one for the upper and one for the lower component, respectively, which are shifted

towards

each other by half the step size.

References 1) 2) 3) 4) 5) 6) 7) 8) 9) 10) 11) 12) 13) 14) 15)

J.D. Walecka, Ann. of Phys. 83 (1974) 491 J. Rafelski, Phys. Rev. D14 (1976) 2358 R. Goldflam and L. Wile&, Phys. Rev. D25 (1982) 1951 R. Friedberg and T.D. Lee, Phys. Rev. D15 (1977) 1694; D16 (1977) 1096; D18 (1977) 2623 J. Boguta, Nucl. Phys. A372 (1981) 386 J. Kunz, D. Masak, U. Post and J. Boguta, Phys. Lett. 169B (1986) 133 Th. Koeppel and M. Harvey, Phys. Rev. D31 (1985) 171 L.R. Dodd and M.A. Lohe, Phys. Rev. D32 (1985) 1816 B.D. Serot and J.D. Walecka, Adv. in Nucl. Phys. 16 (Plenum, New York, 1985) K.T.R. Davies, H. Flocard, S. Krieger and M.S. Weiss, Nucl. Phys A342 (1980) 111 J.D. Bjorken and S.D. Drell, Relativistic quantum mechanics (McGraw-Hill, New York 1964) L. Wilets, private communications P.G. Reinhard, M. Rufa, J. Maruhn, W. Greiner and J. Friederich, 2. Phys. A323 (1986) 13 P.G. Reinhard and R.Y. Cusson, Nucl. Phys. A378 (1982) 418 B. Blattel, J. Kunz, U. Mosel and U. Post, GSI scientific report 1985