text:
# Diffraction-limited speckle interferometry and modeling of the circumstellar envelope of R CrB at maximum and minimum light

## 1 Introduction

The R Coronae Borealis (RCB) stars are a class of unusual objects characterized by sudden declines in their visible light curves as deep as $`\mathrm{\Delta }V8`$. They are extremely hydrogen-deficient and also carbon-rich (e.g. Asplund et al. asplund00 (2000) and references therein). The RCB stars are thought to undergo the formation of dust clouds in random directions, and it is believed that a sudden decline takes place, only when a dust cloud forms in the line of sight (Loreta loreta34 (1934), O’Keefe okeefe39 (1939)). A newly formed dust cloud is expected to be accelerated by radiation pressure and to move away, expanding and dispersing over months, as the object gradually returns to its maximum visual brightness. The mechanism of the dust cloud formation and its temporal evolution are, however, still poorly understood. Especially, the location of dust formation is in dispute: far from the star, $`20R_{}`$ (e.g. Fadeyev fadeyev86 (1986), fadeyev88 (1988), Feast feast96 (1996)), or very close to the photosphere, $`2R_{}`$ (Payne-Gaposchkin PG63 (1963)). Recent photometric and spectroscopic observations as well as theoretical progress on dust formation suggest that the latter scenario may be the case (see, e.g. Clayton clayton96 (1996), Feast feast97ii (1997)), however, no definitive answer is yet available.

The spectral energy distributions (SEDs) of RCB stars exhibit infrared emission peaks around 6 $``$$`\mu `$m. The IR excess, which accounts for typically 30% of the total flux, is constantly present, regardless of the visual brightness of the central star. Therefore, the IR excess originates not from a single newly formed dust cloud, but mainly from a group of dispersed dust clouds with temperatures of approximately 600 – 900 K. For example, Walker et al. (walker96 (1996)) fit the infrared (spectro)photometric data of R CrB with a 650 K blackbody. The study of IRAS observations by Gillett et al. (gillett86 (1986)) led to the detection of an additional, very extended “fossil” shell around R CrB, whose diameter is as large as 18′, with a temperature of $`30`$ K. Walker (walker94 (1994)) found such fossil shells for at least four RCB stars.

Apart from the detection of the fossil shells, most of the observational results on the circumstellar environment around RCB stars were obtained by photometry and spectroscopy. Recently, Clayton & Ayres (clayton01 (2001)) have revealed extended C ii $`\lambda `$1335 emission around two RCB stars, V854~Cen and RY~Sgr, by long-slit spectroscopy. However, such direct information on the spatial distribution of material in the vicinity of the central star has been very rare up to now. In this paper, we present high-resolution speckle interferometry carried out for R CrB at maximum and minimum light. The properties of the warm dust shell will be derived by simultaneous fits of the observed visibilities and SEDs using power-law models. We will also discuss the possible indication of a newly formed hot dust cloud.

## 2 Speckle interferometric observations

The $`K`$-band speckle interferometric observations were carried out with the 6 m telescope at the Special Astrophysical Observatory (SAO) in Russia, using our NICMOS-3 camera for the observation in 1996 and our HAWAII array speckle camera in 1999. The observations are described in Table 1.

Fig. 1 shows the visual light curve of R CrB in the relevant period based on the compiled data of the American Association of Variable Star Observers (AAVSO). The data before 1996 were taken from AAVSO Monograph 4 Supplement 1, while J. A. Mattei kindly provided unpublished data for the rest of the period. On 1996 October 1, R CrB was on its final recovery from the 1995-96 deep minimum and its visual magnitude was around 7 (Mattei mattei00 (2000)). Since $`V`$ is approximately 6 at maximum light, the star was slightly obscured by $`\mathrm{\Delta }V=1`$. On 1999 September 28, the star was just between two sharp deep minima and had a visual magnitude of 10.61 (see Sect. 3). The star was heavily obscured by $`\mathrm{\Delta }V=4.6`$.

Fig. 2 shows the two-dimensional power spectra and azimuthally averaged visibility functions reconstructed with the speckle interferometry method (Labeyrie labeyrie70 (1970)). The error bars include systematic errors caused by seeing variations as well as speckle noise errors. There is no significant deviation from circular symmetry in the two-dimensional power spectra observed at both epochs. The granular features seen in Fig. 2a are the speckle noise resulting from the small number of interferograms. In Fig. 2b, a slight elongation along the vertical axis of the figure can be marginally recognized. However, this is within the error bars, and cannot be regarded as a significant feature. Clayton et al. (clayton97 (1997)) suggest a bipolar geometry for R CrB based on the wavelength dependence of the position angle of polarization. Our speckle observations seem to support a rather symmetric distribution of material. However, a disc or torus may lie very close to the star, and therefore, it is premature to rule out a bipolar geometry completely.

The reconstructed visibility functions shown in Figs. 2c and 2d exhibit no difference larger than the error bars. The source probably consists of the central star (a point source) and an extended dust shell. Visibilities observed for such objects have a plateau at high spatial frequencies resulting from the point source. In the observations presented here, the angular size of the shell is small, and therefore, the plateau is not visible within a cut-off frequency of $``$ 13 cycles/arcsec. The determination of the physical parameters of the shell by simultaneous fits to the observed visibilities and SEDs will be presented in Sect. 5.

## 3 Photometric data

$`UBVRIJHKLM`$ photometry was carried out on 1999 September 29, just one night after the speckle observation. The optical ($`UBVRI`$) data were obtained with the 1.25 m telescope at the Crimean Astrophysical Observatory, and the infrared ($`JHKLM`$) data with the 1.22 m telescope at the Crimean Laboratory of the Sternberg Astronomical Institute.

