text:
# Thermal and electrical transport in single crystalline MgB2

## Abstract

We present the results of measurements of the thermal and electrical conductivity in the basal plane of single-crystalline MgB2 between 2 and 300 K. The analysis of the temperature dependence of the thermal conductivity gives supporting evidence for two different gaps on different sheets of the Fermi surface in the superconducting state. The zero-temperature values of the two gaps are $`\mathrm{\Delta }_1(0)=5.3\pm 0.5`$ meV and $`\mathrm{\Delta }_2(0)=1.65\pm 0.2`$ meV.

The recent discovery of superconductivity in MgB2 below an unexpectedly high critical temperature $`T_c`$ of the order of 40 K initiated intensive studies of its physical properties Nagamatsu01 . Numerous results indicate that the superconducting state of MgB2 is well described by the original BCS theory. Most experiments are compatible with a nodeless superconducting order parameter and imply that the electron-phonon interaction is responsible for the pairing of the electrons Buzea01cm . However, various types of experiments (see Ref. Bouquet01cm  and references therein), mainly using powder or polycrystalline samples, have given evidence for two gaps in the quasiparticle excitation spectrum of this superconductor. Qualitative and especially quantitative experimental checks of this intriguing situation on single crystalline material seem in order.

Below we present results of measurements of the thermal conductivity $`\kappa `$ and the electrical resistivity $`\rho `$ parallel to the basal $`ab`$-plane of the hexagonal crystal lattice of MgB2 as a function of temperature $`T`$ between 2 and 300 K. We demonstrate that our low-temperature results cannot be explained in terms of a single-gap function below the critical temperature, but are well approximated by assuming two different energy gaps in the superconducting state.

The investigated single crystal with dimensions of $`0.5\times 0.17\times 0.035`$ mm3 was grown with a high-pressure cubic anvil technique as described elsewhere Karpinski01 . The resistivity was measured by employing a 4-contact configuration in the $`ab`$-plane. For the thermal conductivity measurements, the standard method with a constant uniaxial heat flow was used. The same contacts were used for measuring the voltage and the temperature difference, respectively. Additional measurements of $`\rho (T)`$ and $`\kappa (T)`$ in constant magnetic fields $`H`$, oriented along the hexagonal $`c`$-axis, were made as well.

The electrical resistivity $`\rho (T)`$ in zero magnetic field is presented in Fig. 1.

The narrow ($`\mathrm{\Delta }T_c=0.1`$ K) superconducting transition occurs at $`T_c=38.1`$ K. In the temperature region between $`T_c`$ and 130 K, $`\rho (T)`$ may be well approximated by $`\rho =\rho _0+A^{}T^3`$, where $`\rho _0`$ and $`A^{}`$ are constants. At temperatures below 50 K, $`\rho \rho _0=2.0`$ $`\mu \mathrm{\Omega }`$ cm. A cubic $`\rho (T)`$ dependence has often been observed in transition metals and is associated with interband electron-phonon scatteringMott35 ; Wilson38 . The corresponding equation

$$\rho (T)=\rho _0+AT^3J_3\left(\mathrm{\Theta }_R/T\right),$$
(1)
where $`J_n(z)=_0^zx^ne^x\left(1e^x\right)^2𝑑x`$, fits our data perfectly. This is demonstrated by the solid line in Fig. 1. A similar fit to data from polycrystalline material is presented in Ref. Putti01cm . The fit parameters $`A`$ and $`\mathrm{\Theta }_R`$ characterize the strength of electron-phonon interaction and the cut-off frequency $`\omega _{\mathrm{max}}`$ of the phonon spectrum, respectively. The fit value of $`\mathrm{\Theta }_R=1226`$ K is in excellent agreement with the value of $`\mathrm{}\omega _{\mathrm{max}}=105`$ meV measured in inelastic neutron scattering experiments Clementyev01 ; Muranaka01\_INS . The magnitude of $`\mathrm{\Theta }_R`$ is considerably higher than the Debye temperature $`\mathrm{\Theta }_D700900`$ K, extracted from low-temperature specific heat measurements Budko01 ; Waelti01 ; Junod01cm . This discrepancy is not unexpected because electrical resistivity and specific heat involve different averages of the lattice frequencies Blackman51 .

The inset of Fig. 1 emphasizes $`\rho (T)`$ close to the superconducting transitions in fields of $`H=0`$ and $`H=50`$ kOe. For $`H=50`$ kOe, the critical temperature, marked by the onset of the deviation from the constant resistivity in the normal state, is shifted to $`T_c=17`$ K.

The thermal conductivity data $`\kappa (T)`$ in zero magnetic field and for $`H=`$ 0.62 and 50 kOe are presented in Fig. 2. The zero-field $`\kappa (T)`$ values are about an order of magnitude higher than previously reported for polycrystalline samples Muranaka01 ; Bauer01 ; Putti01cm ; Schneider01 . Also the temperature dependence of $`\kappa `$ is quite different from those earlier data. Instead of a monotonous increase with temperature we note a distinct maximum of $`\kappa (T)`$ at $`T65`$ K. No anomaly in $`\kappa (T)`$ provides evidence for the transition at $`T_c`$. A pronounced slope change in $`\kappa (T)`$ is observed around $`T6`$ K, however.

The thermal conduction of a superconductor is usually provided by electrons ($`\kappa _e`$) and phonons ($`\kappa _{\mathrm{ph}}`$), such that

$$\kappa =\kappa _e+\kappa _{\mathrm{ph}}.$$
(2)
A separate identification of these two terms at arbitrary temperature is not straightforward. In the normal state, a convenient and often used way to estimate $`\kappa _e`$ is to employ the Wiedemann-Franz law, which relates the electrical resistivity and the electronic contribution to the thermal conductivity via

$$\kappa _e(T)=L_0T/\rho (T),$$
(3)
where $`L_0=2.45\times 10^8`$$`\mathrm{\Omega }`$ K-2 is the Lorenz number. The validity of this law requires an elastic scattering of electrons and it is well established that Eq. (3) is applicable if the scattering of electrons by defects dominates. This is usually true at low temperatures, where $`\rho (T)\rho _0`$. At higher temperatures but below $`\mathrm{\Theta }_D`$, Eq. (3) is usually invalid. In our case, the regime of applicability of Eq. (3) is limited to temperatures below about 50 K (see Fig. 1). At temperatures $`T_cT50`$ K, $`\kappa _e(T)`$ provides about half of the total thermal conductivity.

Upon decreasing the temperature below $`T_c`$, the reduction of the number of unpaired electrons leads to a decrease of $`\kappa _e`$ and an increasing $`\kappa _{\mathrm{ph}}`$. The overall behavior of $`\kappa (T)`$ in the superconducting state depends on the relative magnitudes of $`\kappa _e`$ and $`\kappa _{\mathrm{ph}}`$ and also on the strength of electron-phonon interaction. Using existing theoretical models, $`\kappa (T)`$ data may provide valuable information about the superconducting gap function. In the following, we restrict our analysis of $`\kappa (T)`$ to temperatures below 50 K, where the separation of $`\kappa _e`$ and $`\kappa _{\mathrm{ph}}`$ can reliably be accomplished.

The electronic contribution to $`\kappa (T)`$ in the superconducting state has been calculated by Bardeen, Rickayzen, and Tewordt Bardeen59 . In their model

$$\kappa _e=\kappa _{e,n}f(y),$$
(4)
where

$`f(y)`$ $`=`$ $`{\displaystyle \frac{2F_1(y)+2y\mathrm{ln}(1+e^y)+\frac{y^2}{1+e^y}}{2F_1(0)}},(T<T_c)`$
$`\mathrm{and}`$ (5)
$`f(y)`$ $`=`$ $`1,(TT_c),`$
as well as $`F_n(y)=_0^{\mathrm{}}z^n(1+e^{z+y})^1𝑑z`$, $`y=\mathrm{\Delta }(T)/k_BT`$, and $`\mathrm{\Delta }(T)`$ representing the energy gap. Here $`\kappa _{e,n}`$ is the normal-state electronic thermal conductivity and in our analysis $`\kappa _{e,n}=L_0T/\rho _0`$.

The phonon thermal conductivity was analyzed in terms of the Debye-type relaxation rate approximation. The corresponding equation is BermanBook

$$\kappa _{\mathrm{ph}}=\frac{k_\mathrm{B}}{2\pi ^2v}\left(\frac{k_\mathrm{B}}{\mathrm{}}\right)^3T^3\underset{0}{\overset{\mathrm{\Theta }_D/T}{}}\frac{x^4e^x}{(e^x1)^2}\tau (\omega ,T)𝑑x,$$
(6)
where $`\omega `$ is the frequency of a phonon, $`\tau (\omega ,T)`$ is the corresponding relaxation time, $`v=\mathrm{\Theta }_D\left(k_B\mathrm{/}\mathrm{}\right)(6\pi ^2n)^{1/3}`$ is the average sound velocity, $`n`$ is the number density of atoms, and $`x=\mathrm{}\omega /k_\mathrm{B}T`$. For our analysis, we use $`\mathrm{\Theta }_D=750`$ K, as deduced from specific heat measurements Waelti01 ; Budko01 . The total phonon relaxation rate

$`\tau ^1=L/v+B\omega ^4+CT\omega ^2\mathrm{exp}(\mathrm{\Theta }_D/bT)+D\omega `$
$`+E\omega g(x,y)`$ (7)
can be represented as a sum of terms corresponding to independent scattering mechanisms. The individual terms, dominating in different temperature intervals, introduce phonon scattering by sample boundaries, point defects, phonons, dislocations, and electrons, respectively. The constants $`L`$, $`B`$, $`C`$, $`b`$, $`D`$, and $`E`$ are a measure for the intensity of corresponding phonon relaxation processes. The function $`g(x,y)`$ is given in Ref. Bardeen59 ; we do not reproduce its rather complex form here.

For $`H=50`$ kOe the phonon contribution $`\kappa _{\mathrm{ph}}(T)`$ in the normal state (above 17 K), extracted from $`\kappa (T)`$ and $`\rho (T)`$ using Eqs. (2) and (3), is identical to $`\kappa _{\mathrm{ph}}(T)`$ calculated for $`H=0`$ above $`T_c`$. We thus conclude that magnetic field effects on $`\kappa _{\mathrm{ph}}`$ are negligible. Based on the assumption that phonon scattering by both defects and phonons is the same in the normal and in the superconducting state, we estimate the temperature region where the number of unpaired electrons is too small to produce significant contributions to both phonon scattering and heat transport. Taking into account the lowest previously claimed value of $`\mathrm{\Delta }(0)=1.7`$ meV Tsuda01 , the electronic quasiparticles are effectively excluded from heat transport below 4 K.

