Dynamical Approach to Pair Production from Strong Fields

arXiv:hep-ph/9209212v1 4 Sep 1992 LA-UR-92-2753 DYNAMICAL APPROACH TO PAIR PRODUCTION FROM STRONG FIELDS Fred Cooper Theoretical Division T-8 Los Alamos National Laboratory MS B-285 Los Alamos, NM 87545 USA Lectures given at the NATO ADVANCED STUDY INSTITUTE Particle Production in Highly Excited Matter II Ciocco, Italy July 1992 8/27/92 1. INTRODUCTION In relativistic heavy-ion collisions one is hoping to produce conditions where energy densities are high enough so that a new state of matter– the quark-gluon plasma can be produced. This state of matter lasts for a short period of time following the collision and may or may not be in equilibrium. Following this phase a transition to ordinary hadronic matter takes place and many of the processes which occur during the quark-gluon plasma phase might be masked by processes which occur in the hadronic phase. In order to determine processes which might be signals of the quark gluon plasma one needs to know the dynamical evolution of the plasma. This is because the particles that get produced during that phase have to travel through a time evolving plasma. In order to study this problem one needs a different way of thinking about field theory. Traditionally experiments in Elementary Particle Physics are black box experiments where initial particles enter a region, final particles exit the experimental region and all that is asked is how many particles of what type, energy, etc. enter various detectors. This type of experiment requires only a covariant S-Matrix theory to predict the probabilities to be expected in the detectors. However, if we want to know signatures of the quark gluon plasma,we actually need to follow the time evolution of the plasma and fields produced following the heavy ion collision. This requires a non-covariant real time formalism for the time evolution of the quantum fields. In these talks we would first like to discuss various formalisms for doing real time calculations in quantum field theory and then study in detail a very simplified model of the production of the quark-gluon plasma– Schwinger’s mechanism for Pair-production from strong “classical” gauge fields. The value of doing a “first principles” calculation at this time, even if it is over-simplified, is multifold: (1) We can test the validity of existing semiclassical transport models of lepton production from the quark-gluon plasma. We have already discovered that these models have to be modified to correctly include Pauli blocking and Bose enhancement effects which were ignored. (2) We can determine the effective hydrodynamics and show that certain kinematic assumptions automatically lead to flat rapidity distributions independent of the form of the equation of state. (3) We can determine the dynamical equation of state and in the next order in a systematic calculation in powers of (1/N) we will be able to study whether equilibration will occur and calculate self consistently lepton pair production rates. 1 First I would like to list the various approaches available to study real time processes in Quantum Field Theory. Each of these approaches needs an approximation scheme to reduce the number of degrees of freedom in order to make the problem numerically tractable. Three methods that my collaborators and I have studied in detail are: 1- Functional Schrodinger Equation + variational approximations [1] [2][3] 2- Truncated Heisenberg Equations – Large-N expansion or Mean-field approximations to the Dyson equations [4] [5][6] [7]. 3- Schwinger’s closed time path - Path Integral Formalism in a large-N expansion [8][9][10][11]. There also exists an alternative formalism related to the truncated Heisenberg equations based on the Wigner Distribution function which has been discussed by Rafelski and his collaborators [12]. Once we have chosen a method we have to decide how to specify the initial data at t=0. In these different approaches we have to specify 1- Initial position and width of say a Gaussian Wave Packet at t=0 in the Schrodinger picture 2- Number and pair densities etc. in Heisenberg picture. 3- Initial density matrix in the Path Integral Approach. In classical field theory, such as classical electrodynamics, the theory is finite and any smooth initial configuration of the field is allowed for the initial value problem. When we have a semiclassical field theory for the expectation value of the fields however, the initial data can be reinterpreted in terms of the particle language and even a smooth initial configuration of the field might not be consistent with certain physical constraints such as the initial state having finite number density at t=0 with respect to an adiabatic vacuum. (This requirement is automatic for finite temperature field theory). Thus arbitrary initial data may not be consistent with renormalizability. This is discussed in detail in [3][4]. We also have an additional new problem to face – how to perform renormalization in a non-covariant formulation of the field theory. To do this we isolate the divergences in an adiabatic (WKB) expansion of Green’s functions. This method is similar to the technique of adiabatic regularization used by Parker and Fulling [13] in their study of semiclassical gravity. The problem we will address in detail in these lectures is pair production of either Bosons or Fermions from strong Classical Fields which are either functions of time t, or fluid proper time τ = (t2 − z 2 )1/2 . We will compare the results of the numerical simulation of this problem (for the degradation of the field, the particle spectra, etc.) with a semi- classical transport approach using a Schwinger-inspired source term [14][15][16]. We will also discuss the effective hydrodynamics derived from the expectation value of the energy momentum tensor of the quantum theory. 2 2. SUMMARY OF THE DIFFERENT STRATEGIES IN λϕ4 FIELD THEORY For simplicity let us first study these different approaches to initial value problems in the simplest case- λϕ4 field theory. a) Schrodinger Picture: In the Schrodinger picture the Initial State is described by a wave functional at t=0. For example a Gaussian wave functional is < ϕ|Ψ > = ψ[ϕ, t] = exp[− Z x,y [ϕ(x) − ϕ̂(x)][G−1 (x, y)/4 − iΣ(x, y)][ϕ(y) − ϕ̂(y)]] (2.1) The time evolution is given by the Functional Schrodinger equation [1]: i∂ψ/∂t = Hψ Z 1 1 H = d3 x[− δ 2 /δϕ2 + (∇ϕ)2 + V (ϕ)] 2 2 (2.2) This is a generalization of the usual Schrodinger equation: ψ(x) =< ψ|x >, x → ϕ(x, t); p = −iδ/δx → π = −iδ/δϕ i∂ψ/∂t = Hψ; H = −∂ 2 /∂x2 + V (x) (2.3) ψ(0) = exp[−α(x − xo )2 ]. (2.4) with initial condition: One might imagine solving (2.2) on a computer by introducing a lattice in d dimensions and converting the functional derivatives into partial derivatives. One then quickly realizes that the number of degrees of freedom in equation (2.2) is rather overwhelming. To control this problem one uses variational trial wave functionals which become “exact” in the largeN limit– namely Gaussians. The equations of motion for the variational parameters can be obtained from Dirac’s variational principle [17]: Z Γ = dt < Ψ|i∂/∂t − H|Ψ > (2.5) δΓ = 0 → Schrodinger’s equation: i∂/∂t − H|Ψ >= 0 3 (2.6) In the ϕ representation one can choose a Gaussian trial wave functional: Z < ϕ|Ψv >= ψv [ϕ, t] = exp[− [ϕ(x) − ϕ̂(x, t)] x,y −1 [G (2.7) (x, y, t)/4 − iΣ(x, y, t)][ϕ(y) − ϕ̂(y, t)]] + iπ̂(x, t)[ϕ(x) − ϕ̂(x, t)] where the variational parameters have the meaning: ϕ̂(x, t) =< Ψv |ϕ|Ψv >; π̂(x, t) =< Ψv | − iδ/δϕ|Ψv > G(x, y, t) =< Ψv |ϕ(x)ϕ(y)|Ψv > −ϕ̂(x, t)ϕ̂(y, t) Then the effective action Z Γ(ϕ̂, π̂, G, Σ) = Z = Z − (2.8) for the trial wave functional is dt < Ψv |i∂/∂t − H|Ψv > dtdx[π(x, t)∂ϕ(x, t)/∂t + Z dtdxdyΣ(x, y)∂G(x, y, t)/∂t dt < H > (2.9) where < H >= Z dx{π 2 /2+2ΣGΣ+G−1 /8+1/2(∇ϕ)2 −1/2∇2 G+1/2V ′′ [ϕ]G+1/8V ′′′′ [ϕ]G2 }. < H > is a constant of the motion and is a first integral of the motion. For λϕ4 field theory we get the following equations of motion: π̇(x, t) = ∇2 ϕ − ∂ < V > /∂ϕ; ϕ̇(x, t) = π Ġ(x, t) = 2 Z dz[Σ(x, z)G(z, x) + G(x, z)Σ(z, x)] Z Σ̇(x, t) = −2 dz[Σ(x, z)Σ(z, x) + G−2 /8 1 + [ ∇2x − ∂ < V > /∂G]δ 3 (x − y) 2 (2.10) If there is translational invariance and ϕ=0 we obtain a second order differential equation for G(k,t), the Fourier transform of G(x,t): 2G̈(k, t)G(k, t) − Ġ2 (k, t) + 4Γ(k, t)G2 (k, t) − 1 = 0 Z 1 2 2 2 2 Γ(k, t) = k + m (t); m (t) = µ + λ [dk]G(k, t) 2 4 (2.11) This approximation is called the time-dependent Hartree-Fock Approximation and is equivalent to the leading term in a 1/N expansion of the field theory [3]. To understand this trial wave function let us look at a simple quantum mechanics problem- the harmonic oscillator with a gaussian initial state. Harmonic oscillator: V(x) =1/2 m x2 , Initial conditions: Ψ(x, 0) = [2πG(0)]−1/2 exp{−x2 /[4G(0)]} q(0) =< x >= 0 (2.12) For the harmonic oscillator a Gaussian remains Gaussian as time evolves so that Ψ(x, t) = (2πG(t))−1/2 exp{−x2 [G−1 (t)/4 − iΣ(t)]} (2.13) We find that the conserved Energy can be written in terms of G as follows: E =< H >= Ġ2 /8G + Gm2 /2 + G−1 /8 = Ġ2 /8G + V [g] (2.14) We plot V[g] in fig 1. From fig. 1 we see that the ground state is G = 1/(2m). If at t=0, G0 = 1/(2M) ; m 6 = M then G(t) = 1/2(G0 + G1 ) + 1/2(G0 − G1 ) cos(2m(t − to )) (2.15) Thus the width oscillates with frequency 2m between G0 and G1 . Generalizing to free field theory (which is just independent harmonic oscillators) we have instead for each mode of momentum k: < H(k) >= Ġ2 /8G + (k 2 + m2 )G/2 + G−1 /8 (2.16) This leads to the same result for G(k,t) as for G(t) with m→ ωk = (k 2 +m2 )1/2 . However in field theory, unlike quantum mechanics, an arbitrary initial Gaussian state is not necessarily a physically valid choice since it might correspond to an infinite particle density or energy density when compared to the adiabatic vacuum. Thus the particle interpretation implies that one needs to restrict the large k behavior of G(k) at t=0 to be a physically allowed initial state with finite particle number, energy density etc. Otherwise one gets extra unwanted infinities in loops. b) Heisenberg Picture: Green’s function approach 5 In problems where there is spatial homogeneity one has a Fourier decomposition for a charged field ϕ in terms of mode functions fk (t) which depend only on the time and the usual creation and annihilation operators a and b which satisfy the canonical commutation relations: Φ(x, t) = Z −ikx [dk][fk (t)ak eikx + fk∗ (t)b+ ] ke + 3 3 ′ [ak , a+ k′ ] = [bk , bk′ ] = (2π) δ (k − k ) (2.17) The initial state |i > is totally specified by specifying at t=0 the matrix elements of a and b: d d ′ < i|a+ k ak |i > = (2π) δ (k − k )n+ (k) < i|bk ak |i > = (2π)d δ d (k + k ′ )F (k)etc. (2.18) The equation for the expectation value of the equation of motion is: < i|(− + m2 )ϕ + λ(ϕ+ ϕ)ϕ|i >= 0 (2.19) We see from these equations that we also need to solve the equation of motion for < i|λ(ϕ+ ϕ)ϕ|i >. In general we get a Heirarchy of Green’s function equations- The BBGKY heirarchy. To make practical progress we need a truncation scheme which allows us to solve the lowest order problem and then systematically calculate corrections. In the large N expansion the lowest order approximation leads to a factorization < i|(ϕ+ ϕ)ϕ|i > =< i|(ϕ+ ϕ)|i >< i|ϕ|i > = G(x, x; t) < i|ϕ|i > (2.20) where the fourier transform G(k,t) of G(x-x’; t) obeys the same equation as the width of the Gaussian wave function in the Schrodinger equation in the Hartree approximation. 2G̈(k, t)G(k, t) − Ġ2 (k, t) + 4Γ(k, t)G2 (k, t) − 1 = 0 Z 1 2 2 2 2 Γ(k, t) = k + m (t); m (t) = µ + λ [dk]G(k, t) 2 Z G(x, x; t) = [dk]G(k, t) 6 (2.21) Thus the large-N expansion (Hartree approximation, mean field approximation) truncates the hierarchy of coupled Green’s function equations making it necessary to only solve the coupled one and two-point Green’s function equations. In these mean field equations the problem reduces to an external field problem in that the quantum field ϕ obeys the equation: (− + m2 (t))ϕ = 0 (2.22) Because we have an external field problem with spatial homogeneity: the mode functions f(t) in (2.17) obey: (∂02 + ω 2 )f = 0; ω 2 = k 2 + m2 (t) (2.23) The canonical commutation relations lead to a constraint on the mode functions: fk f˙k∗ − fk∗ f˙k = i (2.24) which is automatically satified by the WKB form ansatz: fk (t) = [2ΩK (t)]−1/2 exp[−iyk (t)] ẏk (t) = Ωk (t) (2.25) which lead to the equation 3 Ω2k (t) + Ω̈k /(2Ωk ) − (Ω̇k /Ωk )2 = ωk2 (t). 4 (2.26) At t=0 one has in general for the initial state: d d ′ < i|a+ k ak |i >= (2π) δ (k − k )n+ (k) < i|bk ak |i >= (2π)d δ d (k − k ′ )F (k) For an adiabatic vacuum: n(k) =F(k) = 0, and the initial conditions on Ω are Ω(k, t = 0) = ω(k, t = 0); Ω̇(k, t = 0) = ω̇(k, t = 0). (2.27) This formalism, however is perfectly general and one could take any initial state with an integrable phase space particle density n(k) and pair density F(k). As a particular 7 choice one could have chosen at t=0 an equilibrium configuration of pions described by a temperature kT = β −1 n(k) = 1/(exp[βE(k)] − 1) (2.28) c) Path Integral Approach: Closed time-path formalism The only formalism that allows a systematic approach to initial value problems is the closed time-path approach of J. Schwinger[8]which was further elaborated by Keldysh[9]and put into a Path Integral framework by Chou, Su, Hao and Yu[10]. This Path Integral approach allows standard Path Integral approximation schemes such as the large N approximation as well as ensuring causality for the Green’s functions for initial value problems [18]. The starting point for determining the Green’s functions of the initial value problem is the generating Functional: + − ∗ Z[J , J , ρ] =< i|T (exp{− Z iJ− ϕ− })|out >< out|T (exp Z iJ+ ϕ+ )|i > (2.29) This can be written as the product of an ordinary Path integral times a complex conjugate one or as a matrix Path integral. Z + − Z[J , J , ρ] = dϕ+ dϕ− < ϕ+ , i|ρ|ϕ− , i > exp i[(S[ϕ+ ] + J+ ϕ+ ) − (S ∗ [ϕ− ] + J− ϕ− )] Z = dϕα exp i(S[ϕα ] + Jα ϕα ) < ϕ1 , i|ρ|ϕ2 , i > (2.30) where < ϕ+ i|ρ|ϕ− , i > is the density matrix defining the initial state. This leads to the following matrix Green’s functions [11]: G++ = δ 2 ln Z/δJ + δJ + |j=0 =< T (ϕ(x1 ), ϕ(x2 ) > G−− = δ 2 ln Z/δJ − δJ − |j=0 =< T ∗ (ϕ(x1 ), ϕ(x2 ) > G+− = δ 2 ln Z/δJ + δJ − |j=0 =< ϕ(x2 ), ϕ(x1 ) > G−+ = δ 2 ln Z/δJ − δJ + |j=0 =< ϕ(x1 ), ϕ(x2 ) > (2.31) The matrix Green’s function structure insures causality. In this approach it is easy to generate a 1/N expansion in analogy with ordinary field theory. The diagrams are the same as in the usual 1/N expansion, except the Green’s functions are the matrix Green’s functions described above.If in lowest order in (1/N) we have an external field problem as described above, one can directly use the mode solutions of the previous methods to determine the lowest order matrix Green’s function of eq. (2.31). This obviates the need to discuss the initial density matrix of the theory, since it is these Green’s functions which then enter the diagrams of the higher order calculations. 8 3. MAIN EXPANSION IDEA: FLAVOR SU(N) In many problems one of the fields can be treated classically to first approximation– pair production in Strong Electric or Gravitational fields. This makes the lowest order problem an external field problem. One way to generate a systematic expansion whose lowest order is an external field problem is by introducing N copies of the original problem and expanding in Flavor SU(N). This is most easily done in the Path Integral formalism. For the initial value problem one would use the matrix Green’s functions discussed above. Having an extra large parameter N allows an evaluation of the Path integral by Laplace’s method (or the method of Steepest Descent). To obtain the large N expansion one realizes that if there are N flavors the loops carry an extra N. Rescaling the fields then display an overall factor of N in the effective action which includes the loops. Examples: λϕ4 : χ = ϕ2 Z Z Z Z = dχ dϕ exp[− (∂µ ϕ)2 + λχϕ2 − λχ2 + µ2 ϕ + Jϕ + Sχ] ϕ → ϕi , i = 1, 2, · · · N ; λ → λ/N ; ϕi → N 1/2 ϕi ; χ → N χ, λχ → λχ. Integrating over ϕ we obtain: Z 1 Z = dχ exp{−N [χ2 + T r ln G−1 − jGj]} 2 Z = dχ exp{−N Sef f (χ)} G−1 = [− + µ2 + λχ]δ(x − y) (3.1) Evaluating the Path Integral at the Saddle point, δSef f (χ)/δχ = 0 leads to the self consistent external field problem [− + µ2 + λχ]ϕ = 0; χ = ϕ2 + G(xx) (3.3) In QED we obtain an external field problem by integrating out the fermions (which have now N flavors to give an extra N to the determinant) and then rescaling the fields to display the overall factor of N in the effective Action: QED: Z Z Z 1 Z = dAµ dΨdΨ exp[ dx{− F 2 + Ψ(iγ∂ − e 6 A + m)Ψ} + Ψη + ηΨ] 4 √ √ Ψ → Ψi ; e → e/ N , A → A N 9 (3.3) Integrate out the N species of fermions Z dAµ exp{−N Sef f (Aµ )} Z 1 Sef f (Aµ ) = dx F 2 + T r ln(S −1 (x, y; A)) + ηS(x, y; A)η] 4 S −1 (x, y; A) = (iγ∂ − e 6 A(x) + m)δ(x − y) (3.4) Evaluating the Path Integral at the saddle point, δSef f (Aµ )/δAµ = 0 leads to the external field problem: (iγ∂ − e 6 A + m)Ψ = 0 (3.6) where A is an external field,Ψ is a quantum field. We also obtain the semiclassical Maxwell Equation: ∂µ F µλ =< j λ >= e < Ψγ λ Ψ > . (3.7) In all these problems one has in leading order in 1/N a straightforward problem of a quantum field theory in a background field which allows a normal mode decomposition in terms of the solutions of the classical field equations. Renormalization can be carried out by an adiabatic expansion of the mode equation[13]. The effect of quantum fluctuations about the semiclassical field can be systematically taken into account by calculating the fluctuations about the leading stationary phase point in the Path Integral order by order in the 1/N expansion. 4. PARTICLE PRODUCTION IN THE CENTRAL RAPIDITY REGION IN HEAVY ION COLLISIONS A popular picture of high-energy heavy ion collisions begins with the creation of a flux tube containing a strong color electric field[19]. The field energy is converted into particles as qq pairs and gluons which are created by tunnelling- the so-called Schwinger mechanism [20][21][22]. The particle production can be modeled as an inside-outside cascade which is symmetric under longitudinal boosts and thus produces a plateau in the particle rapidity distribution. The boost invariant dynamics, in a hydrodynamical picture gets translated into energy densities (such as E2 ) being functions of the proper time. We take this as an initial condition on the fields in an initial value problem based on this pair-production mechanism. First let us look at the case where the electric field is a function of real time t, treating later the more realistic case where E= E(τ ); τ = (t2 − z 2 )1/2 . Thus we first 10 consider particle production from a spatially uniform electric field such as that produced between two parallel plates. This is an idealized model of a flux tube for QCD. The problem of pair production from a constant Electric field (ignoring the back reaction) was studied by J. Schwinger in 1951 [20]. The physics is as follows: One imagines an electron bound by a potential well of order |V0 | ≈ 2m and submitted to an additional electric potential eEx (as shown in fig. 2 ). The ionization probability is proportional to the WKB barrier penetration factor: exp[−2 Z Vo /e o 4 dx{2m(Vo − |eE|x)}1/2] = exp(− m2 /|eE|) 3 (4.1) A direct calculation due to Schwinger from first principles using the effective action in an arbitrary constant electric field (ignoring the back reaction) gives instead w = [αE 2 /(2π 2 )]Σ∞ n=1 (−1)n+1 exp(−nπm2 /|eE|). 2 n (4.2) This equation tells us that pair production is exponentially suppressed unless eE ≥ πm2 . So we expect (as we find in fig. 3) that there is a crossover value of E where the time it takes for E to first reach zero (remember there are plasma oscillations) is relatively short. Schwinger’s result only applies when we can ignore dynamical photons (as well as back reaction)and is related to the lowest order in 1/N calculation where the electric field is treated as a classical object. Schwinger’s analytical result was subsequently used as source term for an approximate transport theory [14], [15], [16]approach to the back reaction connected with pair production which we will later compare with our exact numerical results. We will choose the electric field in the z direction and choose a particularly simple gauge: → → E = E(t)k̂; A = A(t)k̂; E(t) = −dA/dt (4.3) To maintain spatial homogeneity we have from Maxwell’s equation: ∇. E = ρ (4.4) that the plasma of produced particles must be neutral. In scalar QED, the equation for the quantum field ϕ is −(∂α − ieAα )(∂ α − ieAα )Φ + µ2 Φ = 0 11 (4.5) and for the electromagnetic field: ∂α F βα =< C{−ie(Φ+ ∂ β Φ − Φ∂ β Φ+ ) − 2e2 Aβ Φ+ Φ} > (4.6) where C denotes charge symmetrization with respect to Φ+ and Φ. For our constraints on the field E and our choice of gauge we get: −dE/dt =< jZ >= e Z [dk](kZ − eA(t))G(k, t) where G(k, t) = [< ϕ† ϕ + ϕϕ† > −2ϕ∗ ϕ]F T (4.7a) (4.7b) For QED we have instead the field equation: [iγ∂ − e6 A(t) − m]Ψ(x, t) = 0 (4.8) and the semiclassical Maxwell equation: 1 e < i|[Ψ(x, t), γ3Ψ(x, t)]|i > 2 −dE/dt =< jZ >= (4.9) The fact that the external field is independent of space (spatial homogeneity) means that one has a simple normal mode expansion of the fields just as in λϕ4 field theory described earlier. For Scalar QED we have Φ(x, t) = Z −ikx [dk][fk (t)ak eikx + fk∗ (t)b+ ] ke [∂o2 + ωk2 (t)]fk (t) = 0 2 ωk2 (t) = [k − eA(t)]2 + µ2 + k⊥ (4.10) Repeating the arguments of (2.22 - 2.25) we again obtain for the generalized frequency Ωk (t): 3 Ω2k (t) + Ω̈k /(2Ωk ) − (Ω̇k /Ωk )2 = ωk2 (t). 4 (4.11) where now ω is given by (3.10) Spatial homogeneity requires translational invariance, ′ ′ W (x − x , t, t ) = Z ′ [dk]W (k, t, t′)eik(x−x ) . 12 This in turn requires that d d ′ < a+ k ak > = (2π) δ (k − k )n+ (k) d d ′ < b+ k bk > = (2π) δ (k − k )n− (k); < bk ak > = (2π)d δ d (k + k ′ )F (k) (4.12) Thus we obtain for G(k,t) G(k; t) = Ω−1 (k, t){1 + n+ (k) + n− (k) + 2F (k) cos[2yk (t)]} (4.13) This is the most general form of the propagator that one would use in the diagrams of the 1/N expansion, where n and F are the particle and pair phase space densities at t=0. These parameters also totally specify (in leading order in 1/N) the density matrix at t=0. To solve the field theory in leading order in 1/N (ignoring questions of renormalization to be discussed below) one solves the second order differential equation for each mode function Ωk (t), determines G(k,t) and then solves the back reaction equation: −dE/dt = e Z [dk](kZ − eA(t))G(k, t) (4.14) For QED one has to deal with the spinor structure: ψ(x, t) = Z [dk][uks (t)bk eikx + v−ks (t)d†−k e−ikx ] (4.15) If we choose a basis where γ 0 γ 3 is diagonal: γ 0 γ 3 χs = λs χs , s = 1, 2 → λ = 1; s = 3, 4 → λ = −1 χ†r χs = 2δrs (4.16) Then the spinors u and v obey the equation 0 3 ⊥ ⊥ {γ ∂t + iγ π + iγ p + m}  uks vks (t) (t)  =0 (4.17) Squaring the Dirac equation by letting:  uks vks  0 3 ⊥ ⊥ = {−γ ∂t − iγ π − iγ p + m} 13  χs fk+ (t) − χs f−k (t)  (4.18) we find that the mode functions f now obey: [∂02 + ωk2 (t) − iλs π̇]fk (t) = 0, ωk2 (t) = π 2 + p2⊥ + µ2 π = k − eA (4.19) If the operators ak and bk obey the usual anticommutation relations: {aks , a†k′ s } = {bks , a†k′ s } = (2π)3 δ 3 (k − k ′ )δss′ (4.20) the fk are constrained to satisfy ω 2 f ∗α f β + f˙∗α f˙β + iλs π[f ∗α f˙β − f˙∗α f β ] = δ αβ /2 (4.21) Parametrizing the positive and negative frequency solutions: f± (t) = N± exp Z t g± (τ )dτ, (4.22) 0 we find: g + = −[λs π̇ + Ω̇]/2Ω − iΩ (4.23) where the generalized frequencies, Ωk (t) now satisfies the equation: 3 Ω2k (t) + Ω̈k /(2Ωk ) − (Ω̇k /Ωk )2 − π̇ 2 /(4Ω2 ) − λs π̇ Ω̇/Ω2 4 (4.24) Ignoring renormalization, the solution of QED is obtained by simutaneously solving for these modes and also for E(t) which is obtained from the Maxwell equation: dE/dt = 2eΣ4s=1 Z + [dk](p2⊥ + m2 )λs |fks (t)|2 (4.25) 5. RENORMALIZATION The equations of the previous section as they stand are not finite in the continuum since the sum over modes in (4.14) and (4.25) contains a divergence related to the renormalization of the charge (as well as the wave function) resulting from the charged particle loops in the definition of the current. Let us first look at Scalar QED where the back-reaction equation is: 14 −dE/dt =< j >= e Z [dk](kz − eA(t))Ω−1 [1 + N (k)...] (5.1) We first see that N(k) has to fall fast enough at large k to not lead to any further divergences– this is equivalent to the condition that the initial number density ρ is finite. The integral of Ω−1 contains a divergence proportional to dE/dt which renormalizes the charge (as well as the field E). To isolate this divergence one makes an adiabatic expansion of the equation for the generalized frequencies Ωk . That is, we imagine that the time derivatives are small d/dt→ ǫ d/dt : 3 ǫ2 [Ω̈k /(2Ωk ) − (Ω̇k /Ωk )2 ] = ωk2 (t) − Ω2k (t) 4 and we then expand in powers of ǫ 3 1/Ωk = 1/ωk [1 + ǫ2 {ω̈k /4ωk − (ω̇k /ωk )2 } + 0(ǫ4 ωk−4 )] 8 (5.2) (5.3) We see that terms with higher derivatives are associated with more convergence factors of 1/k so that one only has to consider the first two terms in the adiabatic expansion to isolate the divergences which are interpreted as the standard charge renormalization. The log divergence comes from the term ω̈k = e(dE/dt)(k − eA)ω −1 After integrating over k this leads to a term of the form: Z 1 2 2 2 [dk]ωk−3 = π(0) e δe dE/dt; δe = 12 (5.4) (5.5) where π(0) is the usual vacuum polarization at q2 =0. Subtracting this term from both sides of eq. (5.1) we obtain: 2 2 Z −edE/dt(1 + e π(0)) = e [ [dk](kz − eA(t))G − eπ(0)dE/dt]. (5.6) The Ward identity tells us that eE = eR ER ; and the renormalized charge is determined by e2R = e2 /(1 + e2 π(0)) so the explicity mode by mode finite renormalized equation is Z −dER /dt = eR [dk](k − eA(t))[Ω−1 − ω −1 − e2R (k − eA(t))(dE/dt)ω −5/4] 15 (5.7) (5.8) For QED one gets instead after charge renormalization: Z + 2 4 dER /dt = 2eR Σs=1 [dk][(p2⊥ + m2 )λs |fks | − e2R dER /dt ω −3 (5.9) 6. HEAVY ION COLLISIONS AND BOOST INVARIANT DYNAMICS In e+ e− annihilation, hadronic collisions and in heavy-ion collisions particle production in the central rapidity region can be modeled as an inside-outside cascade which is symmetric under longitudinal boosts which leads to a plateau in the particle rapidity distributions. This boost invariance also emerges dynamically in Landau’s hydrodynamical model [23] and forms an essential kinematic ingredient in the analyses of Cooper, Frye and Schonberg [24] as well as Bjorken[25]. It was recognized by Cooper and collaborators and further elaborated by Bjorken that in a hydrodynamical framework scale invariant initial conditions : v = z/t, ǫ(x, t) → ǫ(τ ), τ 2 = t2 − z 2 (6.1) would automatically lead to flat rapidity distributions. In the context of transport or field theory modelling of the heavy ion collision, after an initial time τ0 , energy densities are expected to be functions only of the fluid “proper time” τ . We therefore assume that the kinematics makes the electric field E only a function of the proper time τ . For this kinematical choice it is convenient to introduce new variables τ, η the fluid “proper time” and the fluid rapidity (when v=z/t) via : z = τ sinh η, t = τ cosh η. (6.2) This change of coordinates to (τ, η) from (t,z) can be accommodated by the usual formalism of curved space [26][27] (except the curvature here is zero). One introduces the metric in curved space gαβ = diag(−1, 0, 0, τ 2). (6.3) (−g)−1/2 ∂β [(−g)−1/2 F αβ ] = j α (6.4) Maxwell’s equations becomes for an electric field E(τ ) in the z direction E(τ ) = Fzt = Fητ /τ = −τ −1 ∂τ A(τ ) (6.5) −1/τ ∂τ [1/τ ∂τ A(τ )] =< j η > (6.6) 16 For Scalar Electrodynamics the equation for χ = √ τ ϕ is (∂τ2 + τ −2 [(∂η − ieA(τ ))2 + 1/4] − ∂x2 − ∂y2 + m2 )χ = 0 (6.7) The rescaled field χ has the same Fourier decomposition as φ had in flat space with the mode functions f obeying [∂τ2 + ωk2 (τ )]fk (τ ) = 0 (6.8) 2 ωk2 (τ ) = [k − eA(τ )]2 /τ 2 + k⊥ + µ2 + 1/(4τ 2 ) (6.9) however now so that the longitudinal momenta get suppressed at large τ . For fermions one has the added complication that the covariant derivative now has a spin piece: (denote the Minkowski indices with α, β the curvilinear coordinates with µν) ∇µ = ∂µ + Γµ − ieAµ 1 Γµ = Σαβ Vαυ Vβυ;µ 2 1 Σαβ = [γ α , γ β ] 4 (6.10) and the vierbein represents the transformation to the Minkowski coordinates: ∼µ gµν = Vµα Vνβ ηαβ ; γ = γ α Vαµ (6.11) Maxwell’s equation becomes: 1 ∼η e < i|[Ψ, γ Ψ]|i > 2 1 e < i|[Ψ† , γ 0 γ 3 Ψ]|i > = 2τ −τ −1 dE(τ )/dτ =< j η > = (6.12) and the fermion mode functions now obey [∂τ2 + ωk2 (τ ) − iλs π̇]fk (τ ) = 0 (6.13) ωk2 (τ ) = π 2 + p2⊥ + µ2 ; π = (p − eA(τ ))/τ (6.14) where The divergences in Maxwell’s equation in curved space can be renormalized as before by an adiabatic expansion in the variable τ . The details of this calculation are presented in [28]: 17 7. PARTICLE PRODUCTION RATES AND THE BOGOLIUBOV TRANSFORMATION The wave functions of the first order adiabatic expansion eikx fk0 (t) where fk0 (t) = (2ωk ) −1/2 exp[−i Z t ωk (t′ )dt′ ] (7.1) 0 form an alternative basis for expanding the scalar fields and allows one to define an interpolating number density N(k,t) which becomes the true one as t→ ∞. Expanding the field in terms of f0k (t) we have Φ(x, t) = Z ∗0 −ikx [dk][a0k (t)eikx fk0 (t) + f−k (t)b0† ]where a0k (t → ∞) = aout etc. k k (t)e (7.2) In this expansion the creation and annihiliation operators are time dependent. We also have our previous expansion in terms of the time independent operators a and b related to the initial state: Φ(x, t) = Z ∗ [dk][fk (t)ak eikx + f−k (t)b†k e−ikx ]. (7.3) We recognize that ak and a0k (t) are related by a unitary transformation. The Bogoliubov coefficients are defined by a0k (t) = α(k, t)ak + β ∗ (k, t)b†k b0k (t) = α(k, t)bk + β ∗ (k, t)a†k |α(k, t)|2 + |β(k, t)|2 = 1 (7.4) The number of particles produced per unit volume is just † out † out ak |t = 0 > bk + aout V −1 dN/dk =< t = 0|bout k k (7.5) The interpolating number density is defined in terms of the first order adiabatic operators: V −1 dN (k, t)/d3 k = ht = 0|b0k † (t)b0k (t) + a0k † (t)a0k (t)|t = 0i = (1 + N (k))|β|2 + N (k)|α|2 + 2Re{αβF (k)} 18 (7.6) For N=F=O (the adiabatic vacuum at t=0) 1 V −1 dN (k, t)/dk = |β|2 = (4ωk Ωk )−1 [(Ωk − ωk )2 + (Ω̇k /Ωk − ω̇k /ωk )2 ] 4 (7.7) We see that adiabatic initial conditions (no particle production at t=0) are Ωk (0) = ωk ; Ω̇k (0) = ω̇k (0) (7.8) For fermions we have instead: V −1 dN (k, t)/dk = Σs w(ω + λπ)(2ω)−1 [ω 2 |f + |2 + |f˙+ |2 − iω(f ∗+ ∂0 f + − f + ∂0 f ∗+ )]. Similar expressions exist for the boost invariant problem [28]. 8. TRANSPORT APPROACH TO MULTIPARTICLE PRODUCTION A classical kinetic theory approach to the back-reaction problem as discussed in [14][15][16]introduces a phase space single particle distribution function f (x, p, t) in the presence of a homogeneous electric field and with a phenomenological source term inspired by Schwinger’s solution for the constant external field. df /dt = ∂f /∂t + eE(t)∂f /∂p = dN/dtdzdp = |eE(t)| ln[1 ± exp[−πm2⊥ /|eE(t)|]]δ(p) (8.1) ± stand for boson(+) or fermion case (-). The right hand side of (8.1) is a naive use of Schwinger’s formula (valid when no particles are present and for constant fields with E replaced by E(t) {or E(τ )}. This approach was recently used to predict dilepton production from the quark-gluon plasma [29]. A potential problem with replacing constant E by E(t) is that in the field theory simulations E(t) is rapidly varying in time. A more serious problem is that once particles are produced, Schwinger’s derivation, which was for particle production from the vacuum, is no longer valid. This however can be fixed up by the following arguement. Once particles are present there is an additional quantum mechanical effect due to statistics– Bose enhancement or Pauli Suppression. For the external field 19 problem one always has a normal mode decomposition at each time t. Thus the creation and annihilation operators at different t are again connected by a unitary transformation: b(k, t + ∆t) = α(t + ∆t)b(k, t) + β(t + ∆t)d† (k, t) |α|2 + |β|2 = 1; |b† b| = n+ ; |a† a| = n− ; n+ = n− = n (8.2) n(t + ∆t) = n(t) + 2|β(t + ∆t)|2 {1 ± n} (8.2) ∆n/∆t = 2|β|2 {1 ± n}/∆t (8.3) Therefore or where the +(-) is Bose enhancement (Pauli suppression). The Pauli suppression ensures n(k) ≤ 1 for fermions. Thus to include this effect we will modify the right hand side of (8.1) by multiplying by (1± 2 f (p,t)). This modified transport eqaution, as we will show below gives much better agreement with the field theory calculation. One can solve the Viasov equation using the method of characteristics. From dp/dt = eE and f(p,0)=0 one obtains: f (p, t) = Σi ln[1 ± exp[−πm2 /|eE(ti )|]] (8.4) where the ti are determined from p + eA(t) + eA(ti ) = 0; ti < t (8.5) The back reaction equation is now d2 A/dt2 = jcond + jpol Z jcond = 2e [dp]pf (p, t)/(p2 + m2 )1/2 Z jpol = 2/E [dp](p2 + m2 )1/2 d3 N/dtdxdp (8.7) where dN/dtdxdp = (1 ± 2f (p, t))|eE(t)| ln[1 ± exp[−πm2⊥ /|eE(t)|]]δ(p) A similar expression holds in boost invariant dynamics as discussed in [28]. The transport approach with the enhancement (suppression) factor gives reasonable agreement with the direct numerical solution of the field theory (in lowest order in 1/N) as long as we coarse grain the field theory result in momentum bins. 20 9. HYDRODYNAMIC CONSIDERATIONS: ENERGY FLOW From a hydrodynamical point of view, flat rapidity distributions seen in multiparticle production in p-p as well as A-p and A-A collisions are a result of the hydrodynamics being in a scaling regime for the longitudinal flow. That is for v=z/t (no size scale in the longitudinal dimension) the light cone variables τ, η: z = τ sinh η; t = τ cosh η (9.1) become the fluid proper time τ = t(1 − v 2 )1/2 and fluid rapidity: η = 1/2 ln[(t − x)/(t + x)] ⇒ 1/2 ln[(1 − v)/(1 + v)] = α (9.2) In the rest frame (comoving frame) of a perfect relativistic fluid the stress tensor has the form: Tµν = diagonal (ǫ, p, p, p) (9.3) Boosting by the relativistic fluid velocity four vector uµ (x, t) one has: Tµν = (ǫ + p)uµ uν − pg µν (9.4) Letting u0 = cosh α; u3 = sinh α, we have when v = zt that η = α, the fluid rapidity. If one has an effective equation of state p = p(ǫ) (which happens if both p and ǫ are functions of the single variable (τ ) as well as for the case of local thermal equilibrium) then one can formally define temperature and pressure as follows: Z ǫ + p = T s; dǫ = T ds; lns = dǫ/(ǫ + p) (9.5) Then the equation: uµ ∂ ν Tµν = 0 becomes: ∂ ν (s(τ )uν ) = 0 (9.6) ds/dτ + s/τ = 0 or sτ = constant (9.7) Which in 1 + 1 dimensions becomes The assumption of hydrodynamical models is that the initial energy density for the flow can be related to the center of mass energy and a given volume (say of a Lorentz 21 contracted disk of matter). It is also assumed that the flow of energy is unaffected by the hadronization process and that the fluid rapidity can be identified in the out regime with particle rapidity. Thus after hadronization the number of pions found in a bin of fluid rapidity can be obtained from the energy in a bin of rapidity by dividing by the energy of a single pion having that rapidty. That is one assumes that when the comoving energy density become of the order of one pion/(compton wave length)[23][24]we are in the out regime. This determines a surface defined by ǫ(τf ) = mπ /Vπ (9.8) On that surface of constant τ 0 dN/dη = 1/(mπ u )dE/dη = 1/(mπ cosh α) Z T 0µ dσµ /dη dσµ = 4πa2 (dz, −dt) = 4πa2 τf (cosh η, − sinh η) dN/dη = 4πa2 /mπ [(ǫ + p) cosh α cosh(η − α)p cosh η]/ cosh α = 4πa2 ǫ(τf )/mπ (9.9) which shows that when η = α one gets a flat distribution in fluid rapidity. An extra assumption is needed to identify fluid rapidity α with particle rapidity y = 1/2 ln[(Eπ + pπ )/(Eπ − pπ )], where pπ is the longitudinal momentum of the pion. What I would like to show next is that in a field theory calculation based on the Schwinger mechanism if we make the kinematical assumption that the electric field E is just a function of τ we obtain a flat rapidity distribution. We can also prove that the distribution in fluid rapidity is the same as the distribution in particle rapidity. We will also obtain renormalized expressions for ǫ(τ ) and p(τ ) (non-equilibrium dynamical equation of state). In the pair production problem we have shown that the interpolating phase space number density is given by the Bogoliubov function (7.7) : dN/dηdkη dk⊥ dx⊥ = |β(kη , k⊥ , τ )|2 (9.10) we are interested in transforming from dηdkη to dz dy where y is the particle rapidity y = 1/2 ln[(Eπ + kzπ )/(Eπ − kzπ )]. One can show that the transformation from (η, τ ) to (z,t) is a canonical one (in the sense of Poisson brackets {η, kη } = {τ, Ω} = 1) with canonical momentum kη = −Ez + tp = −τ m⊥ sinh(η − y) Ω = (Et − pz)/τ ) = m⊥ cosh(η − y) 22 (9.11) The phase space is unchanged by this change of variable thus d6 N/(dηdkη dk⊥ dx⊥ ) = d6 N/dzdkz dk⊥ dx⊥ = JdN 6 /dzdydk⊥ dx⊥ (9.12) where J−1 = ∂kη /∂y∂η/∂z. At fixed τ one can show that |J| = dz/dkη which leads to desired result, assumed by Landau that along the breakup surface τ = constant: dN/dy = dN/dη. (9.13) Schwinger’s pair production mechanism leads to an Energy Momentum tensor which is diagonal in the(τ, η, x⊥) coordinate system which is thus a comoving one. In that system one has: T µν = diagonal {ǫ(τ ), pk (τ ), p⊥ (τ ), p⊥ (τ )} (9.14) We see in a 3 dimensional problem, the field theory in this approximation has two separate pressures, one in the longitudinal direction and one in the transverse direction and thus differs from the thermal equilibruim case. However, for a one-dimensional flow we have that the energy in a bin of fluid rapidity is just: dE/dη = Z T 0µ dσµ = A⊥ τ cosh ηǫ(τ ) (9.15) which is just the (1 + 1) dimensional hydrodynamical result of (9.9). This result does not depend on any assumptions of thermalization. In the field theory calculation the expectation value of the stress tensor must be renormalized since the electric field undergoes charge renormalization. We can also determine the two pressures and the energy density as a function of τ . Explicitly we have in the fermion case. ǫ(τ ) =< Tτ τ >= τ Σs Z 2 [dk]Rτ τ (k) + ER /2 (9.16) where Rτ τ (k) = 2(p2⊥ + m2 )(g0+|f + |2 − g0− |f − |2 ) − ω − (p2⊥ + m2 )(π + eȦ)2 /(8ω 5 τ 2 ) Z 1 2 2 2 pk (τ )τ =< Tηη >= τ Σs [dk]λs πRηη (k) − ER τ (9.17) 2 23 where Rηη (k) = 2|f + |2 − (2ω)−1 (ω + λs π)−1 − λs eȦ/8ω 5 τ 2 − λs eĖ/8ω 5 − λs π/4ω 5 τ 2 + 5πλs (π + eȦ)2 /(16ω 7 τ 2 ) and p⊥ (τ ) =< Tyy >=< Txx >= (4τ ) −1 Σs Z 2 [dk]{p2⊥ (p2⊥ + m2 )−1 Rτ τ − 2λπp2⊥ Rηη } + ER /2. (9.18) Thus we are able to numerically determine the dynamical equation of state pi = pi (ǫ) as a function of τ . 10. DISCUSSION OF NUMERICAL RESULTS The physical quantities that we determine numerically are the time evolution of E(t), A(t) , and j(t). We will display in the figures the plasma oscillations and the time scale for field energy to be essentially transferred into pair production. The other quantities of physical interest we determine are the spectra of produced particles dN/dk and the dynamical equation of state. For comparison we have also solved the phenomenological transport theory with and without the quantum correction due to statistics (i.e. Pauli Blocking and Bose Enhancement). In making plots for the spatially homogeneous case we use the dimensionless variables [5]: Ẽ = (eE/m2 ) Ã = eA/m; mt = τ When the Electric field Ẽ is > 1 then it is quite easy for pairs to be produced and in that regime the final result is independent of the initial data. We can see the approach to the tunneling regime by comparing in the regime. 5 < Ẽ0 < 2 the behavior of E(τ ). This is shown in fig. 3 for Ẽ0 = .5, 1, 2. Once Ẽ0 > 2 the behavior of Ẽ(τ ) is only weakly dependent on Ẽ0 . Once the pairs are produced one sees that there are plasma oscillations superimposed on which the electric field degrades. These figures are from early simulations for scalar QED in 1 + 1 dimensions [5]. In fig. 4 we show Ã(t), ẽ(t), < j(t) > for Ẽ0 = 2 for scalar QED in 1+1 dimensions [6]. In fig. 5 we show Ẽ and j̃ for scalar QED in 1+1 dimensions for Ẽ0 = 4. We compare the naive Vlasov approach (dashed line) and the improved Vlasov approach ( dot- dashed line) . We notice that including Bose enhancement corrections is quite important. We also 24 notice that j̃ max = 2 e ρc so that particles continue to be produced when Ẽ is near a maximum. In fig. 6 We show the exact particle spectrum as well as the momentum space smoothed result which is compared to the Vlasov Equation. Here Ẽ0 =1 and we have scalar QED in 1+1 dimensions. In fig. 7 we show the results for E and j for Ẽ0 = 4 in QED in 1+1 dimension compared to the uncorrected Vlasov equation. We notice the dismal agreement. In fig. 8 we see the same curves compared to the improved Vlasov equation. In fig. 9 we show the exact spectrum of produced pairs for QED in 1+1 dimensions for Ẽ0 = 4. We notice that n(k) ≤1 which expresses the Pauli Principle. In fig. 10 we compare the binned version of the field theory result with both the Naive and Improved transport theory. Next we present recent results for Scalar QED in 1+1 dimensions using boost invariant Kinematics. In fig. 11 we plot E,A and j vs. u=log(τ ) for E0 (u0 ) = 4 in the boost invariant case where E is a function of the proper time τ (not to be confused with the previous τ ). In fig. 12 we compare E(u) and j(u) with the boost invariant transport theory with and without Bose enhancement. Finally we present preliminary results[30] for scalar QED in 3+1 dimensions. In fig. 13 we show the time evolution for E(t) and j(t) in 3+1 dimensions and compare with the Boltzmann-Vlasov model with and without Bose-enhancement. 11. ACKNOWLEDGEMENTS The work presented here was done in collaboration with Emil Mottola, Yuval Kluger, So-Young Pi, Ben Svetitsky, Judah Eisenberg, Paul Anderson, Barrett Rogers and M. Samiullah. 25 References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] F. Cooper, S. Y. Pi and P. Stancioff, Phys. Rev. D. 34, 3831 (1986) S. Y. Pi and M. Samiullah, Phys. Rev. D. 36 3128 (1987) Fred Cooper and Emil Mottola, D. 36, 3114(1987) Fred Cooper and Emil Mottola, Phys. Rev. D 40, 456 (1989) F. Cooper, E. Mottola, B. Rogers, and P. Anderson, in intermittency in High Energy Collisions, editied by F. Cooper, R.C. Hwa and I. Sarcevic (World Scientific, Singapore, 1991) p. 399 Y. Kluger, J. M. Eisenberg, B. Svetitsky, F. Cooper, and E. Mottola, Phys. Rev. Lett. 67, 2427(1991) Y. Kluger, J. M. Eisenberg, B. Svetitsky, F. Cooper, and E. Mottola, Phys. Rev. D 45, 4659 (1992) J. Schwinger, J. Math Phys. (N.Y.) 2, 407 (1961) L.V. Keldysh, Sov. Phys. JETP 20 1018 (1965) K. Chou, A. Su, B. Hao, L. Yu, Phys. Rpts. 118 (1985) R. D. Jordon, Phys. Rev. D. 33, 444 (1986) L. Bialynicki-Birula, P. Gornicki, and J. Rafelski, Phys. Rev. D 44, 1825 (1991) L. Parker, and S. A. Fulling, Phys. Rev. D. 7 2357 (1973) K. Kajantie and T. Matsui, Phys. Lett. 164B, 373 (1985) G. Gatoff, A. K. Kerman, and T. Matsui, Phys. Rev. D. 36, 114 (1987) A. Bialas, W. Czyz, A. Dyrek, and W. Florkowski, Nucl. Phys. B296, 611 (1988) P. A. M. Dirac, Proc. Camb. Phil. Soc. 26 (1930) P. Anderson, F. Cooper, S. Habib, E. Mottola, and J. Paz [in preparation] S. Nussinov, Phys. Rev. Lett. 34, 1296 (1975) J. Schwinger, Phys. Rev 82, 664 (1951) C. Itzykson and J. Zuber Quantum Field Theory. McGraw-Hill (1980) A. Casher, H. Neuberger, S. Nussinov, Phys Rev. D. 28, 179 (1979) L. D. Landau, Izv. Akad. Nauyk SSSR 17 (1953) 51 F. Cooper, G. Frye and E. Schonberg, Phys. Rev. D 11(1975) 192 J. D. Bjorken, Phys. Rev. D 27 (1983), 140. N. D. Birell and P.C. W. Davies Quantum Fields in Curved Space (Cambridge University press, Cambridge, England, 1982) S. Weinberg, Gravitation and Cosmology: Principles and Applications of the General Theory of Relativity(Wiley, New York, 1972) F. Cooper, J. M. Eisenberg, Y. Kluger, E. Mottola and B. Svetitsky “Particle Production in the Central Rapidity Region” Tel-Aviv Preprint TAUP 1944-92. M. Asakawa and T. Matsui, Phys. Rev. D43, 2871 (1991) Y. Kluger, PhD. Thesis. Tel-Aviv University (June 1992). (Unpublished) 26