text:
# Are dust shell models well-suited to explain interferometric data of late-type stars in the near-infrared?

## 1 Introduction

Supergiant stars and stars on the Asymptotic Giant Branch (AGB) are surrounded by a dust envelope coupled with gas; for example, see Habing (1996) for a review of these objects, Willson (2000) for a review of mass loss, and Tsuji (2000) for evidence of H2O shells around supergiants.

AGB stars are subject to pulsation and other surface phenomena (shock waves, acoustic waves etc., see Lafon & Berruyer (1991)) that levitate matter above the stellar photosphere. Dust grains appear in the envelope at altitudes where the temperature falls below the condensation temperature. The radiation pressure drives the dust grain overflow far from the star up to distances of a thousand stellar radii causing mass-loss.

The mechanism for mass loss in supergiants is not so clear. The main difference compared to Mira stars is that due to the larger luminosity of supergiants the atmosphere is much more extended and has a much greater pressure scale height. Therefore, the mass loss may have a different driving force although pulsation may also play a role.

Both kinds of stars are characterised by their large diameter and strong luminosity and are consequently very good candidates for infrared interferometry. In the next section, we briefly present the instruments which were used to acquire the interferometric data discussed in this paper, namely FLUOR in the K band ($`2.02.4\mu \mathrm{m}`$) and TISIS in the L band ($`3.44.1\mu \mathrm{m}`$). We selected three particular stars of three different types for which we have a substantial set of data: the supergiant $α$~Ori, the semi-regular SW~Vir, and the Mira-type star R~Leo. Furthermore, their envelopes have quite differing optical thicknesses. It appears that for the objects with large optical depth (SW Vir and R Leo) the observed diameters in L are considerably larger than the respective K band diameters (see Sect. 2.3). Measurements in the L band are more sensitive to cooler material above the stellar photosphere than are the K band data. The question is: what is the physical nature of these layers? One possibility might be astronomical dust. Dust reprocesses the star radiation, in particular in the infrared. Depending on the dust abundance, the star’s L band radiation might be superposed by a dust contribution that probes the inner region of the dust envelope at a few stellar radii away from the star. The measured L band radius might consequently be increased compared to the K band, which is mainly fed by stellar emission.

In Sect. 3, the modelling of the dust envelopes around these stars by a radiative transfer code is explained and the parameters of the model are described. In Sect. 4, we compare the model output with the interferometric and photometric data. Finally, in Sect. 5 we discuss the limits of such models and the questions they raise about the L band observations.

## 2 Observations and data reduction

### 2.1 Instruments used

The stars discussed here were observed with the IOTA (Infrared-Optical Telescope Array) interferometer operating (at that time) with two telescopes (Traub 1998). IOTA is located at the Smithsonian Institution’s Whipple Observatory on Mount Hopkins, Arizona. It offers multiple baseline observations providing visibilities at different spatial frequencies. This enables a model-fit to the visibility data to obtain information on the spatial intensity distribution of the sources. The beams were combined with FLUOR (Fiber Linked Unit for Optical Recombination) in the K band (Coudé du Foresto et al. 1998) and with its extension to the L band TISIS (Thermal Infrared Stellar Interferometric Set-up, Mennesson et al. (1999)). Beam combination with FLUOR is achieved by a single-mode fluoride glass triple coupler. The fibres filter the wavefronts corrugated by the atmospheric turbulence. The phase fluctuations are traded against photometric fluctuations which are monitored for each beam to correct for them a posteriori. The accuracy on visibility estimates measured by FLUOR is usually better than 1% for most sources (Perrin et al. 1998) and can be as good as 0.2% (Perrin 2003).

TISIS is the extension of FLUOR to the thermal infrared ($`3.44.1\mu \mathrm{m}`$). A single coupler is used as beam combiner and photometric signals are not monitored. Approximate correction of turbulence-induced flux fluctuations is achieved with the low-frequency part of the interferometric signal. Accuracy in the L band is therefore not as good as in the K band, but still better than without spatial filtering by fibres. The acquisition protocol is also different from that of the K band, since ground-based interferometric observations encounter a new difficulty at $`3.6\mu \mathrm{m}`$. Thermal background emission produces a fluctuating photometric offset which must be correctly subtracted in order to achieve a good calibration. Hence each source observation must be bracketed by a sky background measurement (chopping) that reduces the instrument efficiency.

### 2.2 Data selection

We chose to include three stars in our study which span different spectral types and stages of evolution, namely the supergiant $`\alpha `$ Orionis (Betelgeuse), the semi-regular variable SW Virginis, and the Mira star R Leonis. Some basic parameters of these objects are compiled in Table 1.

These stars are also some of the most observed in our sample. The L band data are extracted from Chagnon et al. (2002) – we have used the data from the February-March 2000 observations. The K band data for $`\alpha `$ Orionis were presented in Perrin et al. (2003a) and observations cover the period February 1996 to March 1997. The K band data for SW Virginis were collected in May 2000 and are presented in Perrin (2003). Lastly, the K band data for R Leonis were published in Perrin et al. (1999) and were collected in April 1996 and March 1997. Except for SW Virginis, there is a large time difference between the dates of the L band and K band data acquisitions. We believe this to not be a major issue - although these stars are known to be variable, the effects we want to analyse are slower than the periodic changes of stars and larger in amplitude.

In order to get the best constraints on the dust shell model spatial distribution, we have included interferometric data acquired with ISI (Infrared Spatial Interferometer, Hale et al. (2000)) at $`11.15\mu \mathrm{m}`$ as published by Danchi et al. (1994).

A journal of all used data points is given in Table 2.

### 2.3 Uniform disk diameters

The visibilities $`V_i`$, and their standard deviations $`\sigma _i`$, in the K and L bands were fitted by a uniform disk model by minimising the quantity:

$$\chi ^2=\underset{i=1}{\overset{n}{}}\frac{\left(V_iM(\varphi _{\mathrm{UD}},S_i)\right)^2}{\sigma _i^2}$$
(1)
where $`n`$ is the number of available measurements in the respective band. $`M`$ is the uniform disk model:

$$M(\varphi _{\mathrm{UD}},S_i)=\left|\frac{2J_1(\pi \varphi _{\mathrm{UD}}S_i)}{\pi \varphi _{\mathrm{UD}}S_i}\right|,$$
(2)
with $`S_i`$ the spatial frequency and $`\varphi _{\mathrm{UD}}`$ the uniform disk diameter. $`J_1`$ is the Bessel function of first order.

The best fit uniform disk radii, and standard deviations, are listed in Table 3.

The last column is the ratio of the L band and K band diameters. It is very close to 1 in the case of Betelgeuse, but the difference is 20–30% for the other two stars.

In the next sections, we investigate the possibility that dust envelopes around the star may be the cause for the diameter variations from K to L.

## 3 Modelling of the dust envelope

The models used in this paper are based on a radiative transfer code maintained at the Observatoire de la Côte d’Azur in Nice, France. A description of the code is presented by Niccolini et al. (2003). We discuss hereafter the main hypotheses and parameters of this model.

### 3.1 Description of the radiative transfer code

The code is based on the modelling of stellar radiation and circumstellar dust envelope interaction by a Monte Carlo method, assuming thermal equilibrium. The star is assumed to be a blackbody whose emitted energy is split in a number of packets, or “photons”. These photons can either pass the envelope directly, or they are scattered once or more in the envelope, or the photons are absorbed and then thermally re-emitted by the dust grains. A grid over the dust envelope is defined with radial zones and angular sections, with the resolution accounting for expected temperature gradients. From the resulting temperature profile, the code derives the spectral energy distribution (SED) for the three individual photon flux contributions mentioned above. Initially, the stellar radius is left as scaling factor. Assigning a specific stellar radius, the total SED can be compared to photometric data in order to constrain the model. Absorption, emission, and scattering by gaseous components are not taken into account.

Furthermore, the spatial flux distribution at a given wavelength can be derived, i.e., a model image of the object at this wavelength can be created. The Fourier transform of the model image can be compared to visibilities measured at that wavelength which yields a further constraint to the model. For the purpose of producing the image, the code allows a zooming factor while the image dimensions are always fixed at $`300\times 300`$ pixels.

### 3.2 Parameters of the envelope

The code offers several parameters that can be chosen. The size distribution $`s(a)`$ of spherical dust particles having the radius $`a`$ was set as $`s(a)a^{3.5}`$, according to the interstellar medium (Mathis et al. 1977) and independent of the distance from the star. The radius $`a`$ was contained in the range $`0.0050.25\mu \mathrm{m}`$.

The particle density distribution $`n(r)`$ of the dust grains that was used in this study includes acceleration effects on dust by radiation pressure and gas drag (see Appendix A). The resulting density distribution is much more peaked at the inner edge of the envelope than is a $`r^2`$ distribution.

For our study, the dust particle material was constrained to astronomical silicates. Their complex dielectric function (respectively complex refractive index) is given by Draine & Lee (1984). It was sampled at 37 wavelengths between $`0.1\mu \mathrm{m}`$ and $`300\mu \mathrm{m}`$, with the sample points considering local extrema.

The previously mentioned parameters were chosen in consistency with the relevant subset of studies reported by Danchi et al. (1994). They were set and fixed in this study. To fit the model to the data, the following input parameters were adjusted. The stellar radius $`R_{}`$ and the stellar effective temperature $`T_{}`$ govern the flux radiated by the star. Spherical symmetry was assumed for a dust shell around the star with an inner radius $`r_\mathrm{i}`$ and an outer radius $`r_\mathrm{o}`$, in units of $`R_{}`$. Furthermore, the envelope is characterised by its optical depth, $`\tau _\lambda `$, which is defined for overall extinction (absorption and scattering). It can be chosen at a given wavelength, $`\lambda `$, expressed in $`\mu `$m.

