text:
# XMM-Newton observes flaring in the polar UZ For during a low state

## 1 Introduction

UZ For is an eclipsing member of the subclass of cataclysmic variables called AM Her binaries or polars. In these binaries, the strong magnetic field of the white dwarf primary causes it to rotate synchronously with the orbital motion. The magnetic field also prevents the formation of an accretion disk around the white dwarf. The accretion stream from the Roche-lobe filling secondary is funneled along the magnetic field lines and impacts the white dwarf near a magnetic pole. Slightly above the surface, the accretion stream forms a stand-off shock that heats the gas to temperatures in excess of $`10^8`$K. The shock-heated plasma then cools and settles on to the white dwarf while strongly emitting cyclotron radiation (IR to UV) and thermal bremsstrahlung (mostly X-rays). The photosphere below the shock is heated by reprocessing of X-rays and emits blackbody radiation visible at soft X-ray and UV energies. Many polars have been observed to decline in brightness by several magnitudes and remain in a faint state for days to years. The causes of these low states are not known, but, in the absence of an accretion disk, the large brightness variations must be due to changes in the mass transfer rate from the companion star. A comprehensive review of polars is given in Warner (1995).

UZ For (EXO 033319-2554.2) was first detected as a serendipitous X-ray source with EXOSAT (Giommi et al., 1987; Osborne et al., 1988). Subsequent optical spectroscopy and polarimetry established UZ For as an eclipsing polar (Beuermann et al., 1988; Berriman & Smith, 1988). The 126.5–min orbital period is close to the lower edge of the 2–3 hr ”period gap”, a sparsely populated region in the orbital period distribution of cataclysmic variables. UZ For is a high inclination system ($`i80^{}`$), and both accretion regions are eclipsed by the white dwarf for at least half of an orbital cycle. The optical spectrum of UZ For shows strong cyclotron emission lines which indicate magnetic fields of $`53`$MG and $`48`$MG for the two accretion poles (Rousseau et al., 1996). The best available estimates for the distance and the white dwarf mass are $`d=208\pm 40`$pc and $`M_{WD}0.6`$$`0.8M_{}`$ (Ferrario et al., 1989; Bailey & Cropper, 1991).

In this paper we present X-ray and UV data obtained with XMM-Newton while UZ For was in a state of unusually low and irregular accretion. We study the properties of an X-ray/UV flare and show that it was caused by accretion on to the white dwarf. We argue that the intermittent accretion rate increase was due to stellar activity on the companion star. Part of the XMM-Newton data has recently been published in Still & Mukai (2001). The goal of our paper is to utilize all available XMM-Newton data and present additional results not previously published. In particular, we find two more UV flares and clearly identify eclipses in the UV light curves. We show that the beginning of the X-ray/UV flare coincides with the eclipse egress of the main accretion region. This provides strong evidence that the flare was an accretion event on the white dwarf. We estimate the total accreted mass and show that it is consistent with the mass ejected by stellar flares on M dwarfs. We also detect a weak orbital modulation of the X-ray and UV fluxes and demonstrate that it is most likely due to emission from the main accretion region.

## 2 Observations and analysis

UZ For was observed with XMM-Newton on 14 January 2001 as part of the performance verification of the mission. XMM-Newton (Jansen et al., 2001) carried out two separate observations, but during the first observation the EPIC cameras (Strüder et al., 2001; Turner et al., 2001) were used for calibration purposes, and only the Optical Monitor (Mason et al., 2001) was observing the source. During the second observation each of the three EPIC instruments performed one long exposure (16 ks for PN; 18 ks for MOS) with the thin blocking filters and the CCDs in large window mode. The Optical Monitor (OM) performed a total of 11 exposures of 2200 s duration each (3 during the first and 8 during the second observation). All OM exposures were done in fast mode (0.25 s timing resolution) and with the UVW1 filter (240–340 nm).

From the EPIC MOS and PN data we extracted good photon events (FLAG=0) using a circular aperture centred on UZ For and with a 20” radius. Count rates quoted throughout the paper have been corrected for the 75 per cent enclosed energy fraction of this aperture. The background rate was estimated from a larger region on the same CCD as the source image. We included in our analysis events in the energy range 0.10–12 keV for the MOS and 0.15–12 keV for the PN. To maximize the signal-to-noise ratio, we used photon events with patterns 0–12. Only for the spectral analysis of the PN data did we select exclusively events with pattern 0. To create X-ray light curves, we applied a barycentric correction to the photon arrival times and combined the events from all three EPIC instruments. For the extraction of UV light curves from the OM fast mode data, we used a circular aperture of 5.5” radius, which encloses 78 per cent of the detected source photons. Background subtracted UV and X-ray light curves of the available data are shown in Fig. 1. We defined as orbital phase zero the centre of the first eclipse as predicted by the most recent ephemeris (Perryman et al., 2001), which has an accumulated uncertainty of 4 s.

### 2.1 X-ray light curve

Except during a short flare, the X-ray flux from UZ For (Fig. 1) was at an extremely low level never before seen in this object. Even with the EPIC data integrated over the entire observation (excluding the flare), UZ For was detected only at a $`4\sigma `$ level. An improved $`6\sigma `$ detection could be achieved by limiting the integration to the orbital phase range 0.65–0.15. For these orbital phases, the main accretion region is on the side of the white dwarf facing us (Perryman et al., 2001). For orbital phases 0.15–0.65, for which the main accretion region is not visible, the integrated source count rate is consistent with zero at a $`1\sigma `$ level. It appears that the main accretion region was the source of the weak X-ray emission. We cannot rule out the possibility that the X-rays originated from the secondary region, which is visible for orbital phases 0.75–0.35. Yet one might expect that, as for other polars during low states, accretion was limited to the main region (e.g. Beuermann & Schwope, 1989).

The combined EPIC MOS+PN count rate, averaged over the orbital phase range 0.65–0.15 (excluding the flare), is $`0.009(2)`$s-1 (0.15–12 keV). Ramsay et al. (1993) found, based on ROSAT observations, that the X-ray spectrum of UZ For in the high state is best described by an absorbed 28 eV blackbody below 0.5 keV plus a weak thermal bremsstrahlung tail visible at higher energies. Assuming this high-state spectrum, we would predict a XMM-Newton count rate of $`7`$s-1, which is $`800`$ times higher than observed. The EPIC spectrum does not show any evidence of the blackbody component that was seen with ROSAT. However, the spectrum is in good agreement with a thermal bremsstrahlung model. It is possible that, due to a lower temperature of the accretion spot, the blackbody component was shifted below 0.1 keV and not detectable by the EPIC instruments. This would require a blackbody temperature $`15`$eV, provided that the emitting area was the same as during the ROSAT observation. If we consider only the weak bremsstrahlung component found in the ROSAT spectrum, we would predict a XMM-Newton count rate of $`0.26`$s-1. This is still $`30`$ times higher than observed.

We estimated the total bremsstrahlung flux (excluding the flare) by fitting a single temperature thermal bremsstrahlung model to the EPIC spectra. With the temperature fixed at 10 keV, we obtained a bolometric bremsstrahlung flux $`F_{brems}2\times 10^{14}\mathrm{erg}\mathrm{cm}^2\mathrm{s}^1`$. Of this flux, 80 per cent was inside the 0.10–12 keV energy range and directly detected by XMM-Newton. For a distance $`d=200`$pc, the corresponding bolometric luminosity is $`L_{brems}7\times 10^{28}\mathrm{erg}\mathrm{s}^1`$. Here we used a geometric factor of $`3.1\pi `$, which takes into account the scattering of X-rays off the white dwarf’s surface (King & Watson, 1987). Assuming that all of the accretion energy was emitted in the bremsstrahlung component, we can estimate the accretion rate $`\dot{M}`$ using $`L_{brems}GM_{WD}\dot{M}/R_{WD}`$. For a white dwarf with a mass $`M_{WD}=0.7M_{}`$ and a radius $`R_{WD}=8\times 10^8`$cm, we find $`\dot{M}6\times 10^{11}\mathrm{g}\mathrm{s}^1`$. This is likely an underestimate, since at such a low accretion rate most of the energy would be emitted as cyclotron radiation. For the parameters of UZ For, the calculations by Woelk & Beuermann (1996) (fig. 9) predict that cyclotron cooling dominates at specific accretion rates $`\dot{m}<10\mathrm{g}\mathrm{cm}^2\mathrm{s}^1`$. Therefore, if most of the energy were radiated as bremsstrahlung, accretion would have to occur in an unreasonably small region with $`f10^8`$ ($`f`$ is the area as a fraction of the white dwarf’s surface area). For a more typical spot size of $`f=5\times 10^4`$ (Rousseau et al., 1996), we estimate from fig. 9 in Woelk & Beuermann (1996) that the cyclotron luminosity is $`10^2\times L_{brems}`$. The corresponding accretion rate of $`10^{13}`$$`10^{14}\mathrm{g}\mathrm{s}^1`$ is comparable to that determined by Rousseau et al. (1996) for the low state of UZ For.

As shown in Fig. 1, a short flare was detected simultaneously by the X-ray instruments and the Optical Monitor. The flare (Fig. 2), which lasted $`900`$s, increased the X-ray flux by a factor of $`30`$. This is comparable to the flux that we estimated for the bremsstrahlung component during the high state (see above). The beginning of the flare at BJD(TT) 2451924.24401(2) coincides to within a few seconds with the eclipse egress of UZ For (see also Section 2.2). This is strong evidence that the flare emission originates from the white dwarf in UZ For. We will discuss the spectral properties of the flare in Section 2.3.

### 2.2 UV light curve

During most of the observation, the UV flux detected with the Optical Monitor was at a roughly constant level (Fig. 1). Besides the flare that was also detected with the EPIC instruments, two more flaring events are seen near orbital phases 0 and 4. Unfortunately, no X-ray data is available for these time periods. It should be noted that the flaring activity is only seen inside the orbital phase range 0.65–0.15 for which the main accretion region is visible. This is a further indication that the flaring is caused by accretion on to the main pole.

We were able to identify the ingress and egress for several eclipses of the white dwarf by the companion star (see Fig. 1 near integer orbital phases). Coincidentally, each of the eclipses overlapped, at least partially, with a gap in the OM data. We did not detect any UV emission during the eclipses, even after integrating all 530 s of available eclipse data. This non-detection corresponds to an upper count rate limit of $`0.05`$s-1 ($`3\sigma `$) or a flux limit of $`2.5\times 10^{17}\mathrm{erg}\mathrm{cm}^2\mathrm{s}^1`$Å-1 at 270 nm. Besides the accretion stream, the only potential source of emission during eclipse is the M-dwarf companion. Consistent with our limit, the models of Houdebine et al. (1996) predict a flux of $`10^{18}`$$`10^{16}\mathrm{erg}\mathrm{cm}^2\mathrm{s}^1`$Å-1 for the companion in UZ For. Because of the fast rise and decline during eclipse, the white dwarf is probably the sole source of the UV emission. Due to the low accretion rate, a significant contribution from the accretion stream is unlikely.

