text:
# 1 Introduction

## 1 Introduction

Properties and evolution of interfaces is a very interesting topic in soft condensed matter physics, see, e. g.,  for a recent review. There are two main lines of theoretical research in that direction. The first one, computer simulations, provides a relation to microscopic, molecular dynamics level . In the second approach one is satisfied with a coarse-grained description by effective Ginzburg-Landau models.

In the framework of Ginzburg-Landau models the interfaces are described by particular solutions of pertinent nonlinear evolution equations. The nonlinearity is crucial for existence of these solutions, but it is also a serious obstacle in obtaining them in an explicit form. In most cases only numerical solutions are available. A well-known exception is a model which appears in the context of reaction-diffusion chemical processes , and also is related to Landau-deGennes effective model for nematic liquid crystals . It involves a single, non-conserved, real-valued order parameter $`\phi (x^\alpha ,t)`$, and the free energy functional has the form

$$F_0=d^3x\left(\frac{1}{2}\phi \phi +V_0(\phi )\right),$$
(1)
$$V_0(\phi )=A\phi ^2+B\phi ^3+\phi ^4,$$
(2)
where $`=\left(\frac{}{x^\alpha }\right)`$, $`\alpha =1,2,3`$, $`x^\alpha `$ are Cartesian coordinates in the space. For simplicity we have rescaled $`x^\alpha ,\phi `$ and $`F_0`$ so that they are dimensionless and the coefficient in front of the quartic term in $`V_0(\phi )`$ is equal to $`1`$. The constants $`A,B`$ are positive, and $`A`$ is restricted to the interval $`0<A<\frac{9B^2}{32}`$. Then $`V(\phi )`$ has two minima and a maximum between them. In this model exact formulas for profile and velocity of a planar interface are known, see, e.g., .

Exactly soluble models provide natural starting points for perturbative expansions. If the expression (1) for the free energy is changed a little bit, one may hope that the corresponding change of the planar interface can be calculated with the help of a perturbative expansion around the known solution. In this way we would obtain analytical results for velocity and profile of planar interfaces in a class of Ginzburg-Landau models. Moreover, this would open the way for investigations of evolution of curved interfaces in that class of models, because one could use the expansion in curvatures  for which the knowledge of planar interfaces is the only prerequisite.

In the present paper we explore the idea of the perturbative expansion for the planar interfaces. We show that indeed such perturbative scheme can be constructed. Corrections to the profile are obtained in it as solutions of ordinary, linear differential equations, while corrections to the velocity follow from integrability conditions which are due to existence of zero-modes. This type of integrability conditions is well-known in statistical physics . As an example we consider perturbations of the form

$$\delta F=d^3xV_1,$$

where

$$V_1=\frac{1}{2}\epsilon _1\phi ^2\phi \phi +\frac{1}{5}\epsilon _2\phi ^5+\frac{1}{6}\epsilon _3\phi ^6.$$
(3)
Values of the perturbation parameters $`\epsilon _1,\epsilon _2,\epsilon _3`$ are such that $`F_0+\delta F`$ is bounded from below. The $`\epsilon _1`$ correction to the gradient term introduces dependence on $`\phi `$ of the elastic constant $`K=1+\epsilon _1\phi ^2`$. The $`\epsilon _2`$ and $`\epsilon _3`$ terms are corrections to the potential $`V_0(\phi )`$. The perturbations (3) are vanishing in the disordered phase ($`\phi =0`$), while in the ordered phase ($`\phi 0`$) they may give finite contributions.

Our paper is organised as follows. In Section 2 we construct the perturbative expansion. In Section 3 we apply it to the model (1), (2) with the particular perturbation $`\delta F`$ given by formula (3). Summary and remarks are presented in Section 4. In the Appendix we show that the perturbatively calculated velocity of the planar interface coincides with the one calculated from an exact formula, provided that the perturbation is not too strong and that the perturbative series is convergent.

## 2 General structure of the perturbative expansion

We will consider Ginzburg-Landau models in which time evolution of the order parameter $`\phi `$ is governed by equation of dissipative type, namely

$$\frac{\phi }{t}=\mathrm{\Delta }\phi V^{}(\phi ,\phi ),$$
(4)
where $`V=V_0+\epsilon V_1`$, $`t`$ denotes a rescaled dimensionless time, $`\mathrm{\Delta }=_{x^\alpha }_{x^\alpha }`$ is the three-dimensional Laplacian, and

$$V^{}=\frac{V}{\phi }\left(\frac{V}{(\phi )}\right).$$

We have assumed that the perturbation $`V_1`$ does not depend on second and higher derivatives of $`\phi `$. The l.h.s. of Eq. (4) is equal to the variational derivative $`\delta F/\delta \phi `$ multiplied by $`1`$.

The interface is identified with a solution of the evolution equation which interpolates between minima $`\phi _+,\phi _{}`$ of the potential $`V`$. Potential $`V_0`$ given by formula (2) has two minima, and the perturbation $`V_1`$ is by assumption weak enough, so that also $`V`$ has just two minima which smoothly merge with the minima of $`V_0`$ when $`V_1`$ is switched off. In the present paper we concentrate on planar interfaces. Due to translational invariance of Eq.(4), it is sufficient to consider the planar interfaces parallel to the plane $`(x^1,x^2)`$ \- then $`\phi `$ is a function of $`x^3`$ and $`t`$. We require that in the limits $`x^3\pm \mathrm{}`$,

$$\phi (x^3,t)\phi _\pm ,_{x^3}\phi 0.$$
(5)
Let us multiply Eq.(4), in which now $`\phi _{x^3}\phi (x^3,t)`$, by $`_{x^3}\phi `$ and integrate it over $`x^3`$ in the interval $`(\mathrm{},+\mathrm{})`$. With the help of integration by parts, and using a formula for the full derivative of $`V`$ with respect to $`x^3`$, namely

$$\frac{dV}{dx^3}=\frac{V}{\phi }_{x^3}\phi +\frac{V}{(_{x^3}\phi )}_{x^3}^2\phi ,$$

we obtain the following relation

$$_{\mathrm{}}^+\mathrm{}𝑑x^3_t\phi _{x^3}\phi =V(\phi _{})V(\phi _+).$$
(6)
It implies that $`\phi `$ has to be time-dependent if the minima $`\phi _\pm `$ are nondegenerate. The simplest way to ensure that the l.h.s. of formula (6) is constant in time and nonvanishing consists in assuming that $`\phi `$ depends only on the combination $`x^3\upsilon t`$,

$$\phi =\phi (x^3\upsilon t),$$
(7)
where $`\upsilon `$ is a constant equal to the velocity of the interface. This form of dependence on time is also consistent with boundary conditions (5). For such $`\phi `$, relation (6) can be written in the form

$$\upsilon =\frac{V(\phi _+)V(\phi _{})}{_{\mathrm{}}^+\mathrm{}𝑑z(_z\phi )^2},$$
(8)
where $`z=x^3\upsilon t`$. Equation (4) is reduced to

$$_z^2\phi +\upsilon _z\phi V^{}(\phi ,_z\phi )=0,$$
(9)
where $`\phi `$ is a function of $`z`$ only, and

$$\phi (z)\phi _\pm $$
(10)
when $`z\pm \mathrm{}`$, respectively. The asymptotic values $`\phi _\pm `$ are determined from the equations

$$V^{}(\phi _\pm )=0.$$
(11)
From a mathematical viewpoint (8) and (9) should be regarded as a set of equations for $`\upsilon `$ and $`\phi (z)`$, with (10) as the boundary conditions for $`\phi (z)`$.

The perturbative Ansatz for $`\phi (z)`$ has the form

$$\phi (z)=\phi _0(z)+\epsilon \phi _1(z)+\epsilon ^2\phi _2(z)+\mathrm{}.$$
(12)
Here $`\phi _0(z)`$ is the initial, unperturbed interface. It obeys the following equation

$$_z^2\phi _0+\upsilon _0_z\phi _0V_0^{}(\phi _0)=0,$$
(13)
with the boundary conditions at $`z\pm \mathrm{}`$

$$\phi _0(z)a_\pm ,$$

where $`a_\pm `$ are the minima of $`V_0`$. The velocity $`\upsilon _0`$ of the unperturbed interface is given by the formula analogous to (8), namely

$$\upsilon _0=\frac{V_0(a_+)V_0(a_{})}{_{\mathrm{}}^+\mathrm{}𝑑z(_z\phi _0)^2}.$$
(14)
Formula (8) suggests that the velocity of the interface has the following perturbative expansion

$$\upsilon =\upsilon _0+\epsilon \upsilon _1+\epsilon ^2\upsilon _2+\mathrm{}.$$
(15)
Let us stress that in equations (8), (9) and in expansion (12) the independent variable is just $`z`$, and not $`x^3`$ and $`t`$. Therefore, when solving Eqs. (8), (9), the variable $`z`$ is not treated as $`x^3\upsilon t`$, and therefore we do not expand $`z`$ in $`\epsilon `$, contrary to what formula (15) might suggest. Only after $`\phi (z)`$ and $`\upsilon `$ are determined we may substitute $`z=x^3\upsilon t`$ in order to relate the solution of the set of equations (8), (9) with the interface solution of the original evolution equation (4).

Inserting the expansions (12), (15) in Eq.(9) and equating to zero the l.h.s. of (9) order by order in $`\epsilon `$, we obtain the infinite chain of equations for the corrections $`\phi _k,k1`$. These equations have the form

$$\widehat{L}\phi _k=f_k,$$
(16)
where

$$\widehat{L}=_z^2+\upsilon _0_zV_0^{\prime \prime }(\phi _0),$$
(17)
and $`f_k`$ in general depend on $`\upsilon _1,\mathrm{},\upsilon _k`$ and $`\phi _1,\mathrm{},\phi _{k1}`$. Simple calculations give

$$f_1=\upsilon _1_z\phi _0+V_1^{}(\phi _0,_z\phi _0),$$
(18)
and for $`k2`$

$$f_k=\upsilon _k_z\phi _0\underset{j=1}{\overset{k1}{}}\upsilon _j_z\phi _{kj}+h_k(z),$$
(19)
where $`h_k(z)`$ denotes the sum of all terms of the $`k`$-th order in $`\epsilon `$ obtained by inserting formula (12) in $`V_0^{}(\phi )+\epsilon V_1^{}(\phi ,_z\phi )`$ with exception of the term $`V_0^{\prime \prime }(\phi _0)\phi _k`$ which has been included into the l.h.s. of Eq.(16).

It is a well-known fact that inhomogeneous linear equations not always have solutions - certain integrability conditions have to be satisfied. In the case of Eq.(16) such conditions appear because for the operator $`\widehat{L}`$ there exists so called left zero-mode, that is a normalizable function $`\psi _l(z)`$ such that the following identity holds

$$_{\mathrm{}}^+\mathrm{}𝑑z\psi _l(z)\widehat{L}\phi _k=0.$$
(20)
If $`\phi _k`$ obeys Eq.(16) then

$$_{\mathrm{}}^+\mathrm{}𝑑z\psi _l(z)f_k=0.$$
(21)
In our perturbative scheme the integrability conditions (21) serve as equations which determine the perturbative contributions $`\upsilon _k`$ to the velocity of the interface. Using formulas (18), (19) and conditions (21) we find that