Since the code, as used for this study, does not take into account possible changes of the dust properties when heated beyond the dust condensation temperature, the temperature map has to be checked accordingly for each set of parameters. As condensation temperature for silicates, we adopted $`T_{\mathrm{cond}}1000\mathrm{K}`$ (Gail & Sedlmayr 1998).

## 4 Comparison between model and observations

We now study the case of three particular stars. The question is whether a simple spherical dust shell model can account for the large diameter differences observed in the K and L bands. The goal is to find the stellar parameters for which visibility points in the K and L bands and at $`11\mu \mathrm{m}`$ are compatible and simultaneously in agreement with the photometric data.

### 4.1 A supergiant: $`\alpha `$ Orionis

Betelgeuse is a Semi-Regular (SRc) variable star of spectral type M1I. It is a late-type supergiant star with a 2335-day period (see Table 1).

The mass-loss rate $`\dot{M}`$ of $`\alpha `$ Orionis is moderate. Knapp et al. (1998) estimate a mass-loss rate from CO emission measurements of $`3.110^7M_{\mathrm{}}\mathrm{yr}^1`$ at HIPPARCOS distance and with an outflow velocity $`v_\mathrm{o}=14.2\mathrm{km}\mathrm{s}^1`$. Throughout this study, we adopt terminal outflow velocities from their work. For a distance $`d`$ of $`150\mathrm{p}`$ and an outflow velocity $`v_\mathrm{o}=15\mathrm{km}\mathrm{s}^1`$, Danchi et al. (1994) yield a mass-loss rate of $`1710^7M_{\mathrm{}}\mathrm{yr}^1`$. With the HIPPARCOS parallax of $`7.63\mathrm{mas}`$ yielding a distance of $`131\mathrm{pc}`$, and $`\dot{M}v_\mathrm{o}d`$, their mass-loss rate is updated to $`1410^7M_{\mathrm{}}\mathrm{yr}^1`$, which is given in Table 4.

#### Observational data

We used low spectral resolution data from different catalogues available at the SIMBAD (2003) Astronomical Database to constrain the SED over a large range of wavelengths. The TD1 catalogue provides fluxes around $`0.2\mu \mathrm{m}`$. UBV measurements are reported in the GEN and UBV catalogues. Fluxes in the UBVRIJHKL bands are taken from the JP11 Johnson catalogue. The 12, 25, 60, and $`100\mu \mathrm{m}`$ fluxes are from the IRAS catalogue, and Danchi et al. (1994) provide flux curves from 8 to $`22\mu \mathrm{m}`$.

For Betelgeuse, having slight luminosity variations, we can easily compare visibility points from different epochs. This assumption will be fully justified by the fit we obtain.

#### Parameters

Danchi et al. (1994) propose a spherical dust shell model to account for ISI data in the thermal infrared. This shell is supposed to have a small extent ($`r_\mathrm{i}=45.9R_{}`$, $`r_\mathrm{o}=48.2R_{}`$) quite far away from the central star and to be optically very thin ($`\tau _{11}=0.0065`$). It seems to indicate an episodic mass-loss process, leading to an empty area between the star and its envelope. We used these parameters (see Table 5) to reproduce the shorter wavelength visibility points given by FLUOR and TISIS.

The effective stellar temperature of $`3640\mathrm{K}`$ used by Danchi et al. (1994) is fully consistent with the evaluation by Perrin et al. (2003a) giving $`3640\mathrm{K}`$ and $`3690\mathrm{K}`$ depending on models. The stellar radius used in that model is $`21.8\mathrm{mas}`$, satisfactorily accounting for both the photometric data and interferometric data in the three bands. This radius is close to the uniform disk model radii in the K and L bands of $`21.6\mathrm{mas}`$ and $`21.4\mathrm{mas}`$, respectively.

#### Results

The SED of the model is shown in Fig. 1. Apart from the total flux from the object, the plot also contains the contributions of the stellar flux attenuated by the circumstellar shell, the distribution of the flux scattered by the dust grains in the shell (maximum at short wavelength), and the thermal emission of the grains (maximum at longer wavelength). As the SED shows, stellar radiation is prominent over a large wavelength range. Flux emitted by dust becomes prominent only in a small band around $`10\mu \mathrm{m}`$.

Visibility curves are also plotted in Fig. 1. The K and L band visibility curves from the dust shell model resemble each other and show good superposition with the respective visibility curves of the uniform disk model. Apparently, the same structure is seen in these two bands. The SED plot indicates that this should be the stellar photosphere since in these two bands the direct stellar flux is roughly $`10^3`$ times larger than the scattered part and about $`10^5`$ times larger than the dust emission. One can conclude that the dust envelope, having a very small optical depth, has little influence on the intensity distribution at these wavelengths. The dust shell is transparent in K and L and the star is seen through the dust at these wavelengths. As a consequence, the dust shell model derived from the ISI data is fully consistent with the FLUOR and TISIS measurements.

As a further consistency check, yet with minor priority, we computed a mass-loss rate $`\dot{M}`$ from our model parameters by a method similar to Danchi et al. (1994) (see Appendix A for details). With the set of parameters for $`\alpha `$ Ori, we obtain a mass-loss rate of $`1710^7M_{\mathrm{}}\mathrm{yr}^1`$, consistent with Danchi et al. (1994) and roughly a factor of 5 higher than that of Knapp et al. (1998).

Although $`\alpha `$ Ori has the lowest optical depth of the studied objects, its mass-loss rate seems to be the highest. The parameters used here indicate an episodic event in the past. Yet, since this event the dust shell has expanded, with its enlarged inner radius entering the computation. For a further discussion of the mass-loss event see Danchi et al. (1994).

We now study the case of stars with circumstellar envelopes at higher optical depths.

### 4.2 A semi-regular variable: SW Virginis

The pulsation period of this M7III semi-regular (SRb) variable star is 150 days. Its main characteristics are recapitulated in Table 1. The mass-loss rate is estimated by Knapp et al. (1998) to $`1.710^7M_{\mathrm{}}\mathrm{yr}^1`$.

#### Observational data

Photometric data are retrieved from the IRC catalogue for the I and K bands, from the UBV catalogue in the respective filters, from the IRAS catalogue at 12, 60, $`100\mu \mathrm{m}`$, and from Monnier et al. (1998) in the 8 to $`13\mu \mathrm{m}`$ range. The M band flux at $`5\mu \mathrm{m}`$ comes from Wannier & Sahai (1986).

To date, no 11-$`\mu \mathrm{m}`$ interferometric data is available for SW Virginis. Benson et al. (1989) report observations by speckle interferometry in the N band at spatial frequencies below $`210^5\mathrm{rad}^1`$. Since they do not provide measured data directly, but rather give some model parameters, this has not been considered here. It is therefore not possible to assess the consistency of the spatial model between all bands for this object.

#### Parameters

For SW Vir, Perrin et al. (1998) give an effective temperature of $`2921\pm 100\mathrm{K}`$, as derived from the bolometric flux and the limb-darkened diameter. Consistent with this value, the effective temperature was set to $`2800\mathrm{K}`$. For the parameters of the dust shell we chose to adopt $`r_\mathrm{i}=15.0R_{}`$ and $`r_\mathrm{o}=1000R_{}`$, as proposed by van der Veen et al. (1995) on the basis of sub-mm observations. Applying an optical depth $`\tau _{11}=0.045`$ and assuming a stellar radius of $`8.65\mathrm{mas}`$ then accounts for flux data at low spectral resolution.

Danchi et al. (1994) make a distinction between two classes of objects. The first class has a cold circumstellar envelope far from the star (beyond 10 stellar radii), while objects of the second class are surrounded by a warmer dust envelope, near the star (at about 2 or 3 stellar radii). The mass-loss process is different for these two classes. For the second class, the mass-loss process is more continuous than for the first one. SW Virginis has been classified in the past as a Mira star, hence with a somewhat regular mass-loss. Yet, its current classification as a semi-regular variable, possibly a precursor of a Mira star, makes its mass-loss more likely to be irregular, consistently with its irregular photometric variations. It is therefore very tempting to identify it as a member of the first class. The preferred set of parameters obtained for this star is summarised in Table 5. A large optical depth and a large shell extent are necessary to reproduce the photometric contribution of dust at large wavelengths.

#### Results

As Fig. 2 shows, the contribution of dust emission to the SED in the K and L bands of SW Virginis is more important than in the case of $`\alpha `$ Ori, i.e., it is higher in relation to the direct flux. This can be attributed to the larger optical depth. SW Virginis has some infrared excess that requires a dust shell. However, the dust cannot be too warm, otherwise it would produce a much larger infrared excess than is compatible with the SED.

From this model, visibility curves were derived. Although comparatively few visibility measurements are available for this object, it is obvious that the single dust shell model is not simultaneously compatible with the K and L band data. Again, model visibility curves are similar in K and L which lets us assume that we see the same structure, i.e. the stellar photosphere. This is also reflected by the SED for which a similar reasoning applies like in the case of $`\alpha `$ Ori. Indeed, the dust being quite far away from the photosphere, the K band data are well reproduced by the visibility model. However, the dust is too cold and optically too thin to lead to a significant increase of the diameter at larger wavelengths compatible with the L band data. A simple dust shell model, therefore, cannot account for both the interferometric and the photometric data. Unfortunately, no ISI data are available for this star. They would have been very useful to help model the dust shell and place stronger constraints on it.

The mass-loss rate we obtain from our parameters is $`1.410^7M_{\mathrm{}}\mathrm{yr}^1`$, which is roughly consistent with the one of Knapp et al. (1998).

### 4.3 A Mira-type variable: R Leonis