The roughly constant UV count rate of $`0.75`$s-1 that is seen during most of the observation is probably thermal emission from the white dwarf’s photosphere. The corresponding flux of $`3.8\times 10^{16}\mathrm{erg}\mathrm{cm}^2\mathrm{s}^1`$Å-1 at 270 nm is consistent with blackbody radiation from a 11000–K white dwarf ($`R_{WD}=8\times 10^8`$cm, $`d=200`$pc). The temperature is in agreement with that measured by Bailey & Cropper (1991). We find that the UV light curve (excluding the flares) has a sinusoidal orbital modulation of $`30`$ per cent (Fig. 3). The peak of the modulation is seen at orbital phase 0.88(2). Based on the spot longitudes determined by Perryman et al. (2001), we predict that emission from the main accretion region peaks at phase $`0.88`$, while that from the secondary region peaks at $`0.03`$. This strongly suggests that the modulated flux is thermal radiation from the heated photosphere near the main accretion spot. Another possible origin of the modulation is cyclotron radiation from the accretion column. Then, however, one would expect to see a double-peaked light curve with peaks near orbital phases 0.15 and 0.65. Although cyclotron radiation probably dominates the energy output from the accretion column, only a small fraction of it ($`<1`$ per cent) is emitted in the UV. Assuming that the modulation is due to blackbody radiation from the main accretion region, we find that the flux excess of $`1.2\times 10^{16}\mathrm{erg}\mathrm{cm}^2\mathrm{s}^1`$Å-1 at 270 nm is equivalent to a luminosity of $`4.2\times 10^{26}\mathrm{erg}\mathrm{s}^1`$Å-1 ($`d=200`$pc). Because the accretion region is optically thick and inclined with respect to the line of sight, we applied a geometric correction of $`\pi \mathrm{cos}^1(\delta i)9.2`$ (inclination $`i80^{}`$, colatitude $`\delta 150^{}`$; from Perryman et al., 2001). As discussed in Section 2.1, we could not identify a blackbody component in the X-ray spectrum. Using this non-detection, we can place an upper limit on the blackbody temperature. Depending on the hydrogen column density, which is in the range $`0.25`$$`1.3\times 10^{20}`$cm-2 (Osborne et al., 1988), we find a limit of 7–11 eV (80,000–130,000 K). We can also constrain the bolometric blackbody luminosity $`L_{bb}1\times 10^{32}\mathrm{erg}\mathrm{s}^1`$ and the size of the emitting region $`f0.001`$. The limit on $`L_{bb}`$ is fairly high, $`10^3`$ times higher than the bolometric bremsstahlung luminosity $`L_{brems}`$. It is therefore likely that the actual temperature of the emitting region is significantly lower than the upper limit. The lowest possible blackbody luminosity consistent with the observed UV flux is $`1.5\times 10^{30}\mathrm{erg}\mathrm{s}^1`$ (for $`kT_{bb}=1.1`$eV). This is still $`20`$ times higher than $`L_{brems}`$. A large ratio of blackbody to bremsstrahlung luminosity is found in many polars, typically at high accretion rates. This ”soft excess” is commonly explained by blobs in the accretion stream that penetrate deep into the photosphere (Frank et al., 1988). Blackbody radiation from this process is only emitted in the fairly small accretion region. In UZ For, however, the UV flux is varying sinusoidally even when the main accretion spot is not visible (phases 0.15–0.65). This is only possible if the emitting region is significantly larger than the accretion spot and comparable in size to the white dwarf. Apparently, we are seeing a large heated region of the white dwarf’s photosphere around the main accretion spot. For a covering fraction $`f=0.5`$, this region has to be, on average, $`1000`$K hotter than the other half of the white dwarf’s surface. Since $`L_{bb}`$ is much larger than $`L_{brems}`$, reprocessing of X-rays is unlikely to be responsible for the heating of this region. However, the cyclotron luminosity is expected to be much larger than $`L_{brems}`$ and might provide sufficient energy (see Section 2.1). It is also possible that we are seeing the afterglow from an earlier state of high accretion rate, since the time-scale for cooling of a heated pole cap is on the order of tens of days (Gänsicke, 1997).

We determined the times of mid-ingress and mid-egress by fitting the profile expected for a circular emission region to the eclipses in the UV light curve (Fig. 1). Except for the egress near phase 2, we could not obtain useful constraints for the duration of ingress or egress. We therefore assumed a fixed value of 40 s, which is the time it takes to eclipse the entire white dwarf (Bailey & Cropper, 1991). Although not apparent in the UV light curve (Fig. 2), the egress near phase 2, just before the flare, can be clearly identified. About 10 s after the beginning of the new OM exposure, the count rate increased quickly from 0 to $`3`$s-1. This fast rise coincides with the beginning of the flare in the X-ray light curve. The egress duration of $`3\pm 3`$s in the UV is similar to that measured in the optical for the main accretion region ($`2`$s; Bailey & Cropper, 1991). Table 1 shows the fitted times/phases of mid-ingress and mid-egress relative to the eclipse centres. Using high time-resolution optical light curves of the eclipse, Perryman et al. (2001) resolved two distinct accretion regions. They found that the egress of the main region (Spot 1) occurs at phase 0.032, while that of the secondary region (Spot 2) occurs at phase 0.027. The beginning of the flare at phase 0.0307 is in good agreement with the egress of spot 1. This result strongly suggests that the flare was caused by accretion on to the main region and that the onset of the flare was hidden from us by the eclipse. In agreement with this interpretation, the end of the flare (Fig. 2) is seen at the same orbital phase of 0.15 at which the main accretion region disappears behind the white dwarf.

### 2.3 Flare spectrum

As shown in the previous sections, the X-ray/UV flare near orbital phase 2 was most likely caused by a significant increase in the rate of accretion on to the white dwarf. According to the picture of the standard accretion column, the X-ray spectrum of the flare should therefore consist of two spectral components, a bremsstrahlung component from the shock-heated gas and a blackbody component due to reprocessing of hard X-rays in the white dwarf’s atmosphere. We find that the flare spectrum (Fig. 4) is consistent with a pure bremsstrahlung model but does not show statistically significant evidence for a blackbody component. In particular, we do not find the 28 eV blackbody that dominated the ROSAT spectrum during the high state (Ramsay et al., 1993), despite the similar bremsstrahlung fluxes during both observations (Section 2.1). It is plausible that, during the XMM-Newton observation, accretion occurred over a larger area, leading to a blackbody temperature too low for a detection by the EPIC instruments. A lower blackbody temperature is also expected in the absence of blobs in the accretion stream. Such blobs were likely present during the ROSAT observation, as Ramsay et al. (1993) found a large soft X-ray excess. Dense blobs impacting the white dwarf are thought to heat small areas of the photosphere to higher temperatures than the reprocessing of hard X-rays (Hameury & King, 1988). Without blobs in the accretion stream, the blackbody temperature would be much lower than the 28 eV found for the high state. As discussed later, a blackbody component with a temperature below $`12`$eV was possibly detected in the UV. The hardness ratio in Fig. 2 indicates some spectral variability, which is likely due to irregularities in the accretion stream.

The X-ray spectrum of the flare is well fit by a single-temperature thermal bremsstrahlung model ($`\chi _{red}^21`$). We included interstellar absorption in the model fit but could not obtain useful constraints for the hydrogen column density $`N_H`$. An upper limit of $`1\times 10^{20}`$cm-2 follows from the total galactic hydrogen column density in the direction of UZ For (Dickey and Lockman 1990; Stark 1988). Osborne et al. (1988) constrained $`N_H`$ to the range $`2.5\times 10^{19}`$$`1.3\times 10^{20}`$cm-2. Table 2 shows the fitted temperature and unabsorbed bolometric flux of the bremsstahlung model for various fixed column densities. The temperature of $`7`$keV is typical for bremsstrahlung from the accretion column and similar to that found by Rousseau et al. (1996) from the cyclotron spectrum. For a distance $`d=200`$pc we obtain an unabsorbed bolometric luminosity $`L_{brems}2.4\times 10^{30}\mathrm{erg}\mathrm{s}^1`$ (we used a geometric factor of $`3.1\pi `$ as in Section 2.1). This luminosity corresponds to a mass accretion rate $`\dot{M}2\times 10^{13}\mathrm{g}\mathrm{s}^1`$ ($`M_{WD}=0.7M_{}`$, $`R_{WD}=8\times 10^8`$cm). A constant temperature model likely oversimplifies the flare spectrum, which shows considerable temporal variability (see hardness ratio in Fig. 2). However, the low signal-to-noise ratio of the data did not allow us to obtain useful estimates for the temperature variations during the flare. Since $`80`$ per cent of the flux in the constant temperature model is inside the 0.15–12 keV energy range and directly detected by XMM-Newton, a constant temperature fit should provide a good estimate for the bolometric bremsstrahlung flux, despite the spectral variability.

In the Optical Monitor, the flare (Fig. 2) caused an average count rate increase of $`2`$s-1, which corresponds to a flux of $`1\times 10^{15}\mathrm{erg}\mathrm{cm}^2\mathrm{s}^1`$Å-1 at $`270`$nm or an integrated flux in the 240–340 nm band of the UVW1 filter of $`F_{UVW1}1\times 10^{12}\mathrm{erg}\mathrm{cm}^2\mathrm{s}^1`$. The latter value is comparable to the bolometric bremsstrahlung flux $`F_{brems}`$. The likely origin of the UV emission is either cyclotron radiation from the post-shock region or blackbody radiation from the photosphere below the shock (or a combination of the two). The XMM-Newton data do not allow us to distinguish between the two possibilities.

UZ For is known to be a strong emitter of cyclotron radiation (Schwope et al., 1990). Past observations of the cyclotron spectrum covered only optical wavelengths, but model calculations by Rousseau et al. (1996) suggest that cyclotron emission lines may be present in the near UV. If the Optical Monitor did indeed detect cyclotron radiation, only a small fraction of the total cyclotron flux was actually seen. Using the results by Woelk & Beuermann (1996) in a similar way as in Section 2.1, we estimate that the cyclotron luminosity was $`10^1`$$`10^2`$ times higher than $`L_{brems}`$. This luminosity corresponds to an accretion rate of $`10^{14}`$$`10^{15}\mathrm{g}\mathrm{s}^1`$ or a total accreted mass of $`10^{17}`$$`10^{18}`$g for the duration of the flare.

The accretion region is heated by reprocessing of X-rays from the post-shock region or by blobs in the accretion stream that penetrate several scale heights into the atmosphere. This gives rise to blackbody radiation visible at X-ray and UV energies (Warner, 1995). Assuming that all of the UV flare emission is due to this blackbody, we estimate a luminosity of $`4\times 10^{27}\mathrm{erg}\mathrm{s}^1`$Å-1 at 270 nm ($`d=200`$pc). As in Section 2.2 we used a geometric factor of 9.2. Since the blackbody component is not visible in the X-ray spectrum, we can place an upper limit of 7–12 eV on its temperature. The uncertainty derives mostly from the not well known hydrogen column density $`N_H=0.25`$$`1.3\times 10^{20}`$cm-2 (Osborne et al., 1988). (The EPIC MOS spectra in Fig. 4 show a small excess over the bremsstrahlung model at low energies. This might be considered marginal evidence for the blackbody component but could also be due to calibration uncertainties below 0.2 keV.) The limit on the blackbody temperature also puts constraints on the effective emitting area $`f_{eff}0.01`$ and the blackbody luminosity $`L_{bb}0.31.5\times 10^{33}\mathrm{erg}\mathrm{s}^1`$. The actual luminosity is likely close to this limit, since a significantly smaller $`L_{bb}`$ would require an unreasonably large emitting area $`f0.01`$. In particular, it is not possible that $`L_{bb}L_{brems}`$, which rules out reprocessing of X-rays as a major contributor to the UV flux. The large ratio $`L_{bb}/L_{brems}10^2`$ could be due to blobs in the accretion stream (Frank et al., 1988). Then, however, one might expect $`f_{eff}`$ to be smaller than our limit of 0.01. An effective emitting area $`f_{eff}<10^4`$ is required for the 28 eV blackbody found by Ramsay et al. (1993) in the ROSAT spectrum. (As discussed in Hameury & King (1988), $`f_{eff}`$ can be substantially smaller than the size of the accretion region.) Considering these difficulties, it seems unlikely that a large fraction of the UV flare emission is blackbody radiation from the accretion region.

## 3 Discussion