$$\upsilon _1=N^1_{\mathrm{}}^+\mathrm{}𝑑zV_1^{}(\phi _0,_z\phi _0)\psi _l(z),$$
(22)
and for $`k2`$

$$\upsilon _k=N^1_{\mathrm{}}^+\mathrm{}𝑑z\left(h_k(z)\underset{j=1}{\overset{k1}{}}\upsilon _j_z\phi _{kj}\right)\psi _l(z),$$
(23)
where

$$N=_{\mathrm{}}^+\mathrm{}𝑑z\psi _r(z)\psi _l(z).$$

The functions $`\psi _l,\psi _r`$ will be given shortly. In the Appendix we check that these recursive formulas and the expansion (15) give the velocity $`\upsilon `$ which coincides with the one obtained by expanding in $`\epsilon `$ the r.h.s. of formula (8), provided that $`\epsilon `$ is not too large.

The left zero-mode $`\psi _l`$ can be found in the following way. Differentiation of the both sides of Eq.(13) with respect to z gives the identity

$$\widehat{L}\psi _r=0,$$
(24)
where

$$\psi _r=_z\phi _0(z).$$
(25)
The function $`\psi _r`$ is called the right zero-mode. $`\psi _l`$ has the form

$$\psi _l(z)=\mathrm{exp}(\upsilon _0z)\psi _r(z).$$
(26)
Because derivation of identity (20) involves integration by parts, we have to discuss behaviour of $`\psi _l`$ and $`\psi _r`$ for $`z\pm \mathrm{}`$. Notice that Eq.(24) implies that for $`z\pm \mathrm{}`$

$$\psi _r(z)\mathrm{exp}\left(\frac{\upsilon _0}{2}z\right)\mathrm{exp}\left(\sqrt{\frac{\upsilon _{0}^{}{}_{}{}^{2}}{4}+V_0^{\prime \prime }(a_\pm )}z\right),$$

respectively. Therefore $`\psi _l(z)`$ exponentially vanishes for $`z\pm \mathrm{}`$, provided that

$$V_0^{\prime \prime }(a_\pm )>0.$$

For $`V_0`$ given by formula (2) this condition is satisfied if $`0<A<\frac{9B^2}{32}`$, as it has been assumed. When checking identity (20), boundary terms like, e.g., $`(_z\psi _l)\phi _k|_{\mathrm{}}^+\mathrm{}`$, vanish because

$$\phi _k\text{const}\text{for}z\pm \mathrm{},$$
(27)
in accordance with the boundary condition (10).

Equations (16) can be solved in a standard way . Adopting formulas given in  to the case at hand we find that

$$\phi _k=\psi _2(z)_{\mathrm{}}^z𝑑\zeta e^{\upsilon _0\zeta }\psi _r(\zeta )f_k(\zeta )\psi _r(z)_0^z𝑑\zeta e^{\upsilon _0\zeta }\psi _2(\zeta )f_k(\zeta ),$$
(28)
where

$$\psi _2(z)=\psi _r(z)_0^z𝑑\zeta e^{\upsilon _0\zeta }\psi _r^2(\zeta ).$$
(29)
The functions $`\psi _r,\psi _2`$ form a pair of linearly independent solutions of the homogeneous equation $`\widehat{L}\phi =0`$. Notice that the integral $`_{\mathrm{}}^z`$ present in the first term on the r.h.s. of formula (28) vanishes for $`z+\mathrm{}`$ due to the integrability conditions (21). It turns out that in the limits $`z\pm \mathrm{}`$ both terms on the r.h.s. of formula (28) approach constants, and $`\phi _k`$ obey the conditions (27).

General solution of Eqs.(16) is obtained by adding to the r.h.s. of formula (28) the general solution of the corresponding homogeneous equation, that is the function $`c_k\psi _r(z)+d_k\psi _2(z)`$. However, because $`\psi _2(z)`$ exponentially grows for large $`|z|`$ we have to put $`d_k=0`$. In order to determine the constants $`c_k`$ we have to impose a condition on the interface solution $`\phi (z)`$ in addition to the boundary conditions (10). Actually, Eq.(9) and the boundary conditions (10) do not determine $`\phi (z)`$ uniquely - due to invariance with respect to translations in $`z`$ we can take $`\phi (zz_0)`$ with arbitrary $`z_0`$. We will require that

$$\phi (0)=\phi _0(0),$$
(30)
where $`\phi _0(z)`$ is a concrete, explicitly given function describing the unperturbed interface with fixed location on the $`z`$ axis. The condition (30) breaks the translational invariance. It implies that

$$\phi _k(0)=0.$$
(31)
The solutions $`\phi _k(z)`$ given by formula (28) obey this condition, therefore also the constants $`c_k`$ are vanishing. Thus, the perturbative scheme supplemented with the condition (30) yields unique interface solution $`\phi (z)`$. In this solution we can of course shift the variable $`z`$, that is to substitute $`zzz_0`$ simultaneously in all contributions $`\phi _0,\phi _k`$, in order to obtain the other solutions implied by the translational invariance.

## 3 Effects of the perturbations of the form <br>$`V_1=\frac{1}{2}\epsilon _1\phi ^2\phi \phi +\frac{1}{5}\epsilon _2\phi ^5+\frac{1}{6}\epsilon _3\phi ^6`$.

Let us start from a description of the unperturbed interface $`\phi _0(z)`$. The evolution equation has the form (4) with $`V`$ replaced by $`V_0`$ given by formula (2). The interface solution is well-known, see, e.g., . Let us quote the relevant formulas from  and . We use the following abbreviations

$$s=\frac{z}{4l_0},a_+=\frac{1}{2\sqrt{2}l_0},$$

where

$$l_{0}^{}{}_{}{}^{1}=\frac{1}{2\sqrt{2}}\left(3B+\sqrt{9B^232A}\right).$$
(32)
Then

$$\phi _0(z)=\frac{a_+}{2}\left(1+\mathrm{tanh}s\right),$$
(33)
and

$`\psi _r(z)`$ $`=`$ $`{\displaystyle \frac{a_+}{8l_0}}{\displaystyle \frac{1}{\mathrm{cosh}^2s}},`$
$`\psi _2(z)`$ $`=`$ $`{\displaystyle \frac{2l_{0}^{}{}_{}{}^{2}}{a_+}}{\displaystyle \frac{1}{\mathrm{cosh}^2s}}{\displaystyle _0^s}𝑑\sigma \mathrm{exp}[2(2\alpha +1)\sigma ]\left[1+\mathrm{exp}(2\sigma )\right]^4.`$
The integral over $`\sigma `$ can be easily calculated, but the resulting formula for $`\psi _2`$ is quite long. In the absence of perturbations, $`z=x^3\upsilon _0t`$, where

$$\upsilon _0=\frac{3}{4\sqrt{2}}\left(B\sqrt{9B^232A}\right).$$
(34)
The interface $`\phi _0`$ exists and it is stable when $`0<A<\frac{9B^2}{32}`$. Usually, the coefficient $`B`$ is regarded as independent of temperature, while $`A=a(TT_{})`$, where $`a>0`$. The inequality given above fixes the temperature range $`(T_{},T_c)`$ in which the interface exists. Notice that at the temperature $`T_0`$ such that $`A=\frac{B^2}{4}`$ the velocity $`\upsilon _0`$ vanishes. It is clear that $`T_{}<T_0<T_c`$. The left zero-mode has the form

$$\psi _l(z)=\frac{a_+}{8l_0}\mathrm{exp}[2(2\alpha 1)s]\frac{1}{\mathrm{cosh}^2s},$$
(35)
where

$$\alpha =\upsilon _0l_0+\frac{1}{2}.$$

The value of the parameter $`\alpha `$ monotonically increases from 0 to 1 when $`A`$ varies from 0 to $`\frac{9B^2}{32}`$.

The first order correction $`\upsilon _1`$ to the velocity is calculated from formula (22). In the present case

$$V_1^{}(\phi _0)=\epsilon _1\phi _0\psi _{r}^{}{}_{}{}^{2}\epsilon _1\phi _{0}^{}{}_{}{}^{2}_z\psi _r+\epsilon _2\phi _{0}^{}{}_{}{}^{4}+\epsilon _3\phi _{0}^{}{}_{}{}^{5}.$$
(36)
The integrations over $`z`$ ( or equivalently over $`s`$ ) can be related to the Euler $`B`$ function ,

$$B(x,y)=_0^{\mathrm{}}𝑑t\frac{t^{x1}}{(1+t)^{x+y}},$$

by the substitution $`y=\mathrm{exp}s`$. The resulting expressions can be simplified with the help of the well-known formulas

$$B(x,y)=\frac{\mathrm{\Gamma }(x)\mathrm{\Gamma }(y)}{\mathrm{\Gamma }(x+y)},\mathrm{\Gamma }(x+1)=x\mathrm{\Gamma }(x).$$

After straightforward calculations we obtain

$$\upsilon _1=\frac{(\alpha +1)(2\alpha +1)}{160l_{0}^{}{}_{}{}^{3}}\left[(\alpha \frac{1}{2})\epsilon _1\frac{l_0}{\sqrt{2}}\frac{2\alpha +3}{1\alpha }\epsilon _2+\frac{(2\alpha +3)(\alpha +2)}{12(1\alpha )}\epsilon _3\right].$$
(37)
Formula (37) shows that the $`\epsilon _2,\epsilon _3`$ corrections become more and more pronounced as $`\alpha `$ increases towards $`1`$, that is as the temperature increases towards $`T_c`$. The singularity at $`\alpha =1`$ is due to the fact that for $`\alpha =1`$ the left zero-mode $`\psi _l`$ and $`V_1^{}(\phi _0)`$ do not vanish when $`s+\mathrm{}`$, and therefore the integral in formula (22) is divergent.

The $`\epsilon _1`$ term in formula (37) vanishes for $`\alpha =\frac{1}{2}`$. The reason is that the corresponding term in the perturbation potential $`V_1`$ does not contribute to the values of the full potential $`V`$ at the two minima $`\phi _\pm `$. Therefore, it can influence the velocity only by changing the denominator in formula (8). For $`\alpha =\frac{1}{2}`$ the numerator in (8) vanishes if $`\epsilon _2=\epsilon _3=0`$ and the $`\epsilon _1`$ term can not contribute to the velocity.

The facts that $`\epsilon _2`$ term in (37) is negative, and that the $`\epsilon _3`$ term is positive, can be explained by the observation that the $`\epsilon _2`$ term in the potential diminishes the potential energy difference across the interface, while the $`\epsilon _3`$ term in the potential works in the opposite direction.

## 4 Summary and discussion

1. The purpose of this work is to provide a practical tool for analytic, approximate computations of characteristics of the perturbed planar interfaces within the framework of Ginzburg-Landau models. We have shown how one can systematically calculate corrections to the profile and the velocity of the interface. In particular, the scheme can be applied to perturbations of domain walls, which can be regarded as special, static interfaces appearing when the two minima of $`V_0`$ are degenerate. In the example discussed in Section 3 this is the case when $`\alpha =1/2`$. The scheme is relatively simple. The final formulas (22), (23), (28), (29) contain one dimensional integrals which always can be tackled numerically, and in many cases calculated analytically.

