text:
# Kinematics, turbulence and evolution of planetary nebulae Based on observations taken at the ESO

## 1 Introduction

The study of the expansion velocities of planetary nebulae (PN) has two goals. The first goal is to develop our understanding of the hydrodynamical evolution of the nebula. The importance of this is shown by Stasińska & Tylenda (1994), who compare theoretical evolutionary models with the planetary-nebula population of the Galactic Bulge. They show how the comparison of such models with observed properties of the PN population can be used to derive properties such as the mass distribution of the central stars. However, they comment on the lack of expansion velocity data, which is required to “strongly constrain the average properties of Galactic bulge PN and especially the masses of their central stars”.

The second goal is to characterize the parameters of the superwind and study its dependence on the progenitor’s age and metallicity. It is expected that a dust-driven wind is less efficient at low metallicity, leading to heavier white dwarfs (e.g. Zijlstra 1999), but there is no observational confirmation of this. In fact white dwarf masses measured in (metal-poor) globular clusters (Richer et al. 1997) seem low compared to those of central stars of Bulge planetary nebulae (Gesicki & Zijlstra 2000).

In the present paper we discuss a sample of PNe for which we have high-resolution spectra covering \[O iii\], \[N ii\] and H$`\alpha `$, which allow us to model the velocity fields in some detail. In two cases is the 6560Å He ii line used instead of \[O iii\]. Using these lines which trace both the inner and outer regions, we find that the most common velocity profile is ’U’-shaped, with the highest velocities near the outer edge (as predicted by hydrodynamical models) and the inner edge (which is less expected). The mass-averaged expansion velocity is shown to be a well-defined and measurable parameter. This parameter allows us to accurately place the objects on the HR diagram.

## 2 The observations

The observations were performed at the ESO CAT telescope during 1993 and 1994. This was a subsidiary 1.4m telescope feeding into the CES spectrograph located at the neighbouring 3.6m telescope. The long camera was used giving a spectral resolution of 60 000 (corresponding to 5 km s-1). The slit width was 2″, and the CCD was binned to 2″pixels. The spectra used here use only the central row of pixels. (As the nebulae studied here are compact and the CAT was not an imaging-quality telescope, no spatial information was expected.) The spectrum covers one order of the echelle, including H$`\alpha `$ and the \[N ii\] lines at 6548Å and 6853Å. This was combined with earlier published spectra of the \[O iii\] 5007Å line, take with the same telescope but with a resolution of only 30 000.

We assume the nebulae to be spherically symmetric. This is partly necessitated by the data. The objects tend to be resolved by the instrumental setup (as shown by the line splitting which is often seen). But the CAT suffers from field rotation so that the orientation of the slit on the sky changes during the exposure, and may be different between the two separate setups. If the nebulae are not symmetric, the two spectral settings may trace regions in the nebula with different conditions.

For our modeling we will strictly assume spherical symmetry. Assuming symmetry does not make the nebula symmetric—the most common shape for PNe is in fact elliptical. In a previous paper (Gesicki & Zijlstra 2003) we discuss the effect of ellipticity on the data and models, using very compact (extra-galactic) PNe for which integrated spectra could be obtained. We will comment on this problem below, when discussing the interpretation of the deduced velocity profiles.

## 3 Method of analysis

Recently in Gesicki & Zijlstra (2003) we published a very detailed description of the modelling procedure, therefore we provide only a short summary of the method.

The computer model first solves for the photo-ionization equilibrium, where we attempt to reproduce the observed line intensities. Once a satisfactory solution has been obtained, are the line profiles calculated and compared with the observations.

The nebula is approximated as a spherical shell described by the inner and outer radius, the total mass, the radial density distribution, the radial velocity field, and the chemical composition. As far as possible the values for the radii and chemical composition are taken from the literature. The star is assumed to be a black body with a luminosity and effective temperature. The temperature is again taken from the literature if possible, while the luminosity of the star is determined mainly by the H$`\beta `$ flux, once a distance has been chosen.

Once the photoionization model is complete, we calculate the emissivity distribution for the ionic lines. The line profile can now be calculated, taking into account the slit aperture (which does not cover the full nebula) and the seeing (which scatters light rays outside the aperture into the slit). The velocity field is adjusted to fit the line profiles of all three lines simultaneously. O iii and N ii tend to cover the inner and outer regions of the nebula, respectively, while hydrogen is formed throughout the ionized nebula. This set is therefore the minimum required for a realistic description of the full nebula.

## 4 The nebulae

Spectra of 14 nebulae are analyzed in this paper. Twelve are located towards the Galactic Bulge (but only 7 are likely Bulge members), and two are in the Southern Galactic plane. The objects are specifically chosen as having \[O iii\] observations available. More objects were observed at the CAT in H$`\alpha `$ and \[N ii\], but without \[O iii\] the velocity field can be insufficiently constrained. There are however two objects where a He ii line was detected next to H$`\alpha `$, and these are included in the present paper as this line provides a good alternative to \[O iii\].

We have previously analyzed a large set of nebulae using only the \[O iii\] line (Gesicki & Zijlstra 2000). By necessity a simple, linear velocity field was used. The addition of other lines allows us to significantly improve on the earlier results.

The observed nebulae are listed in Table 1, together with some parameters obtained from the fit. The spectra, together with the fitted line profiles, are shown in Figs. 114.

### 4.1 Comments on individual nebulae

002.6+08.1 (H 1-1)  This object, likely in the Bulge, is the only density bounded object in our sample: the evidence for this comes from the weakness of the \[N ii\] (and \[S ii\]) lines. The electron temperature is high. The inner radius is chosen as 0.3 (as opposed to the usual 0.4) times the outer radius, in order to fit the limited line splitting observed in \[O iii\]. The fitted velocity field increases monotonically with radius, with the velocity minimum at the inner radius. Compared to the other nebulae analyzed here, this appears to be unusual. There is no evidence for a turbulent component in the velocity field, in spite of the wels-type (weak-emission-line star). The jump in the \[N ii\] line profile at $`20`$ km s-1 (Fig. 1) is believed to be an artifact of the data reduction process.

003.9$``$02.3 (M 1-35)  The distance adopted for the fit is inconsistent with Bulge membership: at larger distances the ionized mass becomes rather high (0.36 M). A radio image is shown in Zijlstra et al. (1989), but the resolution is limited. There is some indication that the central peak is slightly shifted. The line profiles show some structure. There is a dip near the centre of H$`\alpha `$ which may be an instrumental problem. The \[O iii\] line show a central peak which is absent in the \[N ii\] line. This feature is not reproduced in the fit, but indicates excess high-excitation gas at or near systemic velocity. The best fit is shown in Fig.2, using a ’U’-shaped velocity field, where the minimum velocity is located at an intermediate radius (possibly near the half-mass radius).

006.1+08.3 (M 1-20)  Schwarz et al. (1992) have published an image but it is difficult to measure a diameter from this. The high-resolution VLA radio image (Aaquist & Kwok 1990) shows a compact torus with very small inner radius (0.25′′) and outer radius of roughly 0.8′′. In the perpendicular direction the outer radius is a little over 1′′. The image appears mildly bipolar. The fit uses an inner radius of $`0.2\times R_{\mathrm{out}}`$, as indicated by the radio image. The \[O iii\] line is very narrow. The best fit therefore proposes almost stationary gas (1 km s-1) at the inner radius. The expansion velocity increases rapidly further out.

278.1$``$05.9 (NGC 2867) This object contains a hot \[WC\] central star. We have no \[O iii\] data but the He ii 6560Å line is detected and provides an alternative. Schwarz et al. (1992) present an image showing a dense inner ring and a very large but much fainter halo. Our fit, and probably our data, only relate to the inner ring. The strong \[N ii\] and \[S ii\] suggest that the inner ring is ionization bounded. This raises the question what ionizes the halo. The optical image suggests the ring is clumpy: photon leakage between the clumps could ionize the halo—or the halo may be recombining. The fit suggests considerable turbulence, as appears to be the rule for \[WC\] stars: e.g. Gesicki & Zijlstra (2003) for a discussion of a similar nebula. In addition to the constant turbulence, the expansion velocity increases outward.

