text:
# On the accretion geometry of Cyg X-1 in the low/hard state

## 1 Introduction

The black hole binary (BHB) systems in their low/hard state have spectra dominated by a power law which rolls over at $`200`$ keV. This continuum is generally well fit by thermal Comptonization models, in which low energy photons from the accretion disc are upscattered by energetic electrons (e.g. the review by Zdziarski 2000). To get these hard X-rays, a large fraction of the gravitational energy released by accretion must be dissipated in an optically thin environment i.e. not in the disc itself. However, there is no consensus on how this happens, or on the geometry of this hot region. There are currently two main models, one in which the hot electrons are confined in magnetic flares above a disc which extends down to the last stable orbit around the black hole (magnetic flares), and one in which the electrons form an quasi-spherical hot flow, replacing the inner disc (truncated disc).

The magnetic flare model is motivated by the discovery that the physical mechanism for the disc viscosity is a magneto-hydrodynamic dynamo (e.g. the review by Balbus & Hawley 2002). Buoyancy could cause the magnetic field loops to rise up to the surface of the disc, so they can reconnect in regions of fairly low particle density, forming a patchy corona. Numerical simulations (although these are highly incomplete as in general the simulated discs are not radiative) do show this happening (Hawley 2000), but they do not yet carry enough power to reproduce the observed low/hard state (Miller & Stone 2000).

The truncated disc model has its physical basis in the accretion flow equations. The standard Shakura-Sunyaev disc solution assumes that the accreting material is at one temperature (protons and electrons thermalize) and that the accretion energy released by viscosity is radiated efficiently. At low mass accretion rates neither of these are necessarily true. The thermalization timescale between the electrons and protons can be long, so the flow is intrinsically a two temperature plasma. Where the electrons radiate most of the gravitational energy through Comptonization of photons from the outer disc then the hot inner flow is given by Shapiro, Lightman & Eardley (1976). Alternatively, if the protons carry most of the accretion energy into the black hole then this forms the advection dominated accretion flows (Narayan & Yi 1995). These are related, as in general both advection and radiative cooling are important for a hot accretion flow (Chen et al. 1995; Zdziarski 1998).

These two models of the accretion flow have very different geometries. A potential way to test the geometry is with X-ray reflection. The amount of reflection scales with the solid angle subtended by the optically thick material while the relativistic smearing of the atomic features shows how far this material extends into the gravitational potential (Fabian et al. 2000). With magnetic flares, the disc subtends a solid angle of $`2\pi `$ as seen from the X-ray source, and extends down to the last stable orbit. Its reflected spectrum should be large and strongly smeared. Conversely, a truncated disc illuminated by an inner hot flow subtends a solid angle $`2\pi `$, and the reflected spectrum is small and is only weakly smeared by relativistic effects.

The BHB spectra in the low/hard state show overwhelmingly that the solid angle is significantly less than $`2\pi `$, and that the smearing is less than expected for a disc extending down to the last stable orbit ($`\dot{\mathrm{Z}}`$ycki Done & Smith 1997; 1998; 1999; Gierliński et al. 1997; Done & $`\dot{\mathrm{Z}}`$ycki 1999; Zdziarski et al. 1999; Gilfanov, Churazov & Revnivtsev 1999; 2000). While this is clearly consistent with the idea that the disc is truncated in the low/hard state, the magnetic flare models can be retrieved in several ways.

Firstly, magnetic reconnection on the Sun is known to produce an outflow in the coronal mass ejection events. In the extreme conditions close to the black hole it is possible that this outflow velocity could be large so that the hard X–ray radiation is beamed away from the inner disc (Beloborodov 1999).

An alternative explanation for the lack of reflection and smearing is that the inner disc or top layer of the inner disc is completely ionized. There are then no atomic features, and the disc reflection is unobservable in the 2–20 keV range as it appears instead to be part of the power law continuum (Ross & Fabian 1993; Ross, Fabian & Young 1999). However, these models with passive illumination of the disc require a fairly sharp transition between the extreme ionization and mainly neutral material (Done & $`\dot{\mathrm{Z}}`$ycki 1999; Done, Madejski & $`\dot{\mathrm{Z}}`$ycki 2000; Young et al. 2001). Such a transition can be produced as the disc responds to the intense X–ray illumination. There is a thermal ionization instability which affects X–ray illuminated material in pressure balance, which can lead to a hot, extremely ionized skin forming on top of the rest of the cooler, denser, mainly neutral disc material (Field 1965; Krolik, McKee & Tarter 1981; Kallman & White 1989; Ko & Kallman 1994; Ró$`\dot{\mathrm{z}}`$ańska & Czerny 1996; Nayakshin, Kazanas & Kallman 2000, hereafter NKK; Ró$`\dot{\mathrm{z}}`$ańska & Czerny 2000; Nayakshin & Kallman 2001; Ballantyne, Ross & Fabian 2001). Such X-ray illuminated disc models can fit the 2–20 keV data from BHB with reflection from a disc which subtends a solid angle of $`2\pi `$ and extends down to the last stable orbit around a black hole (Done & Nayakshin 2001b).

While the truncated disc and magnetic flare/X-ray illuminated disc models are indistinguishable with current data in the 2–20 keV range, they are very different at higher energies. In the magnetic flare models, a lot of the 2–10 keV ’continuum’ is actually ionized reflection, so the true continuum level is lower than in the truncated disc models, in which reflection is small. At higher energies, where reflection is negligible, the lower continuum level for the magnetic flare models leads to a smaller predicted flux in the 100-200 keV range than for the truncated discs (Done & Nayakshin 2001a). Here we fit both truncated disc and magnetic flare/X-ray illuminated disc models to the 3–200 keV PCA/HEXTE spectrum of the hard/low state spectra of Cyg X-1. We show that the truncated disc models provide an excellent fit to these data, but that the magnetic flare/X-ray illuminated disc does not, as it dramatically underpredicts the 100–200 keV spectrum.

Similar conclusions were independently reached by Maccarone & Coppi (2002) from fits to the Cyg X-1 broad-band spectrum. However, they approximated the reflection from magnetic flares by highly ionized, single zone reflection models rather than the full complex ionization reflection models used here. Here we are able to show explicitally that the complex ionization magnetic flare models do not fit the high energy spectrum, while the truncated disc models do.

## 2 Models

We use the complex ionization reflection code, xion, described in NKK which computes the self-consistent vertical ionization structure of an X-ray illuminated disc at a given radius. The reflected and diffuse (line and recombination continua) emission from this are smeared by the relativistic effects expected from a disc, and then these spectra are summed over all radii to get the total disc reflected emission for any given source/disc geometry. These models have been tabulated for use in xspec for a central sphere/truncated disc geometry, assuming that the hot inner flow has emissivity $`r^3`$, and starts at the radius where the optically thick disc truncates (Nayakshin et al. 2002, in preparation). They have also been tabulated for a magnetic flare geometry, which results in much stronger illumination of the disc and so gives a deeper ionized skin (Nayakshin 2000). With magnetic flares, the luminosity of the flares is assumed to scale as the local disc flux, so that $`f_x/f_{disc}`$ remains constant with radius. These tabulated models are hereafter referred to as xion-disc (truncated disc) and xion-flares (magnetic flares). The free parameters are the mass accretion rate (this determines the unilluminated disc structure i.e. the starting density as a function of height), $`f_x/f_{disc}`$ which determines the depth of the ionized skin, and the shape of the illuminating spectrum ($`\mathrm{\Gamma }`$ and $`E_{cut}`$ assuming an exponentially cutoff power law). For the truncated disc models there is also the inner radius of the accretion disc, $`R_{in}`$, which controls the amount of relativistic smearing of the iron line features.

The xion code includes Compton up and downscattering by electrons in the disc. Compton upscattering can be very important in smearing the iron features in ionized discs (Ross et al. 1999), while Compton downscattering determines the high energy rollover in the reflected spectrum (Lightman & White 1988). However, the approximations used for Compton scattering in the xion code become progressively less accurate at higher energies (NKK) The Appendix discusses how we merge the more accurate high energy Compton downscattering of the standard xspec pexriv based reflection models (Magdziarz & Zdziarski 1995) with the xion results.

## 3 Spectral fitting

We used the hardest spectra seen from multiple observations of the Cyg X-1 hard/low state (Gilfanov et al. 1999). We extracted the data using the Rex data analysis script with the bright source background from the top layer of PCA detectors 0 and 1. Previous work has shown that this configuration gives a good fit to the Crab data with 0.5 per cent systematic error for RXTE Epoch 3 data (Wilson & Done 2001). We extracted the simultaneous HEXTE data from cluster 0. In fits with both the PCA and HEXTE data we allow a normalisation offset between the two instruments to take into account the cross-calibration uncertainties. We use xspec version 11 (Arnaud et al 1996) and quote all error bars as $`\mathrm{\Delta }\chi ^2=2.7`$, corresponding to 90 per cent confidence limits on 1 parameter, and fix the absorbing column to $`N_H=6\times 10^{21}`$ cm<sup>-2</sup> (Bałucińska-Church et al. 1995).

First we fit the PCA data alone with the old, single ionization parameter reflection models so that we can compare the new reflection models with previous fits. We use the rel-repr model of $`\dot{\mathrm{Z}}`$ycki et al. (1998), which is based on the publically available pexriv reflection but includes the self consistently produced iron line emission and relativistic smearing. We assume a power law continuum model, and also include neutral, unsmeared reflection which can arise from the companion star or outer accretion disc (Ebisawa et al 1996). This gives a $`\chi _\nu ^2=30/39`$ and the amount of reflection from the accretion disc is low, with $`\mathrm{\Omega }/2\pi =0.13\pm 0.08`$, implying a truncated disc.

