text:
# Review on the inclusive rare decays 𝑩→𝑿_𝒔𝜸 and 𝑩→𝑿_𝒅𝜸 in the Standard Model

## 1 Introduction

In the Standard model (SM), rare $`B`$ decays like $`BX_s\gamma `$ or $`BX_s\mathrm{}^+\mathrm{}^{}`$ are induced by one-loop diagrams, where virtual $`W`$ bosons and up-type quarks are exchanged. In many extensions of the SM, there are additional contributions, where the SM particles in the loop are replaced by nonstandard ones, like charged Higgs bosons, gluinos, charginos etc. If the masses of these new particles are not heavier by many orders of magnitude than the heaviest SM particles, the new physics contributions to rare $`B`$ meson decays are expected to be generically large. The sensitivity for nonstandard effects implies the possibility for an indirect observation of new physics, or allows to put limits on the masses and coupling parameters of the new particles.

It is obvious that it is only possible to fully exploit the new physics potential of these decays when both, precise measurements and precise theoretical SM calculations exist.

In the following we mainly concentrate on the decays $`BX_s\gamma `$ and $`BX_d\gamma `$, while the rare semileptonic decay $`BX_s\mathrm{}^+\mathrm{}^{}`$ is reviewed at this conference by T. Hurth Hurth\_conf . There are experimental analyses of the branching ratio $`\mathrm{BR}(BX_s\gamma )`$ by CLEO CLEOold ; CLEOmed ; CLEOnew , ALEPH BRALEPH , BELLE BELLE01 , and BABAR BABAR02  as shown in fig. 1,

leading to the world average Jessop2003

$$\mathrm{BR}(BX_s\gamma )_{\mathrm{exp}}=(3.34\pm 0.38)\times 10^4.$$

In contrast to the exclusive rare decay $`BK^{}\gamma `$, the inclusive counterpart $`BX_s\gamma `$ is theoretically much cleaner as no specific model is needed to describe the hadronic final state. Indeed, nonperturbative effects in the inclusive decay mode are well under control due to the heavy quark expansion technique (HQE), which implies that the decay width $`\mathrm{\Gamma }(BX_s\gamma )`$ is well approximated by the partonic decay rate $`\mathrm{\Gamma }(bX_s\gamma )`$ which can be analyzed in renormalization group improved perturbation theory. The (nonperturbative) power corrections which scale like $`1/m_b^2`$ powermb  and $`1/m_c^2`$ powermc  were estimated to be well below $`10\%`$.

## 2 Theoretical framework

Short distance QCD effects enhance the partonic decay rate $`\mathrm{\Gamma }(bs\gamma )`$ by more than a factor of two. Analytically, these QCD corrections contain large logarithms of the form $`\alpha _s^n(m_b)\mathrm{ln}^m(m_b/M)`$, where $`M=m_t`$ or $`M=m_W`$ and $`mn`$ (with $`n=0,1,2,\mathrm{}`$). In order to get a reasonable prediction for the decay rate, it turns out that one has to resum both, the leading-log (LL) terms ($`m=n`$) as well as the next-to-leading-log (NLL) terms ($`m=n1`$).

To achieve the necessary resummations, one usually constructs in a first step an effective low-energy theory and then resums the large logarithms by renomalization group techniques. The low energy theory is obtained by integrating out the heavy particles which in the SM are the top quark and the $`W`$-boson. The resulting effective Hamiltonian relevant for $`bs\gamma `$ in the SM and many of its extensions reads

$$H_{\mathrm{eff}}(bs\gamma )=\frac{4G_F}{\sqrt{2}}\lambda _t\underset{i=1}{\overset{8}{}}C_i(\mu )O_i(\mu ),$$
(1)
where $`O_i(\mu )`$ are local operators consisting of light fields, $`C_i(\mu )`$ are the corresponding Wilson coefficients, which contain the complete top- and $`W`$\- mass dependence, and $`\lambda _t=V_{tb}V_{ts}^{}`$ with $`V_{ij}`$ being the CKM matrix elements. The CKM dependence globally factorizes, because we work in the approximation $`\lambda _u=0`$.

In the basis introduced by Misiak Mikolaj , the operators read

$`O_1`$ $`=`$ $`\left(\overline{s}_L\gamma ^\mu T^ac_L\right)\left(\overline{c}_L\gamma _\mu T^ab_L\right),`$
$`O_2`$ $`=`$ $`\left(\overline{s}_L\gamma ^\mu c_L\right)\left(\overline{c}_L\gamma _\mu b_L\right),`$
$`O_7`$ $`=`$ $`{\displaystyle \frac{e}{16\pi ^2}}m_b(\mu )\overline{s}\sigma ^{\mu \nu }RbF_{\mu \nu },`$
$`O_8`$ $`=`$ $`{\displaystyle \frac{g_s}{16\pi ^2}}m_b(\mu )\overline{s}\sigma ^{\mu \nu }RT^abG_{\mu \nu }^a.`$ (2)
As the Wilson coefficients of the QCD penguin operators $`O_3,\mathrm{},O_6`$ are small, we do not list them here.

A consistent calculation for $`bs\gamma `$ at NLL precision requires three steps:

* a matching calculation of the full standard model theory with the effective theory at the scale $`\mu =\mu _W`$ to order $`\alpha _s^1`$ for the Wilson coefficients, where $`\mu _W`$ denotes a scale of order $`M_W`$ or $`m_t`$;
* a renormalization group evolution of the Wilson coefficients from the matching scale $`\mu _W`$ down to the low scale $`\mu _b=𝒪(m_b)`$, using the anomalous-dimension matrix to order $`\alpha _s^2`$;
* a calculation of the matrix elements of the operators at the scale $`\mu =\mu _b`$ to order $`\alpha _s^1`$.

As all three steps are rather involved, a common effort of several independent groups was needed in order to calculate the NLL prediction for $`\mathrm{BR}(BX_s\gamma )`$ counterterm ; Adel ; GH ; INFRARED ; Mikolaj ; AG91 ; Pott ; GHW . For a detailed summary of the various steps and intermediate results, we refer to the recent review by T. Hurth Hurth\_rev . However, we would like to point out that the most difficult part, viz. the calculation of three-loop anomalous dimensions performed by Chetyrkin, Misiak and Münz in 1996 Mikolaj , was only confirmed very recently by Gambino, Gorbahn and Haisch Gambino:2003zm . Their paper also contains the three-loop mixing of the four-Fermi operators into $`O_9`$, which is important for the process $`BX_s\mathrm{}^+\mathrm{}^{}`$.

During the completion of the NLL QCD corrections, also calculations of electroweak corrections were started Czarnecki:1998tn ; Kagan:1999ym ; Baranowski:1999tq . At present, the corrections of order $`\alpha _{\mathrm{em}}\mathrm{ln}(\mu _b/M)`$ $`[\alpha _s\mathrm{ln}(\mu _b/M)]^n`$, as well as the subleading terms of order $`\alpha _{\mathrm{em}}[\alpha _s\mathrm{ln}(\mu _b/M)]^n`$ GH00  are systematically available.

## 3 NLL (and partial NNLL) results for BR$`\mathbf{(}B\mathbf{}X_s\gamma \mathbf{)}`$