278.8+04.9 (PB 6)  The photo-ionization model works poorly, as the observed He ii 4686Å (Peña et al. 2001) is far stronger than the model predicts. However, the observed line may arise from the star which is a hot \[WC\] star. As in the previous object, the He ii 6561Å was used instead of the \[O iii\] line. A turbulent solution provides the best fit, with a constant but high expansion velocity.

353.3+06.3 (M 2-6)  The \[O iii\] line is very difficult to fit, with a narrow core but wide and strong wings. The \[N ii\] is very asymmetric. We have assumed that this reflects line splitting where one component is very weak, but there may other interpretations. In fact the radio image (Gathier et al. 1983) is almost unresolved, as is the optical H$`\alpha `$ image of Bedding & Zijlstra (1994), suggesting line splitting is less likely. The fit must be considered as uncertain. However, the extended wings of the \[O iii\] line suggests high-velocity gas at or near the inner radius.

353.5$``$04.9 (H 1-36)  The \[O iii\] has very wide wings. The star is listed as a possible symbiotic (Acker et al. 1992). Aaquist & Kwok (1990) find the radio image to be unresolved, even with the VLA (high resolution) A-array. The line profiles of both H$`\alpha `$ and \[O iii\] are asymmetric, with an extended red wing. The compactness of the nebula indicates that this asymmetry is intrinsic to the nebula. A ’U’-shaped velocity profile is used, with velocity minimum close to the outer radius. The outermost density is also reduced to fit the rather narrow \[N ii\] line.

356.2$``$04.4 (Cn 2-1)  The H$`\alpha `$ image of this object is very compact (Bedding & Zijlstra 1994). But our spectra show significant splitting of the \[N ii\] line, suggesting that the low excitation gas is more extended than the slit width of 2′′. The inner radius was chosen as $`0.3\times R_{\mathrm{out}}`$, to avoid line splitting of \[O iii\]. A good fit could be obtained with a ’U’-shaped velocity field.

357.1+03.6 (M 3-7)  The \[O iii\] line is weak and rather noisy which gives relatively poor constraints on the fit. The velocity of the inner region is therefore not well constrained. The \[N ii\] line shows strong splitting whilst \[O iii\] does not, but otherwise the line widths are similar. The \[N ii\] profiles require a component at low velocity, but the \[O iii\] is too wide for this component to solely dominate the inner nebula. We therefore arrive at higher velocities near the inner and the outer edge and lower in-between. The poorly resolved image in Zijlstra et al. (1989) shows some indication for ellipticity.

357.2+07.4 (M 4-3)  The lack of line splitting suggests an unresolved nebula, in agreement with the radio image (Gathier et al. 1983) and H$`\alpha `$ image (Bedding & Zijlstra 1994). The best fit uses a ’U’-shaped velocity field with the highest velocity at the inner edge. The H$`\alpha `$ line is a little too narrow in the adopted fit, which may indicate that the electron temperature is underestimated.

358.2+03.6 (M 3-10)  A ’U’-shaped velocity field is required with very high velocities (60 km s-1) near the inner and outer edge. The high velocity gas is indicated by the wide wings seen in all three lines. The \[O iii\] line shows some splitting: to reproduce this, the inner radius was chosen as $`0.5\times R_{\mathrm{out}}`$.

358.5$``$04.2 (H 1-46)  The radio image shows this to be a very compact PN (Acker et al. 1992). (The optical image appears to suffer from confusion with a nearby star: Bedding & Zijlstra 1994). A ’U’-shaped velocity field was used. However, the narrow \[O iii\] line is not well fitted. The H-line is also a little too narrow in the fit. Attempts to broaden this with a more complicated velocity field, without affecting the other lines, did not work.

359.1$``$01.7 (M 1-29)  This is the most complicated object in the sample. The \[O iii\] line shows several components including wings at extreme velocities and an unresolved component at the systemic velocity. The \[N ii\] line show far less splitting in the observed profile than in the fit. These facts suggests that a bipolar outflow exists in this object. The radio image confirms that this object is bipolar (Zijlstra et al. 1989). Thus, the extreme wings probably result from the polar flow, where the \[O iii\] extends further into the nebula. The component at the systemic velocity could result from a near-stationary torus (but could also be due to the slit having been centered on the brightest component of the bipolar image). A very small inner radius of $`0.1\times R_{\mathrm{out}}`$ was used, with a density maximum at the inner edge. Such a structure is needed because of the lack of line splitting in the H$`\alpha `$ line, while the nebula is quite extended (7′′). None of the lines are well fitted. A spherically symmetric model is probably a poor approximation to this nebula, and we will not consider this object further. However, as the bipolarity was evident from the spectra, it appears likely that the other PNe in our sample do not show such complicated structures, even if good images are lacking.

359.9$``$04.5 (M 2-27)  The image in Gathier et al. (1983) shows it to be almost unresolved. A ’U’-shaped velocity gives a reasonable fit. There are some residuals in the line cores of \[O iii\] and \[N ii\].

### 4.2 Summary

Of the 14 objects studies, ’U’-shaped velocity fields were required for 10. Of the remaining four, two show significant turbulence. Both nebulae contain \[WC\]-type central stars: turbulence around such stars appears to be common (Acker et al. 2002) although the origin of the turbulence is not clear (Mellema 2001). Three nebulae, including one of the turbulent objects, show monotonically increasing velocities. The prevalence of velocity maxima near the inner edge may be surprising. It was first found from the analysis of three PNe for which many emission line profiles were measured, in deep echelle spectra (Gesicki & Zijlstra 2003). The present result suggests inner acceleration may be quite common. In earlier papers we only looked at monotonic velocity fields. The larger residuals in those fits may also be due to the presence of high-velocity gas near the inner edge.

The data show that the 14 PNe of our sample have on average high $`T_{\mathrm{eff}}`$ values, although within the usual temperature range. This is due to the high stellar temperature needed to obtain \[O iii\] emission: at $`T_{\mathrm{eff}}<\mathrm{4\hspace{0.17em}10}^4\mathrm{K}`$, this line is weak and these objects were not observed at the CAT. This selection effect should be taken into account when comparing the present 14 PNe with the larger sample, it is clearly visible in the figures discussed below, e.g. Fig. 19.

## 5 Discussion

### 5.1 Defining the expansion velocity

The earlier analysis using the \[O iii\] line only assumed a linearly increasing velocity field. This often failed to fit the wings of the \[O iii\] line, but a more complicated velocity field would be too poorly constrained. The present models show that most of the nebulae reveal a ’U’-shape velocity field. Whether this field is more common among PNe is hard to prove. Because of the selection effect described above, the U-shaped velocity field may be a consequence of a stronger wind blowing from the hotter stars in our sample.

We compute the mass-averaged expansion velocity $`V_{\mathrm{av}}`$ as described in Gesicki et al. (1998). This single value was introduced to characterize as simple as possible the nebular velocity field. These values are given in Table 1. In the same column in parentheses we show the $`V_{\mathrm{av}}`$ values obtained in Gesicki & Zijlstra (2000) from the \[O iii\] line only with the assumption of a linear velocity field. It is found that the $`V_{\mathrm{av}}`$ do not differ much, except in two cases where previously we fitted the narrow line core only, neglecting the extended wings. We conclude that the mass averaged expansion velocity $`V_{\mathrm{av}}`$ is a well-defined parameter. Although its value is best derived using several spectral lines, a single line can give a reasonably good approximation to $`V_{\mathrm{av}}`$, if the line shapes are not affected by perturbed velocity fields. It appears that the sharp increases of the expansion velocity towards the inner and outer nebular radius affect the line wings shapes but do not alter much the main flow. Very often the fit to the \[N ii\] line wings requires a density decrease simultaneously with the increase of velocity; this reduced the impact of this region on $`V_{\mathrm{av}}`$. The situation is less clear at the inner nebular radius, mostly because it is not well probed by the emission lines in our sample, but the smaller volume in this region contains a relatively small fraction of the nebular mass.

Our new data rather weaken the tendency found in Gesicki & Zijlstra (2000) that there seems to exist a weak correlation between the nebular outer radius and $`V_{\mathrm{av}}`$. This is because the two lowest $`V_{\mathrm{av}}`$ values found in their paper are now revised upwards.