The following analysis procedure is based on the assumptions discussed above. In a first step, the values of the parameters $`L`$, $`B`$, $`C`$, $`b`$, $`D`$, and $`E`$, related to the processes of phonon relaxation in the normal state, were obtained by fitting the $`\kappa _{\mathrm{ph}}(T)`$ data in $`H=50`$ kOe above 17 K (for this region $`g(x,y)=1`$) to Eqs. (6) and (Thermal and electrical transport in single crystalline MgB2) as well as $`\kappa (T)`$ below 4 K where $`g(x,y)=0`$ and $`\kappa _e=0`$. These fits are shown in the inset of Fig. 3. The values of the parameters, common for both fits, are $`L=2.6\times 10^4`$ m, $`B=1.07\times 10^3`$ s3, $`C=3.0\times 10^{18}`$ s K-1, $`b=6.2`$, $`D=2.1\times 10^5`$, and $`E=7.5\times 10^5`$. We refrain from a detailed discussion of the parameters related to defect, boundary, and phonon-phonon scattering, and in the subsequent analysis we consider them as serving as a background. Our primary interest is the electron-phonon interaction represented by the parameter $`E`$. The value of $`E`$ obviously fixes the shift between the two broken lines in the inset of Fig. 3. We note that the value of $`E`$ is quite robust with respect to the choice of different phonon scattering terms in our fitting procedure.

In the essential step of the analysis, we aimed at establishing the temperature dependence of the gap function $`\mathrm{\Delta }(T)`$, using Eq. (Thermal and electrical transport in single crystalline MgB2). The immediate result of this analysis is that a single gap function with the temperature dependence given by the BCS theory Muehlschlegel59  is not adequate for MgB2. The significant reduction of the low-temperature thermal conductivity by a rather weak magnetic field of 0.62 kOe (see Fig. 2) implies a strong interaction between phonons and the quasiparticles in the cores of the vortices induced by the magnetic field. In case of such a strong electron-phonon interaction and the formation of a single BSC-type gap, a significant change of the slope of $`\kappa (T)`$ at $`T_c`$ is expected, as illustrated in Fig. 3, where the dotted line represents the calculation of $`\kappa (T)`$ for $`\mathrm{\Delta }(0)=1.76k_BT_c`$ and $`\mathrm{\Delta }(T)`$ as tabulated in Ref. Muehlschlegel59 . Hence, for the majority of quasiparticles involved in phonon scattering below $`T_c`$, the gap opens much more slowly with decreasing temperature than expected from the standard BCS prediction. Our fitting procedure is much more successful if we assume the formation of two energy gaps, as has previously been claimed in reports of other experimental studies of MgB2. Thus, we consider two subsystems of quasiparticles with gaps $`\mathrm{\Delta }_1`$ and $`\mathrm{\Delta }_2`$, different parameters $`E_1`$ and $`E_2`$ of phonon-electron scattering, and separate contributions $`\kappa _{e1}`$ and $`\kappa _{e2}`$ to the heat transport. The normal-state parameter $`E`$ is then given by the sum $`E=E_1+E_2`$. In the superconducting state, the last term in Eq. (Thermal and electrical transport in single crystalline MgB2) adopts the form

$$E\omega \left[(1\alpha )g(x,y_1)+\alpha g(x,y_2)\right],$$
(8)
where $`y_1=\mathrm{\Delta }_1(T)/k_BT`$ and $`y_2=\mathrm{\Delta }_2(T)/k_BT`$. The parameter $`\alpha =E_2/(E_1+E_2)`$ characterizes the relative weight of phonon scattering by quasiparticles in the two subsystems. Correspondingly,

$$\kappa _e=\frac{L_0T}{\rho _0}\left[(1\delta )f(y_1)+\delta f(y_2)\right],$$
(9)
where the parameter $`\delta =\kappa _{e2}/\kappa _e`$ and $`\kappa _e=\kappa _{e1}+\kappa _{e2}`$. The fit is shown in Fig. 3 as the solid line.

The resulting values of the fit parameters are $`\mathrm{\Delta }_1(0)=5.3\pm 0.5`$ meV, $`\mathrm{\Delta }_2(0)=1.65\pm 0.2`$ meV, $`\alpha =1\pm 0.02`$ and $`\delta =0.33\pm 0.02`$. The ratio $`2\mathrm{\Delta }_1(0)/(k_BT_c)=3.22`$ is close to the weak-coupling BCS value of 3.52. The second gap, however, is considerably smaller. Our identification of two gaps is in qualitative and quantitative agreement with conclusions from results of specific heat, point-contact spectroscopy, and photoemission experiments (see, e.g., Table I in Ref. Bouquet01cm ). The two gaps in MgB2 are, by most authors, associated with different sheets of the Fermi-surface. According to density-functional calculations of Liu et al. Liu01 , the smaller and the larger gap are associated with 3-dimensional sheets and 2-dimensional tubes, respectively. However, the first-principle FLAPW band calculations of Hase and Yamaji Hase01  claim the opposite. Our result of $`\alpha =1`$ indicates that the dominant part of the interaction between quasiparticles and low-frequency phonons is provided by that part of the electronic excitation spectrum experiencing the smaller gap. The parameter $`E`$ in Eq. (Thermal and electrical transport in single crystalline MgB2) can be expressed as

$$E=\frac{m^2E_{\mathrm{def}}^2}{2\pi v\mathrm{}^3D_0},$$
(10)
where $`m^{}`$ is an effective electron mass, $`E_{\mathrm{def}}`$ is the deformation potential and $`D_0`$ is the mass density. The effective masses of quasiparticles on different parts of the Fermi surface of MgB2 do not differ much Kong01cm . It has been shown An01  that the interaction between the holes of the cylindrical $`\sigma `$-bands and the in-plane bond-stretching optical mode $`E_{2g}`$ with a frequency of 470 cm-1 is provided by an ultrastrong deformation potential. If $`E_{\mathrm{def}}`$ is also strong for low-frequency acoustic phonons of the same polarisation, then our result implies that it is the holes of the 2D sheets of the Fermi surface that experience the smaller energy gap, a counterintuitive result but in agreement with the conclusions of Ref. Hase01 . Considering the contributions of the two subgroups of quasiparticles to the heat transport, our result $`\delta =0.33`$ suggests that the electrons experiencing the larger gap provide approximately 2/3 of total electronic heat conduction.

In conclusion, we present a study of the thermal conductivity $`\kappa `$ in the $`ab`$-plane of a single crystal of MgB2. The analysis of the temperature dependence of $`\kappa `$ provides new convincing evidence for the formation of two superconducting gaps which are associated with different sheets of the Fermi surface.

###### Acknowledgements.

Useful discussions with I. L. Landau, R. Monnier and F. Bouquet are acknowledged. This work was financially supported in part by the Schweizerische Nationalfonds zur Förderung der Wissenschaftlichen Forschung.


format: text

text:
# Simple solutions of fireball hydrodynamics for self-similar, ellipsoidal flows

## Abstract

A new family of simple, analytic solutions of self-similarly expanding fireballs is found for systems with ellipsoidal symmetry and a direction dependent, generalized Hubble flow. Gaussian, shell like or oscillating density profiles emerge for simple choices of an arbitrary scaling function. New, cylindrically or spherically symmetric as well as approximately one dimensional hydrodynamical solutions are obtained for various special choices of the initial conditions.

Introduction — Hydrodynamics is describing the local conservation of matter, momentum and energy. Due to this nature, hydrodynamical solutions are applied to a tremendous range of physical phenomena ranging from the stellar dynamics to the description of high energy collisions of heavy ions as well as collisions of elementary particles. Some of the most famous hydrodynamical solutions, like the Hubble flow of our Universe or the Bjorken flow in ultra-relativistic heavy ion physics have the properties of self-similarity and scale-invariance. Heavy ion collisions are known to create three dimensionally expanding systems. In case of non-central collisions, cylindrical symmetry is violated, but an ellipsoidal symmetry can be well assumed to characterize the final state. The data motivated, spherically or cylindrically symmetric hydrodynamical parameterizations and/or solutions of refs. jnr ; jde ; nr ; 3d ; qm95 ; nrt ; cssol ; cspeter ; ellsol ; csell  are generalized here to the case of such an ellipsoidal symmetry, providing new families of exact analytic hydrodynamical solutions.

The new family of self-similar ellipsoidal solutions — The non-relativistic (NR) hydrodynamical systems are specified by the continuity, Euler and energy equations:

$`_tn+(n𝐯)`$ $`=`$ $`0,`$ (1)
$`_t𝐯+(𝐯)𝐯`$ $`=`$ $`(p)/(mn),`$ (2)
$`_tϵ+(ϵ𝐯)`$ $`=`$ $`p𝐯.`$ (3)
Here $`n`$ denotes the particle number density, $`𝐯`$ stands for the NR flow velocity field, $`ϵ`$ for the energy density, $`p`$ for the pressure and in the following the temperature field is denoted by $`T`$. These fields depend on the time $`t`$ as well as on the coordinates $`𝐫=(r_x,r_y,r_z)`$. We assume, for the sake of simplicity, the following equations of state,

$`p`$ $`=`$ $`nT,ϵ=\kappa p,`$ (4)
which close the set of equations for $`n`$, $`𝐯`$ and $`T`$. The NR ideal gas corresponds to $`\kappa =\frac{3}{2}`$. The new family of exact analytic solutions of the hydrodynamical problem are given for arbitrary values of $`\kappa >0`$, $`m>0`$. The hydro solutions are determined by the choice of a positive function of a non-negative, real variable $`𝒯(s)`$, corresponding to the (dimensionless) scaling function of the temperature. The scaling variable $`s`$ is defined as

$`s`$ $`=`$ $`{\displaystyle \frac{r_x^2}{X^2}}+{\displaystyle \frac{r_y^2}{Y^2}}+{\displaystyle \frac{r_z^2}{Z^2}}.`$ (5)
There the scale parameters depend on time, $`(X,Y,Z)=(X(t),Y(t),Z(t))`$. The ellipsoidal symmetry of the solutions is reflected by the ellipsoidal family of surfaces given by the $`s=s_0=const`$ equation. The temperature and the density field depend on the coordinates $`(r_x,r_y,r_z)`$ only through the scaling variable $`s`$.

The new family of elliptically symmetric solutions of fireball hydrodynamics is given by