No photometric data, except for the visual magnitudes compiled in the AAVSO database, have been published around the date of our speckle observation in 1996 October. Therefore, we use the data published by Shenavrin et al. (shenavrin79 (1979)), who report a series of photometric observations during the 1977 deep minimum. We use the $`UBVRJHK`$ data on 1978 January 22 (JD2443530.6), when the visual magnitude was 7.03, almost the same as that on the date of our speckle observation in 1996 October. For the $`L`$\- and $`M`$-bands, we use photometric data obtained on 1996 June 11 (JD2450246.4) by Shenavrin (unpublished observation), about 4 months before our speckle observation. It is known that the variation of the $`L`$\- and $`M`$-band fluxes has no correlation with the visual light curve, and that they exhibit a semi-periodic variation of $`1.5`$ mag with a period of about 1260 days (e.g. Feast et al. feast97 (1997)). In fact, the $`L`$ and $`M`$ magnitudes on 1996 June 11 are by 1.09 mag and 1.25 mag brighter than those on 1999 September 29, respectively. The $`L`$\- and $`M`$-band fluxes measured 4 months before are expected to represent the values on the date of the speckle observation rather well.

We also use IRAS observations at 12, 25, 60, and 100 $`\mu `$m. The very extended “fossil” shell gives rise to flux excess at 60 and 100 $`\mu `$m, when evaluated with a large aperture (Gillett et al. gillett86 (1986)). In this paper, however, we are primarily interested in the warm dust shell which was not resolved with IRAS, therefore, we adopt the point source processed fluxes for 60 and 100 $`\mu `$m. Long-term variations at wavelengths longer than 10 $`\mu `$m are considered to be smaller than those in the $`L`$\- and $`M`$-bands. Forrest et al. (forrest72 (1972)) show that the amplitude at 11 $`\mu `$m is about half as small as those in the $`L`$\- and $`M`$\- bands. The variation at wavelengths longer than 11 $`\mu `$m is estimated to be $`0.75`$ mag at most.

The interstellar extinction toward R CrB is $`E(BV)=0.05`$ (Asplund et al. asplund97 (1997)). The observed fluxes are dereddened using the method of Savage & Mathis (savage79 (1979)) with $`A_V=3.1E(BV)`$.

## 4 Description of the model

The basics of the model used here are the same as adopted by Gillett et al. (gillett86 (1986)). Fig. 3 illustrates the picture considered in our modeling. It consists of the central star and an optically thin dust shell responsible for the constant IR excess. As mentioned in Sect. 2, the star was obscured by $`\mathrm{\Delta }V=4.6`$ (1999 September 28) and $`\mathrm{\Delta }V=1.0`$ (1996 October 1). Therefore, an additional obscuring dust cloud, which is assumed to be circular (projected onto the sky as seen from the star), is also considered, as shown in Fig. 3.

Adopting $`M_{\mathrm{bol}}`$ and $`T_{\mathrm{eff}}`$ to be $`5.3`$ (Gillett et al. gillett86 (1986)) and 6600 – 6900 K (Asplund et al. asplund00 (2000)), respectively, the radius of the central star is approximately 70 $`R_{\mathrm{}}`$. We adopt 6750 K for the effective temperature, and radiation from the central star is approximately represented by a blackbody. The distance to R CrB adopted here is 1.6 kpc (Gillett et al. gillett86 (1986)), yielding an angular radius of the central star of 0.23 mas.

We consider a spherical shell with an inner radius $`r_{\mathrm{in}}`$ and outer radius $`r_{\mathrm{out}}`$. The shell is assumed to be optically thin not only in the infrared but also in the optical, since there is no evidence for strong circumstellar reddening for R CrB. Asplund et al. (asplund97 (1997)) found that the observed optical flux at maximum light can be well reproduced by their line-blanketed model atmospheres without any circumstellar extinction. In fact, the optical depth of the dust shell in the $`V`$-band is $`0.3`$ in our models presented below. The assumption of spherical symmetry for the shell is justified, given the negative detection of deviation from circular symmetry in the observed two-dimensional power spectra.

The grain number density is assumed to decrease in a power-law of radius,

$$n(r)=Cr^\gamma ,$$
(1)
where $`C`$ is a constant which does not affect the shape of emergent SEDs, and therefore, should be adjusted to fit the flux observed on the earth.

The temperature in an optically thin dust shell is determined by the thermal balance equation

$`{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{L_\nu }{4\pi r^2}}\pi a^2Q_{\mathrm{abs}}(a,\nu )𝑑\nu `$
$`={\displaystyle _0^{\mathrm{}}}4\pi a^2Q_{\mathrm{abs}}(a,\nu )\pi B_\nu (T(r))𝑑\nu ,`$ (2)
where $`L_\nu `$ is the luminosity of the central star at a given frequency, $`a`$ is the radius of a grain, $`Q_{\mathrm{abs}}(a,\nu )`$ is the absorption efficiency of a grain, and $`B_\nu (T(r))`$ is the Planck function. Amorphous carbon is the most probable candidate for the circumstellar dust around RCB stars (e.g. Holm et al. holm82 (1982), Hecht et al. hecht84 (1984)). In Fig. 4, we show temperature distributions predicted by the thermal balance equation, using the extinction of amorphous carbon obtained by Bussoletti et al. (bussoletti87 (1987)) (AC2 sample), Rouleau & Martin (rouleau91 (1991)) (AC1 sample), and Colangeli et al. (colangeli95 (1995)) (ACAR sample). We calculate $`Q_{\mathrm{abs}}`$ from the complex refractive index derived by Rouleau & Martin (rouleau91 (1991)) in the Mie theory, assuming a single grain size, $`a=0.01`$ $`\mu `$m, and using the code published by Bohren & Huffman (bohren83 (1983)). This grain size is based on the result obtained by Hecht et al. (hecht84 (1984)), who analyzed the UV spectra of R CrB and RY Sgr and concluded that the grain size is between 0.005 and 0.06 $`\mu `$m. As Fig. 4 shows, the temperature distributions agree with one another within $`50`$ K.

The flux density observed on the earth from the optically thin shell described above can be calculated by