In the appendix we present a compilation of all PN for which expansion velocities have been determined. We note, however, that in most of these only monotonic velocity fields were considered. As argued above, especially the average expansion velocity (weighted by mass) appears to be relatively robust against such model simplifications. In the next subsections the large sample of 73 objects will be discussed. The smaller sample of 14 will be indicated in the figures.

### 5.2 The effect of metallicity

The thermal broadening of hydrogen amounts to roughly 15 km s-1, which is similar to the typical expansion velocities. For this reason the hydrogen lines are not so sensitive to the details of the velocity fields and in previous papers have been relatively easy to fit. But here we find that the models which best reproduce O and N tend to fail at H$`\alpha `$, making the line either slightly too narrow or too broad.

This may suggest that our photoionization model produces incorrect electron temperatures, as this is the parameter which affects hydrogen but to which the other, heavier atoms are much less sensitive.

The electron temperature is determined partly by the stellar temperature but most strongly by the metallicity. The reason is that most of the cooling of the nebula occurs from the forbidden oxygen lines. A lower abundance of oxygen leads to reduced cooling efficiency which is compensated by a higher electron temperature. We ran some test models, and found this to be the strongest effect in changing the width of the hydrogen line.

We did not attempt to perform our own chemical composition analysis.

### 5.3 Ionized masses

One of the results of photoionization modelling is the amount of ionized gas in the PN. This value is adopted to fit the observed and dereddened $`\mathrm{log}F(\mathrm{H}\beta )`$ at the assumed distance, size and density distribution. It can be adjusted simultaneously with the distance which is often poorly known. Sometimes it is possible to build similarly satisfying models with different masses and distances: the selection is made based on additional data. In Fig. 15 we plot the derived ionized masses against the adopted distance. This figure shows that most of our PNe have ionized masses below 0.2 M. The detailed analysis of Gesicki & Zijlstra (2003) of three PNe resulted in higher ionized masses than generally expected. This finding supports our earlier results of adopting higher masses for the Galactic Bulge nebulae. However it seems likely that some of our PNe can be shifted between small distance & low mass and large distance & high mass regions in the plot.

### 5.4 Dynamical ages and core masses

The original AGB outflow velocity $`V_{\mathrm{AGB}}`$ can be estimated if the luminosity and metallicity of the star is known. Habing et al. (1994) show that $`V_{\mathrm{AGB}}L_{}^{0.33}\delta ^{0.5}`$ where $`\delta `$ is the dust-to-gas ratio, assumed to be proportional to metallicity. We have calculated (see Gesicki & Zijlstra 2000) $`V_{\mathrm{AGB}}`$ for the PNe with known metallicity (using the O/H in Table 2). We assumed that a star with solar metallicity (O/H=8.93) has an AGB outflow velocity of 15 km s-1.

Gesicki & Zijlstra (2000) describe a procedure to estimate the dynamical ages of the PNe and masses of their central stars. We follow their recipe for estimating the nebular dynamical age, by calculating the average of the PN expansion velocity and the original AGB velocity, and apply this value to a radius of 0.8 times the outer radius.

The mass-averaged expansion velocities we plot in Fig. 16 versus the dynamical age. We do not include objects for which the metallicity is not known. There is no apparent correlation which suggests that PNe do not speed up with time. Instead the increasing ionization flux together with increasing stellar wind are used for ionizing and accelerating more mass. The plot may suggest that \[WC\]-type PNe are younger than average PNe.

To find the stellar parameters we applied the computer codes to interpolate between the evolutionary tracks of Blöcker (1995) and Schönberner (1983). These codes were written by S.K. Górny and applied in Górny et al. (1997) for the derivation of PN core masses. The basic idea of these programmes is to produce by cubic spline interpolation a very dense grid of evolutionary tracks. We adopted these routines to extract for each stellar temperature, the luminosity and core mass against the evolutionary age. Then we interpolate the core mass and luminosity for each object, based on its dynamical age and stellar temperature. Because the dynamical age depends on the adopted AGB wind velocity, which value is depending on stellar luminosity, we performed a simple iteration. At first we adopted the luminosities found from photoionization modelling, then having located the objects on their evolutionary tracks we recalculated the AGB velocities with new luminosities and again interpolated stellar parameters.

The average core mass of the 14 PNe modelled in this paper is 0.62 M. This is somewhat higher than the mean value obtained by Gesicki & Zijlstra (2000). This difference can be explained again by the previously mentioned selection effect, because more massive stars evolve faster and can reach higher temperatures before the surrounding nebula disperses.

For all PNe analyzed earlier we present the histograms of core masses in Fig. 17. Both the average and the median values are now 0.61 M. Separately we give the histogram for O-type central stars and for those of \[WC\]-type. Both histograms look similar; however, the \[WC\]-type objects reveal a more pronounced peak at the central position. This finding needs verification on a bigger sample.

Although the relative core masses are well determined, including the narrow range in masses, the possibility of systematic effects cannot be excluded. We assume that the acceleration of the nebula occurred uniformly over time. If more time was spent at the original AGB velocity, the dynamical age would be underestimated and the core mass overestimated. We can quantify this by repeating our calculation using the AGB wind velocity only, applied to the measured outer radius. This maximises the age and returns the lowest possible core mass. We find a systematic uncertainty due to this effect of about $`0.015\mathrm{M}_{}`$. This gives a lower limit to the mean core mass of about 0.6 M. The definition of zero position on evolutionary timescale also introduces a systematic effect.

The evolutionary tracks are affected by the poorly known post-AGB mass-loss rates. The tracks assume a mass loss tapering off after the AGB, until replaced by a Pauldrach et al. (1988) wind. If the mass-loss rate becomes similar to the nuclear burning rate ($`10^7\mathrm{M}_{}\mathrm{yr}^1`$), the evolution in the HR diagram will accelerate and the same time scales will be reproduced for lower core masses. There is at present no evidence that mass-loss rates from central stars of planetary nebulae reach such high values. For the early post-AGB phase, the Blöcker (1995) tracks use much higher mass-loss rates than used by Vassiliadis & Wood (1993); the latter predicts much slower evolution which would give even higher core masses.

The peak of the mass distribution of white dwarfs in the solar neighbourhood is approximately 0.59 M (Napiwotzki et al. 1999). Although nominally discrepant, the difference can be due to systematic effects, especially the effect of mass loss on the post-AGB evolutionary tracks. The mass loss in the early post-AGB phase may be higher than assumed by Blöcker (1995), or significant mass loss event may happen during the PN evolution. The ’U’-shaped velocity fields for which we find evidence could be indicative of a strong wind during the PN phase, but it is not evident that this wind is massive enough to affect the speed of the stellar evolution.

### 5.5 Evidence for photon leakage

In Fig. 18 we plot for each object the ratio of the central star luminosity obtained from photoionization modelling over the dynamical luminosity obtained from timescales and evolutionary tracks interpolation. This ratio varies between 2 and 0.1, with the majority below 1. There is some indication that the mean ratio declines with increasing temperature. The values above unity probably reflect the uncertainty of the ionization models, which for instance are quite sensitive to the distance.

Calculating the core masses for the limiting case of no post-AGB acceleration, which gives lower luminosities, does not remove the discrepancy between the luminosities. Taking larger distances can bring the average values in agreement, but results in very large ionized masses. Also, about half our sample is in the Bulge and has well-determined (within 25%) distances.

The stellar luminosity obtained from our modelling is the one needed to ionize just the amount of mass needed to emit like the observed nebula. The fact that the model luminosity is smaller than the value predicted by evolutionary tracks may suggest that the nebula does not intercept and transform all stellar ionizing flux. In other words, the nebulae are leaking. This phenomenon seems to be stronger for the PNe whose cores evolved to higher temperatures. In the sample of 14 nebulae, only a single object is found to be density bounded. In the full sample, few nebulae are density bounded. The other objects are ionization bounded and are not expected to ’leak’ under the assumption of spherical symmetry. The evidence for leaks points to asymmetries, or clumping.

### 5.6 The objects on the HR diagram