2. We have not delved into the problem of convergence of the perturbative series. Such a mathematical investigation does not belong to the scope of our work. It is clear that if the perturbations shift the two minima of $`V_0`$ only slightly, the corresponding changes of the profile and of the velocity of the interface are also small, hence they can be calculated as small corrections. This does not necessarily mean that the series is convergent – it can belong to a wider class of asymptotic series, but it is acceptable from physical viewpoint.

3. Using the results of the present paper in combination with the curvature expansion mentioned in the Introduction, one can calculate evolution of curved interfaces in the perturbed Ginzburg-Landau models. Comparisons of such theoretical predictions with experimental observations of evolution of the interfaces could help to determine the best formula for the free energy of the system. Let us remind that the formulas for the free energy in the Ginzburg-Landau models rarely can be derived from underlying microscopic theories. In most cases they are postulated on basis of qualitative phenomenological considerations. General discussion of various formulas for the free energy in the context of theory of liquid crystals can be found in .

## 5 Appendix. Resummation of the perturbative contributions to the velocity

The perturbative corrections $`\upsilon _k`$ to the velocity of the interface are given by formulas (22), (23). They have been obtained in somewhat indirect way, namely from the integrability conditions. On the other hand, the velocity is given by the exact formula (8). We would like to check that the sum (15) of the perturbative contributions $`\upsilon _k`$, below denoted by $`\stackrel{~}{\upsilon }`$,

$$\stackrel{~}{\upsilon }=^{df}\upsilon _0+\epsilon \upsilon _1+\underset{k=2}{\overset{\mathrm{}}{}}\epsilon ^k\upsilon _k,$$
(38)
does coincide with the velocity $`\upsilon `$ given by formula (8). Actually, because formula (8) follows directly from Eq.(9), it is sufficient to check that $`\stackrel{~}{\upsilon }`$ coincides with the $`\upsilon `$ present in that equation.

In the first step, we substitute for $`\upsilon _0,\upsilon _1`$ and $`\upsilon _k`$ ( $`k2`$ ) in (38) formulas (14), (22) and (23), respectively. Next, we notice that the definition of $`h_k`$ given below formula (19) is equivalent to the following formula

$$\underset{k=2}{\overset{\mathrm{}}{}}\epsilon ^kh_k=V_0^{}(\phi )+\epsilon V_1^{}(\phi ,_z\phi )V_0^{}(\phi _0)\epsilon V_1^{}(\phi _0,_z\phi _0)V_0^{\prime \prime }(\phi _0)(\phi \phi _0),$$

where

$$V_0^{}(\phi )+\epsilon V_1^{}(\phi ,_z\phi )=_{z}^{}{}_{}{}^{2}\phi +\upsilon _z\phi ,$$

$$V_0^{}(\phi _0)=_z\psi _r+\upsilon _0\psi _r,$$

and

$$V_0^{\prime \prime }(\phi _0)\phi _0=_z\psi _r+\upsilon _0\psi _r\widehat{L}\phi _0,$$

according to Eqs. (9), (13) and definition (17). We also use the formula

$$\underset{k=2}{\overset{\mathrm{}}{}}\epsilon ^k\underset{j=1}{\overset{k1}{}}\upsilon _j_z\phi _{kj}=\left(\underset{i=1}{\overset{\mathrm{}}{}}\epsilon ^i\upsilon _i\right)\left(\underset{j=1}{\overset{\mathrm{}}{}}\epsilon ^j_z\phi _j\right)=(\stackrel{~}{\upsilon }\upsilon _0)_z(\phi \phi _0).$$

After a simple calculation we obtain from formula (38) that

$$(\upsilon \stackrel{~}{\upsilon })_{\mathrm{}}^+\mathrm{}𝑑z\psi _l_z\phi =0.$$
(39)
Thus, indeed

$$\upsilon =\stackrel{~}{\upsilon },$$

provided that

$$_{\mathrm{}}^+\mathrm{}𝑑z\psi _l_z\phi 0.$$

Because

$$_z\phi =\psi _r+\epsilon _z\phi _1+\mathrm{},$$

that last condition is satisfied if $`\epsilon `$ is small enough, that is if the perturbation is not too strong.


format: text

text:
# A novel scenario for the production of antihyperons in relativistic heavy ion collisions

## Abstract

We elaborate on our recent suggestion on antihyperon production in relativistic heavy ion collisions by means of multi-mesonic (fusion-type) reactions. It will be shown that the (rare) antihyperons are driven towards chemical equilibrium with pions, nucleons and kaons on a timescale of 1–3 fm/c in a still moderately baryon-dense hadronic environment.

The original idea behind the collective enhancement of strangeness out of a potentially deconfined state of matter, the quark gluon plasma (QGP), is that the strange (and antistrange) quarks are thought to be produced more easily and hence also more abundantly as compared to the production via highly threshold suppressed inelastic hadronic collisions. In addition, the analysis of measured abundancies of hadronic particles within thermal models supports the idea of having established a (thermodynamically) equilibrated hadronic fireball in some late stage of the reaction (for analyses of Pb+Pb collisions at CERN-SPS see ). In this respect especially a nearly fully chemically equilibrated yield of strange antibaryons, the antihyperons, had originally been advocated as an appropriate QGP signature . Although intriguing, this may not be the correct (or only) interpretation of the observed antihyperon yields: In the following we will elaborate on our recent idea  of rapid antihyperon production by multi-mesonic reactions like $`n_1\pi +n_2K\overline{Y}+p`$ (see fig. 1 for a particular illustration). As we will demonstrate, this might indeed explain the observed excess of antihyperons.

Before, let us first remind that nonequilibrium inelastic hadronic reactions can explain to a very good extent the overall strangeness production seen experimentally . The major amount of produced strange particles (kaons, antikaons and $`\mathrm{\Lambda }`$’s) at SPS-energies can be understood in terms of early and energetic primary, secondary and ternary non-equilibrium interactions. Strange quarks are thought to be produced initially via string excitations. However, applying the usual concept of binary collisions within the transport approaches, it had been the original claim  that chemical equilibrium values for the abundancies of antihyperons can not be explained by successive binary (strangeness exchange) reactions. Only by incorporating more exotic descriptions like color-rope formation  or high-dense droplet formation  in a hadronic transport model a more abundant production of antihyperons can be obtained.

If the hadronic degrees of freedom are in a state of thermodynamical equilibrium, which constitutes the basic concept of thermal model analyses, a dynamically realization to describe such a system has to fulfill the concept of detailed balance in the considered ‘chemical’ reactions. Hence, multi-mesonic (fusion-like) ‘back-reactions’ involving $`n_1`$ pions and $`n_2`$ kaons of the type

$$n_1\pi +n_2K\overline{Y}+N,$$
(1)
corresponding to the inverse of the strong binary baryon-antihyperon annihilation process (similar to the standard baryon annihilation $`\overline{p}+pn\pi `$), have, in principle, to be taken care of in a dynamical simulation. In (1) $`n_2`$ counts the number of anti-strange quarks within the antihyperon $`\overline{Y}`$. $`n_1+n_2`$ is expected to be around $`57`$.

Furthermore, the annihilation rate $`\mathrm{\Gamma }_{\overline{Y}}`$ is expected to be large: It is plausible to assume that the annihilation cross sections are approximately similar in magnitude as for $`p+\overline{p}`$ at the same relative momenta. Hence, with $`\sigma _{p\overline{Y}n_1\pi +n_2K}50`$ mb one has $`(\mathrm{\Gamma }_{\overline{Y}})^11/(\sigma _{\overline{Y}N}v_{\overline{Y}N}\rho _N)13fm/c`$, when adopting for the baryon density $`\rho _N12\rho _0`$. The latter value seems reasonable: In fig. 3 we have depicted the (net) baryon density as a function of time at a space region for particles at midrapidity obtained within the microscopic transport model HSD as employed in . The figure illustrates that a pure hadronic fireball (without any string-like excitations) at a (net) baryon density $`\rho _B2\rho _0`$ has started to be established after approximately 5 fm/c after the onset of the collision.

The annihilation together with the multi-mesonic reactions do effectively lead to the following master equation for the yield of a specific antihyperon specie

$$\frac{d}{dt}N_{\overline{Y}}=\mathrm{\Gamma }_{\overline{Y}}\left\{N_{\overline{Y}}N_{\overline{Y}}^{eq.}\underset{n_1}{}p_{n_1}\left(\frac{N_\pi }{N_\pi ^{eq.}}\right)^{n_1}\left(\frac{N_K}{N_K^{eq.}}\right)^{n_2}\right\}.$$
(2)
Here $`p_{n_1}`$ states the (unknown) relative probability of the reaction (1) to decay into a specific number $`n_1`$ of pions with $`_{n_1}p_{n_1}=1`$. This master equation can be derived from an underlying Boltzmann-type transport equation . The second term describes the gain per unit time in the hyperon yield due to several coalescing pions and kaons. (The microscopic picture of this complicated process is completely hidden; the term is dictated by the principle of detailed balance.) Assuming further that the pions, baryons and kaons do stay in thermal and chemical equilibrium, the master equation becomes simply

$$\frac{d}{dt}N_{\overline{Y}}=\mathrm{\Gamma }_{\overline{Y}}\left\{N_{\overline{Y}}N_{\overline{Y}}^{eq}\right\}.$$
(3)
The annihilation rate thus dictates the timescale for (local) chemical equilibration of the antihyperons together with the pions, nucleons and kaons. This scale competes with the expansion timescale of the late hadronic fireball, which is in the same range or larger. In any case these multi-mesonic, hadronic reactions can explain a sufficiently fast equilibration before or at the (so called) chemical freeze-out (with parameters given by the thermal model analyses, where the (net) baryon density has dropped to around 0.5 - 1 $`\rho _0`$ ). This interpretation does rest on the (conservative) view that before the chemical freeze-out already a hadronic system has established.

To be more quantitative we show in fig. 3 the number of one particular specie of antihyperons (normalized to the net baryon number) as a function of time when solving the above master equation. The initial abundancy is set to zero. The pions, kaons and nucleons are assumed to stay in thermodynamical equilibrium throughout a (late-stage) isentropic expansion with parameters characteristic for SPS energies. The net baryon density $`\rho _B(t)`$ as a function of time has been effectively parametrized from fig. 3. Further details and scenarios will be given in . In addition, the instantaneous equilibrium abundancy of antihyperons, $`N_{\overline{Y}}^{eq.}(T(t),\mu _B(t),\mu _s(t))/N_B`$, is also depicted. This equilibrium number strongly decreases as a function of time (or decreasing temperature). At a time of 8-10 fm/c (3-5 fm/c after the onset of the description) the antihyperons effectively do ‘freeze’ out roughly at their particular equilibrium value (at a temperature around $`150160`$ MeV) at that given time interval. Refering to fig. 3 this happens at a baryon density of around 0.5 - 1 $`\rho _0`$ and which coincides with the parameters of the thermal model analyses. Beyond that ‘point’ at moderately low baryon densities the multi-mesonic creation processes becomes more and more ineffective because of the rapid expansion. This might then also explain the ‘position’ of the chemical freeze-out point for the antihyperons . In this respect it would also be very interesting to adopt the reasoning in the recent work by Rapp and Shuryak  (which, in part, had actually triggered our work), that the pions and also now the kaons might in fact go out of perfect chemical equilibrium below chemical freeze-out, as their total abundancies effectively have to stay constant.