Combining NLL QCD corrections with the electroweak corrections just mentioned and also including the $`1/m_b^2`$ powermb  and $`1/m_c^2`$ powermc  power corrections, the branching ratio reads

$$\mathrm{BR}(BX_s\gamma )=(3.32\pm 0.14\pm 0.26)\times 10^4,$$
(3)
where the first error reflects the dependence on the renormalization scale $`\mu _b`$ varied in the interval $`m_b/2\mu _b2m_b`$, while the second error reflects the error due to the uncertainties in the input parameters.

Among the input parameters the charm quark mass $`m_c`$ plays a crucial role. The charm quark mass dependence only enters the prediction for the decay width at the NLL level, more precisely through the $`𝒪(\alpha _s)`$ correction to the matrix elements $`s\gamma |O_{1,2}|b`$. Until recently, all authors used the pole mass value $`m_c^{\mathrm{pole}}`$ for the charm quark mass in numerical evaluations, leading to a branching ratio as specified in eq. (3).

In 2001, however, Gambino and Misiak Gambino:2001ew  pointed out that the MS-bar mass $`\overline{m}_c`$, normalized at $`\mu m_b/2`$, could be the better choice, because the charm quark appears as an off-shell particle in the loop involved in the above mentioned matrix element with a typical virtuality of $`m_b/2`$. Using this interpretation, $`\overline{m}_c/m_b=0.22\pm 0.04`$ is substantially smaller than the value $`m_c^{\mathrm{pole}}/m_b=0.29\pm 0.02`$ used in eq. (3), leading to a branching ratio of Gambino:2001ew ; Buras:2002tp

$$\mathrm{BR}(BX_s\gamma )=(3.70\pm 0.30)\times 10^4.$$
(4)
We would like to stress here that the above argument in favour of $`\overline{m}_c`$ is an intuitive one. Formally, the difference between using $`m_c^{\mathrm{pole}}`$ or $`\overline{m}_c`$ amounts to a NNLL effect at the level of the branching ratio. This means that a NNLL becomes necessary in order to unambiguously fix this issue.

Before sketching the NNLL program, we would like to stress that settling the $`m_c`$ issue is also important when extracting bounds on new physics, based on NLL calculations Ciuchini:1997xe ; BG98 . For example, in the type-II two-Higgs-doublet model, one obtains a bound from $`bs\gamma `$ on the charged Higgs boson mass of $`m_H>350`$ GeV ($`99\%`$ C.L.) when using $`\overline{m}_c`$. When using on the other hand $`m_c^{\mathrm{pole}}`$, the bound is $`m_H>280`$ GeV ($`99\%`$ C.L.) Gambino:2001ew .

Concerning the NNLL program, it is clear that in order to get a full NNLL QCD result for $`\mathrm{BR}(BX_s\gamma )`$, all the three steps listed above have to be improved by one order in $`\alpha _s`$. This means that three-loop matching calculations are needed, up to four-loop anomalous dimensions have to be worked out and up to three-loop calculations at the level of the matrix elements $`s\gamma |O_i(\mu _b)|b`$ have to be performed. Several groups have been formed in order to attack this ambitious goal. Recently, a calculation of the $`𝒪(\alpha _s^2n_f)`$ corrections to the matrix elements of the operators $`O_1`$, $`O_2`$, $`O_7`$ and $`O_8`$ was published Bieri:2003ue . Diagrammatically, these contributions are generated by inserting quark bubbles ($`n_f`$ denotes the number of light quarks) into the gluon propagators in the diagrams which are involved in the calculations at NLL order. We note that these contributions are not related to the definition problem of $`m_c`$. However, in many other cases they are sources of large corrections. E.g., in the semileptonic decay width $`\mathrm{\Gamma }(BX_c\mathrm{}\nu _{\mathrm{}})`$ these $`𝒪(\alpha _s^2n_f)`$ terms (after replacing $`n_f3\beta _0/2`$, according to the procedure of naive non-abelianization) incorporate more than $`80\%`$ of the complete $`𝒪(\alpha _s^2)`$ corrections Czarnecki:1998kt . The impact of the $`𝒪(\alpha _s^2n_f)`$ corrections to $`\mathrm{BR}(BX_s\gamma )`$ are shown in fig. 2. The dash-dotted curve shows the NLL prediction, while the solid curve incorporates in addition the $`𝒪(\alpha _s^2n_f)`$ terms (after the replacement $`n_f3\beta _0/2`$, according to naive non-abelianization). As one sees from the figure, the $`𝒪(\alpha _s^2n_f)`$ corrections seem to be small. Note, however, that this is a result of a relatively large accidental cancellation between corrections to $`O_2`$ and $`O_7`$. This point is illustrated by the long-dashed and short-dashed curve, which are obtained by switching off the $`𝒪(\alpha _s^2n_f)`$ corrections to $`O_7`$ and $`O_2`$, respectively.

## 4 Partially integrated BR and photon energy spectrum

The photon energy spectrum of the partonic decay $`bs\gamma `$ is a delta function, concentrated at $`(m_b/2)`$, when the $`b`$-quark decays at rest. This delta function gets smeared when considering the inclusive photon energy spectrum from a $`B`$ meson decay. There is a perturbative contribution to this smearing, induced by the Bremsstrahlung process $`bs\gamma g`$ AG91 ; Pott , as well as a nonperturbative one, which is due to the Fermi motion of the decaying $`b`$ quark in the $`B`$ meson.

For small photon energies, the $`\gamma `$-spectrum from $`BX_s\gamma `$ is completely overshadowed by background processes, like $`bc\overline{u}d\gamma `$ and $`bu\overline{u}d\gamma `$. This background falls off very rapidly with increasing photon energy, and becomes small for $`E_\gamma >2`$ GeV AG92 . This implies that only the partial branching ratio

$$\mathrm{BR}(BX_s\gamma )(E_\gamma ^{\mathrm{min}})=_{E_\gamma ^{\mathrm{min}}}^{E_\gamma ^{\mathrm{max}}}\frac{d\mathrm{BR}}{dE_\gamma }𝑑E_\gamma $$
(5)
can be directly measured, with $`E_\gamma ^{\mathrm{min}}=𝒪(2)`$ GeV.

Putting the energy cut at $`E_\gamma ^{\mathrm{min}}=2.0`$ GeV, CLEO used two methods to analyze their data on the photon energy spectrum in their most recent analysis: First, the Ali-Greub model AG91 ; AGmodel , based on the spectator model formulated in ref. ACCMM  and second, methods based on HQET Kagan:1999ym .

The spectator model contains two free parameters, viz. $`p_F`$, the average Fermi momentum of the $`b`$ quark in the $`B`$ meson and the mass of the spectator quark, $`m_{\mathrm{spec}}`$. Equivalently $`(p_F,m_b)`$ can be used as the free parameters, where $`m_b`$ is the average $`b`$ quark mass as defined in ref. AG91 ; AGmodel . In Fig. 3 a comparison between theory and experiment is shown. Using $`p_F=410`$ MeV and $`m_b=4.69`$ GeV the best fit is obtained. We would like to stress that similar values for these parameters are also obtained when fitting the lepton spectra in $`BX_c\mathrm{}\nu `$ and $`BX_u\mathrm{}\nu `$.

