text:
# On coarse-graining by the Inverse Monte Carlo method: Dissipative Particle Dynamics simulations made to a precise tool in soft matter modeling

## I Introduction

Classical molecular dynamics computer simulations with site-site pair potentials can currently be run for molecular systems consisting of the order of $`10^4`$ atoms, corresponding typically to a system size of 50 – 60 Å. Corresponding time scales using the same approach are typically of the order of 10 ns. This may be enough for simulations of isotropic liquids of simple molecules, while for many complex systems the requirements for the length and time scales are much more extensive.

The restrictions of the detailed all-atom description are especially severe in soft matter systems such as biomembranes and polyelectrolytes, although these are just two examples of the diversity of soft materials. This is mainly due to a wide range of time and length scales associated with these systems, and it is obvious that these problems become particularly pronounced in studies of dynamic processes which usually take place under hydrodynamic conditions. The primary aim in coarse-graining is to build a bridge between different time and length scales. Although there is no unique way to do this, the coarse-graining procedure always leads to a reduction of the number of degrees of freedom and thus alleviates the computational burden. The challenge is to establish systematic coarse-graining schemes that allow one to develop simplified model systems in terms of the information extracted from the underlying microscopic systems. Thus far, various approaches have been used to coarse-grain atomic and molecular systems. Below, we will briefly discuss some of them.

As a relatively generic approach, the dissipative particle dynamics (DPD) hoogerbrugge92a ; koelman93a ; espanol95a  has been used to study coarse-grained particle and polymeric systems, see e.g. Ref. warren98a  for a review. In our approach we use DPD as a thermostating method which is discussed in more detail in the following sections. To mention other relatively generic approaches, Flekkøy et al. flekkoy99a ; flekkoy00a  have recently introduced a hybrid method to combine particle and continuum level descriptions. This framework allows a derivation of the DPD model and in practice uses a Voronoi tessellation based technique for simulations of coarse-grained particles at different levels. Another very promising approach was presented by Malevanets and Kapral malevanets99a ; malevanets00a  who proposed a framework to couple a molecular level description with a mesoscale treatment of solvent that conserves hydrodynamics. This approach seems particularly appealing for studying dilute systems with hydrodynamics.

A lot of impetus for multiscale modeling has come from polymer research due to both fundamental and practical reasons. Akkermans and Briels akkermans00a  applied a projector operator formalism based approach to coarse-grain polymer chains using a microscopic model as a starting point. Their approach is appealing and deserves further studies since it accounts for the fluctuating forces at the coarse-grained level. Bolhuis et al. louis00a ; bolhuis01a  used the Ornstein-Zernike equation Hansen86  with a hypernetted-chain closure to invert the pair correlation function in order to derive effective pair potentials for neutral polymers in a solution. The coarse-grained polymers were represented as soft colloidal particles. Their approach is based on a theorem proved by Henderson in 1974 hend74  that stipulates a unique, point-to-point correspondence between pair potentials and radial distribution functions $`g(r)`$. The importance of this theorem becomes clear when one notices that for any $`g(r)`$ under given conditions there is a unique pair potential that includes corrections from the many-body interactions, and that by definition the pair potential then yields the original $`g(r)`$. The problem is thus reduced to inverting the $`g(r)`$ obtained from a microscopic model – how to obtain a coarse-grained model from that information lies at the heart of the problem of coarse-graining. At the representational level the soft ellipsoid model for polymer melts by Murat and Kremer murat98a  is very close to the one of Bolhuis et al. louis00a ; bolhuis01a , but instead of using the pair correlation function Murat and Kremer base their approach on separating the free energy into intra- and interchain parts. In an attempt to model more realistic systems, Reith et al. reith00a  introduced a coarse-graining procedure that starts from a detailed chemical description of a polymer. In order to obtain a coarse-grained approach, a simplex algorithm based optimization is performed to obtain a coarse-grained force field. The structural features produced by this approach are in good agreement with experimentally observed ones. For a more complete overview of the current coarse-graining approaches, see e.g. Ref. mplathe00  and references therein.

Our approach is also based on inverting the radial distribution function but instead of the integral equation based approach of Bolhuis et al. louis00a ; bolhuis01a  we use the inverse Monte Carlo procedure of Lyubartsev and Laaksonen LL95  which will be discussed in the following section. This systematic coarse-graining method is then combined with a momentum conserving dissipative particle dynamics (DPD) thermostat. Our approach consists of three conceptually simple steps. First, we use a detailed atomistic description and perform MD simulations. From these simulations we compute radial distribution functions $`g(r)`$ between different atoms, molecules, or molecular groups. Second, we apply the Inverse Monte Carlo LL95  procedure to invert the radial distribution functions. This yields effective interaction potentials $`V^{\mathrm{eff}}(r)`$ between the selected interaction sites. Finally, we employ these effective, coarse-grained potentials in the DPD algorithm to study the large-scale properties of systems in question. Note that since DPD satisfies momentum conservation, the present approach allows studies of a given system with full hydrodynamics.

Preliminary results of this work were recently reported in Ref. karttunen02 , in which we presented only the most essential ideas of this approach. Our aim in this article is to provide one with a comprehensive discussion of issues related to coarse-graining and mesoscopic simulations in terms of the DPD method. In particular, we discuss in detail how atomistic MD simulations can be coupled to DPD by the Inverse Monte Carlo method. To this end, we use NaCl as our model system. Even in this simple system, and with a very modest degree of coarse-graining, we obtain a computational speed-up of one order of magnitude. Our procedure is well-defined and truly allows easy and controllable tuning of the level of coarse-graining. Finally, we show through coarse-grained DPD simulations that this approach is physically sound and provides results in excellent agreement with the MD simulations and experiments.

This paper is organized as follows. First, in Sect. II, we describe the methods: The Inverse Monte Carlo procedure and the DPD algorithm. In Sect. III we present the results and finally, in Sect. IV, we discuss the results and some general aspects related to the general applicability of the method.

## II Methods and models

As discussed in the Introduction, detailed atomistic simulation techniques are suitable for studies of microscopic features of soft matter systems, e.g. the excited states of a chromophore where the time scales of interest are within a few hundreds of femtoseconds. For studies of systems that are characterized by large time and length scales, such as a protein in water solution over tens or hundreds of nanoseconds, however, atomistic simulation techniques are not very feasible. Our aim in this section is to discuss ways how this problem can be resolved by coarse-graining microscopic systems. In particular, we concentrate on an approach based on solving the “inverse problem” which yields effective interaction potentials which in turn can be coupled to mesoscopic simulation techniques such as the dissipative particle dynamics method discussed below. Full details of the methods can be found in Refs. LL95  and Bes00 ; Vattulainen02 , respectively.

### II.1 The Inverse Problem in statistical mechanics

In the standard statistical mechanical description of soft matter systems one begins from a formulation of a model, which is usually given in terms of a Hamiltonian or interaction potentials between atoms and molecules. When a Hamiltonian is specified, one can use different numerical methods to compute canonical averages such as energies, distribution functions, and response functions which can be compared with experiments. As one has an explicit expression for the canonical averages, their calculation is in principle a fairly straightforward task, although in many cases it may be computationally a very demanding one.

Occasionally one is interested in the solution of the inverse problem, i.e., how to deduce information about molecular interactions from the canonical averages which may be known from experiments. More specifically, the question is how to reconstruct the Hamiltonian based on results for the canonical averages. Many experimental properties can be chosen as a starting point to this end. Among the most important ones for the liquid state are the radial distribution functions $`g(r)`$, which may be known from the structure factor measurements using X-ray or neutron scattering Hansen86 . In 1974, Henderson proved a theorem hend74  stipulating a unique, point-to-point correspondence between pair potentials and radial distribution functions. In essence it states that for a given system under given conditions of temperature and density, two pair potentials which give rise to the same radial distribution function cannot differ more than by an additive constant. As this constant can be defined from the condition that the interaction potential tends to zero at an infinite distance, the potential is uniquely defined. In practice the solution of an inverse problem is not a straightforward calculation, however, since $`g(r)`$ does not provide one with a formal expression for the potential. Therefore, some special techniques are required for the solution of the inverse problem.

In 1988, McGreevy and Pusztai suggested a Reverse Monte Carlo (RMC) method in which the starting point for a simulation was the radial distribution function RMC88 . In this method, Monte Carlo simulations are carried out without any prior knowledge of the interaction potential to fit $`g(r)`$ that serves as an input. The set of configurations obtained in RMC may be used for further structural analysis, for example for calculation of three-dimensional spatial or orientational correlations. However, as the interaction potential is not recovered, the inverse problem is not solved completely, and it is not possible to calculate thermodynamical or dynamical properties using this approach. For the same reason this approach in not suitable for the development of coarse-grained models.

In another approach, Reatto et al. reatto86  used the Hypernetted Chain (HNC) approximation to solve the inverse problem. This and similar approaches, based on some closure of the Ornstein-Zernike equation, have been used in a number of works during the last decade to compute interaction potentials from radial distribution functions rosenfeld97 ; gonzalez98 . It should be noted, however, that computations within the HNC theory are feasible only for relatively simple models. Moreover, although yielding sometimes very accurate results LyMar02 , the HNC theory is not an exact mathematical solution of the statistical-mechanical problem, and its accuracy should be investigated carefully in every specific case. Other works devoted to the inverse problem are listed in Refs. ostheimer89 ; soper96 ; toth00 .

A practical way to solve the inverse problem has been suggested by Lyubartsev and Laaksonen LL95 . This is the method we will concentrate on in the following, and we will use it to compute effective potentials from radial distribution functions obtained from detailed MD simulations. The effective potentials then allow us to to study the large-scale properties of a given model system through coarse-grained DPD simulations.

### II.2 The Inverse Monte Carlo Method – A Simple Case

To illustrate the Inverse Monte Carlo method, we consider the case of a single-component system consisting of identical particles with pairwise interactions. The general case of a multicomponent system is a straightforward extension of it. The Hamiltonian for this system is given as

$$H=\underset{i,j}{}V(r_{ij}),$$
(1)
where $`V(r_{ij})`$ is the pair potential and $`r_{ij}`$ is the distance between the interaction sites $`i`$ and $`j`$. Let us assume that we know the radial distribution function $`g(r_{ij})`$. Our aim is now to find the corresponding interaction potential $`V(r_{ij})`$.

Let us introduce the following grid approximation to the Hamiltonian,

$$\stackrel{~}{V}(r)=V(r_\alpha )V_\alpha $$
(2)
for

$$r_\alpha \frac{1}{2M}<r<r_\alpha +\frac{1}{2M}\text{ and }r_\alpha =(\alpha 0.5)r_{cut}/M,$$
(3)
where $`\alpha =1,\mathrm{},M`$, and $`M`$ is the number of grid points within the interval $`[0,r_{cut}]`$, and $`r_{cut}`$ is a chosen cut-off distance. Then, the Hamiltonian in Eq. (1) can be rewritten as

$$H=\underset{\alpha }{}V_\alpha S_\alpha ,$$
(4)
where $`S_\alpha `$ is the number of pairs with interparticle distances inside the $`\alpha `$-slice. Evidently, $`S_\alpha `$ is an estimator of the radial distribution function: $`S_\alpha =4\pi r^2g(r_\alpha )N^2/(2V)`$. The average values of $`S_\alpha `$ are some functions of the potential $`V_\alpha `$ and can be written as an expansion