UZ For was found in a peculiar state that could have been classified as a regular low state if XMM-Newton had not detected several X-rays/UV flares. Whereas flaring is common during the high state and probably due to dense blobs in the accretion stream (Hameury & King, 1988), flaring during a low state has only been reported for two other polars. EUVE detected two transient events, each lasting $`1`$hr, during a low state of QS Tel (Warren et al., 1993). Three possible explanations for the transients were suggested: magnetic flares on the secondary; dense filaments in the accretion stream impacting the white dwarf; intermittently enhanced mass transfer from the secondary caused by magnetic flares or coronal mass ejections. In AM Her, Shakhovskoy et al. (1993) observed an optical transient that lasted $`20`$min and increased the brightness by 2 mag. The properties of the transient (exponential decay, blue color, no polarization) strongly suggest emission from a stellar flare on the companion star. Also in AM Her, Bonnet-Bidaud et al. (2000) observed several smaller transients ($`0.5`$mag) that were red in color and circularly polarized. These transients were likely caused by irregular accretion on to the white dwarf (e.g. by filaments in the accretion stream).

Our analysis of the XMM-Newton data clearly shows that the flaring in UZ For was caused by accretion on to the white dwarf. This rules out magnetic flares on the secondary as the source of the X-ray and UV emission. It is also unlikely that dense filaments in an otherwise constant accretion stream caused the large increase in X-ray luminosity. The observed bremsstrahlung spectrum and the absence of a blackbody component indicate that the shock was not buried in the photosphere, as is expected for dense filaments impacting the white dwarf. The most likely explanation for the flaring in UZ For is an intermittent increase in the mass transfer rate from the companion star.

The causes of high and low states in polars are still poorly understood. Unlike other CVs, polars do not have an accretion disk, so that changes in the rate of accretion on to the white dwarf must be due to variations of the mass flow rate through the L1-point. Variations in the size of the secondary or of the Roche-lobe cannot be the cause of low states, as those changes occur on time-scales of $`10^4`$$`10^5`$yr (King et al., 1995; Ritter, 1988). However, Howell et al. (2000) showed that for secondaries in short-period CVs the critical Roche-surface is well above the photosphere, so that mass flow rates trough the L1-point are controlled by the chromosphere. Changes in an active chromosphere can occur on time-scale short enough to explain the low states in polars. Alternatively, it has been suggested that star spots may appear at the L1-point and greatly reduce the mass transfer rate (Livio & Pringle, 1994; King & Cannizzo, 1998; Hessman et al., 2000). However, Howell et al. (2000) point out that, for short-period CVs, star spots exist well below the critical Roche surface, and the magnetic activity near the spots might actually increase the accretion rate. In both scenarios, it is possible that coronal mass ejections or solar flares near the L1-point intermittently increase the mass transfer rate, thus causing accretion events as those seen in UZ For.

Little is known about the level of stellar activity in CV secondaries. But because of their rapid rotation, which is synchronized with the binary orbital motion, an active atmosphere can be expected. A number of observations suggest that star spots and magnetic flares do exist on CV secondaries (e.g. Howell et al., 2000; Webb et al., 2002; Shakhovskoy et al., 1993). We showed that the X-ray transient in UZ For was caused by accretion of $`10^{17}`$$`10^{18}`$g of gas on to the white dwarf. It is not unreasonable that this much mass was ejected by a flare on the dM4.5 companion star. From observations of stellar flares on nearby M dwarfs, Pallavicini et al. (1990) estimated a density range of $`10^{12}`$$`10^{13}`$cm-3 and a volume range of $`10^{27}`$$`10^{29}`$cm3, corresponding to flare masses of $`10^{15}`$$`10^{18}`$g. For the largest flares, peak X-ray luminosities up to $`10^{30}\mathrm{erg}\mathrm{s}^1`$ (0.05–2 keV) were measured. A flare with such a high luminosity would have been detected by XMM-Newton with an EPIC count rate of $`0.2`$s-1 over $`10^3`$s. Since no emission from the companion star was seen, the X-ray luminosity of any stellar flare must have been below $`\mathrm{\hspace{0.25em}10}^{29}\mathrm{erg}\mathrm{s}^1`$. Coronal X-ray emission other than that from large stellar flares is unlikely to be detectable, even with the data integrated over the entire observation. For the most active stars of type dM5 or later, Fleming et al. (1993) measured average X-ray luminosities up to $`10^{28}\mathrm{erg}\mathrm{s}^1`$ (0.1–2.4 keV). The corresponding EPIC count rate of $`0.002`$s-1 is slightly below the detection threshold for the UZ For observation. We conclude that, despite the absence of X-ray emission from the companion star, accretion of gas ejected by stellar flares is a viable explanation for the transient events observed in UZ For.

## 4 Conclusion

During the XMM-Newton observation, UZ For was found in an extremely low accretion state with an X-ray luminosity $`800`$ times fainter than during a high state previously observed with ROSAT. Occasional X-ray and UV flaring was detected by the X-ray instruments and the Optical Monitor. The largest flare lasted $`900`$s and increased the X-ray flux by a factor of $`30`$. We found that the beginning of this flare coincided to within a few seconds with the eclipse egress of the main accretion region. This provides strong evidence that the flaring was caused by accretion on to the white dwarf. The X-ray spectrum of the flare is consistent with $`7`$keV thermal bremsstrahlung from the accretion column. A blackbody component, as seen with ROSAT during the high state, was not found. It is plausible that, because of a larger accretion region or the absence of blobs in the accretion stream, the blackbody temperature was too low for a detection by the X-ray instruments. The increase in the UV flux seen during the flare was probably caused by cyclotron radiation from the accretion column. A significant contribution of blackbody radiation to the UV flare emission is unlikely as this would require a very large soft excess. Under the assumption that all accretion energy is emitted as bremsstrahlung, we estimate an accretion rate of $`2\times 10^{13}\mathrm{g}\mathrm{s}^1`$ during the flare. However, since at this low rate most of the energy is emitted as cyclotron radiation, the actual accretion rate was probably $`10^1`$$`10^2`$ times higher. We therefore estimate that during the flare a total of $`10^{17}`$$`10^{18}`$g of gas was accreted on to the white dwarf. The likely cause of the flaring observed in UZ For is stellar activity on the companion star that intermittently increased the mass transfer rate near the L1-point. The mass that was accreted on to the white dwarf during the large transient is consistent with the mass ejected by a stellar flare.

Before and after the X-ray transient, extremely weak X-ray emission, possibly due to the regular low-state accretion on to the main region, was detected. The observed X-ray luminosity corresponds to an accretion rate of $`6\times 10^{11}\mathrm{g}\mathrm{s}^1`$. Since cyclotron radiation dominated the energy output, the actual accretion rate was probably $`10^{13}`$$`10^{14}\mathrm{g}\mathrm{s}^1`$, which is similar to the rates estimated for previous low states. In addition to the flare emission, we detect a roughly constant UV flux consistent with blackbody radiation from a 11000–K white dwarf. A small orbital modulation of the UV flux indicates the presence of a large, heated pole cap around the main accretion region.

Flaring during a low state has only been observed for a few polars. This may be due to insufficient monitoring of polars in low states. Yet low-state observations are essential since, during high and intermediate states, flaring caused by stellar activity on the companion star is likely overlooked and mistakenly attributed to accretion stream instabilities. Future low-state observations will reveal how common flaring due to stellar activity is among polars. An interesting question that might also be answered is, whether flaring occurs preferentially at the beginning or end of low states. Monitoring of irregular accretion in low-state polars may provide a new way to study stellar flares or other types of mass ejections that are too faint to be observed directly.

## acknowledgements

This work is based on observations obtained with XMM-Newton, an ESA science mission with instruments and contributions directly funded by ESA Member States and the USA (NASA). The authors acknowledge support from NASA grant NAG5-7714.


format: text

text:
# Gene networks from DNA microarray data: centrality and lethality

## Abstract

We construct a gene network based on expression data from DNA microarray experiments, by establishing a link between two genes whenever the Pearson’s correlation coefficient between their expression profiles is higher than a certain cutoff. The resulting connectivity distribution is compatible with a power-law decay with exponent $`\gamma 1`$, corrected by an exponential cutoff at large connectivity. The biological relevance of such network is demonstrated by showing that there is a strong statistical correlation between high connectivity number and lethality: in close analogy to what happens for protein interaction networks, essential genes are strongly overerpresented among the hubs of the network, that is the genes connected to many other genes.

Dipartimento di Scienze e Tecnologie Avanzate, Università del Piemonte Orientale, and INFN, gruppo collegato di Alessandria, Italy

and

Dipartimento di Fisica Teorica, Università di Torino, Via P. Giuria 1, 10125 Torino, Italy

e-mail: provero@to.infn.it

DNA microarray experiments are one of the most powerful tools for studying interactions between genes on the scale of the whole genome. It is widely believed that a huge amount of biologically relevant information is encoded in the results of such experiments, and that new analytical methods need to be developed to extract it.

In this work we propose to analyse the expression data obtained in microarray experiments by constructing a network of coregulated genes: the genes are the nodes of the network, and a link is established between two nodes whenever they are similarly expressed across many experimental conditions. On one hand, we show that such network, like many other known networks of self-organizing origin, shows a connectivity distribution that decays with a power law corrected, for large values of the connectivity, by an exponential cutoff . On the other hand, we show that the network encodes biologically relevant information in its topology: exactly as in the case of the protein interaction network in yeast , centrality is strongly correlated to lethality. In other words, among the genes that have the highest connectivity in the network, essential genes, whose deletion produces an inviable mutant, are strongly overrepresented.

We will work on yeast (S. cerevisiae), and use the expression data made publicly available by the authors of Ref. (and that include also the data obtained by the authors of Ref.), who performed a series of microarray experiments with the goal of identifying cell-cycle regulated genes. The data consist in the expression profiles for virtually all of the $`6200`$ yeast genes across a total of 77 timepoints.

The network is constructed by the following procedure:

1. To each gene we associate its expression profile defined as a string of 77 real numbers, representing as it is customary the $`log_2`$ of the ratio between expression (quantity of mRNA) at the given time-point and a reference value of the expression. The data come from different experiments and have been processed by the authors of Ref., to which we refer for details, so as to be comparable to each other across the various experiments.
2. Missing values are replaced with the average expression over the available timepoints. To prevent this manipulation from having a sizable effect on the construction of the network, we retain only those genes for which at least 70 timepoints out of 77 are available. We are thus left with 5293 genes as the nodes of our network.
3. We compute the Pearson’s correlation coefficient $`r`$ for all pairs of nodes in the network.
4. We establish a link between two nodes whenever $`r`$ is larger than a certain cutoff $`C`$.

The only free parameter in the procedure is the cutoff $`C`$. A possible way of choosing it is to compare the probability distribution of $`r`$ for the actual data to the same distribution after the data have been randomized by shuffling the expression values of each gene. For the randomized data, no pair of genes shows a value of $`r`$ greater than $`0.67`$: therefore by choosing $`C=0.67`$ the links we create can be considered of biological origin. A similar procedure was used in Ref., where a network was constructed by establishing a link between two genes whenever the effects of their deletion on the expression of the rest of the genome were linearly correlated.

With this choice of the cutoff $`C`$, 17643 links are established between the genes. Defining the connectivity $`k`$ of a node as the number of links departing from it, we have an average connectivity $`k6.67`$. Defining $`N(k)`$ as the number of nodes with connectivity $`k`$, we see that $`N(k)`$ shows a long tail that reaches up to $`k=173`$: the shape of the distribution is compatible with a power law decay with exponential cutoff:

$$N(k)=ak^\gamma \mathrm{exp}\left(\frac{k}{k_c}\right)$$
(1)
This is shown in Figs. 1-3: Fig. 1 shows $`N(k)`$ in logarithmic scale as a function of $`k`$ (noise in the data has been reduced by logarithmic binning), showing that the decay of $`N(k)`$ is slower than exponential for small to moderate values of $`k`$. Fig. 2 shows the same data with logarithmic scale on the $`k`$ axis too, and demonstrates that the decay is faster, at large $`k`$, than the pure power law characteristic of scale-free networks. Finally Fig. 3 shows the data after correction with the exponential cutoff, with $`k_c38`$. The slope of the straight line is $`\gamma =0.95`$. An analysis of the cumulative distribution confirms these results. Interestingly, a recent study of the transcriptional regulation network in yeast also revealed a scale-free network with $`\gamma 1`$ .

