text:
# Effective String Rope Model for the initial stages of Ultra-Relativistic Heavy Ion Collisions

## 1 Introduction

Fluid dynamical models are widely used to describe heavy ion collisions. Their advantage is that one can vary flexibly the Equation of State (EoS) of the matter and test its consequences on the reaction dynamics and the outcome. This makes fluid dynamical models a very powerful tool to study possible phase transitions in heavy ion collisions - such as the liquid-gas or the Quark-Gluon Plasma (QGP) phase transition. For example, the only models that can handle the supercooled QGP are hydrodynamical models with a corresponding EoS. For highest energies achived nowadays at RHIC hydrodynamic calculations give a good description of the observed radial and elliptic flows , in contrast to microscopic models, like HIJING  and UrQMD .

In energetic collisions of large heavy ions, especially if a QGP is formed in the collision, one-fluid dynamics is a valid and good description for the intermediate stages of the reaction. Here, interactions are strong and frequent, so that other models, (e.g. transport models, string models, etc., that assume binary collisions, with free propagation of constituents between collisions) have limited validity. On the other hand, the initial and final, Freeze Out (FO), stages of the reaction are outside the domain of applicability of the fluid dynamical model.

Thus, the realistic, and detailed description of an energetic heavy ion reaction requires a Multi Module Model, where the different stages of the reaction are each described with a suitable theoretical approach. It is important that these Modules are coupled to each other correctly: on the interface, which is a three dimensional hypersurface in space-time with normal $`d\sigma ^\mu `$, all conservation laws should be satisfied, and entropy should not decrease. These matching conditions were worked out and studied for the matching at FO hypersurfaces in details in Refs. . In Ref  it has been shown that the entropy constraint is very sensitive to the pre-FO EoS, i.e. the arbitrary choice of the pre-FO state does not necessarily lead to the physical FO.

Similar ideas form the base of Ref.  where authors combined two modules - hydro and UrQMD - replacing the hadronic phase of hydrodynamics with hadronic transport model to describe properly chemical and thermal FO.

The initial stages are the most problematic. We will discuss here different models for the initial state in section 2 and will present our modeling of it in section 3.

Our goal is to build a Multi Module Model for ultra-relativistic heavy ion collisions that is valid for RHIC and LHC energies, and maybe for the most energetic SPS collisions. The present work is just the first step, but an important one: the Effective String Rope model is developed for the most problematic module - module describing the initial stages of collisions.

## 2 Initial state of Relativistic Heavy Ion Collision

The situation at the beginning of an energetic heavy ion collision is highly complicated and not easily understandable. None of the theoretical models currently on the physics market can unambiguously describe the initial stages.

As was pointed out in the introduction, at ultra-relativistic energies the one-fluid dynamical models can not be justified for such a situation. Frequently, two- or three-fluid models are used to remedy the difficulties and to model the process of QGP formation and thermalization . Here the problem is transferred to the determination of drag-, friction- and transfer- terms among the fluid components, and a new problem is introduced with the (unjustified) use of an EoS in each component in a nonequilibrated situation where an EoS is not defined. Strictly speaking this approach can only be justified for mixtures of noninteracting ideal gas components.

The kinetic transport models (see for example, Ultra-relativistic Quantum Molecular Dynamics (UrQMD) ) and parton cascades (see for example, the Molnar Parton Cascade (MPC) ) are also frequently applied to describe the high energy collisions from the very beginning. Some of the parton cascades (for example, MPC) are able to handle only $`22`$ reactions so far, what cannot be justified for dense and strongly interacting systems. The parton cascades also have to approximate media-effects, since the parton structure function, measured for a free nucleon, is modified inside the nucleus. This is the so-called shadowing effect . This approximation procedure is also not unambiguous, see for example Ref.  for more details.

Some models, like UrQMD, were initially developed for much lower energies than available nowadays, but even then, when the particle propagates in a hot and dense and nonequilibrated medium of highly excited hadrons, the properties of this particle might change significantly. Thus, properties like effective masses, decay widths, effective momenta, and in-medium cross sections have to be calculated for a given local situation at the point where particle propagates. Unfortunately this is a very complicated and numerically very time consuming task, so in most of the models drastic approximations are made (see for example ). In a relativistic heavy ion collision, especially at the early stages, the excitation energy per particle would be so high that the concept of resonance (we have to deal not only with nucleons, but with many different baryons and mesons) has no sense any more. Here the quark and gluon degrees of freedom must be taken into account. Therefore, some models, including UrQMD, replace the resonances by continuous string excitations for internal excitation energies higher than approximately $`2GeV`$. The properties of such effective “strings” are basically unknown and can be varied within a wide range, almost like free parameters .

Perturbative Quantum Chromo Dynamics (pQCD) would be, in principle, the proper model for describing our systems at very high energies. Unfortunately, pQCD itself is not applicable for heavy ion reactions at RHIC energies. Nevertheless, the pQCD calculations with some extra non-perturbative, phenomenological assumptions, like saturation of a gluon plasma, can be performed. Different models following this scenario have been proposed in . The extension of the pQCD+saturation model for nonzero impact parameters, which is of primary interest for us, leads to serious complications in calculations .

Thus, since we can not unambiguously describe the initial stages of relativistic heavy ion collisions in microscopical models, we have to use some phenomenological models (in fact all the microscopical models discussed above already have some phenomenological elements).

### 2.1 Different models

There are two phenomenological models, which are most frequently discussed in the literature, and which claim to give us at least a qualitative understanding of the initial stages of the relativistic heavy ion collisions.

### 2.2 The Landau model

The Landau model has been developed to generalize the fluid dynamical approach to the energy range where we can not neglect relativistic effects, such as Lorentz contraction of the colliding nuclei. Actually the model was originally developed for $`pp`$ collisions, but the generalization for $`AA`$ is rather straightforward and the fluid dynamical expansion is even more justified for bigger systems.

It is the most appropriate model for relativistic heavy ion collisions in the energy range $`E_{lab}=10100AGeV`$. For energies above $`\gamma ^{CM}310`$, so in the Center of Mass (CM) frame the colliding nuclei are considerably contracted and have disk-like shapes. The Landau fluid dynamical model assumes the initial thermalized state as a static, completely stopped, homogeneous disk, contracted in the beam direction by the factor $`\gamma _{CM}`$ (see Fig. 1). The fluid dynamical expansion initially starts in the direction of the largest pressure gradient - in the beam direction in our case. Only when the system has expanded in the $`z`$ direction to a size comparable to its transverse diameter the transverse expansion also considered. In this second phase the expansion is 2+1 dimensional, while initially it was 1+1 dimensional (all quantities were $`z`$ and $`t`$ dependent).

The inputs to this model are: A) the equation of state; B) the initial energy and baryon densities; C) the FO conditions. Usually the initial energy density, $`e_0`$, is estimated in the following way:

$$e_0=E_{CM}/V_0,\mathrm{where}V_0=V_{rest}/\gamma _{CM}.$$
(1)
### 2.3 The Bjorken model

If we increase the collision energy further to $`E_{CM}=100AGeV`$ or more, the expectation is that the concept of complete stopping does not work anymore. Furthermore, when two nucleons from the target and the projectile nuclei collide at these energies, their momenta are so high that they see each other as compound objects, containing three valence quarks and other partons. The point is that the proper participants at these energies are the partons and not the nucleons. The initial momenta of partons are very high and we expect these two parton clouds to cross each other and keep going with almost their original rapidities. Thus, contrary to complete stopping assumed in the Landau model, we have now the opposite scenario \- nuclei become essentially transparent to each other.

During their interpenetration, however, the target and projectile partons may exchange color charges, what leads to the creation of a chromoelectric field in the middle of the reaction volume. The classical analogy, that one can keep in mind for a better understanding, is the electric field between two condenser plates, moving away from each other. From QCD calculations we know that such a chromoelectric field should be confined in the transverse direction and have the geometry of a string. The energy density of such a string is substantial - already for the strings produced in hadronic collisions the string tension is about $`1GeV/fm`$ , and the theoretical simulations show that even more energetic objects may be produced in heavy ion collisions. Later these strings start to decay by producing quarks, antiquarks and gluons, but because of baryon number conservation, the matter in the middle of reaction volume is baryon free. Theoreticians like very much this last feature of the model, since the baryon free matter is rather pleasant for calculations.

The important input to the model is given by experimental data from ultra-relativistic $`pp`$ collisions. The rapidity distribution of the charged particle multiplicity, $`dN_{ch}/dy`$, can be divided, with reasonable accuracy, into three regions. Two smooth peaks around $`y=\pm y_0`$, where $`y_0`$ is the original rapidity, are believed to be originated by target and projectile partons, interpenetrating through each other almost without rapidity lost. Between those we can see the plateau, - $`dN_{ch}/dy=const`$, - generated by our baryon free flux tube. Assuming the same symmetry, $`dQ/dy=const`$, for all other quantities a simple hydrodynamical model with an analytical solution has been developed by Bjorken in Ref. . The above assumption, $`dQ/dy=const`$, means that our model is boost invariant - all the quantities are functions of the proper time only, $`Q=Q(\tau )`$, where $`\tau =\sqrt{t^2z^2}`$. Obviously the model has a 1+1 dimensional structure.

In order to study transverse expansions, the Bjorken model has to be supplemented with a description of this process separately, see for example . The boost invariance of the model is also very much appreciated by theoreticians, since it is then sufficient to perform simulations for the $`y=0(z=0)`$ plane only, and then boost the result to reproduce all of the phase space. The Bjorken model does not tell us anything about the energy density distribution in the plane perpendicular to the beam direction, say the $`z=0(y=0)`$ plane; therefore additional assumptions are required. If we define a grid in the $`[x,y]`$ plane, then the energy density profiles proportional to the number of binary collisions, $`e(x_i,y_j,z=0)N_{target}(x_i,y_j)N_{target}(x_i,y_j)`$, or proportional to number of participants, $`e(x_i,y_j,z=0)\left(N_{target}(x_i,y_j)+N_{target}(x_i,y_j)\right)`$, are the most popular assumptions in literature (see  for recent overview).

The initial energy density in the Bjorken model can be estimated in the following way. At the $`y=0(z=0)`$ plane the full energy is the transverse energy, $`E_{full}=E_T`$. One may also notice that for the Bjorken boost invariant model, the pseudo-rapidity ($`\eta =\mathrm{ln}\mathrm{tan}\frac{\theta }{2}`$), where $`\theta `$ is a polar angle of emitted particle, is the same as rapidity, $`\eta =y`$. If in the experiment we measure the energy $`E_T`$ emitted in a very small rapidity window, $`[0,y]`$, then we can estimate the initial energy density:

$$e_0=\frac{E_{full}}{V}=\left(\frac{E_T}{Az}\right)_{y=0}.$$

Now using the following expressions: transverse area $`A=\pi R_A^2`$, where $`R_A`$ is the radius of nucleus with A nucleons, and size in the beam direction $`z=\tau \eta =\tau _0y`$, where $`\tau _0`$ is the thermalization time, - we obtain

$$e_0=\frac{1}{\pi R_A^2\tau _0}\left(\frac{dE_T}{dy}\right)_{y=0}.$$
(2)
Here the thermalization time is unknown and should be estimated from some other model for the particular reaction (for example, $`\tau _01fm`$ is usually taken for SPS Pb+Pb $`158GeV/nucl`$ collisions). Another question is - can we really use the post FO final state measurement, $`\left(\frac{dE_T}{dy}\right)_{y=0}`$, directly to estimate the initial energy density. Briefly speaking, this means to neglect the longitudinal expansion caused by pressure during the evolution of our fireball, which would reduce our initial energy density to the one estimated from the data by eq. (2). Therefore, one may expect that the initial energy density is larger than the one given by this simple estimate.

### 2.4 Bjorken model with baryon recoil

Certainly, the assumption of non-stopping partons is an idealization: from basic physical principles we may expect that if two moving quarks stretch a string between themselves with a substantial string tension, they should be decelerated, since during such a process they are converting their kinetic energy into the energy of the string. Thus, the creation of the chromoelectric field should produce some baryon recoil for both colliding nuclei. In the first approximation this was done in Ref. .