$$f_\mathrm{s}(\lambda )=\frac{4\pi m_\mathrm{d}\kappa _\lambda C}{D^2}_{r_{\mathrm{in}}}^{r_{\mathrm{out}}}B_\lambda (T(r))r^{2\gamma }𝑑r,$$
(3)
where $`m_\mathrm{d}`$ and $`\kappa _\lambda `$ is the mass and mass absorption coefficient of a grain, respectively, and $`D`$ is the distance to the star. The constant $`C`$ is adjusted so that the flux predicted from the models can reproduce the observed flux at 25 $`\mu `$m, since its long-term variation mentioned in Sect. 3 is expected to be minimal, and at the same time, it does not exhibit the contribution of the fossil shell. The variation of the 25 $`\mu `$m flux is not available in the literature, but it is estimated to be smaller than the 0.75 mag at 11 $`\mu `$m. $`r_{\mathrm{out}}`$ is set to be $`1.5\times 10^4`$ $`R_{}`$, which gives a good match to the far-infrared part of the observed SEDs, as we will show below.

It should be stressed here that the real circumstellar environment around R CrB is most likely much more complex than depicted by the above spherical shell model. Since dust ejection in RCB stars presumably occurs in clouds, not in a spherical shell, it is likely that the real distribution of material is clumpy or patchy without any clear inner boundary. At large distances, however, it is still plausible that a spherical shell model roughly represents the distribution of material, as long as dust clouds are ejected randomly and frequently. There may be some newly formed clouds inside the inner boundary defined in the above spherical model, if the star ejects dust clouds frequently, for example, every pulsational cycle (40 – 50 days). Without detailed knowledge about the dispersal processes of clouds, however, it is beyond the scope of this paper to construct a more detailed model. In the framework of our models, we consider only one newly formed dust cloud, as Fig. 3 illustrates.

The observed SEDs (see Sect. 5) demonstrate that the contribution of the central star is not negligible in the near-infrared. The flux from the central star was attenuated by a dust cloud in front of the star by $`\mathrm{\Delta }V1`$ on 1996 October 1 and $`\mathrm{\Delta }V4.6`$ on 1999 September 28. In order to estimate the attenuated flux from the star in the near-infrared, the extinction curve of the dust cloud is empirically derived from photometry in the optical, where emission from the dust shell is negligible. The effect of the obscuration due to a dust cloud is expressed as

$$f_{}(\lambda )=f_{}^{\mathrm{max}}(\lambda )\mathrm{exp}(\tau ^{\mathrm{cl}}(\lambda )),$$
(4)
where $`f_{}(\lambda )`$ and $`f_{}^{\mathrm{max}}(\lambda )`$ denote fluxes observed at any given time and at maximum light, respectively, and $`\tau ^{\mathrm{cl}}(\lambda )`$ is the optical depth of the dust cloud in front of the star. Note that the effect of interstellar extinction, though small for R CrB, cancels out, because it affects both $`f_{}(\lambda )`$ and $`f_{}^{\mathrm{max}}(\lambda )`$ by the same amount.

Fig. 5 shows the temporal change of the optical depth of a dust cloud during the 1983 minimum in the $`U`$, $`B`$, $`V`$, $`R`$, and $`I`$ bands. The optical depth in each band was derived using the photometric data in the decline and at maximum light obtained by Goncharova (goncharova92 (1992)). We also plot the extinction curves derived for the photometric data discussed in Sect. 3. The figure illustrates that the extinction is almost independent of the wavelength at the very bottom of the deep minimum ($`V=14.19`$). As the star starts its final recovery, the extinction curve starts to steepen. Pugach (pugach84 (1984)) proposed that the neutral extinction observed during the initial drop to minimum light can be explained by an optically very thick dust cloud whose coverage over the stellar disk varies with time. On the other hand, Hecht et al. (hecht84 (1984)) analyzed the change of the ratio of total to selective extinction ($`A_V/E(BV)`$) during decline events for R CrB and another RCB star, RY Sgr, and suggested that rather large glassy carbon particles with radii from 0.075 to 0.15 $`\mu `$m might be responsible for the nearly neutral extinction seen at the very bottom of minimum light. These sizes are significantly larger than the 0.01 $`\mu `$m we adopt for a single grain size, but it should be noted that we adopt 0.01 $`\mu `$m as the grain size in the optically thin dust shell, while the large grains mentioned above are claimed to be present in a newly formed optically thick dust cloud. For the temporal change of the extinction at the rise phase from minimum light, Hecht et al. (hecht84 (1984)) proposed shattering of large grains as a possible mechanism. Efimov (efimov90 (1990)) also proposed that the observed temporal variation of colors and polarization can be accounted for by the change of grain sizes, but assuming graphite grains instead of amorphous carbon. In any case, the spectral index of the dust cloud approaches 1, as the cloud disperses and becomes part of the optically thin dust shell. Therefore, the dust properties in the optically thin shell may presumably remain constant.

The lack of detailed knowledge about the properties of dust grains formed in R CrB and of their temporal variation forces us to adopt an empirical extinction law. As Fig. 5 shows, the extinction due to a dust cloud can be approximated by $`\tau ^{\mathrm{cl}}(\lambda )(1/\lambda )^p`$. At the bottom of minimum light, $`p`$ is approximately 0, and changes to $`1`$, as the star returns to maximum light. By the least square fit for the extinction curves shown in Fig. 5, we derive $`\mathrm{log}\tau ^{\mathrm{cl}}(\lambda )=0.520.47\mathrm{log}\lambda (\mu \text{m})`$ when the $`V`$ magnitude of the star is 10.61, while $`\mathrm{log}\tau ^{\mathrm{cl}}(\lambda )=0.180.84\mathrm{log}\lambda (\mu \text{m})`$ when the $`V`$ magnitude is 7.03. The extinction at longer wavelengths is then estimated by extrapolation. Obviously this procedure is not optimal, but with no alternatives at hand for disentangling the contributions of the central star and the dust shell for the data studied here, we are forced to adopt this method.

Assuming a blackbody for the flux of the central star, the flux density observed on the earth from the obscured central star is