$$\mathrm{\Delta }S_\alpha =\underset{\gamma }{}\frac{S_\alpha }{V_\gamma }\mathrm{\Delta }V_\gamma +𝒪(\mathrm{\Delta }V^2).$$
(5)
The derivatives $`S_\alpha /V_\gamma `$ can be expressed using exact statistical mechanics relationships LL95

$$\begin{array}{ccc}\frac{S_\alpha }{V_\gamma }\hfill & =& \frac{}{V_\gamma }\frac{\mathrm{dq}S_\alpha (q)\mathrm{exp}(\beta _\lambda K_\lambda S_\lambda (q))}{\mathrm{dq}\mathrm{exp}(\beta _\lambda K_\lambda S_\lambda (q))}\hfill \\ & =& \frac{S_\alpha S_\gamma S_\alpha S_\gamma }{k_BT}.\hfill \end{array}$$
(6)
Equations (5) and (6) allow us to find the interaction potential $`V_\alpha `$ iteratively from the radial distribution functions $`S_\alpha `$. Let $`V_\alpha ^{(0)}`$ be a trial potential for which a natural choice is the potential of mean force (PMF)

$$V_\alpha ^{(0)}=k_BT\mathrm{ln}g(r_\alpha ),$$
(7)
By carrying out standard Monte Carlo simulations, one can evaluate the averages $`S_\alpha `$ and their deviations from the reference values $`S_\alpha ^{}`$ defined from the given radial distribution function as $`\mathrm{\Delta }S_\alpha ^{(0)}=S_\alpha ^{(0)}S_\alpha ^{}`$. By solving the system of linear equations \[Eq. (5)\] with coefficients defined by Eq. (6), and omitting terms $`𝒪(\mathrm{\Delta }V^2)`$, we obtain corrections to the potential $`\mathrm{\Delta }V_\alpha ^{(0)}`$. The procedure is then repeated with the new potential $`V_\alpha ^{(1)}=V_\alpha ^{(0)}+\mathrm{\Delta }V_\alpha ^{(0)}`$ until convergence is reached. The whole procedure resembles a solution of a multidimensional non-linear equation using the Newton-Raphson method.

If the initial approximation of the potential is poor, some regularization of the iteration procedure is needed. In this case we multiply the required change of the radial distribution function by a small factor that is typically between 0 and 1. By doing so, the term $`𝒪(\mathrm{\Delta }V^2)`$ in Eq. (5) can be made small enough to guarantee convergence of the whole procedure, although the number of iterations will increase.

The above algorithm has also the advantage that it provides us with a method to evaluate the uncertainty of the inverse procedure. The radial distribution function has normally some uncertainty. An analysis of the eigenvalues and eigenvectors of the matrix in Eq. (6) allows one to make conclusions of which changes in $`g(r)`$ correspond to which changes in the potential. For example, eigenvectors with eigenvalues close to zero correspond to changes in the potential which have almost negligible effect on the radial distribution function. The presence of these small eigenvalues makes the inverse problem not well-defined, however. In some cases, such as for liquid water, that may pose serious problems in the inversion procedure LL\_CPL00 .

It is instructive to note a relation between the present approach and the renormalization group Monte Carlo method swendsen79 ; pawley84  used to study phase transitions in lattice models (e.g., polymer models and Ising models for ferromagnets) as well as in the quantum field theory. The renormalization procedure introduces a scale change in the system. During this procedure one consecutively obtains a more and more coarse-grained description of the system. Equations (5) and (6) were in fact used by Swendsen et al. swendsen79 ; pawley84  to describe how the parameters of a Hamiltonian change during the renormalization procedure. It now appears that the applications of this method are more general than the original lattice systems near a phase transition, and cover even soft matter problems allowing us to “renormalize” Hamiltonians of molecular systems in such a way that only the degrees of freedom of primary interest are kept in the coarse-grained system.

### II.3 Dissipative Particle Dynamics (DPD)

Dissipative particle dynamics was introduced in 1992 by Hoogerbrugge and Koelman for simulations of hydrodynamic phenomena in complex fluids hoogerbrugge92a ; koelman93a . The original formulation did not obey detailed balance, though, and in 1995 Español and Warren espanol95a  formulated a new DPD algorithm which they showed to be fully consistent with statistical mechanics. This algorithm is nowadays known as DPD. In the following we present the basic DPD formalism and discuss some of its features that are relevant to the present work.

In the basic formulation of DPD, the interactions are pairwise additive and the force exerted by particle $`j`$ on particle $`i`$ is given as a sum of conservative, dissipative, and random forces through $`\stackrel{}{F}_{ij}=\stackrel{}{F}_{ij}^C+\stackrel{}{F}_{ij}^D+\stackrel{}{F}_{ij}^R`$, respectively. These forces are typically given as

$`\stackrel{}{F}_{ij}^C`$ $`=`$ $`F_{ij}^{(c)}(r_{ij})\stackrel{}{e}_{ij},`$ (8)
$`\stackrel{}{F}_{ij}^D`$ $`=`$ $`\gamma \omega ^D(r_{ij})(\stackrel{}{v}_{ij}\stackrel{}{e}_{ij})\stackrel{}{e}_{ij},`$ (9)
$`\stackrel{}{F}_{ij}^R`$ $`=`$ $`\sigma \omega ^R(r_{ij})\xi _{ij}\stackrel{}{e}_{ij},`$ (10)
where $`\stackrel{}{r}_{ij}=\stackrel{}{r}_i\stackrel{}{r}_j`$ is the position vector connecting two particles, $`r_{ij}=|\stackrel{}{r}_{ij}|`$ is the interparticle distance, $`\stackrel{}{e}_{ij}=\stackrel{}{r}_{ij}/r_{ij}`$ is the corresponding unit vector, and $`\stackrel{}{v}_{ij}=\stackrel{}{v}_i\stackrel{}{v}_j`$ is the relative velocity of particles $`i`$ and $`j`$. The terms $`\xi _{ij}`$ are symmetric random variables with zero mean and unit variance, and are independent for different pairs of particles and different times. The conservative force is essentially given by $`F_{ij}^{(c)}(r_{ij})`$. Constants $`\gamma `$ and $`\sigma `$ give the amplitudes of the dissipative and random forces, respectively, and $`\omega ^D`$ and $`\omega ^R`$ are the corresponding weight functions.

Let us now return to the different forces. The DPD formulation does not specify the form of the conservative force. The most common choice in standard DPD simulations is $`F_{ij}^{(c)}(r_{ij})=\alpha _{ij}(1r_{ij})/r_c`$ for $`rr_c`$ and $`F_{ij}^{(c)}(r_{ij})=0`$ otherwise. In other words, the conservative force is derived from a soft pair potential $`U=\alpha _{ij}(1r_{ij}/r_c)^2/2`$, where $`r_c`$ is the cutoff distance. Interactions between different types of coarse-grained particles can then be described by varying the amplitude of the conservative interaction $`\alpha _{ij}`$ groot97a .

The above is the preferred formulation for the conservative force when generic, simple, and soft interactions are adequate. It was shown by Forrest and Suter in 1995 that coarse-graining of molecular representation leads to this type of softening For95 . However, in many cases this formulation is of too generic nature. Although it is possible to use mean field theories such as the Flory-Huggins theory for polymers to find the amplitudes for the conservative forces between different types of particles groot97a , there are cases where this approach does not retain enough information about the actual character of the different atoms, molecules, or interactions (see discussion below).

The weight functions for the dissipative and random forces, $`\omega ^D`$ and $`\omega ^R`$, cannot be chosen arbitrarily. This is easy to understand by the intuitive argument that the thermal heat generated by the random force must be balanced by dissipation. Español and Warren espanol95a  studied this situation analytically and derived a fluctuation-dissipation relation connecting the weight functions and amplitudes of the dissipative and random forces via

$$\omega ^D(r_{ij})=[\omega ^R(r_{ij})]^2\text{ and }\sigma ^2=2\gamma k_BT.$$
(11)
These conditions guarantee convergence to the canonical ensemble as required. It is important to notice that $`\omega ^D`$ and $`\omega ^R`$ are completely decoupled from the conservative part. In other words, we can consider the DPD algorithm as a momentum conserving thermostat for an arbitrary conservative potential. This is in fact the way how we apply the DPD algorithm in our coarse-grained simulations: The conservative potential is defined through the Inverse Monte Carlo procedure and then used in the DPD algorithm.

Finally, to study the dynamics of the system one needs to evolve the system in time. In DPD this is simply done by integrating the Newton’s equations of motion. In our simulations we used a velocity-Verlet based algorithm adapted for DPD (sometimes called DPD–VV) Bes00 ; Vattulainen02 .

### II.4 Combining MD and DPD

With the above definitions we now have all the tools for coarse-graining. We remove fast internal and orientational degrees of freedom from all water molecules, and represent water molecules as one-site particles interacting with other water particles and ions by spherically symmetric potentials. In comparison to the original all-atom model, we have approximately one third of the degrees of freedom left. This can be interpreted as intermediate level of coarse-graining.

Without any loss of generality we can choose the weight functions to be of the standard form $`\omega ^R(r_{ij})=1r_{ij}/r_c`$ as discussed above. This choice has been made invariably in DPD simulations. Since the Inverse Monte Carlo procedure provides us with the effective potentials, thus yielding us $`\stackrel{}{F}_{ij}^C(r)`$, the only task left is to ensure that the MD and DPD systems correspond to each other. In this study, that was done by adjusting the amplitude of the dissipative force $`\gamma `$ which determines the dissipation rate in the DPD system. Here, $`\gamma `$ was determined by requiring the short-time decay of the single-particle velocity autocorrelation function

$$\varphi (t)\stackrel{}{v}_i(t)\stackrel{}{v}_i(0)$$
(12)
to be approximately the same in both MD and DPD systems as shown in Fig. 1. It shows how the early-time decay of $`\varphi (t)`$ can be matched, leading to a value of $`\gamma =0.72`$ for the DPD model. The long-time behavior of the velocity autocorrelation function between MD and DPD is different, however. This is an obvious result since some microscopic degrees of freedom have been coarse grained and thus effects due to hydrogen bonds, for example, are implicitly included in the effective interactions. This leads to an enhanced diffusion rate at intermediate times, and is demonstrated in Fig. 1 as a positive tail for the DPD particles.

This approach is meaningful since (i ) it makes sure that the early-time dynamics is described properly, while (ii ) it does not fix the tracer diffusion coefficient. One should notice that the tracer diffusion coefficient is an integral over the velocity correlation function over long times, and in hydrodynamic systems the long-time tail is important. With this choice, we can assume the diffusion coefficient to be an independent quantity so that its behavior, found by DPD, can be compared with both MD simulations and experiments. Fixing $`\gamma `$ thus determines the transport properties of the system.

## III Results

### III.1 Constructing the potentials

One of the typical simplifications used in molecular simulations is the replacement of solvent molecules by continuum media. For example, in the primitive electrolyte model, ions in water are substituted by charged spheres moving in dielectric media with the dielectric constant set to about 80. This is a serious simplification at small distances (a few ångströms) where it is impossible to define a dielectric constant. Besides this, in the primitive electrolyte model the ion radius is an adjustable parameter without any clear physical meaning. A better model for effective ion-ion interactions in an aqueous solution must take into account the solvation structure around the ions. Practically, effective solvent-mediated ion-ion potentials may be constructed by the Inverse Monte Carlo method from ion-ion radial distribution functions generated in high-quality all-atomic molecular dynamics simulations LL95 ; LyMar02 ; LL97 . This is the approach we used in obtaining the effective potentials.