The authors developed the model for $`pA`$ collisions. The initial momentum of the projectile, $`p`$, was considered to be so huge, that the recoil effect could be disregarded for it, while for the target partons the recoil arises from the acceleration of partons in an effective external field $`F^{\mu \nu }`$, produced in the collision. This baryon recoil was included in a way that guarantees the conservation of baryon current. At the same time the energy-momentum conservation was obviously violated - the effective field produced in the collision was considered to be external, and the projectile matter was not taken into account in the energy-momentum tensor of the system. Thus, the authors ended up with the following system of the modified fluid dynamical equations:

$$_\mu T^{\mu \nu }=F^{\nu \mu }N_\mu +\mathrm{\Sigma }_\pi ^\nu ,$$
(3)
$$_\mu N^\mu =0,$$
(4)
where $`T^{\mu \nu }=e_t\left(\left(1+c_0^2\right)u_t^\mu u_t^\nu c_0^2g^{\mu \nu }\right)`$ is an energy-momentum tensor (only target matter is taken into account), $`c_0^2`$ comes from the EoS $`P=c_0^2e`$, $`\mathrm{\Sigma }_\pi ^\nu `$ is the pion source term , and $`N^\mu =nu^\mu `$. The effective field, $`F^{\mu \nu }`$, is parametrized in the following way:

$$F^{\mu \nu }=\left(\begin{array}{cccc}0& 0& 0& \sigma \\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ \sigma & 0& 0& 0\end{array}\right),$$
(5)
where $`\sigma =const`$ is the field strength, or string tension, of the above described effective external field. The authors have found $`\sigma 34GeV/fm`$.

Without the source term the system (3, 4) can be solved analytically . The source term is necessary at the later stages of the reaction, since it allows the neutralization process, i.e. the mechanism by which the energy stored in the field is converted back into the energy of matter (see  for details). The interplay between the baryon recoil and the energy, coming from the neutralization of the strings - string decay, leads to much higher energy density in the fragmentation region than the original Bjorken model predicts. We will come up with the modification of this nice simple model in section 3.

Experiments have now entered the region where the Bjorken model is expected to be applicable, but there is no clear and unambiguous confirmation that Nature follows this scenario. The preliminary experimental results from RHIC do not show transparency - most particle multiplicities do not show a dip in the rapidity spectra, but rather a plateau around mid-rapidity, which is a sign of strong stopping. Furthermore, a very strong elliptic flow ($`v_2`$ flow component111Nowadays the expansion in Fourier series is usually applied to study the azimuthal distribution of particles :
$$E\frac{d^3N}{d^3p}=\frac{1}{\pi }\frac{d^2N}{dp_t^2dy}\left[1+\underset{n=1}{\overset{\mathrm{}}{}}2v_n\mathrm{cos}(n\varphi )\right],$$
where $`\varphi `$ is the azimuthal angle. The first term in square brackets represents the isotropic radial flow, while the others are referred to anisotropic flow. The first Fourier coefficient $`v_1=\mathrm{cos}\varphi p_x/p_t`$ is called directed flow, and the second one $`v_2=\mathrm{cos}(2\varphi )(p_x/p_t)^2(p_y/p_t)^2`$ is called elliptic flow. ) has been measured, which shows a clear peak around mid rapidity . To build such a strong elliptic flow, strong stopping and momentum equilibration are required. Also the $`\overline{p}/p`$ ratio at mid-rapidity measured at RHIC  (preliminary) is still far from one, which tells us that the middle region is not baryon-free.

### 2.5 Third flow component

As discussed in the introduction, the fluid dynamics is governed by the Equation of State of the matter, and the analysis of the resulting flow patterns turned out to be one of the best tools to extract the EoS from the outcome of a heavy ion collision. The final event shape must carry the information about the pressure development during the collision including the early stages of the collision.

The phase transition to the QGP is connected to a decrease of pressure according to most theoretical estimates, not only in strong first order phase transition models, but even if we have a smooth but rapid, gradual transition. This reduced pressure around the phase transition threshold are known for long (see e.g. ), and it was emphasized as a possible QGP signal: this “soft point” of the EoS might be possible to observe in excitation functions of collective flow data . It was pointed out earlier  that the decrease in the out-of-plane (squeeze out) emission is even more sensitive to the pressure drop and it decreases due to plasma formation stronger than the in-plane collective flow.

In this section we want to discuss another consequence of the same softening in the EoS, which is a recently (1999) identified new distinct flow pattern - the third flow component . This flow pattern can be seen in almost all theoretical fluid dynamics calculations with QGP formation, and also in experimental data, but it was not discussed earlier, see for example .

Two basic flow patterns have been predicted and detected before at lower energies:
A) the directed transverse flow in the reaction plane, or side-splash, or bounce off, which is most frequently presented on the well known $`p_x`$ vs. $`y`$ diagram (or $`v_1`$ flow component vs. $`y`$) and seen at all energies in heavy ion collisions from energies of $`30AMeV`$ to $`165AGeV`$ ;
B) the squeeze-out effect which is an enhanced emission of particles transverse to the reaction plane at center-of-mass rapidities.

At lower energies the directed transverse flow resulted in a smooth, linear $`p_x`$ vs. $`y`$ dependence at CM rapidities. This straight line behavior connecting the maximum at $`y_{proj}`$ and the minimum at $`y_{targ}`$ was so typical that it was used to compare flow data at different beam energies and masses.

At higher energies, say above $`1011AGeV`$ in Lab frame, the deviations from this straight line behavior have been observed . Unfortunately, directed flow has not yet been measured at RHIC.

In theoretical calculations the important feature of this third flow component is that it is clearly seen in fluid dynamical calculations when, and only when, the QGP formation was allowed during the calculation (see  for overview). In sharp contrast, the solutions with a purely hadronic EoS did not show this effect, and the maximum and minimum of the $`p_x`$ curve could be connected with a straight line222Note that all these fluid dynamical calculations were done much before the experiments. The first quantitative flow predictions  preceded the experiments by six years and gave rather good agreement with the data..

### 2.6 Possible source of the third flow component in fluid dynamics

The following possible source of the third flow component has been proposed in Ref. . If QGP is formed in the collision, strong and rapid stopping may take place, and a behavior close to that of a one-fluid model one might be established. I.e., let’s assume that Bjorken scenario is not realized, at least in its form discussed in previous section, but that the Landau fluid dynamical model is applicable for collisions of massive ultra-relativistic heavy ions (see the discussion at the and of section 2.4). The soft and compressible QGP forms a rather flat disk which is at rest in the CM system (see Fig. 3). Then this disk starts to expand in the direction of the largest pressure gradient. At small but finite impact parameters we may assume that this disk is tilted, and the direction of fastest expansion will stay in the reaction plane, but will deviate from the both beam axis and the usual transverse flow direction. Thus, a third flow component might develop from the large pressure gradient and tilted and strongly Lorentz contracted initial state. Now the flow in the reaction plane will have an “elliptic structure”  and must be characterized by two flow components: the usual transverse flow and the one orthogonal to it - the third flow component or antiflow. At the same time as the primary directed flow is weakened by the stronger Lorentz contraction at higher energies, this third flow component is strengthened by increased pressure gradient arising from the Lorentz contraction.

The concept of a tilted “fireball” was proposed initially in the so-called “firestreak” models . The idea behind it is rather simple and has to do with momentum conservation. Without going into mathematical details we can understand it in the following way. First, let us create a grid in the $`[x,y]`$-plane and subdivide the projectile and target into streaks parallel to the beam direction. We try to describe the collision with finite impact parameter. In the CM frame in one streak the pieces of matter from target and projectile will have the same absolute values of rapidities, but different momenta, since they in general contain different amount of matter (see Fig. 4). If the collision happens so fast that there is no interaction between streaks, then the momentum conservation does allow complete stopping, i.e., the final streak - the piece of matter within one tube with the same transverse coordinates after collision of target and projectile fragments \- will have a momentum defined by momentum conservation, and correspondingly its center of mass has to move with some rapidity. Crossing the participant region along $`x`$ axes in the reaction plane we will see that in the middle the final streaks rapidities are small, since the target and projectile fragments were almost equal; while at the edges of the participant region the final streaks move rather fast (see Fig. 4). Thus, after some time our tilted “fireball” will degenerate into a “firestreak”. Please notice the difference between “firestreak” and the tilted Landau disk discussed above: the tilted disk was at rest in CM frame, i.e., the complete stopping has been assumed, while “firestreak” geometry is dynamic. By “dynamic” we mean that different parts of “firestreak” move with different rapidities.

In “firestreak” models such a geometry has been assumed at the freeze-out time, and then, based on energy, charge and baryon number conservations (we already conserve momentum) the authors tried to reproduce the parameters of the post FO distribution . This procedure is somewhat similar to the statistical production model, which is used nowadays to calculate production of different hadronic species at chemical FO hypersurface. Such calculations were not really supported by experiment. Nevertheless, in our opinion, this scenario might be a good qualitative description for the initial stages of the heavy ion collisions until the local thermalization achieved. When the matter is locally thermalized and the pressure is built up, we cannot neglect hydrodynamical expansion, which will smear out this initial distribution producing a more or less spherical fireball before the freeze out process starts. Nevertheless, such a hydrodynamical expansion will initially start in the direction of the largest pressure gradient (see Fig. 4) and, together with the initial velocity field, might produce a two component (elliptic) structure of the flow in the reaction plane, identified in ref. .

## 3 Effective String Rope model

In this section we will discuss so called Effective String Rope model (ESRM) for energy, pressure and flow velocity distributions at the beginning of ultra-relativistic heavy ion collisions. The output of this model can be used as an initial condition for further hydrodynamic calculations.

One important conclusion of heavy ion research in the last decade is that standard ’hadronic’ string models fail to describe heavy ion experiments. All string models had to introduce new, energetic objects: string ropes , quark clusters , or fused strings , in order to describe the abundant formation of massive particles like strange antibaryons. Based on this, we describe the initial moments of the reaction in the framework of classical, or coherent, Yang-Mills theory, following the Bjorken model with baryon recoil Ref.  (see section 2.4) assuming a larger field strength (string tension) than in ordinary hadron-hadron collisions. For example, calculations both in the Quark Gluon String Model (QGSM)  and in the Monte Carlo string fusion model  indicate that the energy density of strings reaches $`810GeV/fm`$ already in SPS reactions. This is nearly $`10`$ times more than the tension used in standard, ’hadronic’, string models where $`\sigma 1GeV/fm`$. In addition we now satisfy all conservation laws exactly, while in Ref. , as it was discussed in section 2.4, energy-momentum conservation was irrelevant and was ignored. Thus, in this approach for the first time the initial transparency/stopping and energy deposited into strings and “string ropes” was determined consistently with each other .

We do not describe the hadronization of our strings. Thus, this approach can be justified only for the very initial stages. The effective way to perform string decay will be discussed in section 3.3.

Recent parton kinetic models indicate that saturation of a gluon plasma takes place in a very short time - $`\tau _{sat}=0.09,0.27fm/c`$ for LHC, SPS energies respectively , and a pressure needs $`\tau _p=12fm/c`$ (SPS) and $`\tau _p=35fm/c`$ (LHC) to be established . More importantly the first experiments at RHIC yield strong elliptic flow, which, as it was mention in introduction, was not reproduced in any other model, except in fluid dynamical models with a QGP EoS . This is a strong experimental indication that transverse pressure builds up early in these reactions, in a few $`fm/c`$, and strong stopping is also necessary to create strong flow before Freeze Out, which usually happens when the radius of the system is not more than $`10fm`$. Our model builds up the locally equilibrated matter \- the initial conditions for further hydrodynamical expansion - at $`t_{lab}=3/5fm/c`$ for central/peripheral regions at RHIC energies, which is in agreement with previous estimations as well as the data.

The presentation in this section is based on Refs. , but we want to stress that the preliminary results were presented at conferences and published in their proceedings, while the model was still being developed and, thus, Ref.  is more up to date than Refs. . And in this work we also improved the presentation of our model, particularly in sections 3.3, 3.4, 3.5. For the continuous description we will repeat in this section some details from Ref. , but we will also underline the differences.

### 3.1 Formulation of model

