text:
# Influence of new reaction rates on 18F production in novae

## 1 Introduction

Classical novae emit gamma-ray radiation at and below 511 keV during the early epochs after the explosion. This emission is produced by electron-positron annihilation in the expanding envelope, and the subsequent Comptonization of the resulting gamma-ray photons, and it shows a line, at 511 keV, and a continuum, between 20 and 511 keV (Gómez-Gomar et al. (1998)). The positrons responsible for this emission come mainly from the disintegration of 18F (Leising & Clayton (1987); Gómez-Gomar et al. (1998)), because its lifetime ($`\tau `$=158 min) is such that positrons are emitted at the “right time”, i.e., when the expanding envelope starts to be transparent to gamma-ray radiation. Therefore, the amount of radiation emitted strongly depends on the 18F content of the nova envelope.

The synthesis of 18F in novae depends largely on some key nuclear reaction rates of 18F destruction and production which are far from being well known. This is the case, in particular, of the 18F(p,$`\gamma )^{19}`$Ne and 18F(p,$`\alpha )^{15}`$O reactions. Recent experimental studies (Graulich et al. (1997); Utku et al. (1998)) drastically improved the knowledge of these reaction rates with respect to previous studies (Wiescher & Kettner (1982)). In a recent paper (Hernanz et al. (1999)), we have analyzed the influence of these rates (Utku et al. (1998)) of 18F(p,$`\gamma )^{19}`$Ne and 18F(p,$`\alpha )^{15}`$O on the final yields of 18F for different models of CO and ONe novae. The effect of the new rates was important, since a factor of 10 reduction in the yields and in the resulting gamma-ray fluxes was obtained for all the models. Therefore, we concluded that a more detailed analysis of the reaction rates was necessary, in order to predict the gamma-ray emission of classical novae.

The rates proposed by Utku et al. (1998) were limited to relatively high resonance energies and temperature domains, more appropriate for temperatures typical of X–ray bursts than for those of classical novae. The reason is that for their determination neither the influence of uncertainties on low energy resonance strengths (based on assumed reduced width) nor the effect of all low energy resonance tails were considered. At higher temperatures the rates are more reliable since their main contributions come from two directly measured resonances. The purpose of this paper is to provide rates for the 18F(p,$`\gamma )^{19}`$Ne and 18F(p,$`\alpha )^{15}`$O reactions valid in the domain of temperature of nova nucleosynthesis, incorporating the latest experimental data. In addition to the nominal rate, we provide upper and lower limits. Since other nuclear reactions also affect 18F synthesis in novae, a global analysis is done, including results from the recent NACRE compilation (Angulo et al. (1999)). With the new rates affecting 18F synthesis, new nova models have been computed, in order to determine the mass of 18F they eject and the impact of the new yields on their early gamma-ray emission.

The organisation of this paper is the following. In section 2 we describe the nucleosynthesis of 18F. In section 3 we discuss in detail the 18F+p rates and the corresponding uncertainties while in the following section (4), we briefly discuss other recently published rates. In section 5 we present new results that show the influence of these new rates on 18F production. Following the conclusion, an appendix gives some analytical approximations to the rates.

## 2  Synthesis of 18F in classical novae

In this Section, we will focus on the nuclear paths leading to 18F synthesis, through a detailed analysis of a representative model of an ONe nova : a 1.25 M white dwarf, accreting solar-like matter at a rate of $`2\times 10^{10}`$ M.yr-1 , assuming a 50% degree of mixing with the outermost core of composition taken from Ritossa, García–Berro & Iben (1996) (mainly 16O and 20Ne with traces of other Ne, Na and Mg isotopes). The evolutionary sequences leading to a nova outburst have been computed with an updated version of the code SHIVA (see José & Hernanz 1998), a one-dimensional, implicit, hydrodynamical code in Lagrangian formulation, which follows the course of the explosion from the onset of accretion up to the expansion and ejection stages. The nuclear reaction network is described in José et al. (1999), and in Hernanz et al. (1999). In particular, we use the reaction rates for proton captures on 18F and 17O based on Utku et al. (1998) and Landré et al. (1989), respectively.

The synthesis of 18F in nova outbursts takes place within the hot-CNO cycle (see fig. 1). Regardless of the nova type (either CO or ONe, according to the composition of the underlying white dwarf), the initial 16O abundance is quite large (up to 25% by mass). Hence, initial 16O is the main source of the formation of 18F, either through the chain 16O(p,$`\gamma `$)17F(p,$`\gamma `$)18Ne($`\beta ^+`$)18F or through 16O(p,$`\gamma `$)17F($`\beta ^+`$)17O(p,$`\gamma `$)18F, which reflects a competition between (p,$`\gamma `$) reactions and $`\beta ^+`$-decays.

Snapshots of the evolution of isotopes relevant to 18F synthesis (i.e., 16,17O, 17,18,19F and 18,19Ne) are shown in Fig. 2. At the early stages of the explosion, when the temperature achieved at the burning shell reaches $`T_{bs}5\times 10^7`$ K (Fig. 2, first panel), the main nuclear activity in this mass region is driven by 16O(p,$`\gamma `$)17F, followed by its $`\beta ^+`$-decay to 17O. However, the temperature is not high enough to burn a significant fraction of 16O. On the contrary, the abundance of 17O rises by several orders of magnitude with respect to its initial content. Both 17,18F are being synthesized at the burning shell by means of proton captures on 16,17O, respectively, whereas the amount of 19F remains essentially unchanged. The minor activity driven by proton capture reactions onto 17,18F results also in a moderate increase of both 18,19Ne (below $`10^{10}`$ by mass).

A similar trend is found when $`T_{bs}7\times 10^7`$ K (Fig. 2, second panel). The temperature rise increases the number of proton captures onto 16O, leading to 17F, which rapidly decays into 17O. The effect of convection, which already extends through most of the envelope (i.e., $`150`$ km), is shown in the smooth distribution of both 17,18F, previously synthesized in the burning shell and efficiently carried away to the outer envelope. 19F is in turn reduced down to 50% due to 19F(p,$`\alpha `$)16O, which dominates 19Ne($`\beta ^+`$)19F.

When $`T_{bs}`$ reaches $`10^8`$ K (Fig. 2, third panel), 17O shows a flat profile along the envelope, with a mean mass fraction of $`3\times 10^3`$, due to 16O(p,$`\gamma `$)17F($`\beta ^+`$)17O, which dominates destruction through 17O(p,$`\gamma `$)18F. With respect to fluorine isotopes, we stress an important synthesis of 17,18F at this stage, driven by 16O(p,$`\gamma `$)17F and 17O(p,$`\gamma `$)18F, whereas 19F is almost fully destroyed. 18,19Ne continue to rise due to proton captures onto 17,18F.

A dramatic change is found as soon as the burning shell reaches $`2\times 10^8`$ K (Fig. 2, fourth panel). At this time, the mass fraction of 16O has been reduced down to $`3\times 10^2`$. For the first time, 17O decreases (by one order of magnitude with respect to the amount shown in previous panel), since destruction through 17O(p,$`\gamma `$)18F dominates 17F($`\beta ^+`$)17O in the vicinity of the burning shell. On the contrary, 17F exhibits a significant rise, up to $`3\times 10^2`$ by mass (i.e., one order of magnitude higher than in previous panel), so that the amount of 17F becomes larger than that of 17O in almost all the envelope. Again, 16O(p,$`\gamma `$)17F dominates both 17O(p,$`\gamma `$)18Ne and 17F($`\beta ^+`$)17O near the burning shell. The evolution of 19F is still dominated by 19F(p,$`\alpha `$)16O, which in fact has no significant influence on the 16O content. Also noticeable is the dramatic rise of both 18,19Ne, which increase by several orders of magnitude, mainly due to the fact that at such temperatures, proton captures onto 17,18F become faster than the corresponding $`\beta ^+`$-decays. It is worth noticing that at this stage, a non-negligible path leading to 18F synthesis is now driven by 18Ne($`\beta ^+`$)18F, basically at the outer envelope shells, which compensates an efficient destruction through proton captures near the burning shell.

The burning shell achieves a maximum temperature of $`2.51\times 10^8`$ K (Fig. 2, fifth panel). 16O is reduced to 0.096 by mass in the burning shell. Whereas 17O is destroyed by (p,$`\alpha `$) and (p,$`\gamma `$) reactions at the burning shell, its mean abundance in the envelope increases due to 17F($`\beta ^+`$)17O which dominates in almost all the envelope. On the contrary, 17F is generated in the burning shell by means of 16O(p,$`\gamma `$)17F, whereas it is being destroyed by $`\beta ^+`$-decays in the outer envelope. The evolution of 18F reflects a competition between destruction and creation modes at different locations of the envelope: whereas its content is reduced by both 18F(p,$`\alpha `$)15O reactions at the burning shell and by 18F($`\beta ^+`$)18O at the outer envelope, a dominant source for 18F synthesis through 17O(p,$`\gamma `$)18F is found at the intermediate shells. Proton captures onto 17,18F and convective transport continue to pump 18,19Ne to the outer envelope shells, which is at the origin of the rise of 19F at this stage (through 19Ne($`\beta ^+`$)19F).

Shortly after, due to the sudden release of energy from the short-lived species 13N, 14,15O and 17F, the envelope begins to expand. As a result of the drop in temperature, (p,$`\gamma `$) and (p,$`\alpha `$) reactions are basically restricted to the vicinity of the burning shell, whereas most of the envelope is dominated by $`\beta ^+`$-decays (Fig. 2, sixth and seventh panels). Hence, whereas 17O is powered by 17F($`\beta ^+`$)17O along the envelope, the amount of 17F and 18,19Ne decreases significantly. The 18F abundance at these late-time stages of the outburst increases due to 17O(p,$`\gamma `$)18F, which is dominant at the intermediate layers of the envelope (whereas in the vicinity of the burning shell 18F is efficiently destroyed by 18F(p,$`\alpha `$)15O, the outer envelope is dominated by 18F($`\beta ^+`$)18O). Moreover, destruction of 19F through 19F(p,$`\alpha `$)16O at the burning shell is dominated by 19Ne($`\beta ^+`$)19F elsewhere.

The final stages of the outburst (Fig. 2, eighth panel), as the envelope expands and cools down, are dominated by the release of nuclear energy by $`\beta ^+`$-decays such as 18F($`\beta ^+`$)18O, or 17F($`\beta ^+`$)17O. The resulting mean abundance of 18F in the ejected shells in this Model is X(18F)= $`2.7\times 10^4`$. Most of the envelope is, however, dominated by the large abundances of 16,17O (X(16O)= $`6.8\times 10^2`$, X(17O)= $`3.9\times 10^2`$). A residual 17F, which is still decaying into 17O, and a non-negligible amount of 19F ($`3\times 10^6`$, by mass) are also present.

In summary, the synthesis of 18F in classical novae is essentially controlled by five proton-capture reactions, 16O(p,$`\gamma `$)17F, 18F(p,$`\gamma `$)19Ne, 18F(p,$`\alpha `$)15O, 17O(p,$`\gamma `$)18F, 17O(p,$`\alpha `$)14N, and several $`\beta ^+`$-decays. Accordingly, the corresponding reaction rates deserve further attention and their uncertainties are discussed in Sect. 3 for the 18F+p rates and Sect. 4 for the other rates.

## 3 The 18F(p,$`\alpha )^{15}`$O and 18F(p,$`\gamma )^{19}`$Ne reaction rates

Uncertainties on the 18F(p,$`\gamma )^{19}`$Ne and 18F(p,$`\alpha )^{15}`$O reaction rates used to reach many orders of magnitude due to the very limited spectroscopic data available for 19Ne in the domain of interest. The rate estimated by Wiescher & Kettner (1982) has become obsolete by the recent measurements of Rehm et al. (1995), Coszach et al. (1995), Rehm et al. (1996, 1997) and in particular by the work of Graulich et al. (1997), Utku et al. (1998) and Butt et al. (1998). We concentrate on the faster 18F(p,$`\alpha )^{15}`$O reaction but most of the discussion also applies to the 18F(p,$`\gamma )^{19}`$Ne one.

Before 1997 and below the excitation energy ($`E_x`$) corresponding to the resonance energy ($`E_r`$) of 0.4 MeV, only two levels were known $`E_x`$ = 6.437 and 6.742 MeV but with unknown or uncertain spins and parity ($`J^\pi `$)111For more detailed explanations of the standard nuclear physics notations used, see e.g. the appendix in José et al. (1999).. On the contrary, various levels were known in the corresponding region of 19F, the conjugate nucleus (see Fig. 3). For instance, in the $`E_r`$=0.–1. MeV region, only 6 levels were known and many spins unknown in 19Ne, while 19 levels had been observed in the corresponding 19F region. Hence, for the estimation of the rates, Wiescher & Kettner (1982) considered three known levels at $`E_x`$ = 6.437, 6.742 and 6.862 MeV and postulated three others. The corresponding rate was highly uncertain since the location of these three postulated levels was approximate and because other levels found in 19F are expected to have counterparts in 19Ne in the region of interest. This is the case, in particular, of the $`E_x,J^\pi `$ = 6.429, 1/2-, 6.497, 3/2+ and 6.527, 3/2+ levels of 19F corresponding to strong expected resonances because of their low centrifugal barrier, i.e. low transfered orbital angular momentum ($`\mathrm{}_p`$) of 1, 0 and 0 respectively.

Even though the recent experiments have not been able to find all the counterparts of the 19F levels (see Fig. 3), they provided direct measurement for the strengths of two resonances which likely dominate the rate in the domain of nova nucleosynthesis and the location of several new 19Ne levels. The strength of the resonance corresponding to the $`E_x`$ = 6.742 MeV level ($`E_r`$=330 keV) has recently been measured directly by Graulich et al. (1997) with a 18F beam provided by the Louvain–La–Neuve facility. This level is thought to be the analog of the 6.787, 3/2- 19F level. The measured strength ($`\omega \gamma `$), $`3.5\pm 1.6`$ eV (Graulich et al. (1997)), is in good agreement with the Wiescher & Kettner (1982) estimate and accordingly does not induce a significant change in the rate. On the contrary, the resonance associated with the $`E_x,J^\pi `$ = 7.067 MeV, 3/2+ level ($`E_r`$=659 keV) strongly modifies the rates. It has been studied by Rehm et al. (1995, 1996), Coszach et al. (1995), Graulich et al. (1997), Utku et al. (1998) and Butt et al. (1998) and its strength has been measured directly (Coszach et al. (1995); Graulich et al. (1997)). It is located well outside of the Gamow peak but due to its total width ($`\mathrm{\Gamma }`$30 keV), the contribution of its tail alone is greater than the Wiescher & Kettner (1982) rate in the domain of nova nucleosynthesis. However, no counterpart of this broad level was known in 19F. The closest 3/2+ 19F level is located at $`E_x`$ = 7.262 but its width is much smaller ($`\mathrm{\Gamma }<`$6 keV). To remove this incompatibility, Rehm et al. (1996) claimed that the 7.067 MeV width could be smaller ($`\mathrm{\Gamma }`$14 keV). A 7.114, 7/2+ level was known in 19F with the possibility that it hides an unresolved $`J^\pi `$ = 3/2+ level as suggested by Smotrich et al. (1961). A recent gamma ray spectroscopy study (Butt et al. (1998)) of 19F produced by the 15N($`\alpha ,\gamma )^{19}`$F reaction has confirmed the presence of a 7.101 MeV, 3/2+, $`\mathrm{\Gamma }`$ = $`28\pm 1`$ keV level. Hence the high width found by Coszach et al. (1995) and Utku et al. (1998) for the $`E_x,J^\pi `$ = 7.067 MeV, 3/2+, 19Ne level can be understood. (Note however that this analog assignment has been questioned very recently by Fortune & Sherr (2000).)

Utku et al. (1998) have also found three new 19Ne levels and have made tentative assignments of analog levels in 19F (see Fig.3). They recognized the 6.437 MeV, 1/2- level ($`E_r`$=26 keV) as the analog of the 6.429, 1/2-, $`\mathrm{\Gamma }`$ = 280 keV one in 19F from its large measured width ($`216\pm 19`$ keV). Even though it is located at a very low energy (in the context of novae) it is so broad that its tail may lead to a significant contribution to the rate (see Fig.4). The 6.449 MeV 3/2+ ($`E_r`$=38 keV) level is one of the three new ones found by Utku et al. (1998). It is also of great importance since the contribution of its tail can dominate the astrophysical factor in the relevant energy range (see Fig.4). In comparison the two other new levels give a smaller (6.698 MeV 5/2+, $`E_r`$=287 keV) or negligible contribution (6.419 MeV 3/2+, $`E_r`$=8 keV).

A new rate has been provided by Utku et al. (1998), using their new data together with available spectroscopic information on 19Ne and 19F supplemented by estimates of radiative and proton widths when missing. This rate is already higher by a factor of up to $`3\times 10^2`$ when compared with Wiescher & Kettner (1982). Nevertheless, even though this rate is now set on a firmer basis than the Wiescher & Kettner (1982) one, it remains quite uncertain as it depends directly on the estimated proton widths. Except for the two resonances whose strengths have been measured directly by Graulich et al. (1997), the proton widths are assumed to be equal to 0.1 or 0.01 Wigner limit, i.e. $`\theta _p^2`$=0.1 or 0.01, according to their parity (Utku et al. (1998)). For the resonances considered here, one has $`\mathrm{\Gamma }_p\mathrm{\Gamma }_{Tot}\mathrm{\Gamma }_\alpha `$ so that their strengths are ($`\omega \gamma \omega \mathrm{\Gamma }_p`$) proportional to chosen $`\theta _p^2`$ values. Accordingly, to obtain lower and upper limits for the rate, we used the two extreme cases $`\theta _p^2`$=0 and $`\theta _p^2`$=1 respectively and regardless of parity but kept Utku et al. (1998) $`\theta _p^2`$ values for the nominal rate (see Fig.6). The lower limit is quite certain as it is given by the contribution of the two directly measured resonances whose parameters are now reliably determined. It remains higher than the Wiescher & Kettner (1982) rate above $``$50$`\times 10^6`$ K. The nominal rate is conditioned by the rather arbitrary choice of $`\theta _p^2`$=0.1 or 0.01 that we adopt following Utku et al. (1998). However, our nominal rate is significantly higher \[up to a factor of $``$10 and 30 around 108 K for (p,$`\alpha `$) and (p,$`\gamma `$) reactions\] than the one from Utku et al. (1998) for typical nova temperatures (see Fig.6). This is due to the inclusion of the contribution of the tail ($`\mathrm{\Gamma }_{Tot}=\mathrm{\Gamma }_\alpha `$=4.3 keV) of the 6.449 MeV level (see Fig. 4) not considered by Utku et al. (1998). The corresponding rates (p,$`\alpha `$) are presented in Fig.6 and have been calculated by numerical integration of the Breit-Wigner formula for all broad resonances. For the (p,$`\gamma `$) rates, we used the same procedure to obtain the low, nominal and high rates (see Figs. 5 and 7). Radiative widths were taken from 19F analog levels as in Utku et al. (1998) except for the important resonance at 659 keV as Bardayan et al. (1999) provided an experimental value ($`\mathrm{\Gamma }_\gamma `$=0.39 eV) for the analog level. We also adopted the direct capture (DC) contribution from Utku et al. (1998) even though the spectroscopic factors used (18O$``$p instead of 18F$``$p or even 18F$``$n) may not be the more appropriates. These rates have been calculated using all the available spectroscopic data collected in Table II of Utku et al. (1998) except for the inclusion of a new gamma width for the 659 keV resonance and for a different proton width for the 26 keV resonance \[1% of the Wigner limits i.e. 2.5$`\times 10^{17}`$ eV, as for other negative parity resonances, instead of 6.6$`\times 10^{20}`$ eV of unknown origin in Utku et al. (1998)\]. These numbers are taken at face value to calculate the nominal rate only. For the low and high rates unknown proton widths are allowed to vary between 0 and the Wigner limit as discussed above. It is clear that this range ($`0<\theta _p^2<1`$) represents the most extreme values as are the corresponding rates. This should be kept in mind when interpreting the astrophysical implications.

However a few points have been neglected due to the lack of experimental or theoretical information. First of all, the identification of all $`E_x`$ $`<`$ 6.6 MeV analog levels is not complete yet (see Fig. 3) and is not certain. In particular, new contributions to the nominal and high rates corresponding to missing levels cannot be ruled out. The tail of the 38 keV resonance has been calculated using the large total width derived from the 19F analog level proposed by Utku et al. (1998). Another 3/2+ level, with smaller width, lies 30 keV below in 19F. It has been assigned (Utku et al. (1998)) to be the analog level corresponding to the 8 keV (Ex=6.419 MeV) resonance in 18F+p (see Fig. 3). It is possible that the analog assignment for these two levels is inversed or that they are mixed but at least one or the other resonance (8 or 38 keV) should be broad enough to dominate the rates at low nova temperature. For the (p,$`\gamma `$) rates, we neglected any difference of radiative widths between analog levels and the uncertainty on the DC contribution in front of the uncertainty on the proton widths. These rates do not either include possible interference effects arising from the interference between the 3/2+ resonances located at 38 and 659 keV. In a favourable case (when $`\theta _p^2`$0.01 for the 38 keV resonance), destructive interference could reduce the low rate by a factor of $``$10 in the sensitive $`T10^8`$ K region. However, at present the rates provided in the appendix are those we recommend.

## 4 Other reactions affecting 18F production

According to the recent compilation of charged-particle induced thermonuclear reaction rates NACRE (Angulo et al. 1999), the rate for 16O(p,$`\gamma `$)17F is rather well known and suffers from a little uncertainty (i.e., a factor of $``$ 2 between the high and low rates). Even though 17F has a much shorter lifetime than 18F, its destruction by proton capture is rather well known at nova temperatures and should not affect the analysis of Sect. 2. Below $``$4$`\times 10^8`$ K, the 17F(p,$`\gamma )^{18}`$Ne reaction rate is dominated by the effect of direct capture on bound 18Ne levels. This contribution is expected to suffer little uncertainty since it has been calculated (Wiescher & Kettner (1982); García et al. (1991)) using experimental neutron spectroscopic factors from the mirror nucleus 18O. The lowest resonances are located at $`E_r`$ = 0.586, 0.655 and 0.600 MeV (Bardayan et al. (1999)) too high to contribute, even by their tails, to the rate below $``$4$`\times 10^8`$ K. Partial widths are reasonably known through experiments, shell model calculations or from the conjugate levels (Wiescher & Kettner (1982); García et al. (1991)) and are such that $`\mathrm{\Gamma }_p\mathrm{\Gamma }_\gamma `$ (i.e. $`\omega \gamma =\omega \mathrm{\Gamma }_\gamma `$). A long standing uncertainty (Wiescher & Kettner (1982); Wiescher, Görres & Thielemann (1988); García et al. (1991)) concerned the location of the $`\mathrm{}`$=0, 0.600 MeV resonance associated with the 4.524 MeV, 3+ level. Fortunately, the resonance energy first measured by García et al. (1991), has been experimentally confirmed very recently by Bardayan et al. (1999). Hence the uncertainty associated with the 17F(p,$`\gamma )^{18}`$Ne rate appears negligible within the domain of nova nucleosynthesis.

According to NACRE (Angulo et al. (1999)), the 17O(p,$`\alpha )^{14}`$N (Fig. 8) and 17O(p,$`\gamma )^{18}`$F (Fig. 9) reaction rates present large uncertainties at temperatures below a few 108 K. The increase in the rates around 50$`\times 10^6`$ (Angulo et al. (1999)) is due to the contribution of a $`E_r`$ = 66 keV resonance (Landré et al. (1989); Blackmon et al. (1995)), not included in Caughlan & Fowler (1988). The large uncertainty around 2$`\times 10^8`$ K (Angulo et al. (1999)), well within the range of temperature reached in nova outbursts, is due to the unknown contribution of a 179.5 keV resonance associated with the Ex = 5786 keV (18F) level. The upper limits for the strengths \[(p,$`\gamma `$) and (p,$`\alpha `$)\] of this resonance are calculated (Angulo et al. (1999)) from the partial widths extracted by Rolfs, Berka & Azuma (1973) and the upper limit for the proton width from Landré et al. (1989). As usual, for the calculation of the contribution of this resonance to the NACRE low, recommended and high rates, its strength has been taken equal to 0, 10% and 100% of the experimental upper limit. Accordingly, the NACRE recommended rate is somewhat arbitrary in the temperature domain affected by the 179.5 keV resonance (Figs. 8 and Figs. 9)

## 5 Results: effect of new rates on 18F synthesis

New models of a 1.25 M ONe nova have been computed with the SHIVA code, in order to analyze the effect of the new nuclear reaction rates. For the 18F+p rates, we have adopted the low, high and nominal prescriptions described in section 3, whereas for the 17O+p ones we have adopted the NACRE rates (see Table LABEL:t:yields). It is important to stress that the energetics of the explosion remains practically unchanged with the new 18F+p and 17O+p reaction rates, since these reactions are not the ones affecting most the evolution (see José & Hernanz (1998) and José et al. (1999) for a deep analysis of the reaction rates which have the largest influence in nova explosions). Therefore, the maximum temperature attained in the burning shell at the base of the accreted envelope, Tmax=2.51$`\times 10^8`$ K, the mean kinetic energy and the ejected mass, $`\mathrm{\Delta }`$$`M_{ejec}`$=1.8$`\times 10^5`$M , are unchanged with respect to previous models (i.e., José et al. (1999)).

Two previously computed evolutionary sequences, with the same white dwarf mass and composition, are adopted as old models for the purpose of comparison: a case with Wiescher & Kettner (1982) rates (see José & Hernanz (1998) for the yields and Gómez-Gomar et al. (1998) for the $`\gamma `$-ray spectra) and a case with Utku et al. (1998) rates (see Hernanz et al. (1999) for the yields and the spectra). In both of them, the Caughlan & Fowler (1988) updated by Landré et al. (1989) rates for 17O+p were adopted. In Table LABEL:t:yields we show the 18F yields obtained with the new rates, as well as those with the old models. We see that yields with the nominal 18F+p rates are smaller by a factor of 30 than those with the older rates, as expected. In fact, the reduction in 18F production is even larger (by a factor of 3) than the one obtained with Utku et al. (1998) rates (see Hernanz et al. (1999)), as a consequence of the inclusion of the tail of the resonance at 38 keV (see section 3).

Concerning the effect of the new 17O+p rates, a reduction by an extra factor of 2 is obtained when the NACRE recommended rates, instead of the old Caughlan & Fowler (1988) and Landré et al. (1989) ones, are adopted (together with 18F+p nominal rates; see Table LABEL:t:yields). In summary, our nominal yield of 18F is 4.84$`\times 10^5`$, which is 60 times smaller than the yield with the old Wiescher & Kettner (1982) rates (José & Hernanz (1998); Gómez-Gomar et al. (1998)) and 6 times smaller than that with Utku et al. (1998) rates (Hernanz et al. (1999)). The consequences of the reduced 18F yields are quite large for the early gamma-ray output, which is directly related to the amount of 18F synthesized. The reduced content of 18F in an expanding envelope with similar physical and dynamical properties, translates directly into a reduced positron annihilation gamma-ray flux, roughly by the same factor as the 18F decrease. This was shown in Hernanz et al. (1999), where complete computations of the gamma-ray spectra were performed with the 18F yields obtained with Utku et al. (1998) rates, for different masses and compositions of the white dwarf. With the rates presented here, the reduction of the predicted fluxes becomes larger by a factor of 6, for the combined nominal 18F+p and recommended 17O+p rates. In summary, if we compare with results in Gómez-Gomar et al. (1998), where 18F yields were computed with Wiescher & Kettner (1982), the 511 keV and continuum fluxes for a 1.25 M ONe nova should be reduced by a factor of 60.

Another important point to stress is the impact of the uncertainty in the rates on the computed yields. The above nominal yields rely upon recommended or nominal rates which contain assumed values ($`\theta _p^2`$=0.01 or 0.1 for 18F+p and 10% of an experimental upper limit for 17O+p). Hence, the nominal yield should not be dissociated from the large remaining uncertainties. From our results (see Table LABEL:t:yields) we conclude that the range between high and low presciptions for the 18F+p rates translates into an uncertainty (maximum versus minimum yield ratio) of a factor of 310, a quite large value, whereas for the 17O+p rates the uncertainty is a factor of 10. (It is worth noticing that the results obtained with the Utku et al. (1998) rates fall within the limits obtained with the new rates; see Table LABEL:t:yields). The corresponding uncertainties in the gamma-ray fluxes are of the same order of magnitude. This points out the interest of more accurate determinations of the 18F+p and 17O+p rates.

## 6 Conclusions

We have investigated the 18F formation and destruction in nova outbursts, identified the key reactions (proton capture on 18F and 17O) and analysed their rates. The proton capture rates on 18F are higher than the Wiescher & Kettner (1982) ones at nova temperatures due to the contribution of the tail of the 659 keV resonance whose large measured width (Coszach et al. (1995); Graulich et al. (1997)) has been indirectly confirmed (Butt et al. (1998)). Another important contribution comes from the tail of the 38 keV resonance which was neglected in previous studies. Its nominal contribution is larger than the 659 keV one but is proportional to its assumed reduced proton width $`\theta _p^2`$. The strengths (proton widths) of the low lying resonances are unknown and induce large uncertainties (factors of 100 to 1000) in the rates at nova temperatures. We have provided updated nominal rates for the two capture reactions together with upper and lower limits. The 17O+p rates also display some large uncertainties at nova temperatures according to the recent compilation of Angulo et al. (1999).

We have used these new nuclear physics results in a fully hydrodynamical nova code to calculate the 18F yields in novae for different rates : our low, nominal and high 18F+p rates, and the low, recommended and high 17O+p rates from Angulo et al. (1999). These results have been compared with models computed with the old Wiescher & Kettner (1982) rates (José & Hernanz (1998)), and with more recent models with the Utku et al. (1998) rates (Hernanz et al. (1999)). Two important results have been obtained. First, there is always a reduction of the amount of 18F synthesized in a nova explosion, with the nominal rates for the 18F+p reactions both alone and combined with the recommended rates for the 17O+p reactions. The nominal 18F yield (nominal 18F+p and recommended for 17O+p) is 4.84$`\times 10^5`$ by mass, which is 60 times smaller than the one obtained with Wiescher & Kettner (1982) rates and 6 times smaller than the one with Utku et al. (1998) rates (and Caughlan & Fowler (1988) and Landré et al. (1989) for 17O+p). The impact on the early gamma-ray spectrum of the nova is a reduction of the flux by the same amount (with respect to Gómez-Gomar et al. (1998) and Hernanz et al. (1999), respectively). Second, the yields are found to be very sensitive to the rates with resulting combined (18F+p and 17O+p) uncertainties of more than three orders of magnitude. This supports the need of new experimental and theoretical studies to improve the knowledge of the 18F+p and 17O+p rates and, consequently, allow for a larger reliability of the predictions of annihilation gamma–ray fluxes from novae, to be observed by current and future instruments.

###### Acknowledgements.

This work was partially supported by PICS 319, PB98-1183-C02, PB98-1183-C03 and ESP98-1348.

Appendix: Updated rates

We give here the tabulated reaction rates for proton capture on 18F resulting from numerical integrations (see § 3 and Figs. 6 and 7). For convenience, we also provide here formulas that approximate the nominal rates (As usual, T9 is the temperature in GK, T9XY=T9\**(X/Y) and T9LN=LOG(T9) in FORTRAN notations.)

The following 18F(p,$`\alpha )^{15}`$O nominal rate formula includes the contribution of the two measured resonances (330 and 659 keV) with additional contributions from the possible resonances at 38 keV, (with high energy tail) and 287 keV.

```
      SVDir = 9.13e10/T923*EXP(-18.052/T913-0.672*T92)*(1.+! Tail
     & 0.0231*T913+6.12*T923+0.988*T9+9.92*T943+4.07*T953) ! 659 keV
     & +5.78e5/T932*EXP(-3.830/T9)+9.91e8/T932*EXP(-7.648/T9)+ ! 330, 659 keV
     & 2.81e-6/T932*EXP(-0.441/T9)+2706.*EXP(3.8319*T9LN-1.3450* ! 38 keV
     & T9LN**2-.0001/T9**3)+4.46e4/T932*EXP(-3.331/T9) ! 287 keV

```

The following 18F(p,$`\gamma )^{19}`$Ne nominal rate formula contains contributions from direct capture and from the 38, 287, 330 and 659 keV resonances.

```
      SVDir = 3.98e7/T923*EXP(-18.052/T913)* ! DC
     & (1.+0.0231*T913+0.0885*T923+0.0143*T9)+
     & 1.34e3/T932*EXP(-3.830/T9)+1.51e4/T932*EXP(-7.648/T9)
     & +EXP(-0.56781+3.6850*T9LN-1.3636*T9LN**2-.0001/T9**3)
     & +8.26e-10/T932*EXP(-0.441/T9)+12.2/T932*EXP(-3.331/T9)

```


format: text

text:
# Phase changes in 38-atom Lennard-Jones clusters. II: A parallel tempering study of equilibrium and dynamic properties in the molecular dynamics and microcanonical ensembles

## I Introduction

The simulation of systems having complex potential energy surfaces (PES) is often difficult owing to the problem of quasi-ergodicity. Quasi-ergodicity can arise in systems having several energy minima separated by high energy barriers. When such situations occur, as for example in proteins, glasses or clusters, the system can become trapped in local basins of the energy landscape, and the ergodic hypothesis fails on the time scale of the simulation. In the canonical ensemble, the high energy regions of the PES are exponentially suppressed and barrier crossings become rare events. Calculations of equilibrium properties when phase space is thus partitioned require methods that overcome quasi-ergodicity by enhanced barrier crossing. Many techniques have been proposed to address this problem, including the use of generalized ensembles such as multicanonical or Tsallisian, simulated tempering, configurational or force bias Monte Carlo, or various versions of the jump-walking algorithm. Most of these techniques have been introduced for Monte Carlo (MC) simulations rather than molecular dynamics (MD) simulations. These techniques have been applied to a variety of sampling and optimization problems, and phase changes in clusters have often been considered as a benchmark to test these methods.

The double-funnel energy landscape of the 38-atom Lennard-Jones (LJ) cluster has been investigated in detail by Doye, Miller and Wales, who recently also estimated the inter-funnel rate constants using master equation dynamics. This landscape is challenging for simulation because of the high free-energy barrier separating the two funnels. In the preceding paper (hereafter referred to as I), we have shown how the parallel tempering algorithm can be used to deal with this particularly complex system for Monte Carlo simulations in the canonical ensemble. Achieving ergodicity in microcanonical simulations is much harder than in the canonical ensemble, because the system is unable to cross any energy barrier higher than the total energy available. The 38-atom Lennard-Jones cluster is fundamentally non-ergodic in a range of energies. This non-ergodicity may not be a serious problem when considering one particular cluster on a short time scale. However, in a statistical sample of such systems it is important to have ergodic results.

To allow MD simulations to cross the high energy barriers, one may think of heating the system (by increasing its kinetic energy), followed by a cooling back to its initial thermodynamic state. Although this process is straightforward, the dynamics becomes biased and non physical during the heating and cooling processes. Moreover, it is difficult to control accurately the heating and cooling rates that should be chosen for any system. This latter aspect is particularly critical for the 38-atom Lennard-Jones cluster where the narrow and deepest funnel is hard to reach even at high temperatures.

Because of the inherent difficulties of molecular dynamics, MC approaches can be especially useful for dealing with the problem of crossing high energy barriers. Monte Carlo methods have been developed in previous work in the microcanonical ensemble. In these approaches the microcanonical sampling is at fixed energy without any additional constraints. Such methods can be contrasted with isoenergetic molecular dynamics where the total, linear and angular momenta are also conserved. These additional constraints must be considered even at zero angular momentum. To differentiate microcanonical simulations, where only the energy is fixed, from molecular dynamics simulations, where additional constraints are imposed, we identify the former to be simulations in the microcanonical ensemble and identify the latter simulations to be in the molecular dynamics ensemble. The differences in the two ensembles can be particularly important when the angular momentum is large enough to induce structural (centrifugal) distortions. Because dynamical properties are calculated using molecular dynamics methods, in this work we find that a combination of Monte Carlo and molecular dynamics methods are most convenient for developing ergodic approaches to dynamics.

In this paper, we adapt the parallel tempering method to both the microcanonical and molecular dynamics ensembles. The application of parallel tempering in the molecular dynamics ensemble requires the incorporation of the conservation of the total linear and angular momenta into the probabilities. In order to extract ergodic dynamical properties, we combine Monte Carlo methods with molecular dynamics to develop a hybrid ergodic MC/MD algorithm. The efficiency of the simulation tools developed in this work is demonstrated by applications to the 38-atom Lennard-Jones cluster, which exhibits a solid-solid transition prior to melting. This transition to an equilibrium phase between truncated octahedral and icosahedral geometries makes this cluster an ideal candidate for investigating how the ergodic hypothesis can influence the dynamical behavior of a complex system.

The contents of the remainder of this paper are as follows. In the next section, we recall the basic principle of Monte Carlo sampling in the microcanonical ensemble, and we present the simple modifications needed to include parallel tempering. We test microcanonical parallel tempering methods on the 38-atom Lennard-Jones cluster, and compare the microcanonical results with those found in I using the canonical ensemble. We focus on equilibrium properties, including the caloric curve and the isomers distributions. In section III we review the characteristics of the molecular dynamics ensemble at fixed total linear and angular momenta and fixed total energy. We extend the parallel tempering Monte Carlo method to the MD ensemble, and we combine microcanonical parallel tempering with molecular dynamics to produce an ergodic MD method. We also apply these methods to several dynamical properties of LJ38; in particular the diffusion constant and the largest Lyapunov exponent. We summarize our findings and discuss our results in section IV.

## II Parallel tempering Monte Carlo in the microcanonical ensemble

The fundamental quantity in the microcanonical ensemble is the density of states $`\mathrm{\Omega }`$. For a system having $`N`$ identical particles, volume $`V`$ and total energy $`E`$, $`\mathrm{\Omega }`$ is defined by

$$\mathrm{\Omega }(N,V,E)=\frac{1}{N!h^{3N}}\delta [H(𝐫,𝐩)E]d^{3N}rd^{3N}p$$
(1)
where $`h`$ is Planck’s constant and where $`H(𝐫,𝐩)`$ denotes the classical Hamiltonian function of the coordinates $`𝐫`$ and momenta $`𝐩`$ of the $`N`$ particles. Knowing the microcanonical density of states $`\mathrm{\Omega }`$, one can calculate the canonical partition function $`Q(N,V,T)`$ by a Laplace transformation. The kinetic part of the Hamiltonian $`H`$ is quadratic in the momenta, and Eq. (1) can be partly integrated to give

$$\mathrm{\Omega }(N,V,E)=\left(\frac{2\pi m}{h^2}\right)^{3N/2}\frac{1}{N!\mathrm{\Gamma }(3N/2)}\mathrm{\Theta }[EU(𝐫)][EU(𝐫)]^{3N/21}d^{3N}r.$$
(2)
In Eq.(2), $`\mathrm{\Gamma }`$ is the Gamma function, $`m`$ is the mass of each particle, $`U(𝐫)=HK`$ is the potential energy and $`\mathrm{\Theta }`$ is the Heaviside step function: $`\mathrm{\Theta }(x)=1`$ if $`x0`$, 0 otherwise. Microcanonical averages of a coordinate-dependent variable $`A(𝐫)`$ can be expressed

$$A(N,V,E)=\frac{{\displaystyle \mathrm{\Theta }[EU(𝐫)][EU(𝐫)]^{3N/21}A(𝐫)d^{3N}r}}{{\displaystyle \mathrm{\Theta }[EU(𝐫)][EU(𝐫)]^{3N/21}d^{3N}r}}.$$
(3)
The microcanonical entropy $`S`$ can be defined by $`S(N,V,E)=k_B\mathrm{ln}\mathrm{\Omega }(N,V,E)`$ with $`k_B`$ the Boltzmann constant. The thermodynamic temperature $`T(N,V,E)`$ is given by the thermodynamic relation $`(S/E)_{N,V}=1/T`$, and can be obtained from a microcanonical average

$$T(N,V,E)=\frac{2}{3N2}\frac{1}{K^1}.$$
(4)
This expression is slightly different from the kinetic temperature $`2K/3N`$, which is a consequence of our choice in the definition of the entropy. As discussed by Pearson and co-workers, it is also possible to define the entropy using the phase space volume

$$\mathrm{\Phi }(N,V,E)=_0^E\mathrm{\Omega }(N,V,E^{})𝑑E^{}.$$
(5)
Definitions of the temperature based on $`\mathrm{\Omega }`$ differ from the temperature based on $`\mathrm{\Phi }`$ to order $`1/N`$, and the two definitions agree only in the thermodynamic limit.

Monte Carlo simulations can be used to explore the microcanonical ensemble by performing a random walk in configuration space. In the standard Metropolis scheme, a trial move from configuration $`𝐫_o`$ to configuration $`𝐫_n`$ is accepted with the probability

$$\mathrm{a}cc(𝐫_o𝐫_n)=\mathrm{min}(1,\frac{\rho _E(𝐫_n)T(𝐫_n𝐫_o)}{\rho _E(𝐫_o)T(𝐫_o𝐫_n)}),$$
(6)
where $`T(𝐫_o𝐫_n)`$ is a trial probability. The acceptance probability expressed in Eq.(6) insures detailed balance so that the random walk visits points in configuration space proportional to the equilibrium distribution $`\rho _E(𝐫)`$ defined by

$$\rho _E(𝐫)=\zeta ^1\mathrm{\Theta }[EU(𝐫)][EU(𝐫)]^{3N/21}$$
(7)
where $`\zeta `$ is the normalization. In practice, $`T(𝐫_o𝐫_n)`$ is a uniform distribution of points of width $`\mathrm{\Delta }`$ centered about $`𝐫_o`$, and $`\mathrm{\Delta }`$ is adjusted as a function of the energy so that not too many trial moves are either accepted or rejected.

Implementation of microcanonical Monte Carlo is as easy as its canonical version. Because Monte Carlo methods are based on random walks in configuration space, in principle the system can cross energy barriers higher than the available energy. However, in difficult cases like LJ38, large atomic displacements having poor acceptance ratios are needed to reach ergodicity.

Parallel tempering has proved to be an important approach to insure ergodicity in canonical Monte Carlo simulations, and parallel tempering can be easily adapted to the microcanonical ensemble by replacing the Boltzmann factors in the acceptance probability by the microcanonical weight $`\rho _E(𝐫)`$. In the parallel tempering scheme, several microcanonical MC simulations are performed simultaneously at different total energies $`\{E_i\}`$. With some predetermined probability, two simulations at energies $`E_i`$ and $`E_j`$ attempt to exchange their current configurations, respectively $`𝐫_i`$ and $`𝐫_j`$, and this exchange is accepted with probability

$`\mathrm{min}(1,{\displaystyle \frac{\rho _{E_i}(𝐫_j)\rho _{E_j}(𝐫_i)}{\rho _{E_i}(𝐫_i)\rho _{E_j}(𝐫_j)}}).`$

The acceptance ratio is analogous to the canonical expression given in I. In microcanonical simulations the potential energies must be smaller than $`\mathrm{min}(E_i,E_j)`$; otherwise the move is rejected. Parallel tempering microcanonical MC works in the same way as in standard canonical MC. As with canonical parallel tempering MC, the gaps between adjacent total energies must be chosen to be small enough so that exchanges between the corresponding trajectories are accepted with a reasonable probability.

By using a histogram reweighting technique, it is possible to extract from the MC simulations the density of states $`\mathrm{\Omega }`$, and then all the thermodynamic quantities in both the microcanonical and the canonical ensembles. The procedure is similar to that described in Ref. , and relies on the calculation of the distribution $`P(U,E)`$ of potential energy $`U`$ at the total energy $`E`$. $`P`$ is fitted to the microcanonical form $`P(U,E)=\mathrm{\Omega }_C(U)(EU)^{3N/21}/\mathrm{\Omega }(E)`$, where $`\mathrm{\Omega }_C`$ stands for the configurational density of states, and $`\mathrm{\Omega }(E)`$ is extracted by convolution of $`\mathrm{\Omega }_C(U)`$ and $`(EU)^{3N/21}`$.

We have tested the parallel tempering Monte Carlo algorithm in the microcanonical ensemble on the 38-atom Lennard-Jones cluster previously investigated. Forty different total energies ranging from $`172.4737\epsilon `$ to $`124\epsilon `$ have been used, and the same simulation conditions have been chosen as in I. In addition to a constraining sphere of radius $`2.25\sigma `$ to prevent evaporation, exchanges have been attempted every 10 passes, with the same method for choosing exchanging trajectories as described in the previous article. The simulations are begun with random configurations of the cluster geometry, and consist of $`1.3\times 10^{10}`$ points accumulated following equilibration moves consisting of $`95\times 10^6`$ Metropolis points(no exchanges) followed by $`190\times 10^6`$ points using parallel tempering. The microcanonical heat capacity calculated in this fashion and shown in Fig. 1, is qualitatively the same as the canonical heat capacity \[see I\]. The melting peak in the microcanonical heat capacity occurs at the same calcuated temperature as the temperature of the melting peak in the canonical heat capacity, and there are slope change regions at temperatures that correspond to equilibrium between the icosahedral basin and the truncated octahedral global minimum structure. The present simulations are also used to obtain structural insight about the cluster as a function of total energy. We have calculated the order parameter $`Q_4`$ as defined in I as a function of temperature, and compared the classification into the three categories of isomers (truncated octahedral, icosahedral or liquid-like) using the energy criterion also outlined in I.

In Fig. 2 we show the caloric curve $`T(E)`$ determined from our parallel tempering microcanonical MC simulations. We also present the canonical curve for comparison. The melting transition near $`T0.166\epsilon /k_B`$ is reflected in the change in slope of the temperature as a function of the energy. The microcanonical curve does not display a van der Waals loop, and remains very close to the canonical curve. The average value of the order parameter $`Q_4`$ is displayed in the lower panel of Fig. 2 as a function of the total energy. As has been discussed in I for the canonical simulation, the order parameter begins to drop at energies where there is the onset of isomerization transitions to the icosahedral basin (near $`E=160\epsilon `$), and the order parameter reaches its lowest value at the melting transition. The isomer distributions have been evaluated either using the parameter $`Q_4`$ or using the energy criterion (see the discussion in paper I). The results have been plotted in Fig. 3 as a function of the total energy. The behavior of the isomer distributions as a function of energy is similar to the canonical distributions as a function of temperature, and the cluster exhibits equilibrium between truncated octahedral and icosahedral geometries in the energy range $`160\epsilon E150\epsilon `$, prior to the solid-like to liquid-like phase change. As in the canonical case, the icosahedral distribution is a symmetric function of the energy when the energy criterion is used rather than the definition based on $`Q_4`$. This difference reflects the differences between the two definitions of icosahedral and liquid basins. The oscillatory structure observed at the peak of $`P_{Q_4}`$ for the icosahedral distribution in the upper panel of Fig. 3 is smaller than the calculated errors (two standard deviations of the mean are shown). Whether the observed structure would persist for a longer simulation is not known to us. Because the definition that assigns configurations to the icosahedral basin is arbitrary, we have chosen not to investigate this structure further.

It is useful to contrast the current results with previous constant energy studies of LJ38. Previous simulations have used molecular dynamics methods where no attempt has been made to insure ergodicity. To contrast these past studies with the molecular dynamics technique discussed in the next section of this paper, we define standard molecular dynamics to represent the usual molecular dynamics method where no special procedure is introduced to insure ergodicity. Simulations of LJ38 using standard molecular dynamics invariably lead to a caloric curve with a clear van der Waals loop and a melting temperature higher than that inferred from Fig. 2. From the results of Ref. , the cluster is trapped in the octahedral basin, and the system does reflect the true dynamical coexistence state between the truncated octahedron and the icosahedral basin. This is the common situation encountered in MD simulations of the LJ38 system; the cluster chooses either to remain trapped in the octahedral basin, or to escape and coexist between the icosahedral solid-like and liquid-like forms. Because the system is unable to return from the octahedral basin, the microcanonical temperature decreases. In the usual case, van der Waals loops arise when there are large energy gaps between the lowest-energy isomers. In the specific case of LJ38, it appears that the presence of extra (icosahedral) isomers only slightly higher in energy than the octahedral structure eliminates this loop in the ergodic microcanonical caloric curve.

In order to extract dynamical quantities, the Monte Carlo method we have presented must be modified to sample the MD ensemble. The modification is the subject of the next section.

## III Ergodic molecular dynamics

The molecular dynamics ensemble differs from the microcanonical ensemble in that two quantities are conserved in addition to the total energy $`E`$, volume $`V`$ and number of particles $`N`$. These two quantities are the total linear momentum $`𝐏`$ and total angular momentum $`𝐋`$. If their values are prescribed, the density of states remains the fundamental quantity of interest, and is now defined by

$`\mathrm{\Omega }(N,V,E,𝐏,𝐋)=`$ (8)
$`{\displaystyle \frac{1}{N!h^{3N}}}{\displaystyle \delta [H(𝐫,𝐩)E]\delta \left(𝐏\underset{i=1}{\overset{N}{}}𝐩_i\right)\delta \left(𝐋\underset{i=1}{\overset{N}{}}𝐫_i\times 𝐩_i\right)d^{3N}rd^{3N}p}.`$ (9)
As is the case in the microcanonical ensemble \[see Eq.(2)\], for atomic systems the momentum integrations in Eq.(9) can be evaluated explicitly. Because the thermodynamic properties are not affected by the translational motion of the center of mass, we can assume that $`𝐏=0`$. We then obtain

$`\mathrm{\Omega }(N,V,E,𝐏=0,𝐋)=`$ (10)
$`\left({\displaystyle \frac{2\pi m}{h^2}}\right)^{3N/23}{\displaystyle \frac{1}{N!\mathrm{\Gamma }(3N/23)}}{\displaystyle \mathrm{\Theta }[EU_𝐋(𝐫)][EU_𝐋(𝐫)]^{3N/24}\frac{d^{3N}r}{\sqrt{\mathrm{d}et𝐈}}},`$ (11)
where $`𝐈`$ is the inertia matrix and $`U_𝐋(𝐫)=U(𝐫)+𝐋^{}𝐈^1𝐋/2`$ is the effective rovibrational energy. This effective potential energy includes the kinetic energy contribution of the rotating cluster considered as a rigid body. The landscape of rotating clusters has been investigated by Miller and Wales in order to study cluster evaporation. Averages in the MD ensemble are now expressed as

$$A=\frac{{\displaystyle \mathrm{\Theta }[EU_𝐋(𝐫)][EU_𝐋(𝐫)]^{3N/24}A(𝐫)\frac{d^{3N}r}{\sqrt{\mathrm{d}et𝐈}}}}{{\displaystyle \mathrm{\Theta }[EU_𝐋(𝐫)][EU_𝐋(𝐫)]^{3N/24}\frac{d^{3N}r}{\sqrt{\mathrm{d}et𝐈}}}}.$$
(12)
As in the microcanonical ensemble, we define the entropy in the molecular dynamics ensemble by $`S=k_B\mathrm{ln}\mathrm{\Omega }`$. The differences between the microcanonical and molecular dynamics ensembles are the exponent $`3N/2`$ which is reduced by 3 owing to the geometrical constraints, the potential energy which now includes the contribution of the centrifugal energy, and the weight $`1/\sqrt{\mathrm{d}et𝐈}`$ which is a consequence of the conservation of the angular momentum. Monte Carlo simulations can sample the MD ensemble by performing a random walk in configuration space. The acceptance probability from configuration $`𝐫_o`$ to configuration $`𝐫_n`$ is

$$\mathrm{a}cc(𝐫_o𝐫_n)=\mathrm{min}(1,\frac{\rho _{E,𝐋}(𝐫_n)T(𝐫_n𝐫_o)}{\rho _{E,𝐋}(𝐫_o)T(𝐫_o𝐫_n)})$$
(13)
in the Metropolis scheme. The microcanonical weight $`\rho _E(𝐫)`$ is now replaced by the MD weight $`\rho _{E,𝐋}`$ given by

$$\rho _{E,𝐋}(𝐫)=\zeta ^1\frac{1}{\sqrt{\mathrm{d}et𝐈}}\mathrm{\Theta }[EU_𝐋(𝐫)][EU_𝐋(𝐫)]^{3N/24},$$
(14)
where $`\zeta `$ is a normalization. The expression for the acceptance probability is similar to Eq. (6), and a practical implementation of Monte Carlo in the MD ensemble is made in the same way as in the true microcanonical ensemble, given the vector $`𝐋`$. Parallel tempering can be also easily combined with the MC simulations. The acceptance probability of exchanging the configurations $`𝐫_i`$ and $`𝐫_j`$ initially at the total energies $`E_i`$ and $`E_j`$ respectively is then

$`\mathrm{min}(1,\left({\displaystyle \frac{[E_iU_𝐋(𝐫_j)][E_jU_𝐋(𝐫_i)]}{[E_iU_𝐋(𝐫_i)][E_jU_𝐋(𝐫_j)]}}\right)^{3N/24})`$

provided that all quantities inside brackets are positive (otherwise the move is rejected). It is remarkable that the geometrical weights have canceled in this expression.

The Monte Carlo method we have just described allows sampling of configuration space rigorously equivalent to the sampling we would obtain using molecular dynamics, but with the additional possibility of crossing the energy barriers higher than the available energy. The method can be used in its present form to extract equilibrium properties only dependent on the energy or geometry, as has been illustrated in the previous section. To compute dynamical quantities, the method can also provide a database of configurations representative of a given total energy. Instead of performing a few very long MD simulations that are in principle unable to reach other parts of the energy surface separated by barriers higher than the available energy, we choose to perform a statistical number of short simulations starting from configurations obtained by parallel tempering Monte Carlo in the MD ensemble with same total energy and angular momentum. By construction, if the MC method is correctly ergodic, then the hybrid MD method we have suggested can be expected to yield ergodic dynamical observables.

We now illustrate this ergodic molecular dynamics method on the LJ38 problem. Two essentially dynamical parameters have been calculated. The first is the self diffusion constant $`D`$, obtained from the derivative of the average mean square atomic displacement

$$D=\frac{1}{6}\frac{d}{dt}[𝐫(t)𝐫(0)]^2,$$
(15)
where the average is taken over all particles of the system and over all short MD simulations. The other parameter is the largest Lyapunov exponent $`\lambda `$, that measures the exponential rate of divergence of the distance between two initially close trajectories in the phase space. If we write the equation describing the Hamiltonian dynamics in condensed form as $`\dot{\psi }(t)=F(\psi )`$ where $`F`$ is a nonlinear function and $`\psi =\{𝐫,𝐩\}`$ the phase space point, then a small perturbation $`\delta \psi `$ evolves according to the simple equation $`d\delta \psi /dt=(F/\psi )\delta \psi `$. The largest Lyapunov exponent $`\lambda `$ is obtained from the time evolution of the vector $`\delta \psi `$:

$$\lambda =\underset{t\mathrm{}}{lim}\underset{\delta \psi (0)0}{lim}\frac{1}{t}\mathrm{ln}\frac{\delta \psi (t)}{\delta \psi (0)}.$$
(16)
In Eq.(16), $``$ is a metric on the phase space. In principle, any metric can be used, and we choose the Euclidian metric including both the momenta and the coordinates. The numerical procedure involves a periodic renormalization of the vector $`\delta \psi `$ to prevent its exponential divergence. The successive lengths are accumulated and contribute to the average value of $`\lambda `$.

In I, the clusters have been defined using a hard sphere constraining potential. Because the angular momentum is not conserved after reflection from such hard wall boundaries, in the molecular dynamics simulations we have chosen a soft repulsive spherical wall $`U_c`$ defined with respect to the center of mass of the cluster for each particle by

$$U_c(𝐫)=\{\begin{array}{cc}0,\hfill & r<R_c\hfill \\ \kappa (rR_c)^4/4,\hfill & rR_c\text{.}\hfill \end{array}$$
(17)
In this equation, the atomic distances $`r`$ are measured with respect to the cluster center of mass. The simulations have been performed setting the angular momentum to zero for simplicity. We stress that even in this case (with $`𝐋=0`$), the weight $`1/\sqrt{\mathrm{d}et𝐈}`$ must be included in the Monte Carlo probabilities so that we effectively sample the MD ensemble. The actual thermodynamic behavior in the MD ensemble at zero angular momentum is nevertheless nearly identical to the microcanonical behavior.

The application to the LJ38 cluster has been made by performing $`10^{10}`$ MC steps following $`10^7`$ equilibration steps in a parallel tempering simulation in the MD ensemble. The same 40 total energies have been chosen as in the previous section, and $`10^5`$ configurations have been stored every $`10^5`$ steps for each simulation. Short molecular dynamics runs of $`10^4`$ time steps following $`10^3`$ equilibration steps have been performed for each of these configurations, with the same total energy as the corresponding MC trajectory of origin, and with zero total linear and angular momenta as well. The parameters used for the constraining wall are respectively $`R_c=2.25\sigma `$ and $`\kappa =100\epsilon `$, for both the MC and MD runs. A simple Verlet algorithm has been used to propagate the MD trajectory with the time step $`\delta t=0.01`$ reduced LJ units. The propagation of the tangent trajectory to calculate the Lyapunov exponent has been determined with a fourth order Runge-Kutta scheme. The final values of $`D`$ and $`\lambda `$ are an average over the $`10^5`$ MD simulations. The variations of $`D`$ and $`\lambda `$ with total energy are depicted in Fig. 4. In both cases, two curves have been plotted, calculated either from standard molecular dynamics (with $`10^8`$ time steps following $`10^7`$ equilibration steps, and starting initially from the lowest-energy structure), or from our hybrid ergodic molecular dynamics method. For both quantities, the two MD schemes clearly yield distinct values in the energy range where equilibrium between truncated octahedral and icosahedral geometries occurs. The thermodynamic temperature, not plotted here, has the same variations as the caloric curve of Fig. 2 when calculated with ergodic MD. Standard molecular dynamics predicts a van der Waals loop centered at $`T0.18\epsilon /k_B`$. For standard MD, the cluster is trapped in the icosahedral basin and is, in practice, unable to reach the octahedral basin. Only the equilibrium between the icosahedral basin and liquid-like structures occurs. As can be seen from the upper panel of Fig. 4, this change in curvature of the temperature is also present for the diffusion constant, which exhibits strong variations at the energy where the octahedral structure vanishes when standard MD is used. In contrast, the variations in ergodic MD are smooth.

The melting temperature implied by the largest Lyapunov exponent is also higher in standard MD than in ergodic MD, even though the variations of the Lyapunov exponent are continuous in both MD schemes. Indeed, using ergodic molecular dynamics we observe a shift of the curve obtained by standard MD toward the lower energies. As shown by Hinde, Berry, and Wales, the Lyapunov exponent and the Kolmogorov entropy are quantities essentially dependent on the local properties of the energy landscapes. One contribution comes from the negative curvature of the landscape, and another contribution is the fluctuation of positive curvature. Both contributions are affected by the cluster being trapped either inside the truncated octahedral basin or inside the icosahedral basin. In this latter case in particular, the different isomers belonging to the icosahedral basin are connected through regions of negative curvature, while only one isomer defines the octahedral funnel.

Because ergodic molecular dynamics allows the cluster to be found in both basins prior to melting, the dynamical behavior is likely to be very different (and more chaotic) with respect to the dynamical behavior of the cluster confined to the octahedral funnel. This difference is precisely what we observe on the lower panel of Fig. 4.

## IV Conclusion

In this paper, we have explored the parallel tempering method in simulations in the microcanonical ensemble. The implementation of the parallel tempering algorithm in this ensemble is straightforward, the Boltzmann factor $`\mathrm{exp}(\beta U)`$ being replaced by the microcanonical weight $`(EU)^{3N/21}`$. Application to the LJ38 cluster has shown the thermodynamic behavior in the microcanonical ensemble to be similar to the behavior in the canonical ensemble. The solid-liquid phase change is preceded by a solid-solid phase change where the cluster is in equilibrium between truncated octahedral and icosahedral geometries. This phase equilibrium is well reproduced in the simulations owing to the power of parallel tempering. The calculated microcanonical caloric curve, which does not display a van der Waals loop, is consistent with the single peaked heat capacity observed in I.

We have extended the parallel tempering microcanonical Monte Carlo algorithm to sample the molecular dynamics ensemble at constant total energy, linear momentum and angular momentum. Combined with standard molecular dynamics, this method circumvents the lack of connectivity between regions of the potential energy surface. The method can ensure ergodicity in microcanonical simulations, which is much more difficult to achieve than in the canonical ensemble. Ironically, this ergodic MD method can be viewed as the counterpart of the techniques developed by Chekmarev and Krivov to study the dynamics of systems confined to only one catchment basin in the energy surface.

We have performed a statistical number of short molecular dynamics runs starting from configurations stored periodically in parallel tempering Monte Carlo simulations. These simulations sample the MD ensemble at the same total energies, linear and angular momenta as the standard molecular dynamics runs. In fact, the length of the MD runs is mainly dictated by the large number of starting configurations. One may think of reducing drastically this number, to allow for the calculation of parameters varying on longer time scales. Unfortunately, if the energy landscape is not known in advance, then it is hard to guess how important are the contributions of the basins not selected as starting configurations. In the case of LJ38 having only 3 main regions on the energy surface, one possibility is to compute a dynamical property as the average value over 3 different simulations starting either from the truncated octahedral geometry, one icosahedral geometry or a low-lying liquid geometry, all carried out at the same total energy. However, as we have seen in Fig. 3, it is not obvious how to choose properly the weights of each basin in this average because of the difficulty in distinguishing between icosahedral and liquid structures in many cases. For this reason, we believe that the first parallel tempering MC step of the hybrid ergodic method is essential in the vicinity of phase changes to capture many starting configurations that are used subsequently in standard molecular dynamics. The enhanced sampling offered by parallel tempering can also act as a statistical representation of the energy surface at a given total energy, and the long time dynamics may be further investigated by using master equations after searching the saddle points.

We have calculated two dynamical quantities with the present hybrid MD/MC method, the diffusion constant and the largest Lyapunov exponent in the 38-atom Lennard-Jones cluster. The variations of both quantities with the total energy are significantly different when evaluated with standard (non-ergodic) molecular dynamics or with our hybrid ergodic MD method. These results emphasize the different contributions of the two funnels of the energy landscape to the average value of the parameters estimated.

The algorithms developed in this investigation allow the calculation of thermodynamic, structural, or dynamical properties of systems such as LJ38 that can be expressed as phase space or time averages. Parallel tempering works using a criterion based on the potential energy but not on the geometry. Consequently permutational isomers can be introduced in the course of the simulation. Quantities such as fluctuations of configuration-dependent properties are much more difficult to extract than actual averages. For instance, the Lindemann index $`\delta `$, which measures the root mean square bond length fluctuation, is often considered to be a reliable parameter for detecting melting in atomic and molecular systems. This quantity cannot be properly estimated with the ergodic MD scheme, and the same difficulty persists for other methods based on the use of different trajectories.

Although the idea of combining Monte Carlo sampling with standard molecular dynamics can be applied to other techniques such as jump-walking, we believe that parallel tempering is the key to the success in the case of LJ38. As in the canonical version, the equilibrium phase between truncated octahedral and icosahedral structures is correctly reproduced in an energy range preceding the melting region, because in this range configurations may be accessed either from higher energy trajectories containing mainly icosahedral geometries, or from lower energy trajectories acting as a reservoir for the octahedral geometry. As noticed by Falcioni and Deem, the parallel tempering algorithm is especially useful at low temperatures, or in our case, at low energy. The long relaxation times inherent in systems like clusters, proteins, critical or glassy liquids, are a serious difficulty for standard simulation methods. We expect the present ergodic method to be particularly useful to deal with the dynamics of such systems.

The method we have presented works at constant total energy. It is possible to improve ergodicity in constant-temperature MD either by using canonical parallel tempering as in the work of Sugita and Okamoto, or by coupling parallel tempering canonical Monte Carlo to short Nosé-Hoover trajectories. In the Nosé-Hoover approach such molecular dynamics simulations do conserve a zero angular momentum, so a rigorous MC sampling should include the geometrical weight $`1/\sqrt{\mathrm{d}et𝐈}`$ in the probabilities also in this case. The present microcanonical scheme can be easily used for rotating bodies, which makes the method suitable for investigating the strong influence of centrifugal effects on phase changes in atomic clusters.

## Acknowledgments

Some of this work has been motivated by the attendance of two of us (DLF and FC) at a recent CECAM meeting on ‘Overcoming broken ergodicity in simulations of condensed matter systems.’ We would like to thank CECAM, J.E. Straub and B. Smit who organized the meeting, and those who attended the workshop for stimulating discussions. This work has been supported in part by the National Science Foundation under grant numbers CHE-9714970 and CDA-9724347. This research has been supported in part by the Phillips Laboratory, Air Force Material Command, USAF, through the use of the MHPCC under cooperative agreement number F29601-93-0001. The views and conclusions contained in this document are those of the authors and should not be interpreted as necessarily representing the official policies or endorsements, either expressed or implied, of Phillips Laboratory or the US Government.