text:
# Epitaxial growth with pulsed deposition: Submonolayer scaling and Villain instability

## I Introduction

Pulsed Laser Depostion (PLD) is an increasingly used growth method which is particularly suited for the fabrication of multilayer thin films ChriseyHubler . The basic mechanism of this technique is that the target material is ablated by a pulsed laser and then deposited in pulses on a substrate so that many atoms arrive at the surface simultaneously. Experimentally each pulse lasts for about a few nanoseconds, and the time between two pulses is of the order of seconds. PLD differs significantly from ordinary molecular beam epitaxy (MBE), where atoms are deposited continuously. Growth experiments show that the surfaces in both cases have different morphologies. In particular it was observed that in some cases PLD leads to smoother surfaces than MBE Jenniches96 . However, as there are only few theoretical studies so far, the mechanism leading to smoother surfaces in PLD is not yet understood so that it is impossible to predict the growth quality for various materials and growth parameters. Therefore the aim of the present work is to improve our understanding of PLD by investigating the growth morphology for different growth parameters in a simple model.

Without going into experimental details it should be mentioned that the physical conditions for PLD are far less well defined than for MBE. The particles deposited may be atoms, clusters or even droplets. They may arrive with energies ranging from $`0.1`$ eV to 1000 eV. If very energetic particles arrive at the surface their kinetic energy is converted into heat at the surface, changing locally the effective mobility of the particles for a short time. In the present work this transient mobility is neglected, i.e., it is assumed that the particles arrive at the substrate surface with energies of about $`0.1`$ eV which are of the same order of magnitude as in the case of thermal deposition.

Recently Narhe et al. investigated the island statistics in the coalescence regime for a PLD process where tin droplets are deposited on a sapphire substrate Narhe . It was observed that the scaling differs significantly from MBE due to the large fraction of multiple droplet coalescence under pulsed vapor delivery. While these results are valid for high deposition energies, the present work investigates a different physical regime, namely pulsed deposition of individual atoms at low energies and the formation of islands in two dimensions.

The model investigated in the present work is defined as a solid-on-solid growth process on a square lattice of $`L\times L`$ sites with integer heights representing the configuration of the adsorbed layer. The particles are deposited in pulses with an intensity $`I`$ which is defined as the number of particles per unit area deposited per pulse. The duration of a pulse is assumed to be zero and the transient enhancement of the mobility of freshly deposited adatoms is neglected. The model is controlled by three parameters, namely, the intensity $`I`$ of the pulses, the diffusion constant $`D`$, and the average flux density of incoming particles $`F`$. One of these parameters can be fixed by choosing the time scale so that we can use $`I`$ and $`D/F`$ as independent parameters. The dynamic rules are defined as follows. (i) In each pulse $`IL^2`$ atoms are instantaneously deposited at random positions on the surface. (ii) Between two pulses a time interval $`\mathrm{\Delta }t=I/F`$ elapses, in which adatoms diffuse on the surface with rate $`D`$. If Ehrlich-Schwoebel barriers are present the rate at which particles hop down the edge of an island is reduced. (iii) If two atoms at the same height occupy neighboring sites they stick together irreversibly, forming the nucleus of a new island or attaching to an already existing island. In our model the influence of edge diffusion is neglected so that islands grow in a fractal manner before they coalesce. We present results from kinetic Monte Carlo simulations on a lattice of $`400\times 400`$ sites. Finite size effects are discussed in the Appendix.

The paper is organized as follows. In the following section we first recall the properties of PLD without Ehrlich-Schwoebel barriers and its crossover to MBE in the limit of very low intensities JensenNiemeyer ; ConferencePaper . In Sec. III we investigate the scaling behavior of the time-dependent nucleation density in PLD, extending our previous analysis in Ref. OurPRL . In Sec. IV we study the influence of Ehrlich-Schwoebel barriers on the morphology of surfaces grown by PLD, finding a parameter range where layerwise growth is improved compared to MBE.

## II From MBE to PLD: <br>Crossover and Scaling

If the intensity $`I`$ is very small PLD behaves essentially in the same way as MBE. In fact, at the lowest possible intensity, where only one atom per pulse is deposited, MBE and PLD are equivalent up to minor statistical differences resulting from finite size effects (see Appendix). Using lattices of $`400^2`$ sites we therefore restrict our analysis to the range $`D/F10^7`$ where these finite size effects are negligible.

At high intensities the growth morphology of PLD differs significantly from MBE. As shown in Fig. 1, there are much more nucleations at an early stage, although the effective flux of incoming particles is the same in both cases. The two regimes are separated by a crossover at a certain intensity $`I_c`$, where the number of deposited atoms per pulse is of the same order of magnitude as the average adatom density in the corresponding MBE process. Obviously, if the pulse intensity is much higher than the MBE adatom density, the adatoms nucleate much faster forming many small islands. For this reason PLD is expected to yield more homogeneous surfaces, which explains the technological interest in this method.

Let us now study the crossover from MBE to PLD in more detail JensenNiemeyer ; ConferencePaper . The adatom density in MBE averaged over space and time is known to scale as $`(D/F)^{2\gamma 1}`$ VPW92 , where