In this paper, we are mainly concerned with establishing the biological relevance of this network, independently of any of its graph-theoretic features. This we will do by showing that the nodes in the network with high connectivity are more likely to be essential genes, whose deletion produces an inviable mutant.

A list of essential yeast genes is publicly available from the Saccharomyces Genome Deletion Project , and comprises 1104 genes, corresponding to 18.7% of the genes deleted in the project. Of the 5293 genes in our network, 964 (18.2%) are included in the list. Fig. 4 shows the fraction $`f(k_{min})`$ of essential genes among the genes having connectivity $`k_{min}`$ or more, as a function of $`k_{min}`$: it grows from $`0.182`$ at $`k_{min}=0`$ (by definition) to $`1`$ for the the 6 most connected genes ($`k155`$).

The figure shows that essential genes are more and more overrepresented as the minimum connectivity $`k_{min}`$ is increased. To evaluate the statistical significance of such overrepresentation, suppose that the number of genes with connectivity $`kk_{min}`$ is $`n`$, and that among these $`m`$ are essential. Then one can evaluate the probability $`P(N,M;n,m)`$ that, out of $`n`$ randomly chosen nodes out of a set of $`N`$, $`m`$ or more are essential genes, when the total number of essential genes is $`M`$. This probability can be computed as the right tail of the appropriate hypergeometric distribution, and reaches very small values: for example the fraction of essential genes reaches 50% for $`k_{min}=37`$, with $`m=127`$ essential genes out of $`n=251`$ nodes, and the probability $`P(N=5293,M=964;n=251,m=127)`$ is about $`410^{33}`$.

In conclusion, we have built a gene network based on expression data obtained with DNA microarray experiment, by joining genes showing similar expression profiles. The resulting network shows a power law decay of the connectivity distribution with an exponential cutoff, and exponent $`\gamma 1`$. Its biological relevance is proved by the strong statistical correlation between centrality and lethality.

Acknowledgements. I would like to thank M. Caselle and C. Herrmann for many useful conversations and suggestions. I am also grateful to M. Barthélémy for many useful comments about a previous version of this manuscript.


text:
# Statistics of Lead Changes in Popularity-Driven Systems

## Abstract

We study statistical properties of the highest degree, or most popular, nodes in growing networks. We show that the number of lead changes increases logarithmically with network size $`N`$, independent of the details of the growth mechanism. The probability that the first node retains the lead approaches a finite constant for popularity-driven growth, and decays as $`N^\varphi (\mathrm{ln}N)^{1/2}`$, with $`\varphi =0.08607\mathrm{}`$, for growth with no popularity bias.

Extremes are vitally important in science and engineering. These quantities are used to determine the likelihood of a rare event, such as the probability of failure of a space shuttle launch or of a dam in flood conditions. The theory of extreme statistics provides a powerful tool to understand such real-world situations Galambos . Extremes are also irresistible in everyday life – we are naturally drawn to compilations of various pinnacles of human endeavor, such as, for example, lists of the most beautiful people, the richest people, the most-cited scientists, athletic records, etc guiness .

This social perspective about extremes raises new questions for which much less is known compared to the magnitude of the extreme value itself MAA02 . For example, how does the identity of the leader – the individual who possesses the extreme value of a particular attribute – change as a function of time? What is the rate at which lead changes occur? What is the probability that a leader retains the lead as a function of time?

We address these questions within the framework of growing networks, where the relevant quantity is the node degree – the number of links that join to each node. We view the degree as quantifying the popularity (or wealth) of the node, and the leader is the node with the highest degree. We focus on generic network models with preferential attachment to already popular nodes yule ; simon ; BA ; KR01 ; DMS2 , and networks with random attachment. The former describe, for example, the distributions of biological genera, word frequencies, publications, urban populations, and income yule ; simon , and contemporary applications to collaboration networks and the World-Wide Web have been developed applics . It has been posited that a hallmark of such systems is “the rich get richer” – that is, more popular nodes tend to remain so simon ; BA . Our basic goal is examine the consequences of preferential and random attachment mechanisms in growing networks and to test whether the adage of the rich get richer really does apply.

The network grows by adding nodes, each of which links to a pre-existing node with an attachment rate $`A_k`$ that depends only on the degree $`k`$ of the target node. We choose $`A_k=k+\lambda `$ with $`\lambda >1`$ yule ; simon . For such networks, the degree distribution has an asymptotic power-law tail, $`N_kN/k^{3+\lambda }`$, where $`N_k`$ is the number of nodes of degree $`k`$ and $`N`$ is the total number of nodes simon ; KR01 . For $`\lambda \mathrm{}`$ the growth mechanism reduces to random attachment and the degree distribution is exponential.

Identity of the leader. We characterize the identity of each node by its index $`J`$. A node of index $`J`$ is the $`J^{\mathrm{th}}`$ one introduced into the network. To start with an unambiguous leader, the initial system contains $`N=3`$ nodes, with the initial leader having degree 2 (and index 1) and the other two nodes having degree 1. A new leader arises when its degree exceeds that of the current leader. For the linear attachment rate, $`A_k=k`$, the average index of the leader $`J_{\mathrm{lead}}(N)`$ saturates to a finite value of approximately 3.4 as $`N\mathrm{}`$ (Fig. 1). With probability $`0.9`$, the leader is from among the 10 earliest nodes, while the probability that the leader is not among the 30 earliest nodes is less than $`0.01`$. Thus only the very earliest nodes have appreciable probabilities to be the leader; the rich really do get richer. Similarly in the general case of $`A_k=k+\lambda `$, the average index of the leader also saturates to a finite value that is a continuously increasing function of $`\lambda `$.

For random attachment ($`A_k=1`$), the average index of the leader grows algebraically, $`J_{\mathrm{lead}}(N)N^\psi `$ with $`\psi 0.41`$ (Fig. 1). The leader is still an early node (since $`\psi <1`$), but not necessarily one of the earliest. For example, simulation indicates that for $`N=10^5`$ a node with index greater than 100 has a probability of approximately $`10^2`$ of being the leader. Thus, in random attachment, the order of node creation plays a significant but not deterministic role in the identity of the leader node.

The identity of the leader can be determined from the joint index-degree distribution. Let $`C_k(J,N)`$ be the average number of nodes of index $`J`$ and degree $`k`$. As shown in KR01 , for constant attachment rate, this joint distribution has the Poisson form,

$$C_k(J,N)=\frac{J}{N}\frac{|\mathrm{ln}(J/N)|^{k1}}{(k1)!}.$$
(1)
From this, the average index of a node of degree $`k`$ is

$$J_k(N)=\frac{\underset{1JN}{}JC_k(J,N)}{_{1JN}C_k(J,N)}=N\left(\frac{2}{3}\right)^k,$$
(2)
implying $`J_{\mathrm{lead}}(N)=N(2/3)^{k_{\mathrm{max}}}`$. We estimate the maximum degree from the extreme value criterion $`_{kk_{\mathrm{max}}}N_k(N)1`$. Using $`N_k(N)=N/2^k`$ KR01 , we find $`2^{k_{\mathrm{max}}}N`$, or $`k_{\mathrm{max}}\mathrm{ln}N/\mathrm{ln}2`$. Therefore

$`J_{\mathrm{lead}}(N)N^\psi ,\mathrm{with}\psi =2{\displaystyle \frac{\mathrm{ln}3}{\mathrm{ln}2}}\mathrm{0.415\hspace{0.17em}037},`$
in excellent agreement with our numerical results.

For the linear attachment rate, the joint index-degree distribution is KR01

$$C_k(J,N)=\sqrt{\frac{J}{N}}\left\{1\sqrt{\frac{J}{N}}\right\}^{k1},$$
(3)
from which the average index of a node of degree $`k`$ is $`J_k(N)=12N/[(k+3)(k+4)]`$. Since $`N_k(N)4N/k^3`$ for the linear attachment rate BA ; KR01 , the extreme statistics criterion $`_{kk_{\mathrm{max}}}N_k(N)1`$ gives $`k_{\mathrm{max}}\sqrt{N}`$. Therefore $`J_{\mathrm{lead}}(N)12N/k_{\mathrm{max}}^2=𝒪(1)`$ indeed saturates to a finite value. A similar result holds in the general case $`A_k=k+\lambda `$. Thus the leader is one of the first few nodes in the network.

Lead changes. We find that the average number of lead changes $`L(N)`$ grows logarithmically in $`N`$ for both the attachment rates $`A_k=1`$ and $`A_k=k`$ (Fig. 2). There is, however, a significant difference in the distribution of the number of lead changes, $`P(L)`$, at fixed $`N`$. For $`A_k=1`$, this distribution is sharply localized, with the average value $`L5.609`$ in a network of $`N=10^5`$ nodes, while the maximum number of lead changes in $`10^5`$ realizations was 16. On the other hand, for $`A_k=k`$, $`P(L)`$ has a significant tail and the maximum number of lead changes is 63. This longer tail in $`P(L)`$ for linear attachment stems from repeated lead changes among the two leading nodes. Even though the distribution is visually broader, the average number of lead changes, $`L5.096`$, is less than that for $`A_k=1`$. Related to lead changes is the number of distinct nodes that enjoy the lead over the history of the network. Simulations indicate that this quantity also grows logarithmically in $`N`$.

This logarithmic behavior can be easily understood for the attachment rate $`A_k=1`$. Here the number of lead changes cannot exceed the upper bound given by the maximal degree $`k_{\mathrm{max}}\mathrm{ln}N/\mathrm{ln}2`$. To establish the logarithmic growth in the general case we first note that when a new node is added, the lead changes if the leadership is currently shared between two (or more) nodes and the new node attaches to a co-leader. The number of co-leader nodes (with degree $`k=k_{\mathrm{max}}`$) is $`N/k_{\mathrm{max}}^{3+\lambda }`$, while the probability of attaching to a co-leader is $`k_{\mathrm{max}}/N`$. Thus the average number of lead changes satisfies

$$\frac{d}{dN}L(N)\frac{k_{\mathrm{max}}}{N}\frac{N}{k_{\mathrm{max}}^{3+\lambda }}.$$
(4)
Since the maximal degree $`k_{\mathrm{max}}`$ grows as $`N^{1/(2+\lambda )}`$, Eq. (4) reduces to $`dL(N)/dNN^1`$ and thus gives the logarithmic growth $`L(N)\mathrm{ln}N`$. This argument can be adapted to networks with arbitrary attachment rates (except those growing faster than linearly with $`k`$ KR01 ), and thus the growth law $`L(N)\mathrm{ln}N`$ is universal. This universality is reminiscent of the radius of random networks which typically are proportional to $`\mathrm{ln}N`$, independent of their construction mechanism.

Fate of the first leader. Figure 3 shows that the degree distribution of the first node depends on the initial conditions for the linear attachment rate; the same is true in the general case $`A_k=k+\lambda `$ while for $`A_k=1`$ the initial condition is asymptotically irrelevant.

We can determine the degree distribution of the first node analytically for the constant and linear attachment rates. (A similar approach is given in Ref. DMS2 ). Let $`P(k,N)`$ be the probability that the first node has degree $`k`$ in a network of $`N`$ links note2 . For $`A_k=k`$, this probability obeys the master equation

$$P(k,N+1)=\frac{k1}{2N}P(k1,N)+\frac{2Nk}{2N}P(k,N).$$
(5)
The first term on the right accounts for the situation when the earliest node has degree $`k1`$. Then a new node attaches to it with probability $`(k1)/2N`$, thereby increasing the probability for the node to have degree $`k`$. Conversely, with probability $`(2Nk)/2N`$ a new node does not attach to the earliest node, thereby giving the second contribution to $`P(k,N+1)`$.