The MD simulations were performed in the NVT ensemble using the flexible SPC water model tourah  and Smith–Dang parameters for Na<sup>+</sup> and Cl<sup>-</sup> ions SD94 , i.e., the Lennard-Jones parameters for the sodium ions were $`\sigma =2.35`$ Å and $`\epsilon =0.544`$ kJ/mol, and for chloride $`\sigma =4.4`$ Å and $`\epsilon =0.419`$ kJ/mol. In the results reported here we set the temperature to 300 K. The electrostatic interactions were taken into account by Ewald summation. Other simulation details are described in Refs. LL95 ; LL97 . From the MD simulations, the radial distribution functions between different pairs of particles were calculated and fed as an input in the Inverse Monte Carlo procedure. The results are shown in Figs. 2 and 3.

Figure 2 shows effective potentials between water molecules, presented as one-center particles, and between other water molecules and ions. For comparison, water-water effective potential calculated from MD simulations for pure water is also displayed. It is interesting to note that the presence of ions has practically no effect on the water-water effective potential.

Figure 3 displays the effective potentials between the ions. They are compared to the effective potentials calculated without water LL95 ; LL97 . The effective potentials are very similar in the two cases both of them having an oscillating character with one or two oscillations. The potential approaches the primitive model potential with dielectric constant close to 80. At distances above 10 Å the effective potential coincides almost perfectly with the Coulombic potential.

Within the inverse Monte Carlo procedure, calculation of effective potentials without water is much easier than in the presence of water. On the other hand, the presence of particles representing water is necessary in DPD simulations. Our test studies have shown that while the use of effective potentials calculated without water may give qualitatively satisfactory results, inclusion of water in the Inverse Monte Carlo simulations essentially leads to an improvement of results, especially in the case of dynamical properties.

It is worth pointing out that the effective potential differs from the potential of mean force defined as $`\psi _{PMF}=k_BT\mathrm{ln}g(r)`$, which corresponds to the Kirkwood approximation for the $`n`$-particle correlation function. The potential of mean force is screened by the other ions in a system and decays as $`(1/r)e^{(r/r_D)}`$, i.e., not as the $`1/r`$ Coulombic potential. The potential of mean force is therefore not a very suitable choice to present ion-ion interactions within a continuum solvent model at finite ion concentrations even in a coarse-grained solvent system. The effective potentials include a contribution from Coulombic forces, however.

### III.2 Coarse-grained simulations

In our first paper using this approach karttunen02 , we showed that the radial distribution functions computed from the coarse-grained DPD simulations for the pairs between the Na<sup>+</sup> and Cl<sup>-</sup> ions agreed with those obtained from the all-atom MD simulations. In addition to the static properties, the diffusion coefficients were shown to be in good agreement between the MD and DPD results at various NaCl concentrations between 0.5 M and 4.1 M. What makes those results very interesting is the fact that, for all the DPD simulations, we used the effective potentials obtained at one single NaCl concentration only. In that case, the effective potentials were obtained from MD simulations at a salt concentration of 0.87 M and then used in DPD simulations for various salt concentrations between 0.5 M and 4.1 M. When compared with MD simulations, the agreement was found to be excellent. Here, we present results for static correlations, coordination numbers, velocity autocorrelation functions, and viscosity.

It is important to notice that due to the self-consistency conditions for the effective potentials, the DPD simulations will always produce the correct pair correlation behavior at the point in phase space where the effective interactions were determined in the first place. As an example, this is demonstrated in Fig. 4 for all different kinds of pairs of particles studied in this work. It is very relevant to ask, however, how the effective interactions change when system variables such as the salt concentration or the temperature of the system change. We have explored this matter through combined MD and DPD simulations and found that for a wide range of salt concentrations in the present system it is not necessary to recalculate the effective potentials. This is demonstrated by the data in Table 1 which compares the coordination numbers obtained from DPD and MD simulations for $`g(r)`$. In the DPD simulations we have again used the effective potentials obtained at 0.87 M. We find that both MD and DPD follow the same qualitative behavior. For sodium the coordination number decreases with increasing salt concentration whereas for chloride we find an opposite trend. Furthermore, the quantitative agreement between the MD and DPD results is remarkably good.

Our finding that the effective interactions depend only weakly on thermodynamic conditions is also supported by the present results for shear viscosity (see below) and earlier studies of the same system for tracer diffusion karttunen02 . In Ref. karttunen02 , we found that DPD simulations with effective potentials yielded tracer diffusion coefficients in good agreement with MD simulations. Further support for this result is given by previous work LyMar02 ; LL97 , where the dependence of effective potentials on the salt concentration was studied in detail. It was found that effective potentials depend very little on the salt concentration: An increase in salt concentration leads to a slow decrease in the dielectric permittivity.

Our results thus suggest that the effective interactions change only slightly when the system parameters are varied, as is here the case for salt concentration. However, while this result holds true in the present system for an aqueous solution of monovalent ions, the generality of this finding remains an open question. Further studies of divalent systems and more complex molecules, among others, are therefore called for.

Next, we determined the shear viscosity using a Green-Kubo relation. As with the other quantities, the shear viscosity should be compared to results from the MD simulations. However, since the shear viscosity coefficient is a collective quantity, meaning that all particles in a system give rise to a single sample, an accurate determination of the shear viscosity coefficient through MD simulations is both difficult and very time-consuming. Thus it was not done in the present case, and as a matter of fact, to the best of our knowledge, there are no published reports of MD simulations of shear viscosity in NaCl solutions. Thus, we compare our DPD data to experimental results crc  instead. The results shown in Fig. 5 indicate that the qualitative behavior of the shear viscosity coefficient obtained through DPD simulations is in good agreement with the experiments. This result is truly remarkable as it shows that both the equilibrium (static) and dynamic behavior of the system are properly described by the coarse-grained approach. This provides us with strong support that the present approach of coarse graining an MD system and applying the obtained effective interactions in DPD simulations is a promising approach indeed.

We would like to note that we have very recently performed simulations using LiCl and CaCl<sub>2</sub> terama02  to study further aspects of this issue, such as the effects of salt concentration, valence, and temperature. Despite the very different natures of these solutions (as compared to NaCl), the conclusions with regard to the coarse-graining method and its general applicability remain the same.

## IV Conclusions

In this article we have described in detail a method that combines the Inverse Monte Carlo method and the momentum conserving dissipative particle dynamics thermostat for coarse-grained simulations. This approach allows systematic coarse-graining of molecular systems and conserves hydrodynamics. Here, the level of coarse-graining has been quite modest and the potentials used did not have particularly soft cores. This is due to the fact that in this system of water and NaCl, we coarse-grained the molecular representation of water molecules only. Thus, the corresponding effective interactions retained features which are typical for systems described on an atomic level. However, it is important to notice two facts. First of all, further coarse-graining is not possible for this kind of a system if one wishes to preserve the identity of individual ions and the essential characteristics of water molecules. One could, of course, continue the coarse-graining procedure to obtain softer potentials by clustering several water molecules together and treating the clusters as new “effective” solvent particles. This would lead to further computational gains but the ions and molecules would then lose essential parts of their identity. As a result of this the interactions would become softer and resemble the quadratic potentials used for conservative interactions in the original DPD formulation. Second, despite the modest level of coarse-graining in this work, the computational gain is still remarkable, around a factor of 20. This speed-up is due to the fact that no explicit Ewald summation was required for the coarse-grained system (since the effective potentials include information about the electrostatics, and the screening length was shorter than the applied cut-off) and that SPC water used in molecular dynamics simulations is replaced by a single coarse-grained particle with the electrostatic interactions included in the effective potentials.

This approach allows one to coarse-grain complex macromolecular systems to a level where a detailed molecular description is replaced with a simplified model of coarse-grained particles. For example, a lipid-acyl chain “$`\mathrm{}`$ – CH<sub>2</sub> – CH<sub>2</sub> – CH<sub>2</sub> – CH<sub>2</sub> – CH<sub>2</sub> – CH<sub>2</sub> $`\mathrm{}`$” can be described as “$`\mathrm{}`$– B – B – $`\mathrm{}`$” where each of the beads “B” represent 2 – 4 methyl groups CH<sub>2</sub>. Since the beads “B” describe the interaction site for a whole molecular group, it is clear that the interactions between different beads will be relatively soft in the spirit of the original DPD formulation, yet retaining their characteristic features well beyond the mean-field description of the original DPD formulation. For macromolecular systems, this coarse-graining procedure thus has the potential to lead to a major reduction in computational effort as well as to softer interaction potentials. The computational gain depends on the type of interactions: While here in the case of NaCl we did not need to use Ewald summation, the situation may be different in the case of divalent ions such as calcium. In addition, one should note that coarse-graining of a macromolecular system requires a lot of insight, and there is no unique way of selecting the units to be coarse-grained. Work is in progress to address these issues and the coarse-graining of macromolecular systems in general.

The approach presented in this article has close connections to other methods as well. The Inverse Monte Carlo procedure has the same foundations as the integral equation methods of Bolhuis et al. louis00a ; bolhuis01a . Second, the coarse-graining procedure used by Reith, Mayer and Müller-Plathe reith00a  has a close resemblance to the Inverse Monte Carlo procedure. The latter group has also applied their method to coarse-grain polymeric systems in the spirit described in the previous paragraph. The novelty here is the application of the effective potentials in coarse-grained simulations by using the dissipative particle dynamics thermostat. The results for the present system indicate that dynamic properties could be studied using this approach. Whether or not this holds true in more complex systems, and what are the limits of this approach are topics of further studies.

In conclusion, the agreement between the DPD simulations, the MD results and experiments strongly indicates that the method presented in this work is both physically sound and provides a controllable method for performing coarse-graining and mesoscale simulations at least to a moderate degree.

## Acknowledgments

This work has been supported by the Swedish Science Council (A.P.L. and A.L.) and the Academy of Finland (M.K.). Further support has been obtained from the Academy of Finland though its Centre of Excellence Program (I.V.) and from the Finnish Academy of Science and Letters (I.V.).


format: text

text:
# A Compact Homogeneous S-space11footnote 1 2000 Mathematics Subject Classification: Primary 54G20, 54D30. Key Words and Phrases: Compact group, pointwise topology.

## 1 Introduction

A space $`X`$ is hereditarily separable ($`\mathrm{𝐻𝑆}`$) iff every subspace is separable. An $`S`$-space is a regular Hausdorff $`\mathrm{𝐻𝑆}`$ space with a non-Lindelöf subspace. A space $`X`$ is homogeneous iff for every $`x,yX`$ there is a homeomorphism $`f`$ of $`X`$ onto $`X`$ with $`f(x)=y`$. Under $`\mathrm{𝐶𝐻}`$, several examples of $`S`$-spaces have been constructed, including topological groups (see ) and compact $`S`$-spaces (see ). It is asked in  whether there are compact homogeneous $`S`$-spaces. As we shall show in Theorem 4.2, there are under $`\mathrm{𝐶𝐻}`$. This cannot be done in $`\mathrm{𝑍𝐹𝐶}`$, since there are no compact $`S`$-spaces under $`\mathrm{𝑀𝐴}+\neg \mathrm{𝐶𝐻}`$ (see ); there are no $`S`$-spaces at all under $`\mathrm{𝑃𝐹𝐴}`$ (see ).

