Radoslaw Ryblewski
email: radoslaw.ryblewski@ifj.edu.pl

Wojciech Florkowski
email: wojciech.florkowski@ifj.edu.pl

Non-boost-invariant motion of dissipative and highly anisotropic fluid

Keywords: numbers: 25.75.-q, 25.75.Ld, 24.10.Nz, 25.75.Nq

The recently formulated framework of anisotropic and dissipative hydrodynamics (ADHYDRO) is used to describe non-boost-invariant motion of the fluid created at the early stages of heavy-ion collisions. Very strong initial asymmetries of pressure are reduced by the entropy production processes. By the appropriate choice of the form of the entropy source we can describe isotropization times of about 1 fm, which agrees with the common expectations that already at such times the perfect-fluid hydrodynamics may be applied. Our previous results are generalized by including the realistic equation of state as the limit of the isotropization processes.

Introduction

In a recent paper we have introduced the new framework of highly-Anisotropic and strongly-Dissipative HYDROdynamics (ADHYDRO) [1]. This framework allows for studies of the systems which are locally highly anisotropic. Such situations happen at the early stages of relativistic heavy-ion collisions where, as suggested by many microscopic calculations [2,3,4], the longitudinal pressure is much lower than the transverse pressure ‡.

In Ref. [1] our approach has been used to analyze the behavior of matter created in the central region of relativistic heavy-ion collisions, where the assumption of boostinvariance is acceptable. We have shown how our equations describe the isotropization of pressure and how the system approaches the regime governed by the perfect-fluid hydrodynamics.

In this paper we relax the assumption of boost-invariance and consider a general one-dimensional expansion along the beam axis (at the early stages of collisions the longitudinal motion dominates). We derive equations describing non-boost-invariant motion of anisotropic fluid and show how the initial anisotropic behavior transforms smoothly into the locally isotropic expansion.

In the second part of the paper we present a scheme implying that the isotropization transition leads to a desired equation of state in the perfect-fluid stage. In this way we construct an effective model which describes several stages of heavy-ion collisions in a very concise way: an early anisotropic phase (Glasma, color-flux tubes), transition to the isotropic phase with the QCD equation of state, perfect-fluid stage (sQGP), and the phase transition to the hadron gas.

We note that similar physical problems connected with the non-boost-invariant motion of the fluid formed at the early stages of relativistic heavy-ion collisions were studied in Refs. [5,6], see also [10,11,12]. In Ref. [5] the effects of the equation of state and the initial conditions were studied, while in Ref. [6] the effects of the shear viscosity were taken additionally into account. In contrast to Refs. [5,6] our approach concentrates on highly asymmetric stages of the evolution, see [7,8]. We also note that our formulation of ADHYDRO was followed by the paper by Martinez and Strickland [9], where a similar framework of anisotropic hydrodynamics was introduced. We discuss the connection between our work and Ref. [9] in the Appendix.

Anisotropic hydrodynamics

In this Section we recapitulate the main assumptions of ADHYDRO, that have been recently introduced in [1]. Our starting point is the following form of the energymomentum tensor,

T µν = (ε + P ⊥ ) U µ U ν -P ⊥ g µν -(P ⊥ -P )V µ V ν . (1)
Here ε, P ⊥ , and P are the energy density, transverse pressure, and longitudinal pressure, respectively. For isotropic fluid, where the two pressures are equal, P ⊥ = P = P , we recover the form of the energy-momentum tensor of the perfectfluid hydrodynamics. The four-vector U µ in (1) represents the hydrodynamic flow,

U µ = γ(1, v x , v y , v z ), γ = (1 -v 2 ) -1/2 , (2)
‡ The directions are defined with respect to the beam axis while V µ defines the direction of the longitudinal axis that plays a special role due to the initial geometry of the collision,

V µ = γ z (v z , 0, 0, 1), γ z = (1 -v 2 z ) -1/2 . (3)
The four-vectors U µ and V µ satisfy the following normalization conditions:

U 2 = 1, V 2 = -1, U • V = 0. (4)
In the local-rest-frame (LRF) of the fluid element we have U µ = (1, 0, 0, 0) and V µ = (0, 0, 0, 1), and the energy-momentum tensor has a simple diagonal structure §,

T µν =     ε 0 0 0 0 P ⊥ 0 0 0 0 P ⊥ 0 0 0 0 P     . (5)
Hence, as expected, the formula (1) allows for different pressures in the longitudinal and transverse directions. Anisotropies of the form (5) are present in the case where the particles are formed by the decays of color strings and in the theory of Color Glass Condensate [13,14,15] . Besides the energy-momentum tensor (1), we introduce the entropy flux

σ µ = σU µ , ( 6
)
where σ is the entropy density. We assume that ε and σ are functions of P ⊥ and P .

In particular, for massless partons the condition T µ µ = 0 gives ε = 2P ⊥ + P .

The form of the energy-momentum tensor (1) resembles the form used in relativistic magnetohydrodynamics [16,17]. In that case the anisotropy is induced by the presence of the magnetic field. At the early stages of heavy-ion collisions we have similar situation -there exist strong color magnetic and electric longitudinal fields (Glasma following CGC [18]) which polarize the medium. In this paper, however, we do not include explicitly the effects of the fields. Such analysis is left for a separate study.

The space-time evolution of the system is governed by the equations expressing the energy-momentum conservation and the entropy growth,