We then replace the single ionization accretion disc reflection with the xion-disc model. The simplest version of this truncated disc geometry assumes the hot inner source is spherical, and starts where the cool disc truncates. This results in a model in which $`\mathrm{\Omega }/2\pi `$ is fixed at $`0.3`$, rather larger than inferred from the simple reflection models above, and so gives a somewhat larger $`\chi _\nu ^2=53/40`$. If instead we allow the amount of reflection to be smaller than this (e.g. if the hot source has a flattened rather than spherical geometry) then we get $`\chi _\nu ^2=43/39`$ for $`f_x/f_d=6`$. This is a worse fit by $`\chi ^2`$ of $`13`$ than that for the single ionization zone reflection models (rel-repr, above). The two models give different fits as even a weakly illuminated disc produces a low optical depth ionized skin. This reflects some small fraction of the incident continuum at low energies (see Fig. A1 in the appendix), changing the curvature of the reflected continuum.

We then use xion-flare to model the magnetic flare geometry in which $`\mathrm{\Omega }/2\pi `$ is fixed at 1 and the inner disc extends down to 3 Schwarzchild radii. This gives $`\chi _\nu ^2=47/41`$ (both the solid angle of reflection and the inner disc radius are fixed, leading to 2 fewer degrees of freedom). This confirms the results of Done & Nayakshin (2001b) that magnetic flares can give as good a fit to the data below 20 keV as a truncated disc. The PCA data are consistent with the presence of a large amount of ionized reflection, i.e. with a disc which subtends a large solid angle to the X-ray source and extends down to the last stable orbit.

However, these two very different reflection models should give very different spectra at higher energies. Reflection must rollover at high energies due to Compton downscattering in the disc. Fig. 1 shows the truncated disc (upper panel) and magnetic flare (lower panel) PCA fits extrapolated to the HEXTE data. The truncated disc model has little reflection, so there is little rollover in the total spectrum. By contrast, the magnetic flare models have a lot of highly ionized reflection, so the contribution of the Compton downscattering rollover is important. Plainly there is no such rollover in the HEXTE data, and the magnetic flare description of the PCA data strongly underestimates the higher energy spectrum. This is made even more marked by the fact that a power law is an overestimate of the high energy continuum: the true Comptonized continuum has a thermal cutoff at 150-200 keV (Gierliński et al. 1997). By contrast, the truncated disc models can qualitatively describe the high energy data.

While this argues strongly against a large ionized reflection component in the data, true spectral fitting is required to show conclusively that this is the case. We do joint fits of the PCA and HEXTE data, using the Comptonization model, compps, of Poutanen & Svensson (1996) rather than a power law for the continuum. With reflection from a spherical source within a truncated disc (and neutral, unsmeared reflection from the companion star) then this gives $`\chi _\nu ^2=114/77`$, while allowing the reflection amount to vary (as in a flattened source) gives $`\chi _\nu ^2=86/76`$. We show the unfolded spectrum and residuals from this model in Fig. 2, where the continuum parameters are $`\tau 0.75`$ and $`kT_e=200`$ keV. By contrast, using reflection from magnetic flares gives a very poor fit, with $`\chi _\nu ^2=259/78`$ (Fig. 3) for a continuum with $`\tau 0.75`$ and $`kT_e210`$ keV.

All these fits assumed isotropic seed photons for the Compton scattering. However, this is unlikely to be the case since the seed photons come from the disc. While there is some range of incident angle for the truncated disc geometry, the magnetic flares are illuminated preferentially from below. For a Comptonizing plasma with optical depth of order unity and temperature of $`200`$ keV (Sunyaev & Titarchuk 1980; Gierliński et al. 1997) then this seed photon anisotropy can have a marked effect on the spectrum (Haardt & Maraschi 1993; Poutanen & Svensson 1996). We use the compps code in a slab geometry, where the hot plasma is illuminated from below (Poutanen & Svensson 1996) as an approximation to the magnetic flare geometry. We refit the data with the magnetic flare reflection model using this new continuum. This gives an even worse fit than before, with $`\chi ^2=531/78`$, where the continuum has $`\tau 1`$ and $`kT_e=85`$ keV. The reason the fit is so poor is that the seed photon anisotropy leads to a break in the spectrum, such that the low energy spectrum is harder than the high energy spectrum (Haardt & Maraschi 1993; Haardt et al. 1993). The low energy spectrum then predicts even less high energy emission than for isotropic seed photons. Fig. 4 shows the unfolded spectrum and residuals from this model.

## 4 Discussion

The data clearly show that the 2–20 keV spectrum from the low/hard state of Cyg X-1 does not contain a large fraction of highly ionized reflection. This rules out models which have static magnetic flares above an untruncated disc unless the flares have a spectrum which is much harder than that predicted by a single temperature Comptonization model.

One way to get some spectral hardening which is expected in a magnetic flare geometry but is neglected in our modelling is for the flares to comptonize some fraction of the reflected photons. The reflected photons are from the disc, so are intercepted by the hot electrons in the flares in the same way as the soft seed photons, and are Compton upscattered to form a hard continuum. With flares covering most of the disc then this can make a $`50100`$ per cent increase to the flux at 200 keV (Petrucci et al. 2001). A covering fraction of unity for the flares (slab corona geometry) is normally ruled out for the low/hard spectra as this produces many non-reflected, thermalized photons from the hard X-ray illumination. These are also Comptonized by the corona, leading to spectral indices which are too soft to explain the low/hard state (Pietrini & Krolik 1995; Stern et al. 1995; Zdziarski et al. 1998). However, a covering fraction of $`0.5`$ might give enough Compton scattering of reflection to produce the required $`50`$ per cent excess emission at 200 keV, while also allowing enough seed photons to escape to produce the required hard continuum.

The problems that the magnetic flare models have in matching the high energy flux are exacerbated by the anisotropy break which should be present in the continuum in this assumed geometry, but which has never been convincingly observed (Gierliński et al. 1999). If the magnetic flares are to fit the high energy data then the anisotropy break must be hidden by having a multiple temperature Comptonised continuum. At some level there must be a distribution of electron temperatures: it is almost inconceivable that a single temperture distribution can be maintained, especially as the sources are variable so the flare spectra probably evolve with time (e.g. Poutanen & Fabian 1999). A multiple temperature continuum gives another way to boost the high energy flux so that we do not necessarily need to strongly Comptonise the reflected continuum. However, it still seems somewhat contrived that a combined complex continuum plus complex ionization reflection spectrum should so precisely mimic a simple single temperature continuum, truncated disc reflection.

By contrast, a truncated disc/hot inner flow geometry at low mass accretion rates is compatible with the observed hard continuum, lack of anisotropy break, low amount of reflection and relativistic smearing. It can also explain the low temperature and luminosity of the direct emission from the disc (e.g. Esin et al. 2000). This geometry can also give a qualitative explanation for a range of observed correlations if the truncation radius decreases with increasing (average) mass accretion rate. The disc penetrates further into the hot flow, increasing the seed photon flux intercepted by the hot inner flow, leading to a softer continuum spectra. This changing geometry gives a larger solid angle subtended by the disc, leading to an increasing amount of reflection (Poutanen Krolik & Ryde 1997; Zdziarski et al. 1999; Gilfanov et al. 1999; 2000), and relativistic smearing ($`\dot{\mathrm{Z}}`$ycki et al. 1999; Gilfanov et al. 2000; Lubiński & Zdziarski 2001). The variability power spectra are also affected as they contain characteristic frequencies which are most probably linked to the inner edge of the disc, so this can explain the correlated increase in break and quasi-periodic oscillation frequency (e.g. the review by van der Klis 2000; Churazov, Gilfanov & Revnivtsev 2001). Lastly, the collapse of an inner hot flow when it becomes optically thick gives a physical mechanism for the state transition (Esin McClintock & Narayan 1997).

Thus, if one were to choose between relatively straightfoward models, then the truncated disk is definitely favored by our analysis while the magnetic flare model is ruled out. However, the straightforward solutions may be too simple to describe the complexity of accretion disk structure near the black hole. The role and magnitude of secondary effects (comptonization of the reflection component; multi-temperature nature of flares) not taken into account in our modelling needs to be clarified in the future with detailed calculations.

## Appendix A A comparison of xion and pexriv reflection

The main objective of the xion code is to accurately calculate the vertical ionisation structure of the X-ray illuminated disc, so the atomic physics/photoionisation calculations are treated in great depth while compton scattering is approximated in the same manner as in Ross & Fabian (1993). These approximations are adequate at lower energies, where the energy shifts from compton scattering are small, but become progressively worse at higher energies where downscattering is large. By contrast, the standard pexriv reflection model in XSPEC treats compton scattering very carefully (Magdziarz & Zdziarski 1995), but has very crude ionisation balance and does not include any vertical structure (Done et al. 1992).

The two codes should give comparable results when the vertical structure is unimportant i.e. when the X-ray illumination is very small. Fig. A1 shows a comparison of the two codes with xion-flares at $`f_x/f_d=0.02`$ and $`\dot{m}=0.001`$ (the lowest tabulated values) for a magnetic flare geometry at inclination of $`30^{}`$. The dotted line shows the comparable pexriv results ($`\mathrm{\Omega }/2\pi =1`$, $`\xi =0`$, inclination of $`30^{}`$). The continuum for both is an illuminating power law of index $`\mathrm{\Gamma }=1.7`$ with exponential rollover at 300 keV. The differences are obvious. At low energies these are expected as the xion code has a residual ionized skin in even a weakly illuminated disc. One would have to extend the grid in xion to even lower values of the “gravity parameter” to render effects of the skin completely negligible in the soft X-ray range. The ionisation effects should nevertheless be very small at high energies where the two codes should yield very simillar results. However the Compton hump extends to much higher energies in pexriv than in xion, and the normalisation is different by a factor of 1.3 at 20 keV.

The difference in normalization is due mainly to the difference in illumination law used. Magdziarz & Zdziarski (1995) have specific intensity $`1/\mathrm{cos}\theta `$, while xion assumes a single ray, incident onto the surface at $`45^{}`$. The first form of the illumination law is appropriate for an optically thin full corona (with a covering fraction equal to unity) which is now known not to work for Cyg X-1 (Gierlinski et al. 1997); the use of a fixed value of $`\theta =45^{}`$ in xion, on the other hand, stems from the necessity to keep the xion runtime to a manageable minimum. Thus both illumination laws are not expected to be strictly correct in reality and are approximations good to some 10-20 %. Also, xion is a table model so linear interpolation from the nearest tabulated parameter values can introduce some inaccuracy.

