text:
# THE DYING-WIND AROUND HD 56126, A POST-AGB CARBON STAR

## 1 INTRODUCTION

Mass loss with rates near $``$ 10-5 M yr-1 during the Asymptotic Giant Branch (AGB) phase is important both in the evolution of the mass-losing star and in replenishing the interstellar medium (see, for example, Habing 1996). After a star evolves off the AGB, its effective temperature rises, its radius shrinks and the wind dies (Vassiliadis & Wood 1993, Blocker 1995, Schroder, Winters & Sedlmayr 1999).

Stars slightly beyond the end of their AGB evolution can be identified by having only relatively cold dust since they are no longer losing much mass. Therefore, they exhibit a large infrared excess but with dust temperatures typically less than 200 K (see, for example, Kwok 1993). The detached dust shell can be extended over the sky and thus provide the opportunity to examine the history of the star’s mass loss.

Meixner et al. (1997, 1999) and Dayal et al. (1998) have acquired mid-infrared images and resolved more than 10 post-AGB stars or “transition objects”. Even in the best-studied cases, the angular diameter of the dust emission is not much larger than the characteristic $`0\stackrel{}{\mathrm{.}}9`$ angular resolution available at the 4m-class telescopes used for these previous studies. With the 10m Keck telescope, it is possible to achieve significantly higher angular resolution (see, for example, Jura & Werner 1999, Jura, Chen & Werner 2000) and thus develop more exacting tests of models for dying winds.

Here, we report the result of using Keck to obtain mid-IR images of HD 56126 (mV = 8.27 mag, spectral type F5I), a post-AGB carbon star with a nebula displaying one of the largest angular sizes in the sample studied by Meixner et al. and extended optical reflection nebulosity as well (Ueta, Meixner & Bobrowsky 2000). This star is also a “21 $`\mu `$m” source (see, for example Kwok, Volk & Hrivnak 1999), a feature attributed to TiC (von Helden et al. 2000). In Section 2 we present our observations while in Section 3, we describe our estimates for the physical properties of the circumstellar envelope and the central star. In section 4, we present a model to explain the gross features of the system while in Section 5, we discuss models to explain deviations from a spherically symmetric wind. In Section 6 we present our conclusions.

## 2 OBSERVATIONS

Our data were obtained on the night 2000 Feb 20 (UT) at the Keck I telescope using the Long Wavelength Spectrometer (LWS) which was built by a team led by B. Jones and is described on the Keck web page. The LWS uses a 128 $`\times `$ 128 SiAs BIB array with a pixel scale at the Keck telescope of 0$`\stackrel{}{\mathrm{.}}`$08 and a total field of view of 10$`\stackrel{}{\mathrm{.}}`$2 $`\times `$ 10$`\stackrel{}{\mathrm{.}}`$2. We used the “chop-nod” mode of observing, and 4 different filters: 7.5 - 8.2 $`\mu `$m, 9.4-10.2 $`\mu `$m, 11.2 - 12.2 $`\mu `$m and 18.4-18.9 $`\mu `$m. We used $`\alpha `$ Boo for flux and point-spread-function calibrations. The data were reduced at UCLA using standard LWS routines.

The image in the 11.7 $`\mu `$m filter is presented in Figure 1; the images at other wavelengths appear similar. Probably because of our higher angular resolution, our measured peak surface brightness of 4.0 Jy arcsec-2 is greater than the peak values of 2.3 Jy arcsec-2 and 3.0 Jy arcsec-2 found by Meixner et al. (1997) and Dayal et al. (1998), respectively. Given the $``$10% uncertainties in our data, our measured flux at 11.7 $`\mu `$m of 28 Jy is consistent with the values for the flux of 30 $`\pm `$ 3 Jy (Dayal et al. 1998), 25.7 $`\pm `$ 2.6 Jy (Meixner et al. 1997) 25 $`\pm `$ 3 Jy (Justtanont et al. 1996) and 25 Jy (IRAS).

Beyond 1$`\stackrel{}{\mathrm{.}}`$5 from the center of the nebula, the isophotes are smooth and elliptically shaped with the major axis lying near Position Angle 45 and being about 7% greater than the minor axis. A hint of ellipticity at this position angle is also seen in the mid-IR images of Meixner et al. (1997) and Dayal et al. (1998) as well as in the optical reflection nebulosity presented by Ueta et al. (2000). In the inner 1$`\stackrel{}{\mathrm{.}}`$5 there is a marked change in the morphology of the dust with two bright blobs lying near the minor axis of the ellipse defined by the outer contours. The bright blob at about 1″ to the southeast of the center of the image is seen in the data presented both by Dayal et al. (1998) and Meixner et al. (1997). The maximum intensity in the bright blob is about a factor of 1.5 greater than the maximum intensity in the faint blob. Our image at 8$`\mu `$m shows that HD 56126 lies within $`\pm `$ 0$`\stackrel{}{\mathrm{.}}`$2 of the center of the dust ring defined by the outer, elliptical isophotes.

## 3 STELLAR AND CIRCUMSTELLAR PARAMETERS

The observed pulsation properties of HD 56126 can be used to place the star in the H-R diagram and to estimate its mass. Since HD 56126 is variable with a period of 36.8 days (Barthes et al. 2000), it is plausible that the mass of the star is near 0.6 M and the time-averaged luminosity is near 6600 L (Tuchman et al. 1993, Jeannin et al. 1997, Barthes et al. 2000). From the observed time-averaged flux (see Meixner et al. 1997), we therefore infer a distance of 2300 pc. The time-averaged effective temperature, $`T_{}`$, and radius, $`R_{}`$, are 5900 K and 5.5 $`\times `$ 1012 cm, respectively (Barthes et al. 2000). Although carbon, nitrogen and oxygen have nearly solar abundances, since the iron group elements have abundances about 0.1 of solar (Parthasarathy, Garcia Lario & Pottasch 1992, Klochkova 1995, Van Winckel & Reyniers 2000), HD 56126 probably is a member of the thick disk and its main sequence progenitor probably had a mass near 1.1 M. Therefore, it is likely that HD 56126 has already experienced a total mass loss of 0.5 M. Evidence that HD 56126 is a post-AGB star is provided by the greater than solar abundance of the s-process elements such as Sr, Y, Zr and Ba.

Below, we propose that $``$1500 years ago, HD 56126 was losing $``$3 $`\times `$ 10-5 M yr-1. In this phase, its effective temperature was probably near 2600 K (Arndt et al. 1997). Therefore, during the past 1500 yr, the average value of $`dT_{}/dt`$ has been +2.2 K yr-1. The star probably maintained a constant luminosity during its post-AGB evolution (see, for example, Vassiliadis & Wood 1993) so that as the temperature increased, the radius decreased from 2.8 $`\times `$ 1013 cm to its current value.

Although, as inferred from the profile of the CO (2-1) circumstellar emission line, most of the gas outflow is at 10.7 km s-1 (Knapp et al. 1998), there are significant deviations from this value. Crawford & Barlow (2000) have obtained optical spectra with a resolution of 860,000, and they found three distinct components in the K I circumstellar shell lines with a net spread in the velocity of the circumstellar gas of 2 km s-1.

## 4 THE SPHERICALLY SYMMETRIC WIND

To compute the surface brightness, $`I_\nu `$, of the nebulosity, we numerically determine:

$$I_\nu =\chi _\nu \rho (D)B_\nu (T_{gr})𝑑s$$
(1)
where $`B_\nu `$ is the Planck function, $`ds`$ is the increment of distance along the line of sight, $`\rho (D)`$ is the dust density at distance, $`D`$, from the star and the dust opacity is $`\chi _\nu `$. To compute the dust temperature as a function of the distance, we follow Sopka et al. (1985) and assume that the dust grains are heated by the light from the central star and re-radiate in the infrared. For carbon-rich particles, we adopt a simple model for the grain opacity, $`\chi _\nu `$ (cm2 g-1) such that

$$\chi _\nu =\chi _{\nu _0}(\nu /\nu _0)^{1.3}$$
(2)
with $`\chi _{\nu _0}`$ = 150 cm2 g-1 at $`\nu _0`$ = 5 $`\times `$ 1012 Hz (60 $`\mu `$m) (Le Bertre 1997). When the circumstellar envelope is optically thin, the dust grain temperature, $`T_{gr}`$, can be written as:

$$T_{gr}=\mathrm{\hspace{0.33em}0.77}T_{}(R_{}/D)^{0.38}$$
(3)
With this prescription, $`T_{gr}`$ at 1$`\mathrm{}`$ from the central star is 170 K while Meixner et al. (1997) and Dayal et al. (1998) inferred values of 160 K and 185 K, respectively. For 11.7 $`\mu `$m photons, $`h\nu /k`$ = 1230 K, and therefore the nebulosity is produced by grains emitting on the Wien portion of the Planck curve.

To compute $`\rho (D)`$, we assume a spherically symmetric envelope where the dust outflow velocity, $`V_{gr}`$, does not change with time. We assume that the mass loss rate of grains, $`\dot{M}_{gr}(t)`$, and $`T_{}`$ were constant until a moment came when the star’s effective temperature started to increase and the wind started to die. During the dying-wind phase, we assume that $`dT_{}/dt`$ was constant. The current boundary between the zone where the wind had a constant mass loss rate and where it started to die, lies at a distance $`D_{die}`$ with $`D_{die}`$ = $`V_{gr}`$ $`t_{die}`$. Here, $`t_{die}`$ denotes the amount of time that has elapsed since that moment when the star’s effective temperature began to increase. From the equation of continuity and the simple mapping between the dust travel time, $`t_{travel}`$, and distance from the star, then:

$$\rho _{gr}(D)=\dot{M}_{gr}(t_{travel})/(4\pi V_{gr}D^2)$$
(4)
In the models for mass loss from carbon-rich AGB stars by Arndt et al. (1997), $`\dot{M}_{gr}`$ varies as $`T_{}^{8.26}`$. Such a steep variation of the mass loss with temperature occurs because the mass loss rate is controlled by radiation pressure on grains, and the dust formation is very sensitive to the gas temperature.

The best fit to the 11.7 $`\mu `$m image shown in Figure 1 is with a model derived from equation (1) where $`D_{die}`$ = 5 $`\times `$ 1016 cm, corresponding to an angular radius of 1$`\stackrel{}{\mathrm{.}}`$5 or $`t_{die}`$ = 1500 yr, and $`dT_{}/dt`$ = +2.2 K yr-1. We show in Figure 2 the results of comparing the results from this model with transverse cuts through the intensity data at 11.7 $`\mu `$m. In the region where $`D`$ $`<`$ $`D_{die}`$, the projected 11.7 $`\mu `$m intensity depends upon the multiplication of two rapidly varying quantities: the mass loss rate and the Wien portion of the Planck function. The qualitative agreement between our data and the observations lends support to the dying wind model, but other possible models for $`\dot{M}(t)`$ are not excluded by our observations.