A modern way - based on first principles - implements the Fermi motion in the framework of the heavy-quark expansion. When probing the spectrum closer to the endpoint, the OPE breaks down, and the leading twist nonperturbative corrections must be resummed into the $`B`$ meson structure function $`f(k_+)`$ shape , where $`k_+`$ is the light-cone momentum of the $`b`$ quark in the $`B`$ meson. The physical spectrum is then obtained by the convolution

$$\frac{d\mathrm{\Gamma }}{dE_\gamma }=_{2E_\gamma m_b}^{\overline{\mathrm{\Lambda }}}𝑑k_+f(k_+)\frac{d\mathrm{\Gamma }_{\mathrm{part}}}{dE_\gamma }(m_b^{}),$$
(6)
where $`(d\mathrm{\Gamma }_{\mathrm{part}}/dE_\gamma )(m_b^{})`$ is the partonic differential rate, written as a function of the “effective mass” $`m_b^{}=m_b+k_+`$. The function $`f(k_+)`$ has support in the range $`\mathrm{}<k_+<\overline{\mathrm{\Lambda }}`$, where $`\overline{\mathrm{\Lambda }}=m_Bm_b`$ in the infinite mass limit. This implies that the addition of the structure function moves the partonic endpoint of the spectrum from $`m_b/2`$ to the physical endpoint $`m_B/2`$. While the shape of the function $`f(k_+)`$ is unknown, the first few moments $`A_n=𝑑k_+k_+^nf(k_+)`$ are known: $`A_0=1`$, $`A_1=0`$ and $`A_2=\lambda _1/3`$. As $`A_n`$ ($`n>2`$) are poorly known, several Ansätze were used for $`f(k_+)`$; e.g. Neubert and Kagan Kagan:1999ym  used $`f(k_+)=N(1x)^ae^{(1+a)x}`$, with $`x=k_+/\overline{\mathrm{\Lambda }}`$. Taking into account the constraints from $`A_0`$, $`A_1`$ and $`A_2`$, the independent parameters in this Ansatz can be chosen to be $`m_b`$ and $`\lambda _1`$. As shown in Kagan:1999ym , the uncertainty of $`m_b`$ dominates the error of the partial branching ratio.

In figure 4 the partial branching ratio is shown for the relevant range of $`m_b`$ as a function of $`E_\gamma ^{\mathrm{min}}`$, keeping $`\lambda _1/\overline{\mathrm{\Lambda }}^2`$ fixed. The data points show three CLEO measurements. In the oldest one the photon energy cut was put at $`E_\gamma ^{\mathrm{min}}=2.2`$ GeV, while in the most recent analysis this cut was lowered to 2.0 GeV, which is very important, because at 2.0 GeV the theoretical error on the partial branching ratio is considerably smaller, as seen from fig. 4.

To determine from the measurement of the partial branching ratio the full BR, one needs from theory the fraction $`R`$ of the $`BX_s\gamma `$ events with photon energies above $`E_\gamma ^{\mathrm{min}}`$. Based on Kagan-Neubert Kagan:1999ym , CLEO CLEOnew  obtained $`R=\left(0.915_{0.055}^{+0.027}\right)`$. A similar result is also obtained when using the spectator model.

It has been shown that up to corrections of $`𝒪(\mathrm{\Lambda }_{\mathrm{QCD}}/m_b)`$, the same shape function also describes $`BX_u\mathrm{}\nu `$ Neubert:1993um . This implies that the photon energy spectrum can be used to predict the fraction of $`BX_u\mathrm{}\nu `$ events with $`E_{\mathrm{lept}}>2.2`$ GeV, where leptons coming from $`BX_c\mathrm{}\nu `$ are absent for kinematical reasons. Taking into account perturbative and $`\mathrm{\Lambda }_{\mathrm{QCD}}/m_b`$ corrections Leibovich:1999xf ; Bauer:2002yu ; Neubert:2002yx , it is possible to extract $`V_{ub}`$ from a measurement of the $`BX_u\mathrm{}\nu `$ decay rate in the region above 2.2 GeV. CLEO used this strategy in ref. CLEO2002  to extract the CKM matrix element $`|V_{ub}|`$, obtaining $`|V_{ub}|=(4.08\pm 0.56_{\mathrm{exp}}\pm 0.29_{\mathrm{th}})\times 10^3`$.

## 5 $`𝑩\mathbf{}𝑿_𝒅𝜸`$ in the SM

The decay $`BX_d\gamma `$ can be treated in a similar way as $`BX_s\gamma `$ AAG98 . The only difference is that $`\lambda _u`$ for $`bd\gamma `$ is not small relative to $`\lambda _t`$ and $`\lambda _c`$; therefore, also the current-current operators $`O_1^u`$ and $`O_2^u`$, weighted by $`\lambda _u`$, contribute.

Unfortunately, these operators induce long-distance contributions to $`BX_d\gamma `$, which at present are not very well understood. To illustrate the problem, we first look at the corresponding charm quark loop, depicted in fig. 5. In this case, one can expand the loop function

$$_0^1𝑑x_0^{1x}𝑑y\frac{xy}{m_c^2\left[1\frac{k_g^2}{m_c^2}x(1x)2xy\frac{k_gk_\gamma }{m_c^2}\right]}$$

in powers of $`t=k_gk_\gamma /m_c^2`$, where $`k_g`$ and $`k_\gamma `$ denote the momentum of the gluon and the photon, respectively. This expansion generates the so-called Voloshin terms powermc , which in $`\mathrm{BR}(BX_s\gamma )`$ are a $`3\%`$ effect. Obviously, there is no such OPE in the case of the $`u`$quark loop. However, Buchalla, Isidori and Rey BIR  argued that an expansion in $`1/t`$ can be done, leading to non-local operators. From naive dimensional counting, the leading contribution is expected to be of order $`\mathrm{\Lambda }_{\mathrm{QCD}}/m_b`$.

In reference AAG98 , where NLL calculations for the process $`BX_d\gamma `$ were presented, the uncertainties due to the long-distance effects were absorbed into the theoretical error. Using $`\mu _b=2.5`$ GeV and the central values of the input parameters, the analysis in reference AAG98  yields a difference between the LL and NLL predictions for $`\mathrm{BR}(BX_d\gamma )`$ of $`10\%`$, increasing the branching ratio in the NLL case. For a fixed value of the CKM-Wolfenstein parameters $`\rho `$ and $`\eta `$, the theoretical uncertainty of the average branching ratio $`\mathrm{BR}(BX_d\gamma )`$ of the decay $`BX_d\gamma `$ and its charge conjugate $`\overline{B}\overline{X_d}\gamma `$ is: $`\mathrm{\Delta }\mathrm{BR}(BX_d\gamma )/\mathrm{BR}(BX_d\gamma )=\pm (610)\%`$. Of particular theoretical interest for constraining $`\rho `$ and $`\eta `$ is the ratio of the branching ratios, defined as

$$R(d\gamma /s\gamma )\frac{\mathrm{BR}(BX_d\gamma )}{\mathrm{BR}(BX_s\gamma )},$$
(7)
in which a good part of the theoretical uncertainties cancels. Varying the CKM-Wolfenstein parameters $`\rho `$ and $`\eta `$ in the range $`0.1\rho 0.4`$ and $`0.2\eta 0.46`$ and taking into account other parametric dependences, the results (without electroweak corrections) are