The deficit of photons at $`200`$ keV is a more serious problem. It is mainly due to the energy grid of xion only extending to 200 keV, again from the need to keep the xion runtime to a managable level. However, it is these high energy continuum photons (which are included in pexriv) which are downscattered to form the reflected continuum above 30 keV. Thus xion as tabulated here gives too few reflected photons at high energies compared with pexriv.

An ideal model would accurately calculate both the ionization structure and compton scattering for a broad range of parameters and do it fast. Since such model is not feasible due to computer limitations, we use the ionization structure of xion at low energies ($`<20`$ keV) and the compton downscattering calculations of pexriv at high energies ($`>30`$ keV), having matched the normalisation of the pexriv reflection to that of xion at 20 keV. The gap in good data between 20–30 keV means that the two different models can be used for the PCA and HEXTE data, respectively, rather than having to interpolate between them.

For a single temperature compton continuum fit to the data above 10 keV in this approach with the xion-disc model gives $`\chi _\nu ^2=62/59`$. By comparison, the reflection model built into the compps code (which is based on pexriv) has $`\chi _\nu ^2=65/59`$ showing that the matching condition is adequate.


format: text

text:
# WeBo 1: A Young Barium Star Surrounded by a Ring-Like Planetary Nebula

## 1 Introduction

In this paper we report the discovery of a new planetary nebula (PN) in Cassiopeia. The PN appears nearly perfectly elliptical, and its central star is a late-type star with enhanced abundances of carbon and s-process elements. The system thus appears to represent the immediate aftermath of the formation of a barium star. In the following sections we present the serendipitous discovery of the nebula, the classification of the nucleus as a Ba II star, a study of the star’s photometric properties and variability, estimates of the distance to the system, and an evolutionary scenario for the origin of this remarkable object. We conclude with suggestions for follow-up studies.

## 2 The Planetary Nebula WeBo 1

### 2.1 Discovery

In 1995, one of us (R.F.W.) discovered an unusual nebulous object while verifying the coordinates of LS I +61303, a well-known stellar X-ray source. The Digitized Sky Survey (DSS)111The DSS was produced at the Space Telescope Science Institute under NASA grant NAG W-2166. The images of these surveys are based on photographic data obtained using the Oschin Schmidt Telescope on Palomar Mountain. was used to generate an image of LS I +61303. During examination of the source and its surrounding field, R.F.W. noticed a faint, elliptical nebula surrounding a 14th-mag star, lying $`5^{}`$ to the southwest of (and unrelated to) LS I +61303.

H.E.B. subsequently obtained narrow-band CCD images of the object, with the 0.9-m telescope at Kitt Peak National Observatory (KPNO) in autumn 1996, and later with the KPNO Mayall 4-m telescope as described below. Filters isolating \[O III\] 5007 Å and H$`\alpha `$+\[N II\] 6584 Å showed the object to be a previously unrecognized222The SIMBAD database reveals, however, that two other groups have independently noticed WeBo 1, but did not recognize it as a PN. Hau et al. (1995) included it in a list of candidate galaxies lying at low galactic latitude. Martí et al. (1998) detected WeBo 1 with the Very Large Array at 6 cm during their search for extended radio emission associated with LS I +61303, and confirmed it optically using the DSS; they suggested that it is a small H II region. PN.

The name WeBo 1 was proposed by Bond, Ciardullo, & Webbink (1996), and will be used here. In the nomenclature of the Strasbourg-ESO Catalogue of Galactic Planetary Nebulae (Acker et al. 1992), which is based on galactic coordinates, the object’s designation would be PN G135.6+01.0. The central star is listed in the USNO-A2.0 catalog at the J2000 coordinates given below in Table 1. Unusually for a PN nucleus, the star appears to be quite red: the USNO-A2.0 approximate photographic $`B`$ and $`R`$ magnitudes are 16.0 and 14.4, respectively.

### 2.2 Nebular Morphology

Our best narrow-band images of WeBo 1 were obtained by H.E.B. on 1999 January 16 with the KPNO 4-m telescope and its Mosaic CCD camera. The \[O III\] and H$`\alpha `$+\[N II\] images are shown in Figs. 1a and 1b, respectively. The frames were taken through cirrus clouds, but under good seeing conditions: the FWHM of stellar images in \[O III\] and in H$`\alpha `$+\[N II\] was $`1\stackrel{}{\mathrm{.}}0`$ and $`0\stackrel{}{\mathrm{.}}8`$, respectively.

In H$`\alpha `$+\[N II\], WeBo 1 has a striking morphology, appearing as a nearly perfect ellipse with major and minor axes of about $`64^{\prime \prime }\times 22^{\prime \prime }`$. The shape strongly suggests that the PN is a thin circular ring with a very low ratio of height to radius, viewed at an inclination of $``$$`69^{}`$. Such a morphology is nearly unique among PNe, matched only by the southern-hemisphere PN SuWt 2, which has a nearly identical appearance (Schuster & West 1976; Bond, Exter, & Pollacco 2001). The H$`\alpha `$+\[N II\] ring has a generally clumpy appearance, a bright, sharp inner rim, and is brightest at the two ends of its major axis (perhaps simply a path-length effect). The interior of the ring is almost hollow, although the brightness level is brighter than the surroundings (which are overlain with diffuse H$`\alpha `$ emission at this low galactic latitude).

In \[O III\] the PN appears more diffuse, but the ring is still apparent. Unlike the H$`\alpha `$+\[N II\] image, the interior of the ring is not hollow. The images suggest a gradient in ionization level, with the interior of the PN filled with high-ionization material radiating in \[O III\], surrounded by cooler material around the periphery of the ring emitting strongly in (presumably) \[N II\]. Images through filters that separate H$`\alpha `$ and \[N II\] would shed further light on the morphology of WeBo 1.

## 3 The Central Star

### 3.1 Classification as a Barium Star

A CCD spectrogram of the central star of WeBo 1 was obtained by D.L.P. at the 2.5-m Isaac Newton Telescope at La Palma, Canary Islands, on 1997 January 12. The Intermediate Dispersion Spectrograph was used at a dispersion of 35.3 Å mm-1. The spectral coverage was 980 Å, at a resolution of 1.6 Å, and the exposure was 1800 s.

Fig. 2 shows the spectrum. To our surprise, the star is cool and chemically peculiar. Molecular bands of C2, CH, and CN are prominent, and the lines of Sr II at 4077 Å and Ba II at 4554 Å are extremely strong. Remarkably, the core of the Ca II H line is filled by very strong emission. (Ca II K is, unfortunately, just outside the spectral range that we covered.)

The star thus has all of the properties of a classical barium star (or “Ba II star”). Ba II stars were first identified as a spectroscopic class by Bidelman & Keenan (1951). They show enhanced abundances of carbon and of heavy elements produced by s-process neutron-capture reactions. The modern view of Ba II stars (see McClure 1984; Jorissen et al. 1998; Bond & Sion 2001) is that they are the binary companions of more massive stars that became asymptotic-giant-branch (AGB) stars, dredged up C and s-process elements from their interiors, and transferred them to the companions. The AGB stars have now become optically invisible white dwarfs, leaving the contaminated companions as the visible Ba II stars.

We classified the spectrum of WeBo 1 by displaying it as a “photographic” spectrogram, and comparing it visually with the atlas of Keenan & McNeil (1976). To estimate the Ba II strength, which classically (Warner 1965) is on a scale from 1 (slightly enhanced over normal) to 5 (extremely strong), we compared our spectrum with the sequence illustrated by Lü et al. (1983, their Fig. 2). The resulting classification is K0 III:p Ba5. The luminosity class is essentially indeterminate from our material (due to the abnormal strength of all of the luminosity indicators at this resolution), and the “p” refers to the presence of C2 bands, which are generally seen only in a subset of the Ba II stars, as well as to the Ca II emission.

### 3.2 BVI Photometry

CCD BVI photometry of the nucleus was obtained by H.E.B. on six occasions on four different photometric nights in 1996 and 1997 September. The KPNO 0.9-m reflector was used with a Tektronix CCD, and the photometry was reduced to the Johnson-Kron-Cousins BVI system through observations of Landolt (1992) standard stars. Two nearby stars were chosen as comparison stars for a variability study (see below) and were also calibrated to the BVI system. Results are given in Table 1. The errors, determined from the internal scatter among the six different measurements, are $`\pm `$0.01 mag in $`V`$ and $`BV`$, and $`\pm `$0.02 mag in $`VI`$.

The central star of WeBo 1, at $`BV=1.72`$ and $`VI=1.77`$, is very red for its spectral type; a normal K0 III star has $`BV0.98`$ and $`VI1.00`$ (Bessell 1979). The star thus appears to have an interstellar reddening of $`E(VI)0.77`$, corresponding to $`E(BV)0.57`$ (using the formula of Dean, Warren, & Cousins 1978). The intrinsic color is $`(BV)_01.15`$, suggesting that the star is somewhat redder in $`BV`$ than a normal K0 giant because of its strong carbon bands as well as the influence of the broad Bond-Neff (1969) flux depression around 4000 Å that is seen in barium stars.

Given that the optically visible star in WeBo 1 is very cool, we can speculate that it is a member of a binary system with a much hotter, but optically inconspicuous, companion that is responsible for the ionization of the PN.

### 3.3 Photometric Variability

In order to search for variability, H.E.B. obtained CCD images in BVI on 40 occasions spread over 6 nights each in 1996 and 1997 September, 4 in 1998 March, and 2 in 1999 January. The 0.9-m KPNO telescope was used, as described above. Since several of these nights were non-photometric, we calculated differential magnitudes between WeBo 1 and the sum of the intensities of the two comparison stars.

The central star is definitely variable, but at a low level. We searched the data for periodic variations, using the Lafler-Kinman (1965) algorithm. The smoothest light curves are obtained for a period of 4.686 days. However, the nearby alias periods at 4.347 and 5.074 days also give plausible light curves, so we cannot make a definitive claim that 4.686 days is the correct period. It is even possible to fit the data with a sub-day alias near 0.82 days, but with significantly larger photometric scatter.