$$f_{}(\lambda )=\left(\frac{R_{}}{D}\right)^2\pi B_\lambda (T_{\mathrm{eff}})\times \mathrm{exp}(\tau ^{\mathrm{cl}}(\lambda )\tau ^\mathrm{s}(\lambda )),$$
(5)
where $`\tau ^\mathrm{s}(\lambda )`$ is the optical depth of the dust shell along the radial direction.

The flux density observed toward R CrB is, therefore, the sum of the flux from the central star obscured by a dust cloud and slightly dimmed by the optically thin dust shell (equation (5)), thermal emission from the optically thin dust shell (equation (3)), and the thermal emission of a newly formed (therefore presumably still optically thick) dust cloud given by

$$f_{\mathrm{cl}}(\lambda )=\left(\frac{R_{\mathrm{cl}}}{D}\right)^2\pi B_\lambda (T_{\mathrm{cl}})\times \mathrm{exp}(\tau ^\mathrm{s}(\lambda )),$$
(6)
where $`R_{\mathrm{cl}}`$ and $`T_{\mathrm{cl}}`$ are the radius and the temperature of the newly formed dust cloud, respectively.

The intensity distribution for the central star and the dust shell can be written as

$`I_\lambda (b)`$ $`=`$ $`B_\lambda (T_{})\mathrm{exp}(\tau ^{\mathrm{cl}}(\lambda )\tau ^\mathrm{s}(\lambda ))\times \mathrm{circ}(b/R_{})`$ (7)
$`+`$ $`2m_\mathrm{d}\kappa _\lambda {\displaystyle _{z_{\mathrm{min}}}^{\sqrt{r_{\mathrm{out}}^2b^2}}}n(r)B_\lambda (T(r))𝑑z`$
$`+`$ $`B_\lambda (T_{\mathrm{cl}})\mathrm{exp}(\tau ^\mathrm{s}(\lambda ))\times \mathrm{circ}(b/R_{\mathrm{cl}}),`$
where $`b`$ is the impact parameter and $`z=\sqrt{r^2b^2}`$. The function $`\mathrm{circ}(b/R_{})`$ takes a value of 1 for $`|b|<R_{}`$ and 0 elsewhere. The lower limit of the integration is $`z_{\mathrm{min}}=0`$ for $`br_{\mathrm{in}}`$, while $`z_{\mathrm{min}}=\sqrt{r_{\mathrm{in}}^2b^2}`$ for $`b<r_{\mathrm{in}}`$. The visibility is calculated by taking the modulus of the Fourier transform of the intensity distribution.

## 5 Simultaneous fit of observed visibilities and SEDs

### 5.1 Model fitting without thermal emission from a dust cloud

We first try to fit the observed SEDs and visibilities using models without thermal emission from a dust cloud, namely, neglecting the term given by equation (6). We adopt $`\gamma =2`$, appropriate for a constant mass loss and expansion velocity. The sensitiveness of the near- and mid-infrared parts of emergent SEDs to $`T_{\mathrm{in}}`$ allows us to determine it by fitting the observed SEDs. The visibility functions, which directly reflect the spatial extent of the shell, enable us to examine the validity of the models described above.

Figs. 6 and 7 show the fit of the observed SEDs and visibility functions at the two epochs of our speckle observations. The SEDs and visibilities are calculated with the opacities of amorphous carbon obtained by Bussoletti et al. (bussoletti87 (1987)) (AC2 sample), Rouleau & Martin (rouleau91 (1991)) (AC1 sample), and Colangeli et al. (colangeli95 (1995)) (ACAR sample). The corresponding temperature distributions shown in Fig. 4 are used in the calculations. Figs. 6a and 6b demonstrate that the SED and visibility for the 1996 data are well reproduced with $`r_{\mathrm{in}}`$ = 76 – 90 $`R_{}`$ and $`T_{\mathrm{in}}`$ = 900 – 950 K. For a given opacity data set, the uncertainty of $`T_{\mathrm{in}}`$ is estimated to be $`\pm 100`$ K, translating to an uncertainty of $`r_{\mathrm{in}}`$ of $`\pm 20`$ $`R_{}`$. Taking the average of the $`T_{\mathrm{in}}`$ and $`r_{\mathrm{in}}`$ derived with three different opacities and adding the uncertainties resulting from the fitting, we derive $`T_{\mathrm{in}}`$ = $`920\pm 103`$ K and $`r_{\mathrm{in}}`$ = $`82\pm 23`$ $`R_{}`$. Fig. 6c illustrates the normalized intensity profile of the best fit model for the 1996 data. It consists of the central star and a ring-like structure characteristic of an optically thin shell. Note that the predicted visibilities become plateau-like at spatial frequencies $`20`$ cycles/arcsec in Fig. 6b. This plateau results from the unresolved central star. It should also be noted that the good match in the wavelength range shorter than 1 $`\mu `$m is simply due to the adoption of the empirical extinction for a dust cloud, as described in Sect. 4.

Figs. 7a and 7b reveal that the visibility observed at minimum light cannot be reproduced simultaneously with the observed SED. The observed SED can be well fitted using the models with $`T_{\mathrm{in}}`$ = 700 – 800 K and $`r_{\mathrm{in}}`$ = 116 – 172 $`R_{}`$, but the visibilities predicted from these models are too low as compared with the observation. Regarding this discrepancy, we first examine the approximations adopted in our models. The line-blanketing effect in the atmosphere of RCB stars is very prominent as Asplund et al. (asplund97 (1997)) show, but the use of line-blanketed atmospheres instead of the blackbody would lower $`T_{\mathrm{in}}`$ by only $`50`$ K. The effect of the uncertainty of the effective temperature is also minor. A decrease of the effective temperature by 500 K leads to a decrease of the dust temperature by $`60`$ K at $`100`$ $`R_{}`$. Combining these effects, the dust temperature in the shell can be by $`100`$ K lower than those used in the fitting above. We tried to fit the observed SED and visibility using such a temperature distribution, but it has turned out that the match to the observations is not much improved.