∂ µ T µν = 0, (8)
∂ µ σ µ = Σ. ( 9
)
The function Σ represents the entropy source. The form of Σ must be treated as an assumption that defines the dynamics of the anisotropic fluid. It is natural to assume that Σ ≥ 0 and Σ = 0 for P ⊥ = P . In this way, in the case where the two pressures are equal, the structure of the perfect-fluid hydrodynamics is recovered. We treat Σ as a function of P ⊥ and P . In this way, Eqs. ( 8) and ( 9) form a closed system of 5 equations for 5 unknown functions: three components of the fluid velocity, P ⊥ , and P . The projections of Eq. ( 8) on U ν and V ν yield § The Lorentz transformation which leads from the center-of-mass frame, where Eqs. (2) and (3) are valid, to LRF, where we have U µ = (1, 0, 0, 0) and V µ = (0, 0, 0, 1), consists of a Lorentz boost along the z-axis, a rotation around the z-axis, and a boost in the direction of the remaining transverse part of v. In this way the initial longitudinal direction defined by the spatial part of V µ is left unchanged.

U µ ∂ µ ε = -(ε + P ⊥ ) ∂ µ U µ + P ⊥ -P U ν V µ ∂ µ V ν , (10)
V µ ∂ µ P = -P -P ⊥ ∂ µ V µ + (ε + P ⊥ ) V ν U µ ∂ µ U ν . (11)

Anisotropy parameter x

Our earlier studies show [1,19,20] that it is useful to switch from P ⊥ and P to the two new variables: the entropy density σ and the anisotropy parameter x. In this case we may write

ε = π 2 σ 4g 0 4/3 R(x), (12)
P ⊥ = π 2 σ 4g 0 4/3 R(x) 3 + xR (x) , (13)
P = π 2 σ 4g 0 4/3 R(x) 3 -2xR (x) . (14)
The function R(x) is defined by the formula [20]

R(x) = 3 g 0 x -1 3 2π 2 1 + x arctan √ x -1 √ x -1 , (15)
and g 0 is the degeneracy factor connected with internal quantum numbers of particles that form the fluid. In our calculations we take g 0 = 16 -the initial system is formed most likely from gluons. The symbol R (x) denotes the derivative of R(x) with respect to x, for x = 1 we have R (1) = 0 and the two pressures are equal. The physical interpretation of x and details describing manipulations that lead to ( 12)-( 15) are given in [1,20]. Here we only note that the ratio P ⊥ /P is a monotonic function of x with P ⊥ > P if x > 1 (to a good approximation P ⊥ /P = x 3/4 , and P ⊥ = P if x = 1).

One-dimensional non-boost-invariant motion

For a purely one dimensional motion we may introduce the following parameterizations:

U µ = (cosh ϑ(τ, η), 0, 0, sinh ϑ(τ, η)) (16)
and

V µ = (sinh ϑ(τ, η), 0, 0, cosh ϑ(τ, η)) , (17)
where ϑ(τ, η) is the fluid rapidity depending on the (longitudinal) proper time

τ = t 2 -z 2 (18)
and the space-time rapidity

η = 1 2 ln t + z t -z . (19)
Equations ( 16) and ( 17) satisfy automatically the normalization conditions (4). We also have

U µ ∂ µ = cosh(ϑ -η) ∂ ∂τ + sinh(ϑ -η) τ ∂ ∂η . (20)
The expressions similar to (20) may be derived for ∂ µ U µ , V µ ∂ µ , and ∂ µ V µ . Using them in Eqs. ( 10) and (11) gives

∂ ∂τ + tanh(ϑ -η) τ ∂ ∂η ε = -(ε + P ) tanh(ϑ -η) ∂ ∂τ + ∂ τ∂η ϑ, (21) tanh
(ϑ -η) ∂ ∂τ + ∂ τ∂η P = -(ε + P ) ∂ ∂τ + tanh(ϑ -η) τ ∂ ∂η ϑ. (22)
Similarly, from Eq. ( 9) one obtains

∂ ∂τ + tanh(ϑ -η) τ ∂ ∂η σ + σ tanh(ϑ -η) ∂ ∂τ + ∂ τ∂η ϑ = Σ cosh(ϑ -η) . ( 23
)
The structure of Eqs. ( 21)-( 23) is very much similar to the equations studied in Refs. [5].

Entropy source

In the boost-invariant case the functions ε, P , P ⊥ , and σ do not depend on η, and ϑ = η. In this case Eq. ( 22) is automatically fulfilled, while Eqs. ( 21) and ( 23) are reduced to

dε dτ = - (ε + P ) τ , (24)
dσ dτ + σ τ = Σ. ( 25
)
Equations ( 24) and ( 25) were studied in [1] where the following form of the entropy source Σ was used

Σ = (1 - √ x) 2 √ x σ τ eq . ( 26
)
Here τ eq is a timescale parameter. The form (26) guarantees that Σ ≥ 0 and Σ(σ, x = 1) = 0. More arguments for such a particular form of Σ are given in [1]. They connect x with the ratio of the transverse and longitudinal temperatures.

It may be shown that for small deviations from the equilibrium Eq. ( 26) leads to quadratic dependence of the entropy source on the variable ξ = 1 -x and this dependence is compatible with the Israel-Stewart theory (where the entropy production depends on the viscous stress squared) and with the Martinez-Strickland model [9].

We have to stress, however, that the structure of the entropy source is an external input for the anisotropic hydrodynamics. It would be interesting to obtain any hints about Σ for large x from the microscopic models of particle production or from the AdS/CFT correspondence.

Isotropization

If ε and P are expressed in terms of σ and x, Eqs. ( 21)-( 23) become three equations for three unknown functions: σ, x, and ϑ:

4 3σ ∂ ∂τ + tanh(ϑ -η) τ ∂ ∂η σ = - R (x) R(x) ∂ ∂τ + tanh(ϑ -η) τ ∂ ∂η x -(1 + h(x)) tanh(ϑ -η) ∂ ∂τ + ∂ τ∂η ϑ, (27)
4 3σ tanh(ϑ -η) ∂ ∂τ + ∂ τ∂η σ = - R (x) R(x) + h (x) h(x) tanh(ϑ -η) ∂ ∂τ + ∂ τ∂η x - 1 + h(x) h(x) ∂ ∂τ + tanh(ϑ -η) τ ∂ ∂η ϑ, (28)
∂ ∂τ + tanh(ϑ -η) τ ∂ ∂η σ + σ tanh(ϑ -η) ∂ ∂τ + ∂ τ∂η ϑ = Σ cosh(ϑ -η) . ( 29
)
Here we have introduced the ratio of the longitudinal pressure and energy density, which equals 1/3 for the isotropic system with x = 1,

h(x) = P ε , h(1) = 1 3 . (30)
Substituting ( 29) into (27) and using the relation

1 + h(x) = 4 3 -2x R (x) R(x) , (31)
we find

∂ ∂τ + tanh(ϑ -η) τ ∂ ∂η x = 2x tanh(ϑ -η) ∂ ∂τ + ∂ τ∂η ϑ - 4H(x) 3τ eq cosh(ϑ -η) .( 32
)
Here we use the notation introduced in [1],

H(x) = R(x) R (x) (1 - √ x) 2 √ x . (33)
In the case where ϑ(τ, η) ≈ η, Eq. (32) may be approximated by the ordinary differential equation derived in [1]

dx dτ = 2x τ - 4H(x) 3τ eq . (34)
Thus, in this case we may immediately use our previous results to argue that x → 1 for τ τ eq . If the condition ϑ(τ, η) ≈ η is not fulfilled, one should study Eqs. ( 27)-(29) numerically.

Initial conditions

The initial conditions for a non-boost-invariant evolution are defined by three functions: σ(τ 0 , η), ϑ(τ 0 , η), and x(τ 0 , η), where τ 0 is the initial proper time. In the numerical calculations we assume τ 0 = 0.2 fm. The initial entropy profile is taken in the form [21,22]

σ(τ 0 , η) = σ 0 exp - (|η| -∆η) 2 2(δη) 2 θ(|η| -∆η) , ( 35
)
where θ is the step function, the parameter ∆η defines the half width of the initial plateau in spacetime rapidity, and δη defines the half width of the Gaussian tails on both sides of the plateau. By the appropriate changes of the parameters ∆η and δη we may vary between the boost-invariant-like and Gaussian-like initial conditions. To match the rapidity distribution measured by BRAHMS [23] we use the values: ∆η = 1 and δη = 1.3. The value of the initial central entropy density is obtained from the condition that the initial energy density is 100 GeV/fm 3 ,

ε 0 = 100 GeV/fm 3 = π 2 σ 0 4g 0 4/3 R(x 0 ). ( 36
)
Here, x 0 is the initial value of the anisotropy parameter at η = 0. In the present calculations we assume that the initial profile of x is constant,

x(τ 0 , η) = x 0 = 100. ( 37
)
This value of x 0 is somewhat arbitrary but it expresses the opinion that the initial longitudinal pressure is much smaller than the initial transverse pressure.

The value of the timescale parameter τ eq has been set equal to 0.25 fm. Finally, in agreement with other hydrodynamic calculations which address the problem of the longitudinal expansion, we take the initial fluid rapidity profile in the form

ϑ(τ 0 , η) = η. ( 38
)
Of course, for the boost-invariant motion we have always ϑ = η. In our case, due to the presence of the longitudinal gradients determined by the entropy profile (35), we expect |ϑ| ≥ |η|.

Results

The results of our numerical calculations with the initial conditions specified in the previous Section are presented in Figs. 1 and2. The upper part of Fig. 1 shows the time evolution of the anisotropy parameter x for two different values of the space-time rapidity: η = 0 (solid line) and η = 6 (dashed line). We observe the initial rapid decrease of x (from 100 to about 5) followed by a longer evolution where x approaches unity. This behavior is very weakly dependent on rapidity and very much similar to that observed previously in the boost-invariant case. The lower part of Fig. 1 shows the time evolution of the ratio P /P ⊥ . The time dependence of this ratio reflects the time dependence of the anisotropy parameter x. Again, we observe a rapid change at the very beginning of the evolution that is followed by a more moderate approach to unity. For τ > 2 fm the system is practically isotropic. Clearly, by changing the value of τ eq we may speed up or slow down the process of isotropization. The most common perfect-fluid description of the data requires that the stage described by perfect hydrodynamics starts at about 1 fm or even earlier. Therefore, we think that the time evolution obtained in our approach may be regarded as realistic if the systems becomes isotropic at about 1 fm. The results shown in Figs. 1 and2 suggest that this condition may be achieved if τ eq ≤ 0.1 fm. Such extremely short time scales indicate very strong interactions in the early anisotropic plasma. The upper part of Fig. 2 shows the difference between the fluid rapidity ϑ and the space-time rapidity η for τ = τ 0 (solid line), τ = 2τ 0 (dashed line), and τ = 10 fm (dotted line). In the boost-invariant case ϑ = η. In the present calculation, due to the extra pressure gradients determined by the entropy profile, the values of the fluid rapidity become larger than the space-time rapidity (for η > 0).

The lower part of Fig. 2 shows the entropy density profiles in rapidity for different values of the proper time. The solid line shows the initial entropy profile (35) with the central entropy density σ 0 determined by the condition (36). The thick dashed line shows the entropy density profile for the proper time τ = 1.2τ 0 . Since the initial dynamics is dissipative, the entropy is produced in the early stage, and the thick dashed line is placed above the solid line. For comparison, the thin dashed line shows the entropy profile expected from the boost-invariant isentropic expansion where σ(τ, η) = σ(τ 0 , η)τ 0 /τ . The thin dashed line is placed below the solid line, since the entropy density decreases in this case due to the longitudinal expansion of the fluid. By comparing the maximal values of the thick and dashed lines we can make the estimate of the relative entropy production due to the dissipative effects present in our approach. From this simple calculation we find that the entropy increases by about 70%.