To summarize, multi-mesonic production of antihyperons is a consequence of detailed balance and, as the rate $`\mathrm{\Gamma }_{\overline{Y}}`$ is indeed very large, seems to be the by far most dominant source in a hadronic gas. This is a remarkable observation, as it clearly demonstrates the importance of hadronic multi-particle channels. At the moment such ‘back-reactions’ are not included in the present transport codes and some new strategy has to be invented  in order to be more competitive for a direct comparison with various experimental findings. Also, our argument should perfectly apply at AGS energies. Indeed a large ratio of $`\overline{\mathrm{\Lambda }}/\overline{p}23`$ for some central rapidity window has been observed by various collaborations, which might indeed also favor a scenario of nearly chemically saturated strange and nonstrange antibaryon populations. More precise experimental data as well as a theoretical analysis along the line of thermal models would be most welcome. Measurements of antihyperon production could also be done at possible future heavy-ion facilities at GSI. This would be a very interesting opportunity to unreveal the mechanism here proposed.
Acknowledgments
This work has been done in part with S. Leupold, with whom the idea of multi-mesonic antihyperon production has been developed. The author is also indebted to W. Cassing for providing fig. 3 and for detailed and valuable comments. He also wants to thank U. Mosel for his continuing interest in this particular subject on the (possible) importance of multi-particle interactions in heavy-ion collisions.


text:
# Centrality Dependence of Multiplicity, Transverse Energy, and Elliptic Flow from Hydrodynamics

## I Introduction

Elliptic flow  is a collective flow pattern which develops in non-central relativistic heavy-ion collisions as a result of the spatial deformation of the initial transverse overlap area. It requires rescattering among the produced particles as a mechanism to map the initial spatial deformation of the reaction zone onto the finally observed hadron momentum distributions. It is quantified by the second harmonic coefficient $`v_2`$ of an azimuthal Fourier decomposition of the measured spectrum $`dN/(dyp_\mathrm{T}dp_\mathrm{T}d\varphi )`$ . Its magnitude $`v_2`$ and its shape $`v_2(p_\mathrm{T})`$ as a function of the hadron transverse momentum are sensitive to the scattering rate among the produced secondaries, especially during the dense early stage of the expansion . The largest elliptic flow signal, especially at high $`p_\mathrm{T}`$, arises in hydrodynamic simulations  which assume local thermal equilibrium at every space-time point, i.e. essentially instantaneous thermalization or infinite scattering rate. Surprisingly, such hydrodynamic simulations  are in very good agreement with first results from $`\sqrt{s}=130A`$ GeV Au+Au collisions at the Relativistic Heavy Ion Collider (RHIC) , up to transverse momenta of $`1.52`$ GeV/$`c`$.

It was recently suggested  that a combined analysis of the full set of hadronic single particle spectra and their elliptic flow as a function of collision centrality should allow to outline the domain of applicability of the hydrodynamic approach. Inside this domain such an analysis would constrain the initial baryon and energy density and the final freeze-out temperature sufficiently well to become sensitive to details of the equation of state (EOS) of the fireball matter . In this context $`v_2`$ provides access to the EOS during the early expansion stage  which, at RHIC energies, is presumably in the quark-gluon plasma (QGP) phase. Direct verification of the phase transition between QGP at high energy density and a hadron resonance gas at lower energy density should then become possible by accurately measuring the excitation function of radial and elliptic flow .

In the hydrodynamic limit the EOS affects the elliptic flow signal through the velocity of sound, $`c_s=\sqrt{dP/de}`$ . However, the sensitivity of $`v_2`$ on $`c_s`$ is not very strong, and even a first order phase transition, where $`c_s`$ vanishes in the mixed phase, affects $`v_2`$ only on the 20$`\%`$ level . This makes accurate measurements and systematic theoretical checks indispensable. One possible source of ambiguity, which has not been systematically investigated in previous studies , is the sensitivity of the radial and elliptic flow pattern on the shape of the initial transverse density profile. The latter depends on the scaling of secondary particle production with the number of colliding nucleons which itself is controlled by the collision centrality. It was recently found  that the momentum-space anisotropy at high $`p_\mathrm{T}`$, where it is not of hydrodynamic origin but due to quark energy loss , is quite sensitive to the initial density profile. A good hydrodynamic baseline for the dependence of $`v_2(p_\mathrm{T})`$ on the initial profiles may thus help to clarify the physics controlling the breakdown of the hydrodynamic model and the transition to the hard physics domain at high $`p_\mathrm{T}`$.

In the present work, we investigate five options which are expected to span the realistic range of possibilities. In the first four, either the initial energy or the initial entropy density are assumed to scale with either the number of wounded nucleons (“soft” or “non-perturbative” particle production) or the number of binary nucleon-nucleon collisions (“hard” or “perturbative” particle production). In the fifth model perturbative particle production is modified by implementing gluon shadowing  in the initial state and by limiting the growth of the production cross section by gluon saturation in the final state (“saturation model” ). This brings in some non-perturbative elements as well. In the first four parametrizations we normalize the initial energy density profile such that for central collisions we reproduce within errors the total charged multiplicity density at midrapidity, $`dN_{\mathrm{ch}}/d\eta |_{|\eta |<1}=\mathrm{\hspace{0.17em}555}\pm 10\%`$, as measured by PHOBOS for Au+Au at $`\sqrt{s}=\mathrm{\hspace{0.17em}130}A`$ GeV . On the other hand, in the saturation model  the initial energy density is fixed by assuming that the transverse energy of produced minijets is entirely converted to thermalized energy density. The different scaling laws implied by the models then translate into different centrality dependencies of $`dN_{\mathrm{ch}}/d\eta `$ which can be tested against the new data presented recently .

The hydrodynamic model describes an adiabatic evolution from one local equilibrium state to another. Our version of the model assumes longitudinal boost invariance which implies conservation of the entropy rapidity density $`dS/dy`$. Using the relation between entropy and particle multiplicity in a thermalized system, the measured centrality dependence of the final multiplicity density $`dN/dy`$ can thus be mapped onto the centrality dependence of the initial parton multiplicity density at the point of thermalization. This allows to constrain models for the initial production of secondary particles, under the assumption of subsequent adiabatic evolution.

We show that different models for the initial energy and entropy production lead to different radial and elliptic flow patterns. At fixed impact parameter, these result in different $`p_\mathrm{T}`$-dependences of the spectra and of the differential elliptic flow $`v_2(p_\mathrm{T})`$. However, the $`p_\mathrm{T}`$-integrated elliptic flow $`v_2`$ as a function of the final charged multiplicity density $`dN_{\mathrm{ch}}/dy`$ and the differential elliptic flow $`v_2(p_\mathrm{T})`$ for minimum bias events show surprisingly little sensitivity to the model used for initializing the hydrodynamic evolution. These two observables have been used to argue for the applicability of the hydrodynamic model at RHIC energies , and they were shown to be sensitive to the EOS of the expanding matter . Their insensitivity to the shape of the initial transverse energy density profile may prove to be crucial for the process of extracting reliable information on the EOS from a detailed analysis of the measured collective flow patterns.

## II Models for the initial transverse energy density profile

### A Soft particle production

In Pb+Pb collisions at the SPS, the rapidity densities at midrapidity of both the total produced transverse energy,***We use $`E_\mathrm{T}=_iE_i\frac{p_{\mathrm{T},i}}{|𝒑_i|}`$ where the sum is over all particles. $`dE_\mathrm{T}/dy`$, and of the charged multiplicity, $`dN_{\mathrm{ch}}/dy`$, scale approximately linearly with the number of participating nucleons, $`N_{\mathrm{part}}`$ , except for very peripheral collisions with $`N_{\mathrm{part}}\stackrel{<}{}100`$. Similar earlier observations in smaller collision systems and/or at lower energies have led to the notion that particle and transverse energy production can be described by the “Wounded Nucleon Model”  in which each participating nucleon contributes to the multiplicity and transverse energy only once in its first collision and not every time it suffers further collisions with other projectile or target nucleons. The validity of this model requires destructive interference effects between subsequent nucleon-nucleon collisions which are thought to be characteristic of non-perturbative or “soft” particle production processes. For two nuclei colliding at impact parameter $`𝒃=b𝒆_x`$, the density of wounded nucleons in the transverse plane (parametrized by $`𝒔=(x,y)`$) can be calculated from the simple geometric formula (Glauber ansatz, see  and references therein):

$`n_{\mathrm{WN}}(𝒔;𝒃)`$ $`=`$ $`T_A\left(𝒔+\frac{1}{2}𝒃\right)\left[1\left(1{\displaystyle \frac{\sigma T_B\left(𝒔\frac{1}{2}𝒃\right)}{B}}\right)^B\right]`$ (1)
$`+`$ $`T_B\left(𝒔\frac{1}{2}𝒃\right)\left[1\left(1{\displaystyle \frac{\sigma T_A(𝒔+\frac{1}{2}𝒃)}{A}}\right)^A\right].`$ (2)
$`T_A`$ is the nuclear thickness function of nucleus $`A`$,

$$T_A(𝒔)=_{\mathrm{}}^+\mathrm{}𝑑z\rho _A(𝒔,z),$$
(3)
with the density $`\rho _A`$ given by a Woods-Saxon profile,

$$\rho _A(𝒓)=\frac{\rho _0}{1+\mathrm{exp}[(rR_0)/\xi ]},$$
(4)
and similarly for nucleus $`B`$. For the Woods-Saxon profile we take standard parameter values, $`R_0=\mathrm{\hspace{0.17em}1.12}A^{1/3}0.86A^{1/3}`$ fm for the radius and $`\xi =\mathrm{\hspace{0.17em}0.54}`$ fm for the surface diffuseness. The integral over (1) gives the number of “wounded” or “participating” nucleons, $`N_{\mathrm{part}}N_{\mathrm{WN}}`$, as a function of impact parameter (see Fig. 1). For the total nucleon-nucleon cross section at $`\sqrt{s}=\mathrm{\hspace{0.17em}130}A`$ GeV we take $`\sigma =\mathrm{\hspace{0.17em}40}`$ mb.

#### 1 Model sWN

Hydrodynamics with boost-invariant longitudinal expansion conserves the entropy per unit rapidity $`dS/dy`$. At fixed freeze-out temperature and chemical potential, $`dS/dy`$ is related one-to-one to the measured charged multiplicity density, $`dN_{\mathrm{ch}}/dy`$. If the latter scales linearly with the number of wounded nucleons, $`N_{\mathrm{WN}}(b)=d^2sn_{\mathrm{WN}}(𝒔;𝒃)`$, barring a strong $`b`$-dependence of the freeze-out conditions it is natural to assume that the initial entropy density in the transverse plane is proportional to that of the wounded nucleons :

$`s(𝒔;\tau _0;𝒃)=K_s(\tau _0)n_{\mathrm{WN}}(𝒔;𝒃).`$ (5)
At SPS and RHIC energies the net baryon density is so small that its influence on the pressure and thus on the developing flow pattern is hardly noticeable. For simplicity we can thus also parametrize the initial net baryon density $`n`$ as being proportional to $`n_{\mathrm{WN}}`$:

$`n(𝒔;\tau _0;𝒃)=K_n(\tau _0)n_{\mathrm{WN}}(𝒔;𝒃).`$ (6)
Equations (5) and (6) together define the first model for the initialization, sWN, to be studied below.

The implementation of the initialization (5) into the hydrodynamic computer code is complicated by the fact that in our formulation the hydrodynamic equations propagate the fields $`n(x)`$ and $`e(x)`$ (and not $`s(x)`$). The initialization (5) thus requires the additional step of computing $`e(𝒔;\tau _0;𝒃)`$ from $`s(𝒔;\tau _0;𝒃)`$ by solving an implicit equation resulting from the thermodynamic identity

$`T(e,n)s=e+p(e,n)\mu (e,n)n,`$ (7)
where $`T(e,n)`$, $`p(e,n)`$, and $`\mu (e,n)`$ are tabulated values for the temperature, pressure and baryon chemical potential corresponding to the selected equation of state.

#### 2 Model eWN

In a hydrodynamic approach with fixed freeze-out conditions, it is strictly speaking not possible that both the transverse energy $`dE_\mathrm{T}/dy`$ and the charged multiplicity $`dN_{\mathrm{ch}}/dy`$ scale linearly with $`N_{\mathrm{WN}}`$. During the hydrodynamic evolution thermal energy is converted into longitudinal and transverse flow energy, and the time available for this conversion increases with the number of wounded nucleons. Since the transverse flow increases the average transverse energy per particle, $`dE_\mathrm{T}/dy`$ should rise more quickly with $`N_{\mathrm{WN}}`$ than $`dN_{\mathrm{ch}}/dy`$. There is some indication for this to happen in the WA98 data from Pb+Pb collisions at the SPS , where $`E_\mathrm{T}/N_{\mathrm{ch}}`$ at midrapidity rises slightly from very peripheral to semiperipheral collisions and then saturates from semiperipheral to central collisions. However, the experimental effect is small and not clearly statistically significant.

Given this unclear experimental situation, and because it simplifies the initialization process, we assumed in our previous work  that it is the initial energy density (and not the entropy density) which scales with the density of wounded nucleons in the transverse plane:

$`e(𝒔;\tau _0;𝒃)=K_e(\tau _0)n_{\mathrm{WN}}(𝒔;𝒃).`$ (8)
Equation (8) together with Eq. (6) defines our second initialization model, eWN. The effect of flow on $`E_\mathrm{T}/N_{\mathrm{ch}}`$ at RHIC and its dependence on the initialization of the hydrodynamic evolution will be discussed in Section III.

### B Hard particle production

At higher and higher collision energies, one expects that hard collisions among quarks and gluons from the nuclear structure functions become more and more important and eventually dominate secondary particle production . In this limit secondary particle production is a result of incoherent parton-parton collisions, and each nucleon-nucleon collision contributes equally to the cross section. Particle and energy production should then be related to the distribution of the number of binary nucleon-nucleon collisions in the transverse plane. It is given in terms of the nuclear thickness functions (3) by

$`{\displaystyle \frac{dN_{AB}^{\mathrm{coll}}(𝒔;𝒃)}{d^2s}}=\sigma T_A\left(𝒔+\frac{1}{2}𝒃\right)T_B\left(𝒔\frac{1}{2}𝒃\right),`$ (9)
where $`\sigma `$ is the nucleon-nucleon cross section. The integral of (9) is the nuclear overlap function, $`T_{AB}(𝒃)=d^2sT_A\left(𝒔+\frac{1}{2}𝒃\right)T_B\left(𝒔\frac{1}{2}𝒃\right)`$, which is normalized to $`AB`$: $`d^2bT_{AB}(𝒃)=AB`$. It gives the total number of binary collisions as a function of impact parameter, $`N_{\mathrm{coll}}(b)`$, which is also shown in Fig. 1.

#### 1 Model sBC

If the system of secondary particles thermalizes quickly via elastic collisions, their number density defines, up to a proportionality constant, the initial entropy density at the beginning of the hydrodynamic expansion:

$`s(𝒔;\tau _0;𝒃)=\stackrel{~}{K}_s(\tau _0)T_A\left(𝒔+\frac{1}{2}𝒃\right)T_B\left(𝒔\frac{1}{2}𝒃\right).`$ (10)
Assuming that the initial net baryon density in the transverse plane can also be calculated perturbatively from the nuclear structure functions , we write

$`n(𝒔;\tau _0;𝒃)=\stackrel{~}{K}_n(\tau _0)T_A\left(𝒔+\frac{1}{2}𝒃\right)T_B\left(𝒔\frac{1}{2}𝒃\right).`$ (11)
As discussed in Section II A 1, the results do not depend on whether we use (11) or (6). The combination of Eqs. (10) and (11) defines model sBC for the initialization.

#### 2 Model eBC

Within the perturbative approach to particle production one can also argue that each nucleon-nucleon collision contributes with equal probability not only to the number of produced secondaries, but also to the energy carried by them. This leads to the ansatz that the initial energy (and not the entropy) density is proportional to the density of binary collisions in the transverse plane:

$`e(𝒔;\tau _0;𝒃)=\stackrel{~}{K}_e(\tau _0)T_A\left(𝒔+\frac{1}{2}𝒃\right)T_B\left(𝒔\frac{1}{2}𝒃\right).`$ (12)
This equation together with Eq. (11) defines model eBC for the initialization.

### C The saturation model

In the saturation model  the production of partons becomes inhibited below a saturation scale $`p_{\mathrm{sat}}`$, determined as the transverse momentum scale where the produced partons start to overlap transversally. The formation time of the QGP at each point $`𝒔`$ is given by $`\tau _{\mathrm{sat}}=1/p_{\mathrm{sat}}(\sqrt{s},A,B,𝒔;𝒃)`$. The local initial energy density profile can then be worked out at the central slice as in ,

$`e(𝒔;𝒃)`$ $`=`$ $`{\displaystyle \frac{dE_\mathrm{T}^{\mathrm{pQCD}}}{d^2sdz}}`$ (13)
$`=`$ $`T_A\left(𝒔+\frac{1}{2}𝒃\right)T_B\left(𝒔\frac{1}{2}𝒃\right)`$ (15)
$`\times \sigma E_\mathrm{T}(\sqrt{s},p_{\mathrm{sat}},\mathrm{\Delta }y,A,B,𝒔;𝒃){\displaystyle \frac{1}{\tau _{\mathrm{sat}}\mathrm{\Delta }y}},`$
where $`dz\tau _{\mathrm{sat}}\mathrm{\Delta }y`$ within the central rapidity unit $`\mathrm{\Delta }y`$. The first $`E_\mathrm{T}`$-moment of the minijet distribution, $`\sigma E_\mathrm{T}`$ , is computed in lowest order pQCD as the first $`p_\mathrm{T}`$-moment of the distributions of minijets with $`p_\mathrm{T}p_{\mathrm{sat}}`$ and $`y`$ within $`\mathrm{\Delta }y`$. The EKS98 shadowing  of the parton distributions is included and the NLO contributions to $`\sigma E_\mathrm{T}`$  are simulated by an approximate factor $`K=\mathrm{\hspace{0.17em}2}`$, as also done in .

In  a fully saturated system was considered by extending the computation down to very low values of $`p_{\mathrm{sat}}=\mathrm{\hspace{0.17em}0.5}`$ GeV, and neglecting the tails of the number (energy) density distributions at large transverse distances. Now, for the hydrodynamic description, we have to consider the tails also. To maintain the spirit of the local saturation model of , we restrict the saturation to the regions where $`p_{\mathrm{sat}}0.75`$ GeV and compute the tail of the energy density profile from Eq. (13) with $`p_{\mathrm{sat}}=0.75`$ GeV. Through this procedure, the multiplicity of  in central collisions is recovered (if computed as in ). For central collisions, the tail contributes only $`7\%`$ to the total multiplicity. We would like to emphasize that we have not tried to fit the saturation model to the centrality data but, rather, to keep the approach as close to the orginal local saturation idea  as possible.

As discussed in , from the energy per particle point of view the system looks thermal already at saturation. The same can be shown to hold also in the locally saturated system . Within the saturation framework it is therefore meaningful to switch on the hydrodynamic evolution already at $`\tau _{\mathrm{sat}}`$ which now is a local quantity.

For the hydrodynamic set-up used in this study, however, the density profiles at a constant initial time are required. To circumvent this problem, we consider two possibilities: first, for each impact parameter $`𝒃`$, we simply use the earliest time, i.e. $`\tau (𝒔=\mathrm{\hspace{0.17em}0})`$, as the initial time. These are shown in Fig. 2. Notice that the maximum initial time is $`(0.75\mathrm{GeV})^1=\mathrm{\hspace{0.17em}0.26}`$ fm/$`c`$, and that a constant initial time is a good approximation for the central region and for small impact parameters. Alternatively, the energy densities computed from Eq. (13) can be evolved at each $`𝒔`$ and at each $`\tau _{\mathrm{sat}}(𝒔;𝒃)`$ to $`\tau _i=\mathrm{\hspace{0.17em}0.26}`$ fm/$`c`$ assuming boost-invariant flow without transverse expansion:We thank P.V. Ruuskanen for helpful discussions related to this point. $`e(𝒔;\tau _i;𝒃)=e(𝒔;\tau _{\mathrm{sat}}(𝒔;𝒃);𝒃)\times \left(\tau _{\mathrm{sat}}(𝒔;𝒃)/\tau _i\right)^{4/3}`$. We have checked, however, that the latter procedure leads only to a few percent increase in the central multiplicities relative to the former one. Therefore, we take in the following the initial time for the saturation model from Fig. 2 and use $`e(𝒔;\tau _{\mathrm{sat}}(b);𝒃)`$ as the initial profile for the transversally expanding hydrodynamics, as discussed next.

### D Initial energy density and spatial anisotropy

Figure 3 shows the initial energy density profiles for the different initialization models. ($`x`$ is the direction of the impact parameter $`𝒃`$, while $`y`$ points perpendicular to the reaction plane.) For models sWN, eWN, sBC and eBC the profiles are normalized such that in each case the total charged multiplicity per unit pseudorapidity at $`b=\mathrm{\hspace{0.17em}0}`$ is $`dN_{\mathrm{ch}}/dy=\mathrm{\hspace{0.17em}670}`$ at $`y=\mathrm{\hspace{0.17em}0}`$. For model eWN  this corresponds to $`dN_{\mathrm{ch}}/d\eta |_{|\eta |<1}=\mathrm{\hspace{0.17em}545}`$ for the $`6\%`$ most central collisions (a bit less for the three other models). Within errors, this is consistent with the first published PHOBOS data , but slightly below the more accurate recent data from PHENIX and PHOBOS . To preserve consistency with our previous publications , we decided against increasing our normalization of $`dN_{\mathrm{ch}}/d\eta `$ to the new data since this would have implied retuning the initial conditions and freeze-out temperature in order to keep the spectra and elliptic flow unchanged.