In Fig. 3 we plot the differential magnitude and color curves for the adopted ephemeris $`T_{\mathrm{min}}=\mathrm{HJD2450346}.6000+4.686E`$. To estimate the uncertainties, we calculated the standard deviations of the comparison star 1 minus comparison star 2 magnitudes, which are 0.011, 0.007, and 0.008 mag in $`B`$, $`V`$, and $`I`$ respectively. The radii of the plotted points in Fig. 3 were then set to these values (with the errors combined in quadrature for $`BV`$ and $`VI`$). The WeBo 1 nucleus shows a nearly sinusoidal variation in all three filters, with a peak-to-peak amplitude of about 0.03 mag. There are no significant changes in color.

We can identify four mechanisms that could produce this variation: (1) pulsation, (2) heating (reflection) effects in a close binary, (3) ellipsoidal variations (at a period of $`2\times 4.7`$ days), and (4) starspots on a rotating star. Pulsation appears unlikely because of the lack of any color variations. For a planetary-nucleus cooling age comparable to the estimated expansion age of the nebula (see below), we estimate its luminosity should exceed $`10^2L_{\mathrm{}}`$, in which case heating of the facing hemisphere of the Ba star in a 4.7-day binary would produce a reflection effect nearly two orders of magnitude larger than the observed photometric amplitude (unless the binary is seen nearly pole-on), with colors decidedly bluer than observed. The reflection effect should also dominate ellipsoidal variability at all inclinations, unless the barium companion is not the ionization source for the PN. We consider starspots the most likely explanation, in which case 4.7 days is the stellar rotation period. The strong Ca II emission indicates that the star is very chromospherically active, supporting the starspot interpretation and suggesting that the activity is dynamo-driven by the relatively rapid rotation.

## 4 Distance, Stellar Radius, and Rotation

We can constrain the distance to WeBo 1 using three different methods. The first one uses the estimated interstellar reddening of $`E(BV)0.57`$ and the three-dimensional extinction model of Arenou, Grenon, & Gómez (1992). The observed reddening is close to the asymptotic limit at large distances in this region of the sky, but the Arenou et al. model allows us to set a 1$`\sigma `$ lower limit on the distance of 0.75 kpc, a result which is in excellent accord with the extinction maps of Neckel & Klare (1980).

An upper limit to the distance of WeBo 1 follows from equating the 4.7-day photometric period with the rotational breakup period of the barium star. Assuming a mass of $`1.5M_{\mathrm{}}`$ for that star (Jorissen et al. 1998), we find that its radius cannot exceed $`13R_{\mathrm{}}`$. The observed magnitude and reddening then imply that a K0 III star of this radius must be no farther than 2.4 kpc.

A third method uses a statistical distance scale for PNe. There are several of these in the literature, but we have chosen that of Cahn, Kaler, & Stanghellini (1992, hereafter CKS). Ciardullo et al. (1999) have shown the CKS formalism to produce reasonably good agreement with PNe of known distance, albeit with the factor-of-two scatter typical of all statistical PN distance scales. Adopting the 6 cm flux density of 4.4 mJy given by Martí et al. (1998), and an angular radius of $`32^{\prime \prime }`$, we use the CKS method to yield a distance estimate of 2.2 kpc.

We adopt a nominal distance of $`1.6\pm 0.85`$ kpc (the approximate mean and range defined by the lower and upper limits). At this distance, the absolute visual magnitude of the central star is $`M_V=+1.3_{+1.6}^{0.9}`$. Thus it has roughly the luminosity of a normal red giant (in agreement with most of the classical Ba II stars, e.g., Mennessier et al. 1997). At this distance, the major axis of the PN is $`0.48\pm 0.26`$ pc, and the radius of the central star is $`9\pm 5R_{\mathrm{}}`$. On the assumption of a nominal $`20\mathrm{km}\mathrm{s}^1`$ expansion velocity, the age of the PN is 12,000 $`\pm `$ 6,000 yr. If 4.7 days is the rotation period of the Ba II star, then with a radius of $`9R_{\mathrm{}}`$ and an inclination $`i=69^{}`$ (the apparent inclination of the PN), it would have a projected rotational velocity $`v\mathrm{sin}i90\mathrm{km}\mathrm{s}^1`$. Since this is less than the $``$$`110\mathrm{km}\mathrm{s}^1`$ resolution of our spectra, we would not expect to see an obvious spectroscopic signature, even though this would be about 2/3 of the rotational breakup velocity.

## 5 Discussion and Future Research

WeBo 1 has several properties in common with the class of “Abell 35”-type planetary nuclei. This class was defined by Bond, Ciardullo, & Meakes (1993), and has been discussed by Bond (1994), Jasniewicz et al. (1996), Jeffries & Stevens (1996), and Gatti et al. (1997). In the three known A 35-type nuclei (A 35, LoTr 1, and LoTr 5), a rapidly rotating late-type giant or subgiant is seen optically, while a hot companion is detected at UV wavelengths. The cool components vary photometrically with periods of a few days, corresponding to their rotation periods. A definitive orbital period has not been found for any of these three objects from radial-velocity studies, suggesting that the orbital periods may be long. This suspicion is confirmed in the case of the field star HD 128220, which lacks a PN but is otherwise similar in all respects to the A 35-type nuclei: its orbital period is 872 days (Howarth & Heber 1990). These systems, then, have almost certainly not undergone common-envelope interactions, which would have decreased their orbital periods by substantial amounts.

As discussed by Jeffries & Stevens (1996), there is a closely related class of wide binaries containing hot white dwarfs and cool, rapidly rotating, magnetically active dwarfs. These authors propose a mechanism in which an AGB star in a wide binary develops a dense stellar wind, part of which is accreted by the companion star. Their calculations suggest that significant spin-up of the companion may occur, along with accretion of chemically enriched material from the AGB star. Although Jeffries & Stevens considered their suggestion somewhat speculative (in the absence of actual 3-D hydrodynamical simulations of the accretion and spin-up), observational support has arisen in the past several years. This includes the finding of mild Ba enhancements in the rapidly rotating dK component of 2RE J0357+283 (Jeffries & Smalley 1996), and in the nuclei of A 35 and LoTr 5 (Thevenin & Jasniewicz 1997). WeBo 1, with its pronounced Ba and C overabundances, now provides further support.

A scenario that emerges to explain the properties of WeBo 1 is thus as follows. The progenitor system was a fairly wide binary whose components had nearly equal masses (initial mass ratio $``$0.98). The more massive star evolved to the AGB stage, at which point the less massive component had also begun to ascend the red-giant branch. The AGB star was constrained to rotate with the orbital period, so that, as it developed a dense wind, the wind was ejected preferentially in the orbital plane, leading to the ring-like nebular morphology. The wind was enriched in Ba and other s-process elements, and had $`\mathrm{C}/\mathrm{O}>1`$. A portion of the wind was accreted by the companion giant, spinning it up to the observed 4.7-day rotation period, and contaminating its photosphere with Ba, C, and other pollutants. At present, the AGB star has completely shed its envelope, exposing its hot core whose UV radiation ionizes the ejected ring.

Several follow-up studies are clearly warranted. (1) Ultraviolet spectra would confirm the expected presence of a hot companion; a determination of its surface gravity from its Ly$`\alpha `$ profile would provide a mass determination, and thus an estimate of the progenitor’s initial mass and constraints on evolutionary scenarios. (2) Radial velocities from ground-based spectra would allow a search for the orbital period, or set a lower limit if, as we suspect, the orbital period is long. (3) High-dispersion spectra of the barium star should be obtained, to determine its rotational velocity and atmospheric elemental abundances. In particular, it would be of great interest to search for lines of technetium, which should be present if the star is really a very recently created barium star. (4) An abundance analysis of the nebula would also be of interest, since it may be possible to detect lines of heavy s-process elements.

The discovery of WeBo 1 would not have been made without the DSS, and we acknowledge the roles of the late B.M. Lasker and the many other persons involved in creating this valuable resource. KPNO is operated by AURA Inc., under contract with the National Science Foundation. We made use of the SIMBAD database, operated at the CDS, Strasbourg, France. R.F.W. thanks B.D. Fields for useful discussions, and acknowledges support from NSF grants AST 92-18074 and AST 96-18462. H.E.B. thanks R. Gallino and D. D. Clayton for useful suggestions, and W. P. Bidelman for introducing him to Ba II stars more than three decades ago.


text:
# High field magnetic resonant properties of 𝛽'-(𝐸𝑇)₂𝑆𝐹₅𝐶𝐹₂𝑆𝑂₃

## I Introduction

The $`(ET)_2X`$ charge transfer salts represent a large number of materials that possess properties ranging from insulator to superconductor ($`ET`$ is bis(ethylenedithio)-tetrathiofulvalene). The same can be said for the sub-class of compounds $`(ET)_2SF_5RSO_3`$, where $`R=CH_2,CF_2,CHF,CH_2CF_2`$wardandschlueter . Modifications of the anion result in fundamental changes in the structural and physical properties, and investigations on these materials have been pursued to determine the relation of the anion to the material behavior. Some interest is focused on the anion relation to the conduction path formed by $`\pi \pi `$ orbital overlap of sulfur in the ET donor moleculesmori -mold . The role of the anion in electrical properties is still unclear. However, it is certain that the anion plays a role in partially determining the structure of the compound and accepting an electron from the ET molecule, forming a potential conduction hole. The hole then promotes the spin character of the system, providing magnetic interaction through band magnetism. It is interesting to investigate such systems that show a transition of a low dimensional magnetism into a three dimensional (3D) AFM state, as well as those which exhibit a transition into a spin-Peierls state.

$`\beta ^{}(ET)_2SF_5CF_2SO_3`$ is characterized by pairs of $`(ET)^{+1/2}`$ molecules which contribute a single electron to the $`SF_5CF_2SO_5^1`$ anion acceptor, forming a charge transfer salt, which stacks with face-to-face ET contacts (Fig.1). The ET dimers form linear spin chains along the *b*-axiswardandrutel .