Transition to quark-gluon plasma

The concept of ADHYDRO is motivated by the microscopic picture of the anisotropic parton distribution function [1,20]. As the system becomes isotropic, the parton distribution becomes also isotropic (locally, in the momentum space). This, however, does not automatically mean that the isotropic system is described by the realistic equation of state. The Boltzmann distribution used in Refs. [1,20] may be inadequate for description of the strongly interacting quark-gluon plasma. To connect the process of isotropization with the process of formation of the equilibrated quark-gluon plasma we may consider the following ansatz for the energy density and pressure

ε(σ, x) = ε qgp (σ) 3

π 2 g 0 R(x), (39)
P ⊥ (σ, x) = P qgp (σ) π 2 g 0 R(x) 3 + xR (x) . (40)
P (σ, x) = P qgp (σ) π 2 g 0 R(x) 3 -2xR (x) . (41)
Here, the functions ε qgp (σ)/3 and P qgp (σ) describe the realistic equation of state constructed in [24]. They appear in Eqs. (39)-(41) in the very much analogous way as the term σ 4/3 in Eqs. ( 12)-( 14). Although ε qgp (σ) and P qgp (σ) were calculated for equilibrium, they are used by us to calculate the energy density and pressure for an arbitrary state characterized by the entropy σ and the anisotropy parameter x.

In order to take into account the effects connected with the anisotropy, we multiply ε qgp (σ)/3 and P qgp (σ) by the same combinations of the function R(x) as in Eqs. ( 12)-( 14). The factor π 2 /g 0 has been introduced to guarantee that ε(σ, x = 1) = ε qgp (σ), and similarly P (σ, x = 1) = P ⊥ (σ, x = 1) = P qgp (σ). We should admit that Eqs. (39)-(41) have no direct microscopic explanation. The main argument for using this form is that it has two attractive limits and may be used as an interpolating formula. At the very early stage, the system consists most likely of massless partons and Eqs. (39)-(41) approach Eqs. ( 12)-( 14), since the realistic equation of state approaches the Stefan-Boltzmann limit (although deviations from this limit are noticeable). On the other hand, if the system becomes isotropic, x → 1, R (x) → 0, and Eqs. (39)-(41) are reduced to the equation of state used in standard hydrodynamics. From the statistical point of view, the initial system of partons is always treated classically and R(x) is given by (15). Inclusion of quantum statistics and consequences of other model assumptions will be discussed in a separate paper.

In Figs. 3 and4 we show the results of the numerical calculations where Eqs. (39)-(41) were substituted in Eqs. ( 21)- (23). By this substitution we have obtained the system of three equations for σ, x, and ϑ, which is quite analogous to Eqs. ( 27)-( 29). The values of the initial parameters and the initial profiles have not been changed except for the initial value of the entropy density, which has been modified to keep the initial energy density unchanged,

ε 0 = 100 GeV/fm 3 = ε qgp (σ 0 ) π 2 g 0 R(x 0 ) 3 . ( 42
)
The overall space-time evolution of the system shown in Figs. 3 and4 resembles very much the evolution shown earlier in Figs. 1 and2. In particular, the time changes of x and the ratio of pressures are practically the same. The small differences are visible in the evolution of the entropy density (mainly due to a different normalization caused by fixing the initial energy density rather than the entropy density) and the fluid rapidity. The difference ϑ -η becomes now a more complicated function of η, the effect noticed also in Ref. [6] and connected with the use of the equation of state from Ref. [24].

Clearly, an attractive feature of using Eqs. ( 39)-(41) instead of Eqs. ( 12)-( 14) is that the system reaching equilibrium is described by the realistic QCD equation of state. In this way, the intermediate hydrodynamic evolution agrees with our expectations and no extra construction (like, e.g., the Landau matching condition used in [25]) is necessary to change from one thermodynamic regime to another one. On the other hand, the connection of Eqs. ( 39)-(41) to any microscopic picture is lacking, since the realistic QCD equation of state cannot be simply interpreted in a quasi-particle picture. Another attractive feature Eqs. ( 39)-( 41) is that they include the phase transition from the quark-gluon plasma to the hadron gas. Thus, using Eqs. (39)-(41) we may describe many stages of heavy-ion collisions in a single approach.

Conclusions

Our results confirm very attractive features of ADHYDRO as an effective model for early stages of heavy-ion collisions. In this paper we have shown that the tendency to approach isotropy is not changed by non-boost-invariant profiles. Very strong initial asymmetries of pressure are reduced by the entropy production processes. By the appropriate choice of the form of the entropy source and the value of the timescale parameter, we have succeeded to reproduce the realistic times for the isotropization expected in heavy-ion collisions. We have also generalized our previous result by including the realistic equation of state as the limit of the isotropization processes.

From the practical point of view, our formalism based on the anisotropic distribution function allows for the determination of the space-time evolution of colorneutral anisotropic distributions, which may be used, for example, as background distributions in the analysis of the plasma instabilities.

The results obtained within ADHYDRO in this paper, and also in our previous publication, open perspectives for further applications of the model. In particular, the model may be applied in the 2+1 case (longitudinal boost-invariance without cylindrical symmetry) or in the general 3+1 case. In the 2+1 case it may be treated as a generalization of the framework where transverse hydrodynamics is combined with the perfect-fluid hydrodynamics by the Landau matching conditions [25]. In ADHYDRO the change from the system where only transverse degrees of freedom are thermalized to the full local equilibrium is described in the continuous way.