At 1$`\stackrel{}{\mathrm{.}}`$5 from the star, the 11.7 $`\mu `$m intensity of the nebula is about 2.0 Jy arcsec-2. Therefore, in the outer region, $`\dot{M}_{gr}`$ = 1.7 $`\times `$ 10-7 M yr-1. Although uncertain, a dust to gas ratio of 5 $`\times `$ 10-3 is plausible and therefore the total mass loss rate before the wind started to die was $``$3 $`\times `$ 10-5 M yr-1. Given the uncertainties in the distance to the star, the dust to gas ratio, the emissivity of the grains and the observations that there is an appreciable amount of carbon in C I as well as CO (Bakker et al. 1996, Bakker & Lambert 1998, Knapp et al. 2000), this estimate of the mass loss rate for the phase just before the winds dies is consistent with previous estimates by Hrivnak et al. (1989), Meixner et al. (1997) and Dayal et al. (1998).

## 5 DISCUSSION

As a first approximation, much of our data can be explained in terms of a nearly-spherical wind of $``$3 $`\times `$ 10-5 M yr-1 which started to die about 1500 yr ago, consistent with models for post-AGB evolution of carbon-rich stars (Vassiliadis & Wood 1993, Blocker 1995, Schroder et al. 1999).

In order to produce a high enough pressure in the outflow to manufacture TiC, von Helden et al. (2000) proposed a short-lived phase with a steady-state mass loss rate of 10-3 M yr-1 which is not supported by our data. Nevertheless, models for shock-driven mass loss from AGB stars may possibly explain the presence of TiC around HD 56126. In calculations by Winters et al. (1997) for a carbon star with a time-averaged mass loss rate of 1.2 $`\times `$ 10-4 M yr-1, as shown in their Figure 2, at 2 stellar radii above the star, the density can be as high as 10-11 g cm-3 when the gas temperature is 1800 K. Therefore, the gas pressure can be $``$1 dyne cm-2 at a temperature nearly low enough to produce solid TiC (Bernatowicz et al. 1996). There may have been phases during the pulsational cycle of HD 56126 while it was on the AGB when the physical conditions were favorable for the production of the apparent carrier of the 21 $`\mu `$m emission feature even if the mass loss rate was as “low” as 3 $`\times `$ 10-5 M yr-1.

The nebula shows significant deviations from spherical symmetry. Because the mass loss rate is extremely sensitive to the photospheric temperature, we tentatively explore the possibility that deviations from spherical symmetry in the dust shell can be explained by temperature variations on the surface of the star while it was losing mass.

### 5.1 The Outer Ellipticity

As can be seen in Figure 1, beyond 1$`\stackrel{}{\mathrm{.}}`$5 from the center, for a given distance from the center, the intensity at position angle 45 is greater than at, say, position angle 135. From equation (1), if the grain opacity is the same throughout the nebula, then these data can be reproduced either by variations in $`\rho `$ or by variations in $`T_{gr}`$. The 18.7 $`\mu `$m image displays the same ellipticity as does the 11.7 $`\mu `$m image so that azimuthal variations in the temperature do not seem to explain the variations in the intensity. Instead, we propose that there are more grains at position angle 45 than at position angle 135 and thus that there was more mass loss directed toward position angle 45.

If the star is rotating and there is enhanced mass loss at the equatorial bulge (Dorfi & Hofner 1996, Soker 1998, Reimers, Dorfi & Hofner 2000), then qualitatively, the data can be explained if the rotation axis is projected along position angle 135. If $`\dot{M}`$ varies as $`T_{}^{8.26}`$ then a 1% variation in the temperature can lead to an 8% variation in the mass loss rate – a change great enough to explain the observed 7% ellipticity in the isophotes.

If the surface of the star is an equipotential, if the star rotates as a solid body with angular velocity $`\omega _{}`$, and if $`\theta `$ denotes the co-latitude, then:

$$(1/2)\omega _{}^2R_{}(\theta )^2sin^2\theta +GM_{}/R_{}(\theta )=constant$$
(5)
When the star was on the AGB with $`T_{}`$ = 2600 K, $`\overline{R}_{}`$ = 2.8 $`\times `$ 1013 cm and $`M`$ = 0.6 M, and if $`R_{eq}/R_{pole}`$ was $``$1.02, consistent with a temperature contrast between the poles and equator of 1.01 and the assumption that $`R^2T_{}^4`$ is constant over the surface of the star, then $`\omega _{}`$ was approximately 1 $`\times `$ 10-8 s-1. The angular momentum of the rotating AGB star was approximately 2/9 $`M_{}R_{}^2\omega _{}`$ (Soker 1998) or 2 $`\times `$ 1051 g cm2 s-1. It is plausible that HD 56126 could have acquired this angular momentum if, when it first swelled onto the AGB, it engulfed a putative 0.1 M companion in a circular orbit at 1 AU which would then have donated 9 $`\times `$ 1051 g cm2 s-1 of angular momentum. An example of a carbon star that has probably engulfed a companion and acquired such a large amount of angular momentum is V Hya (Barnbaum, Morris & Kahane 1995, Kahane et al. 1996). However, we have no evidence in support of the proposed magnitude or orientation of the rotation of HD 56126, and our hypothesis to explain the observed ellipticity is quite tentative.

### 5.2 The Inner Blobs

Meixner et al. (1997) propose that the bright blobs are the consequence of an equatorially-enhanced ejection. Their model implies that the position angle of the rotation axis lies near 45; a very different orientation from what we propose. One difficulty with the model by Meixner et al. (1997) is that the southeastern blob is considerably brighter than the northwestern blob. However, it is possible that the bright blob is a factor of 1.2 closer to the central star than the faint blob, and this could explain the intensity difference between the two blobs. An offset of the central star from the center of the nebula could occur if the system is a wide binary as apparently occurs in TT Cyg (Olofsson et al. 2000).

If, as we propose above, the position angle of the rotation axis is 135, then another model besides that of Meixner et al. (1997) is required to explain the existence of the bright blobs. One possibility is that as the star first evolved away from the AGB, the magnetic field was strongest at the rotational poles with the consequence that the gas was particularly cold in these regions. We speculate that the mass loss continued from the poles for 500 years after it had greatly diminished from the equator. If the field were stronger at one pole than the other, then the two blobs in the nebula could have different intensities. To appreciably change the surface temperature of the star, the energy density in the magnetic field, $`B^2/(8\pi )`$ must be comparable to the thermal pressure in the photosphere (Soker & Clayton 1999). The required value of the magnetic field is therefore $`[(16\pi g)/(3\chi )]^{1/2}`$ where $`\chi `$ is the Kramers opacity in the atmosphere or $``$ 3 $`\times `$ 10-3 cm2 g -1 (Soker 1998) and $`g`$ is the gravitational acceleration which was near 0.1 cm s-2 when the star was just leaving the AGB. Thus a field of $``$20 Gauss is implied. From VLBA maps of polarized SiO masers, magnetic fields approaching this amplitude have been inferred around the oxygen-rich Mira variable TX Cam (Kemball & Diamond 1997), but the true value of the magnetic field is controversial and may be only 0.03 Gauss (Wiebe & Watson 1998, Desmurs et al. 2000).

## 6 CONCLUSIONS

We have obtained high-resolution mid-IR images of the detached dust shell around the carbon star HD 56126 and propose the following:

1. The gross features of the region beyond a radius of 1$`\stackrel{}{\mathrm{.}}5`$ can be reproduced by an approximately spherical wind from a star which had a mass loss rate of $``$3 $`\times `$ 10-5 M yr-1 while it was on the AGB. After the star left the AGB, its effective temperature increased with an average rate of +2.2 K yr-1 and the mass loss rate varied as $`T_{}^{8.26}`$, as expected in models for dust-driven winds.

2. The extreme sensitivity of the mass loss rate to the effective temperature of the star allows us to tentatively propose that both the ellipticity in the outer envelope and the inner blobs can be explained by plausible temperature variations on the surface of the star which might have resulted from rotation and/or large scale magnetic fields.

We thank M. Meixner for a particularly thoughtful referee’s report and M. Morris for his comments. This work has been partly supported both at UCLA by NASA and at the Jet Propulsion Laboratory, California Institute of Technology, under an agreement with NASA.

FIGURE CAPTIONS

Fig 1. The 11.7 $`\mu `$m image of HD 56126. North is to the top and and East is to the left; these directions are displayed on the Figure by lines whose lengths are 1$`\mathrm{}`$ each. The maximum surface brightness that we detected and shown in the red color is 4.0 Jy arcsec-2. The yellow, green and purple colors correspond to 2.4, 1.6 and 0.40 Jy arcsec-2, respectively.
Fig. 2. Two cuts through the 11.7 $`\mu `$m image of HD 56126. Red and green crosses are used for tracing the images at position angles of 135 and 45, respectively, while the results from the spherically symmetric model are plotted as a solid blue line.


format: text

text:
# Fluctuation Effects in High Sheet Resistance Superconducting Films11footnote 1To appear in the Russian review journal Uspekhi

## Abstract

As the normal state sheet resistance, $`R_n`$, of a thin film superconductor increases, its superconducting properties degrade. For $`R_nh/4e^2`$ superconductivity disappears and a transition to a nonsuperconducting state occurs. We present electron tunneling and transport measurements on ultrathin, homogeneously disordered superconducting films in the vicinity of this transition. The data provide strong evidence that fluctuations in the amplitude of the superconducting order parameter dominate the tunneling density of states and the resistive transitions in this regime. We briefly discuss possible sources of these amplitude fluctuation effects. We also describe how the data suggest a novel picture of the superconductor to nonsuperconductor transition in homogeneous 2D systems.

Introduction

It has been known for some time that increasing the normal state sheet resistance, $`R_n`$, of a homogeneously disordered (as opposed to islanded) superconducting film reduces its mean field transition temperature, $`T_{co}`$ and superconducting energy gap, $`\mathrm{\Delta }`$. For $`R_n`$ comparable to $`h/4e^2`$, $`T_{co}`$  and $`\mathrm{\Delta }`$ approach zero and a superconductor to non-superconductor transition (SNST), commonly referred to as a superconductor to insulator transition (SIT), takes place (see Fig. 1). Close to this transition, variations in the density of electronic states on mesoscopic length scales, disorder enhanced inelastic scattering processes and proximity to the quantum critical point separating a superconducting from a non-superconducting phase are expected to cause large, non-mean field like fluctuations in the amplitude of the superconducting order parameter. Finding evidence of these effects is interesting because it can provide insight into the nature of the non-superconducting phase and the physics driving the SNST.