$`n(t,𝐫)`$ $`=`$ $`n_0{\displaystyle \frac{V_0}{V}}\nu (s),`$ (6)
$`𝐯(t,𝐫)`$ $`=`$ $`({\displaystyle \frac{\dot{X}}{X}}r_x,{\displaystyle \frac{\dot{Y}}{Y}}r_y,{\displaystyle \frac{\dot{Z}}{Z}}r_z),`$ (7)
$`T(t,𝐫)`$ $`=`$ $`T_0\left({\displaystyle \frac{V_0}{V}}\right)^{1/\kappa }𝒯(s),`$ (8)
$`\nu (s)`$ $`=`$ $`{\displaystyle \frac{1}{𝒯(s)}}\mathrm{exp}\left({\displaystyle \frac{T_i}{2T_0}}{\displaystyle _0^s}{\displaystyle \frac{du}{𝒯(u)}}\right),`$ (9)
where the constant $`n_0`$ is given by $`n_0=n(t_0,\mathrm{𝟎})`$, the typical volume of the expanding system is $`V=XYZ`$, the initial volume being $`V_0=V(t_0)`$, the dimensionless scaling function of the density profile is denoted by $`\nu (s)`$, the constant $`T_0`$ is defined by $`T_0=T(t_0,\mathrm{𝟎})`$ and a constant of integration is denoted by $`T_i`$. The definitions of $`n_0`$ and $`T_0`$ correspond to the normalization $`\nu (s=0)=1`$ and $`𝒯(s=0)=1`$. Initially, only one of the temperature and density profiles can be chosen as an arbitrary positive function, the equations of state relates the density and temperature profiles, resulting in the matching condition for the profile functions, as expressed by eq. (9).

The equations of motion of the scale parameters are

$$X\ddot{X}=Y\ddot{Y}=Z\ddot{Z}=\frac{T_i}{m}\left(\frac{V_0}{V}\right)^{1/\kappa }.$$
(10)
This time evolution of the radius parameters $`X`$, $`Y`$ and $`Z`$ is equivalent to the classical motion of a particle in a non-central potential, governed by the Hamiltonian

$$H=\frac{1}{2m}\left(P_x^2+P_y^2+P_z^2\right)+\kappa T_i\left(\frac{X_0Y_0Z_0}{XYZ}\right)^{1/\kappa },$$
(11)
where the canonical coordinates are $`(X,Y,Z)`$ and the canonical momenta are $`(P_x,P_y,P_z)=m(\dot{X},\dot{Y},\dot{Z})`$. This generalizes earlier results ellsol ; csell  from $`\kappa =3/2`$ and $`𝒯(s)1`$ to and to arbitrary $`𝒯(s)>0`$ and $`\kappa >0`$.

The conservation of energy by the classical Hamiltonian motion determines the physical meaning of the constant of integration $`T_i`$ as the initial potential energy corresponding to the initial internal energy of the fireball:

$`E_{\mathrm{tot}}`$ $`=`$ $`{\displaystyle \frac{m}{2}}(\dot{X}_0^2+\dot{Y}_0^2+\dot{Z}_0^2)+\kappa T_i,`$ (12)
where the total energy is denoted by $`E_{\mathrm{tot}}`$, and the initial velocities are denoted by $`(\dot{X}_0,\dot{Y}_0,\dot{Z}_0)`$. Due to the repulsive nature of the potential, the coordinates diverge for large values of $`t`$. The asymptotic velocities tend to constants ellsol  of $`(\dot{X}_{\mathrm{as}},\dot{Y}_{\mathrm{as}},\dot{Z}_{\mathrm{as}})`$,

$`E_{\mathrm{tot}}`$ $`=`$ $`{\displaystyle \frac{m}{2}}(\dot{X}_{\mathrm{as}}^2+\dot{Y}_{\mathrm{as}}^2+\dot{Z}_{\mathrm{as}}^2).`$ (13)
This completes the specification of the new family of solutions of fireball hydrodynamics with ellipsoidal symmetry. The form of the dimensionless scaling function $`𝒯(s)`$ can be chosen freely from among the positive functions of a non-negative variable. This freedom corresponds to a freedom in the specification of the initial conditions. Thus (uncountably) infinite new solutions of NR hydrodynamics are found.

Self-similarity of the elliptic hydro of solutions — Each of these new hydrodynamical solutions is scale invariant:

$`𝐫^{}`$ $`=`$ $`(r_x{\displaystyle \frac{X_0}{X}},r_y{\displaystyle \frac{Y_0}{Y}},r_z{\displaystyle \frac{Z_0}{Z}}),`$ (14)
$`n(t,𝐫)`$ $`=`$ $`n(t_0,𝐫^{})\left({\displaystyle \frac{X_0Y_0Z_0}{XYZ}}\right),`$ (15)
$`v_x(t,𝐫)`$ $`=`$ $`v_x(t_0,𝐫^{}){\displaystyle \frac{\dot{X}}{\dot{X}_0}},\mathrm{}`$ (16)
$`T(t,𝐫)`$ $`=`$ $`T(t_0,𝐫^{})\left({\displaystyle \frac{X_0Y_0Z_0}{XYZ}}\right)^{1/\kappa }.`$ (17)
Scale invariance of these solutions is equivalent to their self-similarity. The profile functions depend on time only through the scale parameters $`(X,Y,Z)`$ and on the coordinates only through the scale parameter $`s`$.

Limiting cases — The various physically interesting limiting cases fall into two classes. The first class of limiting cases corresponds to various additional symmetry properties imposed on the scaling variable of eq. (5). The second class of limiting cases corresponds to various choices of $`𝒯(s)`$, the scaling function of the temperature. The Lagrangian equations of motion for the scale parameters, eqs. (10) as well as the shape of the flow velocity field, eq. (7), are the same for all the choices of $`𝒯(s)`$. Trivial prefactors are given in eqs. (6,8). Hence these equations will not be repeated in the forthcoming discussion. We provide some physically interesting examples for the matching pair of density and temperature scaling functions $`(\nu ,𝒯)`$ that satisfy eq. (9).

Spherically symmetric family of solutions — By assuming that initially all the scale parameters as well as their time derivatives are equal, $`X_0=Y_0=Z_0R_0`$ and $`\dot{X}_0=\dot{Y}_0=\dot{Z}_0\dot{R}_0`$, the ellipsoidal family of solutions reduces to the spherical family of solutions of ref. cssol  with a scale parameter $`X=Y=ZR`$. The scaling variable and the equations of motion simplify to

$$s=\frac{𝐫^\mathrm{𝟐}}{R^2},R\ddot{R}=\frac{T_i}{m}\left(\frac{R_0^3}{R^3}\right)^{1/\kappa }.$$
(18)
This family generalizes the Zimányi-Bondorf-Garpman (ZGB) solution jnr , the spherical Gaussian solution nr , and the Buda-Lund type of hydro solutions cssol  to arbitrary scaling functions $`𝒯(s)>0`$ and equation of state parameters $`\kappa >0`$. If $`\dot{R}_0=0`$, and the asymptotic velocity of the expansion is $`\dot{R}_{\mathrm{as}}^2=u^2=T_i/m`$, and if $`\kappa =3/2`$, the equation of motion of the scale parameter is solved in a simple form cssol  as $`R^2=R_0^2+u^2(tt_0)^2`$.

Cylindrically symmetric family of solutions — By imposing cylindrically symmetric initial conditions, $`X_0=Y_0=R_{t0}`$, and $`\dot{X}_0=\dot{Y}_0=\dot{R}_{t0}`$, one finds that $`X=Y=R_t`$, as the equations of motion preserve cylindrical symmetry. Introducing $`r_t=\sqrt{r_x^2+r_y^2}`$, the scale parameter $`s`$ and the equations of motion for the longitudinal and the transverse scales read as

$`s`$ $`=`$ $`{\displaystyle \frac{r_t^2}{R_t^2}}+{\displaystyle \frac{r_z^2}{Z^2}},R_t\ddot{R}_t=Z\ddot{Z}={\displaystyle \frac{T_i}{m}}\left({\displaystyle \frac{R_{t0}^2Z_0}{R_t^2Z_0}}\right)^{1/\kappa }.`$ (19)
These generalize the equations of motion for scales of the cylindrically symmetric, De - Garpman - Sperber - Bondorf \- Zimányi (DGSBZ) solution of ref. jde  to arbitrary $`0<\kappa 3/2`$ and to arbitrary $`𝒯(s)>0`$.

One dimensional expansions — The equations of motion of parameters $`(X,Y,Z)`$ has been studied for the case of $`\kappa =3/2`$ in ref. ellsol . Although the expansion is generally 3 dimensional, a big initial compression in one of the directions ($`r_z`$) was shown to result in an effectively 1 dimensional expansion in this direction, corresponding to Landau type initial conditions. In this case, an analytic solution for the variable $`Z`$ is given in eqs. (23-25) of ref. ellsol , and the conditions of validity of this approximation were given there by eqs. (26-28). These equations can be further simplified during the late stage of the expansion, when acceleration effects are small. In this limiting case, both eqs. (10) and the conservation of energy are satisfied by the asymptotic solution:

$`\dot{X}_a`$ $``$ $`\dot{X}_0,\dot{Y}_a\dot{Y}_0,{\displaystyle \frac{1}{2}}m\dot{Z}_a^2{\displaystyle \frac{3}{2}}T_i,`$ (20)
$`X(t)`$ $``$ $`X_0+\dot{X}_at,`$ (21)
$`Y(t)`$ $``$ $`Y_0+\dot{Y}_at,`$ (22)
$`Z(t)`$ $``$ $`Z_0+\dot{Z}_at.`$ (23)
For simplicity, here we utilized $`t_0=0`$ . This asymptotic approximate solution is valid if the conditions of validity of the 1 dimensional expansion given in ref. ellsol  are satisfied simultaneously with the following constraints: $`Z_0\dot{Z}_at`$, $`X_0\dot{X}_0t`$ and $`Y_0\dot{Y}_0t`$. Alternatively $`\sqrt{m/(3T_i)}Z_0t\mathrm{min}(X_0/\dot{X}_0,Y_0/\dot{Y}_0)`$. Under these conditions all the initial internal energy is converted into kinetic energy in the longitudinal direction, while the kinetic energy in the transverse components is conserved during the time evolution. For $`\dot{X}_0=\dot{Y}_0=0`$, a one dimensional expansion is obtained.

Various choices for the shape of the temperature scaling function $`𝒯(s)`$ generate interesting forms of the hydrodynamical solutions in all the ellipsoidal, cylindrical, spherical or 1 dimensionally expanding classes.

Gaussian solutions — The simplest possible choice for the temperature scaling function is $`𝒯(s)=1`$. After a trivial scale transformation, $`(X,Y,Z)\sqrt{\frac{T_i}{T_0}}(X,Y,Z)`$,