The solution to the master equation (5) for the “dimer” initial condition $``$ $``$ is KR02

$$P(k,N)=\frac{1}{2^{2Nk1}}\frac{(2Nk1)!}{(Nk)!(N1)!}.$$
(6)
For $`N\mathrm{}`$, this simplifies to the Gaussian distribution

$$P(k,N)\frac{1}{\sqrt{\pi N}}e^{k^2/4N}$$
(7)
for finite values of the scaling variable $`k/N^{1/2}`$. Thus the typical degree of the first node is of the order of $`N^{1/2}`$; this is the same scaling behavior as the degree of the leader node. For the trimer initial condition (which we typically used in simulations) we obtained the degree distribution of the first node in the form of a series of ratios of gamma functions KR02 , in which $`P(k,N)`$ has an $`e^{k^2/4N}`$ Gaussian tail, independent of the initial condition. The degree of the first node also approximates that of the leader node MAA02  more and more closely as the degree of the first node in the initial state is increased.

Although $`P(k,N)`$ contains all information about the degree of the first node, the behavior of its moments $`k^a_N=k^aP(k,N)`$ is simpler to appreciate. To determine the moments, it is more convenient to construct their governing recursion relations directly, rather than to calculate the moments from $`P(k,N)`$. Using Eq. (5), the average degree satisfies the recursion relation

$$k_{N+1}=k_N\left(1+\frac{1}{2N}\right),$$
(8)
whose solution is

$$k_N=\mathrm{\Lambda }\frac{\mathrm{\Gamma }\left(N+\frac{1}{2}\right)}{\mathrm{\Gamma }\left(\frac{1}{2}\right)\mathrm{\Gamma }(N)}\frac{\mathrm{\Lambda }}{\sqrt{\pi }}N^{1/2}.$$
(9)
The prefactor $`\mathrm{\Lambda }`$ depends on the initial condition, with $`\mathrm{\Lambda }=2,8/3,16/5,\mathrm{}`$ for the dimer, trimer, tetramer, etc., initial conditions.

This multiplicative dependence on the initial condition means that the first few growth steps substantially affect the average degree of the first node. For example, for the dimer initial condition, the average degree of the first node is, asymptotically, $`k_N2\sqrt{N/\pi }`$. However, if the second link attaches to the first node, an effective trimer initial condition arises and $`k_N(8/3)\sqrt{N/\pi }`$. Thus small initial perturbations lead to huge differences in the degree of the first node.

An intriguing manifestation of the rich get richer phenomenon is the behavior of the survival probability $`S(N)`$ that the first node leads throughout the growth up to size $`N`$ (Fig. 4). For the linear attachment rate, $`S(N)`$ saturates to a finite non-zero value of approximately 0.277 as $`N\mathrm{}`$; saturation also occurs for the general attachment rate $`A_k=k+\lambda `$. Thus for these popularity-driven systems, the rich get richer holds in a strong form – the lead never changes with a positive probability.

For constant attachment rate, $`S(N)`$ decays to zero as $`N\mathrm{}`$, but asymptotic behavior is not apparent even when $`N=10^8`$. A power law $`S(N)N^\varphi `$ is a reasonable fit, but the local exponent is still slowly decreasing at $`N10^8`$ where it has reached $`\varphi (N)0.18`$. To understand the slow approach to asymptotic behavior, we study the degree distribution of the first node. This quantity satisfies the recursion relation

$$P(k,N)=\frac{1}{N}P(k1,N1)+\frac{N1}{N}P(k,N1)$$
(10)
which reduces to the convection-diffusion equation

$$\left(\frac{}{\mathrm{ln}N}+\frac{}{k}\right)P=\frac{1}{2}\frac{^2P}{k^2}$$
(11)
in the continuum limit. The solution is a Gaussian

$$P(k,N)=\frac{1}{\sqrt{2\pi \mathrm{ln}N}}\mathrm{exp}\left\{\frac{(k\mathrm{ln}N)^2}{2\mathrm{ln}N}\right\}.$$
(12)
Therefore the degree of the first node grows as $`\mathrm{ln}N`$, with fluctuations of the order of $`\sqrt{\mathrm{ln}N}`$. On the other hand, the maximal degree grows faster, as $`v\mathrm{ln}N`$ with $`v=1/\mathrm{ln}2`$, and its fluctuations are negligible.

We now estimate the large-$`N`$ behavior of $`S(N)`$ as $`_{kk_{\mathrm{max}}}P(k,N)`$. This approximation gives

$`S(N)`$ $``$ $`{\displaystyle _{v\mathrm{ln}N}^{\mathrm{}}}{\displaystyle \frac{dk}{\sqrt{\mathrm{ln}N}}}\mathrm{exp}\left\{{\displaystyle \frac{(k\mathrm{ln}N)^2}{2\mathrm{ln}N}}\right\}`$ (13)
$``$ $`N^\varphi (\mathrm{ln}N)^{1/2},`$
with $`\varphi =(v1)^2/20.097989\mathrm{}`$. The logarithmic factor leads to a very slow approach to asymptotic behavior.

The above estimate is based on the Gaussian approximate for $`P(k,N)`$ which is not accurate outside the scaling region, namely, for $`k\mathrm{ln}N+\sqrt{\mathrm{ln}N}`$. However, we can determine $`P(k,N)`$ exactly because its defining recursion formula, Eq. (10), is closely related to that of the Stirling numbers of the first kind $`\left[\genfrac{}{}{0pt}{}{N}{k}\right]`$ stir , and the solution for the dimer initial condition is $`P(k,N)=\left[\genfrac{}{}{0pt}{}{N}{k}\right]/N!`$. The corresponding generating function is

$`S_N(x)={\displaystyle \underset{k=1}{\overset{N}{}}}P(k,N)x^k={\displaystyle \frac{x(x+1)\mathrm{}(x+N1)}{N!}}.`$
Using the Cauchy theorem, we express $`P(k,N)`$ in terms of the contour integral $`S_N(x)/x^{k+1}`$. When $`N\mathrm{}`$, this contour integral is easily computed by applying the saddle point technique KR02 . Finally we arrive at Eq. (13) with the same logarithmic prefactor but with a slightly smaller exact exponent $`\varphi =1v+v\mathrm{ln}v0.08607`$.

In summary, lead changes are rare in popularity-driven network growth processes and leadership is restricted to the earliest nodes. With finite probability, the first node remains the leader throughout the evolution. For growth with no popularity bias, leadership is shared among a somewhat larger cadre of nodes. As a consequence the average index of the leader node grows as $`N^\psi `$ with $`\psi =\mathrm{0.415\hspace{0.17em}037}\mathrm{}`$. The possibility of sharing the lead among a larger subset of nodes gives a rich dynamics in which the probability that the first node retains the lead decays as $`N^\varphi (\mathrm{ln}N)^{1/2}`$ with $`\varphi =0.08067\mathrm{}`$.

We are grateful to NSF grant DMR9978902 for partial financial support of this research.


text:
# Charge transfer in molecular conductors — oxidation or reduction?

\[

## Abstract

We discuss the nature of charge transfer in molecular conductors upon connecting to two metallic contacts and imposing a voltage bias across them. The sign of the charge transfer (oxidation vs. reduction) depends on the position of the metal Fermi energy with respect to the molecular levels. In addition, the charge transfer depends on the strength of the coupling (chemisorption vs. physisorption) with the contacts. A convenient way to establish the nature and onset of the charge transfer and the corresponding features in the I-V is to draw an energy level diagram for each spin species. Starting from such a level diagram, we argue that transport in the Tour-Reed switching molecules, which consist of a central phenyl ring with a nitroamine redox center, involves the oxidation of a highest occupied molecular orbital (HOMO)-based level.

\]

In recent years, there have been many reports of conductance characteristics measured for individual molecules or small ensembles of oriented molecules . Understanding transport in molecular conductors is quite challenging because of its special status between solid state physics and molecular chemistry. On the one hand, the system differs from conventional solid-state devices where conduction is understood in terms of alignment of a band structure that is largely unaffected by charge transfer or the nature of the contacts. On the other hand, in contrast to conventional chemistry, the molecule is not an isolated unit in equilibrium; coupling with metallic contacts makes the system an open one, capable of freely exchanging a fractional amount of charge to and from the contacts. In addition in the presence of a voltage bias across the molecule, the system is out of equilibrium and equilibrium statistical mechanical techniques become inapplicable. Predicting I-V characteristics of molecular conductors involves therefore a suitable combination of techniques and insights both from chemistry as well as from mesoscopic physics.

We have recently developed a self-consistent, ab-initio procedure for calculating I-V characteristics of molecular conductors . In this procedure, we self-consistently combine the outputs of a standard quantum chemical software (GAUSSIAN98) with a non-equilibrium Green’s function (NEGF) formalism for transport in an open, voltage-biased system. The method takes into account the hybridization of molecular levels with the contacts exactly, and can be used both for weak and strong contacts. In this paper, we will focus on an intuitive picture that can be used to gain qualitative insight into the nature of conduction through a molecule. This picture is particularly useful when the contacts are weakly coupled, and issues related to level broadening are minor. The first step is to obtain the energy levels of an isolated molecule and locate the Fermi energy $`E_F`$ of the metallic contacts, as shown in Fig. 1 (a). Under bias the electrochemical potentials $`\mu _{1,2}`$ of the contacts separate as shown in Fig. 1(b). We wish to address the following questions in this paper: (a) which molecular levels are involved in the conduction process? (b) Does the molecule get oxidized ($`MM^+`$) or reduced ($`MM^{}`$) at steady state?

Equilibrium energy level diagram ($`V=0`$). We start by obtaining an equilibrium energy level diagram as in Fig. 1(a). The energy levels for a given molecule are obtained using a standard quantum chemical software like Hyperchem, MOPAC or GAUSSIAN98. The equilibrium Fermi energy of gold is known to be $`E_F=5.1`$ eV. It is important to note that semi-empirical programs like Extended Hückel Theory do not give the correct absolute value for molecular energy levels. In that case, the energy levels for gold must also be computed using the same program for consistency .

Which levels conduct?. Fig. 1 shows a schematic description of transport through a mesoscopic device. At equilibrium (Fig. 1(a)) the metal contacts impose their common electrochemical potential on the molecule, equal to the metal fermi energy $`E_F`$. Coupling of the molecule with the contacts broadens each highest occupied (HOMO)/lowest unoccupied (LUMO) molecular orbital-based level, so that the number of electrons below the fermi energy is now fractional (the deviation of this number from the number of electrons for the isolated molecule is a measure of charge transfer between the molecule and the contacts at equilibrium). When a voltage bias is imposed (Fig. 1(b)), the system is driven out of equilibrium and the contact potential splits: $`\mu _{1,2}=E_FeV/2`$. (Note that in contrast to our earlier convention of fixing energy levels and varying the contact potentials with bias by different amounts depending on the ratio of couplings , we now vary the electrochemical potentials as above and let the energy levels adjust accordingly). For adiabatic, reflectionless contacts and ballistic transport, the contact potentials are imposed separately on the right and left-moving electrons . As long as the two electrochemical potentials lie in the gap, all HOMO levels lie below both potentials and are filled at zero temperature, the LUMO levels lie above the potentials and are empty, and there is no current. When the bias is sufficient that a level is crossed by one of the electrochemical potentials, that level is emptied out (oxidized) from one side but filled (reduced) from the other, with a resultant onset of current. Thus only levels which lie between the two electrochemical potentials contribute to conduction at zero temperature (at finite temperature, this range is extended by the thermal energy $`\pm k_BT`$). The onset of current is determined by the proximity of the Fermi energy to the molecular level closest to it.