In this paper, we present tunneling and transport data on homogeneous films in this regime that show evidence of large amplitude fluctuations. The tunneling measurements indicate that the fluctuations increase the quasiparticle density relative to the cooper pair density or equivalently, increase the normal fluid fraction. The transport measurements suggest that the paraconductance or amplitude fluctuation dominated regime of the resistive transitions grows and diverges in size at the SNST (SNST). Taken together, the experiments intimate that amplitude fluctuation effects drive down the superfluid density, by increasing the quasiparticle density to drive the SNST in homogeneous films and the nonsuperconductor consists of weakly localized quasiparticles. This scenario contrasts sharply with the prevalent dirty boson models of a SIT, which attribute the reduction in the superfluid density to the reduction of Josephson coupling between spatially separated film regions.

Experimental Details

The quench condensed films that we have studied are homogeneous Bi/Sb and PbBi/Ge. They were deposited on fire-polished glass substrates held near 4 K on the cold stage of a dilution refrigerator . They consist of a two to three monolayer thick in situ deposited amorphous layer of the Sb or Ge that does not conduct at low temperatures covered by subsequent evaporations of the Bi or PbBi. Experiments  indicate that these films have a uniform morphology. First, they become electrically continuous after deposition of as little as 1-2 bulk monolayers worth of the metal as measured by a quartz microbalance. They superconduct at 3-4 bulk monolayers. Second, their conductance as a function of thickness indicates that the elastic mean free path is on the order of an interatomic spacing and changes very little with metal film thickness. Finally, in situ scanning tunneling microscopy measurements on a similar system imply that the films have a smooth, rather than a granular or islanded, morphology. The transport and tunneling measurements that we present here were obtained using standard, four terminal, low level AC and DC techniques. For a BCS s-wave superconductor, the normalized differential tunnel junction conductance gives a direct measure of the superconducting quasiparticle density of states, $`N_s(E)=E/\sqrt{E^2\mathrm{\Delta }^2}`$ where $`E`$ is measured relative to the Fermi energy and $`\mathrm{\Delta }`$ is the superconducting energy gap. At finite temperatures, the voltage dependent conductance qualitatively resembles the density of states. At low reduced temperatures, $`t<.5`$, the voltage at which the differential conductance crosses 1 corresponds closely to $`\frac{\mathrm{\Delta }}{e}`$.

Results

Tunneling and transport measurements on homogeneous films strongly suggest that the transition out of the superconducting state in these films involves the collapse of the amplitude of the superconducting order parameter and a concomitant rise in the quasiparticle density. Experiments have shown that on approaching the transition region from the superconducting side, the mean field transition temperature decreases and simultaneously, the size of the energy gap decreases. In fact, very close to the transition, the energy gap appears to go to zero . Since the energy gap gives a direct measure of the order parameter amplitude, these earlier tunneling measurements intimate that the amplitude of the superconducting order parameter goes to zero as superconductivity disappears.

In addition, as $`\mathrm{\Delta }`$ decreases, the tunneling density of states assumes an increasingly broadened BCS form as shown by Fig. 2. The differential tunnel junction conductances for these two PbBi/Ge films were obtained at a reduced temperature of $`t`$ 0.2 and normalized by the voltage dependent normal state conductance (see ref for procedure). The voltage dependence in the normal state results from disorder enhanced coulomb interaction effects and is well understood. Qualitatively, the tunnel junction conductances of both films exhibit the expected BCS form including a peak and a gap. Quantitatively, the film with the larger gap fits the BCS form well. However, the conductance of the junction on the film with the smaller energy gap, which is closer to the transition region, is higher in the gap region and the peak is substantially lower than expected for a BCS superconductor. The inset of Fig. 2, in which the same data are plotted on a voltage scale normalized to the peak in the conductance, emphasizes these deviations. For BCS superconductors with different gaps, conductances measured at the same reduced temperature should fall upon one another when the voltage axis is scaled in this fashion. Clearly, the peak is smaller and relatively broader and the conductance in the “gap” is substantially higher than BCS theory predicts. These qualitative deviations from BCS behavior continue to grow upon approaching the transition region .

Fluctuation effects are also apparent in the tunneling and transport in homogeneous films at temperatures in the vicinity of $`T_{co}`$. Fig.3 shows the zero voltage bias conductance (ZBC) of a tunnel junction on and the resistance of thin homogeneous PbBi/Ge films as a function of temperature. To facilitate comparison of the tunneling and transport data among different films, the ZBC’s are normalized to their values at 8K and the resistances are normalized to their maximum value. The ZBC is proportional to the number of quasiparticle states within 4.8kBT of the Fermi energy. In bulk conventional superconductors, the slope of the ZBC as a function of temperature changes discontinuously at $`T_{co}`$ where $`\mathrm{\Delta }`$ appears and the ZBC decreases rapidly with decreasing temperature. The film with the highest $`T_{co}`$ in Fig. 3 exhibits this behavior. Also, the resistive transition for this film is sharp and occurs with the kink in the ZBC. As $`T_{co}`$ decreases the kink in the ZBC rounds and its drop becomes more gradual. Concomitantly the resistive transition, which occurs near the same temperature as the gap opening, becomes broader. While some broadening near $`T_{co}`$ and above is expected, the broadening near the transition region is greater than conventional expectations (see below). It is important to note that for each of the three films, the opening of the energy gap and the resistive transition midpoint, i.e. $`T_{co}`$, roughly coincide. This coincidence implies that the quasiparticle density is high at and above $`T_{co}`$. Correspondingly, the order parameter amplitude is very small and susceptible to thermal flutuations in the region where the broadening of the resistive transitions occurs. In other words, this broadening of the resistive transitions corresponds to a growth of the paraconductance or amplitude fluctuation dominated portion of the superconducting transition.

Analysis of the transport data shown in Fig. 1 suggests that the size of the paraconductance fluctuation dominated regime diverges as $`T_{co}`$ decreases. To accurately measure the size of this regime it is necessary to know the temperature dependence of the normal state conductance. For ”clean” materials this dependence is very weak and often negligible. In these strongly disordered films the temperature dependence is too strong to ignore, but well understood. It stems from weak localization and electron-electron interaction effects and these effects lead to a temperature dependence of the conductance that varies little with film sheet resistance. Fig. 4 demonstrates this weak dependence, showing the conductances of a few Bi/Sb films near the transition between the superconducting and non-superconducting states. The lowest conductance film exhibits no signs of superconductivity while the two upper films both show a conductance that increases with decreasing temperature at the lowest temperatures. We use the temperature dependence of the conductance of the lowest conductance film shown to determine the size of the paraconductance fluctuations, $`\mathrm{\Delta }G_{sc}`$, as shown in Fig. 4. Normalized resistive transitions, obtained using $`\mathrm{\Delta }G_{sc}`$ (see figure caption), are shown in Fig. 5a for films near the transition. Clearly, the relative breadth of the superconducting transitions grows as $`T_{co}`$ decreases. In the case of the film closest to the transition whose resistance drops by less than 10% down to 0.1 K, superconducting fluctuation effects are apparent over a factor of more than 30 in temperature! More quantitatively, the widths corresponding to the temperature interval over which the resistance drops from 90% to 10% of its normal state value, normalized by $`T_{co}`$ are plotted in Fig. 5b. $`T_{co}`$ is estimated by the temperature at which the film sheet resistance has fallen to half its normal state value. As shown, the width appears to diverge logarithmically as $`T_{co}`$ goes to zero.

Discussion

The characteristics of the deviations of the tunneling conductances and resistively measured transitions from the simple mean-field, BCS form provide strong evidence that fluctuations in the amplitude of the superconducting order parameter grow upon approaching the SNST in homogeneous films. Near the transition, the tunneling data cannot be fit to the simple BCS form using a single energy gap: the peaks are too broad and there is a finite number of states in the gap region indicating states at the Fermi energy in the superconducting state. Taken together, these imply that the amplitude of the order parameter, which is proportional to $`\mathrm{\Delta }`$, is not well defined and may be zero in certain regions of a film. The growth and apparent divergence of the width of the resistive transitions over the temperature interval where tunneling indicates that the energy gap has barely opened and hence, the order parameter amplitude is small, intimates that fluctuations in the amplitude of the superconducting order parameter become critical at the SNST.

At the high sheet resistances considered here, $`R_nR_Q`$, there are a few factors that could contribute to the increasing breadth of the density of states, the broadening of the ZBC temperature dependence, and the broadening of the resistive transitions. First, inelastic quasiparticle scattering processes become stronger with increasing sheet resistance and these processes reduce the quasiparticle lifetime. The reduced lifetimes broaden the structure in the quasiparticle density of states and thus, the tunnel junction conductance. They also increase the density of low energy excitations which increases the ZBC and broaden its temperature dependence. While such inelastic scattering effects are in the nonperturbative regime at these sheet resistances, they should still disappear at low temperatures. Thus, they probably cannot explain all of the broadening indicated in Fig. 3 because those data were taken at a relatively low reduced temperature of .2.

There are two other factors that could be responsible for the low temperature broadening. Mesoscopic fluctuations in the density of electronic states near the Fermi energy could give rise to fluctuations in the amplitude of the pairing potential and the energy gap that would persist in the low temperature limit. For $`R_nR_Q`$ these fluctuations are expected to be of order unity. Also, the SNST could be a quantum phase transition with associated critical point amplitude fluctuation effects. Recent theories (see for example ) predict the growth of fluctuations in the superconducting order parameter amplitude that diverge at the transition. The conductance of a large area tunnel junction measures a density of states averaged over these fluctuations. It should be broader than the density of states for a BCS superconductor with a single gap. The resistively measured transition would be sensitive to the accompanying spatial distribution of $`T_{co}`$’s in the film and would also become broader. In two dimensions, the predicted divergence for a superconductor to normal metal transition is logarithmic, in accord with the result shown in Fig. 5b.

The above results and discussion suggest a qualitatively different picture of the SNST for homogeneous films than the “standard” dirty Boson models of the SIT (see for example ). In the latter, increasing the normal state sheet resistance decreases the Josephson coupling between individually superconducting islands in the film. The concomitant reduction in the superfluid density or phase modulus leads to the growth of phase fluctuations between the islands. If the phase fluctuations are large enough, all phase coherence is lost, Cooper pairs become localized to their individual islands, and the system becomes an insulator. In contrast, in the homogeneous films considered here, the growing amplitude fluctuations increase the quasiparticle density. These increases come at the expense of the Cooper pair or superfluid density. At high enough sheet resistances the superfluid density is driven to zero and a nonsuperconducting phase consisting of weakly localized quasiparticles takes over.

Summary

We have presented electron tunneling and transport measurements on homogeneous films that provide evidence that fluctuations in the superconducting to order parameter amplitude grow as the amplitude collapses at the SNST. These fluctuation effects dominate the experimentally accessible portion of the resistive transitions of the films closest to the SNST. Theories suggest that the mesoscopic fluctuation, quasiparticle lifetime, and quantum critical fluctuation effects could each contribute to the growth of the amplitude fluctuations. Taken as a whole, our results suggest that amplitude fluctuations become so strong near $`R_nR_Q`$ that they drive the superfluid density to zero to destroy the superconducting state.