$`6.0\times 10^6`$ $``$ $`\mathrm{BR}(BX_d\gamma )2.6\times 10^5,`$
$`0.017`$ $``$ $`R(d\gamma /s\gamma )0.074.`$
Another observable, which is also sensitive to the CKM parameters $`\rho `$ and $`\eta `$, is the CP rate asymmetry $`a_{\mathrm{CP}}`$, defined as

$$a_{\mathrm{CP}}=\frac{\mathrm{\Gamma }(BX_d\gamma )\mathrm{\Gamma }(\overline{B}\overline{X_d}\gamma )}{\mathrm{\Gamma }(BX_d\gamma )+\mathrm{\Gamma }(\overline{B}\overline{X_d}\gamma )}.$$
(8)
Varying $`\rho `$ and $`\eta `$ in the range specified above, one gets $`7\%a_{\mathrm{CP}}35\%`$ AAG98 . We would like to point out that $`a_{\mathrm{CP}}`$ is at the moment only available to LL precision and therefore suffers from a relatively large renormalization scale dependence.

In summary, this decay mode is very challenging, both in theory and experiment: On the theory side more work is needed concerning the nonperturbative contributions associated with the $`u`$quark loop, while on the experimental side the observation of this decay needs high statistics and a very good discrimination between pions and kaons.

###### Acknowledgements.

This work is partially supported by the Swiss National Foundation and by RTN, BBW-Contract N0. 01.0357 and EC-Contract HPRN-CT-2002-00311 (EURIDICE).


format: text

text:
# Small rare gas clusters in XUV laser pulses

## 1 Introduction

Clusters can absorb laser energy more efficiently than both gas targets and solids. This is because the cluster size is smaller than the skin depth so that all atoms experience the same laser field, despite the locally high particle density. Moreover, unlike solid targets there is no cold bulk that serves as a reservoir of cold electrons. Consequently, the laser-heated electrons cannot as easily escape from the cluster as they can from a solid surface. The electrons that absorb enough energy to leave the cluster as a whole contribute to the so-called “outer ionization.” Other electrons are removed from their parent atom or ion (“inner ionization”) but remain trapped by the positive background of the cluster. These electrons and the ions form a nanoplasma that hydrodynamically expands or Coulomb-explodes. In experiments at wavelengths $`248`$ nm highly energetic electrons shao , ions ditmirenature ; ditmirePRA ; springate , photons mcpherson ; ditmireJPB ; teravet , and neutrons originating from nuclear fusion zweibackDD  were observed (see Refs. review ; reviewII  for reviews).

Recently, Xe cluster experiments in the XUV regime have been performed at the DESY free-electron laser (FEL), Hamburg, Germany wabnitz . Surprisingly, for laser intensities of a few times $`10^{13}`$, $`100`$ fs pulse duration, and $`98`$ nm wavelength multiply charged Xe ions were observed in the laser-cluster experiments even for small clusters whereas isolated Xe atoms loose only one electron under such conditions.

The goal of this work is to shed some light on the mechanism underlying the generation of higher charge states inside clusters for laser parameters close to the FEL experiment. For that purpose molecular dynamics (MD) simulations have been performed. The inner-atomic dynamics is treated explicitly. The classical multi-electron atoms are rendered stable by introducing a momentum-dependent potential. Mechanisms contributing to inner ionization such as field ionization, collisional ionization, and even non-sequential ionization and “shake-off” processes are all self-consistently incorporated on a classical level.

The paper is organized as follows. In Sec. 2 the numerical model is introduced. In Sec. 3 the numerical results are presented. Finally, we conclude in Sec. 4.

Atomic units (a.u.) are used throughout unless noted otherwise.

## 2 MD model of rare gas clusters

The rare gas cluster consisting of $`N_a`$ atoms at the positions $`\stackrel{}{R}_i`$, $`1iN_a`$, and $`Z`$ “active” electrons per atom at the positions $`\stackrel{}{r}_j`$, $`1jZN_a`$ in a laser field $`\stackrel{}{E}(t)`$ is modeled by the Hamiltonian

$`(\stackrel{}{R},\stackrel{}{P};\stackrel{}{r},\stackrel{}{p};t)={\displaystyle \underset{i=1}{\overset{N_a}{}}}{\displaystyle \frac{\stackrel{}{P}_i^2}{2M}}+{\displaystyle \underset{j=1}{\overset{ZN_a}{}}}{\displaystyle \frac{\stackrel{}{p}_j^2}{2}}`$
$`+{\displaystyle \underset{i=1}{\overset{N_a}{}}}{\displaystyle \underset{j=1}{\overset{ZN_a}{}}}\left(V_H(\stackrel{~}{r}_{ij},\stackrel{~}{p}_{ij}){\displaystyle \frac{Z}{|\stackrel{}{R}_i\stackrel{}{r}_j|}}\right)`$
$`+{\displaystyle \underset{i=1}{\overset{N_a}{}}}{\displaystyle \underset{k=1}{\overset{i1}{}}}\left(V_{aa}(|\stackrel{}{R}_i\stackrel{}{R}_k|)+{\displaystyle \frac{Z^2}{|\stackrel{}{R}_i\stackrel{}{R}_k|}}\right)`$
$`+{\displaystyle \underset{j=1}{\overset{ZN_a}{}}}{\displaystyle \underset{l=1}{\overset{j1}{}}}{\displaystyle \frac{1}{|\stackrel{}{r}_j\stackrel{}{r}_l|}}+\stackrel{}{E}(t)\left({\displaystyle \underset{j=1}{\overset{ZN_a}{}}}\stackrel{}{r}_jZ{\displaystyle \underset{i=1}{\overset{N_a}{}}}\stackrel{}{R}_i\right).`$
Apart from the usual terms describing the kinetic energy, the Coulomb interactions, and the interaction with the laser field in dipole approximation, there are the additional potentials $`V_{aa}`$ and $`V_H`$ accounting for the Lennard-Jones interaction between neutral atoms and the Heisenberg uncertainty principle for the electrons inside the atoms, respectively. The purpose of introducing these potentials will become clear when we now describe how the initial, unperturbed cluster configuration is built up.

(I) First, we seek an energetically optimal configuration for the neutral cluster without laser field. To that end we set $`Z=0`$ and

$$V_{aa}(R_{ik})=D\left[\left(\frac{a}{R_{ik}}\right)^{12}\zeta \left(\frac{a}{R_{ik}}\right)^6\right]$$
(2)
where $`\stackrel{}{R}_{ik}=\stackrel{}{R}_i\stackrel{}{R}_k`$, $`R_{ik}=|\stackrel{}{R}_{ik}|`$, and $`\zeta `$ is commonly set to either unity or two in the chemical physics literature. The potential minimum of $`V_{aa}`$ is at $`R_{\mathrm{min}}=(2/\zeta )^{1/6}a`$ and should be chosen close to the known nearest neighbor distance of the cluster under consideration. By varying the parameter $`D`$ the overall strength of the Lennard-Jones potential $`V_{aa}`$ can be adjusted.

Starting from a random atom distribution with the nearest neighbor distances $`>R_{\mathrm{min}}`$ a local minimum in the energy landscape is obtained by propagating the atoms according the equations of motion