The basic idea is to generalize the Bjorken model with baryon recoil, presented in section 2.4, for collisions of two heavy ions and improve it by strictly satisfying conservation laws . Similar ideas were also presented recently in .

First of all, we would like to create a grid in the $`[x,y]`$-plane ($`z`$ is the beam axes, $`[z,x]`$ is the reaction plane). We will describe the nucleus-nucleus collision as a sum of independent streak-by-streak collisions corresponding to the same transverse coordinates, $`\{x_i,y_j\}`$. We assume that baryon recoil for both target and projectile arise from the acceleration of partons in an effective field, $`F^{\mu \nu }`$, produced in the interaction. Of course, the physical picture behind this model is based on chromoelectric flux tube or string models, but for our purpose we consider $`F^{\mu \nu }`$ as an effective Abelian field. The most important consequence of the non-Abelian fields, i.e., its self interaction and the resulting flux tubes of constant cross section, are, nevertheless, reflected in our model: we assume that the field is one-dimensional. The fields generated by the colliding streaks are of constant cross section during the whole evolution, and only their lengths increase with time. As the string tension is constant, the energy of the string increases linearly with increasing length. The single phenomenological parameter we use to describe this field must be fixed from comparison with experimental data.

We describe the streak-streak collision based on the following set of equations:

$$_\mu \underset{i}{}T_i^{\mu \nu }=\underset{i}{}F_i^{\nu \mu }q_iN_{i\mu },$$
(6)
$$_\mu \underset{i}{}N_i^\mu =0,i=1,2,$$
(7)
where $`N_i^\mu `$ is the baryon current of $`i`$th nucleus, $`q_i`$ is the color charge (discussed later in more detail). We are working in the Center of Rapidity Frame (CRF), which is the same for all streaks. We will use the parameterization:

$$N_i^\mu =n_iu_i^\mu ,u_i^\mu =(\mathrm{cosh}y_i,\mathrm{sinh}y_i).$$
(8)
$`T^{\mu \nu }`$ is the energy-momentum tensor. It consists of five parts, corresponding to two nuclei and free field energy (also divided into two parts) and one term defining the QCD perturbative vacuum. (While in the similar equation in the Bjorken model with baryon recoil, eq. (3), $`T^{\mu \nu }`$ contains only the target matter.)

$$T^{\mu \nu }=\underset{i=1,2}{}T_i^{\mu \nu }+T_{pert}^{\mu \nu }=$$

$$\underset{i=1,2}{}\left[e_i\left(\left(1+c_0^2\right)u_i^\mu u_i^\nu c_0^2g^{\mu \nu }\right)+T_{F,i}^{\mu \nu }\right]+Bg^{\mu \nu }.$$
(9)
Here $`B`$ is the bag constant, the equation of state is $`P_i=c_0^2e_i`$, where $`e_i`$ and $`P_i`$ are the energy density and pressure of QGP.

Comparing our equation (6) with the corresponding one in the Bjorken model with baryon recoil, eq. (3), one may notice that we neglected the pion source term. As it was discussed in section 2.4 the source term takes care about the string hadronization process. We are planning to use our model at the beginning of the ultra-relativistic heavy ion collisions only, and the string decay process at later stages will be described separately. Section 3.3 is devoted to this question. The advantage of equations (6, 7) is a simple analytic solution.

Within each streak we form only one flux tube with uniform field strength or field tension, $`\sigma `$, between the target and the projectile streaks. The field extends to the outside edges of the target and projectile streaks, i.e., the field also overlaps with the space-time domains of the matter streaks. For practical and symmetry purposes we, however, divide this field into two spatial domains, a target and a projectile domain, ($`i=1,2`$), separated at a fixed point, $`z_{sep}`$, so that $`\sigma _1=\sigma _2=\sigma `$. The choice of this point will be specified later. (The field is constant and the only change is that it extends with time at its two ends.)

In complete analogy to electro-magnetic field

$$F_i^{\mu \nu }=^\mu A_i^\nu ^\nu A_i^\mu =\left(\begin{array}{cc}0& \sigma _i\\ \sigma _i& 0\end{array}\right),$$
(10)
$$\sigma _i=^3A_i^0^0A_i^3,$$
(11)
$$T_{F,i\mu \nu }=g_{\mu \nu }_{F,i}+\underset{\beta }{}\frac{_{F,i}}{\left(^\mu A_i^\beta \right)}_\nu A_i^\beta ,$$
(12)
$$_{F,i}=\frac{1}{4}F_{i\mu \nu }F_i^{\mu \nu }.$$
(13)
In our case the string tensions, $`\sigma _i`$, will be constant in the space-time region after string creation and before string decay. The creation of fields will be discussed later in more detail.

Since during the interpenetration both nuclei are closely following the light cone, the geometry of the reaction suggests to use the light cone variables

$$(z,t)(x^+,x^{}),x^\pm =t\pm z,$$
(14)
and to assume that $`e_1,y_1,n_1`$ are functions of $`x^{}`$ only (the first nucleus goes from right to left, $`y_1(0)=y_0<0`$) and $`e_2,y_2,n_2`$ depend on $`x^+`$ only (the second nucleus goes from left to right, $`y_2(0)=y_0>0`$).

In terms of lightcone variables (for the transformation matrixes between space-time and lightcone coordinates see Ref. ):

$$N_i^\pm =N_{i,}=n_i(u_i^0\pm u_i^3)=n_ie^{\pm y_i},$$
(15)
$$T_i=\left(\begin{array}{cc}T_i^{++}& T_i^+\\ T_i^+& T_i^{}\end{array}\right)=\frac{1}{2}\left(\begin{array}{cc}h_{i+}e^{2y_i}& h_i\\ h_i& h_{i+}e^{2y_i}\end{array}\right)+T_{F,i},$$
(16)
where

$$h_{i+}=(1+c_0^2)e_i,h_i=(1c_0^2)e_i.$$
(17)
The other tensors in light cone variables are:

$$F_i=\left(\begin{array}{cc}F_i^{++}& F_i^+\\ F_i^+& F_i^{}\end{array}\right)=\left(\begin{array}{cc}0& \sigma _i\\ \sigma _i& 0\end{array}\right).$$
(18)
$$T_{pert}=\left(\begin{array}{cc}0& B\\ B& 0\end{array}\right).$$
(19)
The energy-momentum tensor for free field in the light cone variables is:

$$T_{F,i}=\frac{1}{2}\left(\begin{array}{cc}\sigma _i^2& 0\\ 0& \sigma _i^2\end{array}\right).$$
(20)
At the time of the first contact of the two streaks, $`t=0`$, there is no string tension. We assume that strings are created, i.e., the string tension achieves the value $`\sigma `$ at time $`t=t_0`$, corresponding to complete penetration of streaks through each other (see Fig. 5).

### 3.2 Conservation laws and string creation

In lightcone variables eq. (7) may be rewritten as

$$_{}N_1^{}+_+N_2^+=0.$$
(21)
So, we have a sum of two terms, each depending on different independent variables, and the solution can be found in the following way:

$$\begin{array}{cc}_{}N_1^{}=a,\hfill & _+N_2^+=a,\hfill \\ N_1^{}=ax^{}+(N_1^{})_0,\hfill & N_2^+=ax^++(N_2^+)_0,\hfill \end{array}$$
(22)
where the index $`0`$ indicates the initial proper density, which is the normal nuclear density: $`n_0=0.145fm^3`$. Since both $`N_1^{}`$ and $`N_2^+`$ are positive (and also more or less symmetric) we can conclude that for our case $`a=0`$.

Finally,

$$N_1^{}=n_1e^{y_1}=n_0e^{y_0},N_2^+=n_2e^{y_2}=n_0e^{y_0},$$
(23)
$$n_1=n_0e^{y_0+y_1},n_2=n_0e^{y_0y_2}.$$
(24)
where $`y_0(y_0)`$ is the initial rapidity of nucleus $`2(1)`$ in the CRF, respectively. The other components are given by eq. (15).

Let us make an analogy to the electro-magnetic field, where two charges, $`q_1`$ and $`q_2`$, move in the opposite directions, creating string-like field between them, $`\stackrel{}{E}=(0,0,E)`$, which is constrained transversally into a constant cross section. The Z-axis goes through the charges, $`q_1`$ and $`q_2`$, and directed from $`q_1`$ to $`q_2`$ (let us assume that we have such a field configuration). So, forces acting on our charges, $`q_1E`$ and $`q_2E`$, have opposite sign and both are working against the expansion of the “string”. In our effective model we use color charges, and assume that the vectors of these color charges point in the opposite directions in the color space , so that the forces acting on both target and projectile partons are opposite, both stopping the expansion of the streak. As our field strength (string tension, $`\sigma `$) is not yet defined we normalize the charges to unity:

$$q_1=q_2=1>0,\mathrm{while}\sigma _1=\sigma _2=\sigma .$$
(25)
Then, we have the forces acting in $`z`$ direction: $`q_1\sigma _1=\sigma ,`$ and $`q_2\sigma _2=\sigma `$. Notice again that after string creation fields $`\sigma _1(x)`$ and $`\sigma _2(x)`$ are spatially separated as well as the baryon densities, $`n_1`$ and $`n_2`$; i.e., after complete penetration of the initial streaks through each other (see Fig. 5), $`\sigma _2`$ acts on the partons on the right side of the separating point $`z_{sep}=(l_1l_2)/2`$, while the $`\sigma _1`$ acts on those on the left side. In the absence of matter, in the middle both fields are identical, so the exact position of the separating point does not play any role until it does not enter the target or projectile matter. The fields $`\sigma _i`$ are generated by the corresponding four-potentials $`A_i`$, which are different and spatially separated in the same way (see the dotted separating line, $`z_{sep}=(l_1l_2)/2`$, in Fig. 5).

As described above we do not generate the chromoelectric field self-consistently as a product of color currents, which are also affected by the field. Our effective fields are external with respect to the colliding partons; that’s why we can use the expression (20) for the field energy. On the other hand, if we want to satisfy the conservation laws, we must generate our effective fields in the collision transferring energy from matter to field. It’s possible to define new conserved quantities based on eq. (6). Using the definition of $`F^{\mu \nu }`$, eq. (11), we can rewrite eq. (6) as

$$_\mu T^{\mu \nu }=\underset{i}{}F_i^{\mu \nu }q_iN_{i,\mu }=\underset{i}{}q_i[^\mu \left(A_i^\nu N_{i,\mu }\right)$$

$$A_i^\nu ^\mu N_{i,\mu }^\nu \left(A_i^\mu N_{i,\mu }\right)+A_i^\mu ^\nu N_{i,\mu }].$$
(26)
The solutions for $`N_1^{}`$ and $`N_2^+`$, eq. (23), show that the second term vanishes. The fourth term is a vector $`(A_1^{}_{}N_1^+,A_2^+_+N_2^{})`$ in lightcone coordinates. So, if we impose the conditions

$$A_1^{}=0,A_2^+=0$$
(27)
we can define a new energy-momentum tensor $`\stackrel{~}{T}^{\mu \nu }`$, such that

$$_\mu \stackrel{~}{T}^{\mu \nu }=0,$$
(28)
$$\stackrel{~}{T}^{\mu \nu }=\underset{i}{}\stackrel{~}{T_i}^{\mu \nu }+T_{pert}^{\mu \nu }=$$

$$\underset{i}{}\left(T_i^{\mu \nu }q_iA_i^\nu N_i^\mu +g^{\mu \nu }q_iA_i^\alpha N_{i\alpha }\right)+Bg^{\mu \nu }$$
(29)
One may notice that the above defined new energy-momentum tensor is not symmetric. We know that the energy-momentum tensor is not defined in unique way - if $`T^{\mu \nu }`$ satisfies the eq. $`_\mu T^{\mu \nu }=0`$, then any other tensor $`T^{\mu \nu }+_\kappa \mathrm{\Psi }^{\mu \nu \kappa }`$, satisfies it also, if $`\mathrm{\Psi }^{\mu \nu \kappa }=\mathrm{\Psi }^{\nu \mu \kappa }`$ . How to make it symmetric will be discussed later in this section.

To satisfy the above choice of the fields, (25), and imposed conditions, (27), we take the vector potentials in the following form:

$$\begin{array}{cc}A_1^{}=0,\hfill & A_1^+=\sigma _1x^+=\sigma x^+,\hfill \\ A_2^{}=\sigma _2x^{}=\sigma x^{},\hfill & A_2^+=0.\hfill \end{array}$$
(30)
In principle the vector potentials are not defined in unique way - the general expression would be:

$$\begin{array}{cc}A_1^{}=0,\hfill & A_1^+=\sigma x^++a_1x^{},\hfill \\ A_2^{}=\sigma x^{}+a_2x^+,\hfill & A_2^+=0.\hfill \end{array}$$
(31)
In our published work, Ref. , the choice $`a_1=a_2=0`$ has been used. It appears that this intuitive choice has a deep meaning - with any other $`a_1`$ and $`a_2`$ the symmetrization procedure, which will be discussed below, can not be applied.

How to define the $`\sigma `$ itself? In our calculations we used the parameterization (see Refs. ):

$$\sigma =A\left(\frac{\epsilon _0}{m}\right)^2n_0\sqrt{l_1l_2}=A\left(\frac{\epsilon _0}{m}\right)\frac{\sqrt{𝗇_\mathrm{𝟣}𝗇_\mathrm{𝟤}}}{xy},$$
(32)
where $`m`$ is the nucleon mass, $`\epsilon _0`$ is the initial energy per nucleon, $`l_1`$ and $`l_2`$ are the initial streak lengths in the CRF (see Fig. 5), $`𝗇_1`$ and $`𝗇_2`$ are the numbers of baryons in the initial streaks, $`xy`$ is the cross section of the streaks, and we’ve used the relation $`\gamma =\epsilon _0/m`$. Similar $`(𝗇_\mathrm{𝟣},𝗇_\mathrm{𝟤})`$-dependence of $`\sigma =`$ has been chosen in Ref. . We are working in the system where $`\mathrm{}=c=1`$, so $`\sigma `$ has a dimension of $`length^2=energy/length`$. The typical values of dimensionless parameter $`A`$ are around $`0.060.08`$ (these values of $`A`$ are for $`\sigma `$ measured already in $`GeV/fm`$). Notice, that there is only one free parameter in parameterization (32). The typical values of $`\sigma `$ are $`410GeV/fm`$ for $`\epsilon _0=65GeV`$ per nucleon, and $`\sigma 615GeV/fm`$ for $`\epsilon _0=100GeV`$ per nucleon. These values are consistent with the energy density of the non-hadronized strings in a reaction volume, or “latent energy density”, which is about $`15GeV/fm^3`$ for $`\sqrt{s}=200GeV/nucl`$ . Our $`\sigma `$ given by eq. (32) grows linearly with $`\epsilon _0`$. This may be a too strong assumption, since we predict much bigger field strength for LHC then the Monte Carlo string fusion model . It is easy to modify expression (32) by changing the power $`2`$ to $`1+\alpha `$, where $`\alpha <1`$, and probably this will be done after the comparison with experimental data. In Ref.  $`\alpha =0.3`$ has been chosen based on the low-x behavior of the nuclear structure function or parton density . Since such an improvement is straightforward, we will keep it in mind, but for further discussion we will use the definition (32) as in Ref. .

Using the exact definition of $`A_i^\mu `$, eqs. (30), eqs. (16, 19, 20, 25, 29) and transformation matrixes from Ref.  we obtain

$$\stackrel{~}{T}^{\mu \nu }=\left(\begin{array}{cc}\stackrel{~}{T}^{++}& \stackrel{~}{T}^+\\ \stackrel{~}{T}^+& \stackrel{~}{T}^{}\end{array}\right)=\frac{1}{2}\left(\begin{array}{cc}h_{1+}e^{2y_1}& h_1\\ h_1& h_{1+}e^{2y_1}\end{array}\right)$$

$$+\frac{1}{2}\left(\begin{array}{cc}h_{2+}e^{2y_2}& h_2\\ h_2& h_{2+}e^{2y_2}\end{array}\right)+\frac{1}{2}\left(\begin{array}{cc}\sigma ^2& 2B\\ 2B& \sigma ^2\end{array}\right)$$
(33)
$$+\left(\begin{array}{cc}\sigma x^+N_1^+& 0\\ \sigma x^+N_1^{}& 0\end{array}\right)+\left(\begin{array}{cc}0& \sigma x^{}N_2^+\\ 0& \sigma x^{}N_2^{}\end{array}\right)$$

Notice that the perturbative vacuum term, $`B`$, and the free field energy, $`\frac{\sigma ^2}{2}`$, cover all the interacting volume, while the energy densities of matter and baryon currents are separated in space.

Now we are coming back to the question: how to symmetrize $`\stackrel{~}{T}^{\mu \nu }`$? One may notice that the tensor

$$\stackrel{~}{T}_1^{\mu \nu }=\left(\begin{array}{cc}0& \sigma x^{}N_2^+\\ \sigma x^+N_1^{}& 0\end{array}\right),$$
(34)
satisfies the equation $`_\mu \stackrel{~}{T}_1^{\mu \nu }=0`$ (it’s easy to check using the eqs. (23)). Thus we can redefine $`\stackrel{~}{T}^{\mu \nu }`$ in a symmetric way:

$$\stackrel{~}{T}^{\mu \nu }:=\stackrel{~}{T}^{\mu \nu }\stackrel{~}{T}_1^{\mu \nu }.$$
(35)
This important step has not been done in Refs. . It causes some changes in the calculations of the energy densities at $`t=t_0`$ (see A), but we will see that this does not affect appreciably the final results.

Now the new conserved quantities are

$$Q^0=\stackrel{~}{T}^{00}𝑑V=xy\stackrel{~}{T}^{00}𝑑z,$$
(36)
$$Q^3=\stackrel{~}{T}^{03}𝑑V=xy\stackrel{~}{T}^{03}𝑑z.$$
(37)
We can rewrite the energy-momentum tensor in $`(t,z)`$ coordinates:

$$\stackrel{~}{T}^{\mu \nu }=\left(\begin{array}{cc}\stackrel{~}{T}^{00}& \stackrel{~}{T}^{03}\\ \stackrel{~}{T}^{30}& \stackrel{~}{T}^{33}\end{array}\right)=\left(\begin{array}{cc}\frac{\sigma ^2}{2}+B& 0\\ 0& \frac{\sigma ^2}{2}B\end{array}\right)+$$

$$\left(\begin{array}{cc}(e_1+P_1)\mathrm{cosh}^2y_1P_1& (e_1+P_1)\mathrm{cosh}y_1\mathrm{sinh}y_1\\ (e_1+P_1)\mathrm{cosh}y_1\mathrm{sinh}y_1& (e_1+P_1)\mathrm{sinh}^2y_1+P_1\end{array}\right)+$$

$$\left(\begin{array}{cc}(e_2+P_2)\mathrm{cosh}^2y_2P_2& (e_2+P_2)\mathrm{cosh}y_2\mathrm{sinh}y_2\\ (e_2+P_2)\mathrm{cosh}y_2\mathrm{sinh}y_1& (e_2+P_2)\mathrm{sinh}^2y_2+P_2\end{array}\right)$$
(38)
$$\frac{\sigma x^+}{2}\left(\begin{array}{cc}N_1^+& N_1^+\\ N_1^+& N_1^+\end{array}\right)\frac{\sigma x^{}}{2}\left(\begin{array}{cc}N_2^{}& N_2^{}\\ N_2^{}& N_2^{}\end{array}\right).$$

Obviously, in the absence of the fields, before string creation and after string decay, $`\stackrel{~}{T}^{\mu \nu }T^{\mu \nu }`$ and, thus, $`(Q^0,Q^3)`$ reverts to $`(P^0,P^3)`$ \- components of the four-momenta of the system.

Based on conservation of $`Q^0,Q^3`$ we can calculate energy densities, $`e_1(t_0),e_2(t_0)`$, at the moment $`t=t_0`$, when the string with tension $`\sigma `$ is created - eqs. (89, 90). Another result from A we will need for later discussion is the rapidity of the CM frame:

$$\mathrm{tanh}y_{CM}=Q^3/Q^0=v_0/M,$$
(39)
where $`M=\frac{l_2+l_1}{l_2l_1}.`$

Thus, for $`x^\pm >x_0`$ we should solve eqs. (28), with boundary conditions

$$\begin{array}{cc}N_1^\pm (x^{}=x_0)=n_0e^{y_0}\hfill & N_2^\pm (x^+=x_0)=n_0e^{\pm y_0}\hfill \\ h_{1+}(x^{}=x_0)=e_1(t_0)(1+c_0^2)\hfill & h_{2+}(x^+=x_0)=e_2(t_0)(1+c_0^2)\hfill \\ y_1(x^{}=x_0)=y_0\hfill & y_2(x^+=x_0)=y_0\hfill \\ \sigma _1(x^{}=x_0)=\sigma \hfill & \sigma _2(x^+=x_0)=\sigma \hfill \\ q_1(x^{}=x_0)=1\hfill & q_2(x^+=x_0)=1,\hfill \end{array}$$
(40)
where $`x_0=2t_0|z(0)|`$ defines the string creation surface $`t=t_0`$, for a parton or cell element in the position $`z=z(0)`$ at the time $`t=0`$. The complete analytical solution can be found in Ref. .

### 3.3 Hadronization of the String Ropes

If we let the partons (or cell domains) evolve according to trajectories derived in our model (see eqs. (37, 38) in Ref. ), the partons from the colliding initial streaks will keep going in the initial direction, gradually slowing down up to some time $`t=t_{i,turn}`$; then they will turn and go backwards until the two streaks again penetrate through each other (see Fig. 6). Of course, it is hard to believe that such a contraction of the reaction volume could really happen in heavy ion collisions, because of string decays, string-string interaction, interaction between streaks and other reasons, which would be very difficult to take into account. As we already pointed out earlier we did not install any string hadronization mechanism in our model and, thus, there is no way to decrease the energy supplied in the strings. If we try to describe the evolution of the system with our model after one or both initial streaks turn back, we run into the trouble, since the length of the string starts to decrease, so the field energy should also decrease, but it has no chance to do this. Therefore the above simple analytical model fails to describe the Yo-Yo motion for the initial streaks - it has the realistic character until the length of the final streak has started to decrease, i.e., approximately up to the dash-dotted line in Fig. 6.

In Refs.  we assumed that the final result of the collision of two streaks, after stopping the string’s expansion and after its decay, is one streak of the length $`l_f`$ with homogeneous energy density distribution, $`e_f`$, moving like one object with rapidity $`y_f`$. We assumed that this is due to string-string interactions and string decays, which we are not going to describe in our simple model. As mentioned above the typical values of the string tension, $`\sigma `$, are of the order of $`10GeV/fm`$, and these may be treated as several parallel strings. The string-string interaction will produce a kind of “string rope” between our two streaks, which is responsible for final energy density distribution.

One of the simplest ways to quantitatively take into account string decays is presented in Ref. . One of the complications appearing in this case is that the hadronization of the strings has to be described in the proper time (i.e. all quantities became functions of $`(\tau ,\eta ),\tau =\sqrt{t^2z^2}`$) which is, in general, different for all the pairs of colliding streaks. Then in order to be able to transform all of the final streaks to one frame, say CRF, we would have to know the details of the their dynamics for a long period of the proper time, $`t=\sqrt{\tau ^2+z^2}`$. Unfortunately this is not the case in real modeling, for example in Ref.  authors do not describe the dynamics after all the chromoelectric strings for the particular streak pair decay.

The real situation may be more complicated: when the energy accumulated in the strong color fields will be finally released in a production of $`q\overline{q}`$ pairs and gluons, this process may noticeably change the composition of matter as compared to the chemical equilibrium case . Therefore, the matter created after the mutual stopping of interpenetrating streaks can not, in general, be described by the equilibrium EoS.