The procedure used in calculating the core masses also returns the luminosity of the central star. Together with the central star effective temperatures obtained from the photoionization models, we can plot the PNe in a ’dynamical’ HR diagram. Fig. 19 combines all objects from our earlier publications, the PNe analysed in this paper (framed with squares), and the three PNe of Gesicki & Zijlstra (2003) (framed with circles). Because we have two estimations of stellar luminosities we present the data in two boxes.

The objects show a strong clustering around a single evolutionary track. The fact that the objects agree fairly well with a single track across a range of temperatures suggests that the relative speed of evolution of the tracks is probably correct. If, for instance, the evolution would speed up at a certain temperature, this would appear as a tilted observed sequence compared to the theoretical tracks. Such a tilt is not evident.

We have no objects on the cooling track. The knee in the HR diagram is traversed very quickly (Blöcker 1995) so that few objects would be expected there. Further down on the cooling track the PN become very faint. However, near the knee the procedure to determine the mass and luminosity as a function of temperature and age has two solutions. We select the more luminous one and therefore may have missed one or more objects on the cooling track. The luminosity independently determined in the ionization model provides a check: in only one case (M2-11, a Bulge PN) could this place the object low on the cooling track.

The lower panel of Fig. 19 shows the same diagram but now using the luminosity from the photoionization model. Comparing the two panels shows that the model luminosities do not cluster in the diagram, and are almost always lower than the dynamical luminosities. The lower panel also shows a characteristic decline of luminosity with increasing temperature, in contradiction to the horizontal tracks. This decline has often been seen in nebular luminosities, e.g. Pottasch (1984). Comparison with the dynamical luminosities suggest that measurements of the nebulae (e.g., using the relation $`L_{}=150L_{\mathrm{H}\beta }`$, Gathier & Pottasch 1989, Zijlstra & Pottasch 1989, Magrini et al. 2002), especially with hotter stars, underestimate the stellar luminosity.

### 5.7 \[WC\]-type planetary nebulae

In all figures of section 5 we separated the \[WC\]-type PNe from other types to check how this group of objects behaves. These PNe are characterised by the presence of a hydrogen-deficient central star, which has a much stronger (10–100 times) stellar wind than other PNe cores. Approximately 8% of PNe have \[WC\]-type stars (Górny 2001). The nebulae surrounding \[WC\]-type and O-type stars are, however, very similar. The work of Gesicki & Acker (1996) and Acker et al. (2002) discusses a significant difference between the complete two groups: \[WC\]-type PNe are characterised by strong turbulent motions while the other are not.

Our present work confirms the limited observable differences concerning the turbulent motions. The distribution of stellar temperatures is different, with the \[WC\] stars being on average some 20 000 K cooler–they dominate our sample at the coolest temperatures. However, they were taken from a different sample and this difference can be explained by selection effects. Otherwise, the expansion velocities, ionized masses, the histogram of derived central star masses and the distribution in the HR diagram are the same between the \[WC\] and non-\[WC\] stars. From the presented HR diagram it looks very likely that the \[WC\]-type stars follow the same evolutionary tracks as the other PNe, and are derived from the same progenitor population.

The fact that turbulence is an exclusive property of \[WC\] stars shows that evolution from \[WC\] star to non-\[WC\] star must be rare at best, and that the known \[WC\] stars spend a sufficient amount of time in this phase for their winds to impart sufficient kinetic energy to the nebulae to cause the observed turbulence. The nebular ages are also no different between \[WC\] and non-\[WC\] stars. All these factors suggest that although the two groups arise from the same stellar population, the \[WC\] stars separate relatively early in the post-AGB evolution. The model where the strong wind strips the hydrogen layer of a helium burning star can explain the observed correlation between \[WC\] stars and turbulence, as well as the otherwise very similar nebulae.

## 6 Conclusions

We have presented the expansion velocity fields in a sample of 14 planetary nebulae. The nebulae have high resolution spectra covering for each object three emission lines probing different nebular layers.

Two PNe show significant turbulence and both contain \[WC\]-type central stars. This confirms the existing evidence that turbulence around such stars is common (Acker et al. 2002).

Three nebulae, including one of the turbulent objects, show monotonically increasing velocities.

For ten objects (nonturbulent) we derived a U-shaped velocity field. The outer velocity increase is a well-understood feature of the ionization front. The inner velocity increase was postulated recently by Gesicki & Zijlstra (2003) for two PNe with very hot cores (above 100 000 K). We explained this in terms of a strong wind blowing from a hot central star and sweeping up the inner nebular layers. We thought of it as an exception but the presented results suggests that the inner acceleration may be quite common. It is already present at $`T_{\mathrm{eff}}\mathrm{50\hspace{0.17em}000}\mathrm{K}`$ which is the lowest value in our sample.

We found that the mass-averaged expansion velocity is a reasonably robust parameter. It is well determined even from a single emission line, if this line doesn’t exhibit unusual features. The recently derived U-shape velocity fields result in much improved fits to the line wings, but the $`V_{\mathrm{av}}`$ are almost the same as obtained with simple linearly increasing velocity. We conclude that $`V_{\mathrm{av}}`$ is the proper parameter to describe the main nebular flow.

In the discussion we assembled all previously analysed nebulae. This has not been done before. 24 objects have been discussed by Acker et al. (2002); those PNe were observed in three or more emission lines. 44 objects have been discussed by Gesicki & Zijlstra (2000), observed only in the \[O iii\] 5007Å line. Because we find that the mass-averaged expansion velocity is a robust parameter, we decided to merge the two samples and add three PNe from Gesicki & Zijlstra (2003).

Applying the interpolation routines for evolutionary tracks, we estimate the nebular core masses. The values for our sample of 14 are somewhat larger than the mean value for the full sample. The mean value for the whole sample of 73 PNe is 0.61 M which is higher than the expected average white dwarf mass. This difference can possibly be explained in terms of systematic effects.

For all discussed PNe we also interpolated the predicted by evolutionary calculation stellar luminosities. We find that for the 14 PNe and for the majority of the remaining nebulae, our stellar luminosities are higher than the values obtained from the photo-ionization model of the nebula. We conclude that PNe do not intercept all ionizing radiation even when there is strong evidence of the presence of an ionization boundary.

All analysed PNe show a strong clustering around a single evolutionary track on the HR diagram. This provides support for the evolutionary calculations.

In the discussion we separated the \[WC\]-type PNe from other types. In all plots presented the distributions of the two groups are actually the same. It looks like the \[WC\]-type stars follow the same evolutionary track as the other PNe cores. However, their true position on the HR diagram may be masked by pseudo-photospheres in the strong wind.

## Acknowledgements

The evolutionary tracks interpolation routines were kindly provided by S.K.Górny. This project was financially supported by the CNRS through the JUMELAGE “Astronomie Pologne” programme, by the “Polish State Committee for Scientific Research” through the grant No. 2.P03D.020.17 and by a PPARC visitor grant PPA/G/O/2001/00483.

## Appendix A Full sample of mass-averaged expansion velocities

We present here the data concerning the full sample of PNe with mass-averaged expansion velocities derived. The table contains in consecutive columns: the PNG designation, the usual name, stellar effective temperature, stellar luminosity derived from photoionization models, assumed distance, physical outer radius, model’s ionized mass, mass-averaged expansion velocity, additional turbulent velocity, O/H abundance ratio, nebular dynamical age, interpolated central star mass, remarks.


format: text

text:
# 1 Contour regions of the marginalized likelihood in Ω_𝑚=1-Ω_{𝑡𝑜𝑡}-Ω_Λ versus Ω_Λ plane. Results from CMB data (blue) are superimposed on the contours from SN1a data (orange). Shaded areas correspond to 68% and 95% confidence regions in the Gaussian limit. Degeneracy of the likelihood along the line of the constant angular size of the sound horizon exhibited by CMB data (upper left) is reduced by either imposing a weak prior on the Hubble constant and age (bottom left) or by the Large Scale Structure restrictions (upper right) or both (bottom right).

## Abstract