One concern is the estimation of the $`K`$-band flux of the central star obscured by the dust cloud. As described in Sect. 4, the extinction of the dust cloud in the $`K`$-band is derived by extrapolation from the region shortward of 1 $`\mu `$m. However, the actual extinction curve of the dust cloud could become steeper longward of 1 $`\mu `$m and approach the usual extinction curve of amorphous carbon, which is characterized by a spectral index of $`1.3`$ (Le Bertre lebertre97 (1997)). We find that the observed SED and visibility can be simultaneously fitted by adopting an extinction curve for the obscuring dust cloud such as shown in Fig. 8c. The extinction curve bends at 1.4 $`\mu `$m and the spectral index changes from 0.47 to 1.3. Figs. 8a and 8b reveal that the observed SED and visibility are fairly reproduced. It should be stressed that the physical understanding of such an extinction curve is still unclear, and that the exact shape of the extinction curve cannot be uniquely determined by the fitting presented here. However, the unusual properties of grains in the newly formed dust cloud may be responsible for the discrepancy found for the minimum light data.

### 5.2 Thermal emission from a newly formed dust cloud

We show an alternative to interpret the data obtained at minimum light. We propose that the inconsistency found for the 1999 data set may be attributed to the thermal emission of a newly formed optically thick dust cloud in front of the star. The newly formed dust cloud is supposedly still rather close to the star and its angular size is not yet large enough to be resolved with the 6 m telescope. This assumption may be reasonable, because the 1999 speckle observation was carried out in a decline when the star was obscured by $`\mathrm{\Delta }V=4.6`$. The size and the temperature of such a dust cloud are by no means obvious, but the temperature ($`T_{\mathrm{cl}}`$) should not exceed $`1500`$ K, otherwise the thermal emission of the cloud would be too prominent and would lead to a poorer match to the observed SED. We assume a circular cloud (projected onto the sky as seen from the star) with a radius of $`R_{\mathrm{cl}}`$ between 1 $`R_{}`$ and 10 $`R_{}`$. As shown in Fig. 3, $`R_{\mathrm{cl}}`$ denotes the radius of an individual cloud itself, not the distance between the cloud and the central star. As 1 $`R_{}`$ corresponds to 0.23 mas, even a cloud with a radius of 10 $`R_{}`$ can well be treated as a point source for our speckle interferometry. The thermal emission of the cloud is approximated by the blackbody radiation, as given by equation (6). Fig. 9c shows the normalized intensity profile of a model consisting of the central star, an optically thin dust shell, and a thermally emitting optically thick dust cloud, whose parameters will be derived below.

Figs. 9a and 9b show a simultaneous fit using models with thermal emission from a newly formed dust cloud. The observed SED and visibility are now fairly reproduced with models consisting of the central star, an optically thin dust shell with an inner radius of $`r_{\mathrm{in}}`$ = 150 – 193 $`R_{}`$ (average 172 $`R_{}`$) and $`T_{\mathrm{in}}`$ = 650 – 720 K (average 685 K), and an optically thick dust cloud with $`R_{\mathrm{cl}}`$ = 4.5 $`R_{}`$ and $`T_{\mathrm{cl}}`$ = 1200 K. The fit should be regarded as tentative. The uncertainties of $`T_{\mathrm{in}}`$ and $`r_{\mathrm{in}}`$ are estimated to be $`\pm 100`$ K and $`\pm 20`$ $`R_{}`$, respectively, while those of $`T_{\mathrm{cl}}`$ and $`R_{\mathrm{cl}}`$ are about $`\pm 100`$ K and $`\pm 1`$ $`R_{}`$, respectively.

The difference of $`r_{\mathrm{in}}`$ derived for the 1996 and 1999 data may indicate a variation of the dust ejection frequency. If the star experiences dust ejection less frequently, the inner boundary is expected to become larger, and the IR excess is expected to decrease. In fact, as mentioned in Sect. 3, the $`L`$ and $`M`$ magnitudes decreased monotonically by $`1`$ mag from 1996 to 1999. This is consistent with the trend found for $`r_{\mathrm{in}}`$ derived for the two epochs, but further observations are required to confirm this correlation.

Now we apply the thermally emitting cloud models to the data obtained near maximum light as well for the following reason. Fig. 1 shows that R CrB started its final recovery from the 1995-96 deep minimum about 6 months before our speckle observation on 1996 October 1. If dust ejection occurs rather frequently, for example, every pulsational cycle (40 – 50 days), we can expect that the star underwent dust ejection *out of* the line of sight during these 6 months. It is likely that newly formed dust clouds *out of* the line of sight existed on 1996 October 1, in addition to the dust cloud in front of the star responsible for an obscuration of $`\mathrm{\Delta }V1`$. The latter cloud may already be regarded as part of the optically thin dust shell, because $`\tau _K^{\mathrm{cl}}=0.34`$ is derived for this cloud by the extrapolation described in Sect. 4. Therefore, for near maximum light, we consider a model with the central star dimmed by $`\mathrm{\Delta }V1`$, thermal emission from an optically thin dust shell, and that from a newly formed optically thick dust cloud *out of* the line of sight. Fig. 10 shows the SED and visibility calculated with an optically thick dust cloud with a radius of $`R_{\mathrm{cl}}`$ = 4.5 $`R_{}`$ and $`T_{\mathrm{cl}}`$ = 1200 K, in addition to the dimmed central star and the optically thin dust shell, whose parameters are the same as shown in Fig. 6. The cloud is placed out of the line of sight, 20 $`R_{}`$ (5 mas) offset from the central star in the plane of the sky. Since the central star has already almost regained its brightness near maximum light, the effect of the newly formed dust cloud is minor on the SED and visibility. Therefore, the models including thermal emission from a newly formed dust cloud can also provide good matches to the SED and visibility observed near maximum light.

## 6 Concluding Remarks