Oxidized or reduced?. The coupling to the contacts of an extended molecule is described by a non-hermitian self-energy matrix $`\mathrm{\Sigma }`$ calculated in an appropriately chosen basis . This self-energy effectively partitions the problem into the molecular part and the contact, and leads to a conceptual simplification of the problem. The broadenings of the levels (inverse lifetime) are described by the anti-hermitian components $`\mathrm{\Gamma }_{1,2}`$ of $`\mathrm{\Sigma }_{1,2}`$. For fermi functions $`f_{1,2}`$ corresponding to the two contact potentials $`\mu _{1,2}`$, the steady-state occupancy of a molecular level is given by $`f=(f_1\mathrm{\Gamma }_1+f_2\mathrm{\Gamma }_2)/(\mathrm{\Gamma }_1+\mathrm{\Gamma }_2)`$. If the level lies between $`\mu _{1,2}`$ at zero temperature, $`f_1=0`$ and $`f_2=1`$ so that the occupancy is given by $`f=\mathrm{\Gamma }_2/(\mathrm{\Gamma }_1+\mathrm{\Gamma }_2)`$. One of the contacts is trying to add electrons (reduce), while the other is trying to remove them (oxidize). If the left contact is stronger than the right contact ($`\mathrm{\Gamma }_1\mathrm{\Gamma }_2`$), the level is emptied (oxidized) while in the opposite case the level stays filled (reduced). In either case, a current flows due to the competition between oxidation and reduction. Typically for a self-assembled monolayer, one end of the molecule (usually consisting of a thiol group chemisorbed on a gold surface) has a strong contact, while the other end is physisorbed and has a weak contact. This means that for a positive bias on the physisorbed end, the molecule remains neutral, while for opposite bias, charging occurs from the strongly contacted chemisorbed end .

For a high work function metal, $`E_F`$ usually lies near the HOMO level. This is an issue that clearly needs a lot more attention, since the precise location of $`E_F`$ depends on the model for the contacts, as well as the method of calculating the molecular Hamiltonian. A third terminal (gate)  can help resolve some of these issues. One of the consequences of a self-consistent evaluation of energy levels is that as we raise the fermi energy, the levels float up by an amount equal to the charging energy (for a solid, this is roughly given by the capacitance of a sphere of radius equal to the size of the electronic wave function). For a molecule with strong coupling to at least one contact, or for a molecule in solution, the effective electronic size is infinite, so that the charging energy is small, and classical STM theory holds. For a weakly coupled system, charging effects lead to the levels following the Fermi energy, so that the conduction tends to remain HOMO-based. We verified this for the phenyl-dithiol molecule (PDT) where we couple the outputs of a GAUSSIAN98  evaluation of energy levels using a Becke-3 parameter exchange and Perdew-Wang 91 correlation (B3PW91) approximation , with a non-equilibrium Green’s function (NEGF) evaluation of the density matrix , and iterate to convergence self-consistently . The basis set used was the Los Alamos National Lab set for effective core potentials of the double $`\zeta `$ type (LANL2DZ) . The converged energy levels increase with increasing $`E_F`$ at the rate of $`4`$eV/Volt, which corresponds to the capacitative charging energy of a sphere of size $`5\AA `$.

Level diagram under bias ($`V0`$). As an illustrative example, we have plotted the energy levels of the Tour-Reed molecule 2’-amino- 4,4’-di\[ethynyl\]phenyl-5’-nitro-1-\[thioacetyl\]benzene containing the nitroamine redox center , as a function of voltage applied to the sulphur end (Fig. 2). We consider an extended molecule including three surface gold atoms bonded to sulphur and obtain energy levels analogous to Fig. 1, including the effect of voltage-induced Stark shifts of the molecular orbitals. Exchange and correlation effects are included by performing a B3PW91 calculation with a LANL2DZ basis set using the GAUSSIAN 98 software. In the Tour-Reed experiment, one end of the molecule is self-assembled onto gold using sulphur “alligator clips”, leading to a strongly coupled chemisorbed contact, while the other end is a carbon atom weakly coupled to an evaporated gold contact. Charging occurs when the electrochemical potential $`\mu _1`$ of the sulphur end (dotted line in Fig. 2) crosses a molecular level. From Fig. 2, we see that for positive bias on sulphur, we first cross a HOMO level by the stronger sulphur contact, which oxidizes the molecule (the first LUMO level is actually a nonconducting gold-based level, as we discuss below). For negative bias, the HOMO level is crossed by the weaker carbon contact, so the level stays filled and there is no charging. This means that the I-V should be strongly asymmetric, as is indeed observed. ***For a spatially symmetric molecule, such contact-induced charging asymmetries alone determine the asymmetry in the I-V . However for spatially asymmetric molecules such as the Tour-Reed example here, the I-V asymmetry is much larger, and given predominantly by the unequal Stark shifting of the intrinsic molecular wavefunctions (the asymmetric Stark shifting is seen in Fig. 2). The sulphur contact potential $`\mu _1`$ crosses the HOMO level at around 2 V, whereupon the molecule is oxidized from $`M`$ to $`M^+`$. This is important, since the precise nature of the charging process is an important question raised in the literature . Although cyclic voltammetric curves seem to suggest reduction processes as being operative, one has to remember that in such measurements counterions are present in solution, and substantially compromise any charging effects that may dominate in solid state.

When $`M`$ changes to $`M^+`$, we need a whole new set of energy levels (Fig. 2). Two changes happen: (i) since one electron of a particular spin is pulled out by the contact, the spin degeneracy is removed; (ii) secondly, the removal of one electron of given spin leads to charging effects which cause all the levels of $`M^+`$ obtained from GAUSSIAN98 to float down by about 2 volts. This is unphysical, because GAUSSIAN98 models an isolated molecule with image charges at infinite distance. In the presence of image charges on the contacts, the charging energy will be lower. We therefore raise the energy levels of $`M^+`$ so as to align the energy level of the semi-occupied molecular orbital (SOMO) with the HOMO of the neutral species ($`M`$). Such an alignment is motivated by our assumption that in the absence of correlations and molecular reorganization, the ionization potential of $`M`$ equals the electron affinity of $`M^+`$. This assumption is a point that requires further work, and should come automatically out of a proper ab-initio theory that includes the effect of image charges on the contacts. Note that in order to obtain spin-dependent charging effects that shift all levels except the SOMO substantially, one needs to do a spin-unrestricted (UHF/LSDA) calculation that breaks the spin degeneracy of the energy levels.

It is interesting to note that for an extended molecule, there are some gold-based levels. These can be identified either from their wave functions that are localized strongly on gold, or from their energy levels; the levels run roughly parallel to the electrochemical potential $`\mu _1`$ of the contact on the chemisorbed side (dotted line in Fig. 2). For small molecules these gold levels, as well as metal-induced gap states (MIGS) that arise when a self-energy is included in our calculation, end up contributing to the conductance through direct tunneling processes. For large molecules such states hybridize with molecular wave functions, but don’t contribute directly to tunneling, so they can be ignored. However the gold levels can still influence the positions of the other levels through level-repulsion, so the contacts affect the voltage at which a HOMO level is crossed and oxidation occurs.

We have shown that in self-assembled monolayers such as the Tour-Reed molecules containing the nitroamine redox center, charge transfer occurs due to oxidation from the strongly contacted end involving the level closest to the Fermi energy. We are currently investigating the mechanism of switching thereafter. Understanding the behavior of a molecular switch is important from the point of view of molecular electronics. Several switching molecules have been reported , and various possible mechanisms suggested . A charged species could deform in an applied field, and the electromechanical effect could switch off the current. If the added charge is localized on a sidechain of the molecule, it could rotate the central benzene ring attached to it in the presence of the field, distortiing the conjugated electron wave function enough to switch current off. A more complicated cooperative effect involving the whole self-assembled monolayer could also arguably lead to a current switch. Prior to a full self-consistent ab-initio calculation of transport properties including contact effects, it is important to develop an intuitive understanding of transport in the isolated species. An energy level diagram as illustrated above can serve as a quick intuitive tool to answer some basic questions such as the onset of current, whether an electron is being added or removed, or the existence of an asymmetry in the I-V.

We would like to thank Prashant Damle and Deepak Singh for useful discussions. This work was supported in part by the Defense Advanced Research Projects Agency (DARPA) and the US Army Research Office (ARO) under grants number DAAH04-96-1-0437 and DAAD19-99-1-0198.


text:
# Stability of strangelet at finite temperature

## I Introduction

The study of small lumps of strange quark matter, called strangelets, plays an important role for research the quark gluon plasma (QGP) in recent relativistic heavy ion collision (RHIC) experiments. The reason is that although many signatures of QGP such as $`J/\mathrm{\Psi }`$ suppression, strangeness enhancement, thermal dilepton electro magnetic radiation, etc., have been found , but it is still ambiguous because these signatures can also be explained by hadron gas . To search an unambiguous signature of QGP is the key for RHIC experiments. The strangelet, as was argued by Greiner et. al , is a good candidate which could serve as an unambiguous signature for the QGP.

The essential problem for detectability of strangelet in RHIC experiments is to study its stability during the formation of QGP. Employing MIT bag model, many authors discussed this problem . They came to the same conclusion that the electric charge and the strangeness fraction are of vital important to the experimental searches of strangelets because these two elements affect a change of stability remarkably. But the effect of charge on stability is different for different authors. Jaffe and his co-workers  argued that the strangelet is slightly positive charged. They considered strong and weak decay by nucleon and hyperon emission together and concluded that the stable strangelets will have a low but positive change to mass ratio. Contrary to Jaffe et. al, after considering the initial condition of possible strangelet production in RHIC carefully, Greiner and his co-workers  argued that strangelets are most likely highly negative charged. Here we hope to emphasize that above discussions are limited in the framework of MIT bag model and at zero temperature.

Since the quark deconfinement phase transition can occur at high temperature or/and high density only, it is of interest to extend the investigation for the stability of strangelets to finite temperature. This is the objective of this paper.

However, MIT bag model is a permanent quark confinement model because the confined boundary condition does not change with temperature. In principle, one can not use this model to study the phase transition of QCD directly. The best that we can do is to employ a model which can almost reproduce the properties of strange quark matter obtained by MIT bag model, but in which the quark confinement is not permanent. The quark mass density- and temperature-dependent model (QMDTD)  suggested by us is one of such candidates.

Many years ago, a quark mass density-dependent model (QMDD) was suggested by Fowler, Raha and Weiner  and then it was employed by many authors to discuss the properties of strange quark matter . According to the QMDD model, the masses of $`u,d`$ quarks and strange quarks (and the corresponding anti-quarks) are given by

$`m_q`$ $`=`$ $`{\displaystyle \frac{B}{3n_B}},(q=u,d,\overline{u},\overline{d}),`$ (1)
$`m_{s,\overline{s}}`$ $`=`$ $`m_{s0}+{\displaystyle \frac{B}{3n_B}},`$ (2)
where $`n_B`$ is the baryon number density, $`m_{s0}`$ is the current mass of the strange quark and $`B`$ is the vacuum energy density. As was proved by ref., the properties of strange matter in the QMDD model are nearly the same as those obtained in the MIT bag model. In fact, it is not surprised if one notices that the confinement mechanism of MIT bag model is almost the same as that of QMDD model .

But when we employ the QMDD model to discuss the properties of strangelet, a lot of difficulties emerge . At first, the radius of strangelet decreases as the temperature increases; Secondly, it can not mimic the correct phase diagram of QCD because the temperature $`T`$ tends to infinite when $`n_B0`$. To overcome these difficulties, we suggest a QMDTD model . Instead of a constant $`B`$ in QMDD model, we argue that $`B`$ is a function of temperature and introduced an ansatz

$`B(T)`$ $`=`$ $`B_0\left[1a\left({\displaystyle \frac{T}{T_c}}\right)+b\left({\displaystyle \frac{T}{T_c}}\right)^2\right],0TT_c,`$ (3)
$`B(T)`$ $`=`$ $`0,T>T_c,`$ (4)
where $`B_0`$ is the vacuum energy density inside the bag (bag constant) at zero temperature, $`T_{c\text{ }}=170\text{MeV}`$ is the critical temperature of quark deconfinement phase transition, and $`a`$, $`b`$ are two adjust parameters. Since $`B`$ is zero when $`T=T_c`$, a condition