$`\dot{\stackrel{}{R}}_i`$ $`=`$ $`\stackrel{}{P}_i/M,`$ (3)
$`\dot{\stackrel{}{P}}_i`$ $`=`$ $`{\displaystyle \frac{}{\stackrel{}{R}_i}}\nu \stackrel{}{P}_i`$ (4)
with a non-vanishing but small friction $`\nu `$ and the nuclear mass $`M`$ set to a small value so that relaxation occurs on an acceptable time scale. For small clusters (say, $`N_a<10`$) it is quite likely that also the global minimum is found in this way as long as one starts with a reasonable guess for the cluster structure. Finding the global minimum for bigger clusters is far from trivial. Fortunately, the structures of Lennard-Jones clusters with $`N_a<150`$ are available in the literature wales ; waleswww .

(II) Secondly, we seek an electronic configuration for the single atom with $`Z`$ active electrons and call this the “mother configuration.” Since “classical atoms” with more than one electron are generally unstable, a momentum-dependent potential

$$𝒱(r,p,\xi ,\alpha ,\mu )=\frac{\xi ^2}{4\alpha r^2\mu }\mathrm{exp}\left\{\alpha \left[1\left(\frac{rp}{\xi }\right)^4\right]\right\}$$
(5)
is introduced kw  that enforces approximately the Heisenberg uncertainty relation when applied in the form

$$V_H(\stackrel{~}{r}_{ij},\stackrel{~}{p}_{ij})=𝒱(\stackrel{~}{r}_{ij},\stackrel{~}{p}_{ij},\xi _H,\alpha _H,\mu _{ei})$$
(6)
where

$$\stackrel{~}{r}_{ij}=|\stackrel{}{r}_i\stackrel{}{R}_j|,\stackrel{~}{p}_{ij}=\left|\frac{M\stackrel{}{p}_i\stackrel{}{P}_j}{1+M}\right|.$$
(7)
Here, $`\mu _{ei}=M/(1+M)1`$ is the reduced mass and $`\stackrel{~}{r}_{ij}`$ and $`\stackrel{~}{p}_{ij}`$ are the absolute values of relative distance and momentum, respectively.

The “hardness parameter” $`\alpha _H`$ governs how strictly the uncertainty relation $`\stackrel{~}{r}_{ij}\stackrel{~}{p}_{ij}\xi _H`$ is fulfilled. Big values of $`\alpha _H`$ enforce it severely but also make the differential equations of motion stiff (which is numerically unfavorable). We have chosen $`\alpha _H=2`$. The parameter $`\xi _H`$ may be adjusted in such a way that the known essential ground state features of the system under study (e.g., the total energy or the ionization potentials) are properly mimicked.

Although we neglect the spin in the current study it is worth noticing that the Pauli principle could be modeled in a similar way so that an atomic shell structure is obtained cohen ; cohen2000 .

By solving the Hamilton equations of motion corresponding to (2) with a non-vanishing but small friction $`\nu `$, the “mother configuration” of atoms or molecules with only a few electrons can be easily found. For “bigger” atoms more advanced minimization routines have to be employed cohen ; cohen2000 . Note that contrary to other classical trajectory Monte Carlo (CTMC) methods where the classical atom is modeled by an ensemble of electrons moving on Kepler orbits, in our case the groundstate is stationary, that is $`\dot{\stackrel{}{r}}=\dot{\stackrel{}{p}}=0`$ but, owing to the Heisenberg potential, $`\stackrel{}{p}0`$.

(III) Finally, the entire cluster is build up by taking $`N_a`$ randomly rotated electronic mother configurations and attaching them to the naked ions. The ion positions are known from the first step. Since now all electrons “see” each other (as well as all the other ions), another propagation with non-vanishing friction $`\nu `$ is required. During this step the electrons orient themselves in an energetically favorable way while the ions hardly move due to their huge mass $`M`$. Note that the cluster structure obtained in this way is not completely self-consistent because the ions are still sitting at the positions determined in step (I) where just the Lennard-Jones potential was effective. In a fully self-consistent procedure the semi-classical cluster should be assembled without a Lennard-Jones potential at all. However, this is numerically much more demanding and may lead to an ion distribution which is physically less reasonable than the Lennard-Jones structure after step (I).

We call the cluster configuration obtained after step (III) the “cluster mother configuration.” In order to obtain meaningful results concerning the ionization dynamics of the isolated atom or the cluster in a laser pulse, an ensemble of atoms and clusters had to be simulated. The members of such an ensemble were constructed by randomly rotating the cluster mother configuration. While this did not change the total energy of the cluster or single atom, of course, the orientation with respect to the laser polarization axis changed, and, hence, so did the ionization dynamics. As expected, the ensemble-averaged entities converged more rapidly with respect to the time step for the integration of the equations of motion (a fourth order Runge-Kutta scheme was used) as the individual particle trajectories did. Convergence of the latter may be accelerated by reducing $`\alpha _H`$ or introducing a “soft-core” smoothing parameter for the electron-electron and electron-ion interaction, that is, e.g., replacing $`|\stackrel{}{r}_j\stackrel{}{r}_l|^1`$ by $`[(\stackrel{}{r}_j\stackrel{}{r}_l)^2+a_{ee}]^{1/2}`$ since all this loosens the requirement of small time steps during close encounters of the particles.

## 3 Results

We simulated small Xe clusters in intense laser fields. Each semi-classical Xe atom had $`Z=3`$ active electrons that were supposed to mimic three of the six 5p electrons of Xe. Treating more than three electrons per atom is not necessary for the modest laser intensities we consider in this work.

The other simulation parameters were chosen $`M=131\times 1836`$, $`D=0.04`$, $`\zeta =1`$, $`R_{\mathrm{min}}=8.0`$, $`\alpha _H=2`$, $`\xi _H=1.75`$. The “mother configuration,” i.e., a single Xe atom, consists of the three active electrons sitting on the corners of an equilateral triangle and the ion sitting in the center. The ion-electron distance was $`1.6`$ for all three electrons. The total energy of the configuration was $`_1=2.30`$ while the three ionization potentials were $`I_1=0.32`$, $`I_2=0.80`$, and $`I_3=1.18`$. The first three ionization potentials for the real Xe atom are $`0.45`$, $`0.77`$, and $`1.18`$. The parameter $`\xi _H=1.75`$ was chosen to optimize $`I_2`$ and $`I_3`$ since we are interested in multiple ionization and are particularly interested in the intensity regime where the removal of the first electron is given for granted for both the isolated atom and the atoms inside the cluster but higher charge states only occur in the cluster. Spin was neglected for if we also introduced a Pauli-blocking potential $`V_P`$ of the form (5) into (2) we would have modeled the Li atom (with its closed s-shell plus a loosely bound valence electron) rather than three of the six 5p electrons of Xe.

### 3.1 Isolated atom

In Fig. 1 the charge state of the isolated Xe atom is shown vs the peak laser intensity of a $`T=1764`$ ($`42`$ fs) laser pulse. The envelope of the electric field was linearly ramped up (and down) over $`T_r=331`$ ($`8`$ fs) while it was constant for the rest of the pulse. The two laser frequencies $`\omega _l=0.057`$ (corresponding to $`800`$ nm) and $`\omega _h=8\omega _l`$ (corresponding to $`100`$ nm) were used.