Our 75 mas resolution speckle interferometric observations with the SAO 6 m telescope have spatially resolved the dust shell around R CrB for the first time. Neither the observation near maximum light nor at minimum light shows any clear deviation from circular symmetry.

In order to derive the size of the dust shell, we first considered models consisting of the central star and an optically thin dust shell, neglecting the thermal emission of a newly formed dust cloud. Simultaneous fits of the models to the observed SED and visibility have demonstrated that a model with the central star and an optically thin dust shell with density proportional to $`r^2`$ seems to be appropriate for R CrB near maximum light. The inner boundary is found to be $`82R_{}\text{(19 mas)}\pm 23`$ $`R_{}`$ with a temperature of $`920\pm 103`$ K.

This simple picture fails to simultaneously reproduce the SED and visibility observed at minimum light, which has led us to investigate models with thermal emission from a newly formed optically thick dust cloud whose angular size is not yet large enough to be spatially resolvable. The SED and visibility obtained at minimum light were shown to be well fitted with such models. The presence of a newly formed dust cloud as hot as 1200 K with a radius of 4 – 5 $`R_{}`$ is inferred, together with an optically thin dust shell with $`r_{\mathrm{in}}170`$ $`R_{}`$ and $`T_{\mathrm{in}}690`$ K. Furthermore, the SED and visibility obtained near maximum light were shown to be fitted also using a model with a newly formed dust cloud out of  the line of sight. However, we have also discussed that the discrepancy found for the minimum light data may be attributed to the unusual extinction curve of the obscuring dust cloud. Observations during the very bottom of a deep minimum, when the contribution of the central star is truly negligible, are crucial for investigating the dust shell and dust clouds.

###### Acknowledgements.

We have used, and acknowledge with thanks, data from the AAVSO International Database, based on observations submitted to the AAVSO by variable star observers worldwide. N.R.I. acknowledges the support of the Long-Term Cooperation program of the Alexander von Humboldt Foundation.


format: text

text:
# 1 INTRODUCTION

## 1 INTRODUCTION

When phase transitions are studied in statistical physics, the basic assumption is usually that the dimensions of the macroscopic system are very large compared with that of the constituting elements. However, there are many phenomena in finite systems which resemble phase transitions, for instance, in the physics of clusters of atoms . The main question is how the observed effects in small systems can be related to true phase transitions in macroscopic (or infinite) systems. Attempts in this direction include the exploration of the density of complex zeros of the canonical partition function for finite and small systems by Borrmann et. al.  and the linear behavior for this limiting density . The approach suggested by Janke and Kenna (JK) presents also a new scaling relation to identify the order and strength of a transition from the behavior of small systems .

In this paper we try to evaluate the usefulness and validity of the Borrmann and the JK approaches for the investigation of structural transitions in biomolecules, another important example for “phase transitions” in small systems. For this purpose we consider the helix-coil transition in polyalanine and test whether the two approaches allow a characterization of that transition. In addition, we also investigate a second molecule Ala10-Gly5-Ala10 in the context of Borrmann approach.

## 2 ANALYSIS OF PARTITION FUNCTION ZEROS IN SMALL SYSTEMS

In the canonical ensemble a system is completely described by its partition function

$$Z(\beta )=\underset{E}{}n(E)\mathrm{exp}(\beta E).$$
(1)
Introducing the variable $`u=\mathrm{exp}(k\beta )`$ with conveniently defined constant $`k`$ allows to write the partition function for discrete energy models as a polynominal:

$$Z(u)=\underset{E}{}n(E)u^E.$$
(2)
We label the complex zeros in order of increasing imaginary part as $`u_j=\mathrm{exp}(k\beta _j),j=1,2,\mathrm{}`$.

As established by Yang and Lee  and later by Fisher , the statistical theory of phase transitions can be described by the distribution of complex partition function zeros. In the case of a temperature driven phase transition, we expect the zeros (or at least the ones close to the real axis) condense for large enough system size $`L`$ on a single line

$$u_j=u_c+r_j\mathrm{exp}(i\phi ).$$
(3)
As the system size $`L`$ increases, these zeros will move towards the positive real $`u`$-axis and for large $`L`$ the corresponding value is $`\beta _c`$, the inverse of the physical critical temperature $`T_c`$. Crucial information on phase transitions can be obtained from the way in which the first zero approaches the real $`u`$-axis. However, such an analysis depends on the extrapolation towards the infinite large system and does not allow characterization of the situation in small systems.

One possible extension of the above ideas to “phase transitions” in small systems is the classification scheme by Borrmann et al. . This classification scheme computes density of zeros taking into account numerical estimates of the first four complex zeros. Now, writing the complex zeros as $`u_k=\mathrm{Re}(u_k)+i\tau _k`$, where $`\tau _k`$ stands for $`\mathrm{Im}(u_k)`$, the assumed distribution of zeros on a straight line allows to define two parameters $`\alpha _u`$ and $`\gamma _u`$:

$$\alpha _u=\frac{\mathrm{ln}\varphi (\tau _3)\mathrm{ln}\varphi (\tau _2)}{\mathrm{ln}\tau _3\mathrm{ln}\tau _2},\mathrm{where}\varphi (\tau _k)=\frac{1}{2}\left(\frac{1}{|u_ku_{k1}|}+\frac{1}{|u_{k+1}u_k|}\right),$$
(4)
with $`k`$ labelling the first zeros, and

$$\gamma _u=(\mathrm{Re}(u_2)\mathrm{Re}(u_1))/(\tau _2\tau _1).$$
(5)
Note, that our notation differs from that in Ref.  in that we define the discrete line density $`\varphi `$ in function of the $`u`$-zeros instead of the inverse temperature $`\beta `$. Following the classification scheme by Grossmann and Rosenhauer , phase transitions can now be classified according to the values of these two parameters: for $`\alpha _u=0`$ and $`\gamma _u=0`$ one has a phase transition of first order, it is of second order if $`0<\alpha _u<1`$ and arbitrary $`\gamma _u`$, and for $`\alpha _u>1`$ and arbitrary $`\gamma _u`$ one has a higher order transition. Of course, as expected for a phase transition, the imaginary part of the leading zero $`(\tau _1)`$ should move towards the real axis to obtain in the thermodynamic limit a real temperature. However, for small systems with finite $`\tau _1`$, we may have $`\alpha <0`$ , which translates the condition of first order transition to $`\alpha _u0`$ and $`\gamma _u=0`$. This classification of phase transitions has been tested for the finite Bose-Einstein condensates in a harmonic trap, small magnetic clusters and nuclear multifragmentation .