A measurement of the SQUID susceptibility for $`\beta ^{}(ET)_2SF_5CF_2SO_3`$ (Fig.2) shows a broad maximum at T$``$115 K indicating the low dimensional character of the interactions.

The susceptibility was fit using the Bonner-Fisher model for S=$`\frac{1}{2}`$ Heisenberg antiferromagnet (HAF)bonner . The best fit for $`T>40`$ K was obtained using a coupling constant $`J/k_B`$ $`=`$ 352 K. Thus, the Bonner-Fisher fit together with a crystal structure analysis suggests a $`S=\frac{1}{2}`$, linear magnetic chainrutel .

A kink in the susceptibility is observed at T$``$45 K (which was suggested to be a spin-Peierls transition temperature wardandschlueter ), followed by a pronounced drop at T$``$20 K. As temperature drops below 20 K there is an intrinsic feature which makes it difficult to reveal the nature of the ground state in the low temperature regime. Previous V-band ESR experiments revealed a collapse of the gapless g$``$2 mode below 33 K. The change in the excitation spectrum was interpreted to as possible evidence of the spin-Peierls transitionwardandrutel . A change in the FIR properties was observed below T$``$45 K, which indicate weak lattice distortions in $`\beta ^{}(ET)_2SF_5CF_2SO_3`$ pigos . In order to study the nature of the low temperature transition in $`\beta ^{}(ET)_2SF_5CF_2SO_3`$ , we decided to perform ESR experiments in an extended frequency-field range, using high-field high-resolution submillimeter/millimeter wave ESR spectrometer recently developed at the National High Magnetic Field Laboratory (NHMFL), in Tallahassee, FLzvaygin . In addition, low-temperature properties of $`\beta ^{}(ET)_2SF_5CF_2SO_3`$ are probed using a pulsed NMR spectroscopy technique.

## II Experimental Techniques

The ESR method is a powerful tool in studying the nature of the ground states in solids (see, for instance, refs. dumm1 , dumm2 , and sergei2 ). We investigate the $`\beta ^{}(ET)_2SF_5CF_2SO_3`$ system employing ESR and NMR techniques at high frequencies and magnetic fields. High resolution ESR was carried out in the 25 T W.M. Keck resistive magnet (with magnetic field homogeneity, 12 ppm/cm DSV) at the NHMFL, using backward wave oscillator (BWO) sources. Investigations were performed by holding temperature and frequency constant and sweeping the magnetic field. Details of the high-field millimeter and submillimeter wave facility can be found in Zvyagin, *et al.*zvaygin .

In our experiments we used single-crystalline samples with a typical size of 0.5x0.1x0.1 $`mm^3`$wardandschlueter . To measure such small samples a special sample holder was built. It consists of two metallic horns, conically reducing the microwave radiation down to the sample size (Fig.3). A metallic foil was placed around the sample to reduce noise caused by background radiation. The sample was placed in the Faraday configuration, i.e. the propagation vector of the radiation $`\stackrel{}{k_\omega }`$ is parallel to the external magnetic field H. An optical chopper was used to modulate the microwave power. H was oriented perpendicular to the *c*-axis of the sample. The probe was coupled to the source via circular waveguide sections and directed into the sample space, where the microwave excitation passed through the sample, a polymer platform, and then a marker compound of DPPH. The radiation was then directed out of the probe and into a helium cooled InSb detector.

NMR investigations were performed in a 15 T superconducting magnet. A home-built MagRes2000 spectrometerreyes  and a capacitively tuned and matched probe were employed. The investigation focused on two nuclei, the $`{}_{}{}^{19}F`$ nucleus on the anion, and $`{}_{}{}^{1}H`$ in the ethylene groups on the donor molecule(ET). Remarkably, the NMR results revealed similar behavior for these two sites which are also independent of field (or frequency), suggesting a common mechanism for the spin dynamics at the anion layer and at the ET molecule. For brevity, we present a representative graph showing the results of the 340 MHz (8T) investigation of the 1H nuclei.

## III Results and Discussion

We first present the NMR data performed in the temperature range $``$3-300 K and show evidence indicating that the transition seen at T$``$20 K is not spin-Peierls in nature. The characteristic behavior of the nuclear spin-lattice relaxation rate ($`T_1^1`$) in a prototypical spin-Peierls material $`CuGeO_3`$ is described in Ref.fagot . Above the transition temperature $`T_{sp}`$, where the system would be in a uniform phase, the fluctuations from weakly interacting spins contributes to a large but weakly temperature-dependent relaxation rate. In $`CuGeO_3`$fagot , this behavior persists until the temperature drops below $`T_{sp}`$, where $`T_1^1`$ decreases sharply (by several orders of magnitude). There is also a change in the magnetic hyperfine shift which follows the $`T_1^1`$ curve, decreasing below the transition temperature.

A typical NMR 1H spectrum in $`\beta ^{}(ET)_2SF_5CF_2SO_3`$ is shown in Fig.4a, taken at 8.1 T. The data for $`{}_{}{}^{19}F`$ NMR is not shown. The $`{}_{}{}^{1}H`$ data shows a narrow ET signal and a spurious signal at a lower frequency. This second signal originates from the probe which was difficult to eliminate completely. This signal has a much longer $`T_1`$ than the ET signal and remains well resolved and separated throughout the investigation.

In Fig.5 we show the temperature dependence of $`{}_{}{}^{1}H`$ $`T_1^1`$. At high temperatures, the relaxation rate is slightly temperature-dependent. As the temperature is lowered, instead of the expected drop at $`T_{sp}`$, we observed a sharp peak reaching a maximum at around 20 K followed by a sudden drop at lower temperatures. Concomitantly, the FWHM, shown in Fig.4b, also increases around 20 K, while the resonance frequency (Fig.4b) drops. These data is suggestive of some form of magnetic order. Indeed, this sample shows remarkably similar behavior to the $`LiMn_2O_4`$ system, a material with an AFM ground state.

The data for the 19F NMR (not shown) also reveals the same peak in $`T_1^1`$ at the same temperature. The similarity in these behavior at two different sites, both the $`{}_{}{}^{1}H`$ and $`{}_{}{}^{19}F`$ (on the ET molecule and anion respectively) suggests that the relaxation at both sites is dominated by a single spin dynamicsrutel .

A representative transition to an AFM phase is given by Sugiyama *et al.* in the $`LiMn_2O_4`$ systemsugi . In the case of AFM ordering a peak is found in the $`T_1^1`$ scattering with the Néel temperature ($`T_N`$) indicated by the maximum of the peak. Both the full-width at half maximum (FWHM) and the magnetic hyperfine shift show increases at $`T_N`$.

Generally, a 3D antiferromagnetic ordering produces large hyperfine fields at the nuclear site, thus shifting the resonance and rendering the NMR line unobservable in the detected frequency range. In the present case, the signal remains tractable down to the lowest temperature measured. This may occur if the ordering is not long-ranged in nature or if the sample is not uniform, which is unlikely. It is also important to note that the bulk susceptibility does not show any direct signature of AFM ordering. In addition, although the FWHM follows the correct trend, it is two orders of magnitude smaller than the $`LiMn_2O_4`$ data. Nevertheless and whatever the case maybe, it is clear from the NMR data the absence of any transition at $`T_{sp}`$ ($``$40 K) and that the anomaly at $``$20 K is definitely something other than a spin-Peierls transition.

We now present a discussion of the ESR investigation. In the case of the low-frequency ESR (V-band and W-band) we have a single resonance line which remains relatively constant in width, amplitude and frequency above 20 K and shows a diminishing amplitude, slight broadening, and small shift to lower field below 20 K. The resonance field data, H$``$*c*-axis also show a relatively constant value with a slight change to lower field below 20 Kwardandrutel . The experiment was repeated using the sample investigated in both the NMR and high-frequency ESR (see below). Figure 6 reports data taken at 98 GHz, where similar shifts, broadening, and diminishing amplitude are evidence for the reproducibility of the previous investigationwardandrutel .

This behavior is quite different from that observed at the higher frequency and field. We first discuss the ESR results seen as the raw data in Fig.7, which shows the sample response with changing temperature. At higher temperatures we observe one signal, corresponding to the sample absorption. As the temperature decreases we find the sample absorption line broadens and eventually splits, where the onset of the broadening occurs around 20 K. Following the line to lower temperatures, the broadened line begins to resolve into two, and then three distinct lines (denoted by arrows 1, 2, and 3 in Fig.7).

The high-frequency case (broadening followed by splitting) can be explained by the appearance of three AFM modes at some fixed frequency in the fields: below $`H_{spinflop}`$, above $`H_{spinflop}`$ and in an intermediate critical regime at $`HH_{spinflop}`$turov , which were observed experimentally for instance in $`MnF_2`$hagiwara . The high-frequency-field data indicate the appearance of new excitation modes, that may be understood in terms of the model developed for spin excitations in a uniaxial AFM with H$``$*easy*-axisturov . Of course, tilting of the principal axis (even very small) and a low symmetry (P1) of $`\beta ^{}(ET)_2SF_5CF_2SO_3`$ can cause some deviations from the modelturov .

Both sets of data show the onset of a transition at T$``$20K, where the transition temperature is also verified in the high-frequency data by the onset of splitting and remaining area under the curve seen in Fig.8 and Fig.8(inset). However, several pronounced differences appear between the two sets of data. The first difference is the persistence of excitations below about 10 K in the high-frequency data seen in Fig.8 as opposed to that reported in Ref.wardandrutel . Second, the high-frequency data show the splitting of the line at low temperatures (Fig.8), which may indicate the onset of 3D AFM ordering.

Remarkably, in both cases the behavior is qualitatively similar to that for AFM ordering. The low-frequency case (broadening and shifting) corresponds to similar features reported by Mitsudo *et al.* in $`LaMnO_3`$mitsudo . The broadening is due to 1) the enhancement of the AFM correlations and 2) the change of the frequency-field slope. This leaves us with the conclusion that both the low-frequency and the high-frequency results may suggest an AFM type ordering in $`\beta ^{}(ET)_2SF_5CF_2SO_3`$ at low temperature.