As expected, for the same pulse duration and laser intensity the higher laser frequency ionizes the isolated atom more efficiently than the lower one. For laser intensities $`I<0.001`$ ($`3.5\times 10^{13}`$Wcm-2) it is unlikely to observe charge states higher than the singly ionized Xe. This is in agreement with the results of Wabnitz et al. wabnitz  where at an intensity $`2\times 10^{13}`$Wcm-2 for $`12.7`$-eV photons (corresponding to $`8.2\omega _l`$) and $`100`$ fs pulse durations also higher charge states than Xe+ were found to be absent for isolated atoms.

### 3.2 Small cluster

In Fig. 2 the unperturbed Xe model cluster consisting of $`N_a=27`$ triply charged ions (black) plus the $`ZN_a=3\times 27`$ electrons (gray) is shown. The ion structure is in agreement with Refs. wales ; waleswww . The electrons roughly maintain their triangular setup with respect to their parent ion. However, the triangles are in general not equilateral anymore due to the interaction with the neighboring atoms.

The total energy of the unperturbed cluster amounts to $`_{27}=63.3<N_a_1=62.1`$, i.e., the calculated cluster configuration is indeed energetically more favorable than $`N_a`$ isolated atoms. The binding energy per atom is $`0.04`$ which is a reasonable value for rare gas clusters.

From Fig. 1 one can infer that for laser intensities $`I`$ less than $`0.001`$ no charge states $`>1`$ are expected in the case of isolated Xe atoms. Hence it is interesting to investigate whether higher charge states are observed for clusters. In Fig. 3 the results for the two laser frequencies $`\omega _l=0.057`$, $`\omega _h=8\omega _l`$ and the two laser intensities $`I=4\times 10^4`$ ($`\widehat{E}=0.02`$), $`I=1.6\times 10^3`$ ($`\widehat{E}=0.04`$) are collected for both the isolated atom and the cluster consisting of $`N_a=27`$ atoms. The laser pulse duration and shape was the same as in Sec. 3.1.

The inner ionization of the cluster at the high frequency $`\omega _h`$ (indicated by “inner (hf)” in the plots) is relatively high for both laser intensities. The inner ionization was calculated by counting the number of electrons that are more than $`8`$ a.u. away from their initial position and therefore must have left their parent ion. If recombination plays an important role inner ionization is overestimated by this method. However, by looking at the particle trajectories we inferred that it is unlikely that an electron which left its parent ion is permanently trapped by another ion.

Inner ionization is rather low for $`I_l`$ and $`\omega _l`$. This is because an important prerequisite is missing for these parameters, namely the efficient removal of the outermost electron. It is clear that the very early inner ionization dynamics of clusters is the same as in the isolated atom case since there is not yet a strong influence of the other ions at that stage. In fact, the single atom ionization is small for $`I_l`$ and $`\omega _l`$ (see Fig. 3a, dashed curve close to the bottom) so that the inner ionization of the cluster also remains quite modest. This is different for the higher laser intensity $`I_h`$ (Fig. 3b) where the isolated atom looses the outermost electron. Consequently, in the case of a cluster a space charge builds up that yields further inner ionization (“ionization ignition,” rosepetru ).

Outer ionization (calculated by counting the number of electrons which are farther away than $`r=30`$) is always smaller than the corresponding ionization of the isolated atom for the parameters in Fig. 3 (see dotted curves). This means that a significant fraction of electrons that were removed from their parent ions do not make it to leave the cluster as a whole. This can be attributed to the low quiver energies of the electrons inside the cluster so that no collective electron dynamics sets in that could help to absorb additional laser energy and to overcome the strong backholding field of the cluster ions bauer . Note that the ponderomotive energies for the parameters in Fig. 3 all lie between $`4.8\times 10^4`$ and $`0.12`$ and are thus neither large compared to the photon energy nor large compared to the binding energies.

In Fig. 4 the total energy (per atom) $`(t)`$ during the course of the laser pulse is plotted vs time. The laser and cluster parameters were the same as in Fig. 3. The net absorbed energy (per atom) $`|(T)(0)|`$ is for both frequencies and laser intensities greater for the cluster than for the isolated atom. As already mentioned in the discussion of the ionization degrees shown in Fig. 3, the absorbed energy is low for $`I_l`$ and $`\omega _l`$ owing to the rather small probability for single ionization. For the higher intensity $`I_h`$ instead, slightly more energy is absorbed at the low frequency $`\omega _l`$. This matches with the increased outer ionization at low frequencies (cf. Fig. 3b). Moreover, the kinetic energy of the freed electrons depends on the laser frequency as well.

In order to understand the absorption mechanism at short wavelengths it is worth to look at sample electron trajectories. In Fig. 5 one sees how electrons that at the end of the laser pulse are more than $`30`$ a.u. away from the cluster center (and therefore contribute very likely to outer ionization) manage to escape from the cluster. As long as the distance from the initial position, $`|\stackrel{}{r}\stackrel{}{r}_0|`$, is less than $`8`$ a.u., the electrons can be considered bound with respect to their parent ion. Zig-zag motion between distances from $`8`$ to $`23`$ is due to collisions with other ions while electrons at distances $`>23`$ left the entire cluster.

For long wavelength and low intensity (a) one can see that the few electrons that made it to leave the cluster did so without collisions. This is because of the low probability for the removal of even the outermost electrons: most of the atoms remain neutral so that it is unlikely that the few freed electrons encounter an ion.

At higher frequency (b,d) it is seen from the erratic motion for $`8<|\stackrel{}{r}\stackrel{}{r}_0|<23`$ that the electrons experience several collisions with the cluster ions before they leave the cluster. Moreover, the electrons are temporarily trapped by other cluster ions.

In plot (c), for twice the laser intensity and the low frequency, collisions are observed as well. This is because both the quiver velocity and the thermal velocity are small so that electron-ion collisions are likely (see, e.g., mulser ). The outer electrons are removed efficiently at that laser intensity so that the ionization ignition mechanism starts (see the curve for inner ionization in Fig. 3b). This indicates the eminent importance of the first freed electrons that trigger the entire cluster ionization process.

From the trajectories we deduce that the efficient absorption of laser energy, as it is seen in Fig. 4, is due to collisional heating or, in other words, inverse bremsstrahlung. Similar conclusions were drawn recently by other authors using different methods preprsantra ; preprrost . The increased inner ionization is due to (i) the reduced ionization potential because of the neighboring ions, (ii) the field generated inside the cluster by all other ions, and (iii) collisions. All three mechanisms are included self-consistently in our numerical treatment which, unfortunately, implies that it is difficult to separate the different contributions and to identify their relative importance.

Although high charge states are generated inside the cluster during the course of the laser pulse it is not yet clear which charge states would arrive at the time-of-flight (TOF) detector in a real experiment. Upon the expansion of the cluster the electrons cool and some of them recombine, leading to a lower average charge state. Therefore we simulated the cluster disintegration as well. Since the cluster expansion happens on a picosecond time scale we reduced the ion mass to $`M=50`$ in order to make the simulation feasible remark . While this mass reduction changes the time scale of the cluster expansion it should not affect too strongly the final energy and charge state distribution as long as $`m/M1`$.

The result is shown in Fig. 6 for the high frequency $`\omega _h`$ and the two intensities $`I_l`$ and $`I_h`$. In Fig. 6a it is seen that despite the plasma cooling a significant amount of Xe2+ ions survive the cluster expansion. In the case of the higher intensity $`I_h`$ the doubly charged Xe ion is the most likely charge state and even triply ionized Xe atoms can be observed.