Another extension of partition function zero analysis to small systems is the approach by Janke and Kenna (JK) , which uses that the theoretical average cumulative density of zeros

$$G_L(r_j)=\frac{2j1}{2L^d}$$
(6)
can be written in the thermodynamic limit and for a first order transition as

$$G_{\mathrm{}}(r)=g_{\mathrm{}}(0)r+ar^{w+1}+\mathrm{}.$$
(7)
Here, the slope at the origin is related to the latent heat, $`\mathrm{\Delta }eg_{\mathrm{}}(0)`$. Equations (3) and (7) imply that the distance $`r_j`$ of a zero from its critical point can be written for large enough lattice sizes as $`\mathrm{Im}u_j(L)`$ since $`\mathrm{Re}u_j(L)u_c`$. Hence, in this limit Eq. (6) and (7) lead to the following scaling relation for the cumulative density of zeros as an equation in $`j`$ and $`L`$,

$$\frac{2j1}{2L^d}=a_1(\mathrm{Im}u_j(L))^{a_2}+a_3.$$
(8)
A necessary condition for the existence of a phase transition is that $`a_3`$ is compatible with zero, else it would indicate that the system is in a well-defined phase. The values of the constants $`a_1`$ and $`a_2`$ then characterize the phase transition. For instance, for first order transitions should the constant $`a_2`$ take values $`a_21`$ for small $`r`$, and in this case the slope of this equation is related to the latent heat through the relation

$$\mathrm{\Delta }e=ku_c\mathrm{\hspace{0.17em}2}\pi a_1,$$
(9)
with $`u_c=\mathrm{exp}(k\beta _c)`$. On the other hand, a value of $`a_2`$ larger than $`1`$ indicates a second order transition whose specific heat exponent is given by $`\alpha =2a_2`$.

The above approach was originally developed and tested for systems with well defined first order phase transitions such as the 2D 10-state Potts and 3D 3-state Potts model. The obtained results agree with previous work from numerical simulations and partition function zeros analysis of systems up to $`L=64`$  ($`L=36`$ for the 3D case ). However, some difficulties in the identification of the order of the transition appear when applied to a model with a very weak first order transition .

Here we explore further these two approaches for a non-trivial example of structural transitions in biomolecules, namely the helix-coil transition in polyalanine. The characteristics of this so-called helix-coil transition have been studied extensively . For instance, evidence was presented  that polyalanine exhibits a phase transition between the ordered helical state and the disordered random-coil state when interactions between all atoms in the molecule are taken into account. These previous results, obtained independently by different methods, are used here as benchmarks for the test of the usefulness and to point out the difficulties in applying the classification scheme of Borrmann and Janke/Kenna for research of structural transitions in biomolecules.

Our investigation of the helix-coil transition for polyalanine is based on a detailed, all-atom representation of that homopolymer. Since one can avoid the complications of electrostatic and hydrogen-bond interactions of side chains with the solvent for alanine (a nonpolar amino acid), explicit solvent molecules were neglected. The interaction between the atoms was described by a standard force field, ECEPP/2,. Chains of up to $`N=30`$ monomers were considered, and our results rely on multicanonical simulations  of $`N_{sw}`$ Monte Carlo sweeps starting from a random initial conformation, i.e. without introducing any bias. Our statistics consists of $`N_{sw}`$=400,000, 500,000, 1,000,000, and 3,000,000 sweeps for $`N=10`$, 15, 20, and 30, respectively. Measurements were taken every fourth Monte Carlo sweep. Additional 40,000 sweeps ($`N=10`$) to 500,000 sweeps ($`N=30`$) were needed for the weight factor calculations by the iterative procedure described first in Ref. . In contrast to our first calculation of complex zeros presented in Ref. , where we divided the energy range into intervals of lengths 0.5 kcal/mol in order to make Eq. (2) a polynomial in the variable $`u=e^{\beta /2}`$, we avoided any approximation scheme in the present work. This because the above approximation works very well for the first zero, but not for the next ones. Since we need high precision estimates also for the next zeros we applied the scan method (see Ref.  and references therein).

In Table 1 we present our first four partition function zeros, although the fourth one is less reliable due to the presence of fake zeros . Using these zeros we first calculated the parameters $`\alpha _u`$ and $`\gamma _u`$ which characterize in the Borrmann approach phase transitions in small systems. As one can see in Table 2 the so obtained values for polyalanine are characterized by large fluctuations. It seems that the median of the $`\alpha _u`$ values is $`\alpha _u=0`$ which would indicate a first order transition. However, our data are not good enough to draw such a conclusion. Our error estimates were obtained by means of the Bootstrap method . For this reason, we tried instead the JK scaling relations. Table 3 lists the parameter $`a_3(N)`$ of Eq. (8). Here, the average cumulative density of zeros is replaced by

$$G_N(r_j)=\frac{2j1}{2N},$$
(10)
where we have translated the linear length $`L`$ as $`N^{1/d}`$ . Therefore all finite size scaling relations can be written in terms of the number of monomers $`N`$. In Fig. 1 we show the cumulative distribution of zeros in a log-log plot for chain lengths $`N=10,15,20\mathrm{and}30`$.