Appendix

Our formulation of ADHYDRO was followed by the paper by Martinez and Strickland [9] where similar ideas appeared. It is interesting to show that the two approaches agree for small deviations from equilibrium in the boost-invariant case. In the following, it is useful to use the variable

ξ = x -1. (43)
Equation ( 20) from [9] connects ξ and p hard , Ref. [9]. For small ξ they are both well approximated by the linear form (45/16)ξ.

1 1 + ξ dξ dτ - 2 τ - 6 p hard dp hard dτ = 2Γ 1 -R 3/4 (ξ) 1 + ξ . (44)
Here, p hard is related to the average momentum in the parton distribution function and R(ξ) is defined by the formula

R(ξ) = 1 2 1 1 + ξ + arctan √ ξ √ ξ . (45)
For |ξ| < 1 we expand R(ξ) around zero. Keeping the leading terms in ξ one obtains from (44)

1 1 + ξ dξ dτ - 2 τ - 6 p hard dp hard dτ = - Γξ 2 15 . (46)
In order to connect this result with our approach, we relate the non-equilibrium entropy density σ to p hard [9],

σ = Ap 3 hard x -1/2 . ( 47
)
Here A is an irrelevant constant. Inserting Eq. (47) into Eq. ( 25) and using our ansatz for Σ we find

-
The consistency condition (49) may be further used in Eq. (24a) of Ref. [9] to show that in the limit |ξ| < 1 it agrees with our fundamental expression (34). To do so, we rewrite Eq. (24a) in the analogous form, namely

dx dτ = 2x τ - 4ΓH MS (x) 3 . (50)
In Fig. 5 we show our function H(ξ) and the function H MS (ξ) (ξ = x -1) multiplied by the factor 15/2. The two functions agree for small values of ξ, hence again the consistency is achieved between the two approaches if the condition (49) is applied.

Since the authors of Ref. [9] showed that their approach is reduced to the 2nd order viscous hydrodynamic equations of Israel and Stewart in the leading order in ξ, the comparison showed in this Appendix shows that our approach is also reduced to viscous hydrodynamics in the considered case. We note that H(ξ) and (15/2)H MS (ξ) differ for large values of ξ. In this region there are no direct hints about the behavior of matter from the underlying kinetic theory and the form of the function H(ξ) is not constrained. It may be postulated and eventually verified by making comparisons of the model results with the experimental data. Such analyses represent the most interesting applications of ADHYDRO.

6 Figure 1 .

 61 Figure 1. Upper part: Time evolution of the anisotropy parameter x for two different values of the space-time rapidity: η = 0 (solid line) and η = 6 (dashed line). Lower part: time evolution of the ratio P /P ⊥ for η = 0 (solid line) and η = 6 (dashed line).

Figure 2 .

 2 Figure 2. Upper part: difference between the fluid rapidity ϑ and the space-time rapidity η for τ = τ 0 (solid line), τ = 2τ 0 (dashed line), and τ = 10 fm (dotted line). Lower part: profiles of the entropy density plotted for different values of the proper time. The results of our numerical calculations are represented by the thick lines. The results obtained from the simple scaling σ(τ, η) = σ(τ 0 , η)τ 0 /τ are represented by the thin lines.

6 Figure 3 .

 63 Figure 3. The same as Fig. 1 but with Eqs. (12)-(14) replaced by Eqs. (39)-(41).

Figure 4 .

 4 Figure 4. The same as Fig. 2 but again with Eqs. (12)-(14) replaced by Eqs. (39)-(41).

x 1 Figure 5 .

 15 Figure5. Functions H(ξ) and (15/2)H MS (ξ) used in the present paper and in Ref.[9]. For small ξ they are both well approximated by the linear form (45/16)ξ.

Note that for x < 1 the function (arctan√ x -1)/ √ x -1 should be replaced by (arctanh √ 1 -x)/ √ 1 -x

Acknowledgments

This work was supported in part by the MNiSW grants No. N N202 288638 and N N202 263438.


Andreozzi

Lo Iudice
email: loiudice@na.infn.it

Porrino

A diagonalization algorithm revisited and applied to nuclear shell model

Keywords: numbers: 21, 60, Cs,23, 20, Lv

An importance sampling iterative algorithm for diagonalizing large matrices is upgraded and adopted for large scale nuclear shell model calculations using a spin uncoupled basis. Its numerical implementation shows that the iterative procedure converges rapidly to the exact eigensolutions achieving an effective drastic cut of the sizes of the Hamiltonian matrix.

Introduction

Lanczos algorithm [1,[START_REF] Golub | Matrix Computations[END_REF] is widely adopted in several branches of physics and other disciplines for diagonalizing large size matrices. Its many implementations in nuclear structure are devoted to the shell model (SM) eigenvalue problem.

Among them, some codes adopt a basis of states projected in spin (or spinisospin) (j -j-scheme) [START_REF] Caurier | [END_REF]4,[START_REF] Brown | NUSHELL@MSU MSUNSCL Report[END_REF][START_REF] Novoselsky | [END_REF]] so as to deal with Hamiltonian matrices of relatively small dimensions. These matrices, however, are very dense and require lengthy and cumbersome calculations of fractional parentage and recoupling coefficients.

Other adopt the unprojected m-scheme. In such a basis, the Hamiltonian matrix, though of much larger dimensions, has a far simpler structure and is sparse. Both properties are fully exploited in the powerful Antoine code [START_REF] Caurier | [END_REF]4].

Alternative to Lanczos is the shell model Monte Carlo (SMMC) [7], which avoids the diagonalization of the Hamiltonian matrix by adopting a Monte Carlo technique to evaluate expectation values and strength functions. The SMMC, however, is not suited for computing spectra and has to face the sign problem which has not found a clear-cut solution yet.