Comparison of Fig. 6a with Fig. 1 in wabnitz  (TOF spectrum for $`N=2`$$`20`$) suggests reasonable agreement. However, a precise comparison with the experiment wabnitz  is not straightforward because of the not exactly known experimental laser intensity, the varying sensitivity of the TOF spectrometer with respect to different charge states, and focus effects.

## 4 Conclusions

Molecular dynamics was used to investigate small rare gas clusters in short XUV and 800 nm laser pulses. Inner ionization was treated explicitly. In order to render classical, stable multi-electron atoms possible, a momentum-dependent potential that accounts for the Heisenberg uncertainty principle was introduced. Higher charge states in clusters than in the isolated atoms were observed, in accordance with experimental results. The mechanism underlying the efficient absorption of XUV laser energy was found to be inverse bremsstrahlung.

The study of bigger clusters consisting of several hundreds or even thousands of atoms would be desirable. However, in order to keep the numerical simulations feasible one then has to give up the explicit treatment of the inneratomic dynamics.

## Acknowledgments

The author thanks A. Macchi for proofreading and valuable suggestions. This work was supported by the Deutsche Forschungsgemeinschaft and by the INFM Advanced Research Project CLUSTERS. The permission to run our codes on the Linux cluster at PC2 in Paderborn, Germany, is gratefully acknowledged.


text:
# Experimental Signals of Phase Transition

## 1 INTRODUCTION

In the last 20 years and especially in the most recent ones, considerable progress has been achieved both theoretically and experimentally in the investigation of nuclear reaction dynamics and thermodynamics in the Fermi energy regime .

From a theoretical point of view, strong efforts have been devoted to the understanding of the nuclear equation of state (EoS) either with transport theories or statistical approaches.

From an experimental point of view new generation $`4\pi `$ detectors have been developed and are now operating at different accelerator facilities (Dubna, GANIL, GSI, LNL, LNS, MSU, Texas A-M). They are producing a huge amount of exclusive data and new kind of analyses. Very rich information has already been extracted from experimental studies on intermediate energy heavy ion collisions. Experiments have shown that the final state can be constrained to select the dynamics of the collision and isolate events that populate states closely compatible with equilibrium. Several investigations have demonstrated that excited nuclear systems produced in such collisions undergo in a short time scale (100 fm/c) bulk multifragmentation characterized by final states containing several Intermediate Mass Fragments (IMF,$`Z`$ 3).

A considerable progress has been accomplished on the theoretical as well as on the experimental side in order to define and collect a converging ensemble of signals connecting multifragmentation to the nuclear liquid-gas phase transition  and locating it in the phase diagram of nuclear matter . The opening of the high fragment multiplicity channel, the onset of collective expansion, the flattening of the caloric curves, the fossile signal of spinodal decomposition, a negative branch of the heat capacity, the bimodal distribution of exclusive observables and the finite size and Fisher law scalings have been observed and tentatively related to the Equation of State of the nuclear matter. The possibility of new radioactive beams (RIB) facilities is now prompting exciting theoretical advances of the isospin aspects of the EoS, like the density dependence of the symmetry energy and the modifications to the spinodal instability . Experimentally, first studies on the dependence of the isoscaling parameters on the isospin of the decaying system  already started exploiting stable beams. These works seem to indicate an isospin distillation in asymmetric systems.

All these signals can be considered as circumstantial evidences of a phase transition, but some of them are still controversial and need to be further experimentally investigated before the phase transition can be definitely assessed. This work is a contribution to this aim, but only next generation experiments will allow to reach a comprehensive and detailed understanding of the phase transition and the nuclear EoS.

## 2 THERMODYNAMICS OF NUCLEAR SYSTEMS

### 2.1 EQUILIBRIUM PARTITIONS

In order to perform thermodynamical analyses, one has to collect a data sample which corresponds as closely as possible to an homogeneous population of the phase space. Data must be selected such as to isolate a portion of the cross section where the entire system (or the quasi-projectile) properties keep a negligeable memory of the entrance channel dynamics. This can be experimentally verified checking that for a given source the fragmentation pattern is determined by the size, charge, energy and average freeze-out volume solely, independent of the way the source has been formed, e.g. different impact parameters. In this case the thermodynamics we can access is a microcanonical thermodynamics with energy, number of protons and neutrons, and average volume as state variables. Indeed the excitation energy can be measured on an event-by-event basis by calorimetric techniques. For any shape of the excitation energy distribution the events can thus be sorted in constant energy bins, i.e. in microcanonical ensembles.

In the following we present thermodynamical studies performed on quasi-projectile events from peripheral 35 A MeV $`Au+Au`$ collisions  and central events from 25 A MeV $`Au+C`$, 25 and 35 A MeV $`Au+Cu`$ and 35 A MeV $`Au+Au`$ collisions , measured at the K1200-NSCL Cyclotron of the Michigan State University with the MULTICS-MINIBALL apparatus.

Single source almost complete events have been selected with a constant value for the collected charge11190% of the total charge for central collisions and 90% of the projectile charge for peripheral events in each energy bin . The possible pollution from other sources has been minimized for central collisions through a shape analysis and in the case of $`Au`$ quasi-projectile by substituting the backward light particle emission by the symmetric of the forward emission in the quasi-projectile reference frame. The observed event isotropy (Fig. 1) indicates that the directed flow component coming from a memory of the entrance channel is negligible.

The close similarity between statistical models  and data , together with the isotropy of fragment emission, already suggest that these sets of data are close to a statistical equilibrium.

To progress on this point, we can also compare different data sets. Central events correspond to a narrow distribution of the excitation energy, while the quasi-projectile data supply a widely spread excitation function. Therefore it is possible to find for these data sets common values for the energy deposited in the source and to compare charge distributions and charge partitions. In the left four panels of Fig. 2 the charge distribution of all the reaction products, but the largest one, masured for central $`Au+C`$ events at 25 A MeV, $`Au+Cu`$ at 25, 35 A MeV and $`Au+Au`$ at 35 A MeV (top left, top right, bottom left, bottom right panels, respectively) are compared to the quasi-projectile distributions at the same excitation energy. The calorimetric values of the energy deposited in the systems cover a wide range. Indeed their values are 1.6, 3.1, 4.7 and 7.4 A MeV, respectively. Moreover, the comparison among the charge distributions of the heaviest fragment in each event (Fig. 2) allows to check the scaling of higher moments like the variances, even in more detail.

A very good scaling behaviour is also obtained by comparing the presented sets of data and the charge distributions measured by the FASA collaboration for the reactions $`p,\alpha +Au`$ at 8.1, 4 and 14.6 GeV  and by the Indra collaboration for $`Xe+Sn`$ reaction at 32 A MeV incident energy .

The remarkable scaling between different sets of data means that these data samples can be analyzed, at least in a good first approximation, within statistical microcanonical methods.

### 2.2 PSEUDO-CRITICAL BEHAVIOUR IN NUCLEI