$`\nu (s)`$ $`=`$ $`\mathrm{exp}\left(s/2\right),𝒯(s)=1.`$ (24)
The density profiles are Gaussians and the temperature distribution becomes spatially homogeneous, as follows from eq. (5), and we recover the elliptic Gaussian solutions described in refs. ellsol ; csell . If freeze-out happens at a constant value of the local temperature, $`T(t,𝐫)=T_f`$, these Gaussian hydrodynamical solutions corresponds to a sudden freeze-out at a constant time, $`t=t_f`$, in the whole volume of the fireball. Remarkable features of this model are that i) the slope parameters of the transverse momentum spectra increase linearly with mass, with the coefficient of linearity depending on the relative direction to that of the impact parameter, ii) the parameters of the two-particle Bose-Einstein correlation functions oscillate as a function of the angle between the event plane and the transverse momentum of the pair csell .

De - Garpman - Sperber - Bondorf - Zimányi solution — One can require that the temperature and the density profiles are described as different powers of the same profile function. Such a similarity is achieved if the temperature and the corresponding density profiles are

$`𝒯(s)`$ $`=`$ $`\left(1s\right)\mathrm{\Theta }(1s),`$ (25)
$`\nu (s)`$ $`=`$ $`\left(1s\right)^\alpha \mathrm{\Theta }(1s),\alpha =T_i/(2T_0)1.`$ (26)
These profile functions generalize the spherical ZGB solution jnr ; cssol , and the cylindrically symmetric DGSBZ solution jde  to asymmetric ellipsoids and arbitrary values of $`\kappa `$. The physical meaning of the parameter $`\alpha `$ is determined here, in terms of the total initial internal energy $`T_i`$ and the initial central value $`T_0`$ of the temperature field.

Elliptic Buda - Lund solutions — The Buda - Lund hydro model (BL-H) was developed for the description of single-particle spectra and two-particle Bose-Einstein correlation functions in high energy heavy ion collisions at CERN SPS 3d ; qm95 . The BL-H attempted to characterize the temperature, density and flow fields by their means and variances only, its essential property is that it simultaneously involves a temperature gradient parameter and a Hubble-like flow profile. In its original form, the model was cylindrically symmetric and the (longitudinal) flow profile was relativistic 3d ; qm95 . BL-H type of exact hydro solutions correspond to the scaling functions

$`𝒯(s)`$ $`=`$ $`{\displaystyle \frac{1}{1+bs}},b={\displaystyle \frac{1}{2}}{\displaystyle \frac{\mathrm{\Delta }T}{T}},`$ (27)
$`\nu (s)`$ $`=`$ $`\left(1+bs\right)\mathrm{exp}\left[{\displaystyle \frac{T_i}{2T_0}}(s+bs^2/2)\right].`$ (28)
The dimensionless parameter $`b`$ is interpreted as a measure of the transverse temperature inhomogeneity 3d ; qm99 . Regardless of the symmetry classes, the BL-H density profile has an approximately Gaussian, conventional shape if $`b<T_i/(2T_0)`$. On the other hand, the density profile looks like an ellipsoidal ring of fire, with a density minimum at the center and a density pile-up on the surface, if $`b>T_i/T_0`$. In the spherically symmetric case, similar morphological classes of BL-H solutions were found in ref cssol . Introducing the notation $`T_i=mu^2`$, one finds that $`\frac{mu^2}{T_0}>\frac{\mathrm{\Delta }T}{T}`$ yields ellipsoidal, expanding fireball BL-H solutions, Fig. 1. A detailed analysis of correlations and spectra in Pb + Pb collisions at CERN SPS indicated qm99  such a behavior, corresponding to big and expanding fireballs. On the other hand, the analysis of correlations and spectra of $`h+p`$ reactions at CERN SPS indicated csrev ; krev ; na22  small transverse flow and big transverse temperature inhomogeneity, $`\frac{mu^2}{T_0}<\frac{\mathrm{\Delta }T}{T}`$. This case corresponds to the formation of a shell of fire. Such an example is shown in Fig. 2.

Fireballs with density waves — Here we show that a simple choice of the temperature scaling function can lead to periodically oscillating temperature and density waves in exact solutions of NR hydrodynamics. A pair of oscillating scaling functions is, for example

$`𝒯(s)`$ $`=`$ $`(1+\alpha \mathrm{cos}\beta s)^1,`$ (29)
$`\nu (s)`$ $`=`$ $`\left(1+\alpha \mathrm{cos}\beta s\right)\mathrm{exp}\left[{\displaystyle \frac{T_i}{2T_0}}(s{\displaystyle \frac{\alpha }{\beta }}\mathrm{sin}\beta s)\right].`$
where parameters $`(\alpha ,\beta )`$ correspond to the amplitude and the period of the oscillations, respectively, as shown in Fig. 3. As time evolves, the shapes of the ellipsoids change in coordinate space: the bigger the initial compression, the faster the expansion in that direction, corresponding to eqs. (10).

Connection with other exact solutions of fireball hydrodynamics – Before summarizing the results let us also relate the presented family of exact solutions of hydrodynamics to presently known, exact solutions of relativistic hydrodynamics. One of the most well known relativistic hydrodynamical solution is the 1+1 dimensional Landau-Khalatnikov solution, described in refs. Landau:1953gs , Khalatnikov  and Belenkij:1956cd . The initial condition here is a uniformly heated piece of matter initially at rest, the equation of state is that of an ultrarelativistic gas with three degrees of freedom, but expanding only 1 + 1 dimensions. Landau has shown that such expansions lead to an approximately Gaussian rapidity distribution. This Landau-Khalatnikov solution is presently the only known exact solution of relativistic hydrodynamics that describes an exploding fireball with relativistic acceleration. These solutions are, however, not self-similar, and not explicit, one dimensional solutions, hence they are not easily related to the exact, explicit, selfsimilar and accelerating non-relativistic solutions described in the present manuscript.

However, it is interesting to note, that in the late time limit, the acceleration of the scale parameters vanishes in the non-relativisitic self-similar solutions described here, as given by eqs. (21-23). Hence for very late time, the velocity field becomes spherically symmetric, $`𝐯𝐫/t`$, which corresponds, within the lightcone, to a well known solution of relativistic hydrodynamics: $`u^\mu (x)=x^\mu /\tau `$, which is in the 1+1 dimensional case the reknowned velocity distribution of the Hwa-Bjorken solution of relativistic hydrodynamics, given in refs. Hwa:1974gn ,Bjorken:1982qr , while in 1 + 3 dimensions $`𝐯=𝐫/t`$ or $`u^\mu (x)=x^\mu /\tau `$ is referred to as the Hubble flow, as it corresponds to the flow velocity field of galaxies in an expanding Friedmann universe. These solutions are boost-invariant and both in the 1+1 dimensional Hwa-Bjorken and in the 1+3 dimensional Hubble case, corresponding to a flat rapidity distribution. The initial boundary conditions are given only within a lightcone, on a $`\tau =\tau _0`$ boost-invariant hypersurface, the equation of state is characterized by a constant or piecewise constant speed of sound, and the solutions are explicit, accelerationless and self-similar flows.

In ref. Csorgo:2003rt  these boost-invariant self-similar solutions were generalized to 1+1 dimensional, non-boost invariant, self-similar, accelerationless solutions, using a broad class of equations of state: $`ϵ=mn+\kappa p`$, $`p=nT`$. The initial conditions were given, similarty to the case of the Hwa-Bjorken solution, on a boost-invariant hypersurface with $`\tau =\tau _0`$, but assuming an inhomogeneous initial temperature profile and a corresponding matching initial density profile. In the same paper this solution is extended to axially symmetic, three dimensionally expanding fireballs, that represent the late stages of central heavy ion collisions. However, for non-central collisions, axial symmetry is too restricitve and due to this reason we have generalized these solutions for ellipsoidally expanding relativistic fireballs, and the 1+1 dimensional solution, the 1+3 dimensional axially symmetric solution and the 1+3 dimensional ellipsoidal solutions were written up also in a series of papers for the 3rd Budapest Winter School on Heavy Ion Physics, published in refs. Csorgo:2002ki ,Csorgo:2002bi ,Csorgo:2003ry . These self-similar, relativistic hydrodynamical solutions, in particular, the ellipsoidal expansion described in ref. Csorgo:2003ry , correspond exactly to the late time limit of the non-relativistic solutions described here. See ref. Csorgo:2003ry  for further details on such a correspondence as well as for a more complete list of citations on early papers on the Hwa-Bjorken solution.

It is worthwhile to mention that there are additional, recently found exact solutions of relativistic hydrodynamics that find accelerationless generalized Hubble type of solutions, with direction dependent Hubble constants. The first of these class of solutions has been found by Bíró in refs. Biro:1999eh ,Biro:2000nj , for cylindrically symmetric expansions at the softest point of the equation of state. Sinyukov and Karpenko recently published a solution which can be considered as the generalization of refs. Biro:1999eh ,Biro:2000nj , as well as that of ref. Csorgo:2003ry , for the pre-asymptotic, accelerationless stage of the expansion, when the scales expand already linearly in time, corresponding to eqs. (21-23), but before the time period when the off-sets $`(X_a,Y_a,Z_a)`$ can be considered negligibly small. The Bíró as well as the Sinyukov-Karpenko solutions are also self-similar and accelerationless relativistic solutions, their shortcoming is that they are valid in the medium only when the pressure is a temperature independent constant, or, when the boundary condition is an expansion to the vacuum, they are obtained only for the special ”dust” equation of state, that corresponds to a uniformly vanishing pressure, $`p=p_0=0`$. These solutions can thus be considered as a pre-asymptotic, relativistic solutions corresponding to very late stages of the fireball expansions.

This brief review of known exact solutions of fireball hydrodynamics in the relativistic kinematic domain indicates, that we find a correspondence between the late stages of the presented exact, accelerating, self-similar solutions of non-relativistic, ellipsoidally symmetric fireballs and the known class of relativistic but accelerationless, self-similar, ellipsoidally symmetric solutions of hydrodynamics. The missing link between them is the class of relativistic, accelerating, explicit solutions of hydrodynamics that use a realistic equation of state. The search for this missing class of exact and explicit solutions has been started and the first results will soon be reported elsewhere.

Summary — The non - relativistic hydrodynamical problem has been solved for expanding fireballs with ellipsoidal symmetry for the class of self-similar expansions. The flow velocity distribution is a generalized Hubble field in all the cases. An exact solution is assigned to each positive, integrable function of a non-negative variable. The time evolution of the $`(X,Y,Z)`$ scale parameters corresponds to a Hamiltonian motion of a mass point in a non-central, repulsive potential. The density profiles may be of fireball type, or they may form one or more shells of fire. For initial conditions with higher symmetry, one dimensional, cylindrical and spherical expansions are obtained.