R Leonis is a Mira-type long-period variable with a period of $`310`$ days and of late spectral type M8IIIe (see Table 1). Its mass-loss rate is indicated by Knapp et al. (1998) as $`0.9410^7M_{\mathrm{}}\mathrm{yr}^1`$, whereas the estimate by Danchi et al. (1994) can be updated to $`0.6010^7M_{\mathrm{}}\mathrm{yr}^1`$.

R Leonis is known to have a very strong infrared excess and very steady pulsations which makes it very different from SW Virginis. R Leonis shows large photometric variations, especially at short wavelengths. The problem here is to synthesise a low spectral resolution SED, as well as visibility curves, for a star with large photometric variations. Nevertheless, we decided to include all available data points (both photometric and interferometric) into the plots despite the object variability which will allow a larger range of acceptable model parameters.

#### Observational data

For R Leonis, UBV measurements can be found in the GEN and UBV catalogues. The JP11 Johnson catalogue provides magnitudes in the BVRIJKL bands. Strecker et al. (1978) observed this object between 1.2 and $`4.0\mu \mathrm{m}`$. The fluxes at 12, 20, 60 and $`100\mu \mathrm{m}`$ are from IRAS, and fluxes between 8 and $`20\mu \mathrm{m}`$ come from Danchi et al. (1994).

#### Parameters

The effective temperature that accounts for the photometric data is about $`2500\mathrm{K}`$, a bit more than the temperature of $`2000\mathrm{K}`$ assumed by Danchi et al. (1994) for their model containing only silicate dust around R Leo. The stellar radius is clearly smaller ($`14.9\mathrm{mas}`$ instead of $`19.8\mathrm{mas}`$). However, we have estimated these two parameters with all data combined (interferometric data in the K and N bands and the spectro-photometric data), and they show good consistency (in as much as the star is neither a perfect blackbody nor a uniform disk) with those derived from the K band interferometric data alone by Perrin et al. (1999), who give uniform disk radii of at most $`15\mathrm{mas}`$ and photosphere effective temperatures closer to $`3000\mathrm{K}`$ rather than $`2000\mathrm{K}`$. Our effective temperature is larger than that of Danchi et al. which is not surprising since in their case only photometric data above $`8\mu \mathrm{m}`$ was used, i.e., in a range of wavelengths very sensitive to cool dust emission. This also explains why their stellar radius is larger: for a given bolometric flux, the stellar diameter increases with decreasing temperature.

As in Danchi et al. (1994), we consider a large circumstellar envelope near the central star. But instead of 2 $`R_{}`$ for the inner edge of the dust shell, we chose 3.5 $`R_{}`$ to keep the temperature of the grains close to $`1000\mathrm{K}`$, the assumed temperature of grain condensation. The overall parameters are summarised in Table 5.

#### Results

The fit of the photometric data described in the previous section is displayed in Fig. 3. The contribution of the envelope dominates beyond $`10\mu \mathrm{m}`$ up to about $`30\mu \mathrm{m}`$. With a quite large extent, up to 140 $`R_{}`$, and an optical depth $`\tau _{11}=0.1`$ the dust shell accounts for IRAS and Danchi et al. (1994) observations in the long wavelength range. As can be noticed in the graph, scattered flux and thermal dust emission are even more important than for SW Vir. This is apparently due to the larger optical depth and the smaller inner radius of the envelope. These parameters lead to a mass-loss rate of $`0.7710^7M_{\mathrm{}}\mathrm{yr}^1`$, which is of the same order of and intermediate between the values published by Danchi et al. (1994) and Knapp et al. (1998) (see Table 4).

Visibility curves gained from this model are also shown in Fig. 3 for the K, L, and N bands respectively. The stellar radius derived from our model is $`14.9\mathrm{mas}`$, much smaller than the $`19.8\mathrm{mas}`$ of Danchi et al. (1994). Since our study is based on interferometric data spanning a wider range of wavelengths, the parameters we derive are different from that of Danchi et al. (1994) The ISI data are better suited to characterise the dust, yet the shorter-wavelength data bring more accurate stellar parameters, hence explaining the discrepancies between the two studies. The dust envelope structure may be considered as stationary with little sensitivity to stellar pulsation. FLUOR observations at phase 0.28 (March 1997) in K band give a uniform disk radius of $`15.34\mathrm{mas}`$ (Perrin et al. 1999). In principle, the K band measurements are more sensitive to regions very close to the photosphere. It therefore seems that the radius of R Leonis should be smaller than the evaluations by Danchi et al. (1994) and the difference between our estimate and that of Danchi et al. cannot be attributed to the star pulsation. This difference must be attributed to the difference of effective temperatures.

We now compare interferometric data in the K and L band and that of ISI. Our parameters are still consistent with the ISI visibility points at $`11\mu \mathrm{m}`$ and account for FLUOR observations in the K band at low spatial frequencies. It is to be noticed that, since the ISI visibility data are at low spatial frequencies, they are not as sensitive to the different models as are the visibilities at shorter wavelengths. This explains why they are consistent with both the Danchi et al. (1994) and our model. On the other hand, in the L band our model overestimates the measured visibility values. As in the case of SW Vir, it appears that the large structure seen in the L band is not well accounted for by a simple dust shell model, despite the even larger optical depth. Also the SED shows that contributions from the dust shell are at least a factor $`10^2`$ below the direct stellar flux. A larger radius could explain both $`3.6\mu \mathrm{m}`$ and $`11\mu \mathrm{m}`$ observations, but would disagree with K band measurements and the SED.

## 5 Discussion

### 5.1 Obtained results

From what we have seen so far, the best consistency between the multiple band data and the dust shell model is obtained for $`\alpha `$ Ori. The model described by Danchi et al. (1994) has a stellar radius close to the K and L radii estimated with a uniform disk model. Yet, in this case the influence of the envelope is not very important because of its small thickness and low opacity. Nevertheless, model parameters for ISI data from Danchi et al. are consistent with what we see both for K and L bands.

One may note, however, that Weiner et al. (2000) report further ISI data on $`\alpha `$ Ori which are not fully consistent with the model of Danchi et al. based on observations at low spatial frequencies. Weiner et al. took their measurements at higher spatial frequencies where the dust shell is fully resolved. Its contribution to the visibility is therefore negligible and the observed data are governed by radiation coming directly from the star or from its close environment. Like the L band data of SW Vir and R Leo (see below), the data of Weiner et al. are overestimated by the dust shell model. These findings are further discussed in Perrin et al. (2003a).

For AGB stars with later types like SW Vir and R Leo, for which we measure large diameter differences between $`2.15\mu \mathrm{m}`$ and $`3.77\mu \mathrm{m}`$ (the approximate effective wavelengths of the FLUOR and TISIS instruments), a simple dust shell model cannot account for these differences. However, the shape of near-infrared visibility curves for SW Vir and R Leo seem to show a trend. L band visibility points are systematically overestimated by a simple dust shell model meaning the source is larger than predicted by the model. This shows that the stellar environment mainly seen at this wavelength is much more complex than suggested by a model that takes into account only the star and its dust envelope.