Full suite of the present day Cosmic Microwave background (CMB) data, when combined with weak prior information on the Hubble constant and the age of the Universe, or the Large-Scale structure, provides strong indication for a non-zero density of the vacuum-like dark energy in our universe. This result independently supports the conclusions from Supernovae Ia (SN1a) data. When the model parameter space is extended to allow for the range of the equation of state parameter $`w_Q`$ for the dynamical field $`Q`$ which gives rise to dark energy, the CMB data is found to give a weak upper bound $`w_Q<0.4`$ at $`95\%`$ CL, however combined with SN1a data it strongly favours $`w_Q<0.8`$, consistent with $`\mathrm{\Lambda }`$-term like behaviour.

Limits on the Dark Energy Parameters from Cosmic Microwave Background experiments

1Physics Department, University of Alberta, Edmonton, AB, T6G 2J1, Canada
2CITA, University of Toronto, Toronto, ON, M5S 3H8, Canada

Measurements of anisotropies in the cosmic microwave background (CMB) radiation now tightly constrain the nature and composition of our universe. High signal-to-noise detections of primordial anisotropies have been made at angular scales ranging from the quadrupole  to as small as several arcminutes . Within the context of models with adiabatic initial perturbations, as are generally predicted by inflation, these measurements have been used in combination with various other cosmological constraints to estimate the values of many important cosmological parameters. Combining their CMB data with weak cosmological constraints such as a very loose prior on the Hubble constant $`H_0`$, various teams have made robust determinations of several parameters, including the total energy density of the universe $`\mathrm{\Omega }_{tot}`$, the density of baryons $`\omega _b=\mathrm{\Omega }_bh^2,h=H_0/100\mathrm{k}m/\mathrm{s}/\mathrm{M}pc`$, and $`n_s`$, the spectral index of primordial perturbations .

In Table 1 we present some of our results for the cosmological parameters together with $`68\%`$ CL errors, obtained from the latest suite of CMB observations. Full details can be found in . Our analysis has included Acbar ,Archeops , Boomerang , CBI , DASI , MAXIMA , VSA  and COBE-DMR  data. Combining these highly consistent datasets allows for a very tight determination of the baryon and dark matter physical densities $`\omega _b`$ and $`\omega _{cdm}`$ and the slope of scalar perturbations $`n_s`$, even with no prior restrictions on the parameters.

| Table 1. Cosmological parameters from combination of CMB experiments | | | | | | |
| --- | --- | --- | --- | --- | --- | --- |
| Priors | $`\mathrm{\Omega }_{tot}`$ | $`n_s`$ | $`\omega _b`$ | $`\omega _{cdm}`$ | $`\mathrm{\Omega }_\mathrm{\Lambda }`$ | $`\tau _c`$ |
| no prior | $`1.2_{0.10}^{0.11}`$ | $`0.96_{0.04}^{0.04}`$ | $`0.022_{0.002}^{0.002}`$ | $`0.11_{0.02}^{0.02}`$ | $`0.24_{0.15}^{0.21}`$ | $`<0.33`$ |
| wk | $`1.06_{0.04}^{0.04}`$ | $`0.95_{0.03}^{0.04}`$ | $`0.022_{0.002}^{0.002}`$ | $`0.12_{0.02}^{0.02}`$ | $`0.56_{0.14}^{0.13}`$ | $`<0.31`$ |
| wk+LSS | $`1.05_{0.03}^{0.03}`$ | $`0.96_{0.03}^{0.03}`$ | $`0.022_{0.002}^{0.002}`$ | $`0.11_{0.02}^{0.02}`$ | $`0.66_{0.09}^{0.09}`$ | $`<0.27`$ |
| wk+SN | $`1.03_{0.03}^{0.03}`$ | $`0.96_{0.03}^{0.03}`$ | $`0.023_{0.002}^{0.002}`$ | $`0.10_{0.01}^{0.02}`$ | $`0.70_{0.05}^{0.04}`$ | $`<0.39`$ |
| wk+LSS+SN | $`1.03_{0.03}^{0.03}`$ | $`0.97_{0.03}^{0.03}`$ | $`0.023_{0.002}^{0.002}`$ | $`0.10_{0.01}^{0.02}`$ | $`0.70_{0.04}^{0.04}`$ | $`<0.37`$ |
| $`\mathrm{\Omega }_{tot}=1`$ plus |  |  |  |  |  |  |
| wk | (1.00) | $`0.96_{0.03}^{0.03}`$ | $`0.022_{0.002}^{0.002}`$ | $`0.12_{0.02}^{0.02}`$ | $`0.72_{0.08}^{0.07}`$ | $`<0.17`$ |
| wk+LSS | (1.00) | $`0.97_{0.03}^{0.03}`$ | $`0.022_{0.002}^{0.002}`$ | $`0.12_{0.01}^{0.01}`$ | $`0.70_{0.04}^{0.05}`$ | $`<0.17`$ |
| wk+SN | (1.00) | $`0.96_{0.03}^{0.03}`$ | $`0.022_{0.002}^{0.002}`$ | $`0.12_{0.01}^{0.01}`$ | $`0.71_{0.04}^{0.05}`$ | $`<0.18`$ |
| wk+LSS+SN | (1.00) | $`0.97_{0.03}^{0.02}`$ | $`0.022_{0.002}^{0.002}`$ | $`0.12_{0.01}^{0.01}`$ | $`0.70_{0.04}^{0.04}`$ | $`<0.18`$ |

Here we shall focus on the parameters describing the vacuum (or dark) energy, which may dominate the present day energy balance. In the minimal parameter set of Table 1 the vacuum energy is assumed to behave as $`\mathrm{\Lambda }`$-term being described just by its density parameter $`\mathrm{\Omega }_\mathrm{\Lambda }`$. Supernovae Ia results  provide a strong indication that $`\mathrm{\Omega }_\mathrm{\Lambda }>0`$.

The non-prior CMB analysis exhibits persistent degeneracy in the $`\mathrm{\Omega }_{tot}\mathrm{\Omega }_\mathrm{\Lambda }`$ space (see Figure 1), which extends along the direction of the constant angular size of the sound horizon at recombination — from closed models $`\mathrm{\Omega }_{tot}1.2`$ with little if any vacuum energy contribution, to flat models $`\mathrm{\Omega }_{tot}1,\mathrm{\Omega }_\mathrm{\Lambda }0.7`$. The degeneracy is reduced and the high likelihood region is concentrated near flat models with significant vacuum energy as soon as even weak prior restrictions (“wk” in Table 1) are imposed on the Hubble constant $`0.45<h<0.9`$ and the age of the universe $`t>10`$ Gyr, or when CMB data is combined with the Large-Scale Structure (LSS) information, namely the amplitude $`\sigma _8`$ and the slope of density power spectrum on galaxy cluster scale (see  for the details of the LSS prior). Thus, present day CMB data provides new independent indication for non zero vacuum density. This conclusion is significantly strengthened by the new data obtained this year, in particular, by the Archeops team which has measured CMB temperature fluctuations in the intermediate angular range $`l20150`$, only weakly constrained by previous data.

We note that the data indicates that the optical depth due to late-time reionization in the Universe is not large, with $`95\%`$ upper limit being in the range $`\tau _c<0.30.4`$. The exact value, however, is not yet robust.

The plausible interpretation of the dark energy is that it is the energy of a “quintessence” scalar field $`Q`$, dynamical in origin, thus, $`\mathrm{\Omega }_\mathrm{\Lambda }=\mathrm{\Omega }_Q`$. A popular phenomenology is to describe additional complexity by one more parameter $`w_Q=P_Q/\rho _Q`$, the ratio of the pressure $`P_Q`$ and the density $`\rho _Q`$ of the field which serves as an effective equation of state. As long as we do not have a pure $`\mathrm{\Lambda }`$ term with $`w_Q=1`$, the $`w_Q`$ is expected to vary in time, but the precision of CMB data does not yet warrant such detailed parametrization. We consider $`w_Q`$ constant and restict our attention to inflation motivated flat $`\mathrm{\Omega }_{tot}=1`$ cosmological models (this prior does not affect other parameters, cf. lower part of Table 1). As another simplification we neglect the effect of late-time perturbations in the $`Q`$ field, which may affect CMB signal at the largest scales. This effect is dependent in detail on the exact model for quintessence potential (one should be cautioned not to extend the equation of state analogy too far, and treat perturbations in the $`Q`$ component hydrodynamically).