These results provide analytic insight into the time evolution of expanding fireballs with nontrivial, ellipsoidally symmetric morphology, a non-polynomially hard algorithmic problem that is difficult to solve even numerically on von Neumann type computers.

Acknowledgments: Thanks are due to L. P. Csernai, Y. Hama, F. Grassi, M. Gyulassy, T. Kodama, B. Lukács and J. Zimányi for inspiring discussions. This work was supported by OTKA grants T026435, T034296, T 038406, NWO - OTKA grant N 25487, by an NSF - MTA - OTKA grant INT0089462, the US DOE grant DE - FG02 - 93ER40764 and by the NATO PST.CLG.980086 grant.


text:
# Untitled Document

Abstract: One of the fundamental aspects of statistical behaviour in many-body systems is exponential divergence of neighbouring orbits, which is often discussed in terms of Liapounov exponents. Here we study this topic for the classical gravitational $`N`$-body problem. The application we have in mind is to old stellar systems such as globular star clusters, where $`N10^6`$, and so we concentrate on spherical, centrally concentrated systems with total energy $`E<0.`$ Hitherto no connection has been made between the time scale for divergence (denoted here by $`t_e`$) and the time scale on which the energies of the particles evolve because of two-body encounters (i.e. the two-body relaxation time scale, $`t_r`$), even though both may be calculated by similar considerations.

In this paper we give a simplified model showing that divergence in phase space is initially roughly exponential, on a timescale proportional to the crossing time (defined as a mean time for a star to cross from one side of the system to another). In this phase $`t_et_r`$, if $`N`$ is not too small (i.e. $`N30`$). After several $`e`$-folding times, the model shows that the divergence slows down. Thereafter the divergence (measured by the energies of the bodies) varies with time as $`t^{1/2},`$ on a timescale nearly proportional to the familiar two-body relaxation timescale, i.e. $`t_et_r`$ in this phase. These conclusions are illustrated by numerical results.

Keywords: Gravity, Few-body systems, Relaxation processes, Particle orbits

1. Introduction

The classical gravitational $`N`$-body problem is defined by the equations

$$\ddot{𝐫}_i=G\underset{\begin{array}{c}j=1\\ ji\end{array}}{\overset{N}{}}m_j\frac{𝐫_i𝐫_j}{\left|𝐫_i𝐫_j\right|^3}$$
(1)
where $`𝐫_i`$ is the three-dimensional position vector of the $`i^{th}`$ star, $`m_i`$ is its mass, and $`G`$ is the universal constant of gravitation. We consider applications in which the total energy, $`E,`$ in the barycentric frame is negative and the total angular momentum is negligible. Starting from a rather broad set of initial conditions, such solutions settle down into a roughly spherical distribution of bodies in approximate “dynamic equilibrium” (Fig. 1), i.e. the spatial distribution is nearly time-independent on the time scale of the orbital motions of the particles.

Early numerical integrations with $`N32`$ showed that a small change in initial conditions led to a roughly exponential divergence of solutions (measured in $`6N`$-dimensional phase space), even though the spatial distribution of the bodies in the two solutions might be indistinguishable within statistical fluctuations. The timescale of divergence, $`t_e`$, was of order the crossing time, $`t_{cr},`$ defined in a certain conventional way as the time for a body with a typical speed to move a distance comparable to the size, $`R`$, of the spatial distribution of the particles. Thus

$$t_{cr}\frac{R}{V},$$
$`(2)`$
where $`V`$ is the root mean square speed of the particles. Later work extended numerical results to larger $`N`$, and Goodman et al devised theoretical models confirming that $`t_e/t_{cr}`$ is virtually independent of $`N.`$

One particular statistical specification of the initial conditions which has been studied is the Plummer model, which is often used in stellar dynamics for the study of relaxation and related processes. It is the stellar dynamical analogue of the $`n=5`$ polytrope. For this model it has been found that

$$t_e\frac{0.116t_{cr}}{\mathrm{ln}\left(0.73\mathrm{ln}N\right)}.$$

The functional form is suggested by a theoretical model, and the coefficients are not thought to depend sensitively on the initial conditions. Therefore for large star clusters generally we have

$$t_e0.05t_{cr}.$$
$`(3)`$
Fig.1. Spatial distribution of bodies in a typical simulation. On left is a snapshot, and on the right is the numerically generated space density as a function of radius.

The theoretical models of Goodman et al dealt with the linear divergence of neighbouring solutions, when the separation in position of the $`i^{th}`$ body satisfies the variational equation

$$\ddot{\delta 𝐫_i}=\underset{\begin{array}{c}j=1\\ ji\end{array}}{\overset{N}{}}m_j\left(\frac{\delta 𝐫_i\delta 𝐫_j}{\left|𝐫_i𝐫_j\right|^3}\frac{\left(\delta 𝐫_i\delta 𝐫_j\right).\left(𝐫_i𝐫_j\right)}{\left|𝐫_i𝐫_j\right|^5}\left(𝐫_i𝐫_j\right)\right).$$
$`(4)`$
For practical purposes, however (e.g. for understanding the growth of errors in a numerical integration) the resulting roughly exponential growth quickly leads to separations so large that the linear approximation fails. In this contribution we develop the simplest model of divergence to account for the later, nonlinear growth of the separation between neighbouring solutions. We shall see that the time dependence changes from roughly exponential to roughly power-law, and that the timescale changes from roughly the crossing time to nearly the two-body relaxation timescale, $`t_r.`$ This is the timescale on which the energies of the individual bodies vary significantly. Standard theory shows that

$$t_r\frac{N}{\mathrm{ln}N}t_{cr}$$
$`(5)`$
for systems of the general kind considered here.

2. A model of divergence

2.1 Linear growth of errors

In this section we introduce a toy model for the divergence of neighbouring orbits. Though it gives much insight into the physics of the problem, many details are omitted. In the first instance we apply it to the linear regime in which the approximate eqs.(4) are valid. In this regime more elaborate models have been constructed by Goodman et al.

We make the following assumptions. As in the theory of two-body relaxation we assume that the trajectory of a particle is nearly rectilinear, except for occasional two-body encounters (Fig. 2). We suppose that the important encounters are in the small-angle scattering regime, such that $`pGm/v^2`$ where $`p`$ is the impact parameter and $`v`$ is the relative velocity of the two particles. In computing the effect of one encounter, we suppose we can treat the scatterer as fixed. We also suppose that successive encounters can be treated as if motion takes place on one plane, and that the difference between two orbits is measured by the difference in the impact parameter, $`\delta p`$. We assume that all particles have the same mass $`m`$. Finally, we suppose that the system is in virial equilibrium (see Binney & Tremaine1), which implies that

$$V^2\frac{GmN}{R}.$$
$`(6)`$
Here the symbol $``$ means “is of order”, i.e. that the relation is approximate, and any numerical coefficient is ignored. Thus $`vV`$, for example.

Fig.2. Two successive encounters.

In the small-angle scattering regime the maximum acceleration of the moving particle is of order $`{\displaystyle \frac{Gm}{p^2}}`$ and the duration of the encounter is of order $`{\displaystyle \frac{p}{V}}`$. Thus the change in velocity is of order $`{\displaystyle \frac{Gm}{pV}}`$, and so the angular deflection is of order $`{\displaystyle \frac{Gm}{pV^2}}`$ (Fig.2). After the scattered body has travelled a further distance $`D`$ to its next encounter, its spatial deflection is of order $`{\displaystyle \frac{GmD}{pV^2}}`$.

Now suppose the body had approached the first encounter on a parallel path at a slightly different impact parameter $`p+\delta p.`$ Then, at the time of the second encounter, its position would have been displaced by a distance of order $`\delta p+{\displaystyle \frac{GmD}{p^2V^2}}\delta p.`$ The first term is the displacement that would have occurred even in the absence of the first encounter. The second occurs because, if $`\delta p>0,`$ the body has been deflected less by the first encounter. (The differential approximation used for this term is not valid unless $`|\delta p|p;`$ this is the approximation which restricts the present theory to the linear regime in which eqs.(4) are valid.) The total displacement measures the change in impact parameter at the second encounter. Hence the variation in $`p`$ is multiplied by a factor of order $`\left(1+{\displaystyle \frac{GmD}{p^2V^2}}\right)`$ per encounter111In a fully three-dimensional treatment this becomes a matrix equation..

Now we consider the cumulative effect of several encounters within a restricted range of impact parameters around the value $`p`$, e.g. from $`p/2`$ to $`2p`$, but ignoring other encounters. We start at some time $`t`$ and consider the effect of encounters in a subsequent interval $`\mathrm{\Delta }t`$, chosen sufficiently large that several encounters occur within this interval. The actual number of such encounters is of order $`{\displaystyle \frac{\mathrm{\Delta }tV}{D}},`$ and so the variation in the orbit is given222We ignore two complications which tend to counteract each other: (i) the persistence of effects of early encounters, and (ii) partial cancellation of successive encounters by their vectorial character. by

$$\delta r\left(t+\mathrm{\Delta }t\right)\delta r\left(t\right)\left(1+\frac{GmD}{p^2V^2}\right)^{{\displaystyle \frac{\mathrm{\Delta }tV}{D}}}.$$

Also it is clear that $`p^2Dn1`$, where $`n`$ is the number of particles per unit volume, and so

$$\delta r\left(t+\mathrm{\Delta }t\right)\delta r\left(t\right)\left(1+\frac{Gm}{p^4nV^2}\right)^{\mathrm{\Delta }tVnp^2}.$$
$`(7)`$
It follows from the relation $`n{\displaystyle \frac{N}{R^3}}`$ and eq.(6) that

$$\delta r\left(t+\mathrm{\Delta }t\right)\delta r\left(t\right)\left(1+\frac{R^4}{p^4N^2}\right)^{{\displaystyle \frac{\mathrm{\Delta }t}{t_{cr}}}{\displaystyle \frac{Np^2}{R^2}}},$$
$`(8)`$
where we have used eq.(2).

Encounters take place at a wide range of impact parameters $`p.`$ Writing eq. (8) as

$$\mathrm{ln}\delta r\left(t+\mathrm{\Delta }t\right)\mathrm{ln}\delta r\left(t\right)\frac{\mathrm{\Delta }t}{t_{cr}}N\frac{p^2}{R^2}\mathrm{ln}\left(1+\frac{R^4}{p^4N^2}\right)$$
$`(9)`$
we see that those with $`p<<RN^{1/2}`$ are individually very effective but too rare to dominate, whereas those with $`p>>RN^{1/2}`$ lose out by being individually ineffective, despite being very numerous. Encounters at impact parameter $`pRN^{1/2}`$ are most effective cumulatively, and lead to exponential growth of the deviation $`\delta r,`$ on a timescale of order $`t_{cr}`$.