We should also make mention of the intrinsic nature of the anomaly in the low-T portion of the susceptibility. Assume that the low-T part of the susceptibility is caused by the paramagnetic impurities (chain ends, defects, other impurities). This impurity signal (which appears to be non-negligible in the susceptibility) should manifest itself in the ESR as excitations with g$``$2. We did see the signal from such impurities for T$`<`$ 10 K, where the 98 GHz investigation and Ward *et al.* did not. This suggests that the low-T portion is intrinsic to the system, and may have something to do with the AFM ordering suggested by the ESR.

Let us also discuss the possibility of an incommensurate phase. Is is known that if we apply a high field we can enter into an incommensurate phase, where the magnetic structure is incommensurate with the crystal structure. ESR was studied by W. Palme *et al.* in the incommensurate phase of the spin-Peierls material $`CuGeO_3`$, where hysteresis appears to be a basic property of the systempalme . However, in our investigations we observed no hysteresis. This would again support the non-spin-Peierls nature of the $``$20 K transition.

## IV Conclusion

We investigate $`\beta ^{}(ET)_2SF_5CF_2SO_3`$ , employing both pulsed NMR and high frequency ESR methods. Both the NMR and ESR data provide evidence that the title compound does not exhibit behavior characteristic for a conventional spin-Peierls transition. We submit that both the NMR and the high-frequency ESR investigations provide inconclusive results as to the nature of the ground state in the low-T regime, although there is some evidence of a 3D AFM ordering. It is clear that there is a disparity between the low and high frequency/field dependence, but at present no theory appears to describe the behavior fully. Further investigations using x-ray diffraction and neutron scattering could provide insight into a more exact description of the crystallographic and magnetic structure through the $`T<20`$ K phase transition.

###### Acknowledgements.

The FSU acknowledges support from NSF-DMR 99-71474 and IHRP 500/5031. A portion of this work was performed at the National High Magnetic Field Laboratory, which is supported by NSF Cooperative Agreement No.DMR-0084173 and by the State of Florida. Work at Argonne National Laboratory was supported by the Office of Basic Energy Science, Division of Materials Sciences, US Department of Energy, under Grant W-31-109-ENG-38. Research as Portland State University was supported by NSF Grant CHE-09904316.


text:
# 1 Entropy as a function of 𝑞² and 𝑏.

Recent observations of spectral lines of distant quasars have suggested that the fine structure constant $`\alpha =e^2/\mathrm{}c`$ may have been slightly smaller in the very early universe . Although these claims are still tentative and rather controversial, they have helped rekindle interest in Dirac’s old idea  that the fundamental “constants” of physics may vary in time. In a recent Brief Communication to Nature, Davies, Davis, and Lineweaver have argued that black hole thermodynamics favors theories in which the speed of light $`c`$ decreases in time, and disfavors those in which the fundamental electric charge $`e`$ increases . We show here that when one considers the full thermal environment of a black hole, no such conclusion can be drawn.

The fine structure constant depends on the fundamental charge $`e`$, Planck’s constant $`\mathrm{}`$, and the speed of light $`c`$, and it is natural to ask which of these varies.***Variation of dimensionful constants is inherently ambiguous ; here, we take “varying $`e`$” to mean “suitable variation of all dimensionless quantities that depend on $`e`$.” Davies et al. offer an ingenious argument. A black hole with mass $`M`$ and charge $`Q=ne`$ has an entropy

$$S/k=\frac{\pi G}{\mathrm{}c}\left[M+\sqrt{M^2n^2e^2/G}\right]^2$$
(1)
Evidently an increase in $`e`$ will decrease this entropy, apparently violating the generalized second law of thermodynamics, while a decrease in $`\mathrm{}`$ or $`c`$ will increase the entropy.

As Davies et al. point out, though, such entropic considerations should take into account not just the black hole, but its surroundings as well. An isolated black hole is never in thermal equilibrium: it will always decay by Hawking radiation and, if it is charged, by Schwinger pair production  as well. These processes decrease $`S`$, but they do not violate the second law, since the decrease is compensated by an increase in the entropy of the environment.

To study the full thermodynamics of “varying constants,” one may try to examine the detailed dynamics of heat flow and entropy between a black hole and its environment. A simpler alternative is to consider a proper thermodynamic ensemble, that is, a black hole in a heat bath. To obtain an equilibrium, one must consider a black hole in a “box” of radius $`r_B`$, with fixed temperature $`T`$ and charge $`Q`$ (canonical ensemble) or fixed temperature $`T`$ and electrostatic potential $`\varphi `$ (grand canonical ensemble) at the boundary. For charged black holes, such ensembles were first studied by Braden et al. . In the canonical ensemble, the total entropy takes the form $`S=\pi r_B{}_{}{}^{2}x_{}^{2}`$, where $`x`$ is determined by the seventh-order equation

$$x^5(xq^2)(x1)+b^2(x^2q^2)^2=0$$
(2)
with $`q=\sqrt{G}Q/r_Bc^2`$ and $`b=\mathrm{}c/4\pi r_BkT`$.

While no analytical solution of (2) is possible for generic charge and temperature, a numerical solution is straightforward. Figure 1 shows a graph of $`S_BS/r_B^2`$ against $`q^2`$ and $`b`$, restricted to the highest-entropy configuration. It is apparent—and may be confirmed with more detailed numerical study—that the total entropy increases with increasing $`q^2`$, and thus with increasing $`\alpha `$. For the grand canonical ensemble, exact analytic results can be found, and the conclusion is the same. Black hole thermodynamics thus militates against models in which the fundamental charge $`e`$ decreases, but places no restrictions on models in which $`e`$ increases.

How can one reconcile this result with the argument of Davies et al.? Note first that the Hawking temperature,

$$kT_H=\frac{\mathrm{}c^3}{8\pi G}\frac{\sqrt{M^2n^2e^2/G}}{2M^2n^2e^2/G+2M\sqrt{M^2n^2e^2/G}}$$
(3)
decreases with increasing $`e`$. As $`e`$ increases, a black hole will thus cool below the ambient temperature of the heat bath, and will absorb heat, increasing its mass. By the first law of thermodynamics, the net change in entropy is

$$dS=\frac{1}{T}\left(dE\varphi dQ\right)$$
(4)
and it may be checked explicitly that the increase in quasilocal energy $`E`$ more than compensates for the direct effect of changing $`Q`$.

Of course, thermodynamic arguments of this sort only describe relationships among equilibria, and not the details of the transitions between equilibria. In some ways, this is an advantage: our results are rather insensitive to the details of a theory describing “varying constants,” requiring only that conventional black hole thermodynamics hold at equilibria at which $`\alpha `$ is constant. On the other hand, this generality allows some loopholes: one could, for instance, imagine a scenario in which an abrupt change in $`\alpha `$ led to an initial decrease in entropy, which then grew only slowly as heat was redistributed.

To analyze such possibilities, though, there is no substitute for a detailed dynamical model. In particular, any theory with a variable fine structure constant necessarily contains a new scalar field, $`\alpha `$ itself. The entropy of that field might be negligible at equilibrium, but it surely cannot be ignored during dynamical processes in which $`\alpha `$ is changing. Black hole mass quantization may constrain models with “varying constants” , but in the absence of a detailed dynamical description, it seems that black hole thermodynamics cannot.

Acknowledgments

This work was supported in part by U.S. Department of Energy grant DE-FG03-91ER40674.


text:
# Parallization of Stellar Atmosphere Codes

## 1. Introduction

The last few years saw a shift in high performance computing from the paradigm of vector computing to massively parallel systems which allow to solve increasingly complex, physical problems. We want to concentrate on general concepts to write parallel code which may be helpful to get started and to estimate/choose the right(!) number of CPUs for a given numerical problem. Due to the rapid advances in the hardware, any number for the performance or the hardware mentioned (e.g. Myrinet, Gbit-Ethernet) will be outdated long before this contribution goes into print. To judge the performance of parallel codes, the relevant quantity is the ratio between communication and CPU-speed and, thus, the numbers may still be useful.

Ideally, code for massively parallel computers should be designed and written from scratch, and it should not be based on any legacy code. However, this approach is both costly and very time consuming. Development cycles may be longer than changes in the computational landscape, and the flexibility may be reduced to answer problems in astronomy. Sometimes, it may be beneficial to adopt a current code base because of faster development cycles, higher flexibility and, more important, to use well tested codes. Here, we want to report our attempt of the latter approach to demonstrate problems and possible solutions. Most of our code is written in FORTRAN with some C routines. Our HYDdrodynamical RAdiation transport code is based on a number of individual programs which have been used to carry out many of the analyzes of SNIa and Core Collapse Supernovae (\[Höflich 1988\], …, \[Höflich, 1995\], \[Howell et al. 2002\], …). All components have been written or adopted to a modular form with well defined interfaces which allows an easy coupling (see Fig. 1) and code verification by exchanging modules. The modules consist of physical units to provide a solution for e.g. the nuclear network, the statistical equations to determine the atomic level population, equation of states, the opacities, the hydro or the radiation transport problem. The individual modules are coupled explicitly. Consistency between the solutions is achieved iteratively by perturbation methods (see Höflich 2002, this volume).

The goal of this paper is to provide some help in the transition from existing scalar or vector to parallel codes. After some basics and the general concept, we will discuss various approaches to create parallel code, available tools, and limitations of our approach.

The terms scalability and efficiency characterize how well the execution time decrease with the number of processors N(CPU) and the effective speed per processor compared to N(CPU) one-processor systems. The actual numbers have been obtained using a 20 node PC-cluster with dual-Pentium IIs with 400Mhz, 512MB per node and fast Ethernet interconnections. In addition, we used a 16 node IBM-SP3 system with 4 CPUs and 1GB per node, interconnected by a cross-bar.

## 2. Basics