For the comparison with the other models in Fig. 3, the profile obtained from the saturation model was evolved from the formation time $`\tau _{\mathrm{sat}}(b)`$ (see Fig. 2) to $`\tau _0=\mathrm{\hspace{0.17em}0.6}`$ fm/$`c`$, the time at which the hydrodynamic evolution was started for the other four initializations. Within the line widths of Fig. 3, it did not matter whether this scaling was done by assuming only boost-invariant longitudinal expansion, $`e(𝒔,\tau _0)=e(𝒔,\tau _{\mathrm{sat}})\times (\tau _{\mathrm{sat}}/\tau _0)^{4/3}`$, or by including also the transverse expansion.

Since for a thermalized parton gas at nearly vanishing net baryon density $`es^{4/3}`$, Eq. (10) gives distributions which are more sharply peaked at the origin than Eq. (12). We will see that this results in a different centrality dependence of the total entropy per unit rapidity $`dS/dy`$ (which, for boost-invariant longitudinal expansion, is proportional to the integral of $`s(𝒔)`$ over the transverse plane). Model eBC interpolates between model sWN, where $`dS/dyN_{\mathrm{part}}`$, and model sBC, where $`dS/dyN_{\mathrm{coll}}N_{\mathrm{part}}^{4/3}`$ (the latter proportionality was checked numerically to hold with excellent accuracy over the entire impact parameter range). This is similar to the HIJING model  where the final charged particle rapidity density is a linear superposition of a soft component, scaling with $`N_{\mathrm{part}}`$, and a hard component which scales with $`N_{\mathrm{coll}}N_{\mathrm{part}}^{4/3}`$. In contrast to the models here, however, HIJING has no rescattering among the produced particles and thus no collective flow. Finally, the saturation model is seen to be close to model eWN for central collisions, while for semiperipheral collisions it nearly coincides with model sWN.

Figure 4 shows the initial spatial anisotropy

$$ϵ_x(b)\delta (b)=\frac{y^2x^2}{y^2+x^2}$$
(16)
as a function of impact parameter, evaluated with the initial transverse energy density as weight function, for the five initialization models studied here. At fixed impact parameter, $`ϵ_x`$ varies by up to $`20\%`$. Since the hydrodynamic evolution maps the initial spatial anisotropy onto the final momentum-space, one expects from Fig. 4 that the impact parameter dependence of the elliptic flow $`v_2`$ should show a similar sensitivity to the initialization. We will see that, except for the saturation model, this model sensitivity is almost exactly cancelled by the corresponding variations in the impact parameter dependence of the produced charged particle multiplicity (i.e. entropy).

## III Centrality dependence of multiplicity and transverse energy

Using the initial transverse energy and baryon density profiles from the previous section and the numerical code described in , we solve the hydrodynamic equations for the transverse evolution of the reaction zone, assuming boost-invariant expansion in the longitudinal direction. We study only Au+Au collisions at $`\sqrt{s}=\mathrm{\hspace{0.17em}130}A`$ GeV. For the models eWN, sWN, eBC and sBC we use the same initial time $`\tau _0=\mathrm{\hspace{0.17em}0.6}`$ fm/$`c`$ as in  to start the hydrodynamic expansion. For the saturation model the hydrodynamic expansion is started at $`\tau _{\mathrm{sat}}(b)`$ shown in Fig. 2. All calculations reported here are performed with EOS Q, an equation of state with a first order phase transition from a hadron resonance gas to a non-interacting quark-gluon gas at $`T_\mathrm{c}=\mathrm{\hspace{0.17em}164}`$ MeV , and freeze-out at $`T_\mathrm{f}=\mathrm{\hspace{0.17em}120}`$ MeV. The sensitivity of the spectra and elliptic flow to the EOS and freeze-out temperature were studied in .

### A Charged particle multiplicity

Figure 5 shows the final charged particle rapidity density per participating nucleon pair resulting from the different initializations. In the bottom panel the rapidity density at midrapidity is converted to pseudorapidity density by the transformation

$`{\displaystyle \frac{dN_{\mathrm{ch}}}{d\eta }}={\displaystyle \underset{i\mathrm{charged}}{}}{\displaystyle p_\mathrm{T}𝑑p_\mathrm{T}\frac{p_\mathrm{T}}{\sqrt{m_i^2+𝒑^2}}\frac{dN_i}{dyp_\mathrm{T}dp_\mathrm{T}}},`$ (17)
setting $`y=p_\mathrm{L}=\mathrm{\hspace{0.17em}0}`$. At fixed freeze-out temperature $`dN_{\mathrm{ch}}/dy`$ is a direct measure for the entropy density $`dS/dy`$, and the approximate constancy of the curve for model sWN thus reflects the approximate conservation of entropy by the hydrodynamic evolution. In fact, the slight increase of $`(dN_{\mathrm{ch}}/dy)/N_{\mathrm{part}}`$ for small values of $`N_{\mathrm{part}}`$ (Fig. 5, upper panel) can be traced back to a small amount of entropy production by deflagration shocks which arise during the hydrodynamic expansion stage as a consequence of the first order phase transition  once the initial energy density in the center of the reaction zone increases above the critical value for QGP formation.

Model eWN (solid line in Fig. 5) was used in Refs. , while model sWN was employed in . Fig. 5 shows that both initialization models are disfavoured by the data, model eWN more so than model sWN. The saturation model produces a rather similar centrality dependence of the charged multiplicity as model eWN.

As expected, the saturation model results are close to those obtained in Ref. . The careful reader notices, however, that the values $`(dN_{\mathrm{ch}}/dy)/(0.5N_{\mathrm{part}})`$ for the saturation model shown in Fig. 5 are $`𝒪(10\%)`$ lower than those shown in Fig. 4 of . This is due to a number of partially cancelling effects : (i) Inclusion of the quark degrees of freedom increases the initial entropy by about 20% compared to . (ii) In a realistic hadron resonance gas containing also heavy particles, at a freeze-out temperature of $`T_\mathrm{f}=\mathrm{\hspace{0.17em}120}`$ MeV each particle carries on average 4.87 units of entropy, instead of the 4 units assumed in  (which is a good approximation in the ultrarelativistic limit). This decreases the total multiplicity from that in . (iii) Finally, the decay of unstable resonances changes the fraction of charged particles from 2/3  to $`0.62`$ of the total multiplicity. The net reduction effect for the final charged multiplicity rapidity density is about 10% as mentioned.

The best representation of the shape of the data, even down to rather peripheral collisions (low $`N_{\mathrm{part}}`$-values), seems to be given by model eBC, once appropriately renormalized to the new, slightly higher multiplicity densities in central collisions. Alternatively, one may try a linear combination of the sWM and sBC parametrizations, with a large sWN and a smaller sBC contribution  (the latter is apparently absent at the SPS where the multiplicity scales linearly with $`N_{\mathrm{part}}`$). As discussed in Section II C, the main difference between the model eBC and the saturation model is that the $`𝒔`$-dependence in eBC is contained only in the product $`T_A\left(𝒔+\frac{1}{2}𝒃\right)T_B\left(𝒔\frac{1}{2}𝒃\right)`$ whereas in the saturation model  also the cross sections (and their moments) depend on $`𝒔`$ through the local saturation scale $`p_{\mathrm{sat}}(𝒔)`$. The data do not appear to support this particular implementation  of gluon saturation at the present RHIC energy of $`\sqrt{s}=\mathrm{\hspace{0.17em}130}A`$ GeV for the non-central collisions. In central collisions, however, parton saturation is not excluded as the dominant mechanism. We conclude that obviously the saturation region is now extended too far in transverse direction and that the eBC tails in the initial energy density profile (which were not included at all in ) should be given more weight than they have now. On the other hand, the data strongly indicate that the initial energy deposition process does involve a component which scales with the number of binary nucleon-nucleon collisions. We refrain from a discussion whether or not this indeed proves the onset of “hard” perturbative physics at RHIC .

### B Transverse energy per particle

Hydrodynamic flow is a result of the conversion of a fraction of the thermal energy into collective flow kinetic energy, through work done by the thermodynamic pressure. As a consequence, the thermal energy per particle decreases as a function of time. If the reaction zone were to undergo boost-invariant expansion only along the longitudinal direction, $`dS/dydN_{\mathrm{ch}}/dy`$ would be constant and $`dE_\mathrm{T}/dy`$ would be proportional to the thermal energy density. The reduced thermal energy per particle would thus be directly reflected in the finally observed transverse energy per charged particle, $`(dE_\mathrm{T}/dy)/(dN_{\mathrm{ch}}/dy)`$ . Since in more central collisions the thermodynamic pressure does longitudinal hydrodynamic work for a longer time, the finally observed ratio $`(dE_\mathrm{T}/dy)/(dN_{\mathrm{ch}}/dy)`$ should increase much more slowly with $`N_{\mathrm{part}}`$ than the initial energy per particle  established during the energy deposition process. In reality, the system undergoes not only longitudinal, but also transverse expansion. The transverse collective flow adds a kinetic contribution to $`E_\mathrm{T}`$ which reduces the diluting effect on $`E_\mathrm{T}/N_{\mathrm{ch}}`$ from the longitudinal expansion.

This shows that the centrality dependence of $`E_\mathrm{T}/N_{\mathrm{ch}}`$ is an interesting observable which reflects the interplay of the centrality dependences of initial transverse energy and particle production, longitudinal work done by the hydrodynamic pressure, and radial flow. In Figure 6 we show the corresponding results from the hydrodynamic model for Au+Au collisions at RHIC ($`\sqrt{s}=\mathrm{\hspace{0.17em}130}A`$ GeV) for the five different initialization models. For comparison we also show the centrality dependence for Pb+Pb collisions at the SPS, for model eWN with parameters adjusted to the measured spectra . Comparing the eWN curves for SPS and RHIC, one sees that the larger radial flow at RHIC has very little effect on the average $`E_\mathrm{T}`$ per particle in the final state: the gain in transverse kinetic flow energy is largely compensated by more longitudinal work done at RHIC (smaller $`\tau _0`$, somewhat larger $`\tau _\mathrm{f}`$). The eWN curve at the SPS is consistent in shape with the WA98 data shown in Fig. 14 of , but at the upper edge of their systematic error band.

At RHIC, the models eWN and sWN as well as the saturation model all give similarly flat centrality dependences for the average $`E_\mathrm{T}`$ per charged particle; they differ slightly in normalization and in how steeply this ratio drops for very peripheral collisions. The data in Figure 5 disfavour these initialization models. On the other hand, Model eBC, which works better in Figure 5, exhibits a markedly stronger centrality dependence of $`E_\mathrm{T}/N_{\mathrm{ch}}`$, with a clearly visible slope below $`N_{\mathrm{part}}\mathrm{\hspace{0.17em}250}`$. (This is even more true for the initialization sBC, but for that model the corresponding curve in Fig. 5 is too steep.) This reflects a similar strong centrality dependence of $`(dE_\mathrm{T}/dy)/(dS/dy)`$ in the initial state, which is only partially offset by the increasing amount of longitudinal hydrodynamic work done as the collisions become more central.