Another way of seeing this is to sum the right hand side of eq.(9) over all impact parameters $`p<R`$. Since this term represents the effect of encounters with impact parameters in a range near some value $`p`$, the summation can be accomplished by multiplying by $`dp/p`$ and integrating. The result is that

$$\mathrm{ln}\delta r\left(t+\mathrm{\Delta }t\right)\mathrm{ln}\delta r\left(t\right)\frac{\mathrm{\Delta }t}{t_{cr}}\mathrm{ln}N.$$

Except for the logarithmic factor, this is equivalent to the result obtained by ignoring all encounters except those near $`pRN^{1/2}`$.

Many factors have been omitted from this simple model, including the distribution of velocities and density, and the curved orbits of bodies between encounters. Nevertheless, the results of more detailed models and numerical simulations, already quoted, confirm our basic result, except for a very weak $`N`$-dependence.

2.2 Nonlinear growth of separation

The above theory is valid as long as $`\delta r<<p,`$ and here we may take for $`p`$ the impact parameter for the most effective encounters, i.e. $`pRN^{1/2}.`$ Suppose we are interested in growth of errors in an $`N`$-body integration of eqs. (1), for a system which has been scaled so that $`R1.`$ Then we may have $`\delta r(0)10^{16}`$ for a double precision calculation, and so the linear approximation breaks down after about $`30t_e`$, i.e. between one and two $`t_{cr}`$, by eq.(3).

Thereafter we suppose that encounters with impact parameters $`p<<\delta r`$ are ineffective. Then we may estimate the growth of the separation of neighbouring orbits by substituting $`p\delta r(t)`$ in eq. (9), which gives

$$\mathrm{ln}\delta r(t+\mathrm{\Delta }t)\mathrm{ln}\delta r\left(t\right)\frac{\mathrm{\Delta }t}{t_{cr}}\frac{N\delta r(t)^2}{R^2}\mathrm{ln}\left(1+\frac{R^4}{N^2\delta r(t)^4}\right).$$

We are in a regime where $`\delta r(t)>RN^{1/2},`$ and so we can approximate

$$\mathrm{ln}\delta r(t+\mathrm{\Delta }t)\mathrm{ln}\delta r(t)\frac{\mathrm{\Delta }t}{t_{cr}}\frac{R^2}{N\delta r(t)^2}.$$
$`(10)`$
Since the term on the right depends on $`t`$, we can no longer conclude that $`\mathrm{ln}\delta r(t)`$ increases linearly with $`t.`$ To determine its time dependence we rewrite eq. (10) as a differential equation, i.e.

$$\frac{d}{dt}\mathrm{ln}\delta r(t)\frac{1}{t_{cr}}\frac{R^2}{N\delta r(t)^2}.$$

Ignoring for the moment the distinction between “$``$” and “$`=`$”, we obtain the solution

$$\delta r\left(t\right)=\left(\delta r\left(t_0\right)^2+2\frac{tt_0}{t_{cr}}\frac{R^2}{N}\right)^{1/2},$$

where $`t_0`$ is a constant, which may be interpreted as the time at which the growth of errors enters the nonlinear regime.

Well into the nonlinear regime we now see that $`\delta r\left(t\right)R\left({\displaystyle \frac{t}{Nt_{cr}}}\right)^{1/2}.`$ In order to interpret this result we shall estimate the difference in binding energy, $`\epsilon ,`$ of the body between the two neighbouring solutions. Now $`\epsilon {\displaystyle \frac{GNm}{R}},`$ and we can estimate $`\delta \epsilon {\displaystyle \frac{GNm\delta r}{R^2}}.`$ (We could obtain a similar estimate from consideration of the difference in velocity.) Hence $`{\displaystyle \frac{\delta \epsilon }{\epsilon }}\left({\displaystyle \frac{t}{Nt_{cr}}}\right)^{1/2}.`$ Now the two-body relaxation time, $`t_r`$, may be estimated by eq.(5), and so $`{\displaystyle \frac{\delta \epsilon }{\epsilon }}\left({\displaystyle \frac{t}{t_r}}\right)^{1/2}`$ if we ignore a logarithmic dependence on $`N.`$

3. Discussion

Recall that we are considering two solutions of eq. (1) starting with slightly different initial conditions. Suppose that we measure the separation of the two solutions by the separation in energy, $`\delta \epsilon ,`$ of a typical body. What we have concluded is that, for at most a few crossing times, $`\delta \epsilon (t)`$ grows exponentially, with an $`e`$-folding time comparable with $`t_{cr}`$ itself. Thereafter $`\delta \epsilon (t)`$ approaches a power law dependence, varying as $`t^{1/2},`$ on a timescale of the relaxation time.

The standard theory of relaxation tells us how $`\epsilon `$ (the energy of a given star) evolves on a single solution of the $`N`$-body equation. If we ignore variations of $`\epsilon `$ inside an encounter, $`\epsilon `$ performs a random walk on the timescale $`t_r,`$ and the change in $`\epsilon `$ varies as $`t^{1/2}.`$ (We here ignore the role of “dynamical friction”, which corresponds to the drift term in a Fokker-Planck description of the relaxation.)

Fig.3 illustrates these points using data from numerical $`N`$-body integrations with $`N=256`$. Two systems were integrated simultaneously using identical initial conditions except for a small difference in one coordinate of one particle. The solid curve (a) shows the mean square difference in the energies of the $`N`$ particles333Similar results have been presented by Merritt for motion in the gravitational field of $`N`$ fixed bodies. The corresponding initial conditions were also used for simultaneous integration of the variational equations, and the long-dashed curve (b) shows the corresponding mean square variation of energy. This grows nearly exponentially, but is followed by (a) for only a limited time of order a crossing time. The short-dashed curve (c) shows the mean square difference between the initial energy and the energy at time $`t`$, again averaged over the $`N`$ particles. This is caused by two-body relaxation. Evidently curve (a) departs from curve (b) around the point where the latter crosses curve (c), and then nearly follows (c). In this way we see that the growth of errors, which is exponential only in the linear regime, is consistent with the theory of two-body relaxation.

Fig.3. Mean square energy difference in numerical integrations with $`N=256`$, as a function of time. The meaning of the different curves is stated in the text. The results plotted are the mean of four independent runs. In the adopted units the crossing time is $`2\sqrt{2}`$.

The exponential divergence slows down to a power-law growth because close encounters become increasingly ineffective. There is a geometric way of looking at this. Krylov showed that the divergence could be understood as the behaviour of neighbouring geodesics on a certain manifold. As two neighbouring geodesics deviate further, their deviation is influenced less and less by the fine geometrical structure of the manifold across which they are proceeding.

Acknowledgement We thank the referees for their comments and suggestions.

References

1. Binney J., Tremaine S. (1987). Galactic Dynamics (Princeton University Press, Princeton)

2. Chandrasekhar S., (1942). Principles of Stellar Dynamics (University of Chicago Press, Chicago; also Dover Publications, New York )

3. Goodman J., Heggie D.C., Hut P. (1993). On the Exponential Instability of $`N`$-Body Systems, ApJ, 415, 715-33

4. Hemsendorf M., Merritt D. (2002). Instability of the Gravitational $`N`$-Body Problem in the Large-$`N`$ Limit, astro-ph/0205538

5. Kandrup H.E., Smith H. (1991). On the sensitivity of the $`N`$-body problem to small changes in initial conditions, ApJ, 374, 255-65

6. Krylov N.S. (1979). Works on the Foundations of Statistical Physics (Princeton University Press, Princeton)

7. Merritt D. (2001). Non-integrable Galactic Dynamics, in B.A. Steves, A.J. Maciejewski, eds, The Restless Universe, (SUSSP Publications, Edinburgh; and Institute of Physics Publishing, London)

8. Miller R.H. (1964). Irreversibility in Small Stellar Dynamical Systems, ApJ, 140, 250-6


text:
# Elliptic flow in intermediate energy heavy ion collisions and in-medium effects

## I Introduction

Heavy ion collisions is an area of research which seeks to study physical systems under extreme conditions of density and temperature in experiments occuring in terrestrial laboratories. As such, it represents a rich and challenging field well worthy of intellectual pursuit and its study is important for a deep and a more complete understanding of Nature. At high energies, one goal of this program is to form and study a new state of matter which is a tantalizing prediction of QCD: the quark-gluon plasma . At lower energies, the experimental and theoretical efforts have focused on the need to characterize and quantify the nuclear equation of state . This physics also has an important role to play in the theory of supernovæ and that of neutron star properties . It is that energy regime that we consider in this work.

In order to identify novel many-body features without ambiguity it is imperative to provide a realistic model of the nuclear reaction dynamics. An approach that has proven to be extremely successful is the Boltzmann-Uehling-Uhlenbeck (BUU) model of heavy ion collisions . In BUU simulations, nucleons can suffer hard collisions and can also move on curved trajectories owing to interactions with the self-consistent mean field. The interaction term we use in this work is introduced and described in Refs. . In connection with observables, a considerable amount of information on the nuclear equation of state and on nonequilibrium precursor phenomena is accessible through the measurement of global collective behaviour. A large body of work in intermediate energy nuclear collisions has been devoted to the measurement and to the theoretical calculation of both directed flow , and elliptic flow . For reviews, see , and references therein.

Near the low end of the intermediate energy spectrum, some studies have put forward the possibility of observing experimental signatures of new phenomena. A good example is that of modified in-medium nucleon-nucleon cross sections . Those would signal a departure from vacuum properties. While the confirmation of such manifestations would indeed be extremely interesting, one must keep in mind that such “new physics” issues must be addressed with an approach that incorporates all of the known physics in a computationally tractable model. Our goal in this paper is to first briefly describe such a model, and then to compare its results with experimental data. At the lower beam energies considered in this work, the problem of energy and momentum conservation in transport models is a pressing one. In addition, the momentum dependence of the nuclear mean field is an unavoidable feature, both from the point of view of theory and from that of experiment . With those issues in mind, we use a momentum-dependent lattice Hamiltonian model as a solution of the BUU transport equation. This specific model has previously been used in , and provides for excellent energy and momentum conservation. A similar approach was developed in parallel in . The momentum-independent lattice Hamiltonian method was first used in .