We wish to acknowledge the support of NSF grants DMR-9801983 and DMR-9502920 and helpful conversations with Fei Zhou, Brad Marston, Boris Spivak, Nandini Trivedi, Sean Ling, and Dietrich Belitz.


text:
# Poissonian reducibility and thermal scaling within the lattice gas model and molecular dynamics model

## Abstract

The emission of clusters in the nuclear disassembly is investigated within the framework of isospin dependent lattice gas model and classical molecular dynamics model. As observed in the recent experimental data, it is found that the emission of individual cluster is poissonian and thermal scaling is observed in the linear Arrhenius plots made from the average multiplicity of each cluster. The mass, isotope and charge dependent ”emission barriers” are extracted from the slopes of the Arrhenius plots and their possible physical implications are investigated.

In low-intermediate energy heavy ion collisions (HIC) the hot nuclei with moderate temperature can be formed and they finally deexcite by the different decay modes, such as the emission of multiple intermediate mass fragment ($`IMF`$), $`i.e.`$ multifragmentation. Despite of the extensive studies in experiments and theories, it is still difficult to clarify whether the multifragmentation is statistical or dynamical, sequential or simultaneous. Recently Moretto et al. found that there exists the resilient reducibility and thermal scaling in multiple fragment emission process, which gives one a helpful and clear picture to look and understand the multifragmentation. They observed that the experimental Z-integrated fragment multiplicity distributions $`P_n^m`$ are binomially distributed,

$$P_n^m(p)=\frac{m!}{n!(mn)!}p^n(1p)^{mn}$$
(1)
in each transverse energy ($`E_t`$) window, where $`n`$ is the number of emitted fragments and $`m`$ is interpreted as the number of times the system tries to emit a fragment. The probability of emitting $`n`$ fragments can be reduced to a single-particle emission probability $`p`$ which gives linear Arrhenius plots ($`i.e.`$ excitation functions) when $`ln(1/p)`$ is plotted vs 1/$`\sqrt{(E_t)}`$. By assuming a linear relationship between $`\sqrt{(E_t)}`$ and temperature $`T`$, the linearity of the observed $`ln(1/p)`$ vs 1/$`\sqrt{(E_t)}`$ plot can be explained to a thermal scaling of the multifragment process More ; Tso ; More97 . In this case, these linear Arrhenius plots suggest that $`p`$ has the Boltzman form $`pexp(B/T)`$ with a common fragment barrier $`B`$. However, since the binomial decomposition has been performed on the $`Z`$-integrated multiplicities, typically associated with $`3Z20`$, the Arrhenius plot generated with the resulting one fragment probability $`p`$ is an average over a range of $`Z`$ values. More recently, Beaulieu, Phair and Moretto et al. found that the fit with the binomial distribution can be also replaced with the Poisson distribution in the constraint of charge conservation.

Instead of analyzing for $`Z`$-integrated multiplicities, they analyzed the behavior of individual fragment species of a given $`Z`$ for higher resolution experimental data and noticed that the $`n`$-fragment multiplicities $`P(n)`$ obey a nearly Poisson distribution,

$$P(n)=\frac{<n>^ne^{<n>}}{n!},$$
(2)
where $`n`$ is the number of fragments of a given Z and the average value $`<n>`$ is a function of the total transverse energy $`E_t`$, and were thus reducible to a single-fragment probability proportional to the average value $`<n>`$ for each $`Z`$ Beau . Similarly the $`<n>`$ is found to be proportional to $`exp(B/T)`$ providing that $`T\sqrt{E_t}`$, $`i.e.`$ there exists also a thermal scaling law. As pointed out by Moretto et al. More99 , this kind of ”reducibility” and ”thermal scaling” are empirically pervasive features of nuclear multifragmentation. ”Reducibility” proves nearly stochastic emission process. ”Thermal scaling” gives an indication of thermalization. More recently, Elliott and Moretto et al. discovered that the common features of Poissonian reducibility and thermal scaling can also be revealed in percolation and the Fisher droplet model Elliott . Of course, we should keep in mind that the assumption of $`\sqrt{(E_t)}T`$ may be only valid for compound nuclei formed at low-to-moderate temperatures, but fail at higher temperatures in the experimental data Tsang . In the present work, we will adopt the true temperature to study the reducibility and thermal scaling even though the assumption of $`\sqrt{(E_t)}T`$ is found to be also valid at low-to-moderate temperatures in the present model calculation. By investigating the mean cluster multiplicity as a function of temperature, we will illustrate that the Poissonian reducibility and thermal scaling is also valid for fragment emission in nuclear disassembly via the theoretical reexamination of thermal equilibrium models.

In this Letter, we will analyze the fragment multiplicity distributions for each individual fragment $`Z`$ and $`A`$ value in the framework of isospin dependent lattice gas model (I-LGM) and classical molecular dynamics (I-CMD). We will show that they are Poissonian and the associated mean multiplicities for each $`Z`$ or $`A`$ give linear Arrhenius plots as the experimental data illustrated in Beau . The $`A`$ and $`Z`$ dependent barriers are extracted and investigated as a function of source size. Within our knowledge, this is the first time to explore the Poissonian reducibility and its thermal scaling for the individual fragment in the nuclear disassembly within the lattice gas model and molecular dynamics model.

The lattice gas model was developed to describe the liquid-gas phase transition for atomic system by Lee and Yang Yang52 . The same model has already been applied to nuclear physics for isospin symmetrical systems in the grand canonical ensemble Biro86  with a sampling of the canonical ensemble Jpan95 ; Jpan96 ; Mull97 ; Camp97 ; Jpan98 ; Gulm98 ; Carmona98 ; Ma99 , and also for isospin asymmetrical nuclear matter in the mean field approximation Sray97 . In addition, a classical molecular dynamical model is used to compare with the results of lattice gas model. Here we will make a brief description for the models.

In the lattice gas model, $`A`$ (= $`N+Z`$) nucleons with an occupation number $`s_i`$ which is defined $`s_i`$ = 1 (-1) for a proton (neutron) or $`s_i`$ = 0 for a vacancy, are placed on the $`L`$ sites of lattice. Nucleons in the nearest neighboring sites have interaction with an energy $`ϵ_{s_is_j}`$. The hamiltonian is written as

$$E=\underset{i=1}{\overset{A}{}}\frac{P_i^2}{2m}\underset{i<j}{}ϵ_{s_is_j}s_is_j,$$
(3)
where $`P_i`$ is the momentum of the nucleon and $`m`$ is its mass. The interaction constant $`ϵ_{s_is_j}`$ is chosen to be isospin dependent and be fixed to reproduce the binding energy of the nuclei Jpan98 :

$`ϵ_{nn}`$ $`=`$ $`ϵ_{pp}=0.MeV,`$
$`ϵ_{pn}`$ $`=`$ $`5.33MeV.`$ (4)
Three-dimension cubic lattice with $`L`$ sites is used which results in $`\rho _f`$ = $`\frac{A}{L}\rho _0`$ of an assumed freeze-out density of disassembling system, in which $`\rho _0`$ is the normal nuclear density. The disassembly of the system is to be calculated at $`\rho _f`$, beyond which nucleons are too far apart to interact. Nucleons are put into lattice by Monte Carlo Metropolis sampling. Once the nucleons have been placed we also ascribe to each of them a momentum by Monte Carlo samplings of Maxwell-Boltzmann distribution.

Once this is done the I-LGM immediately gives the cluster distribution using the rule that two nucleons are part of the same cluster if

$$P_r^2/2\mu ϵ_{s_is_j}s_is_j<0,$$
(5)
where $`P_r`$ is the relative momentum of two nucleons and $`\mu `$ is their reduced mass. This prescription is evidenced to be similar to the Coniglio-Klein’s prescription Coni80  in condensed matter physics and be valid in I-LGM Camp97 ; Jpan96 ; Jpan95 ; Gulm98 . To calculate clusters using I-CMD we propagate the particles from the initial configuration for a long time under the influence of the chosen force. The form of the force is chosen to be also isospin dependent in order to compare with the results of I-LGM. The potential for unlike nucleons is

$`v_{\mathrm{np}}(r)({\displaystyle \frac{r}{r_0}}<a)`$ $`=`$ $`A\left[B({\displaystyle \frac{r_0}{r}})^p({\displaystyle \frac{r_0}{r}})^q\right]exp({\displaystyle \frac{1}{\frac{r}{r_0}a}}),`$
$`v_{\mathrm{np}}(r)({\displaystyle \frac{r}{r_0}}>a)`$ $`=`$ $`0.`$ (6)
In the above, $`r_0=1.842fm`$ is the distance between the centers of two adjacent cubes. The parameters of the potentials are $`p`$ = 2, $`q`$ = 1, $`a`$ = 1.3, $`B`$ = 0.924, and $`A`$ = 1966 MeV. With these parameters the potential is minimum at $`r_0`$ with the value -5.33 MeV, is zero when the nucleons are more than 1.3$`r_0`$ apart and becomes stronger repulsive when $`r`$ is significantly less than $`r_0`$. The potential for like nucleons is written as

$`v_{\mathrm{pp}}(r)(r<r_0)`$ $`=`$ $`v_{\mathrm{np}}(r)v_{\mathrm{np}}(r_0),`$
$`v_{\mathrm{pp}}(r)(r>r_0)`$ $`=`$ $`0.`$ (7)
This means there is a repulsive core which goes to zero at $`r_0`$ and is zero afterwards. It is consistent with the fact that we do not put two like nucleons in the same cube. The system evolves for a long time from the initial configuration obtained by the lattice gas model under the influence of the above potential. At asymptotic times the clusters are easily recognized. The cluster distribution and the quantities based on it in the two models can now be compared. In the case of proton-proton interactions, the Coulomb interaction can also be added separately and compared with the cases where the Coulomb effects are ignored.

In this Letter we choose the medium size nuclei 129Xe as a main example to analyze the behavior of individual fragment emission during nuclear disassembly with the helps of I-LGM and I-CMD. In addition, the systems with $`A_{sys}`$ = 80 ($`Z_{sys}`$ = 33) and 274 ($`Z_{sys}`$ = 114) are also studied to investigate the possible source size dependence. In most case, $`\rho _f`$ is chosen to be about 0.38 $`\rho _0`$, since the experimental data can be best fitted by $`\rho _f`$ between 0.3$`\rho _0`$ and 0.4$`\rho _0`$ in the previous LGM calculations Jpan95 ; Beau96 , which corresponds to $`7^3`$ cubic lattice is used for Xe, $`6^3`$ for $`A_{sys}`$ = 80 and $`9^3`$ for $`A_{sys}`$ = 274 . In the condition of the fixed freeze-out density, the only input parameter of the models is the temperature $`T`$. In the I-LGM case, $`\rho _f`$ can be thought as the freeze-out density but in the I-CMD case $`\rho _f`$ is, strictly speaking, not a freeze-out density but merely defines the starting point for time evolution. However since classical evolution of a many particle system is entirely deterministic, the initialization does have in it all the information of the asymptotic cluster distribution, we will continue to call $`\rho _f`$ as the freeze-out density. 1000 events are simulated for each $`T`$ which ensures enough statistics.