The MC technique is used in the quantum Monte Carlo diagonalization (QMCD) method [8,9] just to generate stochastically a truncated basis. This is then used to diagonalize the Hamiltonian. Specific projection methods had to be developed in order to restore the angular momentum and other symmetries, broken by the stochastic sampling. The QMCD was adopted for systematics, especially in the (p, f ) region.

Another technique for truncating the basis is provided by the density matrix renormalization group (DMRG) [10], borrowed from condensed matter [11], based on a partitioning of the full space into smaller subspaces. This method seems to be well suited to schematic Hamiltonians. When applied to shell model, it seems to work better in the coupled scheme [12].

Coming back to Lanczos, its numerical implementation deserves special care in order to avoid numerical errors, whose propagation, during the iterative process, may yield "ghost" states.

Few years ago [13,14], we proposed an iterative algorithm, alternative to Lanczos, which is free of instabilities and of easy implementation. It is also endowed with an importance sampling, which achieves a truncation of the basis. The algorithm was implemented in the coupled j-scheme only, which is too time consuming in the iterative construction of the Hamiltonian matrix. Here, we adopt the m-scheme and elaborate an efficient variant of the method which promises to be considerably more competitive. Its performance will be illustrated through some numerical examples.

The algorithm is as general as Lanczos and, therefore, is not bound to nuclear structure. It may be used in other fields of physics and outside of the realm of physics.

The algorithm

The algorithm has been described in Refs. [13,14]. We give here a brief outline. Let us consider a matrix A = {a ij } representing a self-adjoint operator  in an orthonormal basis {| 1 , | 2 , . . . , |i >, . . . , | N }. The algorithm consists of several iteration loops. The first loop goes through the following steps : a) Pick up the lowest n 0 basis states |i > (n 0 << N), construct and diagonalize the

n 0 × n 0 submatrix A 0 = (a ij ), b)select the lowest v eigenvalues Λ 0 ≡ {λ (0) 1 , . . . , λ (0) k , . . . , λ (0)
v } and the corresponding eigenvectors

| ϕ (0) k = n 0 i=1 c (0) i | i , (1)
c) consider the subspace spanned by these v eigenvectors | ϕ (0) k plus the basis states |j > = |n 0 + 1 >, . . . |n 1 > and construct the new submatrix

A 1 = Λ 0 B 01 B 10 C 1 , (2)
composed of the diagonal block Λ 0 ≡ {λ

(0) k } plus C 1 = {a ij } (i, j = n 0 + 1, n 1 )
. The two blocks are coupled by B 01 and its transpose, composed of the updated matrix elements b

(0) kj = ϕ (0) k | Â | j (k = 1, v), (j = n 0 + 1, n 1 ), d) diagonalize A 1 , extract the new lowest v eigenvalues Λ 1 ≡ {λ (1) k } (k = 1, v) and eigenvectors | ϕ (1)
k , and use them plus the basis states |j > (j = n 1 +1, n 2 ) to construct a new submatrix A 2 just as done in c).

The iterative procedure is now clearly outlined and can be illustrated by the following sequential operations

A 0 =⇒ Λ 0 -→ A 1 = Λ 0 B 01 B 10 C 1 =⇒ Λ 1 -→ . . . -→ A n = Λ (n-1) B (n-1)n B n(n-1) C n =⇒ Λ n = Λ. (3)
The first iteration loop goes on until the basis is exhausted, yielding the approximate v eigenvalues and eigenvectors

E (1) k ≡ λ (n) k , | ψ (1) k ≡| ϕ (n) k = N i=1 c (n) k (i) | i . (4)
These v eigenvectors |ψ > added to the original basis states {| i } are the entry for a new iteration loop. Since these states form a non orthogonal redundant basis, we solve a generalized eigenvalue problem and adopt the Choleski decomposition method to eliminate any redundancy. Apart from this modification, the subsequent iteration loops proceed as the first one and generate a sequence converging to the exact eigensolutions [13]

{E (r) k , ψ (r) k } r→∞ -→ {E k , ψ k }. (5)
As proved and illustrated through typical numerical tests [13], the algorithm is robust, being numerically stable and converging always to the extremal eigenvalues. It yields ghost-free solutions and is also of easy implementation. It requires O( N 2 ) operations. In fact we need N operations to compute the updated matrix element b

(n-1) kj = λ (n-1) k | Â | j (k = 1, v).
For sparse matrices with an average number L of non zero matrix elements, the required number of operations grows linearly with N.

Implementation of the algorithm in the m-scheme

Let  be a SM Hamiltonian H of a nucleus with p valence nucleons. These are supposed to be described by the basis states |i >= |α 1 , . . . α i , . . . , α p >, where α i = a i m i are the quantum numbers of a particle in the shell a i = {n i l i j i }. The states |i > have the total magnetic quantum number

M = m 1 + . . . m i + . . . m p .
It is advantageous to adopt the modified Hamiltonian

H J = H + c[J 2 -J(J + 1)] 2 , ( 6
)
where J is the total spin operator, J a given integer (or half-integer), and c a positive constant.

The full space can be thought to be decomposed into several subspaces

H = H 0 ⊕ H 1 . . . ⊕ H k . . . ⊕ H F . (7)
Each subspace H k is composed of a set of partitions {n i } k = {a n 1 1 , . . . a n i i , . . .} k , where i n i = p. It is, therefore, invariant with respect to J. The partitions in H k differ from those in H k-1 by at most two single particle shells a i .

We proceed according to the following iterative scheme

H 0 ⇒ Λ 0 -→ Λ 0 ⊕ H 1 ⇒ Λ 1 -→ . . . -→ Λ F -1 ⊕ H F ⇒ Λ F . (8)
We, namely, diagonalize the Hamiltonian H J in H 0 obtaining v lowest eigenvalues