Our paper is organized as follows: the next section is a brief discussion of the momentum-dependent lattice Hamiltonian model. The following section introduces the self-consistent in-medium correction to the free-space nucleon-nucleon cross section as dictated by the functional momentum dependence that is observed in the nuclear mean field. Next, a comparison of model results with compiled and recent measurements of directed and elliptic flow in the range of laboratory bombarding energy $`E_k/A:`$25$``$800 MeV is performed. We then summarize and conclude.

## II The model

### A Lattice Hamiltonian solution of the mean field

In this work, we use the BUU transport model which characterizes the time evolution of a system of nucleons during the course of a collision of two heavy ions. The BUU equation reads:

$`{\displaystyle \frac{f(\stackrel{}{r},\stackrel{}{p},t)}{t}}+_\stackrel{}{p}h_\stackrel{}{r}f(\stackrel{}{r},\stackrel{}{p},t)_\stackrel{}{r}h_\stackrel{}{p}f(\stackrel{}{r},\stackrel{}{p},t)`$ (1)
$`={\displaystyle d^3p_1𝑑\mathrm{\Omega }^{}\left(v_{rel}^{}\times \frac{d\sigma ^{}}{d\mathrm{\Omega }^{}}\right)\left(f^{}f_1^{}\overline{f}\overline{f_1}ff_1\overline{f^{}}\overline{f_1^{}}\right)}.`$ (2)
The left-hand side of (1) accounts for the time evolution of the semi-classical nucleon phase space distribution function $`f`$ due to nucleon transport and mean field effects, and the right-hand side accounts for binary nucleon-nucleon collisions which also modify $`f`$. The nucleon single-particle Hamiltonian $`h=t+u`$, contains both the single particle kinetic energy and the mean field single particle potential. The starred quantities in the collision term represent in-medium values of the relative velocity and differential nucleon-nucleon scattering cross section respectively, and $`\overline{f}`$ accounts for Pauli blocking .

For the nucleon-single particle potential, we adopt the momentum-dependent parameterization for nuclear matter used in , supplemented with symmetry and Coulomb terms

$`u(\stackrel{}{r},\stackrel{}{p},\tau _3,t)=A\left({\displaystyle \frac{\rho (\stackrel{}{r},t)}{\rho _0}}\right)+B\left({\displaystyle \frac{\rho (\stackrel{}{r},t)}{\rho _0}}\right)^\sigma +{\displaystyle \frac{2C}{\rho _0}}{\displaystyle d^3p^{}\frac{f(\stackrel{}{r},\stackrel{}{p}^{},t)}{1+\left(\frac{\stackrel{}{p}\stackrel{}{p}^{}}{\mathrm{\Lambda }}\right)^2}}`$ (3)
$`\tau _3{\displaystyle \frac{2D}{\rho _0}}\left(\rho _n(\stackrel{}{r},t)\rho _p(\stackrel{}{r},t)\right)+{\displaystyle \frac{(\tau _3+1/2)}{4\pi ϵ_0}}{\displaystyle d^3r^{}\frac{\rho _p(\stackrel{}{r}^{},t)}{|\stackrel{}{r}\stackrel{}{r}^{}|}}.`$ (4)
This functional is amenable to computational applications and also incorporates nonequilibrium features . The strong interaction part of the mean field is fixed by specifying the value of the parameters $`A,B,\sigma ,C`$ and $`\mathrm{\Lambda }`$, neglecting the isospin and Coulomb terms. The parameter set from , deduced for symmetric nuclear matter, specifies: $`A=322.0`$ MeV, $`B=352.5`$ MeV, $`\sigma =12/11`$, $`C=62.75`$ MeV and $`\mathrm{\Lambda }=1.58p_f^0`$ which gives $`u(p=0)=72.4`$ MeV, $`u(p=p_f)=51.4`$ MeV and $`u(p\mathrm{})=+30.5`$ MeV for a zero temperature equilibrium nuclear matter distribution at saturation density $`\rho _0=0.16`$ fm-3. The zero temperature effective mass at the Fermi surface and nuclear matter compressibility for this parameterization are $`m^{}/m=0.67`$ and $`K=210`$ MeV respectively. The charge-dependent part of the nucleon mean field is taken into account with the isospin and Coulomb contributions to the single particle potential, where $`\tau _3`$ is the third component of isospin of the nucleon and $`\rho _n`$ and $`\rho _p`$ are the local neutron and proton densities, respectively. For the isospin part of the mean field, empirically it is found that $`D=34\pm 4`$ MeV. Furthermore, support for $`D`$ ranging from $`2740`$ MeV is found from various phenomenological investigations \[28, and Refs. therein\]. We adopt the value of $`D=32`$ MeV previously used in the context of simulations of heavy ion collisions . The above parameterization of the mean field provides a good fit to the nuclear equation of state for both symmetric nuclear matter and pure neutron matter as calculated in  for a only slightly larger value of the nuclear compressibility. Note however that the exact high density behavior of the symmetry energy is still an open question that stems partly from the uncertainty in the high density many-body calculations , and partly from the different theoretical paradigms leading to the equation of state .

The nucleon dynamics are realized through the use of the lattice Hamiltonian  equations of motion, whereby the semi-classical nucleon phase space density for $`N_{ens}`$ systems of $`A`$ nucleons is projected onto a configuration space lattice with a finite grid constant $`\delta x`$. At lattice site $`\alpha `$, the semi-classical phase space density reads:

$$f_\alpha (\stackrel{}{p},t)=\frac{1}{N_{ens}}\underset{i}{\overset{A\times N_{ens}}{}}R(\stackrel{}{r}_\alpha \stackrel{}{r}_i(t))\delta (\stackrel{}{p}\stackrel{}{p}_i(t))$$
(5)
With this, the lattice Hamiltonian equations of motion for particle $`i`$ read:

$$\frac{\stackrel{}{r}_i}{t}=\frac{\stackrel{}{p}_i(t)}{(p_i(t)^2+m^2)^{1/2}}+N_{ens}(\delta x)^3\underset{\alpha }{}R(\stackrel{}{r}_\alpha \stackrel{}{r}_i(t))_{\stackrel{}{p}_i}u_\alpha (\stackrel{}{p}_i,\tau _3^i,t)$$
(6)
$$\frac{\stackrel{}{p}_i}{t}=N_{ens}(\delta x)^3\underset{\alpha }{}u_\alpha (\stackrel{}{p}_i,\tau _3^i,t)_{\stackrel{}{r}_i}R(\stackrel{}{r}_\alpha \stackrel{}{r}_i(t)),$$
(7)
where $`u_\alpha (\stackrel{}{p}_i,\tau _3^i,t)`$ is the discretized form of (3), and $`R`$ is the finite width configuration space nucleon form factor. In this work, we use $`N_{ens}=100`$, $`\delta x=1.00`$ fm and a quadratic form factor which provides for good energy and linear momentum conservation as shown in . In one dimension, the form factor is defined as $`R(x)=(3\delta x^2/4x^2)/\delta x^3`$ for $`0x\delta x/2`$, $`R(x)=(3\delta x/2x)^2/2\delta x^3`$ for $`\delta x/2<x3\delta x/2`$ and $`R(x)=0`$ for $`3\delta x/2<x`$. The normalization factor of $`1/N_{ens}`$ in (5) and $`1/\delta x^3`$ and $`1/2\delta x^3`$ in $`R`$ are chosen such that particle number is conserved. That is:

$$A=(\delta x)^3\underset{\alpha }{}d^3pf_\alpha (\stackrel{}{p},t),$$
(8)
at all times $`t`$. Note that the lattice Hamiltonian equations of motion (6) and (7) are derived from the total Hamiltonian $`H`$, and not the single particle Hamiltonian $`h`$, as is the case with the often-used test-particle method . This work constitutes the first application (with Refs.  and ) of the lattice Hamiltonian approach with the interaction of Eq. (3) calculated fully dynamically. See also Ref. .

### B In-medium effects in nucleon collisions

For the calculation of the Uehling-Uhlenbeck collision term, we follow the cascade procedure outlined in  in conjunction with a Pauli blocking algorithm which distinguishes between protons and neutrons. However, the in-medium scattering cross section is calculated such that the functional dependence on momentum that is observed in the nuclear mean field is respected . Fermi’s Golden Rule gives an expression for the total free space nucleon-nucleon cross section in vacuum:

$$\sigma _{free}=\frac{2\pi V}{\mathrm{}}\frac{D_f^{(0)}}{v_{rel}^{(0)}}|t_{fi}|^2,$$
(9)
where $`v_{rel}^{(0)}`$ and $`D_f^{(0)}`$ are the relative velocity of the colliding nucleons and density of final states for the two-nucleon scattering process respectively. Both of these quantities are expressed as the free space value. In matter, we write the total in-medium nucleon-nucleon cross as

$$\sigma ^{}=\frac{2\pi V}{\mathrm{}}\frac{D_f^{}}{v_{rel}^{}}|t_{fi}^{}|^2,$$
(10)
where the starred quantities represent in-medium values. Detailed many-body calculations offer support for $`t_{fi}^{}t_{fi}`$, as expected for low local densities and/or low energy collisions \[38, and Refs. therein\]. We follow this line of thought here. In , the in-medium cross-section is calculated as it is here, however a simplified momentum-dependence in the mean field potential was used. We rewrite the in-medium elastic scattering cross section as:

$$\sigma ^{}=\frac{v_{rel}^{(0)}}{v_{rel}^{}}\frac{D_f^{}}{D_f^{(0)}}\sigma _{free}.$$
(11)
We note that a closed form solution to (11) exists for an equilibrium nuclear matter distribution. In this case, the in-medium nucleon-nucleon total cross section reads:

$$\sigma ^{}=\left(\frac{m^{}(\rho ,p)}{m}\right)^2\times \sigma _{free},$$
(12)
where the nucleon effective mass $`m^{}`$ is defined as:

$$\frac{m^{}(\rho ,p)}{m}=\left(1+\frac{m}{p}\frac{u(\rho ,p)}{p}\right)^1,$$
(13)
and we have replaced the $`\stackrel{}{r}`$ dependence in $`u`$ with a $`\rho `$ dependence which is appropriate for an equilibrium nuclear matter distribution. Note that the density $`\rho (\stackrel{}{r})`$ is a constant in this case and it enters into the above relation self-consistently via the semi-classical phase space distribution function $`f`$ (see, for example, Ref. ). In this work a modern and accurate parameterization of the nucleon-nucleon cross sections is used . We show in figure 1 the ratio of the nucleon effective mass to the free mass as a function of momentum, density and temperature for an equilibrium nuclear matter distribution. As seen in the figure, the nucleon effective mass approaches the free mass in the high momentum, high temperature and low density limit as is expected from the many-body nature of this phenomenon.