As recently pointed out in , a measurement of the ratio plotted in Fig. 6 will be very helpful in complementing the information contained in Fig. 5 and checking for the effects of longitudinal work. We note that non-equilibrium expansion scenarios like those studied in , where the matter performs less longitudinal work than predicted hydrodynamically, should lead to an even stronger dependence of $`E_\mathrm{T}/N_{\mathrm{ch}}`$ on $`N_{\mathrm{part}}`$ than that shown in Fig. 6. In contrast, recent data from the PHENIX Collaboration  exhibit the opposite tendency: they are flatter than the eBC curve in Fig. 6 and, in fact, quantitatively consistent in magnitude and shape with the WA98 data at the SPS . It remains to be seen whether e.g. a linear combination of a large sWN with a small sBC component (which can be made to agree with the multiplicity data in Fig. 5 ) leads to a sufficiently flat behaviour in Fig. 6 to be consistent with the PHENIX data , and whether the magnitude of $`E_\mathrm{T}/N_{\mathrm{ch}}`$ can be lowered by replacing the unrealistic chemical composition from the hydrodynamical model by a more realistic one (our hydrodynamic model assumes that chemical equilibrium is preserved until kinetic freeze-out at $`T_\mathrm{f}=\mathrm{\hspace{0.17em}120}`$ MeV whereas preliminary STAR data indicate a much higher chemical freeze-out temperature of around 170–180 MeV ).

## IV Radial and elliptic flow

### A Radial flow

The single particle spectra for positive pions and antiprotons for the different initialization models are shown in Figure 7, for two representative values of the impact parameter. For central ($`b=\mathrm{\hspace{0.17em}0}`$) collisions, the systematics of the $`m_\mathrm{T}`$-slopes shows that for models eWN, sWN, eBC and sBC the amount of radial flow created in the collision increases in the same order as the central initial energy density shown in Figure 3. The saturation model is an exception: it gives larger radial flow than both eWN and sWN although, at $`\tau _0=\mathrm{\hspace{0.17em}0.6}`$ fm/$`c`$, its central energy density is even below that of model eWN. The reason is that for the saturation model the hydrodynamic evolution starts earlier and the build-up of transverse flow thus begins already at time $`\tau _{\mathrm{sat}}(b)`$ (see Fig. 2).

The effect of the initialization on the single particle spectra is significant, and switching from model eWN used in  to model eBC, which according to Figure 5 provides a better fit to the centrality dependence of the charged multiplicity, considerably flattens the single particle spectra predicted in Ref.  for central collisions. This model sensitivity is weaker for semiperipheral collisions ($`b=\mathrm{\hspace{0.17em}8}`$ fm), again in agreement with Figure 3 which shows that in this range the initial energy density profiles are rather similar for all the models. If hydrodynamics still applies to even more peripheral collisions, models eBC and sBC would there predict steeper spectra than models eWN and sWN.

### B Elliptic flow

Figure 8 shows that for the four models eWN, sWN, eBC, and sBC the differential elliptic flow $`v_2(p_\mathrm{T})`$ follows a similar pattern as the radial flow: the model with the largest initial central energy density produces the strongest elliptic flow. The saturation model does not follow this systematics, which again is presumably due

to the earlier start of the hydrodynamic evolution in this model. The impact parameter dependence of this model sensitivity is now controlled by the initial spatial anisotropy $`ϵ_x`$. According to Figure 4, for most of the impact parameter range it satisfies the hierarchy $`ϵ_x(\mathrm{eWN})<ϵ_x(\mathrm{sWN})<ϵ_x(\mathrm{sat})<ϵ_x(\mathrm{eBC})<ϵ_x(\mathrm{sBC})`$, and this is clearly reflected in the slopes of $`v_2(p_\mathrm{T})`$ for $`p_\mathrm{T}<2`$ GeV/$`c`$, as shown in Figure 8. Again, the effect is significant and can reach up to $`20\%`$. At large impact parameters, $`b>10`$ fm, and for $`p_\mathrm{T}>2`$ GeV/$`c`$ this simple ordering is broken.

In view of these results, it is somewhat of a pleasant surprise to see in Figure 9 that the $`p_\mathrm{T}`$-integrated elliptic flow $`v_2`$ for all charged hadrons, when plotted as a function of the midrapidity density of charged particle multiplicity as an experimentally easily accessible centrality measure, shows very little sensitivity to the initialization model. The models eWN, sWN, eBC, and sBC yield almost identical results, and only for the saturation model $`v_2`$ is slightly larger (except for the most peripheral collisions). Comparing the two extreme models eWN and sBC, one sees from Figure 5 that, at fixed charge multiplicity, the latter favours smaller impact parameters which correspond to smaller $`ϵ_x`$ and $`v_2(p_\mathrm{T})`$ (Figs. 4 and 7), but flatter single particle spectra which give more weight to larger $`v_2`$ at higher $`p_\mathrm{T}`$ (Fig. 6). The net effect of this intricate interplay is an almost complete cancellation of these counteracting tendencies.

A similar cancellation happens for the differential anisotropy $`v_2(p_\mathrm{T})`$ for charged hadrons in minimum bias events (left panel of Figure 10):

$`v_2(p_\mathrm{T})={\displaystyle \frac{b𝑑bv_2(p_\mathrm{T};b)\frac{dN_{\mathrm{ch}}}{dyp_\mathrm{T}dp_\mathrm{T}}(b)}{b𝑑b\frac{dN_{\mathrm{ch}}}{dyp_\mathrm{T}dp_\mathrm{T}}(b)}}.`$ (18)
At fixed $`b`$, the models with smaller elliptic flow are weighted with larger charged multiplicities, but slightly steeper spectra, and again the net result is an almost miraculous cancellation of all sensitivities to the initialization model for $`p_\mathrm{T}\stackrel{<}{}\mathrm{\hspace{0.17em}1.5}`$ GeV/$`c`$. Above $`p_\mathrm{T}=\mathrm{\hspace{0.17em}1}`$ GeV/$`c`$, the saturation model gives the largest $`v_2(p_\mathrm{T})`$.

This cancellation does not carry over to heavier particles. The right panel in Figure 10 shows the differential elliptic flow $`v_2(p_\mathrm{T})`$ for identified protons and antiprotons, which is seen to exhibit more significant variations as the initialization of the hydrodynamic model is changed. Generically, the model which gives the largest radial flow at small impact parameters produces the smallest proton elliptic flow at small values of $`p_\mathrm{T}`$. This is qualitatively consistent with the general analytic discussion of radial flow effects on the elliptic flow for heavy particles presented in Ref. . Again, the saturation model gives the largest elliptic flow of all studied initializations for protons with $`p_\mathrm{T}>1`$ GeV/$`c`$.

## V Conclusions

Within a hydrodynamic model with exact longitudinal boost invariance, we have shown that the centrality dependence of the production of particles, transverse energy and transverse flow is significantly influenced by the shape of the initial energy density profile in the transverse plane. This profile is intimately related to the nature of the primary particle production process which converts beam energy into the matter forming the collision fireball. The available data for the charged multiplicity in Au+Au collisions at RHIC as a function of collision centrality are compatible with an initial energy deposition process involving a component which scales with the number of binary nucleon-nucleon collisions. They strongly disfavour the wounded nucleon parametrization used in earlier hydrodynamic simulations . In non-central collisions they are, at the present RHIC energy, also at variance with predictions from the gluon saturation model . The best agreement with the multiplicity data at $`\sqrt{s}=\mathrm{\hspace{0.17em}130}A`$ GeV  would be obtained by assuming that most of the initial thermalized energy density is proportional to the number of binary nucleon-nucleon collisions, with a smaller “soft” contribution proportional to the number of wounded nucleons. If one instead parametrizes the initial entropy density as a superposition of “hard” and “soft” components, as implied by the approaches in , one needs a larger “soft” and smaller ($`10\%`$) “hard” component . A model with similar properties is obtainable also in the saturation approach considered here by limiting the saturation region to smaller transverse distances and thereby increasing the contribution from the tail which scales with the number of binary collisions. All these models predict an average transverse energy per particle in the final state which rises slowly with increasing number of participating nucleons. A measurement of $`E_\mathrm{T}/N_{\mathrm{ch}}`$ as a function of $`N_{\mathrm{part}}`$ is sensitive to the longitudinal work done during the hydrodynamic evolution and can be used to check the consistency of the hydrodynamic approach.

The centrality dependence of the $`p_\mathrm{T}`$-averaged elliptic flow $`v_2`$ and the differential elliptic flow $`v_2(p_\mathrm{T})`$ for all charged particles from minimum bias collision events exhibit only minor sensitivity to the shape of the initial energy density profile. Thus, the conclusions extracted from earlier studies of elliptic flow using the wounded nucleon parametrization (which is now disfavoured by the new multiplicity data) remain valid. For identified heavy hadrons such as (anti)protons the differential elliptic flow $`v_2(p_\mathrm{T})`$ at low $`p_\mathrm{T}`$ is smaller for the binary collision induced initial conditions than for the wounded nucleon parametrizations. This stems from the somewhat stronger radial flow created by the “hard” initializations, which should also be visible in the single-particle spectra.

In closing we note that the initialization models studied in the present paper assume a one-to-one correspondence between impact parameter, number of participating nucleons, and initial energy density profile in the transverse plane. This neglects the possibility of strong event-by-event fluctuations in the initial density profile at fixed number of charged particles in the final state. As pointed out recently , strong event-by-event fluctuations in the initial energy density profile may result in sizeable event-by-event fluctuations of $`v_2`$ at fixed impact parameter. Our present study shows that small and smooth variations in the initial energy density profile at fixed impact parameter cause effects on the particle multiplicity and on the radial and elliptic flow which compensate each other in such a way that $`v_2(n_{\mathrm{ch}}/n_{\mathrm{max}})`$ is almost unaffected. It would be interesting to study whether this compensation also survives the much larger fluctuations studied in Ref. .

Acknowledgements: We thank P.V. Ruuskanen for several useful discussions. The work of P.K. was supported in part by the Deutsche Forschungsgemeinschaft. P.H. acknowledges financial support by the Director, Office of Science, Office of High Energy and Nuclear Physics, Division of Nuclear Physics, and by the Office of Basic Energy Sciences, Division of Nuclear Sciences, of the U.S. Department of Energy under Contract No. DE-AC03-76SF00098. K.J.E. and K.T. acknowledge financial support from the Academy of Finland.


text:
# Equilibrium and kinetics at the coil–to–globule transition of star and comb heteropolymers in infinitely dilute solutions

## I Introduction

From the polymer design perspective three properties of a polymer are of major importance for the resulting conformational properties and consequent industrial applications. These are the degree of stiffness of the chain, the topology of the polymer, and its composition in terms of monomers involved. In trying to model, for instance, the structure of a natural protein one has to construct the main chain of a given stiffness and to attach various side chains, each of which would have required local stiffness, degree of hydrophobicity and other more specific properties associated with the nature of the corresponding amino acid residue. The challenging nature of the kinetics of protein folding would manifest itself in the shear complexity and diversity of the kinetics of folding of our hypothetical model polymer, which should be rather sensitive on the details of these three types of interactions involved. In the current paper we shall discuss the results from Monte Carlo simulations in continuous space for the kinetics of folding of several sufficiently simple and yet already complex types of polymers based on a rather generic coarse–grained model introduced here. These preliminary studies of the model produce a number of novel transient, and in some cases metastable, kinetic intermediates and rather complex kinetic pathways of folding.