One of the basic characters of the Poisson distribution Eq.(2) is the ratio $`\sigma _{n_i}^2/<n_i>1`$ where $`\sigma _{n_i}^2`$ is the variance of the distribution and $`<n_i>`$ is the mean multiplicity. The first step we will show is this ratio. We give these ratios for clusters classified with different masses ($`A`$), light isotopes ($`ISO`$) and atomic numbers ($`Z`$) for the disassembly of $`{}_{}{}^{129}Xe`$ as a function of temperature in the framework of I-LGM and I-CMD with Coulomb in Figure 1. Obviously, the ratios are close to one except for protons, which indicates that they might belong to the Poisson distributions. This can be further supported by the multiplicity distribution in different temperature. For instance, Fig.2 show the quality of the Poisson fits to the charged particle multiplicity distribution for $`{}_{}{}^{129}Xe`$ in the I-LGM case. These Poisson fits are excellent for all $`Z2`$ over the entire range of $`T`$. The same good Poisson fit is obtained in the cases of I-CMD. Thus we can conclude that Poissonian reducibility is valid in the thermal-equilibrium lattice gas model or molecular dynamics.

To verify thermal scaling in the models, the temperature dependence of the mean yield of clusters is investigated. Consequently, we generate Arrhenius plots by plotting $`ln<n_Z>`$ vs $`1/T`$. Figure 3 give a family of these plots for the disassembly of $`{}_{}{}^{129}Xe`$ within the framework of I-LGM (left panel) and I-CMD with Coulomb interaction (right panel). The observed Arrhenius plots are striking linear for the lower $`T`$ side, and their slopes generally increase with increasing $`Z`$ value. Generally, the thermal scaling is expected when the yields, for a fixed nucleon number system, are dominated by fragment binding. This is the case when the temperatures are low compared to the binding energy per particle. At these temperatures, one can anticipate one large fragment surrounded by many small clusters. The contrary tendency reveals in the high $`T`$ side where $`ln<n_Z>`$ increases with $`1/T`$, $`i.e.`$ decreases with increasing $`T`$. In this case, nuclear Arrhenius plots of $`<n>`$ with $`1/T`$ are not valid but the Poissonian reducibility still remains (see Fig. 2). This behavior of $`<n>`$ at higher $`T`$ is related to the branch of the fall of the multiplicity of $`IMF`$ ($`N_{IMF}`$) with $`T`$ where the disassembling system is in vaporization Ogil ; Tsang93 ; Ma95  and hence only the lightest clusters are dominated and the heavier clusters become fewer and fewer with increasing $`T`$. Afterwards we will focus on the branch of lower temperature to discuss the Arrhenius law. The overall linear trend illustrates that thermal scaling is also present when the individual fragments of a specific $`A`$, $`ISO`$ and $`Z`$ are considered.

From figure 3 the slope parameter can be directly extracted in the lower $`T`$ side as a function of $`Z`$ or $`A`$. In Ref. Beau  Moretto et al. has interpreted these slope parameters as ”emission barriers” of specific individual fragments. Figure 4 gives the emission barrier of individual fragments with different $`A`$, $`ISO`$ and $`Z`$ in the framework of I-LGM, I-CMD with/without Coulomb interaction. The error bar in the figure represents the error in the extraction of the slope parameter. The first indication from this figure is that the emission barrier in the I-LGM case is the nearly same as the I-CMD case without Coulomb force, which supports that I-LGM is equivalent to I-CMD without Coulomb interaction rather well when the nuclear potential parameter is moderately chosen, but I-LGM is a quick model to analyze the behavior of nuclear dissociations. The inclusion of long-range Coulomb interaction makes the emission barrier of individual fragments much lower since the repulsion of Coulomb force reduces the attractive role of potential and hence make clusters escape easily. The second indication is that the emission barriers increase with $`A`$ ($`Z`$) at low $`A`$ ($`Z`$) values and tend to be saturated at high $`A`$ ($`Z`$) ones. Similar experimental results have been observed for individual fragments with different $`Z`$ in Ref.Beau  or different $`A`$ in Ref.Elliott  . However, the middle panel of Fig. 4 shows that bare dependence of emission barrier of $`ISO`$ on $`A`$ in the fixed atomic number Z, which indicates that the Z dependence of barrier is perhaps more intrinsical the A dependence is mostly due to the average effect over the species with the same A but different Z.

On the origin of these barrier, the surface energy and Coulomb energy would play the roles. If the cluster emission is mainly controlled by its surface energy, it would suggest barriers proportional to $`Z^{2/3}(A^{2/3})`$. In the case of I-LGM and I-CMD without Coulomb, we can try to fit the barrier for the particles with different mass number by

$$B_{Coul.off}=c_1\times A_{i}^{}{}_{}{}^{2/3},$$
(8)
or for the particles with different charge number by

$$B_{Coul.off}=c_1\times ((A/Z)_{fit}Z_i)^{2/3},$$
(9)
where $`(A/Z)_{fit}`$ is a fit coefficient of A/Z for emitted particles, and $`A_i`$ ($`Z_i`$) is the mass (charge) of particle. $`c_1`$ is the fit constant for surface energy term. The solid line in the Fig.4a is a function of Eq.(8) with $`c_1`$ = 8.469 and the solid line in the Fig.4c is a function of Eq.(9) with $`c_1`$ = 8.469 and $`(A/Z)_{fit}`$ = 1.866. These excellent fits imply that the surface energy play a major role in controlling the cluster emission when the long range Coulomb force is not considered. However for the cluster emission with the Coulomb field, we can assumed that the barrier is mainly constituted by the surface energy term and an additional Coulomb term as

$$B_{Coul.on}=c_2\times A_i^{2/3}\frac{1.44\times A_i/(A/Z)_{fit}\times Z_{res}}{r_{Coul}(A_{i}^{}{}_{}{}^{1/3}+((A/Z)_{fit}Z_{res})^{1/3})}$$
(10)
for the particles classified with different mass number, or

$$B_{Coul.on}=c_2\times ((A/Z)_{fit}Z_i)^{2/3}\frac{1.44\times Z_i\times Z_{res}}{r_{Coul}((Z_i(A/Z)_{fit})^{1/3}+(Z_{res}(A/Z)_{fit})^{1/3})}$$
(11)
for the particles classified with different charge number, where $`c_2`$ is a fit constant for surface term and $`r_{Coul}`$ is chosen to be 1.22 fm. $`Z_{res}`$ is a fitted average charge number of the residue. $`(A/Z)_{fit}`$ is chosen to be 1.866, as taken from the fits for I-LGM. The overall fits for $`A`$ and $`Z`$ dependent barrier in the case of I-CMD with Coulomb force give $`c_2`$ = 12.921 and $`Z_{res}`$ $``$ 41 with the dot-dashed line in Fig.4a and 4c. The excellent fit supports that the Coulomb energy plays another important role in the cluster emission.

In the case of I-LGM and I-CMD without Coulomb, one would expect the barrier for each $`Z`$ ($`A`$) to be nearly independent of the system studied if only the surface energy is substantial to the emission barrier. The left panel of the figure 5 shows the results for $`B_A`$, $`B_{ISO}`$ and $`B_Z`$ for three different systems in the I-LGM case. The same freeze-out density of 0.38$`\rho _0`$ and the same $`N/Z`$ is chosen for the systems of $`A_{sys}`$ = 80 and $`A_{sys}`$ = 274. Actually, it appears to have no obvious dependence of emission barrier on source size as expected for the role of surface energy. The solid line in the figure is the same as in Fig.4. However, when the long-range Coulomb interaction is considered, the emission barrier reveals a source size dependence. The right panel of figure 5 gives the emission barriers $`B_A`$, $`B_{ISO}`$ and $`B_Z`$ in the case of I-CMD with Coulomb force. It looks that the barrier increase with the decreasing of charge of system, which can be explained with the Eq. (10) and (11) where the decreasing of the residue $`Z_{res}`$ will result in the decreasing of the Coulomb barrier and hence the increasing of the emission barrier. The lines represent the fits with the Eq.(10) and (11) for three different mass systems.

In the above calculations, the freeze-out density of systems is fixed at $``$ 0.38$`\rho _0`$. Considering the freeze-out density is an important debating variable in the latter stage of heavy ion collisions, here we will discuss the possible influence of freeze-out density on the emission barrier of clusters. The calculations at the freeze-out density of 0.177$`\rho _0`$ and 0.597$`\rho _0`$ for $`{}_{}{}^{129}Xe`$, corresponding to $`9^3`$ and $`6^3`$ cubic lattices respectively, are supplemented to compare. Figure 6 gives the results of $`B_A`$, $`B_{ISO}`$ and $`B_Z`$ at different density. It looks that there are no obvious freeze-out density dependence in the both cases of I-LGM and I-CMD. This is also consistent with that assumption that the surface energy is the dominant role in controlling the cluster emission.

In conclusion, the poisson reducibility and thermal scaling of the emitted clusters is explored in the lattice gas model and molecular dynamical model. The calculations are qualitatively consistent with the recent experimental observation by Moretto/Wozniak’s group even though the temperature is supposed to be proportional to the total transverse energy in the latter experiments. A systematic study of the emission barrier on the cluster mass, isotope and charge proves that the cluster emission is mainly controlled by both the surface energy and the Coulomb interaction. In the framework of the lattice gas model and molecular dynamics model without the Coulomb interaction, the emission barrier relies on the cluster charge with the $`Z^{2/3}`$ ($`A^{2/3}`$) law and it does not depend on the the source size and freeze-out density, which indicates that the surface energy play a basic dominant role to control the cluster emission. Conversely, in the framework of molecular dynamics model with the Coulomb force, the emission barrier will decrease strongly according to the Eq.(10) and (11) and it decreases with the increasing of the source size, illustrating that the Coulomb interaction also play another weighty role to control the cluster emission.

###### Acknowledgements.

Author would like to thank Prof. B. Tamain, Prof. S. Das Gupta, Prof W.Q. Shen and Dr. J.C. Pan for helps. This work was supported in part by the NSFC for Distinguished Young Scholar under Grant No. 19725521, the NSFC under Grant No. 19705012, the Special Foundation of the President of Chinese Academy of Sciences, the Major State Basic Research Development Program of China under Contract No. G200077400. It was also supported partly by the IN2P3-CNRS Foundation of France.


text:
# A Diffusion Model for Classical Chaotic Compound Scattering

## 1 Introduction