Figure 2 shows that CMB alone does not yield useful restriction on the $`w_Q`$ ($`w_Q<0.4`$ at 95% CL) but its combination with SN1a data prefers at $`95\%`$ CL the low values $`w_Q<0.8`$, an improvement to the limit reported in .

Acknowledgements.  We are grateful to Boomerang, CBI and Acbar collaborations with whom much of the work reviewed here has been done. The computational facilities at Toronto are funded by the Canadian Fund for Innovation.


text:
# Chapter 1 STACEE Observations of Mrk 421 During the 2001 Flare

## Chapter 1 STACEE Observations of Mrk 421 During the 2001 Flare

### Abstract

STACEE is a ground based $`\gamma `$-ray observatory that uses a heliostat array, built for solar energy research, to detect atmospheric Cherenkov radiation from $`\gamma `$-ray initiated extensive air showers. During the first half of 2001, a prototype detector, STACEE-48, was used to detect the blazar Markarian 421, which was in an extremely active state. Observations from March to May of 2001 yielded an integral flux of $`(8.0\pm 0.7_{stat}\pm 1.5_{sys})\times 10^{10}\mathrm{cm}^2\mathrm{s}^1`$ at energies above $`140\pm 20`$ GeV, and provide some evidence of correlated trends on time scales of a week or more in the GeV and X-ray bands.

### 1. Introduction

The STACEE–48 detector (Covault et al. 2001, Boone 2002) operated from 2000 to 2001 at the National Solar Thermal Test Facility in New Mexico, USA. This prototype of the full STACEE experiment employed 48 steerable mirrors (heliostats), distributed over approximately $`2\times 10^4`$ square meters. This large collection area allows STACEE to operate into the lower energy $`\gamma `$-ray regime left unexplored by previous ground-based techniques.

Science observations of Markarian 421 were conducted from March through May of 2001, constituting 26 nights and roughly 36 hours of exposure time. Of this, 22 hours of source observations were deemed good-quality data and used for this analysis. Simulations of the detector performance indicate that the peak of the STACEE-48 energy response to a source with the approximate spectral index of Markarian 421 occurs at $`140\pm 20`$ GeV, which has been adopted as our energy threshold for these data (see Boone et al. 2002 for details). Thus, these observations represent the lowest energy $`\gamma `$-ray detection of the 2001 flare by a ground-based instrument.

### 2. Spectral Results

The STACEE observations of Markarian 421 yielded an average $`\gamma `$-ray rate of $`7.7\pm 0.7\pm 1.2`$ min-1. Included in this rate is a correction for false triggers due to the presence of a bright star (HD 95934, magnitude 6.16 in the B band) within a few arcminutes of the source. The correction was determined by observing a star of comparable brightness (HIP 80460) with no known associated $`\gamma `$-ray source in the field of view. The excess rate from the HIP 80460 observations was attributed to a constant star effect, and thus subtracted from the Markarian 421 observations. The above rate, assuming a differential spectral index of $`\alpha =2.1`$, corresponds to an average integral $`\gamma `$-ray flux above 140 GeV of

$$\mathrm{\Phi }(E>140\mathrm{GeV})=(8.0\pm 0.7_{stat}\pm 1.5_{sys})\times 10^{10}\mathrm{cm}^2\mathrm{s}^1.$$

If the spectrum of Markarian 421 is assumed to follow a power-law form in the range of energies to which STACEE is sensitive, the normalization of the differential spectrum can be inferred from the measured rate, provided the spectral index is assumed. In Fig. 1a (Boone et al. 2002), the range of possible power-law solutions to the STACEE rate for differential indices between 2.00 and 2.20 is given by the butterfly-shaped region. For comparison, Whipple observations around the same time are plotted for high (I), medium (V), and low (VIII) flux states (Krennrich et al. 2002). EGRET data from cycles 1 through 5 are also plotted for reference (Hartman et al. 1999).

Since STACEE–48 did not record pulse amplitude information for each event, an unbiased determination of the differential flux is difficult. However, because the energy response of STACEE changes significantly as a function of the source elevation and the trigger configuration, it is possible to glean some information about the differential spectrum by using the forward-folding technique. Under this prescription, the Markarian 421 data set was first divided into elevation bands. The simulated detector response to an assumed power-law spectrum at a given elevation was then compared to the appropriate subset of the Markarian 421 observations for several different reimposed trigger criteria. This process was repeated for a number of different elevation bands. The differential index ($`\alpha `$) and the normalization ($`\varphi _0`$) of the assumed spectrum were then adjusted to minimize the difference between the expected and observed responses over all considered combinations.

Figure 2 includes two tests of this method on simulated data sets with integral fluxes approximately ten times that expected from Markarian 421 (contours represent one sigma errors). Note that the normalization energy ($`E_0`$) has been adjusted to the non-standard value of 500 GeV. This transformation was necessary to decrease the strong correlation in the two parameters. The success of the simulated tests indicates that, in principle, the technique is capable of extracting the spectral parameters. However, due to systematic effects in the simulations themselves, these tests do not vouch for the validity of the technique when used on actual data.

This discrepancy is reflected in the application of the technique to the Markarian 421 data in Fig. 1b (contour is the one sigma error). In units of TeV, the reconstructed normalization translates to $`(6.8\pm 1.7)\times 10^{11}`$ cm-2 s-1 TeV-1. This value, in conjunction with the extremely hard spectral index of 1.6 does not agree with the average STACEE flux in Fig. 1a or the STACEE/Whipple flux ratios described in Boone et al. (2002). Whether these inconsistencies arise from inaccuracies in our detector simulations, the inappropriateness of a power-law postulate, insufficient statistics, or the high variability of the source, remains to be seen.

### 3. Temporal Results

Figure 3a (Boone et al. 2002) contains three light curves for Markarian 421 from March through May of 2001. The upper panel, from the RXTE All Sky Monitor is from the 2 to 20 keV energy band (results provided by the RXTE ASM team). The lower panel, from the Whipple detector (Holder et al. 2001), spans energies from 250 GeV to 8 TeV. The STACEE observations, depicted in the central panel, represent the photon flux between 50 and 300 GeV. Points are photon count rates, averaged on a day by day basis.

A correlation analysis was performed on the STACEE and RXTE data sets using the discrete correlation function (DCF) described in Edelson & Krolik (1988). Unbinned discrete correlation coefficients (UDCF) were formed for each pair of points in the respective data sets as follows:

$$\mathrm{UDCF}_{ij}=\frac{(\mathrm{STACEE}_i\overline{\mathrm{STACEE}})(\mathrm{RXTE}_j\overline{\mathrm{RXTE}})}{\sigma _{\mathrm{STACEE}}\sigma _{\mathrm{RXTE}}}.$$

These coefficients were then grouped by their associated time differences ($`\mathrm{\Delta }t=t_{\mathrm{RXTE}}t_{\mathrm{STACEE}}`$), and averaged to form the DCF.

Figure 3b is a plot of the DCF as a function of time lag for the STACEE and RXTE data sets between MJD 51984 and MJD 52052. Error bars are constructed according to the DCF prescription, and can be understood as the error on the mean for each bin. Although the peak value of the DCF occurs at two days (corresponding to the STACEE data leading the RXTE data by two days), the large error bars on this point indicate that it is not statistically significant with respect to a value of zero (only 2.6$`\sigma `$). The most significant point in the DCF is $`\mathrm{\Delta }t=1`$ (4.0$`\sigma `$). However, as most bins have approximately 20 entries, it is expected that even completely uncorrelated data sets could exhibit DCF values of $`0.2`$ (see Edelson & Krolik 1988 for details). Thus, it is difficult to draw any meaningful conclusions from the DCF alone.

### 4. Conclusions

The STACEE–48 instrument was used to successfully detect Markarian 421 during its flare in the early part of 2001. The average flux observed appears consistent with the variable spectral energy distribution (SED) reported by the Whipple experiment. The data are not suitable for an event-by-event energy analysis, and the forward folding technique yields results which are not consistent with the STACEE average flux. Therefore, it appears that the spectral analysis of these data may be limited to the time-averaged flux reported here and elsewhere. However, the SED shown in Fig. 1a suggests that the STACEE sensitivity may include the peak of the TeV emission for Markarian 421. Thus, more detailed spectral information from STACEE is highly desirable. The current incarnation of the experiment, STACEE–64, is instrumented with waveform digitizers on all channels, which should greatly improve the energy resolution.