The values $`a_3`$ are compatible with zero for chains of all length indicating that we have indeed a phase transition. In order to evaluate the kind of transition we also calculate the parameters $`a_1(N)`$ and $`a_2(N)`$ which we also summarize in Table 3. Our analysis procedure differs from JK in that we are calculating $`a_2`$ for a fixed $`N`$before performing a finite-size scaling analysis of that quantity, i.e. we account for the dependence of the average cumulative density on the system size. The parameter $`a_2(N)`$ decreases with system size and a log-log plot of this quantity as a function of chain length, as shown in Fig. 2, suggests a scaling relation

$$a_2(N)=a_2+bN^c.$$
(11)
A numerical fit of our data to this function leads to a value of $`a_2=1.31(4)`$ with goodness of fit $`Q=0.95`$. Using $`\alpha =2a_2`$ we find $`\alpha =0.70(4)`$ which is barely compatible with our previous value of $`\alpha =0.86(10)`$ in Ref. , obtained from the maximum of specific heat. A fit of all four chain lengths can also not exclude a value $`a_2=1`$ (i.e. a first order phase transition) since we can find acceptable fits with $`Q>0.55`$ in the range $`0.92<a_2<1.44`$. Especially, a close examination of Fig. 1 and Fig. 2 shows that the $`N=30`$ data point exhibits a considerable deviation from the trend suggested by the smaller chain lengths. Since the $`N=30`$ data are the least reliable, we also evaluated Eq. (11) omitting the $`N=30`$ chain. This leads to a value of $`a_2=1.16(1)`$ and a critical exponent $`\alpha =0.84(1)`$ which is now compatible with our previous result $`\alpha =0.86(10)`$.

While the JK approach is able to reproduce results for polyalanine obtained in previous work , it does not allow to establish the order of the helix-coil transition from simulation of small chains. Our results for the parameter $`a_2`$ seem to favor a second order transition, but are hampered by large errors (such that a first order transition cannot be excluded) and disputed in Ref.  were indications for a finite latent heat were found.

Borrmann’s approach led to even less decisive results for polyalanine since the fluctuation were too large to draw firm conclusions. However, the advantage of this approach is that it does not require a finite-size scaling. We have therefore used this method to investigate a second molecule, the slightly more complicated Ala10-Gly5-Ala10. A detailed study of this molecule will be published elsewhere .

We describe the interactions between the atoms again by the ECEPP/2  force field (as implemented in the program package SMMP ). Our results rely on multicanonical simulations  of $`4,000,000`$ Monte Carlo sweeps starting from a random initial conformation, i.e. without introducing any bias. Measurements were taken every tenth Monte Carlo sweep. Additional 500,000 sweeps were needed for the calculations of the weight factor.

Analyzing the partition function zeros for this peptide, we find a “phase transition” at two temperatures, each being characterized by a line of complex zeros. The corresponding first four zeros for each characteristic line are listed in Table 4. Both temperatures also correspond with peaks in the specific heat which is displayed in Fig. 3 as a function of temperature. In Ref.  it is shown that the higher temperature $`T_1=480`$ K marks the helix-coil transition, i.e. the temperature where secondary structure elements ($`\alpha `$-helices) are formed. At the lower temperature $`T_2=260`$ K, the peptide then folds into its native structure, a U-turn-like bundle of two (antiparallel) $`alpha`$-helices connected by a turn. Here, we use the Borrmann approach to study these two transitions in more detail. Using the zeros from Table 4 we calculated the parameters $`\alpha _u`$ and $`\gamma _u`$ which characterize in the Borrmann approach phase transitions in small systems. For the first transition, at $`T=480`$ K, which corresponds to the above studied helix-coil transition in polyalanine, we find $`\alpha _u=1.46`$ and $`\gamma _u=0.49`$. While the negative value of $`\alpha _u`$ indicates a first order transition, our results do again not allow us to draw a firm conclusion on the nature of the transition since the values of the two parameters $`\alpha _u`$ and $`\gamma _u`$ are subject to large fluctuations when we vary the number of sweeps considered in the calculation of the zeros. This is not surprising since we were also not able to establish clearly the order of the helix-coil transition in in the previous example of polyalanine. However, our data are more decisive in the case of the second transition, at $`T=260`$ K, which marks the compactification and folding of the peptide. Here we find $`\alpha _u=0.38`$ and $`\gamma _u=1.07`$. These values indicate a second-order transition which is consistent with what one would expect for a transition between extended and compact structures. We hope that by doubling the number of Monte Carlo sweeps we will be able to verify this results and also obtain reliable values for the Borrmann parameter characterizing the helix-coil transition at $`T=480`$ K. These simulations are now under way.

## 3 Conclusion

Let us summarize our results. We have evaluated two recently proposed schemes for characterizing phase transitions in small systems. Simulating polyalanine molecules of chain lengths up $`N=30`$ residues by multicanonical Monte Carlo, we calculated the partition function zeros of these molecules. Analyzing these zeros by the JK approach we were able to reproduce for polyalanine results obtained in previous work  while the Borrmann approach led here to inconclusive results. However, the later method seems more appropriated to study structural transitions in proteins since it does not require a finite size scaling. When applied to Ala10-Gly5-Ala10, the Borrmann approach allowed us to characterize structural transitions in this molecule. However, our results are limited again by the precision with which one can calculate the first four partition function zeros. While these numerical limitations restrict the applicability of the Borrmann and the JK approaches, our results demonstrate that both methods can be useful tools for the study of “phase transitions” in biomolecules.

Acknowledgements
U.H.E. Hansmann gratefully acknowledges support by research grant of the National Science Foundation (CHE-9981874), and N.A. Alves support by CNPq (Brazil).

Figure Captions:

Distribution of zeros for the polyalanine model.

The parameter $`a_2(N)`$ for polyalanine molecules of length $`N`$ in a log-log plot.

The specific heat $`C(T)`$ as a function of temperature for Ala10-Gly5-Ala10.

Table Captions:

Partition function zeros for polyalanine.

Estimates for the parameters $`\alpha _u`$ and $`\gamma _u`$ for polyalanine.

Estimates for the JK parameters for polyalanine.

Partition function zeros corresponding to the two transitions observed for Ala10-Gly5-Ala10.