Three basic approaches to computational resources can be distinguished namely shared memory systems (e.g. Cray T90, NEC), distributed memory systems (Beowulf-clusters) and distributed systems (Grid-computing over the Internet) which represent a sequence from configurations with highly to very loosely coupled CPUs corresponding the trend from expensive to inexpensive computers. Main stream computing uses either shared or distributed memory systems, and we will concentrate on those. Recently, hybrid systems have been introduced which use multiple CPU-nodes to combine the advantages of both.

The basic limitations for multiple-CPUs is the communication namely its latency and bandwidth. The former measures the time it takes to establish a communication link between CPUs, the latter provides a measure for the amount of data per time which can be transfered.

In the following, we want to discuss the concepts and tools. The strategy for parallization depends on the kind of numerical problem. We use either an algorithm based approach, subdivide the computational domain, or use a module based approach.

Before discussing details we want to refer to some WEB-sides which we found useful to get started:

`www.links2go.com/topic/parallel_computing` provides links to a wide variety of resources.

wotug.ukc.ac.uk/parallel and www.cs.rit.edu/ ncs/parallel.html are archives for information on parallel computing.

www.openmp.org and www.crpc.rihce.edu/HPFF/home.html provide a link to communications tools for shared memory systems, namely for OpenMP and HighPerformanceFortran, respectively.

www.mcs.anl.gov and `www.epm.ornl.gov/pvm_home.html` provide a link to communication tools for distributed memory systems, namely MPI and PVM.

www.cs.utk.edu/ browne/perftools-review provides a good overview about tools to analyze MPI and PVM codes.

www.netlib.org provides an excellent resource to math-libraries optimized for parallel machines including ScaLAPACK, BLACS, Linpack, clapack, slatec, sparce, etc.

### 2.1. Shared and distributed memory

The two basic paradigms/concepts of parallel computing are the shared  and distributed  memory model. In shared memory systems, all processors have direct access to a unified memory using high speed buses. This approach avoids the need for a software based data exchange. On the downside, though, the number of processors for a computer system using a shared memory is rather small, and the stringent requirements on the bus-system result in very expensive computers. The distributed memory system allows to connect an unlimited number of nodes to gain large computing power at a low cost using standard network technology or relatively inexpensive and, comparably, slow options such as Ethernet, cross-bars or Myrinet. Each processor has direct access only to the memory on the local node. Data communication has to be handled explicitly by the user. This puts the burden on the user to manage rather sparse resources but gives more control to the user.

### 2.2. Communication tools

In shared memory systems, the communication channels are sufficiently fast to allow parallization for ’moderate’ number of CPUs on a loop-level, very similar to ’classical’ vector computers. Parallization can be achieved by compiler options (auto-tasking, or HPF) or, more efficiently, by directives inserted into the code by the user, e.g. by using OpenMP (OpenMP Review Board, 2000). In distributed memory systems, data communication is achieved by explicit message passing using specialized communication libraries, namely PVM (Parallel Virtual Machine) (Geist et al. 1994), the Message Passing Interface (MPI) (e.g. MPI, 1994, \[Snir et al. 1995\]), or the new standard MPI2 (\[MPI2\]).

PVM has been designed at the Oak Ridge National Laboratory to connect inhomogeneous clusters of workstations. PVM supports fault tolerance and includes a notification scheme, dynamic hosts and it is very portable due to its PD implementation. Dynamical host allocation means that the number of tasks can be changed during a run. This is particular useful in complex calculations to adjust the number of parallel tasks to specific problems, and it may be useful on multi-processor systems used by a large number of users because it allows to free resources 111Of course, the maximum number of tasks per user must be limited. A further advantage of dynamical host allocation is its ”fail-save” feature. If a task terminates due to a failing hardware or a numerical problem, a new task can be started on another CPU or another algorithm to continue the calculation. MPI has been developed to provide fast communication on homogeneous clusters, and it has been adopted as the standard for parallel computing (e.g. Snir et al. 1995). Tasks are static but, as for PVM, the resulting codes are highly portable due to excellent PD implementations (MPICH, 2002; LAM, Byrnes et al. 2002). Static task allocation implies that the number of CPU’s has to be requested at the beginning and cannot be adjusted ’on the run’ and that the entire job terminates if a task fails. MPI is an excellent choice for domain-based parallization and, numerically, very stable algorithm as commonly used in hydro calculations whereas PVM has its virtue if elaborated iteration schemes (e.g. ALI) are used to solve complex systems of non-linear integro-differential equations such as radiation-hydrodynamical problems including statistical equations which, sometimes, fail to converge (see Höflich, this volume) 222Note that dynamical task allocation is not critical if we calculate NLTE-atmospheres for a few snapshots in time (e.g. Hausschildt, this volume).. More precisely, a more aggressive strategy can increase the convergence rates dramatically but, often, on the expense of stability. Even a failor-rate of 1 % of ALI becomes a problem in LC-calculations which require several thousands of time steps (see below). In conclusion, both MPI and PVM are targeted towards distributed memory systems but their design concepts are very different.

To overcome the problems created by two separated communication libraries, the new MPI2 standard has been defined in 1997 by the MPI2-Forum. This new standard combines both the advantages of MPI and PVM, and it includes important enhancements such as single sided communication which greatly decrease the latency for small messages and, thus, boost the overall scalability. Unfortunately, no reference implementation has been released 333The LAM project is being extended to provide a full implementation of MPI2 including dynamical task control which hurts portability of codes. Moreover, most vendor-based implementations are restricted to subset of MPI2 without dynamical task allocation and single sided communication (e.g. IBM).

We use a combination of PVM and OpenMP but hope to migrate to MPI2 in the near future. Using PVM, our code is parallized on the level of large modules, e.g. for the equation of state, statistical equation, or packages of photons for radiation transport reducing the communication overhead. To increase the scalability, OpenMP is employed to distribute the computations on the loop level within time consuming routines and to take advantage of the on-board shared memory of multi-processor nodes.

### 2.3. Latency and Bandwidth

In theory, multi-processing on N processors may reduce the execution time of a code by a factor of N (or more). However, in reality, the typical gain is much lower and, sometimes, multiple CPUs may even cause an increase of the execution time because the speed of non-parallel parts stays the same but communication and administrative overhead increases with the number of processors.

Distributed memory systems have a large latency and low communication bandwidth compared to the processing speed of individual CPUs. It may be useful to ask for the ’break-even-point’ BEP when the communication between two processors requires as much time as the calculations. Some typical examples are given in Table 1. Obviously, parallization is only useful if, per communication, about $`10^3`$ to $`10^5`$ floating point operations are performed and about 10 times more local than global data are accessed. In reality, an increase in speed by parallization can be expected if the workload per task is significantly larger ($`numberofprocessors\times breakevenpoint`$). Parallization is very inefficient on an inner-loop level for all but very large vectors,

### 2.4.  Dynamical Load Balancing

For both the domain and (physical) module based parallization, we distribute the work by dividing the computational load NT-times according to the spatial or frequency coordinates, use groups of photons or individual elements.

In a direct approach, one may distribute packages of size $`N=NT/N(CPU)`$ as following:

do i= 1,N

SEND PACKAGE OF SIZE N TO TASK I

enddo

COLLECT THE RESULTS AND WAIT TILL THE LAST TASK IS FINISHED

This approach works reasonable well for domain-based parallization if the number of cells are about the same in all sub-domains (see section 3.2). However, it becomes very inefficient when a complex system of equations has to be solved. Examples are the solution of implicit equations such as nuclear and atomic networks, or radiation transport by Monte-Carlo type methods for which the actual work depends critical on the optical depth. In our example, all processes have to wait for the slowest task. The direct approach must be modified, and the dynamical range DR must be taken into account. We define the DR as the ratio between maximum to minimum CPU time required for a task. In practice, we use package sizes of $`\stackrel{~}{N}=NT/(N(CPU)DR)`$. DR must be estimated and optimization of the code means improving the estimator for DR 444DR is given as an input parameter but we plan to implement a dynamical estimator.. Small packages increase the overhead in the communication (see table 1). The package sizes and the corresponding work loads are reduced with increasing number of CPUs. Eventually, this limits the scalability of the problem. Often, the performance levels out at 20 to 50 CPUs on our Beowulf cluster. For improvements, we employ auto-tasking or OpenMP for parallel algorithms and secondary masters on multi-processor nodes with share memory (see below).

### 2.5. Concepts for Parallization

Programs may use the Master/Slave approach and/or symmetric multi-processing (Fig. 2). In the master/slave configurations, a master-process is used as hub for the communication between different ’slave’ processes and to balance the load, i.e. to distribute the work. In addition, the master may be used for sequential calculations during its idle time. A centralized approach makes it easy to keep track of the timing and synchronization of tasks, data coherency etc. However, the entire communication load is handled by one master process running on a specific CPU. This produces a bottle-neck which limits the scalability of the code. Direct communication between equal slaves, the symmetric multi-tasking, avoids this problem and allows to take full advantage of the available bandwidth of the system. However, sequential calculations have to be done by each slave, and data coherency may be a problem. In practice, we use a hybrid approach based on a tree scheme (Fig. 2).

## 3. Classes of Computational Problems

The most effective way to parallize code depends largely on the kind of numerical problem. In the following, we want to describe some concepts used in our code HYDRA (Höflich, this volume) and show some results for illustration. Basically we distinguish between three kind of numerical problems: 1) Algorithm, 2) domain and 3) (physical) module based parallization.

### 3.1. Algorithm based parallization

During the last few years, very sophisticated mathematical algorithm have been developed specially designed for parallel systems, in particular for linear algebra problems, Fourier transformations etc. . We use this approach for matrix operations and to solve linear equations. An effective implementation requires a rewrite of the corresponding routines, extensive testing/optimization and, often, is very hardware dependent. Therefore, we use subroutines available on the NET namely ScaLAPACK, a parallel implementation of LAPACK.