In Section 2, we use a slightly modified version of the construction in  to refine the topology of any given second countable space, and turn it into a first countable strong $`S`$-space (i.e., each of its finite powers is an $`S`$-space). In Section 3, we show that if the original space is compact, then there is a natural compactification of the new space which is also a first countable strong S-space. If in addition the original space is zero-dimensional, then the $`\omega ^{\mathrm{th}}`$ power of this compactification will be homogeneous by Motorov , proving Theorem 4.2.

## 2 A Strong S-Space

If $`\tau `$ is a topology on $`X`$, we write $`\tau ^I`$ for the corresponding product topology on $`X^I`$; likewise if $`\tau ^{}\tau `$ is a base we write $`(\tau ^{})^I`$ for the natural corresponding base for $`\tau ^I`$. If $`EX`$, then $`\mathrm{cl}(E,\tau )`$ denotes the closure of $`E`$ with respect to the topology $`\tau `$. This notation will be used when we are discussing two different topologies on the same set $`X`$.

The following two lemmas are well-known; the second is Lemma 7.2 in :

###### Lemma 2.1

If $`X`$ is $`\mathrm{𝐻𝑆}`$ and $`Y`$ is second countable, then $`X\times Y`$ is $`\mathrm{𝐻𝑆}`$

###### Lemma 2.2

$`X^\omega `$ is $`\mathrm{𝐻𝑆}`$ iff $`X^n`$ is $`\mathrm{𝐻𝑆}`$ for all $`n<\omega `$.

The next lemma, an easy exercise, is used in the proof of Theorem 2.4:

###### Lemma 2.3

If $`(x,y)X\times Y`$ and $`SX\times Y`$, then $`(x,y)\mathrm{cl}(S)`$ iff $`y\mathrm{cl}(\pi (S(U\times Y)))`$ for all neighborhoods $`U`$ of $`x`$, where $`\pi :X\times YY`$ is projection.

The following is proved (essentially) in , but our proof below may be a bit simpler:

###### Theorem 2.4

Assume $`\mathrm{𝐶𝐻}`$. Let $`\rho `$ be a second countable $`T_3`$ topology on $`X`$, where $`|X|=\mathrm{}_1`$. Then there is a finer topology $`\tau `$ on $`X`$ such that $`(\omega _1,\tau )`$ is a first countable locally compact strong S-space.

Proof. WLOG, $`X=\omega _1`$. For $`\eta <\omega _1`$ we write $`\rho _\eta `$ for the topology of $`\eta `$ as a subspace of $`(\omega _1,\rho )`$. Applying $`\mathrm{𝐶𝐻}`$, list $`_{0<n<\omega }[(\omega _1)^n]^\omega `$ as $`\{S_\mu :\mu \omega _1\}`$, so that each $`S_\mu \mu ^{n(\mu )}`$ for some $`n(\mu )`$ with $`0<n(\mu )<\omega `$.

For $`\eta \omega _1`$ we construct $`\tau _\eta `$ a topology on $`\eta `$ by induction on $`\eta `$ so as to make the following hold for all $`\xi <\eta \omega _1`$:

* $`\tau _\xi =\tau _\eta 𝒫\left(\xi \right)`$.
* $`\tau _\eta `$ is first countable, locally compact, and $`T_3`$.
* $`\tau _\eta \rho _\eta `$.

Note that (1) implies in particular that $`\xi \tau _\eta `$; that is, $`\xi `$ is open. Thus, if $`\tau =\tau _{\omega _1}`$, then $`(\omega _1,\tau )`$ is not Lindelöf. Also by (1), $`\tau _\eta `$ for limit $`\eta `$ is determined from the $`\tau _\xi `$ for $`\xi <\eta `$. So, we need only specify what happens at successor ordinals.

For $`n1`$ and $`\xi <\omega _1`$, let $`\mathrm{𝐼𝑠𝑒𝑞}(n,\xi )`$ be the set of all $`f(\omega _1)^n`$ which satisfy $`f(0)<f(1)<\mathrm{}<f(n1)=\xi `$. The following condition states our requirement on $`\tau _{\xi +1}`$:

* For each $`\mu <\xi `$ and each $`f\mathrm{𝐼𝑠𝑒𝑞}(n,\xi )`$, where $`n=n_\mu `$:

$$f\mathrm{cl}(S_\mu ,(\tau _{\xi +1})^{n1}\times \rho )f\mathrm{cl}(S_\mu ,(\tau _{\xi +1})^n).$$

If $`n=n_\mu =1`$, then $`(\tau _{\xi +1})^{n1}\times \rho `$ just denotes $`\rho `$. That is, (4) requires

$$\xi \mathrm{cl}(E,\rho )\xi \mathrm{cl}(E,\tau _{\xi +1})$$
$`()`$
for all $`E`$ in the countable family $`\{S_\mu :\mu <\xi \&n(\mu )=1\}`$. It is standard (see ) that one may define $`\tau _{\xi +1}`$ so that this holds. Now, consider (4) in the case $`n=n_\mu 2`$. By (2), $`\tau _\xi `$ is second countable, so let $`\tau _\xi ^{}`$ be a countable base for $`\tau _\xi `$. Applying Lemma 2.3, (4) will hold if whenever $`U=U_0\times \mathrm{}\times U_{n2}(\tau _\xi ^{})^{n1}`$ is a neighborhood of $`f\mathrm{}(n1)`$,

$$\xi \mathrm{cl}(\pi (S_\mu (U\times (\xi +1))),\rho )\xi \mathrm{cl}(\pi (S_\mu (U\times (\xi +1))),\tau _{\xi +1}),$$

where $`\pi :\xi ^{n1}\times (\xi +1)(\xi +1)`$ is projection. But this is just a requirement of the form $`()`$ for countably many more sets $`E`$, so again there is no problem meeting it.

Now, we need to show that $`\tau ^n`$ is $`\mathrm{𝐻𝑆}`$ for each $`0<n<\omega `$. We proceed by induction, so assume that $`\tau ^m`$ is $`\mathrm{𝐻𝑆}`$ for all $`m<n`$. Fix $`A(\omega _1)^n`$; we need to show that $`A`$ is $`\tau ^n`$-separable. Applying the induction hypothesis, we may assume that each $`fA`$ has all coordinates distinct. Also, since permutation of coordinates induces a homeomorphism of $`(\omega _1)^n`$, we may assume that each $`fA`$ is strictly increasing; that is, $`f\mathrm{𝐼𝑠𝑒𝑞}(n,\xi )`$, where $`\xi =f(n1)`$. By the induction hypothesis and Lemma 2.1, $`A`$ is separable in $`(\tau _{\xi +1})^{n1}\times \rho `$. We can then fix $`\mu `$ such that $`n(\mu )=n`$, $`S_\mu A`$, and $`S_\mu `$ is $`(\tau _{\xi +1})^{n1}\times \rho `$-dense in $`A`$. Now, say $`fA`$ with $`\xi =f(n1)>\mu `$. Applying (4), we have $`f\mathrm{cl}(S_\mu ,\tau ^n)`$. Thus, $`A\mathrm{cl}(S_\mu ,\tau ^n)`$ is countable, so $`A`$ is $`\tau ^n`$-separable. $`\mathrm{}`$

## 3 Compactification

###### Definition 3.1

If $`\phi `$ is a continuous map from the $`T_2`$ space $`Y`$ into $`X`$, then $`Y\dot{}_\phi X`$ denotes the disjoint union of $`X`$ and $`Y`$, given the topology which has as a base:

* All open subsets of $`Y`$, together with
* All $`[U,K]:=U(\phi ^1UK)`$, where $`U`$ is open in $`X`$ and $`K`$ is compact in $`Y`$.

Our main interest here is in the case where $`X`$ is compact and $`Y`$ is locally compact. Then, if $`|X|=1`$, we have the 1-point compactification of $`Y`$, and if $`Y`$ is discrete and $`\phi `$ is a bijection we have the Aleksandrov duplicate of $`X`$.

###### Lemma 3.2

Let $`Z=Y\dot{}_\phi X`$, with $`X,Y`$ Hausdorff:

* $`X`$ is closed in $`Z`$, $`Y`$ is open in $`Z`$, and both $`X,Y`$ inherit their original topology as subspaces of $`Z`$.
* If $`Y`$ is locally compact, then $`Z`$ is Hausdorff.
* If $`X`$ is compact, then $`Z`$ is compact.
* If $`X,Y`$ are first countable, $`X`$ is compact, $`Y`$ is locally compact, and each $`\phi ^1(x)`$ is compact, then $`Z`$ is first countable.
* If $`X,Y`$ are zero dimensional, $`X`$ is compact, and $`Y`$ is locally compact, then $`Z`$ is zero dimensional.
* If $`X`$ is second countable and $`Y^\omega `$ is $`\mathrm{𝐻𝑆}`$, then $`Z^\omega `$ is $`\mathrm{𝐻𝑆}`$.

Proof. For (3): If $`𝒰`$ is a basic open cover of $`Z`$, then there are $`n\omega `$ and $`[U_i,K_i]𝒰`$ for $`i<n`$ such that $`_{i<n}U_i=X`$. Thus, $`_{i<n}[U_i,K_i]`$ contains all points of $`Z`$ except for (possibly) the points in the compact set $`_{i<n}K_iY`$.

For (4): $`Z`$ is compact Hausdorff and of countable pseudocharacter.

For (5): $`Z`$ is compact Hausdorff and totally disconnected.

For (6): By Lemma 2.2, it is sufficient to prove that each $`Z^n`$ is $`\mathrm{𝐻𝑆}`$. But $`Z^n`$ is a finite union of subspaces of the form $`X^j\times Y^k`$, which are $`\mathrm{𝐻𝑆}`$ by Lemma 2.1. $`\mathrm{}`$

## 4 Homogeneity

The following was proved by Dow and Pearl :

###### Theorem 4.1

If $`Z`$ is first countable and zero dimensional, then $`Z^\omega `$ is homogeneous.

Actually, we only need here the special case of this result where $`Z`$ is compact and has a dense set of isolated points; this was announced (without proof) earlier by Motorov .

Note that by Šapirovskiĭ , any compact $`\mathrm{𝐻𝑆}`$ space must have countable $`\pi `$-weight (see also , Theorem 7.14), so if it is also homogeneous, it must have size at most $`2^\mathrm{}_0`$ by van Douwen . Under $`\mathrm{𝐶𝐻}`$ this implies, by the Čech – Pospíšil Theorem, that the space must be first countable.

###### Theorem 4.2

$`\left(\mathrm{𝐶𝐻}\right)`$ There is a (necessarily first countable) zero-dimensional compact homogeneous strong $`S`$-space.

Proof. Let $`X`$ be the Cantor set $`2^\omega `$ with its usual topology, let $`Y`$ be $`2^\omega `$ with the topology constructed in Theorem 2.4, let $`\phi `$ be the identity, and let $`Z=Y\dot{}_\phi X`$. By Lemma 3.2, $`Z`$, and hence also $`Z^\omega `$, are zero-dimensional first countable compact strong S-spaces; $`Z^\omega `$ is homogeneous by Theorem 4.1. $`\mathrm{}`$

No compact topological group can be an $`S`$-space or an $`L`$-space. However under $`\mathrm{𝐶𝐻}`$ there are, by , compact $`L`$-spaces which are right topological groups (i.e. they admit a group operation such that multiplication on the right by a fixed element defines a continuous map). We do not know whether there can be compact $`S`$-spaces which are right topological groups.


text:
# The Quark-Gluon Plasma in a Finite Volume

## I Introduction