Conformational transitions of homopolymers with varying degree of flexibility have been of considerable interest for analytical studies  and computer simulations  in recent years. It is well known that folding of a sufficiently stiff chain results in formation of a torus  under certain conditions. Moreover, the transitions coil–to–torus and torus–to–globule are discontinuous, resulting in a rather complex kinetic picture of expansion or collapse, essentially dependent on the quench depth . The phase diagram of a stiff homopolymer obtained from the Gaussian variational method  is found to contain phases of the coil, a spherical (or slightly aspherical) globule in the region of either a low stiffness or a strong monomer attraction, and a number of toroidal phases characterised by distinct winding number in between. The collapse transition changes its behaviour from continuous to discontinuous starting from some value of the stiffness. This distinction is quite dramatic for the kinetics of collapse. While the kinetics of the coil–to–globule transition for a flexible homopolymer proceeds via formation of a necklace of locally collapsed sub–globules or ‘pearls’  rather than via a sausage mechanism, the pearls are not seen in the semi–flexible case where hairpins with abrupt U-turns and indeed sausages, often also corresponding to metastable states, would be typical transient kinetic intermediates .

There is also a great deal of interest in the properties of amphiphilic heteropolymers. Even in the flexible case these exhibit a variety of complex conformational transitions upon a temperature or pH variation . Apart from the stiffness and the sequence of monomers, one can also play with various polymer topologies . Indeed, random branching of polymers is quite common during synthesis. By a carefully controlled polymerisation procedure one can obtain regular stars, which presented a nice object for theoretical study due to their symmetry and simplicity lately . Star polymers are believed to find a number of future applications for coating, as additives and possibly in drug delivery systems related to their low viscosity and other interesting properties.

Speaking of computer simulations, star polymers, for instance, (see e.g. Refs. ) have been studied by a number of techniques, ranging from Monte Carlo, both on and off lattice, to Brownian and Molecular Dynamics . While dynamic Monte Carlo is one of perhaps most efficient approaches in terms of accessible time–scales, on a a lattice it has a number of obvious limitations. Indeed, the core monomer on a lattice has a great difficulty to move at all due to its high coordination number (equal to the number of arms). The resulting rejection of all attempted moves for some kind of motions is called a quasi–nonergodicity, a problem which is hard to deal with even for linear copolymers. If such a situation occurs, it practically means that the results of the simulation may no longer be trusted. A number of complicated non–local moves involving the core monomer have been suggested to alleviate this problem . However, such moves are not permitted in case of Dynamic Monte Carlo needed for study of, for instance, metastable states, moreover they may present problems due to subtle topological effects for truly non–phantom chains .

Referring to our recent paper  for examples of applications to diblock copolymer star, we devote the present work to investigating the kinetics of the coil–to–globule transition for flexible and semi–flexible homopolymer stars and diblock copolymer combs with the hydrophobic backbone (or H-backbone comb briefly) with hydrophilic (polar) side–arms and vice versa (or P-backbone comb).

## II Results

The mathematical model for a generic coarse–grained polymer is based upon the following Hamiltonian (energy functional):

$`H`$ $`=`$ $`{\displaystyle \frac{k_BT}{2\mathrm{}^2}}{\displaystyle \underset{ij}{}}\kappa _{ij}(𝐗_i𝐗_j)^2+{\displaystyle \frac{k_BT}{2\mathrm{}^2}}{\displaystyle \underset{ijk}{}}\lambda _{ijk}(𝐗_i+𝐗_k2𝐗_j)^2`$ (2)
$`+{\displaystyle \frac{1}{2}}{\displaystyle \underset{ij,ij}{}}V_{ij}(|𝐗_i𝐗_j|).`$
Here the first term presents the connectivity structure of the polymer with harmonic springs of a given strength $`\kappa _{ij}`$ introduced between any two connected monomers (which is denoted as $`ij`$). The second term presents the bending energy penalty given by the square of the local curvature with a characteristic stiffness $`\lambda _{ijk}`$ between any three consecutively connected monomers (which is denoted as $`ijk`$) in the form of the Kratky–Harris–Hearst term . Finally, the third term presents pair–wise non–bonded interactions between monomers such as van der Waals forces and so on. In a simple model we can adopt the Lennard–Jones form for the shape of the potential,

$$V_{ij}(r)=\{\begin{array}{cc}+\mathrm{},\hfill & r<d\hfill \\ V_{ij}^0\left(\left(\frac{d}{r}\right)^{12}\left(\frac{d}{r}\right)^6\right),\hfill & r>d\hfill \end{array},$$
(3)
where there is also a hard core part with the monomer diameter $`d`$ (below this is taken to be equal to the statistical length $`\mathrm{}`$, i.e $`d=\mathrm{}`$). Finally, we shall restrict ourselves to a binary heteropolymer with hydrophilic and hydrophobic units, so that the potential coefficients can be further parametrised as ,

$$V_{ij}^0=\frac{V_i+V_j}{2},V_i=V_0(1\sigma _i),\sigma _i=\pm 1.$$
(4)
For details of the adopted Monte Carlo technique we refer to our paper Ref. .

Our first subject is the collapse kinetics of a flexible ($`\lambda =0`$) homopolymer star after a quench from the purely repulsive good solvent condition ($`V_0=0`$) to the poor solvent ($`V_0=5`$). Here and below $`k_BT`$ units are used for the energy and $`\mathrm{}`$ units for the distance. From Fig. 1a one can see that the process involves formation of locally collapsed clusters along the chain similarly to the case of an open or ring chain. However, due to the connectivity structure of the star, the largest cluster is formed around the core monomer and the arms ends seem to belong to fairly large clusters. Also, it is interesting to note that some of the clusters involve monomers from more than one arm (see e.g. on the bottom right of the figure). Further, the ends of arms clusters are being pulled towards the growing core cluster, so that the total number of clusters and free arms rapidly decrease (Fig. 1b). Finally, an almost spherical globule will form and continue to optimise its shape and internal structure until reaching the final equilibrium exponentially slowly .

In Figs. 2, similarly, we examine the kinetic process for a semi–flexible homopolymer star with $`\lambda =5`$. Clearly, as in the case of a semi–flexible rings in Ref.  formation of the localised clusters is no longer possible due to a large energetic penalty. Instead, as in Fig. 2a some of the arms come together forming a kind of fibrils, some of which still have free outstretched arms. This tree–like structure further experiences the process of fibrils growth with their width increasing and the number and lengths decreasing (Figs. 2b and 2c). The overall compaction process is much slower compared to the flexible case and finally a rather long–lived sausage–like object would remain. It is interesting to point out that quite often these kinetically arrested states would have one or a few pieces of single arms sticking out of the the main irregularly shaped body as in Figs. 2c and 2d. These curious looking conformations are rather typical for a semi–flexible chain and should in principle be seen by Electron Microscopy in experiments involving long semi–stiff star polymers. This effect is clearly related to the topology of a star since similar conformations are not present for a ring polymer .

Now let us turn our attention to flexible comb copolymers. We analyse the simple case in which the backbone is constructed of one type of monomers, while the side–arms of another, one of which is hydrophobic (H) and another hydrophilic (P for polar), so that the total numbers of monomers of both kinds are equal. In particular we shall analyse the situation in which in the final state $`V_{PP}^0=0`$ (pure repulsion), $`V_{HH}^0=5`$ and $`V_{HP}^0=2.5`$, corresponding to a near ideal regime for the H-P monomer interaction. In Figs. 3 we present the characteristic equilibrium conformations for the cases of the H-backbone (Fig. 3a) and the P-backbone (Fig. 3b). The former case pretty much reminds the situation for a OUTER-P star , so that in Fig. 3a we see the collapsed hydrophobic backbone in the centre with three outstretched hydrophilic side–groups (the forth one is pointing in the perpendicular to the plane direction and thus is not visible). This object can be viewed as a micellar star. The reversed case of the P-backbone again produces a micellar object with a single collapsed core of hydrophobic side–groups and extended hydrophilic backbone chain looping around the core. This is a much more compact object compared to the H-backbone case and it is in many ways similar to the collapsed OUTER-H star in Ref. . Also, the H-core in Fig. 3b is clearly much less exposed to unfavourable contacts with the solvent molecules than in Fig. 3a.

Next, let us also discuss briefly the transient kinetic states after a quench from the good solvent condition (all $`V_{ij}^0=0`$) to the above discussed final equilibrium values in both cases. In Fig. 4a one can see the formation of locally collapsed clusters along the hydrophobic back backbone, which further grow and merge with each other until producing a single globular core with the hydrophilic side–arms remaining extended and outstretched throughout. This mechanism is akin to the case of an open homopolymer collapse except that perhaps the mobility of the growing clusters is somewhat reduced by the presence of the passive side–arms. Finally, in Fig. 4b we present a typical non–equilibrium conformation for the P-backbone comb. Here, the folding process starts very quickly by a full collapse of all individual hydrophobic side–arms, which remain farther apart from each other along the hydrophilic backbone. This arrangement results in a rather long–lived metastable state, which has a rather distinct structure from the final equilibrium state depicted in Fig. 3b. An example of a similar metastable state in kinetics was also seen for an OUTER-H star in Ref. , where we have also seen that a much deeper quench would render such a state unstable turning the nucleation regime into that of spinodal decomposition.

## III Conclusions

In this work we have presented some preliminary results on the studies of the equilibrium and kinetics of conformational changes in flexible and semi–flexible homo- and block copolymer stars and combs. We have seen a number of unusual conformational states and distinct folding pathways strongly dependent on the macromolecular architecture.

Study of the homopolymer collapse kinetics in stars of varying degree of flexibility has elucidated that the persistent length of the chain does play a major role in determining whether the necklace of clusters or the fibril growth mechanism of folding is realised. This was also true for simpler topological objects such as an open chain or a ring . However, there are some topology induced and previously unseen features such as for instance a ‘monkey and its tail’ conformation in Fig. 2d. Further increase of the stiffness may lead to some novel equilibrium phases and a complex resulting phase diagram in analogy with multiple toroidal states having distinct winding numbers for a ring homopolymer .

Also, it is interesting that a number of similarities with the previously studied case of diblock star copolymers  has been seen here for comb polymers. It is quite likely that many of these features would be also retained for more complex regularly branched dendrimers , which are currently being scrutinised.

It would be important to investigate attendant kinetic laws and estimate characteristic time scales involved as we have done for the homopolymer kinetics . In doing so, it would be tempting to use several alternative simulation techniques such as a coarse–grained Newton and Langevin Dynamics and an analytical technique based on the Gaussian self–consistent treatment  for being able to obtain the most complete picture in this challenging and important problem.

Finally, we believe that modelling chains with distinct types of side–groups is necessary for more realistic description of the kinetics of protein folding. Undoubtedly, a simple heteropolymer model of proteins does miss many of the important interactions and details present in real polypeptides, and a proper account for side–groups and charges of the residues is equally important as the hydrophobic force in driving the folding.

## ACKNOWLEDGMENTS

The authors acknowledge interesting discussions with Professor F. Ganazzoli, Professor H. Orland, Professor T. Garel and Ronan Connolly. This work was supported by grant FR/2000/019 from Enterprise Ireland and UCD President’s research award. We also acknowledge the support of the Centre for High Performance Computing Applications, UCD.