E (0) 1 , . . . , E (0) v and eigenvectors ψ (0) 1 , . . . , ψ (0)
v spanning a subspace Λ 0 . These eigensolutions are exact in this subspace and have all the same spin J if the constant c is chosen so as to push the states with J ′ = J up in energy.

Because of its two-body nature, the Hamiltonian couples the subspace Λ 0 to H 1 only. We therefore diagonalize H J in the subspace Λ 0 ⊕ H 1 to generate new updated eigenvalues E

(1) 1 , . . . , E (1) v and eigenvectors ψ (1) 1 , . . . , ψ (1) v , defining the subspace Λ 1 . We proceed iteratively. Once the updated eigensolutions defining the subspace Λ k are obtained, we diagonalize the Hamiltonian in the upgraded subspace Λ k ⊕ H k+1 . We cover eventually the full space obtaining the exact v eigensolutions {E i , ψ i }.

By exploiting the sparsity of the Hamiltonian matrix, the above procedure allows to reduce the number of operations to O(N), as in Lanczos implementations. This, however, may not be sufficient for spaces of very large dimensions. It is therefore useful to search for a reliable way of cutting the basis. To this purpose, the original sampling procedure [14] has been modified and adapted to the m-scheme.

| < j|H J |ϕ (k-1) i > | 2 a jj -E (k-1) i > ǫ k . (9)
More specifically, in the first step (k = 1), the above condition selects a set of states {|j >} forming a subspace H

(ǫ 1 ) 1 ∈ H 1 . The eigenvalue problem is thus solved in Λ 0 ⊕ H 1 (ǫ 1 ) yielding new v eigensolutions E (1) k (ǫ 1 ), ψ (1)
k (ǫ 1 ) defining the subspace Λ (ǫ 1 ) 1 . We now explore the full subspace complementary to H 0 ⊕ H 1 (ǫ 1 ) and select all the states |j > that fulfill the condition (9) with ǫ 2 replacing ǫ 1 . The states so selected span a subspace H

(ǫ 2 ) 2 ∈ H 0 ⊕ H 1 ⊕ H 2 .
The above procedure is iterated with updated eigensolutions and decreasing sampling values ǫ k until the full space is covered.

A sketch of the sampling procedure is provided by the following sequence

H 0 ⇒ Λ 0 ǫ 1 -→ Λ 0 ⊕ H (ǫ 1 ) 1 ⇒ Λ (ǫ 1 ) 1 ǫ 2 -→ . . . ǫ F -→ Λ (ǫ F -1 ) F -1 ⊕ H (ǫ F ) F ⇒ Λ (ǫ F )
F .(10) Clearly, we get the exact solutions for ǫ F → 0, namely lim ǫ F →0 Λ

(ǫ F ) F = Λ F , or, more explicitly, lim ǫ F →0 {E i (ǫ F ), ψ i (ǫ F )} = {E i , ψ i }. (11)
Few remarks are mandatory. Like the algorithm sketched in Section 2, the sampling procedure is extremely stable. It yields, indeed, orthonormal eigenfunctions. These, at each step, result from the diagonalization of a symmetric submatrix. The subspaces selected by the sampling are not strictly invariant with respect to J. The invariance, however, is restored as the sampling value ǫ becomes sufficiently small. We will show, in fact, that all sampled eigenstates reach soon a good J.

The just outlined sampling process reminds of the DMRG [10]. In both approaches, in fact, one starts with a subspace of small dimensions n 0 and then goes through several steps consisting of enlarging the space and, then, going back to a restricted space of the original dimensions n 0 spanned by updated eigensolutions. There are crucial differences, however. Ours is based on a new diagonalization algorithm. Moreover, in our case the space truncation is achieved by updating both energies and wavefunctions, while in the DMRG it is based on prescriptions imposed on the wave functions (density matrix).

Numerical implementation and results

We applied the sampling algorithm to 132 Xe and 130 Xe. These nuclei are supposed to be near the transition point from a vibrational to a γ-soft structure and are being object of experimental investigations [15]. The model space includes the shells {2d5/2, 1g7/2, 2d3/2, 3s1/2, 1h11/2} for both valence neutrons and protons. The dimensions of the resulting Hamiltonian matrix are N ≃ 3.7 × 10 7 for 132 Xe and N ≃ 0.8 × 10 9 for 130 Xe.

It is advantageous to describe the above isotopes in terms of valence proton particles external to a Z = 50 core and neutron holes referredl to a N = 82 core. The single particle energies are shown in Table 1. The neutron hole energies were deduced from the spectrum of 131 Sn, while the proton particle energies were taken from Ref. [14].

A renormalized G matrix [16] derived from the CD-Bonn potential [17] was taken as a two-body interaction. The E2 transitions were computed using the effective charges e p = 1.5 and e n = 0.7 for protons and neutrons, respectively.

The first task consisted in choosing the initial subspace H 0 . Its dimensions n 0 increase with the number of eigenstates of good J we intend to generate. To yield up to ten eigenstates, the space dimensions came out to be of the order n 0 ∼ 100. Having chosen H 0 , we applied the iterative sampling procedure (10) with decreasing values of ǫ k . Each ǫ k determines uniquely the dimension n k of the Hamiltonian matrix to be In the case of 132 Xe, the iterative procedure was carried out so as to include up to 40% ÷ 50% of the basis states. As shown in Fig. 1, about 10% of the basis states is enough to lead the J π = 0 + and J π = 2 + eigenvalues to convergence. Indeed, the J π = 0 + absolute energies E k decrease by less than 80 keV in going from n/N ∼ 0.08 to n/N ∼ 0.46 ( Table 2). The spectrum converges even faster. Indeed, the energy difference E(0 + 2 ) -E(0 + 1 ) changes by ∼ 10 keV when we move from n/N ∼ 0.08 to n/N ∼ 0.46.