The study of the quark gluon plasma in the laboratory through hadronic collisions represents one of the largest current initiatives in nuclear physics. Although these volumes are manifestly finite, modeling of the collisions seldom addresses the impact of the finiteness as most treatments are based on the the grand canonical ensemble where the conservation constraints of baryon number, strangeness, isospin and color are ignored. Indeed, one expects that these constraints would become meaningless for larger volumes. The effective volume for a $`pp`$ collision might be only a few cubic fm, while the volume of a RHIC collision might be a few thousand cubic fm. But even in a RHIC collision, the duration of the quark gluon phase might be only a few fm/c, meaning that the effective volume for charge conservation may be a few dozen cubic fm.

Recently, the effects of local charge conservation have gained more attention due to the relationship of charge conservation to charge fluctuations and charge balance functionsbassdanpratt ; jeonkoch ; asakawaprl ; asakawa ; gavinqm2002 ; gavinkapusta ; pruneaugavinvoloshin ; shuryakstephanov ; stephanovqm2001 . These observables are intimately related to the delay of hadronization expected should a quark-gluon plasma be created. Recent papers have focused on the importance of conserving electric charge, strangeness and baryon number, as well as the dynamics of charge equilibration kochlin . The statistical mechanics of a finite-volume hadron gas have been studied in the context of pion flavor distributionssenisofluc .

The approaches for solving the canonical ensemble can be divided into two classes. In the first set of approaches, one can sum over all partition functions in the grand canonical ensemble assuming the chemical potentials are imaginary, $`Z\mathrm{Tr}e^{\beta H+i\beta \mu Q}`$. By integrating over all $`\mu `$, the phase factor becomes a delta function which restricts the phase space to states with $`Q=0`$ becattini ; turko ; redlichturko . This approach can even be applied to non-Abelian symmetries such as SU(N)elzeplb ; elzepra . By associating $`N1`$ charges, $`Q_1Q_{N1}`$, with the $`N1`$-fold Cartan subgroup, one can restrict the ensemble to states with specific charges $`Q`$, which can then be associated with the SU(3) multiplet labels, $`(p,q)`$, by a transformation. These techniques work well for simple systems, e.g. massless partons, where the integral over complex phases can be be performed analytically. Multiple charges can be readily incorporated, including conservation of overall momentum.

The second class of approaches for the canonical ensemble centers about recursion relations which were proposed by Chase and Mekjianmekjian ; dasguptamekjian  for the study of nuclear fragmentation. These methods have been extended to include quantum statisticsprattfermi  and non-additive chargessenisofluc , such as angular momentum and isospin. In addition to partition functions, these methods can also generate multiplicity distributionsprattdasgupta ; senisofluc . The advantage of these methods is that, given the one-particle partition functions, they provide exact answers. Since these methods involve sums, rather than integration over phases, they provide robust answers for arbitrary energy levels. Furthermore, symmetrization is included to all orders. If the number of particles to be considered remains below 100, numerical calculations tend to take only seconds or minutes at most.

We review recursive techniques for canonical ensembles and present extensions to incorporate conservation of color in the next section. In the subsequent section a simple example of a non-interacting quark-gluon gas is explored. Canonical partition functions are generated which enforce conservation of baryon number, strangeness and charge, while requiring the system to be in an isosinglet and in a color singlet. Bose and Fermi statistics are included to all order. Energy, entropy, and particle number densities are all shown to be significantly suppressed for small volumes, less than 5 fm3. Arguments are made that these effects are unlikely to be important in a relativistic heavy ion collision, but might be quite important for $`pp`$ collisions.

## II Recursive techniques for quantum colored particles

For additive charges, canonical partition functions can be calculated in a straight-forward manner. Neglecting symmetrization, conservation of a set of charges $`\stackrel{}{Q}`$ can be enforced through the recursion relation,

$$Z_{A,\stackrel{}{Q}}(T)=\underset{k}{}\frac{a_k\omega _k(T)}{A}Z_{Aa_k,\stackrel{}{Q}\stackrel{}{q}_k}(T).$$
(1)
Here, $`k`$ refers to all particles species of charge $`\stackrel{}{q}_k`$ whose single-particle partition functions are:

$$\omega _k=\underset{i}{}d_ie^{E_i/T},$$
(2)
where $`i`$ denotes the individual energy levels of energy $`E_i`$ and degeneracy $`d_i`$. The number $`a_k`$ must be positive for every species. For low energy nuclear physics applications, $`a_k`$ could be the baryon number and $`A`$ would represent one of the conserved charges. For a high-energy system where anti-particles must be included and there are no positive-definite conserved charges, $`A`$ can refer to the number of hadrons. This number is positive-definite, but as the number of hadrons is not conserved, one must consider all $`A`$ to obtain the partition function, $`Z=_AZ_A`$.

Symmetrization can be incorporated by considering all permutations of identical particles prattfermi ; prattmultibose ,

$$Z_{A,\stackrel{}{Q}}=\frac{1}{A!}\underset{j_1\mathrm{}j_A,𝒫(i)}{}j_1\mathrm{}j_A|e^{H/T}|𝒫(j_1\mathrm{}j_A)(1)^{N_P(𝒫)},$$
(3)
where a sum over all $`A!`$ permutations of the single particle levels accounts for the symmetrization and $`N_P(𝒫)`$ counts the number of pair-wise interchanges of fermions involved in the permutation $`𝒫`$.

By considering the cyclic permutation,

$$\begin{array}{cc}\hfill C^{k,\mathrm{}}(T)& \underset{j_1\mathrm{}j_{\mathrm{}}}{}j_1,j_2\mathrm{}j_{\mathrm{}}|e^{H/T}|j_2,\mathrm{}j_{\mathrm{}},j_1\hfill \\ \hfill =& \underset{i}{}e^{\mathrm{}E_i/T}=\omega _k(T/\mathrm{}),\hfill \end{array}$$
(4)
one can derive a modified recursion relation,

$$Z_{A,\stackrel{}{Q}}(T)=\underset{k}{}\underset{\mathrm{}}{}\frac{a_kC^{k,\mathrm{}}(T)}{A}(\pm 1)^{\mathrm{}+1}Z_{A\mathrm{}a_k,\stackrel{}{Q\mathrm{}q_k}}(T).$$
(5)
By summing over all possible values of $`\mathrm{}`$, symmetrization is included to all orders.

Before addressing SU(3) symmetries, we review the techniques used for SU(2) as described in senisofluc . Extending the techniques described above to incorporate an SU(2) symmetry such as isospin or angular momentum can be accomplished in two ways. First, one can consider the projection $`M`$ as an additive charge and solve for $`Z_M(T)`$ with the method outlined above. Given the $`(2I+1)`$ degeneracy of all multiplets of size $`I`$, the partition function for states with fixed $`I`$ is easily obtained from the partition function for fixed $`M`$, $`Z_I=Z_MZ_{M+1}`$. The second technique requires one to understand the decomposition of the cycle diagrams into SU(2) multiplets. In terms of the projection operator, $`𝒪_M^I`$, which projects states with total isospin $`I`$ and projection $`M`$, one can define the isospin decomposition of the cycle diagram.

$$\begin{array}{cc}\hfill C_i^{\kappa ,\mathrm{}}(T)& \underset{j_1\mathrm{}j_{\mathrm{}}}{}j_1,j_2\mathrm{}j_{\mathrm{}}|e^{H/T}𝒪_m^i|j_2,\mathrm{}j_{\mathrm{}},j_1\hfill \\ \hfill =& \omega _\kappa (T/\mathrm{})\chi _i^{\kappa ,\mathrm{}},\hfill \end{array}$$
(6)
where $`\chi _i^{\kappa ,\mathrm{}}`$ is the operator $`C_k^{\mathrm{},i}`$ evaluated with the particles being combined to one quantum state of zero energy. The subscript $`\kappa `$ refers to a specific multiplet, not the individual states. For instance, $`\kappa `$ might refer to pions and $`\omega _\kappa (T)`$ would be the partition function for a given species of pion. The details about the energy levels are absorbed in $`\omega _\kappa `$ while $`\chi `$ is independent of the energy levels, but does depend on the isospin of the particular species, $`i_\kappa `$.

Once $`\chi _i^{\kappa ,\mathrm{}}`$ is known, one can calculate the partition function for fixed $`I`$ recursively.

$$Z_{A,\stackrel{}{Q},I}(T)=\underset{\kappa }{}\underset{\mathrm{}}{}\frac{a_k\omega _k(T/\mathrm{})}{A}(\pm 1)^{\mathrm{}+1}Z_{A\mathrm{}a_k,\stackrel{}{Q}\mathrm{}\stackrel{}{q}_k,I^{}}(T)\chi _i^{\kappa ,\mathrm{}}𝒩(I^{},i;I).$$
(7)
One could have added indices $`M`$, $`M^{}`$ and $`m`$ onto $`Z_{A,I}`$, $`Z_{A^{},I^{}}`$ and $`\chi _i^{\kappa ,\mathrm{}}`$ respectively. However, as each of these functions involve a trace over all states, and the Hamiltonian has no dependence on the isospin projection, none of these quantities has any dependence on the projection and the indices can be suppressed. The function $`𝒩(I^{},i,I)`$ counts the number of multiplets of size $`I`$ that one obtains by coupling $`I^{}`$ with $`i`$. For SU(2), these obey a simple form.