In figure 2 we show the ratio of the in-medium to that of the free space nucleon-nucleon total cross section as a function of both density and momentum for an equilibrium nuclear matter distribution at $`T=0`$. Note that for momentum-independent mean field potentials, the effective mass is equal to the free mass and consequently, the in-medium cross section is equal to the free space cross section.

### C Dynamical effects in heavy ion collisions

We now wish to apply our model developed thus far to simulate the collision of two heavy ions. In this scenario, dynamical effects will generate a nuclear matter distribution that is in general not that of a zero temperature system. Furthermore, one also expects a nuclear matter distribution that is different from the equilibrium situation. In this case, one cannot use the zero temperature equilibrium closed form solution of (3) as given in . Instead, (3), (6), (7) and (11) must be explicitly and dynamically calculated at all times during the collision process. In this work, we present the first results of the self-consistent calculation of (11) as applied to the dynamical collision of heavy ions.

In the initial stages of the collision of two heavy ions, large deviations from equilibrium will result in binary nucleon-nucleon centre of mass collision energies higher than that at later stages in the collision process. From figure 2 one can infer that initially, owing to the density and momentum-dependence in (12), the in-medium cross section will be larger at earlier stages in the collision compared to later times. In addition, one can also infer that higher centre of mass nucleus-nucleus collision energies will result in an in-medium cross section that is closer to the free space value. We illustrate this for a simulation of the head-on collision of two 209Bi nuclei at laboratory bombarding energies of 25, 150, 500 and 1000 MeV/A in the left panel of figure 3. In this figure, we show the mean value of the in-medium nucleon-nucleon cross section for unblocked collisions as a function of time. We note that previous attempts to model the in-medium nucleon-nucleon cross section have used a constant overall scaling factor . Others have used a phenomenological density-dependent in-medium cross section of the form

$$\sigma ^{}=(1\alpha \rho /\rho _0)\sigma _{free},$$
(14)
where $`\alpha `$ is an adjustable parameter. We show in the right panel of figure 3 that adoption of this simple ansatz results in an in-medium cross section that has considerably different behaviour from that of the self-consistent calculation. The self-consistent modification to the free space cross section in (11) effectively probes the shape of the momentum-dependent part of the single particle potential at all momenta. Such a detailed modeling is not possible with neither a constant scaling factor nor the density-dependent relation in (14). We will address the consequences of these differences in terms of elliptic flow in section III. Before we leave this section, we show in figure 4 a sample (250 events at each energy) of the values of the self-consistent in-medium cross section as a function of density for the $`t:1020`$ fm/c time slice from figure 3 for all energies.

Owing to the additional momentum-dependence in the in-medium cross section, we see that for a given density, the former can take on drastically different values. The solid lines in that figure indicate the value of the in-medium cross section that would be obtained by using the relation given in (14) for two different values of $`\alpha `$. Again, we see that (14) clearly oversimplifies the in-medium result from (11).

To summarize, we note that the calculation of the self-consistent momentum (and density) dependent in-medium cross section indicates that the latter approaches the free space cross section for high centre of mass nucleus-nucleus collision energies and steadily decreases with centre of mass nucleus-nucleus collision energies. Also, at early stages in the nucleus-nucleus collision, non-equilibrium effects results in an in-medium cross section that is closer to the free space value then compared to later times in the nucleus-nucleus collision process. This implementation of a self-consistent, momentum-dependent, BUU transport model represents a considerable numerical challenge. Efforts to develop algorithms in a parallel architecture are under way.

## III Collective Effects

### A Directed Flow

One of the often discussed signatures of collective motion in heavy ion collisions is the directed flow. Typically, the projection into the reaction plane of the transverse momenta of the reaction products is plotted against the rapidity, resulting in the well known “S-shaped” plot. The directed flow is obtained by taking the slope of this plot at mid-rapidity. Specifically, one fits the flow plot with:

$$\frac{d<p_x>}{dy}=a_1+a_2y+a_3y^3,$$
(15)
where the rapidity fit interval is taken to be $`|y|<y_{proj}`$, and $`y_{proj}`$ is the initial projectile rapidity. In this case, the flow $`F`$ is given by the value of $`a_2`$ in (15). Often, a narrow mid-rapidity region is selected, and (15) is used to fit the flow, with $`a_3=0`$. Some recent flow data from the FOPI collaboration  is shown in figure 5.

In this figure, the directed flow for Au+Au collisions at 400 MeV/A is shown as a function of impact parameter. The experimental data in this plot was obtained by using a coalescence model which weights all detected fragments by its measured charge. We also show in this figure the results of our BUU simulations with a free space nucleon-nucleon cross section, and the in-medium cross section given in (11). Here, the normalized flow $`F^{(0)}=F\times (y_{proj}/p_{proj})_{CM}`$ is used. The coalescence-invariant prescription was developed to enable direct comparisons with transport models such as the BUU. All the BUU results were obtained by using (15) with and without $`a_3=0`$. For the $`|y_n|=|y/y_{proj}|<0.45`$ cut, we find the two methods to be within one standard deviation of each other. Qualitatively, we find that both the free space and in-medium simulations predict the observed maximum in the flow magnitude for intermediate impact parameters. More quantitatively, we find the free space cross section over-predicts the flow magnitude for all impact parameters. The calculations with the in-medium cross section however, are generally in good agreement with the data at all impact parameters.

The directed flow excitation function also provides a means by which comparisons of measured signals to model calculations may be performed. With this in mind, we show in figure 6 a compilation of data near the balance energy  for central impact parameters and charge $`z=2`$ fragments and for semi-central impact parameters at $`E_k/A:150800`$ MeV  for charge $`z=1,2`$ fragments. The BUU results for the central impact parameters are shown for $`|y_n|<0.5`$ and $`\widehat{b}=b/b_{max}<0.39`$ fm and for $`0.19<y_n<+0.29`$ and ($`3.5<b<6.5`$) fm for the semi-central impact parameters. To compare our BUU results with the measured fragment flow, free nucleons are omitted from the flow analysis; see Ref.  for a discussion.

For the entire energy range probed here, we find that the calculations employing the in-medium nucleon-nucleon cross section result in a reduced value of the directed flow as compared with the free space cross section. From this figure, we find the in-medium result to be in better agreement with the data. A more stringent test of the model in terms of the directed flow is shown in figure 7. In this figure, the data for free protons  is shown along with the model predictions for free protons. The figure indicates that the predicted directed flow over-estimates the observed value. We point out however, that the coalescence-invariant comparison in figure 5 with the more recent data is in better agreement with the model calculations.

### B Elliptic Flow

We now turn to comparisons of the measured elliptic flow with our BUU model. The elliptic flow is a measure that quantifies the azimuthal anisotropy of the momentum distribution. Specifically, we fit the azimuthal distribution of nucleons about the reaction plane with a Fourier expansion of the form:

$$\frac{dN}{d\varphi }=p_0\left(1+2p_1cos(\varphi )+2p_2cos(2\varphi )\right),$$
(16)
where the ellipticity coefficient $`p_2`$ depends on the in-plane and out-of-plane flow amplitudes. In addition, the anisotropy ratio is defined by: $`R_n=(12p_2)/(1+2p_2)`$. A ratio $`R_n`$ larger than unity signals a preferred out-of-plane emission.

In  it was shown that BUU calculations of the elliptic flow for a Bi+Bi system with a momentum-dependent nuclear mean field favoured a nucleon effective mass of $`m^{}/m0.60.7`$. We show in figure 8, the results that we obtain for BUU calculations with the mean field in (3) of effective mass $`m^{}/m=0.67`$, for both a free space cross section and separately, the in-medium cross sections obtained with equations (11) and (14). The figure indicates that adoption of the self-consistent correction to the free space cross section is significant for the anisotropy ratio at the lower energies only (200 and 400 MeV/A). For the higher bombarding energies, the self-consistent result does not differ significantly from the free space result. Since the in-medium effect from (11) probes the shape of the momentum-dependence in the mean field it is more effective for lower bombarding energies, with many nucleon-nucleon collisions occurring near the Fermi surface where the momentum-dependence is steepest. At high bombarding energies, the nucleon-nucleon collisions occur far from the Fermi surface where the momentum-dependence in the mean field is relatively flat, and as such, the correction to the free space cross section is small. The experimental results are from . Unfortunately there are no data at the lower energy, where the differences in theoretical results are the largest.

We now turn to the elliptic flow excitation function for Au+Au recently measured by the FOPI collaboration . In figure 9, we show the ellipticity coefficient measured for fragments with $`A=1,2,3`$ along with the appropriate BUU result. We find that both parameterizations of the nucleon-nucleon cross section (free space and in-medium result from (11)) are able to qualitatively reproduce the experimental trend. Quantitatively, the calculation with the in-medium cross-section is better able to reproduce the measured signal.

We note that the ellipticity parameter is calculated in a frame where the principle axis of the momentum ellipsoid coincides with the beam axis. A more stringent test of the model is shown in figure 10 where the ellipticity coefficient for free protons with transverse momenta $`p_{}^{(0)}>`$0.8 only, is plotted as a function of incident bombarding energy. Note that $`p_{}^{(0)}p_{}/p_{beam}`$. From this figure, we find that the elliptic flow in the energy range ($`150<E_k/A<600`$) MeV is over-predicted by the model. This agrees with figure 8, where we see that the magnitude of the elliptic flow is over-predicted at high transverse momentum. This identifies the phase space areas where the BUU seems to perform perhaps less satisfactorily. However there, the phase space distribution functions are small and therefore those regions do not contribute importantly to the momentum-integrated observables. This is seen in our results as well as in the experimental data. This observation however represents valuable information for future theoretical developments.

## IV Conclusion

In summary, we have shown for the first time the effect of the parameter-free self-consistent calculation of the nucleon-nucleon in-medium scattering cross section in BUU calculations of both directed and elliptic flow from $`E_k/A:`$25$``$800 MeV. Our results indicate that, globally, the in-medium cross section better describes the flow data than with the free space cross section. In addition, the self-consistent calculation of the in-medium cross section results in a value of the latter that is less than the free space cross section for low bombarding energies and approaches the free space cross section at high bombarding energies. While in some sense transport models utilizing vacuum cross sections are a well-defined low-density approach, there is still some work to be done to formulate and use a totally self-consistent many-body transport theory. For example, one needs to generalize to the coupled set of transport equations that include all additional degrees of freedom . The work in this direction is continuing.

###### Acknowledgements.

This work is supported in part by the Natural Sciences and Engineering Research Council of Canada and in part by the Fonds FCAR of the Québec government. One of us (D.P.) is happy to acknowledge the financial support of McGill University through the Alexander McFee Memorial Fellowship.