The same convergence rate is found for the levels of 130 Xe (Figs. 2 and Table 2). Indeed, the fraction of states considered (7% ÷ 8%) is enough to bring the energy eigenvalues to the plateau. On the ground of the analysis made for 132 Xe, the energies obtained here may differ from the asymptotic values (n/N = 1) by at most 100 keV or by 10 ÷ 20 keV, when referred to the ground state.

All the corresponding eigenstates have good J. As shown in Table 2, even with 7%÷ 8% of the basis states, the J values coincide with the exact ones up to the second or third decimal digit. Thus, although the subspaces selected according to the sampling criterion (9) are not strictly J-invariant, the invariance is eventually restored by a relatively small fraction of basis states. Apparently, the components left out by the sampling are small and of little relevance to all observables. This statement is supported by the good convergence of the E2 transition strengths. In 132 Xe, the E2 strengths reach their saturated values for n/N ∼ 0.1. In going from n/N ∼ 0.08 to n/N ∼ 0.46, the strengths increases from B(E2) ∼ 694 e 2 fm 4 to B(E2) ∼ 708 e 2 fm 4 , a ∼2% increment. The same convergence rate is obtained for the B(E2) values in 130 Xe. The values obtained for n/N ∼ 0.08 may be smaller than the asymptotic one by ∼ 2% ÷ 3%.

The good convergence of the angular momentum J and the B(E2) probes the accuracy of the sampled wavefunctions.

Being interested in illustrating the method, we have refrained ourselves from making any comparison with experiments. We may only say that the computed levels and transition strengths are not far from the corresponding measured quantities [15]. By few small changes, like a better tuning of the single particle energies, the calculation should be able to provide a satisfactory description of the low-energy spectroscopic properties of these two important Xe isotopes.

Conclusions

From the present analysis we may conclude that the new m-scheme implementation of the algorithm allows to take advantage of the sparsity of the Hamiltonian. The importance sampling, combined with sparsity, enhances the convergence of the iterative process leading to an effective truncation of the basis.

The algorithm is able to generate an almost arbitrary number of extremal eigenstates for a given J π and, therefore, can describe spectra of great complexity. All these eigenstates are automatically orthonormal and have good J values.

Its stability is to be pointed out once again. Indeed, the iterative procedure reaches convergence, independently of the subspace we start with.

As the sizes of the Hamiltonian matrix become too large, the time required to perform the iterations becomes too long, at least on a desktop. On the other hand, the smooth monotonic behavior of energies and transition strengths versus n/N allows to make reliable and straightforward extrapolations to their asymptotic values (n/N = 1).

If we stick on a desktop, the extrapolation may allow to treat nuclei, like 128 Xe, at most. In order to go further, we need to develop a parallel version of the code and explore if and to what extent its performance is enhanced.

Figure 1 .

 1 Figure 1. Convergence rate of the lowest J π = 0 + and J π = 2 + eigenvalues in 132 Xe.

Figure 2 .

 2 Figure 2. Convergence rate of the lowest J π = 0 + and J π = 2 + eigenvalues in 130 Xe.

Figure 3 .

 3 Figure 3. Convergence properties of the E2 reduced strengths.

Table 1 .

 1 Single proton particle and single neutron hole energies (in MeV).Let us fix a sequence of positive small numbers of decreasing values ǫ 1 > . . . ǫ k > . . . > ǫ F . Having generated the lowest v eigenvectors ψ , we proceed following the same scheme (8) of the exact case with one constraint: In going from H k-1 to H k , we pick up only the basis states |j > that fulfill the condition

	(nlj) p	2d 5/2	1g 7/2	2d 3/2	3s1/2	1h 11/2
	ǫ p	0.00	0.2	2.2	2.3	2.9
	(nlj) -1 n	(2d 3/2 ) -1	(1h 11/2 ) -1	(3s 1/2 ) -1	(2d 5/2 ) -1	(1g 7/2 ) -1
	ǫ n	0.00	0.14	0.33	1.65	2.43
					(0) 1 , . . . , ψ (0) v	in H 0

Table 2 .

 2 Convergence rate of the absolute energies (in MeV) in 130 Xe and 132 Xe.

	132 Xe J π = 0 + J samp			(n/N) k
			0.0 81	0.219	0.290	0.372	0.463
	0 + 1	0.0005 E k = -41.138 -41.201 -41.209 -41.213 -41.215
	0 + 2	0.0004 E k = -39.857 -39.912 -39.918 -39.922 -39.924
	J π = 2 + J samp			(n/N) k
			0.086	0.166	0.223	0.293	0.374
	2 + 1	2.0002 E k = -40.273 -40.329 -40.343 -40.351 -40.356
	2 + 2	2.0001 E k = -39.720 -39.774 -39.786 -39.794 -39.799
	130 Xe J π = 0 + J samp			(n/N) k
			0.011	0.018	0.029	0.045	0.067
	0 + 1	0.012	E k = -59.141 -59.273 -59.369 -59.435 -59.481
	0 + 2	0.012	E k = -57.942 -58.074 -58.167 -58.229 -58.272
	J π = 2 + J samp			(n/N) k
			0.013	0.021	0.034	0.053	0.079
	2 + 1	2.002	E k = -58.365 -58.517 -58.621 -58.692 -58.740
	2 + 2	2.002	E k = -57.699 -57.866 -57.994 -58.082 -58.144
	diagonalized.

Acknowledgments

We thank M. Hjorth-Jensen and N. Pietralla for useful information and stimulating discussions.