A visual inspection of the contemporaneous light curves from RXTE, STACEE, and Whipple suggests correlated trends across all three energy bands. And though a more formal correlation analysis with the DCF does not yet provide compelling statistical evidence for a correlation on day time scales, it does not preclude it either. We suspect that this ambiguity may be due to the rather sparse sampling in the STACEE data set. In fact, it may be that the traditional DCF analysis is generally less suited to GeV-TeV observations, which tend to redefine the notion of “sparsely sampled” data. Alternative timing analysis methods are currently under investigation.

### 5. References

1. Boone, L. M., et al. 2002, ApJ, 579, L5

2. Boone, L. M. 2002, PhD. thesis, Univ. California, Santa Cruz

3. Covault, C. E., et al. 2001, in Proc 27th ICRC (Hamburg), 1, 2810

4. Edelson, R. A. & Krolik, J. H., 1988, ApJ, 333, 646

5. Hartman, R. C., et al. 1999, ApJS, 123, 79

6. Holder, J., et al. 2001, in Proc 27th ICRC (Hamburg), 1, 2613

7. Krennrich, F., et al. 2002, ApJ, 575, L9


text:
# Homogeneity and Size Effects on the Liquid-Gas Coexistence Curve

## Abstract

The effects of (in)homogeneity and size on the phase diagram of Lennard-Jones fluids are investigated. It is shown that standard multifragmentation scenarios (finite equilibrated systems with conserved center of mass position and momentum) are implying a strong radial inhomogeneity of the system strongly affecting the phase diagram. The homogeneity constraint is therefore necessary for finite systems in order to align to the “meaning” of infinite systems phase diagrams. In this respect, a method which deduces the equation of state of homogeneous finite systems from the one corresponding to bulk matter is designed. The resultant phase diagrams show a strong dependence on the system’s size.

Phase transitions have been studied from a long time in the limit of infinite matter. However, there are many physical situations away from the thermodynamic limit at both microscopic and macroscopic scales. Such systems are presently of special interest as so far there is little knowledge about their thermodynamical behavior. There are thus fundamental questions related to these systems such as whether they do present phase transitions and, if yes, how can these be identified. (Recent attempts to extend the standard thermodynamics towards the “small” systems can be found in Refs. gross ; bor ; chomaz ; prl .)

Highly excited nuclear systems are good “laboratories” for thermodynamics studies. Due to the van der Waals type of the nucleon-nucleon interaction these systems are supposed siemens  to exhibit a liquid-gas phase transition as the classical fluids. The connection is however not straight forward due to their (relatively) small number of constituents and the presence of the (long-range) Coulomb force. The effect of the Coulomb interaction on the nuclear liquid-gas phase transitions was studied previously (see e.g. Refs. lee-mek ; prl ). The effects of other features specific to multifragmentation such as finite size (for attempts to address the finitness of the system within Hartree-Fock theories see e.g. Refs. jaqaman ; satpathy ) and degree of homogeneity are still mostly unknown. These aspects are addressed in the present paper.

The classical Lennard-Jones 6-12 (LJ) fluid is used for illustrating the paper’s ideas. While neglecting for the quantum effects acting in nuclear matter, it has the advantages of generality and tractability. The LJ potential writes:

$$v_0(r)=4ϵ\left[\left(\frac{\sigma }{r}\right)^{12}\left(\frac{\sigma }{r}\right)^6\right]$$
(1)
Two versions of the above potential are considered herein. The first is the truncated and shifted (TS) LJ potential: $`v(r)=v_0(r)v_0(r_c)`$ when $`r<r_c`$; $`v(r)=0`$ when $`rr_c`$. The second is the truncated and long range corrected (TLRC) LJ potential: $`v(r)=v_0(r)`$ when $`r<r_c`$; $`v(r)=0`$ when $`rr_c`$, corrections being subsequently included in order to account for the effect of the neglected tail 111 Contributions from the neglected tail to the system’s potential energy and virial energy term \[see eq. (4)\] are taken respectively as:
$`\mathrm{\Delta }v=\frac{1}{2}\rho A_{r_c}^{\mathrm{}}d𝐫v(r)=8\pi \rho Aϵ[\sigma ^{12}/(9r_c^9)\sigma ^6/(3r_c^3)]`$; $`\mathrm{\Delta }𝒱=\frac{1}{2}\rho A_{r_c}^{\mathrm{}}d𝐫(rv(r)/r)=8\pi \rho Aϵ[4\sigma ^{12}/(9r_c^9)2\sigma ^6/(3r_c^3)]`$.. In the present work we use $`r_c=2.5\sigma `$. Phase diagrams of such fluids are subsequently constructed for various situations. Phase transitions in finite systems being here addressed, an adequate representaion is necessary in order to identify them. For example, the $`P(V)|_T`$ curves corresponding to an isochore canonical ensemble will prompt the (first order) phase transition through backbendings even in the case of small systems prl ; npa\_c . Pressure can be easily evaluated starting from its canonical definition: $`P=T\mathrm{ln}Z(\beta ,V)/V`$, where $`Z(\beta ,V)`$ is the system’s canonical partition function:

$$Z(\beta ,V)=\frac{1}{A!\lambda _T^{3A}}\underset{i=1}{\overset{A}{}}d𝐫_i\mathrm{exp}\left(\beta \underset{i<j}{\overset{A}{}}v(r_{ij})\right),$$
(2)
where $`\beta 1/T`$ and $`\lambda _T`$ is the thermal wavelength. The result is the following virial expression:

$$P=\frac{AT}{V}\frac{1}{3V}\underset{i<j}{\overset{A}{}}r_{ij}\frac{v(r_{ij})}{r_{ij}},$$
(3)
where $``$ has the meaning of canonical average. When conservation of the system’s center of mass (c.m.) and c.m. momentum is also considered in $`Z(\beta ,V)`$, one simply has to change $`A`$ into $`A1`$ in eqs. (2) \[the power of $`\lambda _T`$\] and (3) \[the first term\]. While rather redundant for very large systems ($`A\mathrm{}`$), these constraints give important effects when dealing with small systems. Expressions like (3) can be easily evaluated by means of Metropolis simulations. Knowing the formal expression of $`Z(\beta ,V)`$, (2), the statistical simulation of the corresponding canonical ensemble at constant volume is straight-forward: Pick a randomly chosen particle and make a random variation of its initial position (in a fixed interval $`\mathrm{\Delta }V`$). Then, consider this move according to the acceptance $`\mathrm{exp}(\beta \mathrm{\Delta }v_t)`$, where $`\mathrm{\Delta }v_t`$ is the change in the system’s total potential energy. The simulation can be easily adapted to the case of the conserved system’s c.m. and c.m. momentum: Two randomly chosen particles are moved simultaneously with $`\mathrm{\Delta }𝐫`$ and $`\mathrm{\Delta }𝐫`$ such that the system’s c.m. remains fixed ($`\mathrm{\Delta }𝐫`$ being a variation randomly chosen in the volume interval $`\mathrm{\Delta }V`$). Then, the move is considered according to the same acceptance.