Indeed, the temperature of this area is close to the temperature at which grains form. Consequently, we must also take the presence of gas into account in modelling. Doing this will radically change the radiation transfer behaviour. Recently, spectroscopic detection of SiO, OH, and H2O masers in late-type star envelopes have revealed the necessity of modelling the gas to explain large diameters seen in L band. ISO/SWS and other studies (Tsuji et al. (1997); Yamamura et al. (1999); Tsuji (2000) brought up the existence of molecular layers of warm and cold H2O, SiO, and CO2 around AGB stars. Such layers in AGB stars certainly exist and have an impact on the stellar radiation which is observed by infrared interferometers. Strong molecular scattering by H2O and CO were suggested to account for K band visibility measurements of R Leo by FLUOR in 1997 (Perrin et al. 1999) and seem to explain more easily the near-infrared measurements (Mennesson et al. 2002; Perrin et al. 2003b). It is very likely that dust needs to be coupled with gas to produce a more reliable model of these sources and that dust alone cannot be the only cause of the measurements taken in a wide range of wavelengths. The inclusion of a detailed photospheric model with line opacity due to gas would improve the agreement of the SED, particularly in the blue, and would improve the discrimination between stellar and shell components.

### 5.2 Data used

Certain aspects of our study can be further weighed. Simultaneous interferometric data would certainly help constrain models more accurately. Same epoch measurements are needed because of the large diameter variations of AGB stars. At least for the position angles it can be noted, though, that they are fairly close in range, which leaves little possibility that the discrepancy we see may be explained by an elliptical shape (see below). Furthermore, 11.15-$`\mu \mathrm{m}`$ data would have been very useful to see whether our SW Vir dust model parameters are correct since the envelope has more impact on the model at this wavelength. Such data might be provided in the near future by single-dish instruments on 10-m class telescopes or by upcoming mid-infrared interferometers like MIDI at the VLTI.

### 5.3 Model code parameters

A limitation of the model as used for this study is that it does not consider heterogeneous grain size distributions $`s(a)`$ depending upon the distance from the star or temporal effects (like grain formation, growth, or destruction). Also, different limits of the grain radius $`a`$ might have affected the fit of the SED in the UV and the far-infrared. Yet, varying the grain size range from $`0.0050.25\mu \mathrm{m}`$ to $`0.0052.5\mu \mathrm{m}`$ and $`0.052.5\mu \mathrm{m}`$ for R Leo as a study case, the resulting diagrams showed only small changes. It appeared that, in the modelled visibility curves for the K and L bands, the contribution of spatial frequencies below $`2010^5\mathrm{rad}^1`$ were somewhat increased. Above all, neither the star nor the envelope are ideal blackbodies which adds to these deviations.

Concerning the particle density distribution $`n(r)`$, Suh (1999) questions the necessity of assuming a distribution which is peaked at the inner edge. Indeed, comparing results of calculations with two sets of parameters that differ only in $`n(r)`$ showed only minor changes in K and L bands. For $`11.15\mu \mathrm{m}`$, changes were somewhat clearer, whereby for $`\alpha `$ Ori and R Leo the accelerated version fitted ISI data better. Other issues are the possible existence of smaller structures in the envelope like clumps or aspherical dust distributions, as they were investigated, for example, by Lopez et al. (1997). These have not been tested here and require more intense visibility modelling. This would be definitely required if sources showed systematic asymmetric effects, yet which are likely to produce ringing features on top of a visibility curve compared to the more basic effect of radiative transfer in spherical shells. Nevertheless, closure phase observations could help test for asymmetries.

As for the dust particle material, slightly differing optical properties were proposed by several authors. See, for example, Suh (1999) for a comparison of optical properties of silicate dust grains in the envelopes around AGB stars. We chose the ones of Draine & Lee (1984) which are widely used. Using optical properties suggested by Suh (1999) for comparison made it necessary to choose different dust shell parameters in order to achieve reasonable consistency with ISI data. Yet, L band data were still overestimated, similar to Fig. 3. Another issue might be the dust composition. For example, for R Leo, Danchi et al. (1994) tested also a mixture of silicates and graphites. Yet, this led to a required stellar radius in the dust shell model almost twice as large as observed in the K band, which is supposedly sensitive to the region close to the stellar photosphere. Therefore, this has not been further investigated here. A third aspect might be the condensation temperature $`T_{\mathrm{cond}}`$ that we adopted. Some authors chose higher $`T_{\mathrm{cond}}`$, see, for example, Willson (2000) for an overview. Taking R Leo as study case, we set $`r_\mathrm{i}=2R_{}`$ and chose $`r_\mathrm{o}`$ and an optical depth $`\tau _{11}`$ so that also ISI data were well fitted. The output of the model showed a somewhat increased IR excess around $`11\mu \mathrm{m}`$, and the temperature at the inner radius turned out to be about $`1400\mathrm{K}`$. Still, while the SED, K band, and ISI data were reasonably well fitted, the modelled visibility curve in the L band was clearly higher than the measured data, similar to Fig. 3.

In summary, we do not believe the aspects mentioned above to affect our conclusions significantly.

## 6 Conclusions

High angular resolution techniques give important constraints on spatial intensity distribution of late-type stars. Using diameter measurements in the K and L bands, we have tried to explain both K and L observations for three stars surrounded by circumstellar envelopes with different optical depths. We applied a model with a single dust shell and consistent with low resolution photometric data as well as thermal infrared interferometry data from ISI when available.

For the supergiant Betelgeuse, the optically very thin dust shell has little influence on what is seen and is transparent in the K and L bands. The observed diameters are almost equal. In contrast, the AGB stars SW Vir and R Leo show bigger diameters in the L band than in the K band. Here, the optical depths of the shells are larger and near-infrared model visibilities react quite sensitively on changes of the dust shell parameters. Yet, the models we find in agreement with the SED, K band interferometry, and the specially dust-sensitive $`11\mu \mathrm{m}`$ interferometry cannot reproduce the observed diameter differences.

Our study leads to the conclusion of the insufficiency of a simple shell model with dust only or without more complex structures. New observations at several wavelengths from visible to far-infrared would help understand better external envelopes of AGB stars. A denser coverage of the spatial frequency plane will also be very useful to strongly constrain models. Interferometric observations with higher spectral resolution could enable to localise molecular components of these envelopes and certainly help explain differences between K and L band.

## Appendix A Evaluation of the mass-loss rate

Let $`\dot{M}_{\mathrm{dust}}`$ be the constant mass-loss rate of dust ejected from the star, $`m_\mathrm{d}`$ the average single dust grain mass, and $`n(r)`$ the dust particle density, then at the inner edge $`r_\mathrm{i}`$ of the dust shell it holds that

$$\dot{M}_{\mathrm{dust}}=\frac{4\pi r_\mathrm{i}^2\mathrm{d}rm_\mathrm{d}n(r_\mathrm{i})}{\mathrm{d}t}=4\pi r_\mathrm{i}^2v_\mathrm{i}m_\mathrm{d}n_\mathrm{i},$$
(3)
where $`v_\mathrm{i}`$ is the dust outflow velocity and $`n_\mathrm{i}=n(r_\mathrm{i})`$ the particle density.

In the dust shell, the particle density of grains is assumed as (Schutte & Tielens 1989; Danchi et al. 1994)

$$n(r)=n_\mathrm{i}f(r)=n_\mathrm{i}\frac{r_\mathrm{i}^2}{r^2}\frac{1}{(1Ar_\mathrm{i}/r)^{1/2}}$$
(4)
where the parameter

$$A=1(v_\mathrm{i}/v_\mathrm{o})^2$$
(5)
can be adapted to different cases. It is zero, if $`v_\mathrm{i}`$ is assumed to equal the terminal outflow velocity $`v_\mathrm{o}`$ of the dust far away from the star. In this case, the particle density simply follows $`n(r)r^2`$. For taking into account acceleration of the dust by radiation pressure and coupling with the gas, $`v_\mathrm{o}`$ will be higher than $`v_\mathrm{i}`$. The parameter $`A`$ will then be larger than in the first case and the particle density much more pronounced at $`r_\mathrm{i}`$. In the present study we chose $`A=0.99`$, which applies for typical conditions (Danchi et al. 1994).

The optical depth $`\tau _\lambda `$ of the envelope is determined by

$`\tau _\lambda `$ $`=`$ $`{\displaystyle _{r_\mathrm{i}}^{r_\mathrm{o}}}C_{\mathrm{ext}}(\lambda )n(r)dr`$ (6)
$`=`$ $`C_{\mathrm{ext}}(\lambda )n_\mathrm{i}{\displaystyle _{r_\mathrm{i}}^{r_\mathrm{o}}}f(r)dr,`$
where $`n(r)`$ was replaced by (4). $`C_{\mathrm{ext}}`$ is the average extinction mean section of the dust grains. Replacing $`n_\mathrm{i}`$ in (3) by expression (6) yields a relation for the mass-loss rate by dust:

$$\dot{M}_{\mathrm{dust}}=\frac{4\pi r_\mathrm{i}^2v_\mathrm{i}m_\mathrm{d}\tau _\lambda }{C_{\mathrm{ext}}(\lambda )_{r_\mathrm{i}}^{r_\mathrm{o}}f(r)dr}.$$
(7)
In the previous two expressions, $`r_\mathrm{i}`$, $`r_\mathrm{o}`$, and $`\tau _\lambda `$ are the chosen model parameters. $`C_{\mathrm{ext}}`$ is calculated in the modelling code for a given grain size distribution from the optical properties of the dust by Mie theory (Niccolini et al. 2003). Likewise, the integral in (7) is performed by the code. The dust grain mass $`m_\mathrm{d}`$ has to be estimated by averaging the size distribution $`s(a)`$ and defining the mass density. For the presently considered silicate dust, we chose as $`\rho _\mathrm{d}=3.3\mathrm{g}\mathrm{cm}^3`$ (also used by Draine & Lee (1984) and Danchi et al. (1994)). The terminal velocity of the shell particles is taken from high-spectral resolution observations of CO lines (Knapp et al. 1998).

Furthermore, for the total mass-loss rate also the gas has to be taken into account, i.e.,

$$\dot{M}_{\mathrm{total}}=\dot{M}_{\mathrm{dust}}+\dot{M}_{\mathrm{gas}}.$$
(8)
Gas-to-dust mass ratios are not well known and differ quite a lot. We adopted here, as typical value and assuming dust and gas to have the same terminal velocities (Schutte & Tielens 1989; Danchi et al. 1994; Habing 1996),

$$\dot{M}_{\mathrm{gas}}250\dot{M}_{\mathrm{dust}}.$$
(9)
This means that most of the mass-loss is carried by gas. For $`\dot{M}_{\mathrm{gas}}`$, an expression analogous to (3) can be given where $`v_\mathrm{o}`$ enters for $`v_\mathrm{i}`$. Considering (9), and therefore neglecting the mass-loss carried by dust, the total mass-loss rate is eventually

$$\dot{M}_{\mathrm{total}}=250\frac{4\pi r_\mathrm{i}^2v_\mathrm{o}m_\mathrm{d}\tau _\lambda }{C_{\mathrm{ext}}(\lambda )_{r_\mathrm{i}}^{r_\mathrm{o}}f(r)dr}.$$
(10)
###### Acknowledgements.

P. Sch. appreciated support by the European Community through a Marie Curie Training Fellowship for an extended stay at Observatoire Paris-Meudon. It also allowed the collaboration with the Observatoire de la Côte d’Azur. – The authors are solely responsible for the information communicated. This publication made use of NASA’s Astrophysics Data System (ADS) Bibliographic Services (http://cdsads.u-strasbg.fr/) and of the SIMBAD database, operated at the Centre de Donné astronomiques de Strasbourg (CDS), Strasbourg, France (http://simbad.u-strasbg.fr/). We thank the referee and Lee Anne Willson for helpful comments.


format: text

text:
# Fcc-bcc transition for Yukawa interactions determined by applied strain deformation

## I Introduction

The screened-Coulomb or Yukawa pair potential, $`U(r)=\mathrm{\Phi }e^{\kappa r}/r`$, has been the focus of great theoretical interest for two reasons. One is that it describes a wide range of interactions, changing continuously from a pure Coulomb potential to an effective hard-sphere potential as the inverse screening length $`\kappa `$ increases. The second is that it provides an approximate description of the effective interactions between large ions that are screened by more mobile counterions. In this context it has been used to describe the interactions between ions surrounded by electrons in metals Ashcroft and Mermin (1976), dust grains surrounded by electrons in dusty plasmas Whipple et al. (1985); Hamaguchi and Farouki (1994); Rosenfeld (1994), and macroions surrounded by counterions in charge-stabilized colloidal suspensions Rosenfeld (1994); Verwey and Overbeek (1948); Alexander et al. (1984); Stevens et al. (1996); Palberg et al. (1995).

The phase diagram of systems of particles interacting with a Yukawa potential has been studied with both analytic Hone et al. (1983); Shih and Stroud (1983); Vaulina et al. (2002); Vaulina (2002); Tejero et al. (1992); Rosenberg and Thirumalai (1987) and numerical Rosenberg and Thirumalai (1987); Meijer and Frenkel (1991); Zahorchak et al. (1992); Stevens and Robbins (1993); Farouki and Hamaguchi (1994); Hamaguchi et al. (1996); Kremer et al. (1986); Robbins et al. (1988); Dupont et al. (1993); Hamaguchi et al. (1997); Miller and Reinhardt (2000) techniques and compared to experiments on dusty plasmas Chu and I (1994); et.al (2003) and colloidal suspensions Monovoukas and Gast (1989); Sirota et al. (1989); Schope et al. (1998). The high temperature phase is a fluid. There is no liquid-gas transition because the interactions are purely repulsive. The stable crystalline phase at zero temperature changes from bcc to fcc as $`\kappa `$ increases. The higher entropy of the bcc phase leads to a greater range of stability as temperature increases until the melting line is reached. Previous results for the fcc-bcc transition line Hone et al. (1983); Kremer et al. (1986); Robbins et al. (1988); Dupont et al. (1993); Hamaguchi et al. (1997); Miller and Reinhardt (2000) vary substantially and the most recent detailed calculation Hamaguchi et al. (1997) quotes an uncertainty of about 10% roughly halfway between the zero-temperature transition point and the triple point.

In this paper we use a different approach to obtain the bcc-fcc phase boundary with an uncertainty of only about 1%. Bounds on the free energy difference between the two phases are obtained by calculating the work done during a continuous deformation between them. The effect of deformation rate, truncation of the potential, and system size and geometry are all analyzed to determine systematic errors. The resulting bcc-fcc transition line is in qualitative agreement with recent simulation results, and quantitative differences are comparable to the larger error bars quoted by previous studies. We estimate the location of the bcc-fcc-fluid triple point using previously published melting-line results Stevens and Robbins (1993); Hamaguchi et al. (1997), and find that it lies at higher inverse screening lengths than previously reported.

The role of anharmonicity in stabilizing the fcc phase is analyzed in detail. While anharmonicity increases the energy of the fcc phase relative to that of the bcc, there is an even larger increase in the relative entropic contribution to the free energy that increases the range of stability of the fcc phase. This appears to reflect an increase in the frequency of the long wavelength shear modes that dominate the bcc entropy in the harmonic approximation Robbins et al. (1988).

Our results are also compared to phenomenological criteria proposed by Vaulina et. al. Vaulina et al. (2002). These authors predict a transition at a critical value of the mean-squared displacement about lattice sites, and calculate the displacement using a simple Einstein-like model. We find that the actual displacement from MD simulations on our transition line is in reasonable agreement with their phenomenological criterion, but substantially larger than predicted by their Einstein model.

The details of our calculations are presented in the following section. Section III provides a detailed analysis of systematic errors and presents our results for the phase boundary. In Section IV, we compare our results to previous transition lines, and Section V provides a summary and conclusions.

## II Method

### II.1 Free energy difference calculations

NVT ensembles are most natural for the study of Yukawa systems for two reasons. First, since the Yukawa potential is purely repulsive, the macroions in an experiment will expand to fill the container. Second, the inverse screening length $`\kappa `$ is density dependent in charged colloidal suspensions and dusty plasmas 111$`\kappa `$ typically depends explicitly on the densities of the mobile counterions that screen interactions between macroions. These densities change with macroion density to maintain charge neutrality.. This density dependence is system-specific, and affects the pressures and bulk moduli. Thus any calculation of coexistence regions will be non-universal. For this reason we focus on finding the Helmholtz free energy difference $`\mathrm{\Delta }F=F_{fcc}F_{bcc}`$ at fixed volume. A brief discussion of coexistence is given in Section II.2.

Postma, Reinhardt, and others Straatsma et al. (1986); Reinhardt and Hunter (1992); Hunter et al. (1993) have shown that the free energy difference between two phases of a system may be calculated in numerical simulations by evaluating the external work done on the system along a thermodynamic path connecting the phases. From elementary thermodynamics, the mechanical work $`W_{AB}`$ done on a system on an isothermal path from state A to state B gives an upper bound on the change $`\mathrm{\Delta }F_{AB}=F_BF_A`$ in the system’s Helmholtz free energy Reif (1965). The work $`W_{BA}`$ done on the system during the reverse process $`BA`$ is an upper bound on $`\mathrm{\Delta }F_{BA}`$, and hence $`(W_{BA})`$ is a lower bound on $`\mathrm{\Delta }F_{AB}`$.

Bounds on $`\mathrm{\Delta }F`$ for Yukawa systems can be obtained using a continuous constant-volume Bain deformation path (Fig. 1) connecting the bcc and fcc lattices. An initially bcc lattice deformed such that its three cubic symmetry directions are scaled respectively by $`(\zeta ^{1/2},\zeta ^{1/2},\zeta )`$ is transformed into an fcc lattice of the same density as $`\zeta `$ varies continuously from 1 to $`2^{1/3}`$ Milstein et al. (1994). We calculate the work done along this path in the forward and reverse directions using strain-controlled molecular dynamics simulations Hoover et al. (1980).

Assuming that the systems traverse these paths homogeneously, we can calculate the work done from the global stresses and strains. We define

$$W_{bf}=V_{bcc}^{fcc}\overline{\sigma }𝑑\overline{ϵ},$$
(1)
$$W_{fb}=V_{fcc}^{bcc}\overline{\sigma }𝑑\overline{ϵ},$$
(2)
where $`\overline{\sigma }`$ and $`\overline{ϵ}`$ are the stress and true strain tensors. $`W_{bf}`$ and $`(W_{fb})`$ are upper and lower bounds on $`\mathrm{\Delta }F`$. For these bounds to be narrow, the intermediate configurations of our systems must remain statistically representative of the $`\zeta `$-dependent equilibrium distributions as $`\zeta `$ is varied Straatsma et al. (1986). In particular, the stress tensor $`\overline{\sigma }(\zeta )`$ must remain near its equilibrium value.

### II.2 Potential parameters

The phase behavior of Yukawa systems is most conveniently expressed in terms of dimensionless screening and temperature parameters. One natural, phase- independent length scale is $`a=n^{1/3}`$, where $`n=N/V`$ is the macroion number density. The Yukawa potential may then be expressed as

$$U(r)=\frac{\mathrm{\Phi }}{a}\frac{e^{\lambda r/a}}{r/a},$$
(3)
where $`\lambda =\kappa a`$ is the dimensionless screening parameter. The limits $`\lambda 0`$ and $`\lambda \mathrm{}`$ correspond to the exhaustively studied one-component plasma and hard-sphere systems.

A natural time scale is provided by $`\tau _E`$, the period of an Einstein oscillator in a crystal. The Einstein periods for the fcc and bcc phases change by an order of magnitude over the range of $`\lambda `$ studied here $`(3\lambda 8)`$, yet differ from each other by less than 1.2% at any given $`\lambda `$ within this range. To obtain consistent results across a wide range of screening lengths, we normalize all time scales in this study to $`\tau _E(\lambda )`$, using the fcc values given in Ref. Robbins et al. (1988).

A natural energy scale is given by the Einstein phonon energies $`m\omega _E^2a^2`$, where m is the macroion mass and $`\omega _E=2\pi /\tau _E`$ is the Einstein frequency. Following Kremer, Robbins, and Grest Kremer et al. (1986), we define the dimensionless temperature

$$\stackrel{~}{T}=\frac{k_BT}{m\omega _E^2a^2}$$
(4)
using the fcc phonon energies, and plot our phase diagram in $`(\lambda ,\stackrel{~}{T})`$ space. A dimensionless inverse temperature $`\mathrm{\Gamma }=(\mathrm{\Phi }/a)/k_BT`$ called the coupling parameter is used in many studies of dusty plasmas. The advantage of using $`\stackrel{~}{T}`$ rather than $`\mathrm{\Gamma }`$ in Yukawa phase diagrams is that the transition lines are approximately linear in $`\lambda `$.

The bcc and fcc phases coexist in equilibrium over a part of the phase diagram. Following previous authors Kremer et al. (1986); Robbins et al. (1988); Dupont et al. (1993); Hamaguchi et al. (1997), we define the bcc-fcc transition line as the curve $`\stackrel{~}{T}_{trans}(\lambda )`$ on which $`\mathrm{\Delta }F=\mathrm{\Delta }F(\lambda ,\stackrel{~}{T})=0`$. This transition line will certainly lie within the coexistence region, regardless of the thermodynamic state dependence of $`\kappa `$ and $`\mathrm{\Phi }`$. We find $`\stackrel{~}{T}_{trans}(\lambda )`$ by calculating $`\mathrm{\Delta }F`$ at many points $`(\lambda _i,\stackrel{~}{T}_i)`$ on the phase diagram.

### II.3 MD simulation details

We simulate NVT ensembles of identical particles using a velocity-Verlet Swope et al. (1982) algorithm to integrate the particle trajectories. The temperature is maintained with a Langevin thermostat Schneider and Stoll (1978). Periodic boundary conditions are used to maintain the density. The equations of motion for the position $`\stackrel{}{q}_i`$ and peculiar momentum $`\stackrel{}{p}_i`$ of the ith particle are

$$\begin{array}{ccc}\dot{\stackrel{}{q}}_i& =& \stackrel{}{p}_i/m+\dot{\overline{ϵ}}\stackrel{}{q}_i,\hfill \\ \dot{\stackrel{}{p}}_i& =& \stackrel{}{F}_i\dot{\overline{ϵ}}\stackrel{}{p}_i+\stackrel{}{\eta }_i\stackrel{}{p}_i/\tau _{Lang},\hfill \end{array}$$
(5)
where $`\dot{\overline{ϵ}}`$ is the true strain rate tensor, $`\stackrel{}{F}_i`$ is the force due to Yukawa interactions, $`\stackrel{}{\eta }_i`$ is a random noise term, and $`\tau _{Lang}`$ is the characteristic relaxation time of the thermostat. We use a timestep $`\delta t=.01\tau _E`$ to insure proper integration of Eqs. (5) and set $`\tau _{Lang}=10\tau _E`$. Changing $`\delta t`$ and $`\tau _{Lang}`$ by a factor of two in either direction had no effect on the phase diagram.

For numerical efficiency we truncate interactions at a cutoff radius $`r_c`$. Due to the presence of long range order in Yukawa crystals, care must be taken in choosing this cutoff radius. We present the details of our determination of $`r_c(\lambda )`$ in Section III.2.

In most of our simulations, we impose the bcc$``$fcc Bain transformation as follows. We start with a lattice of 3456 particles ($`12^3`$ bcc unit cells) in a cubic simulation cell with edges of length $`L_x=L_y=L_z=L_0`$ aligned with the $`<100>`$ directions of the lattice. The system is equilibrated for 200 Einstein periods. We then fix $`\dot{\zeta }=\dot{L_z}/L_0`$ for a time $`\mathrm{\Delta }t`$ sufficient to reach the fcc structure: $`\mathrm{\Delta }t=(2^{\frac{1}{3}}1)\dot{\zeta }^1`$. The other cell edges $`L_x`$ and $`L_y`$ are varied to maintain constant volume and tetragonality $`(L_x=L_y=\sqrt{L_0^3/L_z})`$. The true strain rate tensor $`\dot{\overline{ϵ}}`$ is then given by

$$\begin{array}{c}\dot{ϵ}_{zz}=\dot{L_z}/L_z=\dot{\zeta }/\zeta ,\hfill \\ \dot{ϵ}_{xx}=\dot{ϵ}_{yy}=\dot{L_x}/L_x=\dot{\zeta }/2\zeta ,\hfill \\ \dot{ϵ}_{xy}=\dot{ϵ}_{xz}=\dot{ϵ}_{yz}=0.\hfill \end{array}$$
(6)
We compute the diagonal elements $`(P_x,P_y,P_z)`$ of the pressure tensor using standard methods Allen and Tildesly (1987). Equation (1) then takes on the more physically familiar form

$$W_{bf}=_0^{\mathrm{\Delta }t}\left(P_xL_yL_z\dot{L}_x+P_yL_xL_z\dot{L}_y+P_zL_xL_y\dot{L}_z\right)𝑑t.$$
(7)
After the system has reached the fcc structure, the deformation process is reversed by changing the sign of $`\dot{\zeta }`$. As the system returns to bcc, $`W_{fb}`$ is calculated using the analogue of Eq. (7).

To minimize uncertainties in $`\stackrel{~}{T}_{trans}(\lambda )`$, $`\dot{\zeta }`$ must be small enough for the system to remain near equilibrium. One requirement is that the strain-rate components $`(\dot{\overline{ϵ}}\stackrel{}{q}_i)`$ of the velocities must be small compared to the thermal velocity. The Bain transformation time $`\mathrm{\Delta }t`$ (which is proportional to $`\dot{\zeta }^1`$) must also be large compared to $`\tau _{Lang}`$ to allow the thermostat to transfer heat to or away from the system as necessary to maintain constant temperature. Since $`\tau _{Lang}`$ sets the time over which the system samples the canonical ensemble, the thermodynamic sampling improves as $`\mathrm{\Delta }t/\tau _{Lang}`$ increases.

The precision of the calculated transition line depends on the difference $`(\mathrm{\Delta }F_{max}\mathrm{\Delta }F_{min})W_{bf}+W_{fb}=W_{cycle}`$ between the bounds on $`\mathrm{\Delta }F`$. These bounds converge to each other in the reversible thermodynamic (zero strain rate) limit. In this limit, the average work $`<W_{cycle}>`$ done on the system over a full deformation cycle (bcc$``$fcc$``$bcc or vice versa) should vanish. In simulations at finite strain rate, however, there is a positive systematic error in $`W_{cycle}`$ due to energy dissipation Jarzynski (1997). This can be physically interpreted as arising from viscosity. Each applied strain increment takes the system slightly out of equilibrium. When $`\dot{\zeta }`$ is small, one expects the stresses to deviate from their equilibrium values by an amount $`\overline{\sigma }_{visc}\dot{\zeta }`$ Larson (1999). Sources of viscous dissipation include the intrinsic viscosity and the drag forces $`\stackrel{}{p}_i/\tau _{Lang}`$ on the particles applied by the Langevin thermostat. The viscous dissipation rate is given by $`\overline{\sigma }_{visc}\dot{\overline{ϵ}}`$, so one expects the dissipated power to be proportional to $`\dot{\zeta }^2`$. Since the total simulation time scales as $`\dot{\zeta }^1`$, the total dissipated energy, and hence the deviation of $`<W_{cycle}>`$ from zero, should be linearly proportional to $`\dot{\zeta }`$. We present the $`\dot{\zeta }`$-dependence of our results in Section III.1.

## III Results

### III.1 Strain rate dependence

At temperatures near the transition line, calculations of the geometrical structure factor and pair correlation function verify that our systems traverse the Bain transformations homogeneously. This homogeneity allows us to use Eqs. (1,2) for calculating $`W_{bf}`$ and $`W_{fb}`$ and leads to tight bounds on $`\mathrm{\Delta }F`$.

We obtained results similar to those shown in Figure 2 over the entire investigated range of $`\lambda `$ and for three different values of $`\dot{\zeta }`$. Near the transition line, $`W_{bf}(\lambda ,\stackrel{~}{T})`$ and $`W_{fb}(\lambda ,\stackrel{~}{T})`$ vary linearly with $`\stackrel{~}{T}`$ and have nearly opposite ($`\lambda `$-dependent) slopes. The scatter about linear fits to $`W_{bf}(\lambda ,\stackrel{~}{T})`$ and $`W_{fb}(\lambda ,\stackrel{~}{T})`$ is consistent with fluctuations in $`W_{bf}`$ and $`W_{fb}`$ at fixed $`(\lambda ,\stackrel{~}{T})`$. The intersections of these fits with $`W=0`$ give two estimates, $`\stackrel{~}{T}_{bf}`$ and $`\stackrel{~}{T}_{fb}`$, for $`\stackrel{~}{T}_{trans}`$. These are obtained using data at ten evenly spaced $`\stackrel{~}{T}`$ within about 5% of the transition line.

For a given system, $`\stackrel{~}{T}_{fb}`$ and $`\stackrel{~}{T}_{bf}`$ provide upper and lower bounds on $`\stackrel{~}{T}_{trans}`$ since $`W_{bf}>\mathrm{\Delta }F>W_{fb}`$ and $`(\mathrm{\Delta }F/\stackrel{~}{T})>0`$. We define the fractional uncertainty due to dissipative hysteresis as

$$\delta _{hyst}=(\stackrel{~}{T}_{fb}\stackrel{~}{T}_{bf})/(\stackrel{~}{T}_{fb}+\stackrel{~}{T}_{bf}).$$
(8)
The conditions $`\delta _{hyst}\dot{\zeta }`$ and $`<W_{cycle}>\dot{\zeta }`$ are equivalent due to the linear dependence of $`W_{bf}`$ and $`W_{fb}`$ on $`\stackrel{~}{T}`$. Figure 3 shows the $`\dot{\zeta }`$-dependence of $`\delta _{hyst}`$ for $`\lambda =4`$ and $`\lambda =5`$. The results are consistent with our hypothesis that the energy dissipated is linear in $`\dot{\zeta }`$.

We identify $`\stackrel{~}{T}_{trans}=(\stackrel{~}{T}_{bf}+\stackrel{~}{T}_{fb})/2`$ as the best estimate for the transition temperature for a given system size, system geometry, and potential cutoff radius. Table 1 shows that the fractional variation of $`\stackrel{~}{T}_{trans}`$ with $`\dot{\zeta }`$ is much smaller than $`\delta _{hyst}`$ 222These values are for the standard system size and geometry described in Section II.3 and cutoff radii reported in Section III.2..

In the following, we present results for $`|\dot{\zeta }|=10^4/\tau _E`$. Based on Table 1, for this value of $`|\dot{\zeta }|`$ the random and finite strain rate uncertainties in $`\stackrel{~}{T}_{trans}`$ are comparable, both about 0.2%. The combined error is estimated to be less than 0.4%. The uncertainties given in subsequent tables include only statistical uncertainties from the linear fits used to calculate $`\stackrel{~}{T}_{bf}`$ and $`\stackrel{~}{T}_{fb}`$.

### III.2 Potential cutoff dependence

We estimate the errors introduced by truncating the force at $`r_c`$ by calculating the error $`\delta E(r_c)`$ in the potential energy difference. If the error in $`\mathrm{\Delta }F`$ is of the same order, then the fractional error in the transition temperature is

$$\delta _{cut}\frac{\delta \stackrel{~}{T}(r_c)}{\stackrel{~}{T}}\frac{1}{\stackrel{~}{T}}\frac{\stackrel{~}{T}}{\mathrm{\Delta }F}\delta E(r_c).$$
(9)
Here $`(\stackrel{~}{T}/\mathrm{\Delta }F)`$ is known near the transition line from the work calculations.

The cutoff-induced error in the potential energy difference can be written in terms of the pair correlation functions $`g_{fcc}(r)`$ and $`g_{bcc}(r)`$ of the fcc and bcc crystals. For N particles

$$\delta E(r_c)=\frac{N}{2}_{r_c}^{\mathrm{}}U(r)(g_{fcc}(r)g_{bcc}(r))4\pi r^2𝑑r/a^3$$
(10)
If $`\mathrm{\Omega }\text{max}\{|g_{fcc}(r)g_{bcc}(r)|;rr_c\}`$, then

$$\frac{\left|\delta E\left(r_c\right)\right|}{N\left(\mathrm{\Phi }/a\right)}<\frac{\mathrm{\Omega }}{\mathrm{\Phi }a^2}_{r_c}^{\mathrm{}}U\left(r\right)2\pi r^2𝑑r=\frac{2\pi \mathrm{\Omega }\left(1+\lambda r_c/a\right)e^{\lambda r_c/a}}{\lambda ^2}$$
(11)
One expects $`\mathrm{\Omega }`$ to be of order 1 at finite temperature.

For $`\lambda `$ = 3, 4, and 5, we estimated $`\delta E(r_c)`$ at $`\stackrel{~}{T}\stackrel{~}{T}_{trans}`$ by calculating the pair correlation functions in large systems using large cutoff radii and long integration times. Due to the exponential falloff of U(r) and finite temperature smoothing of g(r), the infinite upper bound in Eq. 10) can be replaced by a finite value $`r_l`$ without introducing significant errors. We found $`\lambda r_l=30a`$ to be sufficiently large.

Figure 4 shows our estimate of $`|\delta E(r_c)/(N\mathrm{\Phi }/a)|`$ from Eq. (10) and the value of $`|\delta E/(N\mathrm{\Phi }/a)|`$ corresponding to $`\delta _{cut}=0.01`$ for $`\lambda =3`$, the longest-range potential considered. We found that $`\mathrm{\Omega }`$ decreases from 1.4 to .62 as $`r_c`$ increases from 3.5a to 6.5a. The actual error is always smaller than the bound given by $`\mathrm{\Omega }`$ because $`g_{fcc}g_{bcc}`$ oscillates in sign. The envelope shown corresponds to $`\mathrm{\Omega }=1/3`$. Because of the sharp variation of $`\delta E(r_c)`$, we use the envelope to estimate $`\delta _{cut}`$.

To test Eq. (9) we calculated $`\stackrel{~}{T}_{trans}`$ as a function of $`r_c`$ for $`\lambda =3`$. Results are shown in Table 2. The fractional changes in $`\stackrel{~}{T}_{trans}`$ from $`r_c=3.5a`$ and $`r_c=4.667a`$ to $`r_c=6.667a`$ are 19% and 0.9%, respectively. Both changes are about one-fifth of the estimates for $`\delta _{cut}`$ from Eqs. (9,10). No statistically significant changes are expected or observed for $`r_c5.833a`$. We conclude that errors estimated from the envelopes of curves like Fig. 4 give a conservative estimate of cutoff errors.

To ensure that the fractional systematic errors were no larger than our random and rate errors, we chose $`r_c`$ slightly above the values corresponding to $`|\delta _{cut}|=0.002`$. For $`\lambda `$ = 3, 4, and 5, we used cutoff radii of 5.833a, 4.375a, and 3.5a in the simulations used to determine $`\stackrel{~}{T}_{trans}`$. Smaller $`r_c`$ can be used at higher $`\lambda `$ both because the interactions weaken and $`\stackrel{~}{T}_{trans}/\stackrel{~}{T}_{melt}`$ increases, leading to a smaller $`\mathrm{\Omega }`$. For $`\lambda 5`$ we fixed the cutoff radius at $`r_c=3.5a`$.

### III.3 System size and geometry dependence

To examine finite size effects we also considered a 432-particle system (initial state $`6^3`$ bcc unit cells). Because the corresponding fcc state has transverse length $`6.73a`$, the minimum image convention requires $`r_c<3.367a`$, and we used $`r_c=3.3a`$. To separate out $`r_c`$-dependence from system size dependence, we also recalculated the transition line for $`N=3456`$ for $`5\lambda 8`$ with $`r_c=3.3a`$.

Table 3 shows a comparison of our calculated transition temperatures. The $`N=432`$ values were systematically lower, but the effect was small. From theoretical considerations one expects the leading finite size corrections to $`\mathrm{\Delta }F/N`$ to be proportional to $`1/N`$ Polson et al. (2000). This should produce a corresponding error in $`\stackrel{~}{T}_{trans}`$. As shown in Table 3, the changes in $`\stackrel{~}{T}_{trans}`$ from $`N=432`$ to $`N=3456`$ were all about 1%. The changes in $`\stackrel{~}{T}_{trans}`$ from $`N=3456`$ to $`N=\mathrm{}`$ for this system geometry should be about 8 times smaller.

Another test indicates that finite size effects are larger than the above estimate. The geometry was changed so that the fcc state has equal cell edges and the bcc state has $`L_x=L_y=\sqrt{2}L_z`$. These simulations contained $`10^3`$ fcc unit cells (4000 particles) with the $`<100>`$ directions parallel to the simulation cell edges. After Bain transformation, the bcc state has two $`<110>`$ directions parallel to the simulation cell edges. As shown in Table 4, the values of $`\stackrel{~}{T}_{trans}`$ obtained for both $`\lambda =4`$ and $`\lambda =7`$ were 0.6% lower than those obtained with the standard system geometry 333Both sets of simulations used the cutoff radii indicated in Section III.2.. Other simulations verified that this was due solely to the change in boundary conditions. We conclude that our dominant source of uncertainty is finite size and is less than 1%.

We attribute the observed sensitivity to geometry to the change in allowed low frequency modes. These modes play a disproportionate role in determining the entropy in lattice dynamics calculations Robbins et al. (1990) and drive the fcc$``$bcc transition with increasing temperature Robbins et al. (1988). Since the shear velocity is highly anisotropic in the bcc phase, changing the boundaries affects the sampling of these low frequency modes and thus $`\mathrm{\Delta }F`$.

### III.4 Transition line

Table 5 shows our calculated $`\stackrel{~}{T}_{trans}(\lambda )`$ with statistical uncertainties. As described above, the combined systematic errors due to finite strain rate, system size, and potential cutoff are estimated to be less than 1%. Results of a cubic polynomial fit to the data are also given:

$$\begin{array}{ccc}10^4\stackrel{~}{T}_{trans}^{fit}(\lambda )\hfill & =& 6.46678(\lambda \lambda _0)+0.43001(\lambda \lambda _0)^2\hfill \\ & & 0.06806(\lambda \lambda _0)^3,\hfill \end{array}$$
(12)
where $`\lambda _0=1.718`$ is the zero-temperature transition point obtained from lattice statics calculations e Silva and Mockross (1980). Lower-order polynomials fail to adequately fit the data within our uncertainties.

Figure 5 shows the polynomial fit and two previously published solid-fluid coexistence lines Stevens and Robbins (1993); Hamaguchi et al. (1997). The intersections of these lines give estimated values of the bcc-fcc-fluid triple point. Using results from Ref. Stevens and Robbins (1993) we find $`(\lambda _{tp}=7.45,\stackrel{~}{T}_{tp}=0.00384)`$. Those from Ref. Hamaguchi et al. (1997) yield $`(\lambda _{tp}=7.84,\stackrel{~}{T}_{tp}=0.00401)`$.

If we assume that the parameters $`\kappa `$ and $`\mathrm{\Phi }`$ are density-independent, we can calculate the width of the bcc-fcc coexistence region from the pressures and bulk moduli of the two phases on the line where $`\mathrm{\Delta }F=0`$. The bcc pressure is larger than the fcc pressure by only about 0.04% for $`\lambda =4`$, and by 0.65% for $`\lambda =7`$. This results in a higher density in the fcc phase at coexistence, but only by about 0.015% at $`\lambda =4`$ and 0.2% at $`\lambda =7`$. The corresponding changes in $`\lambda `$ and $`\stackrel{~}{T}`$ are much smaller than the uncertainties in our calculated transition line. The coexistence region in experimental systems may be much larger due to variations in $`\kappa `$ and $`\mathrm{\Phi }`$ with density Schope et al. (1998); et.al (2003). As noted above, these variations are system specific and a more complete treatment is beyond the scope of this paper.

### III.5 Anharmonic effects

The dotted line in Figure 5 shows lattice dynamics results for the bcc-fcc transition line Robbins et al. (1988). In this approximation the energy and entropy differences, $`\mathrm{\Delta }E_{LD}`$ and $`\mathrm{\Delta }S_{LD}`$, are independent of T. The fcc-bcc transition line is given by $`T_{LD}=\mathrm{\Delta }E_{LD}/\mathrm{\Delta }S_{LD}`$. The resulting curve lies below $`\stackrel{~}{T}_{trans}`$, indicating that the fcc phase is stabilized by anharmonic effects Robbins et al. (1988); Hamaguchi et al. (1997). This implies that the anharmonic component of the free energy difference,

$$\mathrm{\Delta }F_{an}\mathrm{\Delta }E_{an}T\mathrm{\Delta }S_{an}=\mathrm{\Delta }F\mathrm{\Delta }F_{LD},$$
(13)
is negative on the transition line. The relative signs and magnitudes of $`\mathrm{\Delta }E_{an}`$ and $`\mathrm{\Delta }S_{an}`$ may be calculated by comparing our accurate measurements of free and total energy differences with the lattice-dynamics results.

Table 6 shows results for anharmonic contributions to the free and total energy differences on the fit transition line 444Results for $`\lambda =8`$ are not presented in Tables 6 and 7 respectively because at $`\stackrel{~}{T}_{trans}`$ solid diffusion prohibits accurate calculations of the rms displacements, and lattice dynamics calculations of $`\mathrm{\Delta }S`$ are unavailable.. The values of $`\mathrm{\Delta }F_{an}`$ are known from the work calculations, while the values of $`\mathrm{\Delta }E_{an}`$ were found from separate equilibrium simulations. The anharmonic corrections to the total energy favor the bcc phase for all $`\lambda `$, i.e., $`E_{fcc}E_{bcc}`$ on the transition line has increased relative to its zero-temperature value. The anharmonic contributions to the free energy difference, however, are larger in magnitude and opposite in sign, implying that anharmonic entropic contributions to $`\mathrm{\Delta }F`$ favor the fcc phase at all $`\lambda `$ and overwhelm energetic contributions.

In lattice calculations the larger entropy of the bcc phase comes mainly from the lower frequency of its shear modes. Some of these modes have negative energy for $`\lambda >7.67`$, causing the bcc phase to become linearly unstable at low temperatures Robbins et al. (1988). It is interesting that the onset of this low temperature instability is close to $`\lambda _{tp}`$. However, we have performed runs near the melting line for $`\lambda `$ as large as 10 and find that the bcc phase remains metastable. This implies that anharmonic effects have increased the frequency of long wavelength shear modes, providing an explanation for the decreased entropy advantage of the bcc phase.

## IV Comparison to Previous Results

Miller and Reinhardt were the first authors to use Bain deformation paths to obtain bounds on $`\mathrm{\Delta }F`$ for Yukawa systems Miller and Reinhardt (2000). They calculated the work by integrating the change in the Hamiltonian rather than from the stresses and strains. The large discrepancy between their $`\lambda =7`$ transition temperature and our result is likely due to their extremely small system size $`(N=108)`$, which was just used to illustrate their method.

The earliest MD calculations Kremer et al. (1986); Robbins et al. (1988) of the transition line also deviate substantially from ours, particularly at large $`\lambda `$. The line shown in Figure 5 is a fit between points where the fcc and bcc phases were found to be stable. The gap between points was about 20% and the final shape was strongly influenced by a bcc-stable point above the melting line. Other points where the bcc phase was stable lie close to our $`\stackrel{~}{T}_{trans}`$ but are shifted up due to the smaller $`r_c`$ used.

Our transition line is in qualitative agreement with more recent MD and Monte Carlo results Hamaguchi et al. (1997); Dupont et al. (1993). Dupont et. al. calculated a fcc-bcc coexistence point and the fcc-bcc-fluid triple point using small systems $`(N250)`$. Although their triple point $`(\lambda _{tp}=6.75,\stackrel{~}{T}_{tp}=0.0034)`$ lies well below ours $`(\lambda _{tp}=7.7\pm 0.3,\stackrel{~}{T}_{tp}=0.0039\pm 0.0001)`$, it lies only about 2% below our fcc-bcc transition line, and well below recently published melting lines Stevens and Robbins (1993); Hamaguchi et al. (1997).

Hamaguchi, Farouki, and Dubin also obtain a lower triple point $`(\lambda _{tp}=6.90,\stackrel{~}{T}_{tp}=0.0038)`$ because their bcc-fcc transition temperatures are systematically (6-10%) higher than ours Hamaguchi et al. (1997). One possible explanation is that their equilibration times were too short. They used the $`\lambda `$-independent time unit $`\tau =\sqrt{3}\omega _p^1`$, where $`\omega _p=\sqrt{4\pi n\mathrm{\Phi }/m}`$ is the plasma frequency. Starting with perfect bcc and fcc lattices as their initial conditions, they equilibrated their systems for a maximum of 300$`\tau `$ before beginning their free energy measurements. This corresponds to about 27$`\tau _E`$ for $`\lambda =3`$ and only 4$`\tau _E`$ for $`\lambda =8`$. Since the latter is only about four times the velocity autocorrelation time, and comparable to the time for sound to propagate across their simulation cells, it is doubtful that their systems had equilibrated sufficiently at high $`\lambda `$. Too short an equilibration time could cause overestimation of the stability of the phase with lower entropy, the fcc phase, which is consistent with their findings.

Vaulina and colleagues have proposed phenomenological criteria for the bcc-fcc transition Vaulina et al. (2002). They assume that $`\kappa ^1`$ is an effective hard-sphere radius and predict that the value of the rms displacement at the fcc transition, $`\mathrm{\Delta }_{trans}`$, satisfies

$$2(1\pi \sqrt{2}/6)^{1/3}\mathrm{\Delta }_{trans}=R_{WS}\kappa ^1$$
(14)
where $`R_{WS}=(4\pi /3)^{1/3}a`$ is the Wigner-Seitz radius. They then use an approximate formula for the effective frequency in an Einstein-like model to determine $`\mathrm{\Delta }`$ for the fcc and bcc structures. These values of $`\mathrm{\Delta }`$ give two predictions for the transition line. As shown in Figure 6, their predictions are qualitatively correct but lie roughly 10-40% above our $`\stackrel{~}{T}_{trans}`$.

The discrepancy in Fig. 6 could be due to a failure either of Eq. (14) or of the approximations used to find $`\mathrm{\Delta }`$. To test this we performed equilibrium simulations at $`\stackrel{~}{T}_{trans}^{fit}`$ in both bcc $`(N=3456)`$ and fcc $`(N=4000)`$ systems. Our results for the rms displacements, $`\mathrm{\Delta }_{fcc}`$ and $`\mathrm{\Delta }_{bcc}`$, are compared to the predictions of Eq. (14) in Table 7. The rms displacements for the fcc structure lie quite close to the prediction for small $`\lambda `$, and about 13% above it at $`\lambda =7`$. Finite size effects decrease $`\mathrm{\Delta }_{fcc}`$ relative to the N = $`\mathrm{}`$ value Robbins et al. (1990)555Equilibrium simulations of larger (N $``$ 100000) systems for $`\lambda =5`$ showed that both $`\mathrm{\Delta }_{fcc}`$ and $`\mathrm{\Delta }_{bcc}`$ increase with N, worsening the agreement with $`\mathrm{\Delta }_{trans}`$, but that the finite-size corrections are small compared to the discrepancies shown in Table VII. . These results indicate that most of the error in Vaulina et. al.’s transition lines comes from substantial underestimation of the rms displacements. They calculate $`\mathrm{\Delta }`$ in the harmonic approximation, and anharmonic corrections increase $`\mathrm{\Delta }`$ for these $`\lambda `$ Robbins et al. (1988). Note that the measured bcc displacements in Table VI are larger due to the bcc lattice’s softer shear modes Robbins et al. (1988).

## V Summary and Conclusions

We calculated the bcc-fcc coexistence line of Yukawa systems to an uncertainty of approximately 1% through integration of the mechanical work along Bain transformation paths. The range of bcc stability was found to be slightly greater than that found in previous comprehensive studies Robbins et al. (1988); Dupont et al. (1993); Hamaguchi et al. (1997), and the triple point shown to lie at higher inverse screening length. The large changes in $`\stackrel{~}{T}_{trans}`$ with $`r_c`$ for small $`\lambda `$ indicate that the relative stability of fcc and bcc phases depends sensitively on long-range correlations, and calls into question the use of local nearest-neighbor arguments to calculate the transition line. Nevertheless, we found that one such phenomenological criterion Vaulina et al. (2002), derived from the idea that the fcc phase is stable when interparticle interactions are hard-sphere-like Hoover and Ree (1968), predicts the transition line remarkably well when combined with MD results for the mean-squared displacement.

Comparison with lattice-dynamics results shows that anharmonic terms in the total energy favor the bcc phase for all $`\lambda `$, but that these corrections are overbalanced by anharmonic contributions to the entropy. The change in entropy appears to reflect an increase in the frequency of long-wavelength shear modes in the bcc phase. This increase also stabilizes the bcc phase against a linear shear instability observed for $`\lambda >7.67`$ at low $`\stackrel{~}{T}`$.

We found that shifts in the transition line due to finite size effects are less than 1% if N$``$3000-4000, but that the presence of long-range order at temperatures near the transition line in weakly screened systems requires a cutoff radius larger than that used in some previous studies Kremer et al. (1986); Robbins et al. (1988). Accurate simulations of weakly screened ($`\lambda <3`$) systems in this temperature range require either larger system sizes or an Ewald-like summation over periodic images Hamaguchi et al. (1997). However, we have also shown that a reasonably small potential cutoff need not introduce large errors in a transition line calculation in the moderate-screening regime, provided the cutoff is chosen with some care.

It is known that the phase behavior of real systems such as charge-stabilized colloidal suspensions is not fully described by pointlike Yukawa interactions. Recent simulations of charged macroions in a dynamic neutralizing background have shown that the repulsive interactions between macroions are truncated by many-body effects, destabilizing the crystalline phases in the weak screening limit Dobnikar et al. (2003a, b); Hynninen and Dijkstra (2003a). Independently, hard-core repulsions significantly alter the phase diagram when the volume fraction is more than a few percent Meijer and Azhar (1997); Azhar et al. (2000); Hynninen and Dijkstra (2003b). However, we hope that our high-precision calculation of the point Yukawa fcc-bcc transition line may serve as a benchmark for further studies of more sophisticated models.

## VI Acknowledgements

The simulations in this paper were carried out using the LAMMPS molecular dynamics software. Support from NSF Grant DMR-0083286 is gratefully acknowledged.