$$\gamma =\frac{1}{d_f+4}=\{\begin{array}{cc}1/6\hfill & \text{for compact islands}\hfill \\ 0.18\hfill & \text{for DLA fractals }\text{Meakin88}\hfill \end{array}.$$
(1)
Thus the crossover takes place at the critical intensity

$$I_c(D/F)^{2\gamma 1}.$$
(2)
A quantity which distinguishes the two different growth modes shown in Fig. 1 is the average island distance. Performing numerical simulations we realized that the scaling regime of this quantity is not only restricted by finite-size effects but also by lattice effects, which become relevant for high intensities as well as for low values of $`D/F`$, where the average distance of islands is of the order of a few lattice constants. Combining these bounds we find that for a system with $`400^2`$ sites finite-size and lattice effects are negligible in the parameter range $`10^4D/F10^7`$ and $`I10^2`$.

Fig. 2 shows the average island distance as a function of the intensity for various values of $`D/F`$ at a fixed coverage of $`0.2`$ monolayers (ML). For low intensities $`I<I_c`$ the island distance depends only on $`D/F`$. Plotting these saturation values versus $`D/F`$ (shown in the inset of the figure) one recovers the well-known power law for the island distance in MBE

$$\mathrm{}_D(D/F)^\gamma $$
(3)
with $`\gamma =0.17(1)`$. This estimate lies between the values for compact growth and diffusion-limited aggregation (see Eq. (1)), supporting the assumption that the islands are characterized by an effective fractal dimension $`1.6<d_f<2`$.

For intensities $`I>I_c`$ the island distance is independent of $`D/F`$ and can be described by a power law

$$\mathrm{}_DI^\nu $$
(4)
with an exponent $`\nu =0.26(1)`$. The independence of the island distance of $`D/F`$ is a result of the high density of adatoms which nucleate so quickly that they do not make use of their full mobility given by $`D/F`$.

Combining Eqs. (3) and (4) we obtain the scaling form ConferencePaper

$$\mathrm{}_D(D/F)^\gamma h(I/I_c),$$
(5)
where $`h`$ is a scaling function with the asymptotic behavior

$$h(y)\{\begin{array}{cc}\text{const.}\hfill & \text{for }y1\hfill \\ y^\nu \hfill & \text{for }y1.\hfill \end{array}$$
(6)
As shown in Fig. 3, this scaling form leads to a convincing data collapse. Moreover, the independence of $`D/F`$ at high intensities together with Eq. (2) implies that

$$\nu =\frac{\gamma }{12\gamma }\{\begin{array}{cc}0.25\hfill & \text{for compact islands}\hfill \\ 0.28\hfill & \text{for DLA fractals}\hfill \end{array}.$$
(7)
confirming our numerical estimate $`\nu =0.26(1)`$.

## III Time-dependent nucleation density

In a recent paper OurPRL  we investigated the time-dependent nucleation density in PLD, reporting an unusual type of scaling behavior. In order to avoid the crossover to MBE, we considered the limit $`D/F\mathrm{}`$ where $`I_c=0`$. In this limit all adatoms nucleate or attach to existing islands before the next pulse arrives. Fig. 4 shows a log-log plot of the the nucleation density at the bottom layer $`n(I,\mathrm{\Theta })`$ as a function of the coverage $`\mathrm{\Theta }=Ft`$ for various intensities $`I`$. Obviously this quantity does not display ordinary power law scaling since it is impossible to collapse the curves by shifting them horizontally and vertically. However, in Ref. OurPRL  we observed that the normalized nucleation density

$$M(I,\mathrm{\Theta })=n(I,\mathrm{\Theta })/n(I,1)$$
(8)
obeys an unusual logarithmic scaling law of the form

$$\mathrm{ln}M(I,\mathrm{\Theta })\left(\mathrm{ln}I\right)g(\mathrm{ln}\mathrm{\Theta }/\mathrm{ln}I).$$
(9)
This scaling form was also proposed by Kadanoff et al. and Tang in the context of multiscaling in self-organized criticality Kadanoff88 ; Tang89 . As shown in Fig. 5 this scaling form leads to a convincing data collapse. More recently the same type of scaling has also been observed in one-dimensional systems, so that we can rule out logarithmic corrections at the marginal dimension of random walks as the origin for this type of non-conventional scaling Lee ; Krapivsky .

In Fig. 5 the left terminal points of the curves at $`g(1)=0.44(2)`$ can be used to determine $`\gamma `$ as follows. It is known that the nucleation density after the first pulse grows linearly with the intensity OurPRL  while the nucleation density after completion of one monolayer grows as $`n(I,1)I^{2\nu }`$ (see Eq. (4)). Therefore, the normalized nucleation density after the first pulse scales as $`M(I,I)I^{12\nu }`$ so that

$$g(1)=12\nu =\frac{14\gamma }{12\gamma }.$$
(10)
Solving this equation we obtain $`\gamma =0.179(4)`$ in agreement with the estimate in Fig. 2.

In order to understand how this unusual scaling behavior for PLD crosses over to the ordinary power-law scaling of MBE let us now turn to the case of finite $`D/F`$, where the system is characterized by a typical length scale $`\mathrm{}_0(D/F)^{1/4}`$. Generalizing the results of Ref. OurPRL  we consider islands with an arbitrary fractal dimension, i.e., we regard $`\gamma `$ as a free parameter.

In MBE the nucleation density is known to exhibit ordinary power law scaling of the form

$$n(\mathrm{}_0,\mathrm{\Theta })=\mathrm{}_0^2f(\mathrm{\Theta }\mathrm{}_0^2),$$
(11)
where $`f`$ is a scaling function with the asymptotic behavior Tang93

$$f(z)\{\begin{array}{cc}z^3\hfill & \text{ for }0z1\hfill \\ z^{1/3}\hfill & \text{ for }1zz_{\mathrm{max}}.\hfill \end{array}$$
(12)
The upper bound $`z_{\mathrm{max}}`$ is determined by the condition that the whole surface is covered by islands so that no further nucleation in the respective layer is possible VPW92 . Because of $`I_c(D/F)^{2\gamma 1}`$ the length scale $`\mathrm{}_0`$ is related to the crossover intensity by

$$\mathrm{}_0I_c^{1/(48\gamma )}.$$
(13)
Using this scaling form we obtain the expression

$$M(I_c,\mathrm{\Theta })=\{\begin{array}{cc}\mathrm{\Theta }^3I_c^{\frac{31/3}{24\gamma }}\hfill & \text{ for }0z1\hfill \\ \mathrm{\Theta }^{1/3}\hfill & \text{ for }1zz_{\mathrm{max}},\hfill \end{array}$$
(14)
where $`z=\mathrm{\Theta }I_c^{1/(24\gamma )}`$. Taking the logarithm and extrapolating to $`z=\mathrm{}_0^2`$ we arrive at

$$\frac{\mathrm{ln}M(I_c,\mathrm{\Theta })}{\mathrm{ln}I_c}=\{\begin{array}{cc}3\frac{\mathrm{ln}\mathrm{\Theta }}{\mathrm{ln}I_c}+\frac{4}{6\gamma 3}\hfill & \text{ for }\frac{1}{24\gamma }\frac{\mathrm{ln}\mathrm{\Theta }}{\mathrm{ln}I_c}\hfill \\ \frac{\mathrm{ln}\mathrm{\Theta }}{3\mathrm{ln}I_c}\hfill & \text{ for }0\frac{\mathrm{ln}\mathrm{\Theta }}{\mathrm{ln}I_c}\frac{1}{24\gamma }.\hfill \end{array}$$
(15)
In the limit $`II_c0`$ (i.e., $`D/F\mathrm{})`$ the crossover between both regimes becomes sharper and converges to a piecewise linear curve, which is shown in Fig. 5 as a dashed line. The crossover point is located at

$`\mathrm{ln}\mathrm{\Theta }_c/\mathrm{ln}I_c`$ $`=`$ $`1/(24\gamma ),`$
$`\mathrm{ln}M(I_c,\mathrm{\Theta })/\mathrm{ln}I_c`$ $`=`$ $`1/(612\gamma ).`$ (16)
Surprisingly, the crossover point lies on the collapsed curves for PLD within numerical errors. This is plausible for the following reasons. On the one hand the PLD-curve must be an upper bound for the MBE-curve since the island density in PLD is always larger. On the other hand, if the gap between the two extrapolated curves did not close at the crossover between PLD- and MBE-behavior, it would imply that there is an additional characteristic length in the system, for which we have no evidence.

Furthermore, we note that the scaling function itself roughly follows a power law

$$g(z)g(1)z^\beta .$$
(17)
As shown in the double-logarithmic inset of Fig. 5, the effective exponent $`\beta `$ varies between 2.4 and 2.0. In order to verify this estimate we derive the exponent $`\beta `$ from Eq. (16) by assuming that the crossover point lies exactly on the collapsed curves for PLD. This leads to an expression for $`\beta `$ in terms of the exponent $`\gamma `$, namely

$$\beta =\frac{\mathrm{ln}(624\gamma )}{\mathrm{ln}(24\gamma )}\{\begin{array}{cc}2.4\hfill & \text{ for compact islands}\hfill \\ 2.15\hfill & \text{ for DLA fractals }\hfill \end{array}.$$
(18)
Since this is just the range in which the numerical values for $`\beta `$ vary, we are led to the conclusion that the effective fractal dimensions of the islands for low and high coverages are different. For low coverages the islands are spaced relatively far apart so that the growth is DLA-like while for high coverages the islands coalesce and become more and more compact.

## IV Influence of Ehrlich-Schwoebel barriers

In most experimental situations interlayer transport is reduced by an additional energy barrier $`E_{\mathrm{ES}}`$ which the atoms have to overcome when hopping down an edge of an island Schwoebel . This barrier is called Ehrlich-Schwoebel (ES) barrier and is of the order of about $`0.1`$ eV for metals SmilauerSchwoebel . A useful measure for the influence of this barrier is the Schwoebel length $`\mathrm{}_{\mathrm{ES}}`$, which is defined as KrugSchwoebel

$$\mathrm{}_{\mathrm{ES}}=a\text{exp}\left[\frac{E_{\mathrm{ES}}}{k_\mathrm{B}T}\right],$$
(19)
where $`a=1`$ denotes the lattice constant. For MBE it is known that Ehrlich-Schwoebel barriers impede interlayer transport. This leads to a growth instability which was predicted and first theoretically investigated by Villain Villain91  followed by others Siegert  and has also been observed experimentally SchwoebelObservations . In contrast to MBE, not much is known about the influence of Ehrlich-Schwoebel barriers on PLD. However, it has been observed experimentally that PLD leads to better growth results than MBE in certain situations where Ehrlich-Schwoebel barriers are present Jenniches96 .

The influence of ES barriers on PLD is twofold. On the one hand, for high intensities many adatoms are deposited on the same island, which should increase the influence of ES barriers. On the other hand, in PLD the islands are much smaller so that adatoms tend to leave an island very quickly, thereby reducing the influence of ES barriers. Therefore, the question arises whether the experimentally observed improved quality of layer-by-layer growth in PLD can be related to a reduced influence of ES barriers and whether it is possible to choose the parameters of the model in such a way that PLD produces smoother surfaces than MBE.

This question has been addressed previously by Schinzer et al. Biehl . They studied the special case where one intense pulse (0.23 ML) is deposited at the beginning of each monolayer, while the remaining atoms are deposited with a continuous flux. The aim was to check an idea of Rosenfeld et al. Rosenfeld , that it should improve the smoothness of the surface, if the islands are forced to be much smaller than the diffusion length of the atoms deposited after the pulse. Surprisingly, Schinzer et al. Biehl  found that this is not the case in their simulation. On the contrary, in the presence of the pulse the surface became rougher.

A key concept for describing the influence of ES-barriers is the time it takes before the first nucleation in the second layer takes place. Layer-by-layer growth requires that nucleations in the second layer do not start significantly before the first layer is completed Villain91 . In the case of MBE the study of the second-layer nucleation time turned out to be very useful in order to predict the growth mode for a given set of parameters Tersoff94 ; Krug00 . In the following we demonstrate that this concept can be successfully applied to PLD as well. To this end we first investigate the growth behavior of PLD in the limit of $`D/F\mathrm{}`$. It will be shown that for any finite Ehrlich-Schwoebel barrier MBE produces a better layer-by-layer growth than PLD in agreement with the finding of Schinzer et al.Biehl . Then we compare MBE and PLD for different values of $`I`$ and $`D/F`$ over a wide range of Ehrlich-Schwoebel barriers. It turns out that there is indeed a regime where PLD produces better growth results than MBE.

### IV.1 Infinite $`D/F`$

We start by explaining why in the limit $`D/F\mathrm{}`$ PLD does not delay the Villain instability in spite of the island size reduction compared to MBE, provided that the ES-barrier is sufficiently high and not leaky e.g. at kink sites. In this limit the time scales are separated as follows. Because of the high ES barrier the nucleation time of adatoms on top of islands is much smaller than their residence time, while for $`D/F\mathrm{}`$ the residence time is in turn much smaller than the time interval between two pulses. This means that we can restrict our analysis to a single pulse.

The second layer nucleation probability $`p_{\mathrm{nuc}}(I,A)`$ is defined as the probability that there is at least one nucleation event on islands with area $`A`$ after deposition of a pulse with intensity $`I`$. For high ES barrier this means that during a single pulse at least two atoms have to be deposited on the same island. Obviously, the probability of depositing $`k`$ atoms on an island with area $`A`$ during one pulse is given by a Poisson distribution with average $`IA`$,

$$p_k=\frac{\left(IA\right)^k}{k!}e^{IA}.$$
(20)
Hence the nucleation probability can be expressed as

$`p_{\mathrm{nuc}}(I,A)`$ $`=`$ $`1p_0p_1=1e^{IA}\left(1+IA\right)`$ (21)
$``$ $`I^2A^2+O(I^3A^3).`$
Fig. 6a shows the numerical results for $`p_{\mathrm{nuc}}`$ for $`I=0.01`$ and different ES-barriers. As can be seen, for increasing ES barrier the measured curves approach the predicted one, (21), which is shown as a dashed line. The agreement is good for $`\mathrm{}_{ES}10^3`$, where the nucleation time is much shorter than the residence time. However, the agreement is not convincing for island areas $`A<10`$, where the discrete lattice spacing starts to play a role. Moreover, the upper curves deviate for $`A>40`$ due to coalescence of large islands. This explanation is supported by the results shown in Fig. 6b, where two different intensities are compared. Indeed, for the lower intensity $`I=0.005`$, where coalescence starts at larger island sizes, the measured and the predicted curves agree much better. Thus we conclude that the second layer nucleation probability $`p_{\mathrm{nuc}}`$ is adequately described by (21).

With Eq. (21) we can now answer which of the opposing trends in PLD – island size reduction and increasing nucleation probability on top of islands – will dominate. In the limit $`D/F\mathrm{}`$, where the critical intensity $`I_c`$ tends to zero, the island area scales as $`A\mathrm{}_D^2I^{2\nu }`$ for all intensities. Together with Eq. (21) one obtains to leading order

$$p_{\mathrm{nuc}}I^{2(12\nu )}.$$
(22)
Using the previous estimate $`\nu 0.26`$ the exponent is given by $`2(12\nu )1`$. Therefore, the second layer nucleation probability $`p_{\mathrm{nuc}}`$ grows with increasing intensity, enhancing the Villain instability. Contrarily, MBE with infinite $`D/F`$ always shows perfect layer-by-layer growth, even for high but finite ES-barriers. Thus, we conclude that in the limit $`D/F\mathrm{}`$ PLD cannot improve layer-by-layer growth. In order to apply this result to practical situations, one has to find the lower boundary for $`D/F`$, below which the behavior is different. This question will be addressed now.

### IV.2 Comparison of PLD and MBE for finite $`D/F`$

For perfect layer-by-layer growth the squared surface width $`w^2=(h^2h^2)`$ is known to oscillate between zero for completed monolayers and $`1/4`$ for half monolayers Wolf95 . Without ES barriers such oscillations can be seen in MBE as well as in PLD. As shown in Fig. 7a, they are most pronounced in MBE while in PLD they become more and more damped as the intensity increases. Moreover, it can be seen that in PLD without ES barriers the roughness always increases with increasing pulse intensity.

Fig. 7b shows the corresponding result for a very high Ehrlich-Schwoebel length $`\mathrm{}_{\mathrm{ES}}=10^4`$. There are no oscillations since the surface roughens very quickly due to the Villain instability Villain91 . Surprisingly, the width is now maximal for MBE and decreases with increasing pulse intensity, i.e., the trend is reversed. This reversal can even be observed visually by monitoring the interface at different intensities, as shown in Fig 8.

Let us now investigate this reversal in more detail. Since the oscillations in PLD are extremely weak or not present at all, a numerical measurement of the damping time (such as in Ref. Kallabis97  for MBE) is not feasible. Instead we monitored the squared surface width after deposition of a two monolayers – a coverage which is also relevant for experimental applications Jenniches96 . Although the choice of two monolayers is arbitrary and does not permit a rigorous quantitative analysis, this criterion is very simple and tells us for which parameters PLD produces smoother surfaces than MBE.

The roughness after deposition of two monolayers is plotted in Fig. 9. In agreement with previous results for MBE Brendel , all curves increase monotonously, i.e., the roughness increases with increasing Ehrlich-Schwoebel barrier. However, the curves for MBE and PLD cross each other. The typical ES barrier, where this crossing takes place, varies roughly as $`\mathrm{}_{\mathrm{ES}}(D/F)^{1/2}`$. This observation is in agreement with our previous result that for $`D/F\mathrm{}`$ surfaces grown by MBE are always smoother compared to PLD.

Thus there is indeed a range of parameters where PLD with a high intensity produces smoother surfaces than MBE, even if the atoms are deposited with thermal energy. However, the Ehrlich-Schwoebel barriers are unphysically large (typical experimental values are $`\mathrm{}_{ES}1\mathrm{}10`$). Therefore we believe that nonthermal energy deposition effects are important for explaining the experimentally observed growth improvement in PLD.

Nevertheless it is interesting to address the question, what causes the crossing of the curves in Fig. 9. As we have explained, second layer nucleation becomes more likely the larger the pulse intensity is in spite of the decreasing island size. Therefore one would always expect the roughness to increase with pulse intensity, as it is indeed the case for sufficiently small Ehrlich-Schwoebel barriers. For $`\mathrm{}_{\mathrm{ES}}\mathrm{}`$, however, there is no interlayer transport and the roughness is the same as in random deposition, irrespective of the pulse intensity. Thus, the curves in Fig. 9 are expected to saturate eventually at the same roughness. Thus the observed reversal results from an interplay of a large but finite ES barrier and a finite value of $`D/F`$.

To understand the reversal let us consider how the Villain instability unfolds for large but finite ES barriers. Depending on the intensity, the first pulse leads to the formation of many small islands. After the first pulse there are two temporal regimes of different roughening behavior. Because of the Villain instability effective uphill currents lead to a quick formation of large mounds. As new nucleations preferentially take place on top of these mounds, the nucleations are vertically aligned Somfai , as shown on the r.h.s. of Fig. IV.2. Therefore, the number and the lateral size of the mounds is essentially determined by the initial configuration of small islands after the first pulse. This determines at least for transient times also the height of the mounds and therefore the roughness of the surface. In our simulation model we do not observe slope selection. In a more realistic case one would expect that the mounds grow until their edges reach a critical slope where the uphill current becomes zero. Then the process enters a different temporal regime where the mounds compete with one another, leading to an effective coarsening process where the roughness increases only slowly.

By increasing the pulse intensity second layer nucleations become more likely so that the Villain instability is accelerated. However, the typical roughness from where on the process enters the second temporal regime of slow coarsening depends mainly on the initial density of islands after the first pulse. Therefore, this typical roughness decreases with increasing intensity. After sufficiently long time this effect dominates and leads to the observed crossings of the curves in Fig. 9. Obviously, this mechanism works in any dimension.

## V Conclusion

In the present paper we have studied a model for pulsed laser deposition with and without Ehrlich-Schwoebel barriers. The model assumes that the particles are deposited at thermal energies where the transient mobility of adatoms as well as implantation effects can be neglected. First we investigated the case without Ehrlich-Schwoebel barriers. For low pulse intensities PLD displays essentially the same behavior as MBE. Increasing the intensity it crosses over to a different behavior characterized by the nucleation of many small islands after the first pulse. The behavior in this regime can be analyzed by studying the time-dependent nucleation density. Extending previous results we have shown that in the limit $`D/F\mathrm{}`$ this quantity displays an unusual type of logarithmic scaling behavior.

Turning to PLD with Ehrlich-Schwoebel barriers without leakages we first showed that in the limit $`D/F\mathrm{}`$ the second layer nucleation probability grows with increasing intensity, enhancing the Villain instability. This means that in this limit PLD cannot improve layer-by-layer growth. For finite $`D/F`$, however, the situation is different. Studying the surface width for the Schwoebel lengths between $`1`$ and $`10^4`$ we found that PLD produces a smoother surface than MBE if the barrier is strong enough. This reversal can be explained by the influence of the initial nucleation density after the first pulse on the roughening due to the Villain instability. However, we believe that this mechanism alone cannot explain the experimentally observed improvement of layerwise growth in PLD where the barriers are much lower. In order to describe this effect more accurately the transient mobility of adatoms would have to be taken into account.

As we pointed out, Schinzer et al. Biehl  had reached a similar conclusion that pulsed deposition leads to increased roughness. However, this statement, which contradicts an idea by Rosenfeld et al. Rosenfeld  cannot be universally valid, as our more detailed investigation showed. In our opinion the negative effect of pulses is mainly due to the fact that the probability of multiple deposition on top of islands is increased. This indicates, that the very high intensity of the single pulse per monolayer was responsible for the roughness in the simulations of Schinzer et al. Biehl .

However, Schinzer et al. Biehl  also discovered a way to improve MBE with pulsed deposition: If desorption cannot be neglected, atoms will preferentially evaporate from the top terraces. Thereby the negative effect of multiple deposition can be compensated, so that one gains the benefit of enforcing small islands by the pulse in the beginning of the monolayer.

Our model describes only a limiting case of pulsed laser deposition. It has the virtue of making new scaling concepts clear. In practice the atoms are seldomly deposited with thermal energy in pulsed laser deposition. In particular it would be interesting to work out the effect of transient mobility due to heating by the pulse.

Acknowledgements:
We would like to thank J. Krug for pointing out Refs. Kadanoff88 ; Tang89 . We also thank the Deutsche Forschungsgemeinschaft for support within SFB 491.

## Appendix A Finite-Size effects in the limit of low intensity

As discussed in Sec. II, for small intensities PLD displays essentially the same growth behavior as MBE. However, on finite lattices both processes are not exactly equivalent. This can be seen in Fig. 10, where the average island distances for MBE and for PLD with an intensity of one particle per pulse $`I=L^2`$ are compared over a wide range of $`D/F`$. For very low values of $`D/F`$ the curves coincide but do not follow a power law because of lattice effects. Only in the range $`10^4D/F10^7`$ the curves follow approximately the expected power law $`\mathrm{}_D(D/F)^\gamma `$. Finally, for $`D/F10^7`$ the island distances of MBE and PLD differ increasingly from each other.

In order to demonstrate that this discrepancy between PLD and MBE is a finite size effect we determine the values of $`(D/F)_{\mathrm{dev}}`$ from where on the island distances in both models differ by a certain factor. As shown in Fig. 11 these values increase algebraically with the system size as

$$(D/F)_{\mathrm{dev}}L^\chi ,$$
(23)
where the exponent $`\chi `$ is found to be close to $`3`$. This power law behavior can be explained as follows. In both cases single atoms are deposited and the average time between two depositions $`\tau =1/L^2`$ is the same. However, in PLD the deposition takes place in constant time intervals, whereas in MBE atoms are randomly deposited so that the time intervals between deposition events obey a Poissonian distribution $`P(\tau )=L^2e^{L^2\tau }`$ with the variance $`\sigma =\sqrt{\tau ^2\tau ^2}=1/L^2`$. If such a fluctuation leads to a $`\tau `$ which is smaller than the average time interval between two depositions, more nucleations will be produced. As the formation of nucleations is irreversible this enhancement is not compensated by fluctuations in opposite direction where $`\tau `$ is large. Therefore, the influence of fluctuations increases the number of nucleations, leading to a smaller average island distance than in PLD. However, this effect can only be seen if the fluctuations are strong enough, i.e., they have to be at least of the same order as $`\mathrm{}_{D}^{}{}_{}{}^{2}/(D/F)`$, which is the average diffusion time (in ML) before an adatom reaches the edge of an island or another adatom. Since $`\mathrm{}_D(D/F)^\gamma `$ this argument leads to the scaling relation

$$(D/F)_{\mathrm{dev}}L^\chi ,\chi =\frac{2}{12\gamma }.$$
(24)
For $`\gamma =1/6`$ the exponent is $`\chi =3`$ while for DLA it is given by $`\chi 3.1`$. This result is in fair agreement with the numerically determined exponents in Fig. 11.


format: text

text:
# Theory of amplification of sound (acoustic phonons) by absorption of laser radiation in quantum wires with parabolic potential.

## 1 Introduction

It is well known that when a laser radiation is applied to a material, the number of acoustic phonons inside is varied with time. Studied this phenomenon can lead to new knowledge about the electron-phonon interaction mechanism, especially in low dimensional structures.
The theory of amplification of sound (acoustic phonons) by absorption of laser radiation in bulk semiconductors and low dimensional structures was studied in a number of papers . In  the problem was restricted to bulk semiconductor with monophoton absorption and multiphton absorption processes. In  the problem was considered for quantum wells with monophoton and multiphoton absorption processes, or in a quantizing magnetic field that do not have the restriction conditions of .
In this paper we study theoretically the amplification of sound by absorption of laser radiation in quantum wires with the elliptic cross section using the assymmetric parabolic confining potential:

$$V(x,y)=\frac{m^{}}{2}(\mathrm{\Omega }_x^2x^2+\mathrm{\Omega }_y^2y^2),$$

Here $`\mathrm{\Omega }_x,\mathrm{\Omega }_y`$ are the effective frequencies of the potential, $`m^{}`$ is the effective electron mass. The energy of electron in this model has the form :

$$\epsilon _{n,l}(p_z)=\sqrt{\frac{\mathrm{\Omega }_x}{2}}(n+\frac{1}{2})+\sqrt{\frac{\mathrm{\Omega }_y}{2}}(l+\frac{1}{2})+\frac{p_z^2}{2m^{}}.$$

Effect of scattering on the impurities and the reflection of electron mode from boundary between the wire and the electron reservoir are not considered.

## 2 The general formula

The Hamiltonian for electron-phonon system in a quantum wire in external field can be written as:

$$\begin{array}{c}H(t)=\underset{n,l,\stackrel{}{k}}{}\epsilon _{n,l}\left(\stackrel{}{k}\frac{e}{c}\stackrel{}{A}(t)\right)+\underset{\stackrel{}{q}}{}\omega _\stackrel{}{q}b_\stackrel{}{q}^+b_\stackrel{}{q}+\hfill \\ \hfill +\underset{n,l,n^{},l^{},\stackrel{}{k},\stackrel{}{q}}{}C_{n,l,n^{},l^{}}(\stackrel{}{q})a_{n,l,\stackrel{}{k}+\stackrel{}{q}}^+a_{n^{},l^{},\stackrel{}{k}}(b_\stackrel{}{q}^++b_\stackrel{}{q}),\end{array}$$
(1)
where $`a_{n,l,\stackrel{}{k}}^+`$ and $`a_{n,l,\stackrel{}{k}}`$ ($`b_\stackrel{}{q}^+`$ and $`b_\stackrel{}{q}`$ ) are the creation and annihilation operators of electron (phonon); $`\stackrel{}{k}`$ is the electron wave vector (along the wire’s axis: z axis); $`\stackrel{}{q}`$ is the phonon wave vector; $`C_{n,l,n^{},l^{}}`$is the interaction constant of electron-acoustic phonon scattering; $`\stackrel{}{A}(t)=\frac{c}{\mathrm{\Omega }}E_0cos(\mathrm{\Omega }t)`$ is the potential vector, depend on the external field. Process using the method in , we obtain the quantum transport equation for acoustic phonons in quantum wires:

$$\begin{array}{c}\frac{}{t}b_\stackrel{}{q}_t+i\omega _\stackrel{}{q}b_\stackrel{}{q}_t=\underset{n,l,n^{},l^{}}{}|C_{n,l,n^{},l^{}}(\stackrel{}{q})|^2\underset{\stackrel{}{k}}{}(n_{n,l}(\stackrel{}{k}\stackrel{}{q})n_{n^{},l^{}}(\stackrel{}{k}))\times \hfill \\ \hfill \times _{\mathrm{}}^tb_\stackrel{}{q}_{t_1}\underset{\nu ,s=\mathrm{}}{\overset{\mathrm{}}{}}J_\nu \left(\frac{\lambda }{\mathrm{\Omega }}\right)J_s\left(\frac{\lambda }{\mathrm{\Omega }}\right)\times \\ \hfill \times exp\left(i[\epsilon _{n^{},l^{}}(\stackrel{}{k})\epsilon _{n,l}(\stackrel{}{k}\stackrel{}{q})](t_1t)i\nu \mathrm{\Omega }t_1+is\mathrm{\Omega }t\right)dt_1\end{array}$$
(2)
Where the symbol $`<x>_t`$ means the usual thermodynamics average of operator x; $`n_{n,l}(\stackrel{}{k})`$ is the distribution function of electron; $`\lambda =\frac{e\stackrel{}{q}\stackrel{}{E_0}}{m^{}\mathrm{\Omega }}`$. All formulae are written in units where $`\mathrm{}=1`$
Using Fourier transformation, with $`\delta +0`$, we have the dispersion equation:

$$\begin{array}{c}\omega \omega _\stackrel{}{q}=\pi \underset{n,l,n^{},l^{}}{}|C_{n,l,n^{},l^{}}(\stackrel{}{q})|^2\underset{\stackrel{}{k}}{}[n_{n^{},l^{}}(\stackrel{}{k})n_{n,l}(\stackrel{}{k}\stackrel{}{q})]\times \hfill \\ \hfill \times \underset{\nu =\mathrm{}}{\overset{\mathrm{}}{}}J_\nu ^2\left(\frac{\lambda }{\mathrm{\Omega }}\right)\frac{1}{\epsilon _{n^{},l^{}}(\stackrel{}{k})\epsilon _{n,l}(\stackrel{}{k}\stackrel{}{q})\omega _\stackrel{}{q}\nu \mathrm{\Omega }i\delta }.\end{array}$$
(3)
From this we obtain the general formula for absorption coefficient of acoustic phonons in quantum wires:

$$\begin{array}{c}\alpha (\stackrel{}{q})=\pi \underset{n,l,n^{},l^{}}{}|C_{n,l,n^{},l^{}}(\stackrel{}{q})|^2\underset{\stackrel{}{k}}{}[n_{n^{},l^{}}(\stackrel{}{k})n_{n,l}(\stackrel{}{k}\stackrel{}{q})]\times \hfill \\ \hfill \times \underset{\nu =\mathrm{}}{\overset{\mathrm{}}{}}J_\nu ^2\left(\frac{\lambda }{\mathrm{\Omega }}\right)\delta \{\epsilon _{n^{},l^{}}(\stackrel{}{k})\epsilon _{n,l}(\stackrel{}{k}\stackrel{}{q})\omega _\stackrel{}{q}\nu \mathrm{\Omega }\},\end{array}$$
(4)
where $`\delta (x)`$ is the Dirac equation.

## 3 The amplification of sound

In the following, we will calculate the absorption coefficient of sound by absorption of laser radiation in quantum wires in the two cases: monophoton absorption process and multiphoton absorption process.

### 3.1 Monophoton absorption process

From Eq.(4), assuming the parameter in Bessel function is small enough: $`\lambda \mathrm{\Omega }`$ and the electron is non-degenerate, we obtain the expression for the sound absorption coefficient:

$$\begin{array}{c}\alpha (\stackrel{}{q})=\frac{m^{}n_0\lambda ^2}{4q\mathrm{\Omega }^2}\underset{n,l,n^{},l^{}}{}|C_{n,l,n^{},l^{}}(\stackrel{}{q})|^2\mathrm{\Lambda }_{n,l,n^{},l^{}}(\stackrel{}{q},\mathrm{\Omega })\times \hfill \\ \hfill \times exp\left(\frac{\beta }{2}\left(\sqrt{\frac{\mathrm{\Omega }_x}{2}}\left(n^{}+\frac{1}{2}\right)+\sqrt{\frac{\mathrm{\Omega }_y}{2}}\left(l^{}+\frac{1}{2}\right)\right)\frac{m^{}\beta }{2q^2}(a^2+\mathrm{\Omega }^2)\right),\end{array}$$
(5)
where $`n_0`$ is the electron density,

$$a=\sqrt{\frac{\mathrm{\Omega }_x}{2}}(nn^{})+\sqrt{\frac{\mathrm{\Omega }_y}{2}}(ll^{})+\omega _\stackrel{}{q}+\frac{q^2}{2m^{}},$$

and

$$\mathrm{\Lambda }_{n,l,n^{},l^{}}(\stackrel{}{q},\mathrm{\Omega })=e^{\frac{\beta \omega _\stackrel{}{q}}{2}}\left(e^{\frac{\beta m^{}a\mathrm{\Omega }}{q^2}+\frac{\beta \mathrm{\Omega }}{2}}sh\left(\beta \frac{\omega _\stackrel{}{q}+\mathrm{\Omega }}{2}\right)+e^{\frac{\beta m^{}a\mathrm{\Omega }}{q^2}\frac{\beta \mathrm{\Omega }}{2}}sh\left(\beta \frac{\omega _\stackrel{}{q}\mathrm{\Omega }}{2}\right)\right).$$

Due to the $`\delta `$ function in (4), the momentum must satisfy:

$$\stackrel{}{k}\frac{\stackrel{}{q}}{2}+\frac{m^{}}{\stackrel{}{q}}\left(\sqrt{\frac{\mathrm{\Omega }_x}{2}}(nn^{})+\sqrt{\frac{\mathrm{\Omega }_y}{2}}(ll^{})+\omega _\stackrel{}{q}\right)+\frac{\mathrm{\Omega }m^{}}{\stackrel{}{q}}.$$
(6)
From Eq.(5), when $`\omega _\stackrel{}{q}\mathrm{\Omega }`$, we see that the absorption coefficent is negative, or it is the amplification of sound: $`\alpha (\stackrel{}{q})<0`$:

$$\begin{array}{c}\alpha (\stackrel{}{q})=\frac{m^{}n_0\lambda ^2}{2\stackrel{}{q}\mathrm{\Omega }^2}\underset{n,l,n^{},l^{}}{}|C_{n,l,n^{},l^{}}(\stackrel{}{q})|^2exp(\frac{\beta \omega _\stackrel{}{q}}{2}\frac{m^{}\beta \mathrm{\Omega }^2}{2q^2})\times \hfill \\ \hfill \times exp[\frac{\beta }{2}(\sqrt{\frac{\mathrm{\Omega }_x}{2}}(n^{}+\frac{1}{2})+\sqrt{\frac{\mathrm{\Omega }_y}{2}}(l^{}+\frac{1}{2}))\frac{m^{}\beta a^2}{2q^2}]\times \\ \hfill \times sh\left(\frac{\beta \mathrm{\Omega }}{2}\right)sh\left(\frac{\beta m^{}\mathrm{\Omega }}{q^2}\left(\omega _\stackrel{}{q}+\sqrt{\frac{\mathrm{\Omega }_x}{2}}(nn^{})+\sqrt{\frac{\mathrm{\Omega }_y}{2}}(ll^{})\right)\right).\end{array}$$
(7)
### 3.2 Multiphoton absorption process

Assume the parameter in Bessel function is large: $`\lambda \mathrm{\Omega }`$, use the approximate formula :

$$\underset{\nu }{}J_\nu ^2\left(\frac{\lambda }{\mathrm{\Omega }}\right)\delta (E\nu \mathrm{\Omega })=\frac{\theta (\lambda ^2E^2)}{\pi \sqrt{\lambda ^2E^2}},$$

with

$$\theta =\{\begin{array}{cc}1,\hfill & \text{if }x>0\hfill \\ 0,\hfill & \text{if }x<0\hfill \end{array}.$$

From the general formula for absorption coefficient of sound (4), we have:

$$\begin{array}{c}\alpha (\stackrel{}{q})=\frac{\pi ^{\frac{3}{2}}n_0m^{}}{2q}e^{\frac{\beta m^{}\lambda ^2}{2q^2}}\underset{n,l,n^{},l^{}}{}|C_{n,l,n^{},l^{}}(\stackrel{}{q})|^2\hfill \\ \hfill \underset{\nu =0}{\overset{\mathrm{}}{}}\frac{\mathrm{\Gamma }\left(\nu +\frac{1}{2}\right)}{\nu !}\left(\chi _{n,l}\left(\frac{q^2}{2m^{}}\right)\chi _{n^{},l^{}}\left(\frac{q^2}{2m^{}}\right)\right),\end{array}$$
(8)
where:

$$\begin{array}{c}\chi _{n,l}(x)=exp\left[\frac{\beta }{2}(\sqrt{\frac{\mathrm{\Omega }_x}{2}}(n^{}+\frac{1}{2})+\sqrt{\frac{\mathrm{\Omega }_y}{2}}(l^{}+\frac{1}{2}))\right]\times \hfill \\ \hfill \times exp[\frac{\beta m^{}}{2q^2}(\sqrt{\frac{\mathrm{\Omega }_x}{2}}(nn^{})+\sqrt{\frac{\mathrm{\Omega }_y}{2}}(ll^{})+x+\omega _\stackrel{}{q})^2]\times \\ \hfill \times \left(\frac{\lambda }{\sqrt{\frac{\mathrm{\Omega }_x}{2}}(nn^{})+\sqrt{\frac{\mathrm{\Omega }_y}{2}}(ll^{})+x+\omega _\stackrel{}{q}}\right)^\nu I_\nu \left(\frac{\beta m^{}\lambda }{q^2}\right),\end{array}$$

$`I_\nu (x)`$ is the complex Bessel function of order $`\nu `$
. From the $`\theta `$ function above, we derive the momentum condition:

$$k<\pm q+\frac{m^{}}{q}\left(\sqrt{\frac{\mathrm{\Omega }x}{2}}(nn^{})+\sqrt{\frac{\mathrm{\Omega }y}{2}}(ll^{})+\omega _\stackrel{}{q}\lambda \right).$$
(9)
Note that if $`\chi _{n,l}\left(\frac{q^2}{2m^{}}\right)<\chi _{n^{},l^{}}\left(\frac{q^2}{2m^{}}\right)`$, then $`\alpha (\stackrel{}{q})<0`$, which mean we have the amplification of sound.

## 4 Discussion and numerical results

Let us begin by discussing the results Eq.(5) and Eq.(8), which give the absorption coefficient of sound (acoustic phonons) in quantum wires with parabolic potential in the two cases: monophoton absorption and multiphoton absorption respectively. These expressions are carried out only when the conditions (6) and (9), which demonstrate the conservation of energy-momentum, are satisfied. The results are much different from previous work on bulk semiconductors  due to the effect of confinement of electron in low dimensional structures. In quantum wires, the electron’s mobility is enhanced, what leads to an unusual behaviour under external stimuli (in this case is the laser radiation). Furthermore, they are also different from results with quantum wells , and quantum wires with cylindrical potential  due to the complexity of the electron energy and wave function. The formula (5) demonstrates the relation of the absorption coefficient on laser radiation with order two (the $`\lambda =\frac{e\stackrel{}{q}\stackrel{}{E_0}}{m^{}\mathrm{\Omega }}`$ parameter), whereas Eq.(8) shows the same relation with order greater than two (since the complex Bessel function $`I_\nu (x)`$ contain $`\lambda `$). Note that in both cases, in proper conditions, the absorption coefficient is negative, which means the number of phonon increases with time.

From the obtained results, we plot for the dependence of the absorption of sound $`\alpha `$ on wave vector $`\stackrel{}{q}`$ and laser frequency $`\mathrm{\Omega }`$ for the case of the monophoton absorption process. The plot is for the case $`nn^{}`$, $`ll^{}`$, that is for interband absorption, $`\beta =k_BT=0.05meV,r_0=10^8m,m^{}=0.0067m_0,m_0`$ is the electron mass. The figure shows that the amplification of sound is maximum at $`q=0.0005`$ (the minimum on the graph). Compare this plot with those of cylindrical quantum wires, we see that the maximum and the minimum of $`\alpha `$ happen at the same values of $`\stackrel{}{q}`$, however, the variation of $`\alpha `$ is slightly difference: it decreases more slowly in the interval of the wavevector from 0.0015 meV to 0.0025 meV.

## 5 Conclusion

In this paper, the amplification of sound (acoustic phonons) by absorption of laser radiation is investigated in the two cases: monophoton absorption process and multiphoton absorption process. The analytical results are difference from previous work, however, the numerical results are similar to those in the case of cylindrical quantum wires. Thus a change in potential does not affect much the electron-phonon scattering in quantum wires.


text:
# The High 𝐸_𝑇 Drop of 𝐽/𝜓 to Drell-Yan Ratio from the Statistical 𝑐𝑐̄ Coalescence Model

## Abstract

The dependence of the $`J/\psi `$ yield on the transverse energy $`E_T`$ in heavy ion collisions is considered within the statistical $`c\overline{c}`$ coalescence model. The model fits the NA50 data for Pb+Pb collisions at the CERN SPS even in the high-$`E_T`$ region ($`E_T100`$ GeV). Here $`E_T`$-fluctuations and $`E_T`$-losses in the dimuon event sample naturally create the celebrated drop in the $`J/\psi `$ to Drell-Yan ratio.

Measurements of the $`J/\psi `$ suppression pattern in heavy ion collisions had been proposed as a diagnostic tool for the deconfinement of strongly interacting matter at the early stage of the reaction  (see also  for a modern review). However, it was found that the $`J/\psi `$ suppression in proton-nucleus and nucleus-nucleus collisions with light projectiles (up to S+U) can be explained by inelastic collisions with the nucleons of the incident nuclei (the so-called ‘normal nuclear suppression’) . In contrast, the NA50 experiments with a heavy projectile and a heavy target (Pb+Pb) revealed ‘anomalous’ $`J/\psi `$ suppression , i.e. an essentially stronger drop of the $`J/\psi `$ over Drell-Yan ratio $`R`$ than it could be expected from a simple extrapolation of the light projectile data. The dependence of $`R`$ on the neutral transverse energy $`E_T`$, which is used as an estimator of the collision centrality, shows a two step behavior. The onset  of the anomalous $`J/\psi `$ suppression (deviation of the data downwards from the ‘normal nuclear suppression’ curve) supposedly takes place at $`E_T35`$ GeV. The second drop has been noted at $`E_T100`$ GeV. The initial interpretation  attributed these two drops to the successive melting of the charmonium $`\chi _c`$ states (which contribute about 45% to the measured $`J/\psi `$ number) and primary $`J/\psi `$’s in the quark-gluon plasma.

Other interpretations are, however, possible. Because of relatively high statistical errors of the data at low $`E_T`$, the behavior around $`E_T40`$ GeV can be equally well interpreted as smooth decreasing of the measured ratio $`R(E_T)`$ due to charmonium suppression by comoving hadrons***Later it was found that even the sweeping nucleons of the colliding nuclei might, under certain assumptions, produce the observed suppression . . Discussing the second drop of the $`J/\psi `$ to Drell-Yan ratio, one has to take into account  that the neutral transverse energy $`E_T`$ estimates different quantities in the regions below and above $`E_T100`$ GeV. At $`E_T100`$ GeV, it is proportional with good accuracy to the number of nucleon participants $`N_p`$ and provides a reliable measure of the collision centrality. The large $`E_T`$ ‘tail’, however, is produced mainly due to fluctuations of the transverse energy at almost fixed (nearly maximal) centrality. The influence of these fluctuations on the $`J/\psi `$ suppression can partially explain the second drop. Better agreement of the co-mover model with the data in the region $`E_T100`$ is achieved when a loss of transverse energy in the dimuon data sample with respect to the minimum bias one is also taken into account .

The $`E_T`$ fluctuations appear to be important in the deconfinement scenario as well. It was found  that while the first drop can be explained by melting of $`\chi _c`$’s, the onset of primary $`J/\psi `$ melting does not produce any visible structure on the curve. The second drop appears due to the fluctuations of $`E_T`$.

All mentioned approaches , despite of their differences, have one common feature: charmonia are assumed to be created exclusively at the initial stage of the reaction in primary nucleon collisions. Any subsequent interactions (with sweeping nucleons, co-moving hadrons or quark-gluon medium) may only destroy them. Despite of rather successful agreement with the $`J/\psi `$ data, such a scenario seems to be in trouble explaining the $`\psi ^{}`$ yield. The recent lattice simulations  suggest that the temperature of $`\psi ^{}`$ dissociation $`T_d(\psi ^{})`$ lies far below the deconfinement point $`T_c`$: $`T_d(\psi ^{})0.1`$$`0.2T_c`$ . Therefore, not only the quark-gluon plasma, but also a hadronic co-mover medium should completely eliminate $`\psi ^{}`$ charmonia in central Pb+Pb collisions at SPS. However, the experiment revealed a sizable $`\psi ^{}`$ yield (see, for instance, ). It was observed  that $`\psi ^{}`$ to $`J/\psi `$ ratio decreases with centrality only in peripheral lead-lead collisions, but remains approximately constant and equal to its thermal value at sufficiently large number of participants $`N_p100`$.

A natural way to explain the $`\psi ^{}`$ yield in nucleus-nucleus collisions without running into contradiction with the lattice data is to assume that a charmed quark and an antiquark that have been created at the initial (‘hard’) stage of the reaction can coalesce to form a hidden charm meson at a later stage, close to the chemical freeze-out point. This possibility has been considered in the literature  within various models.

The statistical $`c\overline{c}`$ coalescence model of $`J/\psi `$ production suggested in  will be considered in the present paper. The charmed quark-antiquark ($`c\overline{c}`$) pairs are assumed to be created at the initial stage of the nucleus-nucleus reaction in hard parton collisions. Creation of $`c\overline{c}`$ pairs after the hard initial stage as well as their possible annihilation is neglected. Therefore, the number of charmed quark-antiquark pairs remains approximately unchanged during the subsequent evolution. They are then distributed over open and hidden charm particles at the final stage of the reaction in accordance with laws of statistical mechanics. It is postulated that the deconfined medium prevents formation of charmonia at early stage of the reaction. In fact, it was found  that $`J/\psi `$ becomes unbound slightly above the deconfinement point $`T_d(J/\psi )1.1T_c`$. Therefore, all charmonia observed in central Pb+Pb collisions supposed to be created at the final stage of the reaction.

It was shown in our previous paper  that the statistical $`c\overline{c}`$ coalescence model is in excellent agreement with the NA50 data for (semi)central ($`N_p100`$) Pb+Pb collisions. However, the analysis of  was restricted to the region $`E_T100`$ GeV, where $`E_T`$ estimates the collision centrality. Fluctuations of $`E_T`$ were not taken into account. The aim of the present paper is to extend the applicability domain of the model to $`E_T100`$ GeV and describe the second drop of the anomalous $`J/\psi `$ suppression. Two effects appear to be important in this region:

* Influence of fluctuations of thermodynamic parameters of the system on the $`J/\psi `$ production.
* A loss of neutral transverse energy, $`E_T`$, in the event sample collected with the dimuon trigger relative to the minimum bias events.

As far as charmed (anti)quarks are created exclusively at the initial ‘hard’ stage of the reaction in collisions of primary nucleons, the average number of $`c\overline{c}`$ pairs $`c\overline{c}_{AB(b)}`$ created in a nucleus-nucleus (A+B) collision at fixed impact parameter $`b`$ is proportional to the number of primary nucleon-nucleon collisions and can be found in Glauber’s approach:

$$c\overline{c}_{AB(b)}=AB\sigma _{c\overline{c}}^{NN}T_{AB}(b).$$
(1)
Here $`T_{AB}(b)`$ is the nuclear overlap function calculated assuming a Woods-Saxon distribution of nucleons in nuclei  (see e.g. Appendix of ). $`\sigma _{c\overline{c}}^{NN}`$ is the production cross section of a $`c\overline{c}`$ pair in a nucleon-nucleon collision. There are indirect experimental indications  and a phenomenological model  suggesting an enhancement of the open charm in heavy ion collisions with respect to a direct extrapolation of nucleon-nucleon data. Therefore, $`\sigma _{c\overline{c}}^{NN}`$ in nucleus-nucleus collisions is not necessarily equal to its value measured in a proton-proton experiment. We consider $`\sigma _{c\overline{c}}^{NN}`$ as a free parameter. Its value is fixed by fitting the NA50 data.

The primary nucleon-nucleon collisions in a nucleus-nucleus reaction are independent and the probability to produce a $`c\overline{c}`$ pair in a single collision is small. Therefore, event-by-event fluctuations of the number of $`c\overline{c}`$ pairs in nucleus-nucleus collisions at fixed impact parameter $`b`$ are Poisson distributed. Assuming exact conservation of this number one obtains the following statistical coalescence model result for the average number of produced $`J/\psi `$’s per A+B collision :

$$J/\psi _{AB(b)}^{(T,\mu _B)}=c\overline{c}_{AB(b)}\left[1+c\overline{c}_{AB(b)}\right]\frac{N_{J/\psi }^{tot}(T,\mu _B)}{(N_O(T,\mu _B)/2)^2}+\mathrm{},$$
(2)
Where $`T`$ and $`\mu _B`$ are, respectively, the system temperature and baryonic chemical potential. The dots in (2) state for higher order terms with respect to the ratio $`\frac{N_{J/\psi }^{tot}(T,\mu _B)}{(N_O(T,\mu _B)/2)^2}`$, which can be safely neglected at $`N_p>100`$. Here $`N_O`$ is the total open charm multiplicity calculated within the grand canonical ensemble of the equilibrium hadron gas model:

$$N_O=V\underset{j=D,\overline{D},D^{},\overline{D}^{},\mathrm{}}{}n_j(T,\mu _B)Vn_O(T,\mu _B).$$
(3)
The sum in (3) runs over all known (anti)charmed particle species . The total $`J/\psi `$ multiplicity $`N_{J/\psi }^{tot}(T,\mu _B)`$ includes the contribution of excited charmonium states decaying into $`J/\psi `$:

$$N_{J/\psi }^{tot}=V\underset{j=J/\psi ,\chi _1,\chi _2,\psi ^{}}{}R(jJ/\psi )n_j(T,\mu _B)Vn_{J/\psi }^{tot}(T,\mu _B).$$
(4)
Here $`R(jJ/\psi )`$ is the decay branching ratio of the charmonium $`j`$ into $`J/\psi `$: $`R(J/\psi J/\psi )1`$, $`R(\chi _1J/\psi )0.27`$, $`R(\chi _2J/\psi )0.14`$ and $`R(\psi ^{}J/\psi )0.54`$. The volume parameter $`V`$ is fixed by the condition of the baryon number conservation:

$$N_p(b)=V\underset{j=N,\overline{N},\mathrm{\Delta },\overline{\mathrm{\Delta }},\mathrm{}}{}b_jn_j(T,\mu _B)Vn_B(T,\mu _B).$$
(5)
Here $`b_j`$ is the baryon number of particle type $`j`$ and the sum in (5) runs over all known (anti)baryon species . The number $`N_p(b)`$ of nucleon participants at fixed impact parameter $`b`$ is found in Glauber’s approach (see e.g. Appendix of )In what follows we neglect fluctuations of the number of nucleon participants at fixed impact parameter..

The particle densities in (35) are found from:

$$n_j(T,\mu _B)=\frac{d_j}{2\pi ^2}_0^{\mathrm{}}k^2𝑑k\left[\mathrm{exp}\left(\frac{\sqrt{m_j^2+k^2}\mu _j}{T}\right)\pm 1\right]^1,$$
(6)
where $`m_j`$ and $`d_j`$ are, respectively, the mass and degeneracy factor of particle type $`j`$; $`\mu _j=b_j\mu _B+s_j\mu _S+c_j\mu _C`$ is the chemical potential with $`s_j`$ and $`c_j`$ being, respectively, the strangeness and charm of particle type $`j`$. The strange and charm chemical potentials $`\mu _S`$ and $`\mu _C`$ are fixed by the requirement of zero net strangeness and charm in the system.

As was already discussed above, the analysis  of the recent lattice data  has shown that $`\psi ^{}`$ cannot exist in hot hadronic medium. To explain the observed $`\psi ^{}`$ multiplicity one has to assume that charmonia are formed close to the chemical freeze-out point. (For (semi)central Pb+Pb collisions at SPS the chemical freeze-out (almost) coincides with the hadronization.) The thermodynamical parameters $`T`$ and $`\mu _B`$ at chemical freeze-out can be estimated by fitting the hadron gas model to the hadron yield data in Pb+Pb collisions at SPS . This procedure, however, provides only average (in a sense that will be defined later) values $`\overline{T}`$ and $`\overline{\mu }_B`$. In fact, the thermodynamic parameters are not the same for all events, but rather they are subject to fluctuations of dynamical origin.

Only a fraction of the initial energy of the incident nuclei is converted to thermal energy $`E_{th}`$ (including the rest masses) of the produced hadron gas. The rest is redistributed to kinetic energy of the collective flow. The contribution of each nucleon participant to $`E_{th}`$ is not fixed. It fluctuates in accord with a hitherto unknown dynamical mechanism.

Suppose one is able to measure both $`N_p`$ and $`E_{th}`$. Then, events can be grouped in bins with different thermal energy per participant, $`E_{th}/N_p`$. Different thermal energies per net baryon number cannot be obtained at the same values of the temperature and baryonic chemical potential of the hadron gas. Therefore, fitting the hadron gas model to the hadron yield data in each bin, one would obtain different values of $`T`$ and $`\mu _B`$. In reality, however, such a selection of events is not feasible, therefore no experimental information about fluctuations of $`T`$ and $`\mu _B`$ is available. The procedure of  has to mix events with different $`E_{th}/N_p`$ and estimates average values $`\overline{T}`$ and $`\overline{\mu }_B`$. One can, however, make a reasonable assumption about the fluctuations. It was found that at different initial energies of colliding nuclei, the produced hadron gas chemically freezes-out at constant total thermal energy per hadron :

$$E_{th}(T,\mu _B)/N_h(T,\mu _B)=\text{const}1\text{ GeV}.$$
(7)
Different initial energies result in different thermal energy per participant. Therefore, it is natural to assume that the same criterion is valid when $`E_{th}/N_p`$ fluctuates at fixed initial energy.

At given $`T`$ and $`\mu _B`$, the thermal energy of the hadron gas can be found from

$$E_{th}(T,\mu _B)=V\underset{j}{}\epsilon _j(T,\mu _j),$$
(8)
where $`\epsilon _j(T,\mu _B)`$ is the contribution of the hadron species $`j`$ to the thermal energy density:

$$\epsilon _j(T,\mu _B)=\frac{d_j}{2\pi ^2}_0^{\mathrm{}}k^2𝑑k\sqrt{m_j^2+k^2}\left[\mathrm{exp}\left(\frac{\sqrt{m_j^2+k^2}\mu _j}{T}\right)\pm 1\right]^1.$$
(9)
The number of hadrons $`N_h(T,\mu _B)`$ is given by

$$N_h(T,\mu _B)=V\underset{j}{}n_j(T,\mu _j).$$
(10)
The sum in equations (8) and (10) runs over all known stable hadron and resonance species .

Equation (7) defines the curve in the $`(T,\mu _B)`$ plane along which the freeze-out parameters fluctuate. To compare our model with the experiment the relation of each point on this curve to the transverse energy $`E_T`$ of neutral hadrons at given impact parameter $`b`$ is needed. The point $`(\overline{T},\overline{\mu }_B)`$ corresponds to the average value $`\overline{E}_T`$, which in the framework of wounded nucleon model  is proportional to the number of participant nucleons:

$$\overline{E}_T=qN_p(b).$$
(11)
The fluctuations of $`E_T`$ are known to be Gaussian-distributed :

$$P_{MB}(E_T|b)=\frac{1}{\sqrt{2\pi q^2aN_p(b)}}\mathrm{exp}\left(\frac{[E_TqN_p(b)]^2}{2q^2aN_p(b)}\right).$$
(12)
The subscript ‘MB’ states for ‘minimum bias’. This probability distribution is valid for the event sample collected with the ‘minimum bias’ trigger, which includes allIn contrast, the dimuon trigger requires production of a dimuon pair. events having the given value of $`E_T`$. The parameter values $`q=0.274`$ GeV and $`a=1.27`$  are fixed from the minimum bias transverse energy distribution . The deviations of $`E_T`$ from its central value are mostly due to change of the number of neutral hadrons rather than due to change of average energy per hadron. Therefore, $`E_T`$ is approximately proportional to the total multiplicity of final hadrons: $`E_TN_h^{tot}(T,\mu _B)`$. Due to resonance decay, $`N_h^{tot}(T,\mu _B)`$ is larger than the number of particles in the hadron gas at freeze-out $`N_h(T,\mu _B)`$. Namely, every term in the sum (10) should be multiplied by a factor $`g_j1`$:

$$g_j=\{\begin{array}{cc}1\hfill & \text{for stable and weakly decaying particles;}\hfill \\ _kkR(jkh)\hfill & \text{for resonances.}\hfill \end{array}$$
(13)
Here $`R(jkh)`$ is the probability that the resonance $`j`$ produces $`k`$ hadrons in the final state.

Hence, the average number of $`J/\psi `$ mesons, $`J/\psi _{AB(b)}^{(E_T)}`$, produced in an $`A+B`$ collision at fixed impact parameter $`b`$ and transverse energy $`E_T`$ can be found in the following way:

1. $`T`$ and $`\mu _B`$ are found as a solution of coupled transcendental equations

$`E_{th}(T,\mu _B)/N_h(T,\mu _B)`$ $`=`$ $`E_{th}(\overline{T},\overline{\mu }_B)/N_h(\overline{T},\overline{\mu }_B),`$ (14)
$`N_h^{tot}(T,\mu _B)/N_h^{tot}(\overline{T},\overline{\mu }_B)`$ $`=`$ $`E_T/\overline{E}_T.`$ (15)
The first of the above equations expresses the freeze-out criterion (7): constant thermal energy per hadron. The second one appears due to the proportionality between the multiplicity of final particles and the measured neutral transverse energy.
2. The values of $`T`$ and $`\mu _B`$ found from (14) and (15) are substituted into the formula (2) giving the desired value of $`J/\psi _{AB(b)}^{(E_T)}`$.

In real experiment, the impact parameter is not fixed. The relevant quantity, the differential cross section of the $`J/\psi `$ production in collisions of nuclei A and B with respect to the transverse energy $`E_T`$, is obtained by integrating $`J/\psi _{AB(b)}^{(E_T)}`$ over the impact parameter $`b`$:

$$\frac{d\sigma _{ABJ/\psi }}{dE_T}=2\pi _0^{\mathrm{}}𝑑bbJ/\psi _{AB(b)}^{(E_T)}P_{J/\psi }(E_T|b).$$
(16)
Here $`P_{J/\psi }(E_T|b)`$ is the probability to produce neutral transverse energy $`E_T`$ at fixed impact parameter $`b`$, provided that a $`J/\psi `$ particle is produced in the same event.

Formula (16) refers to the total number of produced $`J/\psi `$’s, while the measured quantity is the number of $`\mu ^+\mu ^{}`$ pairs originating from their decay and satisfying certain kinematics conditions. Therefore, the cross-section (16) should be multiplied by the probability of $`J/\psi `$ decay into a dimuon pair $`B_{\mu \mu }^{J/\psi }=(5.88\pm 0.10)\%`$  and by a factor $`\eta `$, the fraction of the dimuons satisfying the kinematical conditions of the NA50 spectrometer. Because the considered $`J/\psi `$ production mechanism is completely different from the “standard” production in hard collisions, the rapidity distribution of $`J/\psi `$’s may be very different from that in N+N collisions. Consequently, $`\eta `$ should not be necessarily equal to $`\eta _{NN}0.24`$ which is obtained from Schuler’s parameterization  of proton-proton data. Therefore, $`\eta `$ will be treated as one more free parameter.

Finally, to make a comparison with the NA50 data, the $`J/\psi `$ cross section should be divided by the Drell-Yan cross section. Similarly to (16), it is found from

$$\frac{d\sigma _{ABDY^{}}}{dE_T}=d^2bDY^{}_{AB(b)}P_{DY}(E_T|b).$$
(17)
Here quantity $`P_{DY}(E_T|b)`$ is similar to $`P_{J/\psi }(E_T|b)`$, $`DY^{}_{AB(b)}`$ is the average number of Drell-Yan pairs produced in an A+B collision at fixed impact parameter $`b`$. The prime means that only the pairs satisfying the kinematical conditions of the NA50 spectrometer are taken into account. Similarly to $`c\overline{c}`$ pairs, the number of Drell-Yan pairs is proportional to the number of primary nucleon-nucleon collisions:

$$DY^{}_{AB(b)}=AB\sigma _{DY^{}}^{NN}T_{AB}(b),$$
(18)
where $`\sigma _{DY^{}}^{NN}`$ is the production cross section of Drell-Yan pairs in nucleon-nucleon collisions. This quantity is isospin dependent, therefore the average value should be used:

$$\sigma _{DY^{}}^{NN}=\frac{\sigma _{DY^{}}^{AB}}{AB}.$$
(19)
For the case of Pb+Pb collisions, $`A=B=208`$ and $`\sigma _{DY^{}}^{PbPb}=1.49\pm 0.13`$ $`\mu `$b .

Hence, the following formula will be used to fit the NA50 data:

$$R(E_T)=\eta B_{\mu \mu }^{J/\psi }\left(\frac{d\sigma _{ABJ/\psi }}{dE_T}\right)/\left(\frac{d\sigma _{ABDY^{}}}{dE_T}\right).$$
(20)
As the first step, we assume that both $`P_{J/\psi }(E_T|b)`$ and $`P_{DY}(E_T|b)`$ in (16) and (17) coincide with the ’minimum bias’ probability (12). The result of the calculations is shown in Fig. 1 with the dotted line. Calculations with $`T\overline{T}`$ and $`\mu _B\overline{\mu }_B`$, i.e. ignoring the fluctuations, were done for comparison (the thin dashed line). As is seen, the influence of fluctuations becomes essential only at $`E_T>100`$. This can be explained as follows. If $`E_T`$ is not very large, $`E_T<<qN_p(0)`$, the integration in (16) runs over the region where $`E_T/N_p>q`$, as well as over the one with $`E_T/N_p<q`$. The central point of the Gaussian (12), $`E_T/N_p=q`$, corresponds to $`T=\overline{T}`$ and $`\mu _B=\overline{\mu }_B`$. Deviations of $`T`$ and $`\mu _B`$ in both directions from their central values almost cancel each other. Therefore, the result does not differ essentially from that obtained at $`T\overline{T}`$ and $`\mu _B\overline{\mu }_B`$.

In contrast, when $`E_T`$ approaches the value $`qN_p(0)`$ or exceeds it, the cancellation does not take place any more. At large $`E_T`$, the region $`E_T/N_pq`$ falls into nonphysical domain $`N_p>N_p(0)`$. Only the tail of the Gaussian (12) with $`E_T/N_p>q`$ contributes to the integral (16). This corresponds to fluctuational tail with $`T>\overline{T}`$ and $`\mu _B<\overline{\mu }_B`$.

Although both the thermal ’density’ of $`J/\psi `$, $`n_{J/\psi }^{tot}(T,\mu _B)`$, as well as that of the open charm, $`n_O(T,\mu _B)`$, are very sensitive to the temperature, their ratio in the formula (2),

$$\frac{N_{J/\psi }^{tot}(T,\mu _B)}{(N_O(T,\mu _B)/2)^2}=\frac{1}{V}\frac{n_{J/\psi }^{tot}(T,\mu _B)}{(n_O(T,\mu _B)/2)^2},$$
(21)
is much less sensitive. The additional $`J/\psi `$ suppression takes place mostly due to increasing of the system volume $`V`$. Smaller $`\mu _B`$ leads to smaller baryon density $`n_B(T,\mu _B`$) in (5). At fixed $`N_p`$ this results in a larger volume. Decreasing the number of $`J/\psi `$ with increasing volume at fixed number of $`c\overline{c}`$ pairs is intuitively clear: the larger is the volume the smaller is the probability that $`c`$ and $`\overline{c}`$ meet each other and form a quarkonium state.

As is seen from figure 1, the suppression due to the $`T`$ and $`\mu _B`$ fluctuations is not strong enough to explain the experimental data at $`E_T>100`$ GeV. The reason is that we have assumed that $`P_{J/\psi }(E_T|b)P_{MB}(E_T|b)`$. In fact, however, production of a $`J/\psi `$ particle costs a fraction of energy of the colliding nuclei. Therefore, the multiplicity of ordinary hadrons (and consequently $`E_T`$) in a $`J/\psi `$ event is, in average, smaller than that in a ’minimum bias’ event (i.e. without registering $`J/\psi `$)§§§The ’minimum bias’ probability distribution $`P_{MB}(E_T|b)`$ is basically related to the events without $`J/\psi `$ or a Drell-Yan pair in the final state. Only a tiny fraction of Pb+Pb collisions results in production of charmonia or Drell-Yan pairs, therefore their influence on the shape of the minimum bias probability distribution is negligible. at the same impact parameter. Hence, as was pointed out in , the probability $`P_{J/\psi }(E_T|b)`$ is slightly different from $`P_{MB}(E_T|b)`$. Similarly to , we assume that the nucleon pair that have produced a $`J/\psi `$ meson does not contribute to the transverse energy. Therefore, the probability $`P_{J/\psi }(E_T|b)`$ is obtained from formula (12) by the replacement $`N_p(b)N_p(b)2`$. This induces a tiny (about $`0.5`$%) shift of the $`E_T`$ probability distribution.

The transverse energy distribution of the Drell-Yan event sample $`P_{DY}(E_T|b)`$ should also differ from $`P_{MB}(E_T|b)`$ because of the energy loss for the Drell-Yan pair production. However, in the large $`E_T`$ region, where the mentioned difference is important, our result is to be compared with the data obtained from the minimum bias analysis. This means that the so called ”theoretical Drell-Yan” instead of the physical one was used as a reference for $`J/\psi `$ suppression pattern. This ”theoretical Drell-Yan” was obtained neglecting the $`E_T`$ loss in the dimuon data sample. Therefore, to be consistent with the data presented by the NA50, we put $`P_{DY}(E_T|b)P_{MB}(E_T|b)`$.

The result of the calculations taking into account the modification of $`P_{J/\psi }(E_T|b)`$ (ignoring, for the moment, fluctuations of $`T`$ and $`\mu _B`$) is shown in figure 1 with the dash-dotted line. Again, a sizable effect (compare with the thin dashed line) is observed only at $`E_T>100`$. This is explained by the following properties of the Gaussian distribution. Let us consider a ratio of two Gaussians. One of them is shifted with respect to another one:

$$\frac{\mathrm{exp}\left\{\left[x(x_0\mathrm{\Delta })\right]^2\right\}}{\mathrm{exp}\left\{\left[xx_0\right]^2\right\}}=\mathrm{exp}\left\{2\mathrm{\Delta }\left[x(x_0\mathrm{\Delta })\right]\right\}.$$
(22)
If the shift $`\mathrm{\Delta }`$ is small, the ratio (22) is very close to $`1`$ near the central point $`x_0`$. But, due to exponential behavior, it differs essentially from unity at large distances from $`x_0`$, which corresponds to tails of the Gaussians. As was mentioned above, the shift of $`P_{J/\psi }(E_T|b)`$ with respect to $`P_{MB}(E_T|b)`$ is about $`0.5`$%. Therefore, when the main contribution to the integral (16) comes from the vicinity of the central value of the Gaussian (12), the effect is small. In contrast, at large $`E_T`$, when the integration runs only over the tail of the Gaussian, the effect becomes essential.

Taking into account both effects, $`T`$ and $`\mu _B`$ fluctuations and $`E_T`$ loss in the $`J/\psi `$ event sample, one finds excellent agreement with the fitted data including high $`E_T`$ ($`E_T>100`$ GeV) region (the thick solid line in figure 1). The free parameters $`\sigma _{c\overline{c}}^{NN}`$ and $`\eta `$ were fixed by fitting the model to the NA50 dataOnly the complete (i.e. including $`T`$ and $`\mu _B`$ fluctuations and $`E_T`$ loss) model has been fitted to the data. The calculations shown by the thin dashed, dotted and dash-dotted lines in figure 1 has been done with the parameter set fixed from the complete model fit..

As was mentioned above, the statistical coalescence model is not expected to describe small systems. This can be seen from the $`\psi ^{}`$ data . In the framework of this model the multiplicity of $`\psi ^{}`$ is given by formula (2) with the replacement $`N_{J/\psi }^{tot}N_\psi ^{}`$. Consequently, the $`\psi ^{}`$ to $`J/\psi `$ ratio as a function of the centrality should be constant and equal to its thermal equilibrium value. The experimental data  (see also a compilation in ) are consistent with this picture only at rather large numbers of participants . Therefore our fit is restricted to the region

$$E_T>27\text{ GeV,}$$
(23)
which corresponds to $`N_p100`$. Note that all NA50 data except the two leftmost points from the 1996 standard analysis set lie in the region (23).

To check the robustness of our model with respect to possible uncertainties we use two independent freeze-out parameter sets. The first one  has been obtained assuming strangeness and antistrangeness suppression by a factor $`\gamma _s`$:

$$\overline{T}=158\mathrm{MeV},\overline{\mu }_B=238\mathrm{MeV}\gamma _s=0.79.$$
(24)
In the second one  the complete strangeness equilibration has been supposed:

$$\overline{T}=168\mathrm{MeV},\overline{\mu }_B=266\mathrm{MeV}\gamma _s1.$$
(25)
The predicted effective $`c\overline{c}`$ production cross section as well as the fit quality appeared to be insensitive to the freeze-out parameters:

$$\sigma _{c\overline{c}}^{NN}=35.7\begin{array}{c}+10.8\hfill \\ 8.8\hfill \end{array}\mu \text{b},\chi ^2/\text{dof}=1.06$$
(26)
for both sets (24) and (25). The parameter $`\eta `$ only slightly depends on $`\overline{T}`$ and $`\overline{\mu }_B`$:

$$\eta =0.134\begin{array}{c}+0.061\hfill \\ 0.042\hfill \end{array}$$
(27)
for the set (24), and

$$\eta =0.126\begin{array}{c}+0.059\hfill \\ 0.039\hfill \end{array}$$
(28)
for the set (25).

As it was already pointed out in our previous papers , the model predicts a rather large enhancement of the total charm (by a factor of about $`3.5`$ within the rapidity window of the NA50 spectrometer). A direct measurement of the open charm in Pb+Pb collisions at the CERN SPS would allow for a test of this prediction.

Extrapolation of the model to low $`E_T`$ region reveals strong disagreement with the experimental data (about 3.5 standard deviations at the leftmost experimental point). This suggests that the charmonium production mechanism in central and peripheral Pb+Pb collisions may be of a very different nature. The ‘standard’ mechanism, formation of pre-resonant charmonium states in primary nucleon collisions and their subsequent suppression in nuclear medium, seems to dominate peripheral collisions. In contrast, the present statistical $`c\overline{c}`$ coalescence model assumes the formation of a deconfined quark-gluon medium in central collisions. This medium destroys all A different possibility has been considered in . It has been assumed that both direct and statistical coalescence mechanisms of $`J/\psi `$ production are present at all energies and at all collision centralities. No open charm enhancement has been permitted in . Therefore the contribution of statistical coalescence mechanism has been found to be small at the SPS. Then, the primordial $`J/\psi `$’s that survived suppression in the quark-gluon plasma dominate even in the most central Pb+Pb collisions. Such an approach is, however, hardly able to describe the centrality dependence of the $`J/\psi `$ suppression pattern. primary $`c\overline{c}`$ bound states. The final state charmonia are formed at a later stage of the reaction from $`c`$ and $`\overline{c}`$ quarks available in the medium. The distribution of $`c`$ and $`\overline{c}`$ over charmonia and open charm hadrons follows from the laws of statistical mechanics.

In conclusion, we have shown that the statistical $`c\overline{c}`$ coalescence model provides a quantitative description of the ‘anomalous $`J/\psi `$ suppression’ in Pb+Pb collisions at CERN SPS. In particular, the ‘second drop’ of the $`J/\psi `$ suppression pattern, at $`E_T100`$ GeV, can be naturally explained within the model as the result of two simultaneous effects: the fluctuations of the freeze-out thermodynamic parameters and transverse energy losses in the dimuon events with respect to the minimum bias ones. The model predicts a rather strong enhancement of the open charm. Therefore a direct measurement of the open charm would be a crucial test of the statistical $`c\overline{c}`$ coalescence model at the SPS.

###### Acknowledgements.

The authors are thankful to K.A. Bugaev, P. Bordalo, M. Gaździcki and L. McLerran for comments and discussions. We acknowledge the financial support of the Alexander von Humboldt Foundation, DFG, GSI and BMBF, Germany. The research described in this publication was made possible in part by Award # UP1-2119 of the U.S. Civilian Research and Development Foundation for the Independent States of the Former Soviet Union (CRDF) and INTAS grant 00-00366 and INTAS grant 00-00366.