$$1a+b=0$$
(5)
is imposed and only one parameter $`a`$ can be adjusted. As pointed out by ref., in order to satisfy two physical conditions of strangelet at finite temperature, namely, (1), the radius of the strangelet must increases when the temperature rises, (2), the energy of the strangelet must increases when the temperature rises, the parameter $`a`$ is restricted in a small range:

$$0.65a0.8.$$
(6)
In this paper, we fix the value of two parameters $`a`$, $`b`$ in this suitable range

$$a=0.65,b=0.35.$$
(7)
With these parameters set, we use the QMDTD model to discuss the thermodynamical properties of strangelet, especially, to investigate the stability of strangelet via strong hadron emission and weak hadronic decay. We will study the effects of temperature, charge and strangeness fraction on the stability and hope that our study can have impact to the detectability of the strangelet.

The organization of this paper is as follows. In the following section, we give the formulae of thermodynamical calculations for QMDTD model. The results of thermodynamical properties of strangelet are presented in section 3. In section 4, we will discuss the stability of strangelet via possible strong and weak decays. The last section is a summary.

## II Thermodynamical formulae

To calculate the dynamical and thermodynamical quantities of the strangelet, we must look for the density of states first. The density of states of a spherical cavity in which the free particles be contained can be expressed as

$$\rho (k)=\frac{dN(k)}{dk},$$
(8)
where $`N(k)`$ is the total number of particle states and can be written in terms of dimensionless variable $`kR`$ as

$$N(k)=A(kR)^3+B(kR)^2+C(kR),$$
(9)
where $`R`$ is the radius of the bag. The three terms of the right hand side of Eq.(9) refer to the contributions of the volume, surface and curvature, respectively. The coefficients $`A`$, $`B`$ and $`C`$ are expected to be very slow varying functions of $`kR`$, and their expressions are model dependent. For MIT bag model, these coefficients have been obtained by numerical calculations in our previous paper . The volume term $`A`$ is a constant, and it has the value of

$$A=\frac{2g}{9\pi },$$
(10)
where $`g`$ is the total degeneracies. For example, it is the total number of spin and color degrees of freedom for a quark with flavor treated separately. The surface term $`B`$ is

$$B\left(\frac{m}{k}\right)=\frac{g}{2\pi }\left\{\left[1+\left(\frac{m}{k}\right)^2\right]\mathrm{tan}^1\left(\frac{k}{m}\right)\left(\frac{m}{k}\right)\frac{\pi }{2}\right\},$$
(11)
where $`m`$ is the mass of quark. Equations (10) and (11) are in good agreement with the ones given by multi-reflection theory . The curvature term $`C`$ can not be evaluated by this theory except for the two limiting case: $`m0`$ and $`m\mathrm{}`$. Madsen proposed that

$$\stackrel{~}{C}(\frac{m}{k})=\frac{g}{2\pi }\left\{\frac{1}{3}+\left(\frac{k}{m}+\frac{m}{k}\right)\mathrm{tan}^1\frac{k}{m}\frac{\pi k}{2m}\right\}.$$
(12)
But as was pointed out by ref., the beat fit of numerical data for curvature term is

$$C(\frac{m}{k})=\stackrel{~}{C}(\frac{m}{k})+\left(\frac{m}{k}\right)^{1.45}\frac{g}{3.42\left(\frac{m}{k}6.5\right)^2+100}.$$
(13)
Now we are in the position to calculate the thermodynamical quantities of strangelets for QMDTD model. The thermodynamical potential $`\mathrm{\Omega }`$ is

$$\mathrm{\Omega }=\underset{i}{}\mathrm{\Omega }_i=\underset{i}{}\frac{g_iT}{(2\pi )^3}_0^{\mathrm{}}𝑑k\frac{dN_i}{dk}\mathrm{ln}\left(1+e^{\beta (\epsilon _i(k)\mu _i)}\right),$$
(14)
where $`i`$ stands for $`u,d,s`$ (or $`\overline{u},\overline{d},\overline{s}`$ ) quarks, $`g_i=6`$ for quarks and antiquarks. $`{\displaystyle \frac{dN_i}{dk}}`$ is the density of states for various flavor quarks, it is given by eqs.(8)-(13). $`\mu _i`$ is the corresponding chemical potential (for antiparticle $`\mu _{\overline{i}}=\mu _i`$).

$$\epsilon _i(k)=\sqrt{m_i^2+k^2}$$
(15)
is the single particle energy and $`m_{i\text{ }}`$is mass for quarks and antiquarks.

According to the QMDTD model , the masses of quarks are

$`m_{u,\overline{u},d,\overline{d}}`$ $`=`$ $`{\displaystyle \frac{B_0}{3n_B}}\left[1a\left({\displaystyle \frac{T}{T_c}}\right)+b\left({\displaystyle \frac{T}{T_c}}\right)^2\right],0TT_c,`$ (16)
$`m_{u,\overline{u},d,\overline{d}}`$ $`=`$ $`0,TT_c,`$ (17)
$`m_{s,\overline{s}}`$ $`=`$ $`m_{s0}+{\displaystyle \frac{B_0}{3n_B}}\left[1a\left({\displaystyle \frac{T}{T_c}}\right)+b\left({\displaystyle \frac{T}{T_c}}\right)^2\right],0TT_c,`$ (18)
$`m_{s,\overline{s}}`$ $`=`$ $`0,TT_c,`$ (19)
where $`m_{s0}`$ is the current mass of the strange quark matter, $`n_B`$ is the baryon number density

$$n_B=A/V,$$
(20)
and $`A`$ is the baryon number of the strangelet, $`V={\displaystyle \frac{4}{3}}\pi R^3`$ is the volume of the strangelet. Using the standard statistical treatment, and noticing that $`\mathrm{\Omega }`$ is not only a function of temperature, volume and chemical potential, but also of density, it can be proved that the total pressure $`p`$ and the total energy density $`\epsilon `$ are given by

$$p=\frac{1}{V}\frac{(\mathrm{\Omega }/n_B)}{(1/n_B)}|_{T,\mu _i}=\frac{\mathrm{\Omega }}{V}+\frac{n_B}{V}\frac{\mathrm{\Omega }}{n_B}|_{T,\mu _i},$$
(21)
$$\epsilon =\frac{\mathrm{\Omega }}{V}+\underset{i}{}\mu _in_i\frac{T}{V}\frac{\mathrm{\Omega }}{T}|_{\mu _i,n_B}.$$
(22)
The number density of each particle can be obtained by means of

$$n_i=\frac{1}{V}\frac{\mathrm{\Omega }}{\mu _i}|_{T,n_B}.$$
(23)
At finite temperature, we must include the contributions of the anti-particles, therefore, the baryon number for $`i`$ quark is given by

$$\mathrm{\Delta }N_i=(n_in_{\overline{i}})\times V=\frac{g_i}{(2\pi )^3}_0^{\mathrm{}}𝑑k\frac{dN_i}{dk}\left(\frac{1}{\mathrm{exp}[\beta (\epsilon _i\mu _i)]+1}\frac{1}{\mathrm{exp}[\beta (\epsilon _i+\mu _i)]+1}\right).$$
(24)
The strangeness number $`S`$ of the strangelet reads

$$S=\mathrm{\Delta }N_s,$$
(25)
the baryon number $`A`$ of the strangelet satisfies

$$A=\frac{1}{3}(\mathrm{\Delta }N_u+\mathrm{\Delta }N_d+\mathrm{\Delta }N_s).$$
(26)
The electric charge $`Z`$ of the strangelet is

$$Z=\frac{2}{3}\mathrm{\Delta }N_u\frac{1}{3}\mathrm{\Delta }N_d\frac{1}{3}\mathrm{\Delta }N_s.$$
(27)
At finite temperature, the stability condition of strangelets for the radius reads

$$\frac{\delta F}{\delta R}=0.$$
(28)
where the free energy $`F`$ of strangelet is,

$$F=ET\stackrel{~}{S},$$
(29)
$`E=\epsilon V`$ is the total energy, and

$$\stackrel{~}{S}=\underset{i}{}\stackrel{~}{S}_i=\underset{i}{}\frac{\mathrm{\Omega }}{T}|_{\mu _i,n_B}$$
(30)
is the entropy.

Given strangeness number, baryon number and electric charge, for any strangelet, we can calculate chemical potentials for $`u,d,s`$ quarks self-consistently by the equations (25),(26),(27) and (28). Then the thermodynamic potential, free energy and the stable radius of the strangelet at finite temperature can be obtained self-consistently.

## III Thermodynamical properties of strangelet

The numerical calculations have been done with the parameters set:

$$B_0=170\text{MeV}\text{fm}^3,m_{s0}=150\text{MeV},T_c=170\text{MeV},$$
(31)
and the possible area for strangelet in our calculations is chosen as

$`S`$ $`>`$ $`0,`$ (32)
$`Z`$ $``$ $`A,`$ (33)
$`S+Z`$ $``$ $`2A.`$ (34)
We calculate the free energy of the strangelet first. Figure 1 shows the free energy per baryon $`F/A`$ in a $`SZ`$ plane for all possible strangelets with $`A=5`$ and $`T=20\text{MeV}`$. This is a set of downward protruding curves. The corresponding curves but for $`T=50\text{MeV}`$ are shown in figure 2. Comparing these two figures, we find that the positions of dots appeared in figure 2 are lower than the corresponding positions (with the same strangeness number $`S`$ and charge $`Z`$) in figure 1. For example, for $`S=3`$ and $`Z=2`$, the free energy per baryon of the dot is $`F/A=1064.9\text{MeV}`$ in figure 1, but reads $`F/A=898.9\text{MeV}`$ in figure 2. For fixed $`S`$ and $`Z`$, the free energy of the strangelet decreases when temperature increases. This result is reasonable because the entropy increases with the temperature, the term $`T\stackrel{~}{S}`$ in eq.(29) increases considerably and $`F`$ decreases. We will see in the next section this result affects on the decay of the strangelet remarkably.

To illustrate our result transparently, we draw the $`F/A`$ vs $`S`$ curves for fixed $`A=5,Z=5`$ but different temperatures $`T`$ $`=20\text{MeV}`$ and $`50\text{MeV}`$ in figure 3, respectively. We find that the minimum of these two curves are located at the same strangeness $`S=4`$ but with different free energies.. The whole curve for $`T=20\text{MeV}`$ is located on the upper position of the curve for $`T=50\text{MeV}`$. We also draw the $`F/A`$ vs $`S`$ curve for $`A=10,Z=10`$ and $`T=50\text{MeV}`$ in figure 4. Compared with figure 3, the shape of the curve in figure 4 almost does not changes except the position of the minimum is changed to $`S=8`$. These results are similar to that one given by .

Now we turn to study the charge and strangeness dependences of the radius of strangelet. We draw $`F/A`$ vs $`R`$ curves with fixed $`A=5`$ ,$`Z=1`$ and $`T=50\text{MeV}`$ but different $`S=2,6`$ and $`9`$ in figure 5, respectively. The stable radiuses given by Eq.(28) for different strangeness are different. Figure 5 shows that the stable radius changes from $`1.64\text{fm}`$ to $`1.615\text{fm}`$ when strangeness number changes from $`2`$ to $`9`$. The same curves for $`A=5`$ ,$`Z=4`$ and $`T=50\text{MeV}`$ but different $`S=2,14`$ are shown in figure 6. We see that the stable radius become $`1.67\text{fm}`$ for $`S=2`$ and $`1.62\text{fm}`$ for $`S=14`$. Therefore, we come to a conclusion that the stable radius of strangelet decreases when $`S`$ increases.