The Multichannel Quantum Defect Theory (MQDT) is a well established quantum approximation framework to accurate spectroscopic studies of atoms  and molecules  near (i.e. both below and above) an ionization limit. More recently the MQDT theory for molecules has been used as an interesting tool to test ideas about the quantum signatures of classical chaos . Since this correspondence is valid only near the semi classical limit, i.e. for large quantum numbers, some of these these quantum numbers have been taken much larger than their usual experimental values. To simplify the mathematical analysis, an approximation (the conservation of the angular momentum $`L`$ of the outer (Rydberg) electron) valid in the case of usual small molecules only for low $`L`$ has been kept for much higher values of $`L`$. Indeed we apply a low-$`L`$ – approximation in the semi-classical limit, i.e. for large $`L`$. Other approximations better suited to actual small molecules have been used in recent calculations . A possible experimental way out is to study much larger molecules, and this is not out of reach since recent experimental studies have found Rydberg states on molecules as large as Benzene Argon complexes, Diazacyclooctane (DABCO) and bis (benzene) chromium (BBC). These experimental studies however involve other issues, namely the vibrational transfers between the very dense set of vibrational levels of such large molecules. While all these experimental issues are very interesting to study, we will not consider them in this paper.

We will use the model only for its interesting theoretical properties with respect to the problem at hand, namely the study of chaotic scattering. Indeed one of us (M.L.) with co workers proposed a classical map, the quantization of which is exactly solved by the MQDT . This approximate system deserves to be studied for its own sake precisely because we can compute the corresponding quantum results without approximation. The study of the correspondence between classical and quantum results, however, will be deferred to future work. In the present paper we establish only interesting classical results. With this proviso we shall in this paper refer to the approximate classical system as the Rydberg molecule.

The Rydberg molecule provides a remarkably simple example of a two dimensional system that displays compound behavior and positive energy bound states typical of many-body systems. It will therefore be of interest to study its scattering functions.

Among the classical scattering functions the time delays are commonly used as they display the singularity structure of such functions very clearly and simultaneously have an intuitive meaning. The number of iterations of an internal Poincaré map are a convenient alternative to time delays that does not reflect distortions by far flung trajectories that are particularly important in the Coulomb case. As we are precisely in this situation we shall mainly use the latter, though we shall also display the mentioned distortions. The results for the number of iterations display a very surprising structure of successive decay laws. We cannot hope to obtain an understanding by reconstructing the chaotic saddle. Indeed we have not been able to resolve the fractal structure of the scattering functions and therefore no branching tree was obtained. The reason was precisely the complexity induced by the far flung trajectories, which produce much too narrow structures (see the discussion at the beginning of sec. 3). We therefore develop a diffusive model that reflects this situation and can be qualitatively related to the dynamical behavior of our system.

## 2 The Classical Model

In ref.  the Rydberg molecule is described classically as a map of the unit sphere onto itself. This sphere corresponds to the orientations of the orbital angular momentum $`\stackrel{}{L}`$ of the electron (whose absolute value is approximated to be constant) after each collision with the molecular core. The construction of the map and its description turn out to be somewhat involved. Conceptually, the simplest approach seems to be that of looking at the system in the laboratory frame. The final results will be stated in a form that is essentially frame independent.

We define a unit vector $`\stackrel{}{P}`$ in the direction of the core axis, the electron angular momentum $`\stackrel{}{L}`$, the core angular momentum $`\stackrel{}{N}`$ and the electronic energy $`E_e`$. This yields a total of seven quantities. (We do not consider the quantities $`\stackrel{}{P}`$, since they are determined up to an initial condition by the time evolution of $`\stackrel{}{N}`$ which is perpendicular to $`\stackrel{}{P}`$.) Between these seven quantities, there exist five relations given by total energy ($`E`$) conservation, the conservation of total angular momentum

$$\stackrel{}{J}=\stackrel{}{L}+\stackrel{}{N}$$
(1)
and finally, the approximate conservation of the magnitude of the electronic angular momentum $`L`$. From this follows that one only needs two numbers to specify fully the configuration of the system once total energy, total angular momentum and $`L`$ are known. For these we choose the two polar angles $`\theta `$ and $`\varphi `$ determining the relative orientation of $`\stackrel{}{L}`$ with respect to $`\stackrel{}{P}`$ as polar axis and $`\stackrel{}{N}`$ as reference axis for $`\varphi `$ or, alternatively, the angles $`\alpha `$ and $`\beta `$ determining the orientation of $`\stackrel{}{L}`$ with respect to $`\stackrel{}{N}`$ and $`Y`$ (see fig. 1). As a matter of notation, we define $`R_\stackrel{}{a}(\psi )`$ to be a rotation of angle $`\psi `$ around the axis $`\stackrel{}{a}`$ in the positive direction. We are now ready to define the two basic components of the map: All quantities with an index $`n`$ refer to a given time step, whereas primed quantities are intermediate quantities for computation only. We therefore aim to express such quantities as $`\stackrel{}{L}_{n+1}`$ in terms of $`\stackrel{}{L}_n`$.

* Free Keplerian motion: When the electron is far from the core, it feels only the isotropic $`1/r`$ potential of the core and thus follows a free Coulomb trajectory, so that its angular momentum $`\stackrel{}{L}`$ is fixed. The core on the other hand, feels no torque from the electron and thus performs a free rotation. This implies for $`\stackrel{}{L}`$, $`\stackrel{}{N}`$ and $`\stackrel{}{P}`$ the following evolution:

$`\stackrel{}{L}^{}`$ $`=`$ $`\stackrel{}{L}_n`$
$`\stackrel{}{N}^{}`$ $`=`$ $`\stackrel{}{N}_n`$
$`\stackrel{}{P}^{}`$ $`=`$ $`R_{\stackrel{}{N}_n}\left(2\pi {\displaystyle \frac{T_e}{T_N}}\right)\stackrel{}{P}_n`$ (2)
$`T_e`$ $`=`$ $`2\pi \left(2E_e\right)^{3/2}`$
$`T_N`$ $`=`$ $`2\pi \left(2BN_n\right)^1`$
Here $`T_e`$ and $`T_N`$ represent the Keplerian period of the electron and the free rotation period of the core respectively, and $`2B`$ is the reciprocal of the moment of inertia of the core, according to usual molecular spectroscopy notation. Atomic units ($`e=\mathrm{}=m=1`$) are used throughout and the equations given simply follow from the elementary mechanical properties of Kepler motion and of free rotation.
* After the free motion, one has a “collision” of the electron with the core, when it is near enough to feel the anisotropic part of its potential. Due to cylindrical symmetry of this potential, the projection $`\mathrm{\Lambda }=\stackrel{}{L}\stackrel{}{P}^{}`$ of $`\stackrel{}{L}`$ on the molecular axis $`\stackrel{}{P}^{}`$ is conserved. Again, as pointed out above, we assume that the magnitude of $`\stackrel{}{L}`$ is conserved. This collision is thus described by a precession of $`\stackrel{}{L}`$ around the molecular axis $`\stackrel{}{P}^{}`$ by an angle $`\delta \varphi `$. Due to invariance of the core potential in any plane containing the molecular axis $`\stackrel{}{P}^{}`$, $`\delta \varphi `$ must be odd in $`\theta \pi \theta `$. We take in the present paper the simplest choice $`\delta \varphi =K\mathrm{cos}\theta =K\stackrel{}{L}\stackrel{}{P}^{}/L`$, which defines the coupling strength $`K`$. In equations, this leads to

$`\stackrel{}{P}_{n+1}`$ $`=`$ $`\stackrel{}{P}^{}`$
$`\stackrel{}{L}_{n+1}`$ $`=`$ $`R_\stackrel{}{P}^{}(\delta \varphi )\stackrel{}{L}^{}`$ (3)
$`\stackrel{}{N}_{n+1}`$ $`=`$ $`\stackrel{}{J}\stackrel{}{L}_{n+1}`$
It should be noted that the effects of the change in $`\stackrel{}{N}`$ in the last equation are rather profound. First, and most important, there is an energetic consequence. This change involves a change in the magnitude $`N`$ of $`\stackrel{}{N}`$ and hence in the core kinetic energy $`E_N=BN^2`$. This in turn reflects itself in a change in $`E_e=EE_N`$, which in turn, will affect both $`T_N`$ and $`T_e`$ in the following free motion. The system is therefore quite far from being a simple kicked top, in which the motion of $`\stackrel{}{L}`$ would be decoupled from the rest of the dynamics. Second there is a kinematical consequence: since $`\stackrel{}{L}`$ has rotated by $`\delta \varphi `$ around $`\stackrel{}{P}^{}`$, and thus moved with respect to $`\stackrel{}{J}`$, this entails a reorientation of $`\stackrel{}{N}`$ and, since the angle $`\varphi `$ is defined with $`\stackrel{}{N}`$ as reference axis (see fig. 1), this produces an extra $`\delta \varphi ^{}`$, which was called frame recoil in ref.

In the Appendix B of this ref.  are given explicit expressions for the map established in molecular frame, which, even if less obvious to visualize, are easier to write down explicitly.

In any case what we shall be always presenting are the values of $`(\theta _n,\varphi _n)`$ describing the orientation of $`\stackrel{}{L}`$ with respect to $`\stackrel{}{P},\stackrel{}{N}`$, which are equally well defined in the laboratory frame or in a fixed molecular frame.

The most important consequence of this exchange of energy between electron and core is that if the total energy $`E`$ is just slightly above the core energy for the lowest core angular momenta allowed, we will find points on this sphere, for which the electron energy $`E_e`$ is positive and others for which it is negative. The positive electron energy region is thus given by $`E>E_N=BN^2`$, which using $`\stackrel{}{J}^2=\left(\stackrel{}{L}+\stackrel{}{N}\right)^2`$ gives

$$BN^2=B\left(L\mathrm{cos}\alpha +\sqrt{J^2L^2\mathrm{sin}^2\alpha }\right)^2<E$$
(4)
This condition depends only on the angle $`\alpha `$ between $`\stackrel{}{L}`$ and the $`\stackrel{}{N}`$ axis, so that the positive electron energy region is a cap about the positive $`\stackrel{}{N}`$\- axis.

Therefore we can describe a scattering process by initiating the map with a point in the positive electron energy region and terminating it whenever we hit this region again. Indeed the restriction of the MQDT map to positive electronic energy has been identified  with the Jung scattering map . Ionization processes would correspond to starting at some internal point and stopping when the positive energy region is hit for the first time.

## 3 Numerical Results

For the part of the sphere corresponding to negative electron energy, one observes the usual KAM scenario as the coupling constant $`K`$ increases. The size of the ordered regions decreases as $`K`$ grows. To get an idea of the general KAM structure in Fig. 2 we show views of the map at total energies $`E`$ just above the ionization energy for two different coupling constants. Note that, for strong coupling, elliptic regions are (nearly) invisible. For the weaker constant these regions are large at the back of the sphere and display the typical features of a twist map. These elliptic regions correspond to positive (total) energy bound states,which are stable in the classical regime. In quantum mechanics we will expect them to decay slowly, and it will be interesting, in some future work, to study their decay law, which will not obviously be exponential, due to the fractal structure of islands surrounding the main elliptic islands.