$$\begin{array}{cc}\hfill 𝒩(I^{},i;I)=& \underset{M^{},m}{}I,M|I^{},M^{},i,m^2\hfill \\ \hfill =& \{\begin{array}{cc}1,& |I^{}i|II^{}+i\\ 0,& \mathrm{otherwise}\end{array}\hfill \end{array}$$
(8)
Before performing the recursive calculations described in Eq. (7), one needs a simple expression for $`\chi _i^{\kappa ,\mathrm{}}`$. This can be accomplished by realizing that by summing over $`i`$ one obtains the projection of states with fixed $`m`$.

$$\underset{i}{}\chi _{i,m}^{\kappa ,\mathrm{}}=\underset{j_1\mathrm{}j_{\mathrm{}}}{}j_1,j_2\mathrm{}j_{\mathrm{}}|𝒪_m|j_2\mathrm{}j_{\mathrm{}},j_1.$$
(9)
Since the states $`|j_1\mathrm{}j_{\mathrm{}}`$ are eigenstates of $`𝒪_m`$, the bra and ket must be identical, i.e., all particles must be identical and have the same isospin projection.

$$\underset{i}{}\chi _{\kappa ,m}^{\mathrm{},i}=\{\begin{array}{cc}1,& m=0,\pm \mathrm{},\pm 2\mathrm{},\mathrm{}\pm \mathrm{}i_\kappa \\ 0,& \mathrm{otherwise}\end{array}$$
(10)
Since $`\chi _{i,m}^{\kappa ,\mathrm{}}`$ is independent of $`m`$ when $`m`$ is a member of the isomultiplet $`i`$, the sum over $`i`$ in Eq. (10) can be limited to $`im`$. It is then straight-forward to see that

$$\chi _i^{\kappa ,\mathrm{}}=\{\begin{array}{cc}\hfill 1,& i=0,\mathrm{},2\mathrm{},\mathrm{}\mathrm{}i_\kappa \\ \hfill 1,& i=\mathrm{}1,2\mathrm{}1,\mathrm{}\mathrm{}i_\kappa 1\\ \hfill 0,& \mathrm{otherwise}\end{array}$$
(11)
Here, the index $`m`$ is suppressed since $`\chi `$ is independent of $`m`$ as long as $`m`$ is a member of the isomultiplet $`i`$. For the $`i_\kappa =1`$ case, e.g., pions or $`\rho `$ mesons, $`\chi _i^{\mathrm{}}=1`$ for $`i=\mathrm{}`$ and $`i=0`$, and $`\chi _i^{\mathrm{}}=1`$ for $`i=\mathrm{}1`$.

Restricting SU(3) color can be accomplished with the same procedure. Following the same steps as the SU(2) case, one can derive analogous expressions,

$$Z_{A,\stackrel{}{Q},(P,Q)}(T)=\underset{\kappa ,\mathrm{},(P^{},Q),(p,q)}{}\frac{a_\kappa \omega _\kappa (T/\mathrm{})}{A}(\pm 1)^{\mathrm{}+1}Z_{A\mathrm{}a_\kappa ,\stackrel{}{Q}\mathrm{}\stackrel{}{q}_\kappa ,(P^{},Q^{})}(T)\chi _{(p,q)}^{\kappa ,\mathrm{}}𝒩((P^{},Q^{}),(p,q);(P,Q)).$$
(12)
Here, $`(P,Q)`$ denotes an SU(3) multiplet, e.g., the gluon color octet is represented by $`(P=1,Q=1)`$. For an explanation fo the notation see youngtableauxref . Again, the challenge in making Eq. (12) tractable is in finding expressions for $`\chi _{(p,q)}^{\kappa ,\mathrm{}}`$ and $`𝒩((P^{},Q^{}),(p,q);(P,Q))`$. Coupling SU(3) multiplets, i.e., finding expressions for $`𝒩`$, follows rules based on manipulating Young tableaus youngtableauxref . For the calculations in this paper, these rules were programmed numerically.

Finding an expression for $`\chi _{(p,q)}^{\kappa ,\mathrm{}}`$ requires performing a color decomposition of the cycle diagrams and can be done analogously as was shown above for SU(2). First, one must project states of a given $`\mu `$, where $`\mu `$ represents the eigenvalues of the projection operators in SU(3), e.g., hyper-charge and $`I_3`$ for SU(3) flavor. As with the SU(2) example above, this projection is realized by summing $`\chi _{(p,q),\mu }^{\kappa ,\mathrm{}}`$ over all $`(p,q)`$ multiplets which include the projection $`\mu `$.

$$\underset{(p,q)}{}\chi _{(p,q),\mu }^{\kappa ,\mathrm{}}=\underset{j_1\mathrm{}j_{\mathrm{}}}{}j_1,j_2\mathrm{}j_{\mathrm{}}|𝒪_\mu |j_2\mathrm{}j_{\mathrm{}},j_1.$$
(13)
Again, the $`\mu `$ dependence in $`\chi `$ represents the known degeneracy of $`\mu `$ within a given $`(p,q)`$ multiplet.

At this point, we proceed by considering gluons as an example, and drop the index $`\kappa `$. For $`\mathrm{}`$ gluons in a given quantum state the trace of $`j_1,j_2\mathrm{}j_{\mathrm{}}|𝒪_\mu |j_2\mathrm{}j_{\mathrm{}},j_1`$ will be zero unless $`\mu `$ corresponds to $`\mathrm{}`$ identical gluons. The upper-left panel of Fig. 1 displays the projections $`\mu `$ that result from this trace applied for five gluons. The upper-right and lower-left panels display the $`(\mathrm{},\mathrm{})`$ and $`(\mathrm{}2,\mathrm{}+1)`$ multiplets. By careful inspection, one can see that the states along the diagonal, as shown in the lower-right panel of Fig. 1, are represented by the combination, $`(\mathrm{},\mathrm{})(\mathrm{}+1,\mathrm{}2)(\mathrm{}2,\mathrm{}+1)+(\mathrm{}1,\mathrm{}1)`$. By subtracting the same object evaluated with $`\mathrm{}^{}=\mathrm{}1`$ and adding two color singlets, one can see that

$$\chi _{(p,q)}^{\mathrm{gluons},\mathrm{}}=\{\begin{array}{cc}\hfill 1,& (p,q)=(\mathrm{},\mathrm{})\mathrm{or}(\mathrm{}3,\mathrm{})\mathrm{or}(\mathrm{},\mathrm{}3)\hfill \\ \hfill 1,& (p,q)=(\mathrm{}2,\mathrm{}+1)\mathrm{or}(\mathrm{}+1,\mathrm{}2)\mathrm{or}(\mathrm{}2,\mathrm{}2)\hfill \\ \hfill 2,& (p=0,q=0)\hfill \\ \hfill 0,& \mathrm{otherwise}\hfill \end{array}$$
(14)
Inserting this expression into Eq. (12) allows one to account for Bose effects in the gluonic partition function.

Quarks are represented by the multiplet $`(p=1,q=0)`$ while anti-quarks are represented by the multiplet $`(0,1)`$. Following the same ideas as were illustrated for the gluons, one can show

$$\chi _{(p,q)}^{\mathrm{quarks},\mathrm{}}=\{\begin{array}{cc}\hfill 1,& (p,q)=(\mathrm{},0)\mathrm{or}(\mathrm{}3,0)\hfill \\ \hfill 1,& (p,q)=(\mathrm{}2,1)\hfill \\ \hfill 0,& \mathrm{otherwise}\hfill \end{array}$$
(15)
By combining the recursion relations for additive charges, which are used to account for baryon and strangeness conservation, with the recursive techniques for SU(2), which are used to enforce conservation of isospin, and the methods for SU(3) which account for color, one can construct a recursive prescription which account for all the charges in a parton gas. One might also choose to write partition functions for subsystems, e.g., the strange quarks or the gluons, then convolute the partition functions together to find the partition function of the combined system. For instance, partitions for the two subsystems $`a`$ and $`b`$ can be combined to find the combined partition function,

$$Z_{(P,Q)}(T)=\underset{(p_a,q_a),(p_b,q_b)}{}Z_{(p_a,q_a)}^a(T)Z_{(p_b,q_b)}^b(T)𝒩((p_a,q_a),(p_b,q_b);(P,Q)).$$
(16)
## III Results and conclusions

The methods described in the last section were applied to the example of a parton gas in a volume $`V`$ at a temperature, $`T=250`$ MeV. The gluons, up-quarks and down-quarks were assumed to be massless, while the strange quark was assumed to have mass of 150 MeV. The single-particle partition functions were calculated by integrating over the momenta,

$$\begin{array}{cc}\hfill \omega (T)=& \frac{2V}{(2\pi )^3}d^3pe^{\beta \sqrt{p^2+m^2}}\hfill \\ \hfill =& \frac{mV}{2\beta ^2\pi ^2}\left(\beta mK_0(\beta m)+2K_1(\beta m)\right).\hfill \end{array}$$
(17)
The two helicities were accounted for by the factor of two preceding the expression. One could easily account for shell effects by replacing the integral over momentum with a discrete sum. Although discrete states are more consistent given the Bose and Fermi effects described above, the continuous form for $`\omega `$ is used here so that charge conservation effects can be viewed separately.

The calculations were performed according to the following prescription:

1. The gluon partition function, $`Z_{(p,q)}^{\mathrm{gluon}}`$, was calculated recursively by first calculating for all numbers of gluons $`A`$, then summing over $`A`$.
2. The partition function for strange quarks, $`Z_{A,(p,q)}^s`$, was generated. The partition function for anti-strange quarks, $`Z_{A,(p,q)}^{\overline{s}}`$, was then generated by switching $`(p,q)`$ with $`(q,p)`$. The partition function of the strange/anti-strange quark system was found by summing over all partition functions with equal numbers of $`s`$ and $`\overline{s}`$ quarks.

$$Z_{(p,q)}^{s\overline{s}}(T)=\underset{A,(p_a,q_q)(p_b,q_a)}{}Z_{A,(p_a,q_a)}^s(T)Z_{A,(p_b,q_b)}^{\overline{s}}(T)𝒩((p_a,q_a),(p_b,q_b);(p,q)).$$
(18)
3. The partition function for up and down quarks were calculated separately. In terms of the isospin projection, $`m`$, and the net number of up/down quarks, $`A`$, the partition function for the up and down quarks is:

$$Z_{A,m,(p,q)}^{ud}(T)=\underset{(p_a,q_a)(p_b,q_b)}{}Z_{m+A/2,(p_a,q_a)}^u(T)Z_{m+A/2,(p_b,q_b)}^d(T)𝒩((p_a,q_a),(p_b,q_b);(p,q)).$$
(19)
4. The partition function for $`\overline{u}\overline{d}`$ quarks is calculated by interchanging $`p`$ and $`q`$. In the same manner as the $`s\overline{s}`$ quarks were convoluted, the $`ud`$ and $`\overline{u}\overline{d}`$ partition functions were convoluted to find the partition function for $`u`$, $`d`$, $`\overline{u}`$ and $`\overline{d}`$ quarks with the constraint of zero baryon number, $`Z_{M,(p,q)}^{ud\overline{u}\overline{d}}`$, where $`M`$ is the isospin projection.
5. The partition function for the $`ud\overline{u}\overline{d}`$ system constrained to an isosinglet is calculated by taking the difference of the $`M=0`$ and $`M=1`$ partition functions.
6. Using Eq. (16), the partition function for the $`uds\overline{u}\overline{d}\overline{s}`$ system were generated. This was then convoluted with the partition function for gluons to find the partition function for the entire system.

Energy densities can be calculated with the well-known formula, $`E=(/\beta )\mathrm{ln}(Z)`$. This requires performing the calculation at two adjacent temperatures, doubling the CPU time of calculating the partition function. Thus, an alternative approach was developed where recursion relations were derived for $`(/\beta )Z`$. Recursion relations were also generated for the trace of $`Ne^{\beta H}`$, where $`N`$ could be one of several number operators, $`N_{\mathrm{gluons}}`$, $`N_s`$, $`N_u`$ or $`N_d`$. The recursion relations for these functions were calculated simultaneously with the recursion relations for the partition functions with little penalty in CPU time, because all calculations involve the same calculations of $`𝒩((p,q),(P^{},Q^{});(P,Q))`$.

Figure 2 shows the energy and entropy densities as a function of the volume. Both the entropy, $`S=\mathrm{ln}(Z)+\beta E`$, and the energy are lowered by the reduction in the available states due to the conservation constraints. The grand canonical limit, represented by the dashed lines, is approached at large volumes. As shown in elzeplb , restricting the matter to a color singlet reduces the partition function by a factor which scales as $`V^4`$. The other charge conservation constraints result in additional reductions to the quark sectors. A reduction of $`V^4`$ to the partition functions corresponds to a reduction of $`4\mathrm{ln}(V)`$ to the entropy and energy densities. Since the bulk contributions scale as $`V`$ the logarithmic conservation penalties become irrelevant at large volumes. The chemical decomposition of the various species behave similarly as shown in Fig. 3. The color penalty for gluons is more severe than the color penalties for quarks since gluons effectively have larger color charges and are less likely to couple to a singlet.

From Fig.s 2 and 3, it is clear that conservation rules are important for volumes of 5 fm3 or less. It is therefore important to take such considerations into account in $`pp`$ or $`pA`$ collisions. Heavy ion collisions at RHIC occupy thousands of cubic fm. However, the local nature of charge conservation results in an effective volume which is determined by the initial conditions, and the diffusion of the various charges. Given that the color of a projectile nucleon is spread over several units of rapidity by the initial stopping process, and given the spread of charge over a few fm in the transverse direction, it would be difficult to make a case that the effective volumes are less than a few dozen cubic fm. It should be emphasized that the determining factor for the importance of conservation constraints is the number of partons in the effective volume. For massless partons, the density of partons is expected to exceed 10 fm-3 in a during the first one or two fm/c of a RHIC collision. If the degrees of freedom are restricted by another means, e.g., a large effective mass for partons, the relative penalty for color conservation would increase.

###### Acknowledgements.

This work was supported by the National Science Foundation, Grant No. PHY-00-70818. J.R. acknowledges support by the Studienstiftung des Deutschen Volkes (German National Merit Foundation). The authors wish to thank D.D. Dietrich for sharing his profound understanding of SU(3).


text:
# Statistical hadronization of charm in heavy-ion collisions at SPS, RHIC and LHC

## 1 Introduction

The concept of a possible thermal origin of charmonia was initially introduced  to explain the measurement of the J/$`\psi `$/hadron ratio in nuclear collisions. Because of their large mass and small production cross section at thermal energies charm quarks are, however, not likely to be thermally produced. On the other hand, significant production of charm quark pairs takes place in initial hard collisions. This led to the idea of statistical hadronization of charm quarks , which has sparked an intense activity in this field . At about the same time an independent effort  based on a kinetic model has been developed to study charm production in heavy ion collisions at collider energies. Initial interest focussed on the available SPS data on $`J/\psi `$ production in Pb–Pb collisions. As we stress below, these data can indeed be described within the statistical approach, but only when assuming the charm production cross section to be enhanced beyond the perturbative QCD (pQCD) predictions.

The largest differences between the results of the statistical coalescence scenario and earlier models are expected at collider energies. For example, for central Au-Au collision at RHIC energy, the Satz–Matsui approach  predicts a very strong (up to a factor of 20 ) suppression of $`J/\psi `$ yields as compared to a direct production. In the present approach this suppression is overcome by statistical recombination of $`J/\psi `$ mesons from the initially produced $`c\overline{c}`$ pairs, so that much larger yields are expected.111 The absence of $`J/\psi `$ suppression at RHIC energy is also predicted by the kinetic model . We therefore focus in this letter on the predictions of the model for open and hidden charm meson yields at RHIC (Au–Au) and LHC (Pb–Pb) energies, with particular emphasis on the centrality dependence of rapidity densities at RHIC. We, furthermore, present predictions for the energy dependence of the production of hadrons with open and hidden charm.

## 2 Model description and input parameters

The model assumes that all charm quarks are produced in primary hard collisions and equilibrate222This implies thermal, but not chemical equilibrium for charm quarks. in the quark-gluon plasma (QGP). In particular, no $`J/\psi `$ mesons are preformed in the QGP, implying complete color screening . There are indications from recent lattice QCD calculations that, already at SPS energies, J/$`\psi `$ mesons can be dissociated in a deconfined medium via collisions with thermal gluons . As noted earlier , the cross sections of charmed hadrons are too small to allow for their chemical equilibration in a hadronic gas. We assume that open and hidden charm hadrons are formed at chemical freeze-out according to statistical laws. Consistent with the fact that, at SPS energy and beyond, chemical freeze-out appears to be at the phase boundary , our model implies that a QGP phase was a stage in the evolution of the fireball. The analysis of $`J/\psi `$ spectra at SPS  lends further support to the statistical hadronization picture where $`J/\psi `$ decouples at chemical freeze-out. A recent analysis of electron spectra at RHIC  also strenghtens the case for an early thermalization of heavy quarks. However, in that analysis it was pointed out that both the hydrodynamical approach and PYTHIA reproduce the measured single-electron spectra. As suggested meanwhile , the measurement of the elliptic flow ($`v_2`$) of open charm hadrons and $`J/\psi `$ may help to disentangle the two scenarios. Recent work  in extracting the quark $`v_2`$ from the one of hadrons may help towards this goal. Most important though, in our view, is a high-precision direct measurement of open charm. As expected, and as seen in ref. , the $`p_t`$ spectra predicted by PYTHIA and hydro are different in detail at low $`p_t`$ and differ manifestly at high $`p_t`$ ($`p_t`$ mass of charm quark). On the other hand, a recent theoretical analysis  indicates that charm quarks might not thermalize quickly because of their large mass.

In statistical models charm production needs generally to be treated within the framework of canonical thermodynamics . Thus, the charm balance equation required  during hadronization is expressed as:

$$N_{c\overline{c}}^{dir}=\frac{1}{2}g_cN_{oc}^{th}\frac{I_1(g_cN_{oc}^{th})}{I_0(g_cN_{oc}^{th})}+g_c^2N_{c\overline{c}}^{th}.$$
(1)
Here $`N_{c\overline{c}}^{dir}`$ is the number of directly produced $`c\overline{c}`$ pairs and $`I_n`$ are modified Bessel functions. In the fireball of volume $`V`$ the total number of open $`N_{oc}^{th}=n_{oc}^{th}V`$ and hidden $`N_{c\overline{c}}^{th}=n_{c\overline{c}}^{th}V`$ charm hadrons are computed from their grand-canonical densities $`n_{oc}^{th}`$ and $`n_{c\overline{c}}^{th}`$, respectively. The densities of different particle species in the grand canonical ensemble are calculated following the statistical model introduced in . All known charmed mesons and hyperons and their decays  have to be included in the calculations.

The balance equation (1) defines a fugacity parameter $`g_c`$ that accounts for deviations of charm multiplicity from the value that is expected in complete chemical equilibrium. The yield of open charm mesons and hyperons $`i`$ and of charmonia $`j`$ is obtained from:

$$N_i=g_cN_i^{th}\frac{I_1(g_cN_{oc}^{th})}{I_0(g_cN_{oc}^{th})}\mathrm{and}N_j=g_c^2N_j^{th}.$$
(2)
The above model for charm production and hadronization can be only used  if the number of participating nucleons $`N_{part}`$ is sufficiently large. Taking into account the measured dependence of the relative yield of $`\psi ^{}`$ to $`J/\psi `$ on centrality  in Pb–Pb collisions at SPS energy the model appears appropriate  if $`N_{part}>`$100.

To calculate the yields of open and hidden charm hadrons for a given centrality and collision energy one needs to fix a set of parameters in Eq.(1) and (2):

i) A constant temperature of 170 MeV and a baryonic chemical potential $`\mu _b`$ according to the parametrization $`\mu _b=1270[\mathrm{MeV}]/(1+\sqrt{s_{NN}}[\mathrm{GeV}]/4.3)`$  are used for our calculations. For SPS and RHIC energies these thermal parameters are consistent, within statistical errors, with those required to describe experimental data on different hadron yields . In Table 1 we present the values of these input parameters and the resulting chemical potentials for isospin, strangeness and charmness for the canonical ensemble calculations at SPS, RHIC and LHC.

ii) The volume of the fireball corresponding to a slice of one rapidity unit at midrapidity, $`V_{\mathrm{\Delta }y=1}`$, is obtained from the charged particle rapidity density $`\mathrm{d}N_{ch}/\mathrm{d}y`$, via the relation $`\mathrm{d}N_{ch}/\mathrm{d}y=n_{ch}^{th}V_{\mathrm{\Delta }y=1}`$, where $`n_{ch}^{th}`$ is the charged particle density computed within the thermal model. The charged particle rapidity densities (and total yields in case of SPS) are taken from experiments at SPS  and RHIC  and extrapolated to LHC energy. In Fig. 1 we show a compilation of experimental data on $`\mathrm{d}N_{ch}/\mathrm{d}y`$ at midrapidity for central collisions in the energy range from AGS  up to RHIC333A constant Jacobian of 1.1 has been used to convert the $`\mathrm{d}N_{ch}/\mathrm{d}\eta `$ data to $`\mathrm{d}N_{ch}/\mathrm{d}y`$.. The continuous line in Fig. 1 is the prediction  of the saturation model, the dashed line is a flatter power law dependence, both arbitrarily normalized. We use the $`(\sqrt{s_{NN}})^{0.3}`$ dependence of $`\mathrm{d}N_{ch}/\mathrm{d}y`$ to get the value of $`\mathrm{d}N_{ch}/\mathrm{d}y`$2000 for LHC energy. Throughout this paper we define central collisions by the average value $`N_{part}`$=350, which roughly corresponds to the topmost 5% of the nuclear cross section. For the centrality dependences we assume that the volume of the fireball is proportional to $`N_{part}`$.

iii) The yield of open charm $`N_{c\overline{c}}^{dir}`$ is taken from next-to-leading order (NLO) pQCD calculations for p–p collisions  and scaled to nucleus–nucleus collision via the nuclear overlap function  $`T_{AA}`$. For a given centrality: $`N_{c\overline{c}}^{dir}(N_{part})=\sigma (ppc\overline{c})T_{AA}(N_{part})`$. The pQCD calculations with the MRST HO parton distribution function (PDF) are used here . We note that the pQCD results, in particular at LHC energy, are sensitive on the choice of PDF and/or charm quark mass value .

The input values $`\mathrm{d}N_{ch}/\mathrm{d}y`$ and $`\mathrm{d}N_{c\overline{c}}^{dir}/\mathrm{d}y`$ and the corresponding volume at midrapidity and enhancement factor are summarized in Table 2 for model calculations for different collision energies.

## 3 Model results and predictions

We first compare predictions of the model to 4$`\pi `$-integrated $`J/\psi `$ data at the SPS energy. The NA50 data  are replotted in Fig. 2 as outlined in . The model results shown in Fig. 2 are going beyond that of  since we have included a complete set of charmed mesons and baryons as well as updated the values of the charm production cross section and the volume of the fireball. Recent results  of the NA49 Collaboration on total charged particle multiplicity in central collisions yield $`N_{ch}=1533`$, leading to a fireball total volume V=3070 fm3 and to a corresponding total yield of thermal open charm pairs of $`N_{oc}^{th}`$=0.98. This is to be contrasted with $`N_{c\overline{c}}^{dir}`$=0.137 from NLO calculations , leading to a value of $`g_c`$=0.78. Although $`g_c`$ is here close to unity, this obviously does not indicate that charm production appears at chemical equilibrium, as the suppression factor is a strongly varying function of the collision energy. We have already indicated that, within the time scales available in heavy ion collisions, the chemical equilibration of charm is very unlikely both in confined and deconfined media .

In Fig. 2 we show the comparison between the results of our model and NA50 data for two different values of $`N_{c\overline{c}}^{dir}`$: from NLO calculations  and scaled up by a factor of 2.8. Using the NLO cross sections for charm production scaled by the nuclear overlap function, the model understimates the measured yield. To explain the overall magnitude of the data, one needs to increase the $`N_{c\overline{c}}^{dir}`$ yield by a factor of 2.8 as compared to NLO calculations. We mention in this context that the observed  enhancement of the di-muon yield at intermediate masses has been interpreted as a possible indication for an anomalous increase of the charm production cross section . A third calculation (resulting in the dash-dotted line in Fig. 2) is using the NLO cross section scaled-up by 1.6, which is the ratio of the open charm cross section estimated by NA50 for p–p collisions at 450 GeV/c  and the present NLO values. For this case the $`N_{part}`$ scaling is not the overlap function, but is taken according to the measured di-muon enhancement as a function of $`N_{part}`$ . The resulting $`J/\psi `$ yields from the statistical model are on average in agreement with the data, albeit with a flatter centrality dependence than by using the nuclear overlap function. Thus our charm enhancement factor of 2.8 needed to explain the $`J/\psi `$ data is very similar to the factor needed to explain the intermediate mass dilepton enhancement assuming that it arises exclusively from charm enhancement . We note, however, that other plausible explanations exist of the observed enhancement in terms of thermal radiation .

The drop of the $`J/\psi `$ yield per participant observed in the data for $`N_{part}>`$350 (see Fig. 2) is currently understood in terms of further $`J/\psi `$ dissociation due to energy density fluctuations for a given overlap geometry . Alternatively, the drop has been connected with a possible trigger bias in the data . Note that this second drop is only visible in the minimum bias analysis . These aspects are not included in the current model, however when being incorporated the $`J/\psi `$ yield for very central collisions can be as well reproduced .

We turn now to discuss our model predictions for charmonia and open charm production at collider energies and compare them with the results obtained at SPS. Notice that from now on we focus on rapidity densities, which are the relevant observables at the colliders. In Table 3 we summarize the yields for a selection of hadrons with open and hidden charm. All predicted yields increase strongly with beam energy, reflecting the increased charm cross section and the concomitant importance of statistical recombination. Also, ratios of open charm hadrons evolve with increasing energy, reflecting the corresponding decrease in the charm chemical potential.

Model predictions for the centrality dependence of $`J/\psi `$ and D+ rapidity densities normalized to $`N_{part}`$ are shown in Fig. 3. The results for $`J/\psi `$ mesons exhibit, in addition to the dramatic change in magnitude, a striking change in the shape of the centrality dependence. In terms of the model this change is a consequence of the transition from the canonical to the grand-canonical regime . For D+-mesons, the expected approximate scaling of the ratio D+/$`N_{part}`$ $`N_{part}^{1/3}`$ (dashed lines in Fig. 3) is only roughly fulfilled due to departures of the nuclear overlap function from the simple $`N_{part}^{4/3}`$ dependence.

The results summarized in Table 3 and shown in Fig. 3 obviously depend on two input parameters, $`\mathrm{d}N_{ch}/\mathrm{d}y`$ and $`\mathrm{d}N_{c\overline{c}}^{dir}/\mathrm{d}y`$. For LHC energy, neither one of these parameters is well known. An increase of charged particle multiplicities by up to a factor of three beyond our “nominal” value $`\mathrm{d}N_{ch}/\mathrm{d}y`$=2000 for central collisions is conceivable. The saturation model, e.g., leads to a prediction of  $`\mathrm{d}N_{ch}/\mathrm{d}y`$2300. However, due to quite large uncertainties on the amount of shadowing at LHC energy, these results may be still modified. The yield of $`\mathrm{d}N_{c\overline{c}}^{dir}/\mathrm{d}y`$ is also not well known at LHC energy. Although these uncertainties affect considerably the magnitude of the predicted yields, their centrality dependence remains qualitatively unchanged: the yields per participant are increasing functions of $`N_{part}`$. We also note here that, while detailed predictions differ significantly, qualitatively similar results (see ref. ) have been obtained for a kinetic model study of $`J/\psi `$ production at the LHC.

In Fig. 4 we present the predicted centrality dependence of the $`J/\psi `$ rapidity density normalized to $`N_{part}`$ for RHIC energy ($`\sqrt{s_{NN}}`$=200 GeV). The three panels show its sensitivity on $`\mathrm{d}N_{ch}/\mathrm{d}y`$, $`\mathrm{d}N_{c\overline{c}}^{dir}/\mathrm{d}y`$, and (freeze-out) temperature $`T`$. The calculations are compared to experimental results of the PHENIX Collaboration . Note that the point for the central collisions is the upper limit extracted by PHENIX for 90% C.L. . The experimental data have been rescaled according to our procedure  to calculate $`N_{part}`$ and the number of binary collisions, $`N_{coll}`$. Using the centrality intervals of PHENIX, we obtain $`N_{coll}`$ = 275 and 782 (the differences to the PHENIX values of 296 and 779, respectively, are small). We use $`\sigma _{NN}`$=42 mb and a Woods-Saxon density profile.

Within the still large experimental error bars, the measurements agree with our model predictions. In Fig. 4 only the statistical errors of the mid-central data point are plotted. The systematic errors are also large . peripheral collisions. A stringent test of the present model can only be made when high statistics $`J/\psi `$ data are available. In any case, very large $`J/\psi `$ suppression factors as predicted, e.g., by  seem to be not supported by the data.

We turn now to a more detailed discussion of the sensitivity of our calculations to the various input parameters as quantified in Fig. 4. First we consider the influence of a $`10\%`$ variation of $`\mathrm{d}N_{ch}/\mathrm{d}y`$ on the centrality dependence of $`J/\psi `$ yield. Note that the total experimental uncertainty of $`\mathrm{d}N_{ch}/\mathrm{d}\eta `$ (which is for the moment the measured observable for most experiments) at RHIC is below 10% . The sensitivity on the $`\mathrm{d}N_{ch}/\mathrm{d}y`$ values stems from the volume into which the (fixed) initial number of charm quarks is distributed. The smaller the particle multiplicities and thus also the fireball volume, the more probable it is for charm quarks and antiquarks to combine and form quarkonia. That is why one sees, in the top panel in Fig. 4, that the $`J/\psi `$ yield is increasing with decreasing charge particle multiplicity.

The sensitivity of the predicted $`J/\psi `$ yields on $`\mathrm{d}N_{c\overline{c}}^{dir}/\mathrm{d}y`$ is also straightforward. The larger this number is in a fixed volume the larger is the yield of charmed hadrons. In case of charmonia the dependence on $`\mathrm{d}N_{c\overline{c}}^{dir}/\mathrm{d}y`$ is non-linear due to their double charm quark content, as reflected by the factor $`g_c^2`$ in equation (2). To illustrate the sensitivity of the model predictions on $`\mathrm{d}N_{c\overline{c}}^{dir}/\mathrm{d}y`$, we exhibit the results of a 20% variation with respect to the value given in Table 2. The open charm cross section is not yet measured at RHIC. However, some indirect measurements can be well reproduced, within the experimental errors, by PYTHIA calculations using a p–p charm cross sections scaled with the number of collisions $`N_{coll}`$ of 650 $`\mu `$b . The corresponding value at $`\sqrt{s_{NN}}`$=130 GeV is 330 $`\mu `$b . For comparison, the NLO pQCD values we are using  are 390 and 235 $`\mu `$b, respectively. Despite the still large experimental uncertainties, this discrepancy, recognized earlier in ref. , needs to be understood. We note that, dependent on the input parameters used in the NLO calculations , possible variations of the open charm production cross section for the RHIC energy are of the order of $`\pm `$20%. In terms of our model this variation corresponds to about a $`\pm `$30% change in $`J/\psi `$ yield which is also centrality dependent (see middle panel in Fig. 4). If we use the PHENIX p–p cross section of 650 $`\mu `$b, the calculated yield is a factor 2.5 larger for $`N_{part}`$=350 and increases somewhat stronger with centrality, in disagreement with the data. As apparent in Fig. 4, the predictive power of our model or of any similar model relies heavily on the accurate knowledge of the overall charm production cross section. A simultaneous description of the centrality dependence of open charm together with $`J/\psi `$ production is, in this respect, mandatory to test the concept of the statistical origin of open and hidden charm hadrons in heavy ion collisions at relativistic energies.

The apparent weak dependence of $`J/\psi `$ yield on freeze-out temperature, seen in Fig. 4, may be surprising. This is particularly the case since in the thermal model considered in  this dependence was shown to be very strong. In our model this result is a consequence of the charm balance equation (2). The temperature variation leads, obviously, to a different number of thermally produced charmed hadrons as in , but this is compensated by the $`g_c`$ factor. The approximate temperature dependence of $`g_c`$ and the $`J/\psi `$ yield are:

$$g_c(T)1/N_D^{th}e^{\frac{m_D}{T}},N_{J/\psi }(T)=g_c^2N_{J/\psi }^{th}e^{\frac{2m_Dm_{J/\psi }}{T}}$$
(3)
As a result of the small mass difference in the exponent the $`J/\psi `$ yield exhibits only a weak sensitivity on $`T`$. This is in contrast to the purely thermal case where the yield scales with $`\mathrm{exp}(m_{J/\psi }/T)`$. The results are summarized in Table 4, where we present the number of $`J/\psi `$ and $`D^+`$ mesons produced thermally ($`N_{J/\psi }^{th}`$ and $`N_{D^+}^{th}`$, respectively) and within the frame of our statistical hadronization model. Also shown are the values of $`g_c`$ and the yield ratio $`\psi ^{}/J/\psi `$ for three different values of $`T`$. The compensations are evident in the last column of Table 4, where ratios for the two extreme temperatures are calculated. The only exception is the ratio $`\psi ^{}/J/\psi `$, which is obviously identical in the statistical hadronization scenario and in the thermal model and coincides, for T$``$170 MeV, with the measured value at SPS .

Most of our results presented above are obtained considering a one unit rapidity window at midrapidity, while results for the full volume were presented only for the SPS. Unlike the kinetic model of Thews et al. , our model does not contain dynamical aspects of the coalescence process. However, in our approach, the width of the rapidity window does influence the results in the canonical regime. For the grand-canonical case, attained only at LHC energy, there is no dependence on the width of the rapidity window, due to a simple cancellation between the variation of the volume, proportional to the rapidity slice in case of a flat rapidity distribution, and the variation of $`N_{c\overline{c}}^{dir}`$, also proportional to the width of the rapidity slice. In Fig. 5 we present the centrality dependence of $`J/\psi `$ rapidity densities for RHIC energy and for different rapidity windows $`\mathrm{\Delta }`$y from 0.5 to 3. The dependence on $`\mathrm{\Delta }`$y resembles that of the kinetic model , but is less pronounced. The available data are not yet precise enough to rule out any of the scenarios considered. However, for the kinetic model , the cases of small $`\mathrm{\Delta }`$y seem to be ruled out by the present PHENIX data.

We stress in this context that the size of $`\mathrm{\Delta }`$y window has a potentially large impact on the results at SPS energy. It is conceivable that no charm enhancement is needed to explain the data if one considers a sufficiently narrow rapidity window for the statistical hadronization.

The results presented above were obtained under the assumption of statistical hadronization of quarks and gluons. In addition, for charm quarks the yields were constrained by the charm balance equation that was formulated in the canonical ensemble. We have further assumed that charm quarks are entirely produced via primary hard scattering and thermalized in the QGP. No secondary production of charm in the initial and final state was included in our calculations. Final state effects like nuclear absorption of $`J/\psi `$  are also neglected.

## 4 Conclusions

We have demonstrated that the statistical coalescence approach yields a good description of the measured centrality dependence of $`J/\psi `$ production at SPS energy, albeit with a charm cross section increased by a factor of 2.8 compared to current NLO perturbative QCD calculations. Rapidity densities for open and hidden charm mesons are predicted to increase strongly with energy, with striking changes in centrality dependence. First RHIC data on $`J/\psi `$ production support the current predictions, although the experimental errors are for the moment too large to allow firm conclusions. We emphasize that the predictive power of our model relies heavily on the accuracy of the charm cross section, not yet measured in heavy-ion collisions. The statistical coalescence implies travel of charm quarks over significant distances e.g. in a QGP. If the model predictions will describe consistently precision data this would be a clear signal for the presence of a deconfined phase.

## Acknowledgments

We acknowledge stimulating discussions with R.L. Thews. K.R. acknowledges the support of the Alexander von Humboldt Foundation (AvH).