To illustrate the charge dependence of the stable radius, we draw $`F/A`$ vs $`R`$ curves with fixed $`A=5`$ ,$`S=2`$ and $`T=50\text{MeV}`$ but different electric charge $`Z=1,5`$ and $`8`$ in figure 7, respectively. The stable radius increases from $`1.64\text{fm}`$ to $`1.69\text{fm}`$ when the electric charge rises from $`1`$ to $`8`$. The same curves for $`A=5`$ ,$`S=9`$ and $`T=50\text{MeV}`$ but different charge $`Z=1,4`$ are shown in figure 8. We see that the stable radius changes from $`1.615\text{fm}`$ to $`1.605\text{fm}`$ when $`Z`$ decreases from $`1`$ to $`4`$. We find that the stable radius of strangelet increases with electric charge .

## IV Stability of the strangelet

In this section, we follow the line of ref. to investigate the stability of the strangelet and extend their study to finite temperature by using the QMDTD model.

### A Strong decay and unstable strangelet

As was pointed out by refs., small clusters of strange matter are most favoured for detection. As two examples, hereafter we study two cases with $`A=5`$ and $`A=10`$, respectively.

Instead of the binding energy at zero temperature, we calculate the free energy per baryon of the possible strangelet at finite temperature first. Figure 9 and 10 show the free energy per baryon of all possible strangelets as a function of the strangeness $`S`$ at the temperature $`T=50\text{MeV}`$ but for $`A=5`$ and $`A=10`$, respectively. The solid lines in these two figures connect the masses of the nucleon, $`\mathrm{\Lambda },\mathrm{\Xi }`$ and $`\mathrm{\Omega }`$. As a first cut for potential candidates of stable strangelets, these ones lying above this line can (or probably will) completely decay to the pure hadron state via strong processes, and only these ones beneath the line will be possible for metastable or stable strangelets .

We consider the influence of temperature now. As shown in last section, the free energy per baryon increases when temperature decreases. The positions of the dots will raise and many dots will across the line and become unstable when temperature decreases. Figure 11 shows this result clearly. The figure 11 is the same as figure 9 except for temperature $`T=20\text{MeV}`$. Comparing these two figures, we find many dots across the line when temperature decreases from $`50\text{MeV}`$ to $`20\text{MeV}`$. At low temperature, more strangelets can decay to the pure hadronic state by strong process.

Here we must emphasize that although any strangelet initially formed under this line can not decay to a pure hadron state, it is still possibly unstable because it can decay to a hadron and another strangelet changing baryon number, strangeness number and charge . We will look for possible strong decays, i.e. single baryon {n, p, $`\mathrm{\Lambda },\mathrm{\Sigma }^+,\mathrm{\Sigma }^{},\mathrm{\Xi }^0,\mathrm{\Xi }^{}`$ and $`\mathrm{\Omega }`$} emission and mesonic decays at the finite size configuration at finite temperature.

The baryon number $`A`$, strangeness number $`S`$ and electric charge $`Z`$ are conserved in the strong process. A general expression of a strong baryon decay for a strangelet $`Q(A,S,Z)`$ can be written as

$$Q(A,S,Z)Q(A1,SS_x,ZZ_x)+x(1,S_x,Z_x).$$
(35)
And this process is allowed if the free energy balance of the corresponding reaction is

$$F(A,S,Z)>F(A1,SS_x,ZZ_x)+m_{x,}$$
(36)
where $`F`$ stands for the total free energy of the strangelet and $`x`$ stands for a baryon with strangeness number $`S_x`$ and electric charge $`Z_x`$. Our results of numerical calculation are shown in figures 12, 13 and 14. For $`A=5`$ and $`T=50\text{MeV},`$ all strangelets including the ones lying above the lines of figure 9 and the ones which satisfy the inequality (36) are drawn by the circles in figure 12. These strangelets can strong decay and are unstable. The unstable strangelets for $`A=10`$ and $`T=50\text{MeV}`$ are drawn by circle in figure 13. Figure 12 and figure 13 show that the strangelets with small strangeness number $`S`$ and situated in the left side of the figure are unstable.

To study the temperature effect, the same figure for $`A=5`$ and $`T=20\text{MeV}`$ are shown in figure 14. Comparing figure 12 and figure 14, we find that many strangelets become unstable when temperature decreases. Lower temperature is in favour of the strong decay of the strangelet.

### B Weak decay and metastable strangelet

According to the definition of Greiner et. al , ”a strangelet is called metastable in the following if its energy lies under the correspond (free) hadronic matter of the same baryon number, charge and strangeness, and if it can not emit a single hadron or multiple hadrons by strong processes.” At finite temperature, instead of energy, we use free energy. A metastable strangelet can then only decay via weak process like the nonleptonic (hadronic) decays. In this sub-section, we study all possible weak hadronic decay for metastable strangelets.

In the weak process, the baryon number $`A`$ and the electric charge $`Z`$ are conserved, but the strangeness number $`S`$ is not conserved. For the weak decay, $`\mathrm{\Delta }S=\pm 1`$. Therefore, a general expression of a weak baryon decay for a strangelet $`Q(A,S,Z)`$ can be written as

$$Q(A,S,Z)Q(A1,SS_x1,ZZ_x)+x(1,S_x,Z_x).$$
(37)
This process is allowed if the free energy balance of the corresponding reaction is

$$F(A,S,Z)>F(A1,SS_x1,ZZ_x)+m_{x.}$$
(38)
The metastable strangelets which satisfy the inequality (38) are shown in figure 12, 13 and 14 with filled circle. We find that they are almost situated in the area with negative charge (or slightly positive charge) and high strangeness number. In particular, high strangeness is in favour of the stability of strangelet. As shown in figures 12, 13 and 14, a minimal strangeness $`S_c`$ above which ($`S>S_c`$) strangelets are metastable or stable exists. The values of the minimal strangeness $`S_c`$ are $`5`$ in figure 12 ($`A=5,T=50\text{MeV}`$) and $`10`$ in figure 13 ($`A=10,T=50\text{MeV}`$). Defining strangeness fraction $`f_{s\text{ }}`$ as

$$f_{s\text{ }}=\frac{S}{A},$$
(39)
Greiner and his co-workers had predicted that there exists a critical $`f_{sc\text{ }}`$at zero temperature and pointed that the metastable or stable strangelets could only be found above this value $`f_{s\text{ }}>f_{sc\text{}}`$Our results strongly support their prediction and extend it to the finite temperature. We obtained $`f_{sc}=1`$ at $`T=50\text{MeV}`$. This value does not depend on the baryon numbers. But as shown in figure 14, $`f_{sc}`$ depends on the temperature. It becomes $`f_{sc}={\displaystyle \frac{10}{5}}=2`$ at $`T=20\text{MeV}`$.

Now we turn to discuss the charge dependence of the stability of strangelet. As shown in figures 12, 13 and 14, higher negative charge is in favour of the metastable and/or stable strangelets. A maximal charge, which is $`Z_m=1`$ in figure 12 ($`A=5,T=50\text{MeV}`$), $`Z_m=6`$ in figure 13 ($`A=10,T=50\text{MeV}`$) and $`Z_m=3`$ in figure 13 ($`A=5,T=20\text{MeV}`$), respectively, exists. The strangelets would be stable or metastable when $`ZZ_m`$. Defining charge fraction $`f_z`$ as

$$f_{z\text{ }}=\frac{Z}{A},$$
(40)
We find that the maximal charge fraction depends on not only the temperature but also the baryon number $`A`$. For example, $`f_{zm}`$ equals to $`{\displaystyle \frac{1}{5}}`$ for $`A=5`$ and $`{\displaystyle \frac{6}{10}}`$ for $`A=10`$ when $`T=50\text{MeV}`$.

### C Stable strangelet

The strangelet which is stable against both strong and weak decay is called stable strangelet. The stable strangelets are shown in figure 12, 13 and 14 by filled squares. We find that only a few strangelets can not decay via strong and weak reactions and be completely stable. For example, as shown in figure 12, for $`A=5`$ and $`T=50\text{MeV}`$, in total $`128`$ strangelets, only $`4`$ strangelets are stable ($`Z=3,S=12`$), ($`Z=4,S=13`$), ($`Z=4,S=14`$) and ($`Z=5,S=15`$).

Finally, we hope to emphasize that the high negative charge and the high strangeness are in favour of the stable strangelets obviously at finite temperature. The stable strangelets in figures 12, 13 and 14 are all highly negative charged. The conclusion given by Greiner et. al at zero temperature  is still correct at finite temperature.

## V Summary

In summary, by using the QMDTD model, we have studied the thermodynamical properties and the stability of strangelets at finite temperature. We obtain:

(1), For fixed strangeness and charge, free energy per baryon decreases as the temperature increases. The stable radius of the strangelet decreases when the strangeness increases or the charge decreases.

(2), The higher temperature is in favour of the stability of the strangelets. Comparing figure 12 with figure 14, we see clearly that the total area of metastable and stable strangelets expands when the temperature increases.

(3), The higher negative charge and higher strangeness number are in favour of the stability of the strangelets. The stable strangelets are highly negative charged and have high strangeness.

## VI Acknowledgment

This work is supported in part by the NNSF of China under contract Nos. 19975010, 10047005 and 19947001.

## VII Figure Captions

Figure 1. The dots stand for the free energy per baryon $`F/A`$ in a $`SZ`$ plane for all possible strangelets with $`A=5`$ and $`T=20\text{MeV}`$.

Figure 2. The same as figure 1, but for the temperature $`T=50\text{MeV}`$.

Figure 3. The free energy per baryon $`F/A`$ as functions of strangeness number $`S`$ with $`A=5`$ and $`Z=5`$, the two lines represent the temperature $`T=20\text{MeV}`$ and $`50\text{MeV}`$, respectively.

Figure 4. The free energy per baryon $`F/A`$ as a function of strangeness number $`S`$ with $`A=10`$ , $`Z=10`$ and $`T=50\text{MeV}`$.

Figure 5. The free energy per baryon $`F/A`$ as functions of radius $`R`$ with $`A=5`$ , $`Z=1`$ and $`T=50\text{MeV}`$, the three lines represent the strangeness number $`S=2,6`$ and $`9`$, respectively.

Figure 6. The free energy per baryon $`F/A`$ as functions of radius $`R`$ with $`A=5`$ , $`Z=4`$ and $`T=50\text{MeV}`$, the two lines represent the strangeness number $`S=2`$ and $`14`$, respectively.

Figure 7. The free energy per baryon $`F/A`$ as functions of radius $`R`$ with $`A=5`$ , $`S=2`$ and $`T=50\text{MeV}`$, the three lines represent the electric charge $`Z=1,5`$ and $`8`$, respectively.

Figure 8. The free energy per baryon $`F/A`$ as functions of radius $`R`$ with $`A=5`$ , $`S=9`$ and $`T=50\text{MeV}`$, the two lines represent the electric charge $`Z=4`$ and $`1`$, respectively.

Figure 9. The dots stand for the free energy per baryon $`F/A`$ with various strangeness numbers $`S`$ for all possible strangelets with baryon number $`A=5`$ at temperature $`T=50\text{MeV}`$. The masses of nucleon, $`\mathrm{\Lambda },\mathrm{\Xi }`$ and $`\mathrm{\Omega }`$ are represented by the filled squares, respectively.

Figure 10. The same as figure 9, but for baryon number $`A=10`$.

Figure 11. The same as figure 9, but for temperature $`T=20\text{MeV}`$.

Figure 12. The electric charge $`Z`$ as a function of the strangeness number $`S`$ for unstable strangelets (open circles), metastable strangelets (filled circles) and stable strangelets (filled squares) with baryon number $`A=5`$ at temperature $`T=50\text{MeV}`$.

Figure 13. The same as figure 12, but for baryon number $`A=10`$.

Figure 14. The same as figure 12, but for temperature $`T=20\text{MeV}`$.