For the stopped and equilibrated matter a set of homogeneous final streaks is the simplest assumption. Its advantages are: simple expressions for $`e_f,y_f`$ and a very simple way to generate an initial state for further hydrodynamical description. We have to give only four numbers (coordinates of the ends of the final streak, $`z_{left}`$ and $`z_{right}`$, $`e_f`$ and $`y_f`$) for each point of a transverse grid, $`\{x_i,y_j\}`$ . Later in section 3.5 we present an analytical solution for the expanding final streaks which seems to be more realistic. Remember - we describe the initial state which is not directly observable in experiments. Thus, even a flat initial rapidity distribution may end up after the hydrodynamical evolution in both a forward-backward peaked or a centrally peaked distribution depending on several circumstances.

The final energy density and rapidity, $`e_f`$ and $`y_f`$, should be determined from conservation laws. The conservation of the energy and momentum gives for the final rapidity:

$$\mathrm{cosh}^2y_f=\begin{array}{c}(M^2(1+c_0^2)2c_0^2v_0^2)+\\ \underset{¯}{\sqrt{(M^2(1+c_0^2)2c_0^2v_0^2)^2+4c_0^4v_0^2(M^2v_0^2)}}\\ 2(1+c_0^2)(M^2v_0^2)\end{array},$$
(41)
where we neglected $`Bl_f`$ next to $`Q^0/xy`$ and use notation $`M`$, introduced in eq. (84); $`v_0=\mathrm{tanh}y_0`$ is the initial velocity333Please notice the misprint in equation for $`\mathrm{cosh}^2y_f`$ in Refs. : $`(M^2(1+c_0^2)+2c_0^2v_0^2)`$ instead of $`(M^2(1+c_0^2)2c_0^2v_0^2)`$.. Then the expressions for $`e_f`$ is:

$$e_f=\frac{\frac{Q^0}{xy}}{((1+c_0^2)\mathrm{cosh}^2y_fc_0^2)l_f}.$$
(42)
The typical trajectory of the streak ends is presented in Fig. 7 (A). From $`t=t_0`$ they move according to trajectories derived in Ref.  (see eqs. (37, 38)) until they reach the rapidity $`y_i=y_f`$. Later, the final streak starts to move like one object with uniform rapidity, $`y_f`$, until we reach the time when the fluid dynamical calculation starts. The time and position of final streak formation can be found from the condition $`y_i=y_f`$ (see eq. (43) in Ref. ). We shall call this scenario the “nonexpanding final streak assumption”.

Unfortunately this simple model of the final state causes a principle problem. Comparing expressions (41) and (39) one may notice a strange thing: the final streak moves like one object with a certain rapidity calculated from energy-momentum conservation, but this rapidity, $`y_f`$, may differ, in general, from the rapidity of the Center of Mass of the system given by eq. (39). Obviously, if all the elements of the system move with the same rapidity, then the Center of Mass of the system also moves with the same rapidity. This problem have been pointed out in the Refs. , but the above oversimplified assumptions have been kept because of the simple analytical solution for the $`e_f`$, $`y_f`$. We think that now we have a deeper understanding of the source of the problem and we would like to update the discussion in Ref. . This is illustrated and elucidated in next section.

### 3.4 Flow velocity field generated by pressure in relativistic fluid mechanics

Let us propose to the reader the following problem:

Imagine the homogeneous piece of matter with energy density $`e`$ and Equation of State $`P=P(e)`$ moving with velocity $`v_0`$ in the Lab frame (see Fig. 8). The question is to find the energy-momentum tensor, $`T^{\mu \nu }`$, and four-momentum, $`P^\mu `$, in the Center of Mass frame and in the Lab frame and show that they are consistent with each other.

#### CM frame:

In the CM frame, which is in our case the rest frame of matter the energy momentum tensor and four-momentum of the system are:

$$T^{\mu \nu }=\left(\begin{array}{cccc}e& 0& 0& 0\\ 0& P& 0& 0\\ 0& 0& P& 0\\ 0& 0& 0& P\end{array}\right),$$
(43)
$$P^\mu =(E,0,0,0),$$
(44)
where $`E=eV_0`$, $`V_0`$ is the volume of our matter in the CM frame, $`\mu ,\nu =0,1,2,3`$.

The four-momentum is connected to the energy-momentum tensor in the following way :

$$P^\mu =T^{\mu \nu }𝑑S_\nu ,$$
(45)
where the integral is over any 3-dimensional hypersurface covering all the volume, i.e., intersecting all the world lines of the matter in space-time. If we integrate over the hyperplane $`x^0=t=const`$, then

$$P^\mu =T^{\mu 0}𝑑V.$$
(46)
It’s easy to check that equations (43) and (44) are indeed connected by the relation (46) for any EoS.

#### Lab frame:

In the Lab frame our system is moving with velocity $`\stackrel{}{v}=(0,0,v_0)`$. The four-momentum undergoes the Lorentz transformation:

$$P_L^\mu =(E\gamma _0,0,0,E\gamma _0v_0).$$
(47)
The invariant expression for the energy-momentum tensor in any reference frame is

$$T^{\mu \nu }=(e+P)u^\mu u^\nu Pg^{\mu \nu }.$$
(48)
In our case $`u^\mu =\gamma _0(1,0,0,v_0)`$. Thus,

$$T_L^{\mu \nu }=\left(\begin{array}{cccc}(e+P)\gamma _0^2P& 0& 0& (e+P)\gamma _0^2v_0\\ 0& P& 0& 0\\ 0& 0& P& 0\\ (e+P)\gamma _0^2v_0& 0& 0& (e+P)\gamma _0^2v_0^2+P\end{array}\right),$$
(49)
Now we have to use eq. (46) to find $`P_L^\mu `$ from $`T_L^{\mu \nu }`$.

$$P_L^0=𝑑V\left[(e+P)\gamma _0^2P\right]=\left[(e+P)\gamma _0^2P\right]V_0/\gamma _0,$$
(50)
$$P_L^3=𝑑V(e+P)\gamma _0^2v_0=(e+P)\gamma _0^2v_0V_0/\gamma _0,$$
(51)
where $`V_0/\gamma _0`$ is the volume of our system in Lab frame. Comparing eqs. (50, 51) with eq. (47) we obtain the following system of equations:

$$\{\begin{array}{c}\frac{(e+P)\gamma _0^2P}{\gamma _0}=e\gamma _0\hfill \\ (e+P)\gamma _0v_0=e\gamma _0v_0\hfill \end{array}$$
(52)
and we see that for a consistent description we have to put $`P=0`$. Why do we get such a strong restriction on our EoS?!

Certainly the problem is that we made wrong assumptions. We assumed that all the parts of our system are moving with the same velocity. This is true for the dust, where the pressure is exactly zero. If $`P>0`$ such a situation may happen for some particular $`t`$ in the particular reference frame (remember we integrate $`T^{\mu \nu }`$ over the hypersurface $`t=const`$). If we boost our system the parts of matter on the new integration hypersurface $`t_1=const`$ will not move with the same velocity, as illustrated in Fig. 9: in the Lab frame all the matter on the hypersurface $`t=t_0`$ has the same velocity, but in the CM frame we have to either integrate over hypersurface $`t_1=t_1^0`$, where different parts of matter have different velocities or integrate over hypersurface $`t=t_0`$ homogeneous in velocity, but in this case $`dS_\mu (dV,\stackrel{}{0})`$ in eq. (45). Briefly speaking the presence of pressure requires, in general, the flow velocity field in CM frame, $`\stackrel{}{v}_{CM}(t,\stackrel{}{r})`$, for real physical situation, and our invariant energy-momentum tensor, $`T^{\mu \nu }`$, defined by eq. (48) seems to be clever enough to know that.

Thus, as we see, the assumption of the final streak moving like one object oversimplifies the real physical situation. The way to improve the model is to let the final streak expand with the rapidity field which will allow us to satisfy exactly all the conservation laws. We can describe such a streak expansion based on the analytical solution for one–dimensional expansion of matter into the vacuum.

### 3.5 Expanding final streaks

Apart of the principal reasons discussed above, our test results  indicated that the uniform final streaks without expansion lead to somewhat unrealistic expansion patterns. Thus we had to improve the model and include the expansion of the final streaks.

In this section we discuss the one-dimensional expansion of the finite streak into the vacuum (generalizing the description in ). The initial condition is

$`e(z,t_0)`$ $`=`$ $`\{\begin{array}{cc}0\hfill & \mathrm{}<z<R_{},\hfill \\ e_0\hfill & ,R_{}zR_+\hfill \\ 0\hfill & ,R_+<z<\mathrm{},\hfill \end{array}`$ (56)
$`v(z,t_0)`$ $`=`$ $`\{\begin{array}{cc}1\hfill & ,\mathrm{}<z<R_{}\hfill \\ \mathrm{tanh}y_0\hfill & ,R_{}zR_+\hfill \\ 1\hfill & ,R_+<z<\mathrm{},\hfill \end{array},`$ (60)
where $`R_{}`$, $`R_+`$ are the borders of the system. The choice $`v=\pm 1`$ in the vacuum is purely conventional, but it guarantees a continuous hydrodynamic solution at the boundary to the vacuum, since in the limit of infinite dilution the velocity of matter approaches unity.

Now consider a simple wave, i.e., continuous hydrodynamic flow. This flow respects entropy conservation . For isentropic flow, the $`t`$ and $`z`$component of the equation of motion, $`_\mu T^{\mu \nu }=0`$, can be combined to yield

$$\left(\frac{}{t}+\frac{v\pm c_0}{1\pm vc_0}\frac{}{z}\right)_\pm =0,$$
(61)
where the so-called Riemann invariants are

$$_\pm yy_0\pm _{e_0}^e\frac{c_0\mathrm{d}e^{}}{e^{}+p(e^{})}.$$
(62)
Obviously, the Riemann invariants are constant,

$$\frac{\mathrm{d}_\pm }{\mathrm{d}t}|_{𝒞_\pm }=0,$$
(63)
along the so-called characteristic curves $`𝒞_\pm (x,t)`$, the positions $`x_\pm (t)`$ of which are determined by integrating

$$\frac{\mathrm{d}x_\pm }{\mathrm{d}t}=\frac{v\pm c_0}{1\pm vc_0}.$$
(64)
The characteristics are the world-lines of sonic disturbances on the hydrodynamic flow pattern. For simple waves moving to the right (i.e., $`v>0`$), one can prove  that $`_+=const`$, and it suffices to consider the $`𝒞_{}`$–characteristic. For simple waves moving to the left ($`v<0`$), $`_{}=const`$, and only the $`𝒞_+`$–characteristics have to be considered.

For the initial conditions (56, 60) we determine the value of $`_\pm `$ at $`t=t_0,z=R_\pm `$ to be $`_\pm =0`$ and thus

$$y_\pm =y_0_{e_0}^e\frac{c_0\mathrm{d}e^{}}{e^{}+p(e^{})}.$$
(65)
If we consider the QGP EoS, $`p(e)=c_0^2e,c_0^2=const`$, this gives

$$v(e)_\pm =\mathrm{tanh}y(e)_pm=\mathrm{tanh}\left[y_0\frac{c_0}{1+c_0^2}\mathrm{ln}\left\{\frac{e}{e_0}\right\}\right]=$$

$$=\frac{\mathrm{tanh}y_0\frac{1(e/e_0)^{2c_0/(1+c_0^2)}}{1+(e/e_0)^{2c_0/(1+c_0^2)}}}{1\mathrm{tanh}y_0\frac{1(e/e_0)^{2c_0/(1+c_0^2)}}{1+(e/e_0)^{2c_0/(1+c_0^2)}}}.$$
(66)
This equation nicely shows the defining property of a simple wave, namely that there is a unique relationship between the value of the fluid velocity and its thermodynamic state . Furthermore, with the help of (64) for $`𝒞_{}`$, we can now calculate for any given $`e`$, $`0ee_0`$, the position $`z(t;e)`$ at which this value for the energy density occurs for given $`t`$,

$$z_\pm (t;e)R_{}=\frac{v(e)\pm c_0}{1\pm v(e)c_0}(tt_0).$$
(67)
Eq. (67) has similarity form, i.e., the profile of the rarefaction wave does not change with time when plotted as a function of the similarity variable $`\zeta _\pm \frac{zR_\pm }{tt_0}`$.

Let us now complete the solution of the hydrodynamic problem. Causality requires that the initial conditions (56, 60) remain unchanged for $`|\zeta _\pm |>1`$. Thus, we only have to determine the solution in the range $`1\zeta _\pm 1`$, i.e. in the forward light cone (remember that $`\zeta _+`$ ($`\zeta _{}`$), obtained from the $`𝒞_{}`$ ($`𝒞_+`$)–characteristic, eq. (67), corresponds to $`y_+`$ ($`y_{}`$), obtained from $`_+=const`$ ($`_{}=const`$) condition).

For the forward going rarefaction wave, $`v>0`$, generated on the right end of the initial streak, from eqs. (66, 67) we infer that the head of this wave (the point where the rarefaction of matter starts, i.e., where the energy density $`e`$ starts to fall below $`e_0`$) travels with the velocity $`\frac{\mathrm{tanh}y_0c_0}{1c_0\mathrm{tanh}y_0}`$ to the left. On the other hand, the base of the rarefaction wave (the point where the vacuum ends, i.e., where $`e`$ starts to acquire non-vanishing values) travels with light velocity $`v=1`$ to the right. Thus, the energy density as a function of $`\zeta `$ can be written as

$$e(\zeta _+)=e_0\{\begin{array}{c}1,\hfill \\ 1\zeta _+\frac{\mathrm{tanh}y_0c_0}{1c_0\mathrm{tanh}y_0}\hfill \\ \left[\frac{1+\mathrm{tanh}y_0}{1\mathrm{tanh}y_0}\frac{1c_0}{1+c_0}\frac{1\zeta _+}{1+\zeta _+}\right]^{(1+c_0^2)/2c_0},\hfill \\ \frac{\mathrm{tanh}y_0c_0}{1c_0\mathrm{tanh}y_0}<\zeta _+1.\hfill \end{array}$$
(68)
The velocity can then be inferred from (66), or simply from (67):

$$\mathrm{tanh}y(\zeta _+)=\{\begin{array}{cc}\mathrm{tanh}y_0\hfill & ,1\zeta _+\frac{\mathrm{tanh}y_0c_0}{1c_0\mathrm{tanh}y_0}\hfill \\ \frac{c_0+\zeta _+}{1+\zeta _+c_0}\hfill & ,\frac{\mathrm{tanh}y_0c_0}{1c_0\mathrm{tanh}y_0}<\zeta _+1.\hfill \end{array}$$
(69)
The similar analyzes for the backward going rarefaction wave, $`v<0`$, generated on the left end of the initial streak, gives:

$$e(\zeta _{})=e_0\{\begin{array}{c}1,\hfill \\ \frac{\mathrm{tanh}y_0+c_0}{1+c_0\mathrm{tanh}y_0}\zeta _{}1\hfill \\ \left[\frac{1\mathrm{tanh}y_0}{1+\mathrm{tanh}y_0}\frac{1c_0}{1+c_0}\frac{1+\zeta _{}}{1\zeta _{}}\right]^{(1+c_0^2)/2c_0},\hfill \\ 1<\zeta _{}\frac{\mathrm{tanh}y_0+c_0}{1+c_0\mathrm{tanh}y_0},\hfill \end{array}$$
(70)
$$\mathrm{tanh}y(\zeta _{})=\{\begin{array}{cc}\mathrm{tanh}y_0\hfill & ,\frac{\mathrm{tanh}y_0+c_0}{1+c_0\mathrm{tanh}y_0}\zeta _{}1\hfill \\ \frac{\zeta _{}c_0}{1\zeta _{}c_0}\hfill & ,1<\zeta _{}\frac{\mathrm{tanh}y_0+c_0}{1+c_0\mathrm{tanh}y_0}.\hfill \end{array}$$
(71)
This simple analytical solution is valid as long as the two rarefaction waves did not overlap in the middle of the system. Further evolution becomes more complicated and doesn’t have a similarity form any more.

The point where two rarefaction waves start to overlap, $`z_{over},t_{over}`$ can be found from the system

$$\{\begin{array}{c}\zeta _+(z,t)=\frac{zR_+}{tt_0}=\frac{\mathrm{tanh}y_0c_0}{1c_0\mathrm{tanh}y_0}\hfill \\ \zeta _{}(z,t)=\frac{zR_{}}{tt_0}=\frac{\mathrm{tanh}y_0+c_0}{1+c_0\mathrm{tanh}y_0},\hfill \end{array}$$
(72)
which gives

$$z_{over}=\begin{array}{c}R_+(1\mathrm{tanh}y_0c_0)(\mathrm{tanh}y_0+c_0)\\ \underset{¯}{R_{}(1+\mathrm{tanh}y_0c_0)(\mathrm{tanh}y_0c_0)}\\ 2c_0(1\mathrm{tanh}^2y_0)\end{array},$$
(73)
$$t_{over}=t_0+\frac{(z_{over}R_+)(1\mathrm{tanh}y_0c_0)}{(\mathrm{tanh}y_0c_0)}.$$
(74)
Thus the analytical solution (68-71) valid for $`t_0tt_{over}`$.

For the symmetric initial state, $`R_+=R_{}=R`$, eqs. (73,74) became

$$z_{over}=R\frac{\mathrm{tanh}y_0(1c_0^2)}{c_0(1\mathrm{tanh}^2y_0)},$$
(75)
$$t_{over}=t_0+R\frac{1\mathrm{tanh}y_0c_0}{c_0(1\mathrm{tanh}^2y_0)}.$$
(76)
Thus, based on the solution presented above we can have a more advanced description of the final streaks than the one presented in section 3.3. Let us assume that the homogenous final streaks, with some $`e_f`$, $`y_f`$, are formed in the CRF when the larger of initial streaks reach the rapidity $`y_f`$. Up to this point the fluid cell trajectories are the same as those discussed in section 3.3, but after the homogeneous final streak is formed, it starts to expand according to the simple rarefaction wave solutions presented above. Figure 7 (B) shows new trajectories of the streak ends, and Figures 7 (C,D) present energy density and rapidity profile for this expanding streak.

Such an initial state with expanding streaks will also help to avoid the problem which may cause the development of numerical artifacts, namely a step-like in beam direction initial energy density distribution (output of ESRM): it has a jump from $`e(x_i,y_j)=const`$ inside the matter to $`0`$ in the outside vacuum (of course, where is also a jump of $`e`$ as a function of $`x,y`$, but it is much smoother and this is not the direction of the initial expansion). In order to avoid (or at least to suppress) the effect of this, it was proposed in Ref.  to smooth over initial energy density distribution, for example by a Gaussian shape. Our simple analytic solution smooths over this jump in a natural way.

### 3.6 Initial conditions for hydrodynamical calculations

In this section we present the results of our calculations. We are interested in the shape of the QGP formed when the expansion of the combined target-projectile streaks stop and their matter is locally equilibrated. This will be the initial state for further hydrodynamical calculations. The time, $`t_h`$, at which we start the hydrodynamical description is a second (in addition to $`A`$) free parameter of our model. Of course, $`t_h`$ should be larger than the time of final streak formation, at least in the central, most hot and dense region. For the peripheral streaks the string tension is low, and the transparency is large, but peripheral matter does not play a leading role in further hydrodynamic expansion. Therefore, we will also build the final streaks ($`y_f,e_f`$) for peripheral streak-streak collisions, with lengths, $`l_f`$, corresponding to the lengths of the interacting region at the moment $`t=t_h`$, even if the final rapidity, $`y_f`$, was not yet achieved for this particular collision.444 The calculations of the initial state presented in this section are done for colliding nuclei with a homogeneous density distribution. The modification is pretty straightforward - we can use the Wood-Saxon density distribution for the Au nucleus in our case. This will lead to the increasing of the density in the middle of the colliding nuclei, therefore the middle steaks will contain more partons, will generate a stronger field and, thus, will form the final streaks faster. While the peripheral streaks will contain less matter, therefore their dynamics, which is oversimplified, will have less influence on the final state. On the other hand the time, $`t_h`$, cannot be longer than the time for which the simple analytical solution presented in the previous section is valid. I.e., $`t_hMinimum\{t_{form,i}+t_{over,i}\}`$, where $`t_{form,i}`$ is a time of formation of the final streak for particular streak-streak collisions $`i`$: when the larger of two initial streak gets rapidity $`y_f`$ (given by eq. (43) in Ref. ); $`t_{over,i}`$ is given by eq. (74). In this section we present all the results at $`t_h=5fm/c`$, while for $`\epsilon _0=100GeV/nucl`$, $`b=0.252R_{Au}`$, $`A=0.08`$ the maximal possible $`t_h=Minimum\{t_{form,i}+t_{over,i}\}5.04fm/c`$.

The reader may notice that all of the above discussion about final streaks is for the energy density distribution. What about the baryon charge distribution in the final streaks? Our knowledge here is very poor and, in principle, we are free to assume anything between two extreme cases: A) the complete transparency picture, where charges would remain at the ends of the final streak; this is the qualitative picture in Ref. ; B) homogeneous in baryon charge final streaks, the simplified situation discussed in Refs.

First, we present the initial rapidities for final streaks for different $`x`$ coordinates in the reaction plane – Fig. 10. Fig. 11 shows the energy density distribution in the reaction plane with non-expanding final streaks (Au+Au collision at $`\epsilon _0=65GeV/nucl`$, $`b=0.252R_{Au}`$, $`A=0.08`$). Fig. 12 shows the energy density distribution in the reaction plane for the same collision, but with expanding streaks. The QGP forms a tilted disk for $`b0`$. Thus, the direction of fastest expansion, the same as the largest pressure gradient, will be in the reaction plane, but will deviate from both the beam axis and the usual transverse flow direction. So, the new flow component, called “antiflow” or “third flow component” (see section 2.5), appears in addition to the usual transverse flow component in the reaction plane (for non-expanding streaks this was shown in ). With increasing beam energy the usual transverse flow is getting weaker, while this new flow component is strengthened.

More results for full RHIC energy ($`\epsilon _0=100GeV/nucl`$) for different string tensions and impact parameters are presented in Figs. 1314.

Note that the initial state with non-expanding streaks (Figs. 11) looks pretty much like the one discussed in “firestreak” model (see Fig. 4), since the final streaks are moving with rapidities presented in Fig. 10. If we let final streaks expand, this smooths over picture, but most of the matter, nevertheless, keeps similar energy density profile and velocity distributions.

From Figs. 11 and 13-14 we may see that for the central collisions at $`\epsilon _0=100GeV/nucl`$ the maximum energy density $`E=T^{00}`$ reaches as high as $`E_{max}100GeV/fm^3`$. It seems to be larger than what one should expect from the Bjorken estimate (see section 2.3) - $`e_0=4.3GeV/fm^3`$ for $`\sqrt{s}=130AGeV`$ . This is natural since our initial energy density distribution is much more peaked as a function of $`z`$ than is the one obtained in the Bjorken picture (our energy density distribution is also peaked as a function of transverse coordinate $`(x,y)`$, but here the difference seems to be not that large, see for example ). One should, nevertheless, keep in mind that we have also a huge initial pressure, $`P=c_0^2e`$, which will cause a fast hydrodynamical expansion, and thus, as discussed in section 2.3, will reduce $`\left(\frac{dE_T}{dy}\right)_{y=0}`$. This is a usual situation for hydrodynamical calculation, for example, in Ref.  $`e_0=23GeV/fm^3`$ has been used to reproduce the elliptic flow at RHIC for $`\epsilon _0=65GeV/nucl`$, although the Bjorken initial geometry had been assumed.

In a recent work  another scenario of the creation of the third flow component based on the shadowing picture was proposed. The calculations were performed in the framework of QGSM. The shadowing effect manifests itself in the following way - hadrons emitted at small rapidities in the antiflow direction can propagate freely, while hadrons emitted in the normal directed flow direction still remain within the expanding subsystem of interacting particles. We would like to point out how one can distinguish this scenario from the one discussed above in this paper. The tilted initial state has the maximal effect for non-zero, but not very large, impact parameters, approximately up to $`b=0.5(R_A+R_B)`$, i.e., in semicentral reactions. For semiperipheral and peripheral collisions our string rope tension, $`\sigma `$, would be small, and we would have to wait for a long time for the final streak formation, i.e., at the initial stages we would have one-dimensional Bjorken transparency. On the other hand, shadowing just becomes stronger with increasing impact parameter and gives the most drastic effect for semiperipheral and peripheral collisions.

## 4 Conclusions

A realistic and detailed description of an energetic heavy ion reaction requires a Multi Module Model, wherein the different stages of the reaction are each described with a suitable theoretical approach. The initial stages are the most problematic at RHIC energies. We tried to show that for the qualitative understanding of the basic features of the initial stages a phenomenological “firestreak” model can be used. For the detailed calculations of the initial state the Effective String Rope model has been developed. This model is based on earlier Coherent Yang-Mills field theoretical models and introduces an effective string tension based on Monte-Carlo string cascade and parton cascade model results. Contrary to earlier expectations based on standard string tensions of $`1GeV/fm`$, which lead to the Bjorken model type of initial state, the effective string tension introduced in our model for collisions of energetic massive heavy ions causes limited transparency. The increased string tension is a consequence of collective effects related to QGP formation. These collective effects in central and semicentral collisions lead to an effective string tension of the order of $`10GeV/fm`$ and consequently cause much less transparency than earlier estimates. The resulting initial locally equilibrated state of matter in semicentral collisions takes a rather unusual form (see Figs. 11, 12, 13-14) and the resulting flow velocity field (see Fig. 10, 7 (D)) is in agreement with the “firestreak” phenomenological picture. Such an initial state will manifest itself by the asymmetry of the produced collective flow, even after the subsequent hydrodynamical expansion .

## Acknowledgment

Two of us (V.K.M and D.D.S.) thank the support of the Bergen Computational Physics Laboratory in the framework of the European Community - Access to Research Infrastructure action of the Improving Human Potential Programme and the Humboldt Foundation.

## Appendix A Initial conditions after string creation

Our conserved quantities are (36, 37)

$$Q^0=\stackrel{~}{T}^{00}𝑑V=xy\stackrel{~}{T}^{00}𝑑z,$$
(77)
$$Q^3=\stackrel{~}{T}^{03}𝑑V=xy\stackrel{~}{T}^{03}𝑑z,$$
(78)
where $`\stackrel{~}{T}^{00}`$ and $`\stackrel{~}{T}^{03}`$ are given by eq. (38). Before string creation the initial values of the modified energy-momentum tensor, $`\stackrel{~}{T}^{\mu \nu }`$, are

$$\stackrel{~}{T}_1^{00}=\stackrel{~}{T}_2^{00}=e_0\mathrm{cosh}^2y_0=\left(\frac{\epsilon _0}{m}\right)^2n_0m,$$
(79)
$$\stackrel{~}{T}_2^{03}=\stackrel{~}{T}_1^{03}=e_0\mathrm{tanh}y_0\mathrm{cosh}^2y_0=\left(\frac{\epsilon _0}{m}\right)^2n_0mv_0$$
(80)
where $`m`$ is the nucleon mass, $`\epsilon _0`$ is the initial energy per nucleon, and we have used $`\mathrm{cosh}^2y_0=\gamma _0^2=\left(\frac{\epsilon _0}{m}\right)^2`$. $`v_0=\mathrm{tanh}y_0`$ is the initial velocity, $`v_0=1`$ is a good approximation for ultra-relativistic heavy ion collisions. Thus,

$$Q^0=xy\left(\frac{\epsilon _0}{m}\right)^2n_0m(l_1+l_2),$$
(81)
$$Q^3=xy\left(\frac{\epsilon _0}{m}\right)^2n_0m(l_2l_1)v_0,$$
(82)
where $`l_1`$ and $`l_2`$ are the initial lengths of streaks (see Fig. 5), $`x`$, $`y`$ are the grid sizes in $`x`$ and $`y`$ directions. Thus the rapidity of CM frame is

$$\mathrm{tanh}y=\frac{Q^3}{Q^0}=\frac{v_0}{M},$$
(83)
where

$$M=\frac{l_2+l_1}{l_2l_1}.$$
(84)
After string creation –

$$\stackrel{~}{T}^{00}=e_1\mathrm{cosh}^2y_1+c_0^2e_1\mathrm{sinh}^2y_1+e_2\mathrm{cosh}^2y_2+c_0^2e_2\mathrm{sinh}^2y_2$$

$$+\frac{1}{2}\sigma ^2+B\frac{\sigma x^+}{2}n_0e^{y_0}e^{2y_1}\frac{\sigma x^{}}{2}n_0e^{y_0}e^{2y_2},$$
(85)
$$\stackrel{~}{T}^{03}=e_1(1+c_0^2)\mathrm{cosh}y_1\mathrm{sinh}y_1+e_2(1+c_0^2)\mathrm{cosh}y_2\mathrm{sinh}y_2$$

$$\frac{\sigma x^+}{2}n_0e^{y_0}e^{2y_1}+\frac{\sigma x^{}}{2}n_0e^{y_0}e^{2y_2},$$
(86)
At the point of complete penetration of streaks, $`t=t_0=(l_1+l_2)/2`$ (see Fig 5), we introduced energy densities $`e_1(t_0)`$ and $`e_2(t_0)`$. We assumed transparency, i.e., that complete penetration happened so fast, that the field itself that was created during this time did not have time to stop partons. So, the rapidities $`y_{1(2)}(t_0)=y_0(y_0)`$ correspondingly, and the proper baryon densities did not change, and thus, the baryon current is conserved. Then the energy and momentum conservation laws can be written in the form:

$$\frac{Q^0}{xy}=\left[(1+c_0^2)\mathrm{cosh}^2y_0c_0^2\right]\left(e_1(t_0)l_1+e_2(t_0)l_2\right)$$

$$+\left(\frac{\sigma ^2}{2}+B\right)(l_1+l_2)+\frac{\sigma n_0e^{y_0}}{4}\left(l_1^2+l_2^2\right),$$
(87)
$$\frac{Q^3}{xy}=\left[(1+c_0^2)\mathrm{cosh}^2y_0\right]\left(e_1(t_0)l_1+e_2(t_0)l_2\right)$$

$$\frac{\sigma n_0e^{y_0}}{4}\left(l_1^2l_2^2\right).$$
(88)
We neglect $`c_0^2`$ next to $`(1+c_0^2)\mathrm{cosh}^2y_0`$ in eq. (87), then eqs. (87, 88) may be solved

$$e_1(t_0)=\frac{n_0m}{1+c_0^2}\frac{\frac{\sigma ^2}{2}+B}{\left(\frac{\epsilon _0}{m}\right)^2(1+c_0^2)}\frac{l_1+l_2}{2l_1}\frac{\sigma n_0e^{y_0}}{4\left(\frac{\epsilon _0}{m}\right)^2(1+c_0^2)}l_1,$$
(89)
$$e_2(t_0)=\frac{n_0m}{1+c_0^2}\frac{\frac{\sigma ^2}{2}+B}{\left(\frac{\epsilon _0}{m}\right)^2(1+c_0^2)}\frac{l_1+l_2}{2l_2}\frac{\sigma n_0e^{y_0}}{4\left(\frac{\epsilon _0}{m}\right)^2(1+c_0^2)}l_2.$$
(90)


format: text

text:
# Kondo impurities in nanoscopic systems: new confinement-induced regimes

\[

## Abstract

We present results for Kondo impurities in nanoscopic systems. Using exact diagonalization in small clusters and Wilson’s renormalization group we analyze an isolated system and a nanoscopic system weakly coupled to a macroscopic reservoir. In the latter case, new regimes not observed in macroscopic homogeneous systems are induced by the confinement of conduction electrons. These new confinement-induced regimes are very different depending on whether the Fermi energy coincides with the energy of a resonant state or lies between two quasi-bound states.

\]

The physics of nanoscopic and mesoscopic systems became of grate interest due to the possibility of making experiments in extremely small samples. Some examples of nano and mesoscopic systems recently studied are small metallic and superconducting islands, quantum dots, nanotubes, and quantum corrals. Depending on the type of system under study, different parameters can be controlled.

The advances in nanotechnologies revived the interest in the Kondo effect, one of the paradigms of strongly correlated systems. On one hand, Scanning Tunneling Microscopy (STM) allowed the direct measurement of local spectroscopic properties of Kondo impurities. On the other hand it has been shown that quantum dots and single walled carbon nanotubes weakly coupled to contacts may behave as Kondo impurities generating new alternatives to study the phenomena.

In this work we report results for Kondo impurities in nanoscopic systems. Some experiments show that the effects of confinement are important and recently Affleck and Simon proposed to use a mesoscopic system to measure the Kondo screening length. Although the problem has been addressed theoretically, there are still many open questions. What is the temperature dependence of the Kondo screening in nanoscopic systems? What is the effect adding leads or coupling the system to a macroscopic reservoir?

In the case of systems weakly coupled to a macroscopic reservoir, the one electron states acquire a finite lifetime and the spectrum consists of resonant states. This situation is the most relevant to compare with experiments, however due to the intrinsic difficulties to formulate the problem is has not been analyzed in detail. As an example we may mention the case of quantum corrals, were the electrons are not well confined in the corral due to leaking of the wave function and hybridization with the bulk states. Consequently the number of electrons in the corral is not a good quantum number and the local density of states as a function of the energy presents well defined maxima corresponding to the quasi-confined states of the corral. We show that in systems like the quantum corrals new regimes are induced by the confinement of conduction electrons. These new confinement-induced regimes are very different if the Fermi energy is at resonance or between two quasi-bound states. We also discuss the case of quantum dots in Aharonov-Bohm rings.

We have performed numerical diagonalization in small clusters and implemented Wilson’s numerical renormalization group (NRG). Our starting point is the Anderson model for magnetic impurities. The model Hamiltonian reads:

$`H_{AM}`$ $`=`$ $`{\displaystyle \underset{\sigma }{}}\epsilon _dd_\sigma ^{}d_\sigma +Ud_{}^{}d_{}d_{}^{}d_{}+{\displaystyle \underset{\nu ,\sigma }{}}\epsilon _\nu c_{\nu \sigma }^{}c_{\nu \sigma }`$ (2)
$`+{\displaystyle \underset{\nu ,\sigma }{}}V_\nu c_{\nu \sigma }^{}d_\sigma +V_\nu ^{}d_\sigma ^{}c_{\nu \sigma }\mu _iBS_{iz},`$
here the operator $`d_\sigma ^{}`$ creates an electron with spin $`\sigma `$ at the impurity orbital with energy $`\epsilon _d`$ and Coulomb repulsion $`U`$, $`c_{\nu \sigma }^{}`$ creates an electron in an extended state with quantum numbers $`\nu `$ and $`\sigma `$ and energy $`\epsilon _\nu `$. The last term represents the effect of an external magnetic field along the $`z`$-direction coupled to the impurity spin $`𝐒_i=d_\sigma ^{}\tau _{\sigma ,\sigma ^{}}d_\sigma ^{}`$ where $`\tau `$ are the Pauli matrices.

We start by analyzing the case of an isolated system. The properties of the system depend crucially on a few parameters, in particular we show below that in the Kondo limit a finite system presents quite a marked even-odd electron number asymmetry.

We first consider the case of Hamiltonian (1) describing a linear chain of $`N`$ equivalent sites with the impurity at one end. In this case $`ϵ_\nu =2t\mathrm{cos}(\nu \pi /(N+1))`$, where $`t`$ is the hopping matrix element in the chain and $`V_\nu =V_0\mathrm{sin}(\nu \pi /(N+1))`$ with $`\nu =1,2,\mathrm{},N`$. We performed exact diagonalization in these finite systems using a Lanczos algorithm.

For an even number of electrons in the system, the ground state is a singlet indicated as $`|0`$ and the expectation value of the impurity spin $`0|S_{iz}|0`$ is zero. The impurity spin is screened by the free electrons and the zero temperature susceptibility is finite. We can estimate the characteristic energy scale $`k_BT_K`$ as the energy involved in the screening of the impurity spin. The magnetization as a function of the external field $`B`$, shown in the inset of Fig.1(a), gives an estimation of $`T_K`$. We define $`k_BT_K=\mu _iB_c`$ where $`B_c`$ is the crossover field indicated in the inset. As $`t`$ increases the mean density of conduction states decreases and consequently $`T_K`$ decreases. As indicated in Fig.1(a) two regimes are obtained: $`k_BT_K<\mathrm{\Delta }`$ and $`k_BT_K`$ $`>\mathrm{\Delta }`$ where $`\mathrm{\Delta }4t/N`$ is the level spacing. For the parameters used the crossover between these two regimes takes place at a hopping $`t_{cross}0.27`$.

For an odd number of electrons, the ground state is a Kramers spin-$`1/2`$ doublet, $`|`$ and $`|`$. The expectation value $`s|S_{iz}|=|S_{iz}|`$ is different from zero and, in the low temperature limit the impurity susceptibility diverges as $`\chi =\mu _i^2s^2/k_BT`$. In Fig.1(b) the effective magnetic moment $`\mu =\mu _is`$ (we take $`\mu _i=1`$) is plotted as function of $`t`$. For large $`\mathrm{\Delta }`$ (large $`t`$) the impurity spin is unscreened and $`s1/2`$ as expected for a gaped system; for small $`\mathrm{\Delta }`$ the impurity spin tends to be completely screened. The crossover between these two regimes takes place at $`t_{cross}`$.

Although these results were obtained for a very small system, they give the correct physical picture for zero temperature. While for $`\mathrm{\Delta }k_BT_K`$ the impurity spin is essentially screened regardless of the parity of the electron number, for $`\mathrm{\Delta }k_BT_K`$ an important even–odd asymmetry is obtained. In order to analyze larger systems we have to resort to some approximate scheme. We have adapted the NRG to the present problem. The idea of the NRG is to obtain the quantum many-body states of the system on all energy scales, or length scales, in a sequence of steps. To do that, Wilson defined a basis of concentric wavefunctions in which the Hamiltonian takes the form of a linear chain. A truncated chain with $`n`$ sites, described by an effective Hamiltonian $`H_n`$, gives the correct physics on an energy scale $`w_n`$. The NRG transformation $`𝐓`$ relates the Hamiltonians describing successive lower energy scales, $`H_{n+1}=𝐓[H_n]`$, and leads to a systematic way of calculating the thermodynamic properties at successive lower temperatures.

Going back to our problem and for the sake of simplicity, consider a spherical metallic cluster of radius $`R_c`$ with the impurity at the center. The cluster is embedded in a bulk material with which it is weakly coupled through a large surface barrier. The Hamiltonian can be put in the form $`H=H_{AM}+V(r,R_c)`$ where the last term is a potential barrier placed at a distance $`R_c`$ from the impurity. The Hamiltonian $`H_{AM}`$ is rewritten in Wilson’s basis as schematically shown in Fig.2. The potential barrier is described by including a higher diagonal energy to the $`(n_c+1)`$–shell centered around $`R_c`$. Alternatively the barrier can be simulated with a smaller hopping matrix element in the region where $`V(r,R_c)`$ is different from zero. We adopted this last description reducing the hoppings by a factor $`\alpha `$.

If the barrier is impenetrable ($`\alpha =0`$), the central cluster is isolated. The process of renormalization ends up after $`n_c`$ steps and with the obtained low energy states with a definite number of particles $`N_e`$ the partition function $`Z`$ is evaluated and the thermodynamic properties of the system are obtained. This procedure shows that the high temperature properties, evaluated with a number of shells $`n<n_c`$, the impurity behaves as in a bulk material. Only when the temperature is lowered and the discrete nature of the one-electron states becomes evident the behavior of the system deviates from that of the macroscopic system. The impurity susceptibility can be evaluated for even and odd number of particles. As we discuss below, to describe even or odd number of particles it is convenient to take even or odd number of shells $`n_c`$. The susceptibility is calculated as:

$`\chi =\mu _i^2({\displaystyle \underset{\nu }{}}{\displaystyle \frac{P_\nu |\nu |S_{iz}|\nu |^2}{k_BT}}+{\displaystyle \underset{\nu \xi }{}}{\displaystyle \frac{2P_\nu |\nu |S_{iz}|\xi |^2}{E_\xi E_\nu }})`$

where the summation is done over the low energy states $`|\nu `$ with energies $`E_\nu `$ and $`P_\nu =\mathrm{exp}(E_\nu /k_BT)/Z`$. The matrix elements $`\nu |S_{iz}|\xi `$ have to be evaluated in a recursive way at each renormalization step. The susceptibility reflects the thermodynamic properties of the system and we use it as an indicator of the degree of screening of the impurity spin. In Fig.3 the impurity susceptibility is shown for different system sizes. For comparison the same quantity evaluated in the thermodynamic limit ($`n_c\mathrm{}`$) is shown. As $`n_c`$ decreases the characteristic energy separation between the one–electron states $`\mathrm{\Delta }`$ increases. In the NRG this energy separation at the Fermi level is $`\mathrm{\Delta }D\mathrm{\Lambda }^{n_c/2}`$ where $`D`$ is half the band width and $`\mathrm{\Lambda }2`$ is Wilson’s discretization parameter.

In Fig.3(a) the effective impurity magnetic moment squared $`\mu ^2=k_BT\chi `$ for an even number of electrons is shown as a function of the temperature. When the temperature is of the order of $`\mathrm{\Delta }`$ the behavior deviates from that of the macroscopic sample and $`\mu ^2`$ decreases exponentially due to the gaps induced by confinement. For an odd number of electrons the behavior is illustrated in Fig.3(b). In this case, as $`k_BT`$ approaches $`\mathrm{\Delta }`$ the screening of the impurity spin ends up and a finite magnetic moment survives down to zero temperature, the ground state is a spin-1/2 doublet. In the inset the low temperature magnetic moment versus $`\mathrm{log}(R_c/\xi _K^{\mathrm{}})`$ is shown with $`\xi _K^{\mathrm{}}`$ being the Kondo screening length (see below), the behavior shows the same tendency of that of Fig.1(b) corresponding to exact results of the impurity in a linear chain. Here the crossover between an almost unscreened impurity spin to an almost completely screened spin occurs for $`\mathrm{\Delta }`$ $`k_BT_K^{\mathrm{}}`$.

The behavior of the system with a finite barrier is richer and more interesting. Now the central cluster is in contact with a reservoir and the number of electrons in the cluster is not a good quantum number. In the numerical NRG approach the Fermi energy is set to zero. The one particle spectrum in the absence of the impurity has a state at zero energy for even iterations $`n`$ and two states at approximately $`\pm D\mathrm{\Lambda }^{n/2}`$ for odd $`n`$. Then by taking the barrier enclosing an even or odd number of shells, $`n_c`$, we describe a system with the Fermi level at the energy of a resonant state (at resonance) or between two of them (off resonance) respectively. These two situations correspond approximately to an even and odd average number of electrons in the cluster respectively.

In Fig.4(a) the impurity magnetic moment $`\mu ^2`$ is shown for the a system at resonance. Only when the temperature is lowered and reaches the value $`k_BT\mathrm{\Delta }`$ the structure in the local density of states becomes important. At low temperatures the universal Kondo behavior is recovered with a Kondo temperature $`T_K^R`$ higher than $`T_K^{\mathrm{}}`$, corresponding to a higher density of states. The crossover between the high and low temperature regimes occurs in a small temperature interval leading to a rapid decrease of the magnetic moment $`\mu ^2`$.

For the system off resonance, as $`k_BT`$ approaches $`\mathrm{\Delta }`$ the magnetic moment saturates leading to a plateau in the curve $`\mu ^2`$ versus $`T`$ as shown in Fig.4(b). Only at lower temperatures the universal Kondo behavior is recovered with a Kondo temperature $`T_K^{OR}`$ that is lower than $`T_K^{\mathrm{}}`$. The plateau corresponds to a new regime of partially screened magnetic moment. In this regime, as the temperature is lowered, the thermodynamic properties do not change and are similar to those of a finite system with an even number of particles. The central cluster has a magnetic moment and it is the whole cluster that behaves as a Kondo impurity or a “Kondo grain” when the temperature is lowered. For a given potential barrier height the size of the grain can be changed: systems with small $`R_c`$, a large $`\mathrm{\Delta }`$, present a wide plateau with a partially-screened magnetic moment $`\mu _{ps}^21/4`$. As $`R_c`$ increases, both the width $`\delta T`$ of the plateau where the partially-screened magnetic moment $`\mu _{ps}`$ is stable and the value of $`\mu _{ps}^2`$ decrease. For a given size, the $`\delta T`$ increases with the barrier height as illustrated in Fig.4(c).

In summary, for $`k_BT_K^{\mathrm{}}\mathrm{\Delta }`$ the fine structure (on the scale of $`T_K^{\mathrm{}}`$) of the density of states does not change the properties of the system and the universal Kondo behavior is recovered. For $`k_BT_K^{\mathrm{}}\mathrm{\Delta }`$ the local density of states, on the scale of $`T_K^{\mathrm{}}`$, is a smooth function and the universal Kondo behavior is obtained at low temperatures. In this case the system at resonance has a Kondo temperature $`T_K^R`$ much larger than the corresponding to the system out of resonance $`T_K^{OR}`$, however in both cases the high temperature susceptibility is non universal. For $`k_BT_K^{\mathrm{}}\mathrm{\Delta }`$ new confinement induced regimes are observed: for the system at resonance, as the temperature is lowered, there is a rapid decrease in the magnetic moment $`\mu ^2`$ when $`k_BT\mathrm{\Delta }`$; for the system off resonance there is a new partially-screened regime that leads to a plateau in the temperature dependence of $`\mu ^2`$. This plateau can be interpreted as the free-spin regime of the whole Kondo grain that has magnetic moment partially localized at the impurity. At very low temperatures the electrons that are outside the cluster complete the screening. The temperature at which this occurs depends on the coupling between nanoscopic system and the bath. The condition for the existence of these new regimes $`k_BT_K^{\mathrm{}}\mathrm{\Delta }`$ can be put in the form $`\xi _K^{\mathrm{}}R_c`$ where $`\xi _K^{\mathrm{}}=\mathrm{}v_F/k_BT_K^{\mathrm{}}`$ is the Kondo screening length.

For the case of Co impurities in quantum corrals built on the Cu (111) surface, taking $`T_K^{\mathrm{}}=50K`$ and using the parameters of the Cu (111) surface band we estimate $`\xi _K^{\mathrm{}}400\AA `$ which is larger than the typical corral radius. Consequently if the Fermi level is out of resonance we expect the Kondo screening to occur only at very low temperatures. Larger corrals could be built or the effective barrier height decreased by removing some atoms from the corral fence and in this way control the effective Kondo temperature of the impurity.

Our results may have important consequences for the interpretation of experiments with QD in Aharonov-Bohm (AB) rings. In the AB configuration, the local density of states inside the ring has the structure of resonant states. A direct consequence of this structure are the oscillations in the transmittance of a ring as a magnetic field threading the circuit is varied. For rings built on GaAs as in some of the experimental setups the energy separation between resonances $`\mathrm{\Delta }`$ could be larger than one Kelvin, in fact at 1K oscillations are still observed indicating that at this temperature thermodynamic fluctuations do not wash out the density of states structure. If the Kondo temperature of a quantum dot inserted in one of the ring arms is of the order or smaller than $`\mathrm{\Delta }`$, new effects are to be taken into account. In this case the quantum dot behaves quite differently if the Fermi level lies at resonance or off resonance, corresponding to a maximum or a minimum in the transmittance respectively. The magnetic field can continuously shift the resonances and if the system is in the temperature interval $`T_K^{OR}`$ $`<`$$`T`$ $`<`$$`T_K^R`$ , even at constant temperature the QD continuously evolves from a low temperature regime ($`T`$ $`<`$$`T_K^R`$) to a high temperature regime ($`T_K^{OR}`$ $`<`$$`T`$ ). Moreover, the energy of the resonances depend not only on the external field but also on the phase shift introduced by the QD that is selfconsistently adjusted. On the light of our analysis, in this regime the AB oscillations can not be interpreted directly in terms of a Kondo impurity at a fixed temperature ratio $`T/T_K`$.

Acknowledgment: We thank I. Affleck and N. Andrei for stimulating discussions. This work was partially supported by the CONICET and ANPCYT, grants N. 02151 and 99 3-6343.