As basic building blocks, ScaLAPACK uses BasicLinearAlgebraSubroutines (BLAS). In Fig. 3, the execution speed is given as a function of the rank of a matrix. Graphs are given for a single CPU, 2,4 and 8 processor PC-clusters with fast Ethernet interconnections using TCP-stacks or directly addressing the memory on each node (GAMMA). Obviously, there is no gain in speed for a rank $`1000`$. 2 CPUs are faster than 4 CPUs for $`N<2500`$. Our tests on an IBM-SP3 showed similar results for a matrix of rank 2000. By increasing the number of processors by a factor of 16 ($`464`$ processors) the speedup was very moderate (factor 2.2). Therefore, we do not use algorithm based parallization on PC-clusters or distributed memory systems but restrict this approach to multi-processor, shared memory nodes.

When using the distributed memory model (PVM or MPI) on a dual-PII-node and matrices of rank 100 and 1000, we achieve a speedup of 0.9 and 1.6, respectively. Using OpenMP, the corresponding numbers are 1.3 and 1.5. Therefore, OpenMP is employed for ranks $`1000`$ and, for larger matrices, we use PVM or MPI.

As a general rule for linear algebra problems with of moderate rank ($`5000`$), parallization should be restricted to shared memory systems/nodes. Optimized algorithm and machine specific libraries should be used. For small matrices ($`N1000`$), OpenMP seems to give better performance than explicit message parsing.

### 3.2. Domain-based parallization

Here, the computational domain is subdivided and the workload is distributed over several CPUs. We use this approach for 3-D hydro (\[Khokhlov 1998\]) and for the spherical radiation transport in the comoving frame (Mihalas, Kunacz & Hummer 1975). For the hydro, each processor solves the set of equations in a specific spatial domain. The hydro code solves the compressible reactive flow equations in Eulerian frame using an explicit Piecewise Parabolic Method (PPM) (\[Collela & Woodward 1984\]). For Adaptive Mesh Refinement (AMR), it uses a Fully Threaded Tree (FTT). Communication between neighboring domains is needed to exchange the results at the boundary. The scalability is limited by the latency. For typical problems with an effective resolution of $`1000^3`$ to $`2000^3`$ (see Fig. 4), the code scales well up to 20 and 100 processors with an efficiency of 0.5 to 0.2 on systems with a low latency (e.g. IBM-SP3, table 1). On our Beowulf cluster with fast-Ethernet, the absolute peak performance is already reached on 8 processors with an efficiency of $`0.15`$ due to the large latency of the cluster (see Table 1). Another example is the spherical radiation transport module for expanding atmospheres. The equations are discretized by ND depth and NFR frequency points. Angular resolution is achieved by integrating along ND rays parallel to a main axis with impact parameters according to the radial grid. Typically, we use between 90 and 900 depth points and $`10^4`$ to $`10^5`$ frequency points. Because the velocity field is monotonic, the intensity $`I`$ at the frequency IFR depends on the results at the previous frequency IFR-1 making, on a first glance, the frequency space not suitable for parallization despite the large workload which includes the integration over all rays. In principle, parallization could be done along individual rays. However, the number of operations depends on the number of shells touched by a ray. As a consequence, the workload of individual tasks will vary by a factor of ND and some tasks would have very little work per communication. Both results in a very low efficiency. Instead, we distribute the workload using frequency groups, and employ load balancing as described above. As pay-off, the boundary conditions of each frequency group depends on the neighbors and the solution has to be iterated to obtain the correct intensity $`I`$ at the boundary. In practice, the solution for $`I`$ is pretty well known from the model at a previous time step or the NLTE iteration. Typically the efficiency is 30 to 50 % even on our Beowulf cluster with a fast-Ethernet connection.

### 3.3. Physical Modules

We found it very useful to employ parallization which is based on physical modules, e.g. to solve the statistical equations for all elements at a specific grid point (see Fig. 1). Each physical module is designed to solve a specific problem such as the equation of state, nuclear networks, the radiation transport problem etc.. In figure 1, we indicate the basic variable used to distribute the computational load.

The advantages are obvious. Firstly, each of the modules consists of complex codes (with several thousands of lines) and the solution of the problem requires a fair amount of computations. Secondly, only the results have to be communicated although the amount of local data can be significant keeping small the the communication overhead. E.g. for the nuclear network or the rate equations, a huge amount of atomic and nuclear data are required but only the resulting abundances and level populations have to be communicated. Both these properties are the basic recipe for good scalability and efficiency. On the down-side, this approach requires relatively large memory per node, typically 512MB to 1GB per node.

For the user, this approach is beneficial too. Moderate changes are required in course of the transition from scalar or vector code to massively parallel systems. This is particular true for the master/slave approach. The modular approach with well defined interfaces allows easy code verification (see above). The modules can be based on well tested tested scalar and vector codes, and modules can be parallized and tested one at a time. By using a distributed memory model, the memory space is well separated. Data are only seen globally if they are communicated, all but avoiding memory leaks. Of course, all programmers keep control over all local and global variables and their names. However, by mistake, the same variable name may be used e.g. in the hydro and the nuclear reaction network which, historically, are based on separate codes.

For the data communication between various modules, we use dedicated subroutines. E.g., for the nuclear network, all data relevant are collected by the master and communicated to a dedicated subroutine. Subsequently, this subroutine uses calls to the communication library (PVM or MPI2) to transfer the data to the dedicated routine of the Nuclear Network which makes the data visible to the module e.g. via COMMON blocks.

### 3.4. Some Examples

#### NLTE-LC calculations for Supernovae:

As example, we want to discuss the budget for a spherical NLTE-calculation of SNe Ia light curves. For results, see Höflich (2002, this volume). Not all parts of our code have been parallized and, consequently, scaling is far from ideal. In particular, the 1-D hydro module has not yet been adopted for shared memory computers. However, the solution of NLTE-rate equations and the radiation transport equations are fully parallel. In radiation-hydro problems, the solution of the Boltzmann equation for the radiation transport and rate equations for the level population is more CPU intensive than the hydro by factors of 50 to 100.

The calculations have been done on our 40 processor Beowulf cluster using spherical hydro in Lagrangian coordinates, and detailed equation of state for the atomic level population and expansion opacities. The radiation transport is solved using the moment equations. The Eddington factors are calculated by solving the time-independent radiation transport equation in comoving frame. $`\gamma `$-ray transport is done using our Monte Carlo scheme (Fig. 1). For the atomic level populations, we consider a total of 59 super-levels. Radiation transport has been performed using 2500 frequency groups. The envelope has been discretize by 272 depth points. To calculate the Eddington factors, the hydro-grid has been re-mapped to properly resolve the photosphere by 90 radial points. a) The rate equations and the radiation transport equations are fully parallized. The dynamical range DR has been set to 3 for both. Depending on the phase of the light curve, the efficiency $`E(rate)`$ was between 30 to 50 %. b) As described above, the comoving frame equations are parallized using the frequency grid. Formally, the problem scales almost perfect. However, we need about 2 to 5 iteration for the boundary conditions, reducing the effective efficiency to $`25\%`$. c) The $`\gamma `$-ray transport by MC is used to determine the energy deposition functions. The energy deposition varies very slowly and hardly increases the total CPU time. d) Finally, the 1-D hydro, moment and energy equations are tuned for shared memory systems (i.e. nodes). The CPU-time requirements for a), b) and d) are $``$ 200, 800 and 3.8 CPU-seconds per NLTE-iteration, respectively. At each time step, we need about 5 to 10 NLTE-iteration. 2000 time steps are required for an entire light curve up to day 60. Our LC computations take about 8 to 16 days which is rather favorable to the 110…220 days on a single processor system.

#### 3-D gamma ray transport:

The final example has been chosen as a warning not to extrapolate scalability of a code and to demonstrate the importance of testing and fine tuning. Monte Carlo methods are well known for excellent scalability which, often, is called ’trivial’. In Fig. 5, the increase of performance and efficiency is shows as a function of CPUs. Apparently, scaling is nearly perfect for up to 10 processors in all cases which include photon packages of $`10^5`$ and $`10^6`$ photons, and the use of homogeneous and heterogeneous cluster of workstations. However, with increasing number of CPUs, the efficency drops rapidly for heterogeneous clusters of workstations. The translation of data into processor independent formats for the communication is extremely costly and heterogeneous systems should be avoided. In addition, this example clearly shows the importance to use large work loads and the influence of the ’latency’. Often the latency of a system is more important for the scalability of a problem than bandwidth.

## 4. Final Discussions and Conclusions

Parallelization may allow to solve problems previously not feasible. Stellar atmospheric programs and radiation hydro problems are most suitable because they use very complex physics which, in general allow a (physical) module based approach. This enables to the use of legacy code. Nevertheless, the efforts to change from scalar and vector-computers are non-trivial and time-consuming.

Currently, three architectures are used for parallel computing: 1) shared memory and 2) distributed memory systems, and 3) distributed systems. For most problems, only the first two are a valuable option. Inhomogeneous clusters of workstations (e.g. PCs + Suns) should be avoided. Often, the scalability of a parallel code can be increased significantly by employing both the shared and distributed memory model. Nowadays, systems become rather common which bundle clusters of multi-processor nodes.

Parallization tools on shared memory systems are directive based (e.g. Open MP) or or done by the compiler. Typically, communication within multi-processor nodes is sufficient fast for task distribution on a loop level. These systems are highly preferred for algorithm based problems such a linear equations, matrix inversions etc.

The domain and physical module based approach is very suitable for distributed memory systems, and load balancing can be implemented in a ’straight forward’ way. In particular, parallization on a module basis shows a very good scalability and it is easy to implement. Explicit message passing may be a bonus in large program packages because the memory space is well splitted which greatly reduces the problem of memory leaks between physical modules. Dynamical task allocation as a crucial feature in radiation-hydro problems which include statistical equations because it provides a ’fail-save’ mode (see sect. 3.3). Fortunately, it is realized in both PVM and the new standard MPI2. New projects should use the communication libraries based on the MPI2 standard as soon as full implementations become common place.

For high performance, it is critical to use the ’right’ number of processors. Communication overhead may actually decrease the overall performance significantly. Optimization and dynamical load balancing is critical for good performance.

##### Acknowledgments.

This research is supported in part by NASA Grant LSTA-98-022. The calculations have been performed on the Beowulf cluster of the Department for Astronomy and the HighPerformanceCenter at the University of Texas at Austin.