We next proceed to analyze the time delays in scattering both for strong and weak coupling. Here we have to address the question whether we wish to look at the real time or at the number of iterations of the map. The real time is certainly more attractive at first sight, but the excursions of the electron on a Rydberg orbit take times that approach infinity as the electron energy approaches zero from below. Not only do these long times distort the entire picture but they are actually also unrealistic because any electron will feel some perturbation on a sufficiently long excursion. We shall therefore first consider the delay in terms of iterations of the map in detail and at the very end show the true time evolution.

The most common idea nowadays to study these scattering delays would be to launch a set of trajectories on the circle corresponding to the scattering threshold, and to display the fractal structure of these delays as a function of the angle along this circle. As said in the introduction we saw no structure this way. To understand what is going on we use the following trick. Suppose that for some reason (e.g. external electric field) the electron is ionized for some negative electron energy $`E_e`$: this would correspond to suppress far flung trajectories. We can then launch a set of trajectories on this new threshold circle and study the new delays. We have done this for the same parameters as for all other cases in this paper (see the caption of table 1), which correspond to a very narrow cap of angle $`\alpha _0=0.075\pi `$ of positive electron energies, but selecting artificial ionizing thresholds electron energies which correspond to half sphere ($`\alpha _0=0.5\pi `$) and quarter sphere ($`\alpha _0=0.25\pi `$). The corresponding structure (delay versus angle $`\beta `$ on the ring (see fig. 1)) is displayed on fig. 3. The half sphere case displays clearly the expected fractal delay structure (we have checked that it is indeed fractal by expanding some parts of this structure). For the narrower quarter sphere case the fractal structure becomes very narrow, and for the true zero energy ionization threshold it is no more visible.

We will thus now study the diffusion process, for which it will turn out as we will see that these far flung trajectories instead of being a drawback are an advantage because they are responsible of the interesting structure we will see.

Let us first consider the strong coupling case $`K=1`$: Fig. 4 shows the time evolution of a set of 20000 incoming particles arranged on a ring concentric to the circle that defines the scattering threshold. The sequence shows the original distribution and its evolution after 1, 3, 10 and 100 iterations. Note that the total number of particles diminishes due to scattered particles. Thus the density diminishes but after a small number of iterations this density tends to be very uniform, i.e. the system has equilibrated. Note that the escaping particles are drawn as empty squares that appear on the cap of the sphere that corresponds to positive electron energies.

In Fig. 5 top, we show the number of particles escaping at a given iteration in a semi-logarithmic plot. We see the straight line characteristic of exponential decay dominating the figure, and an anomaly in the first channel. To get a better feeling for this short time behavior, we next look at a log-log plot of the same data in Fig. 5 bottom. Here we see that the short time ($`t10`$) behavior is a straight line i.e. a power-law. Whereas the exponential behavior is standard for chaotic scattering, the power law for the transient is puzzling at first sight.

To gain better understanding we pass to the weak coupling case given by $`K=0.1`$. Again we first look in Fig. 6 at a sequence of pictures after 1, 5, 1000, 10 000, iterations. Here we see that the particles diffuse very slowly towards the back of the sphere and they never penetrate the elliptic islands. At very large numbers of iterations some get quite close, and thus may penetrate the fractal tangle of small islands.

To get a quantitative understanding of the observed behavior, we again plot the number of escaping particles versus iteration in a log-log plot (fig. 7). We see an initial power law as in strong coupling for now roughly 1000 channels. This power law can now be accurately fitted to determine its power $`3/2`$ (cf table 1). This result is compatible with the findings of fig. 5, -bottom

Furthermore in this weak coupling case an other distinct deviation from exponential behavior appears at very long times ($`t\mathrm{2\hspace{0.17em}10}^5`$). This last region is compatible with a power law between -1 and -3. Let us treat first this second power law. This is what we expect for particles caught in the chaotic tangle . That this is indeed the effect we are seeing, can easily be shown as follows. In fig. 8 we see that some particles show regions in their time-evolution plot where the angle $`\alpha `$ between $`\stackrel{}{L}`$ and $`\stackrel{}{N}`$ does not change over many iterations indicating trapping in the fractal tangle around an island which has this behavior according to fig. 2. We may actually subtract these iterations and make a new plot of delays where these are absent as shown in Fig. 9. The final power law tail disappears, but the lower times region is not affected.

## 4 Diffusion Model

Let us focus now on the interpretation of this first power law, which is the most interesting part of our results. The law

$$At^{3/2}\mathrm{exp}(t/T)$$
(5)
corresponds, as is well-known, to the probability that a 1D random walk starting from the origin in an interval of length $`L\sqrt{T}`$ returns to the origin for the first time at time $`t`$. We could readily explain the connection to the Rydberg model by the following argument: Our map consists of two parts: during the long Rydberg excursion of the electron the core rotates around the $`\stackrel{}{N}`$ axis very often and we have to take the total rotation modulo $`2\pi `$ which makes for a reasonable randomizer of the rotation angle $`\beta `$. Yet the kick $`\delta \varphi =K\mathrm{cos}\theta `$ produced in the orthogonal direction (around $`\stackrel{}{P}`$) depends on the final value of this angle $`\beta `$. In particular it is towards the front or back of the sphere depending whether this $`\stackrel{}{N}`$ rotation left it on the upper $`z>0`$ or lower $`z<0`$ part of the sphere. It thus becomes a random walk towards the front or back of the sphere. This could explain our findings. However there is a distinct qualitative difference between this 1D random walk law and our Rydberg results. It is displayed by comparing the previous fig. 9 with the theoretical fig. 10. Indeed, the law (5) always leads to a value that is lower than the corresponding power law $`At^{3/2}`$, whereas in our Rydberg case there is a long, nearly horizontal, shoulder between the power law and the exponential decay. Mathematically the two histograms, (7) and (9) for the Rydberg case can be fitted by a formula

$$\left(A_1t^{\alpha _1}+A_2\right)\mathrm{exp}(t/T)+A_3t^{\alpha _3}$$
(6)
with $`A_3=0`$ for the latter case. Parameters are given in table 1. There are negligible residuals: the sluggish transition between the two parts is an effect of non additivity in a log-log plot, and there is no statistically significant possibility to fit two independent time constants for $`A_1`$ and $`A_2`$. This introduces a second characteristic time,

$$t_c=(A_1/A_2)^{1/\alpha _1},$$
(7)
which is the transition time between power law and exponential regimes, besides the time constant of the exponential decay $`T`$, which is much longer. The time $`t_c`$ has no equivalent in the 1D random walk model. In both strong and weak coupling cases $`t_c`$ corresponds roughly according to figs. 4 and 6 to the time needed for the particles to fill two thirds of the sphere.

To check these ideas we replace the true rotation $`\delta \beta `$ by a random one, equidistributed on $`[0,2\pi )`$, while keeping the exact Rydberg formula for the collisional kick $`\delta \varphi +\delta \varphi ^{}`$ and resulting $`\delta \alpha `$, which is the significant parameter to monitor distance to ionization. The result (fig. 11 and table 1) is virtually indistinguishable from the case where trapped parts of trajectories are eliminated (fig. 9). Since in this case the regular part at the back of the sphere has disappeared (see fig. 12), this result first corroborates our previous conclusion that the long time tail is due to trajectories trapped in the chaotic tangle. Second it shows that the model of random $`\beta `$ is correct but that the right Rydberg value of the kick $`\delta \varphi +\delta \varphi ^{}`$ introduces a qualitative difference with respect to the simple model of a constant amplitude kick. To gain further insight, since the formula for $`\delta \varphi +\delta \varphi ^{}`$ is too complex to enable to visualize easily its effects, we compute a first order development of $`\delta \alpha `$ as a function of the small parameter $`K`$:

$$\delta \alpha =\left(1\frac{L}{\sqrt{J^2L^2\mathrm{sin}^2\alpha }}\mathrm{cos}\alpha \right)^1\frac{1}{2}K\mathrm{sin}\alpha \mathrm{sin}2\beta $$
(8)
Neglecting the recoil $`\delta \varphi ^{}`$, which is valid if $`LJ`$, it simplifies further as:

$$\delta \alpha =\frac{1}{2}K\mathrm{sin}\alpha \mathrm{sin}2\beta $$
(9)
With any of these first order formulas, the resulting histogram (fig. 13 and table 1) corresponds to the 1D random walk prediction of Eq. (5) and fig. 10, with no exponential shoulder. What is characteristic in these equations is that for an equipartition of $`\beta `$ on $`[0,2\pi )`$ the average value of $`\delta \alpha `$ is zero. The exact value of the law, especially the way it varies with $`\alpha `$ is of little importance. This is expected since a variation of $`|\delta \alpha |`$ with $`\alpha `$ can be compensated for by a renormalisation of $`\alpha `$ to keep its variation constant on average. The existence of the exponential shoulder is thus related to the fact that the true $`\delta \alpha `$ is not zero on average. This is shown by writing $`\delta \alpha `$ up to second order in $`K`$ (neglecting $`\delta \varphi ^{}`$ to simplify):

$$\delta \alpha =\frac{1}{2}K\mathrm{sin}\alpha \mathrm{sin}2\beta +\frac{1}{4}K^2\mathrm{sin}2\alpha \mathrm{sin}^4\beta .$$
(10)
The resulting histogram (fig. 14) displays indeed the exponential shoulder. The explanation is as follows. The average value of $`\delta \alpha `$ is positive for the front part of the sphere ($`0<\alpha <\pi /2`$), thus pushing it backwards towards the middle (increasing $`\alpha `$), and conversely negative for the back part of the sphere, pushing it also towards the middle. This creates a kind of trap in the middle of the sphere, and the exponential law results from the need to escape that trap.

To justify this explanation let’s put it quantitatively. We first renormalize $`\alpha `$ to have constant average $`|\delta \alpha |`$ by redefining a new coordinate

$$x=\frac{1}{2}\mathrm{ln}\frac{1+\mathrm{cos}\alpha }{1\mathrm{cos}\alpha }$$
(11)
which varies between $`\mathrm{}`$ and $`+\mathrm{}`$ when $`\alpha `$ varies between $`\pi `$ (back) and $`0`$ (front of the sphere). $`x`$ has been computed so that to first order $`\delta x`$ is independent of $`x`$, and then to second order in $`K`$

$`\delta x`$ $`=`$ $`\frac{1}{2}K\mathrm{sin}2\beta \frac{1}{2}K^2\mathrm{cos}\alpha \left(\mathrm{sin}^4\beta \frac{1}{4}\mathrm{sin}^22\beta \right)`$ (12)
$`=`$ $`\frac{1}{2}K\mathrm{sin}2\beta \frac{1}{2}K^2\mathrm{tanh}x\left(\mathrm{sin}^4\beta \frac{1}{4}\mathrm{sin}^22\beta \right)`$
We obtain thus a random motion with differing average amplitude in the two directions. This is equivalent to a thermal motion within a potential $`V(x)`$. The average step is $`\pm K/2\pi `$ and the difference between the two directions comes from a Boltzmann factor:

$$\mathrm{exp}\left(V^{}(x)\times K/\pi \right)1V^{}(x)\times K/\pi =\frac{\delta x^+}{\delta x^{}}=1\frac{\pi }{4}K\mathrm{tanh}x,$$
(13)
where we have normalized temperature and potential so that $`kT=1`$. Thus

$$V(x)=\frac{\pi ^2}{4}\mathrm{ln}\mathrm{cosh}x$$
(14)
This potential is plotted in fig. 15. The potential curve is limited on the right at the actual value of $`x`$ corresponding to the ring of ionization in all preceding figures. The depth of the trap is $`3.6`$, justifying a long time exponential to escape it.

In order to compare the results given here with the times reported in table 1, one should note the following: these times are given in number of iterations. However, at each iteration, $`x`$ changes by an amount of order $`K`$, so that the natural diffusive time scale, over which $`x`$ varies by a quantity of order one, is of order $`K^2`$. We therefore expect the times given in table 1 to contain an overall factor of order $`K^2`$ apart from the first passage probability. That this is indeed so, is verified in the case $`K=0.2`$, for which all times are approximately four times shorter than for $`K=0.1`$.

The physical origin of this potential is obvious. The electron can both gain or loose energy at each collision with the core. The most probable situation is equipartition of energy, i.e. the middle of the sphere, and ionization which needs nearly all energy in the electron is a rather improbable situation. This potential takes thus into account the amount of available phase space.

We now give a theory of this phenomenon. We wish to solve the following problem: if we introduce a particle at the point $`a`$ and this particle diffuses under the influence of a binding potential $`V(x)`$, what is the probability that it returns to $`a`$ for the first time at time $`t`$, if the potential is of such a nature that being at the position $`a`$ has quite a small probability in equilibrium? We shall show that in these circumstances it is to be expected on quite general grounds that this probability has the behavior outlined above, namely a $`t^{3/2}`$ initial decay, followed by a plateau and finally ending in an exponential decay. At the intuitive level, this can be seen as follows: Near $`a`$, the process can be viewed as a diffusion with a drift away from $`a`$. This means that the distribution of first return times is of the type $`t^{3/2}e^{t/t_c}`$. This distribution, however, is not normalized to one, and the missing probability corresponds to particles that diffuse away to infinity. These particles, however, will not really diffuse to infinity, but rather reach the equilibrium distribution $`\mathrm{exp}(V(x))`$. Under these circumstances, it will from then on reach $`a`$ at a rate $`\mathrm{exp}(V(a))`$. This therefore yields an exponential decay in the first passage probability on a time scale $`T`$ much larger than $`t_c`$, hence yielding the expected plateau. A more formal derivation is found in the Appendix.

## 5 Final Remarks

We finally display the true time delays in the two cases of strong and weak coupling, that is, we compute the actual time of each Rydberg excursion. We note in Fig. 16 that the different regimes detected in the iteration plots become quite blurred. This is due to the fact that two very similar trajectories near threshold have greatly different excursion times as the latter diverge at this threshold. The first consequence is that time curves span only two orders of magnitude whereas iteration curves span eight orders of magnitude. The first channel in the strong case and the first isolated peak in the weak coupling case correspond to electrons which are stabilized in the first collision and ejected exactly in the second. These features depend very critically on the exact position with respect of threshold of the incident ring of particles, whereas the remaining of the picture as well as iteration curves are stable against such minor variations. The $`t^{3/2}`$ behavior in the weak coupling case corresponds only to a longer decreasing plateau, the slope of which is no more 1.5. The delay between $`t^{3/2}`$ and exponential decay is no more visible, And finally the last $`t^n`$ decay has also disappeared. If we want to obtain information about the process, we conclude that the true time delays do not yield this information in a clear cut fashion.

We may conclude that our model displays typical properties of compound behavior in many respects. Yet the remarkably complicated structure of the time delay indicates that more information can be obtained than is commonly assumed e.g. in nuclear compound reactions. The diffusive model seems promising, as it can readily describe situations where the regimes are not as clearly separated, which may be more common.

As our analysis is limited to classical mechanics it remains an interesting open question, how this situation translates to quantum mechanics. Also one will have to ask what features are very model dependent and what features may be found in more general context. In particular we think about the possibility of a diffusive description of pre compound situations.

## Appendix A Appendix

In this Appendix, we shall show at a more formal level that the behavior of the type described in the text is quite general for first return probabilities in a potential that goes to infinity when $`|x|\mathrm{}`$. The master equation satisfied by the probability $`P(x,t)`$ of reaching $`x`$ at time $`t`$ is given by

$$\frac{P(x,t)}{t}=\frac{}{x}\left[\frac{P}{x}+V^{}(x)P\right]LP,$$
(15)
where the last equality defines the operator $`L`$. Here the initial condition is given by

$$P(x,0)=\delta (xa).$$
(16)
From this one obtains

$$P(a,t)=\underset{m=0}{\overset{\mathrm{}}{}}\frac{|\varphi _m(a)|^2}{\varphi _0(a)}e^{E_mt},$$
(17)
where $`E_n`$ and $`\varphi _n(x)`$ are the eigenvalues and eigenfunctions of $`L`$ respectively. The division by $`\varphi _0(a)`$ is related to the fact that $`L`$ is not self-adjoint with respect to the ordinary scalar product, so that the $`\varphi _n(x)`$ are not orthogonal. However, replacing $`\varphi _n(x)`$ by $`\varphi _n(x)/\sqrt{\varphi _0(x)}`$ restores self-adjointness.

As is well-known, the Laplace transform $`Q(s)`$ of the first passage probability $`q(t)`$ is given by

$$Q(s)=\frac{P(s)}{1+P(s)},$$
(18)
where $`P(s)`$ is the Laplace transform of $`P(a,t)`$ given by (17).

Now we need both a large and a small time estimation for the expression (17). The large time estimate is trivial: Since (15) always has an equilibrium solution, $`L`$ always has a zero eigenvalue. At large times, this eigenvalue always dominates, and one has

$$P(a,t)=\varphi _0(a)=e^{V(a)}.$$
(19)
For short times, on the other hand, one may simply assume that the diffusing particle does not go far from $`a`$, so that we replace $`V(x)`$ by its Taylor development to first order. Then (15) can be solved exactly to yield

$$P_d(s)=\frac{1}{\sqrt{4s+V^{}(a)^2}},$$
(20)
where $`P_d(s)`$ is the approximation to $`P(s)`$ in this approximation. From (17) follows, however, that $`P(a,t)`$ is never less than $`\varphi _0(a)`$, which we may incorporate approximately by setting

$$P(s)=P_d(s)+\frac{\varphi _0(a)}{s}.$$
(21)
From this follows

$$Q(s)=\frac{\varphi _0(a)\sqrt{4s+V^{}(a)^2}+s}{(s+\varphi _0(a))\sqrt{4s+V^{}(a)^2}+s}.$$
(22)
If $`sV^{}(a)^2`$, this reduces to a rational function in $`s`$, which corresponds to a double exponential decay, with decay times $`\tau _1`$ and $`\tau _2`$ approximately given by

$$\tau _1=\frac{2}{V^{}(a)\left(1+V^{}(a)\right)}\tau _2=\frac{1+V^{}(a)}{\varphi _0(a)V^{}(a)}$$
(23)
The $`\tau _2`$ decay clearly dominates the picture at those times. Since $`\tau _2`$ is much larger than the time at which this approximation starts to be valid, namely approximately $`V^{}(a)^2`$, this implies the existence of the plateau followed by an exponential decay. On the other hand, in the large $`s`$ approximation we may altogether neglect the correction made in (21) and use $`P_d(s)`$. This yields

$$Q(s)=\frac{1}{\sqrt{4s+V^{}(a)^2}+1},$$
(24)
which yields a $`t^{3/2}`$ decay cut off by an exponential at a time scale of order $`V^{}(a)^2`$. This meshes therefore quite well with the large time asymptotics, which just then begins to set in.

## List of Tables

## List of Figures


text:
# The Quest for the Dominant Stellar Population in the Giant Elliptical NGC 5018

## 1. Introduction

NGC 5018, a relatively near, giant E3 galaxy, deserved quite a large attention in the past. Evidence of past interaction with other systems comes mainly from the extended pattern of dust lanes/patches clearly recognizable from ground-based observations. What is more, this object captures immediately the attention of the observer for (i) its anomalous location on the Mg2-luminosity relation of Schweizer et al. (1990), having one of the weakest absorption-line strengths (Mg2=0.218) relative to its absolute magnitude and central velocity dispersion among giant ellipticals and (ii) the extremely weak far-UV emission Bertola et al. (1993) find in its IUE spectrum (strictly resembling the UV energy distribution of the dwarf elliptical M32). Interesting enough, though one expects ongoing star formation plays a major role in a recent merger such as NGC 5018, both the observed low value of its optical Mg2 line strength index and its weak far-UV flux can fit into scenarios not invoking recent star formation activity. According to this latter view, NGC 5018 has been alternatively described either as a sort of “Celacanthus galaxy”, i.e. a rare, chemically unevolved giant elliptical hosting a dominating M32-like population (Bertola et al. 1993), or a system fully dominated by an intermediate-age (2.8 Gyr) population with a roughly solar composition as put forward by Leonardi & Worthey (2000). On the other side, since other groups (e.g. Carollo & Danziger 1994) are prone to conclude that a large amount of truly young population and strong dust absorption might conspire to cause the above UV/optical properties, strict constraints on the amount of internal extinction are urgently needed in order to properly address the issue of the dominant population age and metallicity.

## 2. Data Analysis and Discussion

The above goal can now be achieved thanks to the recent set of Schweizer’s WFPC2 F555W and F814W deep images (program ID 7468) which allows us to finally get rid of the persistent suspicion of a heavy dust absorption affecting the UV emission of this galaxy. In particular, starting from the high S/N E(V$``$I) map shown in Fig. 1, one derives—after removing the underlying V$``$I color gradient—a luminosity weighted average value E(B$``$V)$``$0.02 due to structured dust within the IUE large aperture (10”$`\times `$20”). This, in turn, assures that the collected UV flux is only marginally affected by the internal extinction due to dust (namely $``$0.06 dex at $`\lambda `$ 1500 Å).

## References

Bertola, F., Burstein, D., and Buson, L.M. 1993, ApJ, 403, 573

Carollo, C.M, and Danziger, I.J. 1994, MNRAS, 270, 743

Leonardi, A.J., and Worthey, G. 2000, ApJ, 534, 650

Schweizer, F., Seitzer, P., Faber, S.M., et al. 1990, ApJ, 364, L33