The present study is started by considering the case of a system at fixed temperature $`T`$, composed of $`A`$ particles interacting via TS LJ potential, contained into a spherical recipient of volume $`V`$, having the c.m. constrained to coincide with the center of the recipient. The corresponding canonical ensemble is simulated as described earlier. Then, isothermal $`P(V)`$ curves can be evaluated by means of eq. (3), adapted to the considered conservation laws. Subsequently, the borders of the liquid-gas coexistence region are evaluated by performing Maxwell constructions on all $`P(V)|_T`$ curves bending backwards. The resulting phase diagrams (in temperature vs. density representation) corresponding to two systems of different sizes $`A=20`$ and $`A=50`$ are represented in the upper part of Fig. 1. The increase of the critical temperature with the system’s size can be observed. A striking feature of these phase diagrams is the small densities corresponding to the borders of the liquid-gas coexistence regions. In this respect, note that the liquid border is situated at densities smaller than 0.2 $`\sigma ^3`$ which differs a lot from the LJ phase diagrams corresponding to infinite homogeneous systems hansen-verlet  where the liquid border goes up 0.8 $`\sigma ^3`$. The reason for this discrepancy can be easily understood from the lower part of Fig. 1. There, the radial density profiles corresponding to three sample points chosen from the coexistence line of the A=50 system are represented. In all three cases the system appears to be inhomogeneous, its density varying from higher values (towards the center of the recipient) to very small ones (towards the recipient’s walls). The low density tails of these radial profiles are therefore inducing the above mentioned effect on the global density of the system (calculated as $`A/V`$). In particular, note that at small values of the distance from the recipient’s center the densities of the considered sample points tend to be consistent with the ones from Ref. hansen-verlet . Obviously, this inhomogeneity effect is dictated by the c.m. conservation constraint which forces the larger clusters to stay towards the center of the recipient and the smaller ones towards the borders. This example is particularly important in the context of multifragmentation where the system’s c.m. is naturally conserved from event to event and similar density profiles are expected to occur.

Therefore, in order to meaningfully compare the phase diagram calculated for finite systems with the ones corresponding to infinite systems an extra constraint has to be imposed: the homogeneity. The infinite and homogeneous system can be simply approached within the (above-described) Metropolis simulation by implementing a cubic recipient with periodic boundary conditions (PBC). (No c.m. conservation constraint has to be imposed.) A number of 200 particles interacting via TLRC or TS LJ potentials are placed in a cubic box with PBC. The very good agreement between our calculated phase diagram corresponding to the TLRC case (see Fig. 3, the “$`\mathrm{}`$” curve) and the one from Ref. hansen-verlet  is confirming the accuracy of the method here employed. The PBC simulations are used for evaluating the bulk “virial energy” per particle, defined as:

$$𝒱=\frac{1}{3A}\underset{i<j}{\overset{A}{}}r_{ij}\frac{v(r_{ij})}{r_{ij}}$$
(4)
Eq. (3) can be now translated into:

$$P=\rho \left(T𝒱\right),$$
(5)
where $`\rho `$ is the system’s density. Therefore, in an infinite and homogeneous system, the virial energy per unity of volume can be written as: $`\stackrel{~}{𝒱}=\rho 𝒱`$. Given the above definitions one may perform surface corrections for evaluating the $`𝒱`$ term corresponding to finite systems. Since the system is supposed to be homogeneous, its total virial energy can be expressed as the difference between a bulk term and a surface one:

$$𝒱_t=𝒱_b(AA_sf)=\stackrel{~}{𝒱}_b(VV_sf),$$
(6)
where the index $`b`$ specifies that the respective term is a bulk one, $`A_s`$ and $`V_s`$ are respectively the number of particles from the surface and the “volume” of the surface, and finally, $`1f`$ is the ratio between the virial energy of a particle from the surface and the virial energy of a bulk particle. While for very large values of the recipient’s radius ($`R`$) $`1f`$ is rigorously equal to 1/2 (i.e. the particles on the surface have half the number of nearest neighbors they have in bulk), for small values of $`R`$ (the case of small systems) curvature corrections have to be applied to this factor. In the spirit of the previous definition, $`f`$ can be fairly approached by the ratio between the surface of spherical cap situated outside the recipient, corresponding to a sphere of radius $`\sigma `$ having the center on the surface of the recipient and $`4\pi \sigma ^2`$. After some algebraic manipulation one gets:

$$f=\frac{1}{2}\left(1+\frac{\sigma }{2R}\right).$$
(7)
Note that when $`R\mathrm{}`$ then $`f1/2`$. Considering the surface “width” equal to $`\sigma `$ (i.e. no two particles from that region, situated on the same radius, attract each other) one can re-express eq. (6) as:

$$𝒱_t=𝒱_bA\left(1\frac{V_s}{V}f\right),$$
(8)
so that:

$$𝒱=\frac{𝒱_t}{A}=𝒱_b\left\{1\frac{1}{2}\left(1+\frac{\sigma }{2R}\right)\left[1\left(1\frac{\sigma }{R}\right)^3\right]\right\}.$$
(9)
Eq. (9) can be further expressed in terms of $`\rho `$ using the identity $`\sigma /R=[4\pi \rho \sigma ^3/(3A)]^{1/3}`$. Subsequently, eq. (5) \[with $`𝒱`$ given by eq. (9)\] can be applied for calculating the pressure of a finite system of size $`A`$ at various values of $`\rho `$.

An illustration of the results of the method is given in Fig. 2 where $`P(\rho )|_T`$ curves, corresponding to the temperature $`T=1.05ϵ`$ are represented for various sizes of the system. One can observe that the depth of the backbending of the curve is diminishing as the size of the source is decreasing such that at system sizes as small as $`A=20`$ the backbending completely disappears.

Maxwell constructions have been further performed on the $`P(V)|_T`$ curves corresponding to various values of $`T`$ and $`A`$ allowing thus the construction of phase diagrams for various sizes of the system. The results corresponding to the TLRC and TS LJ forms of potential are given in Figs. 3 and 4 respectively. The points from the coexistence region borders obtained via Maxwell constructions can be further interpolated in order to get an estimation of the critical point and a clearer view on the systems’ phase diagrams by means of the Guggenheim scaling relations for the coexistence curve:

$$\rho _\pm =\rho _c\left(1+aϵ\pm bϵ^\beta \right),$$
(10)
where $`\rho _+`$ corresponds to the liquid branch, $`\rho _{}`$ to the vapor branch of the coexistence curve, $`\rho _c`$ is the critical density, $`ϵ(T_cT)/T_c`$ ($`T_c`$ being the critical temperature) and $`\beta `$ is the critical exponent of the coexistence curve. Eq. (10) is fitted on the points obtained via the Maxwell construction method by adjusting the parameters: $`T_c`$, $`\rho _c`$, $`\beta `$, $`a`$ and $`b`$. As observed in Figs. 3 and 4 the quality of the obtained fits is very good. (It is worth mentioning that for the TS potential, infinite system case, $`\beta =0.34`$ which sharply corresponds to the liquid-gas universality class.) For both considered potentials the critical temperature is drastically decreasing and the critical density is increasing with decreasing the size of the source. In particular, note that important deviations from the “$`\mathrm{}`$” phase diagram can be observed even for the phase diagram of a system as large as $`A=10^4`$ (see Fig. 4). This result is particularly important since it proves that thermodynamically the bulk limit is not (even by far) attainable in nuclear multifragmentation experiments (where the equilibrated systems formed are usually smaller then $`A=300`$).

Summarizing, homogeneity and size effects on the system’s phase diagram have been discussed in the framework of the classical LJ fluid. Phase diagrams are constructed using a method adequate for revealing the coexistence region even for small systems: Maxwell constructions on canonical $`P(V)|_T`$ curves. It was shown that standard scenarios like finite systems with conserved c.m., physically consistent with the multifragmentation phenomenon imply a strong radial inhomogeneity in the systems importantly affecting their phase diagrams. Thus, in order to align to the meaning of the phase diagrams corresponding to the infinite fluids ggnhm , an extra constraint has to be introduced for the case of the finite systems: homogeneity. To this aim, a method for constructing the $`P(V)|_T`$ curves corresponding to homogeneous finite systems based on making surface corrections on the bulk virial energy term was designed. This method makes possible the deduction of the $`P(V)|_T`$ curves corresponding to any size of the system, $`A`$, by using the information embedded in the bulk $`P(V)|_T`$ curves. The resulting phase diagrams corresponding to homogeneous systems of various sizes show a strong dependence on the size of the system. In this respect, the critical temperature is drastically decreasing when the size of the system is reduced. For example, important deviations from the bulk phase diagram can be noticed even for systems as large as $`A=10^4`$. This means that the largest equilibrated systems formed in nuclear multifragmentation experiments are not even close to the bulk limit. To this effect one should, of course, couple the inhomogeneity (discussed earlier) and the Coulomb ones. And as shown in Ref. prl  Coulomb is independently bringing a very important contribution towards lowering the system’s critical point.

###### Acknowledgements.

This work was supported by the Alexander von Humboldt Foundation.