Since the early 80’s, size distributions have been fitted with power laws , and more sophisticated critical analyses have been performed following theoretical concepts coming from percolation theory. More recently, an astonishing good scaling behaviour has been observed in the EOS data  and tentatively associated to the critical point of the nuclear liquid-gas phase transition expected to occur in nuclear matter in the framework of the Fisher droplet model . The debate on the order of the transition has been further animated by a very recent analysis of the EOS and Isis data  which shows a high quality scaling of the fragment size distribution over a wide range of charges and deposited energies with an ansatz for the scaling function taken from the Fisher droplet model. The Fermi gas ”critical” temperature (about 8 MeV) extracted in these papers is identified as the temperature of the thermodynamical critical point and the whole coexistence line of the nuclear phase diagram is reconstructed under the assumption that the Fisher model gives a good description of the multifragmentation phenomenon . In this interpretation, multifragmentation would correspond to the critical point of the nuclear matter EoS (with a lower temperature due to finite size and Coulomb effects), i.e. to a second order phase transition.

However, the experimental observation of a flattening of the caloric curve  and recent studies  in the framework of the SMM model  (in which the resulting critical temperature is about 20 MeV) point rather to a first order phase transition, and this is also suggested by other thermodynamical statistical multifragmentation models .

One may also wonder what physical meaning can be attributed to a Fermi gas estimator of the temperature at excitation energies of several A MeV. Indeed for the $`Au`$ quasi-projectile data we are discussing, the temperature, evaluated from an isotope and a kinetic energy thermometer , strongly deviates from the Fermi gas ansatz (left part of Fig. 3).

Surprisingly enough, a comparable quality of scaling (Fig. 3) and a consistent set of critical exponents can also be extracted from the $`Au`$ QP data samples, for which the heat capacity evaluation points to a first order phase transition , and the quality of the scaling does not depend the ansatz adopted for the $`T`$ parameter. The same ansatz for the scaling function applied to central events  gives similar exponents and points to the same critical-like energy.

However, critical exponents and scale invariance are compatible with many different physical phenomena and are not necessarily linked to a thermodynamic second order phase transition . In particular the observed signals of critical behavior can be compatible with a first order phase transition, since in different statistical models size distributions that mimic a scale invariant behavior are observed inside the coexistence zone of small systems . This means that scaling per sé does not demonstrate the existence of a phase transition, and even less defines its order or allows to localize the system on the phase diagram.

On the other hand, the fact that all the analyzed reactions behave as a universal multifragmentation process independent of the entrance channel and directly correlated with the available energy only, is a strong indication of a microcanonical equilibrium.

## 3 MULTIFRAGMENTATION AND SIGNALS OF PHASE TRANSITION

It has been shown  that for a given total energy the average partial energy stored in a subsystem of the microcanonical ensemble is a good thermometer while the fluctuations associated to the partial energy can be used to evaluate the heat capacity. An example of such a decomposition is given by the kinetic $`E_k`$ and the interaction $`E_I`$ energies. In particular first order phase transitions are marked by singularities and negative heat capacities , corresponding to fluctuations anomalously larger than the canonical expectation. If the system is in statistical equilibrium, a measurement of anomalous fluctuations at a given energy is an unambiguous proof of a thermal first order phase transition.

Experimentally, the total energy $`E^{}`$ deposited in the system can be evaluated event by event by calorimetry. If one is also able to reasonably estimate for multifragmentation data the relation between the measured charge and mass of the reaction products and the value of the freeze-out volume, one can obtain the event-by-event Q-value and the Coulomb energy, and thus the interaction energy, according to:

$$E_I=\underset{i=1}{\overset{M}{}}m_i+E_{coul}^{FO}m_0(i=1,M)$$
(1)
where M is the total multiplicity, $`m_i`$ ($`m_0`$) are the mass excess of the primary products (of the source). The interaction energy fluctuation can then be studied as a function of the total energy and the heat capacity can be evaluated according to:

$$C=\frac{C_k^2T^2}{C_kT^2\sigma _k^2}$$
(2)
where $`\sigma _k^2=\sigma _I^2`$ is the fluctuation of the interaction energy $`E_I`$ from Eq.(1), $`T`$ is the temperature, and $`C_k`$ is the kinetic heat capacity that can be evaluated by taking the numerical derivative of $`E_k=E^{}E_I`$ with respect to $`T`$. Eq.(2) shows that a negative heat capacity corresponds to partial energy fluctuations in the microcanonical ensemble that exceed the corresponding fluctuations in the canonical ensemble ($`\sigma _{can}=C_kT^2`$).

The average freeze-out volume can be determined  from the detected fragment kinetic energies through a many-body Coulomb trajectory calculation (Fig. 4). For the sets of data here presented, a volume close to three times the normal volume reproduces the measured kinetic energies. Only for the central $`Au+Au`$ collisions, the interplay between the Coulomb energy and the radial flow (about 1 A MeV) gives an uncertainty in the determination of the freeze-out volume. Indeed, the fragment kinetic energies would also be compatible with zero radial flow and an increased Coulomb repulsion from a more compact configuration. This ambiguity, however, does not affect our main conclusions, as is evident from the results shown in Fig.5.

An estimator of the microcanonical temperature $`T`$ can be obtained by inverting the kinetic equation of state:

$$E_k=\underset{i=1}{\overset{M}{}}\frac{A_i}{a(T)}T^2+\frac{3}{2}(M1)T$$
(3)
The unknown parameters of Eq.(3) are the average side feeding correction $`\mathrm{\Delta }A`$ on the fragment masses $`A_i`$, and the level density parameter $`a`$ of primary fragments. Since the Coulomb energy is positively correlated with the charged products multiplicity, the value obtained for the freeze-out volume in each excitation energy bin depends on the side feeding correction. This means that $`V_{FO},a_i`$ and the percentage of evaporated particles have to be fixed consistently with an iterative procedure. When this is done , we find that the microcanonical temperature, evaluated through Eq.(3) results in agreement with the measured isotope temperature (left panel of Fig. 3).

In order to minimize spurious fluctuations due to unmeasured quantities, a constant side feeding correction is applied in each excitation energy bin, and no volume fluctuations are allowed . The partial energy fluctuation overcomes the canonical expectation approximately at the same values of $`E^{}`$ for peripheral and central collisions data. The same thing is true for central $`Xe+Sn`$ collisions from 32 to 50 A MeV, measured with the Indra device . Therefore, the different data sets are fully compatible, the results do not depend neither on the detector, nor on the data selection.

It is however important to stress that the quantitative study of nuclear thermodynamics is still at the first stage. Only systematic studies of correlated observables, performed with

sophisticated experimental devices, with low energy threshold, high granularity and isotopic resolution will bring a detailed understanding of the phase transition and will allow to locate the position of the observed multifragmentation in the phase diagram of nuclear matter . Recent investigations on the opening of the multifragmentation channel, on the caloric curve and on the limiting temperature for $`A100`$ systems started at moderate incident energies ($``$ 15 A MeV) at INFN Laboratori Nazionali di Legnaro with the Garfield apparatus. First results  indicate that already at an excitation energy about 3 A MeV, multifragmentation events with nearly equal size fragments represent a non negligible fraction of the decay channels. The very good mass resolution of the Garfield Si-CsI telescopes (Fig.6) allows an isotope analysis. The dependence on the isospin of the entrance channel is apparent. Isoscaling analyses, now in progress, will bring information on the isospin distillation at the onset of the multifragmentation .

The authors would like to thank Garfield collaboration for kindly providing preliminary results.