text:
# Gas Rich Dwarfs from the PSS-II III. H I Profiles and Dynamical Masses

## 1 INTRODUCTION

The study of low surface brightness (LSB) galaxies presents an enormous observational challenge. LSB galaxies often too faint for more than limited optical study by broadband imaging. Spectroscopy and near-IR are practically impossible due to the fact that too few galaxy photons are emitted per square arcsec compared to the brightness of the night sky. Even a nominal number of HST orbits requires rebinning to obtain reliable S/N, which then loses the high resolution advantage of space imaging. However, one intrinsic characteristic has provided an avenue to study this class of galaxies, a majority of LSB galaxies contain copious amounts of neutral hydrogen (Schombert et al. 1992, Impey et al. 1996). This provides both a quick method to confirm the reality of objects, barely visible in the optical, while simultaneously providing their distance and crude kinematic information.

Dwarf galaxies are the most numerous class of galaxies in the Hubble sequence forming a large component of the number density of galaxies. For example, over 80% of the Local Group are classed as dwarfs and their space densities are higher even in cluster environments. They are found in all types of environments, from dense clusters to the field. They span the range of morphological types, from elliptical to highly irregular. And they cover a full range of current star formation rates, from quiescent dE’s to intense BCD’s (Sandage and Binggeli 1984). Although dwarf galaxies are low in metallicity (Skillman and Kennicutt 1993), they are by no means zero-metallicity protogalaxies. And, while dwarf galaxies are interesting for their internal properties and star formation histories, they also serve as probes to the characteristics of dark matter (Ashman 1992) and are tracers of large scale structure (Eder et al. 1989, Zwaan et al. 1995).

Despite numerous studies on cluster and field dwarfs, it is still unclear whether dwarf galaxies are a continuation of the spiral sequence or whether they form a distinct population, possibly fossil remnants from the era of galaxy formation. Information about the gas content of dwarf galaxies is therefore critical to an understanding of their place in the hierarchy of galaxies. Unfortunately, our knowledge of the properties of dwarf irregulars is hampered by small catalogs because of difficulty in detecting them due to their LSB nature and small angular size. In order to improve the sample of dwarf irregular galaxies available for study, we have identified a new sample of 458 candidate dwarfs from inspection of 49 fields of the new Palomar Sky Survey (PSS-II) directed towards the North and South Galactic Caps. Previous papers have presented the catalog and the optical data (Eder et al. 1989, Schombert, Pildis, and Eder 1997, Pildis, Schombert and Eder 1997). This paper presents the single-dish Arecibo 21 cm data for this sample and determines the dynamical mass for dwarfs with optical measurements. The impact of these parameters on the star formation history of LSB dwarfs is explored in a latter paper (Schombert, McGaugh and Eder 2000).

## 2 OBSERVATIONS

In this paper, we report on the 21 cm observations of three samples of dwarf galaxy candidates selected visually from photographic plates of the Second Palomar Sky Survey (PSS-II). All of the samples were chosen so as to optimize their observation with the Arecibo radio telescope, meaning they lie between $`5`$ and $`+35`$ in declination. The first sample (Sample 1) of 178 candidates was obtained from 15 fields located on regions of the fall sky occupied by a nearby void. The plates used were high quality preproduction plates. Of these candidates, 132 were observed at 21 cm. The void is centered at RA = 0h45m, Dec = +20, v = 3500 km sec-1, and is of roughly spherical shape with a diameter of 1500 km sec-1. The second sample (Sample 2) of 278 was identified using methods similar to those for the first sample from 35 random fields between RA’s of 7h to 16h and Dec’s of 10 to 25, irrespective of galaxy density. Of these, all were observed at 21 cm. The third much smaller sample of 27 (Sample 3) consisted of new dwarf candidates identified from a visual reexamination of one of the fields of Sample 1 (Field 409), using the final PSS-II production plates. Sample 3 also included observations of some of the dwarf candidates from Sample 1 which were not originally observed with 21 cm due to a lack of telescope time.

All objects on the plates with angular diameters greater than 20 arcsecs, irregular appearance (Sm, Im, Irr, dI) and low surface brightness were selected. Low mass implies instability, at least with respect to the pattern of star formation. Thus, irregular structure is nominally associated with dwarf galaxies (unless the star formation is completed, as in a dwarf elliptical). Basically, late-type galaxies can be broken down into three morphological types; those with bulges and LSB disks (Sc and Sd class, those with apparent axial symmetry but no bulges (Sm class) and those with no symmetry (Im class). However, morphology is only weakly related to the mass of a galaxy. There are many irregular galaxies that are massive (Hunter and Gallagher 1986) and many spirals which are dwarf-like in luminosity and size (Schombert et al. 1995). In addition, many high luminosity events, such as galaxy-galaxy interactions, are also associated with irregular morphology, although usually of a nature that identifies them as tidal in origin. Irregular morphology, by itself, is not a sufficient criteria for dwarf classification, but is sufficient for the development of a candidate list. The coordinates of each candidate dwarf were taken directly from the plate material using a fine ruler and the position of SAO stars. Accuracy of the coordinates were typically $`\pm `$15 arcsecs based on the centering of CCD imaging, which is sufficient for detection with the Arecibo beam of 3 arcmin diameter. Further details of the selection of the three samples can be found in our previous papers: Eder et al. 1989 (Sample 1) and Schombert, Pildis, and Eder 1997 (Samples 2 and 3).

Follow-up CCD images in V and I for most of those galaxies detected in H I from Sample 2 were obtained on the Hiltner 2.4m telescope located at Michigan-Dartmouth-M.I.T. (MDM) Observatory. We obtained images using either a Thomson 400$`\times `$576 pixel CCD (0.25 arcsec pixel-1) or a Ford-Loral 2048$`\times `$2048 pixel CCD binned 3$`\times `$3 (0.51 arcsec pixel-1), with minimal exposure times of 25 minutes in Johnson $`V`$ and 15 minutes in Johnson $`I`$. Analysis of these images reveal that the candidates are gas-rich dwarf galaxies with total luminosities of $`M_V>18`$ mags, central surface brightnesses between 21 and 23 I mags arcsec-2, exponential scale lengths of typically less than 3 kpc and mean isophotal diameters (of approximately the Holmberg diameter) of 8 kpc (as compared with a mean of 30 kpc for the UGC galaxies, Pildis, Schombert and Eder 1997). A sample of comparable size to ours would be the Virgo Cluster dwarfs (Hoffman et al. 1997) but these dwarfs reside in a cluster environment while ours inhabit a range of environmental densities.

We observed the dwarf candidates for the H I line at 21 cm with the Arecibo 305m telescope during the 1988, 1992 and 1993 observing season. All observations were made with the 21 cm dual-circular feed positioned to provide a maximum gain (8 K Jy-1) at 1400 MHz. The 2048 channel autocorrelator was used and the independent, opposite polarized signals were each divided into two subcorrelators of 512 channels. In order to search a larger velocity space, the secondary local oscillators of each polarization set of subcorrelators were offset on either side of the standard local oscillator frequency of 260 MHz by 8.75 MHz, allowing a total velocity coverage of 8000 km sec-1, a velocity resolution of 8.6 km sec-1, and some overlap at the band edges. The observations were centered on 4000 km sec-1, which avoided detection of the strong Galactic hydrogen signal on the low-velocity end, and extended to 8120 km sec-1. Observations were made in the total power mode with 5 minute ON- source and OFF-source observations. In most cases, only one 5-minute ON-source integration was required for detection. Wherever possible, the zenith angle was kept less than 14 degrees to minimize the degradation of the gain.

Efforts were made to assure that the signal was due to the candidate galaxy. The field of each position was examined on the original plates and also with NED to check for possible sources of confusion. Because many of the detections were very narrow, the possibility existed that they were actually caused by radio frequency interference (RFI). However, RFI confusion can be distinguished by their polarization properties and/or, with 3-level sampling, by the ringing they cause in the band which is eliminated with hanning smoothing. Those narrow detections that were still suspicious were reobserved to confirm their authenticity. Higher resolution observations of several narrow line detections revealed double-horned features.

## 3 ANALYSIS

### 3.1 Galaxy Data

All three samples were analyzed and reduced in the same fashion. Calibration of the spectra for atmospheric correction and noise diode variation with frequency was performed during the phase of translation into export format for compatibility with the Arecibo Observatory ANALYZ-GALPAC analysis system. Baseline fitting, removal, and parameter extraction were all performed using the latter package. The width, referred to hereafter as $`W_{50}`$, is defined as the full width across the profile measured at a level of 50% of the peak flux.

In Table 1 we present the results of the analysis for the 250 galaxies whose digital spectra have been examined. In addition, 143 galaxies which were searched for, but not detected, are also listed in Table 1. The nominal search range was 1000 to 8,000 km sec-1, with occasional extensions to 12,000 km sec-1. The typical RMS values were 3 mJy. All distance related values in this paper use values of $`H_o=75`$ km sec-1 Mpc-1. All the data included in this paper can also be found in ASCII format at zebu.uoregon.edu/$``$js.

Details of the entries in Table 1 are as follows:

* Column (1): Object ID, based on the PSS-II field number and a running number for the dwarfs found on that field.
* Columns (2) and (3): Right ascension and declination in the 1950.0 epoch. These are measured off the PSS-II plates themselves and are accurate to 15 arcsecs.
* Column (4): The heliocentric velocity, $`V_h`$ of the H I line, taken at the midpoint of the profile (50% of peak level).
* Column (5): The velocity corrected for the motion of the Sun relative to the center of mass of the Local Group, according to the precepts of de Vaucouleurs, de Vaucouleurs, and Corwin (1976), $`V_o`$, in km sec-1.
* Column (6): The velocity width, $`W_{50}`$, in km sec-1. No corrections have been made for inclination or turbulent broadening.
* Column (7): The profile shape where S means single-horned and D means double-horned. A semi-colon indicates uncertain classification.
* Column (8): The observed integrated 21-cm H I line flux, $`S=S_V𝑑V`$ in Jy km sec-1.
* Column (9): The H I gas mass of the object, in $`M_{\mathrm{}}`$, calculated using the standard conversion from flux to gas mass.
* Column (10): Observing run number from which the data was taken, where run #1 is Sample 1 observed in 1986, #2 is Sample 2 and #3 is Sample 3, both observed in 1992 and 1993.

### 3.2 Notes on Individual Galaxies

409-4,476-5 – Velocity scale in the spectrum is not correct because the LO synthesizer was set incorrectly. The velocities in the Table are the correct values.

409-7 – 1.5 arcmin from SRGb 062.12 (v=7333 km sec-1)

409-11 – 8.5 arcmin from 409-10

409-16 – 6.1 arcmin from N7773 (v=8486 km sec-1)

409-17 – 3.0 arcmin from KUG 2356+303

409-18 – 1.6 arcmin from U152 (v=4863 km sec-1)

474-9 – 8.1 arcmin from N304 (v=4991 km sec-1

476-2 – also U948

476-4 – also U1073

476-5 – also U1084

495-1 – 5.8 arcmin from U4361 (v=3741 km sec-1)

495-3 – detected as F495-V1 (LSB2)

508-1 – The dwarf HI signal could be embedded in that of NGC 4961, detected here off center, but is buried in flux from primary galaxy. NGC 4961 detected by Shostak. The velocities and widths listed in the Table are measured at 50% of the mean flux

508-2 – near A1656

512-1 – lots of RFI

512-4 – 8.8 arcmin from KUG 1430+264

512-5 – observed in velocity range 0 - 12,000 km sec-1

512-10 – 6.5 arcmin from IC4475

514-2 – observed in velocity range 0 - 12,000 km sec-1

514-3 – paired with 514-4

514-5 – The spectrum is missing from the Figure

514-6 – observed in velocity range 0 - 12,000 km sec-1

516-1 – observed in velocity range 0 - 12,000 km sec-1

538-1 – 4.9 arcmin from U3 (v=7882 km sec-1)

538-2 – 7.8 arcmin from N7817 (v=2309 km sec-1)

538-15 – also U12843

538-16 – also U12846

538-18 – 8.3 arcmin from U12916 (v=6352 km sec-1)

539-2 – also U159

539-5 – also IC1542

539-6 – also U200

540-2 – 7.9 near U425 (v=5841 km sec-1)

540-3 – also U477

540-8 – 1.3 arcmin from F474-10

559-3 – Malin object, 9.3 arcmin from U3920

563-2 – detected as F563-V2 (LSB2)

563-3 – 7.7 arcmin from U4588 (v=4264 km sec-1), detected as F563-V1 (LSB2)

563-4 – dwarf spiral, detected as F563-1 (LSB1)

564-2 – measured parameters are uncertain because of overlapping RFI

564-8 – detected as F564-V3 (LSB2)

564-12 – 7.4 arcmin from U4729 (v=3900 km sec-1)

564-15 – dwarf spiral

565-2 – UGC 5005, detected by Bicay and Giovanelli

565-3 – possibly two merged signals - several nearby galaxies are within 1 arcmin on the CCD field

565-5 – signal is blend of UGC 5086 and NGC 2903 in the beam and sidelobes

565-10 – KARA 68-056, not detected as F565-V4 (LSB2)

568-3 – REIZ 363, 3.4 arcmin from U5629 (v=1236 km sec-1)

568-4 – ZWG 094.005

568-5 – observed in velocity range 0 - 12,000 km sec-1 but had strong RFI (rms=2.3mJy)

570-3 – detected as F570-V1 (LSB2)

570-4 – detected as F570-7 (LSB2), 4.5 arcmin from IC2703

570-6 – poor S/N, companion to D571-5

571-3 – 8.7 arcmin from AGC210538 (v=6196 km sec-1)

571-5 – two galaxy signals appear in the spectrum, due to two galaxies within 20 arcsec on the CCD image. The parameters listed in the table are from the face-on diffuse dwarf galaxy. The other signal is double peaked and therefore emanates from the inclined spiral ($`v`$=5985 km sec-1, W=86 km sec-1, FI = 157.2 mJy km sec-1 )

571-6 – also not detected as F571-10 (LSB2)

572-5 – 8.3 arcmin from U6770

575-1 – I3840, detected as F575-2 (LSB2)

575-2 – UGC 8011

575-3 – I4107

575-4 – pair to 575-3

575-5 – KARA 68.215, in same beam as F575-3 detected in LSB2. The velocities and widths listed in the Table are measured at 50% of the mean flux

576-2 – strong continuum

576-6 – ZWG 101.005

576-7 – measured parameters uncertain because of interference

576-11 – MCG +03-34-034, 2.7 arcmin from U8448 (v=7151 km sec-1), probably U8448 off band

577-5 – dwarf spiral

577-6 – measured parameters uncertain because the galaxy spectrum is near the edge of the band, companion to D648-4

582-1 – probably off-beam detection of F582-2 (LSB2)

582-4 – 9.5 arcmin from F583-2

584-2 – also U10140, 6.1 arcmin from U411A (v=2525 km sec-1)

609-2 – also U385

609-3 – companion to D609-2

609-4 – UGC 472, detected by Bothun et al. 1985

609-5 – also U560

611-5 – detected as F611-1 (LSB2)

611-6 – also U883

611-7 – also U891

611-11 – 6.1 arcmin from F611-V2

611-13 – also U1026

611-14 – also U1056

611-18 – 8.7 arcmin from U1087 (v=4485 km sec-1)

631-3 – questionable signal at v= 168 km sec-1, W = 33, FI = 123 mJy km sec-1

631-4 – questionable signal at v= 352 km sec-1, W = 50, FI = 147 mJy km sec-1 631-7 UGC 4115, detected by Bicay and Giovanelli

631-7 – also U4115

631-8 – dwarf spiral

631-9 – questionable signal at v= 2280 km sec-1, W = 59, FI = 178 mJy km sec-1

631-11 – questionable signal at v= 6130 km sec-1, W = 60, FI = 123 mJy km sec-1

637-1 – very noisy spectrum

637-5 – 7.8 arcmin from IC581

637-7 – 7.5 arcmin from U5342 (v=4560 km sec-1) (pair?)

637-10 – observed in velocity range 0 - 12,000 km sec-1

637-21 – observed in velocity range 0 - 12,000 km sec-1

640-2 – 640-2,3 and 4 are all in the same beam (UGC 6006, UGC 6007). Signal is due to U6006 also detected by Schneider et al.

640-7 – UGC 5948, detected by Schneider et al.

640-8 – near NGC 3489

640-10 – UGC 5944

640-11 – measured parameters uncertain due to overlapping RFI, 5.4 arcmin from N3412 (v=865 km sec-1)

640-14 – near ZWG 066.072

640-15 – NGC 3559, detected by Bicay and Giovanelli

646-4 – Rea 66

646-5 – 646-6 is also in the beam, signal could be from either one

646-7 – UGC 8091, detected by Lewis, Helou, and Salpeter

646-11 – UGC 8061 - detected by Schneider et al.

648-1 – solar interference

648-2 – very bad baseline and noise

648-4 – companion to D577-6

648-6 – solar interference

651-6 – comp to ZWG 076.018

656-1 – U10398 - detected by Schneider et al.

656-2 – U10281 - detected by Schneider et al.

656-4 – observed in velocity range 0 - 12,000 km sec-1

656-5 – 9.2 arcmin from U10218 (v=1080 km sec-1)

685-6 – also U1323

685-18 – also U1515

685-19 – remeasured parameters since Eder et al. , 4.4 arcmin from IC1772

709-5 – KARA 68.060

709-7 – 9.5 arcmin from N3130 (v=8206 km sec-1)

709-9 – 9.2 arcmin from U5304 (v=12308 km sec-1)

709-10 – 5.3 arcmin from U5304 (v=12308 km sec-1)

721-1 – ZWG074.037 also in beam

721-5 – dwarf spiral, detected as F721-V4 (LSB2)

721-7 – 4.1 arcmin from N5463 (v=7235 km sec-1) (pair)

721-8 – companion to D721-9

721-10 – limits of the measured galaxy spectrum are uncertain

721-14 – signal is from U8995 at 140219+090224 detected by Bothun et al. 1985

721-15 – broad wings - parameters uncertain - NGC 5511 interaction

723-3 – signal on edge of spectrum, parameters uncertain

723-4 – 6.3 arcmin from U9356 (v=2225 km sec-1)

723-9 – also U9500

749-3 – also U12416

749-6 – also U12480

749-7 – 5.9 arcmin from N7587 (v=8917 km sec-1)

749-8 – 5.0 arcmin from N7593 (v=4108 km sec-1)

749-9 – 5.8 arcmin from U12497 (v=3761 km sec-1)

749-10 – also U12497

749-12 – also U12553

749-14 – also U12562, 8.2 arcmin from N7641 (v=7872 km sec-1)

749-17 – 8.4 arcmin from U12687 (v=6158 km sec-1)

774-1 – dwarf spiral

774-4 – PK 215+11.1

822-4 – remeasured parameters since Eder et al. , 3.7 arcmin from D822-5

822-6 – also F750-2

822-14 – 8.8 arcmin from N7757 (v=2956 km sec-1) (pair)

### 3.3 Mass Estimates

To determine the dynamical masses from the H I profiles, inclinations and metric sizes are needed. To this end, CCD images were obtained for 116 galaxies at the 2.4m MDM telescope. The isophotometric ellipse fitting was done as described in Pildis, Schombert and Eder 1997. Scale lengths were derived from least-squares fitting of an exponential function to the data in the surface brightness-radius plane.

Due to their LSB nature, the inclination measurements of dwarf galaxies are problematic (Tully et al. 1978, Lo et al. 1993). Since H I rotation velocities are being used, strictly speaking H I diameters should also be used. However, numerous studies have shown that the H I surface densities follow the optical surface densities (see McGaugh 1992). Thus, we adopt the optical parameters for inclination as based on the eccentricity of the ellipse at the 25 mag arcsec-2 isophote. The $`R_{25}`$ isophotal radius, axial ratio and $`I`$ magnitude are listed in Table 2.

The dynamical mass is determined from the virial theorem using the prescription of Staveley-Smith, Davies and Kinman (1992, hereafter SDK). However, to complicate matters, ellipticity studies indicate that dwarf galaxies are primarily triaxal with particular signatures in their structural parameters that confirm their non-oblate nature (Schombert, McGaugh and Eder 2000). If dwarf galaxies are not oblate rotators, then as discussed in SDK, the use of a pure rotation term to calculate the dynamical mass will underestimate this value by as much as 30%. Hence, we adopt SDK’s virial mass estimate, $`M_{dyn}`$, as given by

$$M_{dyn}=2.3\times 10^5(V_{rot}^2+3\sigma ^2)R_{HI}$$

where $`V_{rot}`$ is the rotation velocity, $`\sigma `$ is the isotropic velocity dispersion and $`R_{HI}`$ is the outer radius of the H I emission. The velocity dispersion is assumed to be 10 km sec-1, the standard value for late-type galaxies (Shostak and van der Kruit 1984). High resolution H I maps of gas-rich dwarf galaxies also indicate a typical velocity dispersion of 10 km sec-1 (Skillman et al. 1988, Lo et al. 1993). Note that, for dwarfs without rotation, $`V_{rot}`$ could be replaced with $`\sigma `$ and the dynamical masses is underestimated by a factor of 3, but we have assumed rotation in our discussion.

We calculate $`V_{rot}`$ from the measured H I profile width, correcting for inclination and velocity dispersion such that

$$W_{50}=2(V_{rot}\mathrm{sin}i+1.18\sigma )$$

where the inclination, $`i`$, is determined from the axial ratio of the outer isophotes of the dwarf galaxy. The remaining variable is $`R_{HI}`$, the radius of the edge of the H I emission. We can measure $`R_{25}`$, the radius of the 25 $`B`$ mag arcsec-2 isophote directly from the optical images. H I mapping studies indicate that the H I diameter of dwarf galaxies is much larger than the optical radius. Van Zee, Skillman and Salzer (1998) found that the H I diameter is approximately 4 times the optical radius for a mapping study of 5 dwarfs. However, these dwarfs were selected as having an overabundance of H I. Salpeter and Hoffman (1996) discovered a correlation between optical and H I radius, such that $`R_{HI}=2.3R_{opt}^{0.11}`$ for a sample of 109 dwarf irregular galaxies. Based on this evidence, we have set $`R_{HI}=2.3R_{25}`$ for the dwarf sample and determined all scale length values based on this conversion. The calculated values $`M_{dyn}`$ are listed in Table 2 along with the ratio of $`M_{HI}/M_{dyn}`$.

## 4 DISCUSSION

### 4.1 What is a dwarf galaxy?

Our understanding of dwarf galaxies, their properties, contents and star formation histories, has changed sharply, mostly because our catalogs of dwarf galaxies have expanded in the last decade. Our initial view of dwarf galaxies were based on our studies of the LMC, SMC and other Irr class galaxies in the Local Group (Holmberg 1957) and low luminosity dwarf spheriodals such as Leo A. These objects, although interesting in their own right, do not properly represent the range in mass, appearance and star formation history that define what we consider today to be the class of dwarf galaxies. For example, in the late 1950’s, de Vaucouleurs (1959) extended the late-type Hubble classes into the Sd, Sm and Im types to delineate the gradual loss of spiral structure. Most of these extreme late-type galaxies would later be defined as dwarfs based on their optical magnitudes (see Sandage and Binggeli 1984); however, they do not define the dwarf galaxy population.

Inspection of nearby clusters (Virgo and Fornax) revealed an extensive population of gas-poor elliptical dwarfs and gas-rich irregular dwarfs, subsequently called dE and dI to indicate their morphologically link to ordinary ellipticals and irregulars. The similarities between the dwarf ellipticals and irregulars to their giant counterparts (such as mass-metallicity relations, fundamental plane parameters and star formation) did appear to simplify the Hubble sequence as one divided by morphological class on one axis and by luminosity (and therefore mass) on the other. However, later work on field dwarfs (Eder et al. 1989, Schombert et al. 1997) showed that there is very little difference between the cluster dI class and the de Vaucouleurs Im class.

With the addition of a wide variety of morphological types, the definition of a dwarf galaxy has become blurred in the last decade. If asked, a theorist would define a dwarf galaxy by mass, then consider various formation scenarios that produce a mass spectrum. Unfortunately, dark matter becomes fractionally more significant in galaxies with decreasing mass (Ashman 1992), and this means that galaxy mass does not map, in a direct fashion, into observables such as luminosity, size or even rotational velocity. In addition, to make a difficult situation even worse, lower luminosity, density and size make the necessary observables laborious or impossible to obtain with our current technology and/or limited telescope time.

Within the observables, there are several methods in which to define a dwarf galaxy. For example, one can define a dwarf galaxy by optical luminosity (Tammann 1980). Since a galaxy’s optical luminosity is the primarily due to the photons emitted by stellar photospheres, then the total luminosity should reflect the stellar mass. With a stable IMF, the connection between luminosity and mass is parameterized by $`\mathrm{{\rm Y}}_{}=M_{}/L`$, the stellar mass, $`M_{}`$ divided by the total luminosity of a galaxy, $`L`$. Therefore, a dwarf galaxy could be defined by some boundary in optical luminosity as it reflects into stellar mass. However, since the fraction of dark matter (baryonic and non-baryonic) increases with later type galaxies (Ashman 1992), the use of luminosity to trace mass begins to decouple or has too many variables.

Dwarf galaxies are typically low in density, so division by central surface brightness (based on fits to the galaxy’s profile) can be used as an indicator of dwarfness. There is some physical reasoning to support the idea that all LSB galaxies should be dwarf systems. To produce low stellar densities would require very low past star formation rates and, hence, low total stellar mass and luminosity. Low surface brightness does not necessarily mean low stellar densities, but the colors of LSB galaxies rule out any highly unusual $`\mathrm{{\rm Y}}_{}`$ (McGaugh and de Blok 1997) implying that low optical density is indeed low stellar density. Or, conversely, we might expect a dwarf galaxy of low mass to have low mean density which, in turn, inhibits star formation to produce a LSB object. One could argue that initial conditions, such as the amplitude of the density perturbation in the early Universe from which the galaxy will form, produces a necessary connection between a dwarf galaxy and a LSB galaxy. Unfortunately, the discovery of large Malin galaxies destroyed the notion that all LSB galaxies are dwarfs. In fact, most LSB disk systems are comparable in size to the ordinary spirals that make up the Hubble sequence (McGaugh, Schombert and Bothun 1995).

Low mass implies instability, at least with respect to the pattern of star formation. So, irregular structure is nominally associated with dwarf galaxies (unless the star formation is completed, as in a dwarf elliptical). Basically, late-type galaxies can be broken down into three morphological types; those with bulges and LSB disks (Sc and Sd class, those with apparent axial symmetry but no bulges (Sm class) and those with no symmetry (Im class). However, morphology is only weakly related to the mass of a galaxy. There are many irregular galaxies that are massive (Hunter and Gallagher 1986) and many spirals which are dwarf-like in luminosity and size (Schombert et al. 1995). In addition, many high luminosity events, such as galaxy-galaxy interactions, are also associated with irregular morphology, although usually of a nature that identifies them as tidal in origin. Irregular morphology, by itself, is not a sufficient criteria for dwarf classification.

The most promising approach is to use the size of a galaxy as a measure of its dwarfness. A combination of angular size and irregular morphology were the criteria used in the initial construction of the PSS-II LSB dwarf catalog. Our previous LSB surveys (Schombert and Bothun 1988) was rich in LSB disks with scale lengths greater than 4 kpc, but few dwarf galaxies. That survey used a minimal angular size of one arcmin, the same as the UGC, to determine the bias against LSB galaxies in our catalogs. This PSS-II LSB dwarf survey set the criteria to be between 30 arcsecs and one arcmin, with a dramatically different metric size distribution as seen in Figure 3 of Schombert et al. (1997). With the acquisition of surface photometry and a H I redshift, a proper comparison of scale length sizes for the sample can be made (shown in the top panel of Figure 3). Most of the objects (90% of the sample) turn out to be quite small, less that the typical disk galaxy (3.2 kpc, de Jong 1996), testifying to the power of using morphology to select dwarf galaxies. We suspect that many of the undetected dwarf-like objects are probably background Sc’s outside the velocity range of the observing set-up for the Arecibo telescope. We conclude from Figure 3 that a robust method of collecting a sample of dwarf galaxies is to select objects by irregular morphology over a limited angular size (in our case, 30 to 60 arcsecs).

### 4.2 H I Profiles

A sample of 458 candidates were found by visual inspection of the PSS-II plates. Each candidate object was assigned a quality index as a measure of its probability to be a dwarf galaxy based on its morphological appearance. This is a purely subjective system, based on an objects irregular shape, lack of coherent spiral patterns, small to non-existent H II regions and a general LSB appearance. Most of these galaxies would be classified as extreme late-type (Sm, Im, Irr or dI) and a more complete discussion of their morphological and properties can be found in Schombert, McGaugh and Eder (2000). From our initial sample, 432 with high quality indices were observed on the Arecibo 300m telescope at 21-cm out to a velocity of 8,000 km sec-1.

Of the 432 galaxies observed, 259 (60%) were detected. The objects not detected are probably background LSB spirals (with spiral patterns that were not visible on the plates), gas-poor dwarfs (dE’s) or dwarfs whose H I masses and distance fall below the sensitive limit of the telescope. A small sample of the nondetections were reobserved in a velocity range of 4,000 to 12,000 km sec-1 in order to determine how many of the nondetected candidates were actually distant spiral galaxies. Of 18 candidates, 10 were detected at the higher velocities. However, only 4 of these have velocity widths which suggest a normal spiral galaxy (one being the Malin object D559-3 with a velocity width of 576 km sec-1). The loss of gas-poor dwarfs from the sample is unfortunate since the distribution of dE’s, outside of a cluster environment, is not known. In fact, the existence of a true dE, independent of a high local density or nearby companion, has not been confirmed (Binggeli, Tarenghi and Sandage 1990).

The spectra of all the detections from the three samples are presented in Figure 1. A hanning smoothing function has been applied, followed by a boxcar of 3. Each panel displays the raw data as well as the fitted baselines. A majority of the sample requires only a 2nd or 3rd order polynomial to flatten the baseline. As seen by other studies, most dwarfs have the classic double-horned profile indicative of rotation, although the data is too low in resolution to discern features such as asymmetries, filled-in cores or extended line wings (Hoffman et al. 1996, Matthews, van Driel and Gallagher 1998, Staveley-Smith, Davies and Kinman 1992). Skewed profiles may be due to confusion with nearby sources or pointing errors.

The H I profiles can be divided into two types, gaussian and double-horned. A galaxy with a flat or rising rotation curve, and a distribution of gas that declines with radius in the form of a power law, will exhibit a H I profile that has two peaks with a flat plateau region between them (see Giovanelli and Haynes 1988). The 246 spectra in Figure 1 were classified as either double horned, single horned or uncertain (listed in Table 1). Double horned profiles (such as D475-4 or D704-2) have distinguishable peaks or a plateau shape. Single horned profiles display a smooth gaussian shape with no indication of flattening. In many cases (e.g. D709-5), the profile is broad, but the double horn shape is unclear due to noise. In other cases (e.g. D702-1), there is an asymmetry that distorts the profile. Profiles with noisy plateaus are classed as uncertain. For the 246 spectra, 140 are double horned, 61 are single horned and 45 are uncertain. Assuming the uncertain profiles are evenly divided between single and double (most are broad and probably double horned), then double horned profiles are found in at least 70% of the sample.

The distribution of H I fluxes and H I profile widths are shown in Figure 2. All the detections are well above the 3$`\sigma `$ level given the mean fluxes and velocity widths. The mean velocity width, even without inclination corrections, is low compared to a typical late-type galaxy sample (see Schombert, Pildis and Eder 1997) confirming our intent to acquire a sample of low mass objects. A handful of objects with $`W_{50}>150`$ km sec-1 represent a few Malin type objects, irregular in their morphology, that mimic a dwarf galaxy’s appearance yet are actually large, distant spirals.

The metric scale length ($`\alpha `$) distribution, discussed in §4.1, is shown in the top panel of Figure 3. Based on exponential fits to the $`I`$ band surface brightness profiles, $`\alpha `$ provides the most fundamental property of a galaxy, its characteristic size. While the isophotal size of a galaxy, combined with its mean surface brightness, determines its visibility and detection probability. The typical scale length for a disk galaxy is indicated.

The redshift distribution of the detected dwarfs is shown in the bottom panel of Figure 3. The spatial distribution was shown in Figure 6 of Schombert, Pildis and Eder (1997). To repeat that paper’s conclusion, LSB dwarfs trace the same large scale structure as brighter galaxies. Figure 3 confirms that a majority of the dwarfs are located between 2000 and 6000 km sec-1, although there is a strong mass dependence. This dependence on mass is shown in the lower panel of the Figure 3 as a plot of H I mass versus redshift. The solid line shows the telescope 3$`\sigma `$ detection limit for a dwarf galaxy with a line width ($`W_{50}`$) of 75 km sec-1. A majority of the detections lie above this line. Objects below the limit all have line widths smaller than the assumed 75 km sec-1 and, thus, a higher S/N per channel. Note that the lowest mass object ($`M_{HI}=2.5\times 10^5M_{\mathrm{}}`$) is located at less than 500 km sec-1. This plot demonstrates that the search for very low mass H I-rich galaxies is extremely difficult as a combination of confusion with Galactic H I emission and small angular size makes their inclusion in any catalog problematic (see Schneider and Schombert 1999). It is also apparent from Figure 3 that a search for objects with $`M_{HI}<10^7M_{\mathrm{}}`$ will require a smaller angular limit in order to catalog candidates and deeper H I observations to obtain sufficient S/N (see also Zwaan et al. 1997) and that contrary to expectations (Briggs 1997) our galaxy catalogs are incomplete even in the local region of the Universe.

### 4.3 H I Masses and $`M_{HI}/L`$

The distribution of H I masses, calculated using the standard formula of distance and H I flux (Giovanelli and Haynes 1988) is shown in the top panel of Figure 4. The sample has a mean value of $`M_{HI}=7.9\times 10^8M_{\mathrm{}}`$ and a sharp cutoff at $`M_{HI}=4.0\times 10^9M_{\mathrm{}}`$. This can be compared to the UGC H I sample (see Figure 4, Schombert, Pildis and Eder 1997) where the mean H I mass of a UGC selected galaxy is $`M_{HI}=3.0\times 10^9M_{\mathrm{}}`$. Figure 4 shows that, whereas LSB galaxies are gas-rich, this does not mean that LSB dwarfs have high gas masses. The distribution of H I masses for this sample still defines the low-end of a gas mass distribution from any galaxy catalog and a galaxy selected to be dwarf-like by H I mass is as valid as classification by total mass or luminosity.

The bottom panel of Figure 4 displays mean H I surface density ($`\mathrm{\Sigma }_{HI}`$) for each dwarf with available optical imaging. The H I radius has been used to determine the average surface density (assuming the same axial ratio as the outer optical isophote). The $`\mathrm{\Sigma }_{HI}`$ values range from 0.6 to 20 $`M_{\mathrm{}}`$ pc-2 with a mean value of 5. The importance of $`\mathrm{\Sigma }_{HI}`$ to star formation was outlined in Kennicutt (1989) where he found that star formation follows a power-law in ordinary galaxies (the so-called Schmidt law, SFR $`\mathrm{\Sigma }_{HI}^\alpha `$). Key to Kennicutt’s findings was that there exists a threshold for star formation based on the Toomre (1964) instability criterion, the balance between gravity and rotation/thermal pressures. In a study of five dwarf galaxies, van Zee et al. (1997) found that the critical threshold was between $`9\times 10^{20}`$ and $`6\times 10^{21}`$ atoms cm-2, which corresponds to the range of 5 and 30 $`M_{\mathrm{}}`$ pc-2 in H I mass after corrections for neutral helium and other metals. This indicates that the mean H I surface densities for most of the LSB dwarfs in this sample are below the critical threshold density for star formation. Of course, there can exist local density enhancements that could form sites of star formation. However, the average surface density for the galaxy is below the threshold and this would explain the LSB nature to the dwarf sample, where star formation has occurred in the past (since there is some measurable stellar luminosity), but that this star formation has been inhibited and the current star formation levels are very low (McGaugh and de Blok 1997).

The fraction of a galaxy’s mass in the form of gas increases with later Hubble types and can dominate the kinematics of a dwarf galaxy (Meurer et al. 1996). In addition, the amount of gas in dwarf galaxies increases with decreasing surface brightness to the point where many of these dwarfs have gas fractions between 0.7 and 0.9 whereas the typical spiral have $`f_g0.4`$ (Schombert, McGaugh and Eder 2000). Thus, the Hubble sequence can be described as a declining ratio of stellar to gaseous material, where early-type galaxies are dominated by stars, later types by gas. This ratio is typically expressed as $`M_{HI}/L`$ and is shown in Figure 5 plotted against galaxy central surface brightness ($`\mu _o`$) and galaxy color ($`VI`$). The LSB dwarfs have a much higher value of $`M_{HI}/L`$ then the ordinary spirals, taken from the de Jong (1996). The distribution of $`M_{HI}/L`$ herein is similar to that of van Zee, Haynes and Giovanelli (1995), a study of extreme $`M_{HI}/L`$ galaxies from the Haynes and Giovanelli H I survey with the $`M_{HI}/L`$ values ranging from 0.3 to 12 and a mean value of 2.0. This average $`M_{HI}/L`$ found here is the same as the value found by van Zee, Haynes and Giovanelli (after adjustment for the $`I`$ band luminosities) and we also find very few galaxies with $`M_{HI}/L>10`$ even though this sample was selected by optical morphology rather than $`M_{HI}/L`$ ratio.

By itself, the H I properties of LSB dwarfs are poorly correlated with optical surface brightness or color. But, when compared with the de Jong sample, it is clear that the dwarfs continue the trend of increasing $`M_{HI}/L`$ with decreasing surface brightness and bluer colors. Both these correlations provide insight into the star formation history of dwarf galaxies. Higher $`M_{HI}/L`$’s with lower central surface brightness reflects the obvious connection between a galaxy’s gas supply and the density of its stellar population. LSB dwarf galaxies have converted very little of their gas into stars (i.e. their past star formation has been inefficient and/or at a very low rate). Thus, they maintain a high gas mass value in comparison to an even lower luminosity.

The correlation with galaxy color follows from the surface brightness correlation. Though star formation in LSB galaxies occurs at a low rate (McGaugh 1992), even the smallest amount of star formation has a large impact on the total galaxy color due to the low surface brightness of the underlying population. Galaxies with highest $`M_{HI}/L`$ values have the bluest $`VI`$ colors, meaning that the portion of the stellar population that contributes a majority of the optical luminosity in a LSB dwarf galaxy is quite young (see also Schombert, McGaugh and Eder 2000).

### 4.4 Dynamical Masses and Baryonic Matter

The dynamical mass values are listed in Table 2 for the dwarf galaxies with optical imaging and the histograms of $`M_{dyn}`$ and $`M_{dyn}/L`$ values are found in the of Figure 6. Note that this ratio is determined using $`I`$ band luminosities. For the typical dwarf color of $`BI=1.5`$, the $`M_{dyn}/L`$ values are a factor of 4 higher than similar $`B`$ values (Staveley-Smith, Davies and Kinman 1992). In general, the dynamical masses are a factor of 10 higher than the H I masses, but still much lower than a typical spiral sample (de Jong 1996). Again, as we have seen with indirect mass indicators (luminosity and size), galaxies selected by dwarf-like morphology define a low mass population. As shown in Figure 6, LSB dwarfs range in $`M_{dyn}/L`$ from 1.5 to 50 with a mean of 20. For comparison, the mean $`M/L`$ of spiral galaxy pairs is 40 at radii of 100 kpc (Honma 1999).

A comparison of dynamical mass to stellar mass (as represented by the total $`I`$ band magnitude of the dwarf) is found in Figure 7. Magnitudes in the $`I`$ band are preferred over other bandpasses because stellar population models (Worthey 1994) show that the stellar mass to luminosity ratio ($`M_{}/L`$) vary with time and star formation rate as a function of wavelength, but is most stable in the far-red. The use of the $`I`$ band observations minimizes these stellar populations effects due to its distance in wavelength from the region around the 4000Å break. Thus, $`I`$ band measurements 1) provide a more accurate estimate of the stellar mass of a galaxy, 2) obtain a luminosity measure that vary little with recent star formation and 3) determine structural parameters (such as scale length) which are undistorted by recent star formation events.

Although there is a great deal of scatter in the $`M_T^I`$ versus $`M_{dyn}`$ diagram, the LSB dwarf sample defines a linear correlation with a constant $`M/L`$ slope. Immediately obvious is that none of the dwarfs have their total masses composed solely of stars. Nor is the total mass explained by the sum of the stellar and gas masses (see below). In other words, LSB dwarfs are dark matter dominated. Interestingly, the linear correlation means that the dark matter follows the stellar light, even for this sample of extreme $`M/L`$ objects. One could interpret this one-to-one correspondence to indicate that dark matter is baryonic (i.e. baryons follow baryons). However, it is equally probable that the dark matter component is non-baryonic and simply accretes baryon material in a fashion that is proportional to its own mass. Certainly, with the discovery of massive neutrinos (see Turner 1997), it is now known that at least some fraction of the dark matter component in galaxies is non-baryonic.

The remaining question, regarding the dynamical mass estimates, is how much of a dwarf galaxies mass is in the form of observable baryonic material. This question is answered, graphically, in Figure 8 where we plot the ratios of stellar mass to dynamical mass and gas mass to dynamical mass. The stellar mass is calculated assuming a $`M_{}/L`$ of 1.2, based on the optical colors of LSB galaxies (McGaugh and de Blok 1997). Again, the advantage of $`I`$ band magnitudes is that this value is very insensitive to recent star formation, varying only 10% for the range of observed optical colors in LSB dwarfs. The gas mass is calculated from the H I mass and multiplying for a factor of 1.4 to account for helium and metals. There has been no observed CO emission in LSB galaxies (Schombert et al. 1990) and the current star formation signature by H$`\alpha `$ emission is weak. Thus, we take the amount of molecular and ionized gas to be negligible. Although there has been discussion of ionization by extragalactic UV (Corbelli and Salpeter 1993), none has been detected in LSB galaxies. The sum of the stellar and gas masses produces the mass in baryons ($`M_{baryons}`$) for each galaxy shown in Figure 8.

Unlike the rest of the Hubble sequence, LSB dwarfs tend to have equal amounts of mass in the form of stars and gas. However, both these components are a small fraction of the total mass. The mean for $`M_{stars}/M_{dyn}`$ is 0.11 and the mean for $`M_{gas}/M_{dyn}`$ is 0.08. This is similar to the value found by Burlak (1996) of $`M_{gas}/M_{stars}=0.7`$. Combining the gas and stellar mass numbers, we find that the mass in baryons for these dwarf galaxies is less than 20% the total mass, on average. Again, in agreement with a value of 26% from Burlak.

## 5 CONCLUSIONS

The paper is the third in a series concerning the properties of newly discovered LSB dwarfs. These galaxies were selected by irregular morphology from deep photographic plates (PSS-II). While the emphasis was on their visual appearance, this method to find new, uncataloged galaxies is automatically biased the search to objects with very low central surface brightnesses (i.e. not visible on previous sky surveys). The combination of low surface brightness and dwarf-like morphology was highly effective at producing a sample of low mass, low luminosity and small metric size galaxies. The results from the H I data presented herein can be summarized as the following:

* Low mass dwarfs exists, but are rare. Contrary to claims (Briggs 1997), our current galaxy catalogs are incomplete for objects with H I masses less than $`10^7M_{\mathrm{}}`$ since these objects are small and low in surface brightness and would be confused with Galactic emission in H I surveys.
* The H I mass distribution has a mean of $`M_{HI}=7.9\times 10^8M_{\mathrm{}}`$ with a sharp cutoff above $`M_{HI}=4.0\times 10^9M_{\mathrm{}}`$ and a long tail to low H I masses. This distribution is much lower than a typical spiral sample and demonstrates that, although LSB galaxies are H I-rich, the gas component does not make them high mass. Mean H I surface densities are below the critical density for star formation based on the Toomre criteria, justifying statements that LSB dwarfs are inefficient at star formation and have quiescent past histories in agreement with their LSB nature.
* LSB dwarfs continue the trend from late-type spirals of high $`M_{HI}/L`$, increasing with lower central surface brightnesses and bluer colors. The explanation of this correlation lies in a quiescent star formation history which leaves a large gas reservoir on top of a faint stellar population. Any recent star formation, even at a very low rate, has a sharp blueward effect on the underlying LSB structure.
* The $`M_T^I`$,$`M_{dyn}`$ correlation (Figure 7) demonstrates that the LSB dwarf sample that the dark matter (whether baryonic or non-baryonic) follows the baryonic matter. The distribution of $`M_{dyn}/L`$ ranges from 1.5 to 50 (note that these are $`I`$ band numbers, a factor of 4 higher than equivalent values in the $`B`$ band).

###### Acknowledgements.

The authors wish to thank the generous support of the Arecibo Observatory for the allocation of time to search for HI emission from the candidate dwarf galaxies and Michigan- Dartmouth-M.I.T Observatory in carrying out the photometry portion of this program This work is based on photographic plates obtained at the Palomar Observatory 48-inch Oschin Telescope for the Second Palomar Observatory Sky Survey which was funded by the Eastman Kodak Company, the National Geographic Society, the Samuel Oschin Foundation, the Alfred Sloan Foundation, the National Science Foundation and the National Aeronautics and Space Administration.


format: text

text:
# On the convergence of Regge calculus to general relativity

## 1 Introduction

Regge calculus  is a discrete theory of gravity which replaces the smoothly curved spacetime of general relativity with a lattice. The curvature of the lattice spacetime is concentrated entirely on the two-dimensional hinges of the four-dimensional lattice cells.

Regge calculus holds much promise for the numerical investigation of both classical and quantum gravity. Although the lattice approach appears well suited to numerical applications, progress in the field has been slow. Only recently have the first completely generic four-dimensional numerical simulations been performed .

After the recent papers by Brewin  and M. Miller , and despite the proven track record of the Regge calculus, a lively debate arose  as to whether or not solutions of the Regge equations would converge to solutions of the Einstein equations in some suitable limit. Neither Brewin or M. Miller directly computed solutions of the Regge equations. Instead, they took the somewhat easier approach of evaluating the Regge equations on an exact solution of Einstein’s equations. They did this using sophisticated interpolation schemes, based on geodesics, to map a range of Einstein solutions onto simplicial lattices. The residual of the resultant Regge equations, calculated using these interpolated lattice edge lengths, was then examined in the limit of very fine lattice discretisations. Brewin observed that the residual scaled as $`𝒪(1)`$ as the lattice was refined, and he inferred that solutions of the Regge equations, in generic spacetimes, would not converge to solutions of the Einstein equations in the limit of fine discretisations. M. Miller chose a slightly different measure for the residual of the Regge equations, but when expressed in terms of Brewin’s measure, M. Miller’s results are consistent with those obtained by Brewin. Despite this M. Miller does draw sharply different conclusions, based upon further calculations associated with averages of the Regge equations, namely that solutions of the standard Regge equations will converge as the fourth order in the lattice spacing to solutions of the Einstein equations. Similar behaviour, in which the residual converges only over averages of the equations, has been observed by Sorkin  in the context of massless scalar fields on a 2-dimensional simplicial space.

The situation became even more confused after the work of Gentle and W. Miller , who demonstrated explicit quadratic convergence of particular solutions of Regge calculus to a solution, the Kasner spacetime, of general relativity. This spacetime was one of the test cases used by Brewin and M. Miller. This seems most odd – we appear to have a convergent set of solutions from an apparently non-convergent set of equations. How can this be?

In this paper we explore a possible explanation for this behaviour, as proposed by Brewin , which is consistent with all previous numerical studies. It is important to note that the direct solution of the Regge equations has always generated solutions which converge to the corresponding solutions of general relativity. Although the observed rate of convergence is dependent upon the particular lattice and symmetry restrictions imposed, when the lattice construction allows the full expression of the gravitational degrees of freedom, second order convergence of the lattice solutions to the continuum has always been observed .

We begin with a more complete description of the problem in section 2, and then discuss a possible explanation in section 3. Finally, in section 4, we present numerical evidence to support our proposal.

## 2 Investigating the convergence of Regge calculus

There are various ways in which the convergence properties of a numerical technique can be explored. The most rigorous approach involves the direct comparison of the approximate equations and the full system, concentrating on the leading order discrepancies. The underlying assumption is that the solutions of the approximate equations will converge at the same rate as the approximate equations themselves converge to the original system.

The Regge and Einstein equations are too complex for such a direct approach to be beneficial. Not only are the equations inherently complex non-linear systems, but before such an analysis can proceed one must choose a lattice upon which to express the Regge equations, together with a coordinate system in which to write out the Einstein equations. Once these decisions have been made, we are faced with a more problematic choice: how are the two sets of equations to be compared? In general there are more Regge equations (one per lattice edge) than Einstein equations (ten per spacetime event).

There is no clear prescription for averaging the lattice equations to obtain the correct number of Einstein equations in the continuum limit. Brewin  and M. Miller  have both developed their own schemes for averaging the Regge equations. A direct assault on the convergence properties of Regge calculus is therefore likely to be applicable only on a specific lattice, with a particular choice of averaging in the continuum limit. That is, the results would indicate that, in some averaged sense, the Regge equations converge to the Einstein equations. However, the result would not constitute a proof valid for all possible lattice choices.

Given these limitations, it is natural to ask if similar results can be obtained by directly solving the Regge equations and comparing the lattice solutions to corresponding Einstein spacetimes. Brewin and M. Miller have investigated the inverse problem: given a known (analytic) solution of general relativity, and interpolating that solution onto a lattice, what can be said about the convergence of the Regge equations?

The approach taken by Brewin and M. Miller was to introduce some discretisation process in the smooth manifold on which the Einstein equations are defined. This provides a way to map the smooth metric solutions to a lattice, giving a set of lattice edge lengths $`𝐋_E`$ derived from the continuum metric. Both authors used geodesic lengths to map continuum information onto the lattice spacetime. These new “continuum” lattice edge lengths do not in general satisfy the Regge equations, but we can evaluate the residual

$$r=|𝐑(𝐋_E)|$$
(1)
which is an indication of how well the interpolated Einstein solution satisfies the Regge equations. It is this residual, using an appropriately chosen norm, which both Brewin and M. Miller considered .

Brewin observed that the residual of the simplicial equations remained roughly constant as the lattice was refined on a fixed region of spacetime. This result led Brewin to question the validity of Regge calculus, since any useful numerical scheme must converge to the underlying solution of the partial differential equations as the resolution is improved. M. Miller observed second order convergence of the residual, for any smooth metric, whether or not they were solutions of Einstein’s equations. Had he used the same norm for the residual as used by Brewin we believe he too would have observed $`𝒪(1)`$ convergence of the residuals.

It is important to note that the experiments of Brewin and M. Miller do not directly evaluate the convergence of the numerical solutions. Rather, they investigate the convergence of the lattice equations to the Einstein equations. As the lattice is refined, it is reasonable to expect that the interpolated Einstein solution will satisfy the Regge equations increasingly accurately. This, however, is not what Brewin observed numerically.

## 3 A possible explanation

The observations of Brewin are particularly puzzling in light of the many previous applications of Regge calculus. Almost every numerical application of the method has displayed convergence towards the corresponding continuum solution, with most studies indicating that numerical Regge calculus is a second order accurate approximation to general relativity.

The edge lengths measured in the Regge lattice ($`L`$) and the corresponding interpolated Einstein edges ($`L_E`$, obtained by assigning geodesic lengths calculated in the continuum) are related as

$$L=L_E+𝒪(\delta ^{p+1})$$
(2)
with previous numerical studies suggesting that $`p=2`$ (see Gentle  for case studies and a general review). Throughout this paper we assume that $`\delta `$ is a typical length scale in the lattice, and noting that the edges themselves are of this magnitude, we say that edges which satisfy equation (2) are $`p`$-order accurate approximations to the continuum solution.

Brewin  has proposed a mechanism whereby the solutions of the Regge equations converge to the corresponding Einstein solutions, while at the same time the residual of the Regge equations evaluated on the interpolated Einstein lattice edges does not converge. The key to Brewin’s proposal is allowing the functional form of the error terms (the discrepancy between the Einstein and Regge solutions) to depend on the discretisation scale $`\delta `$. This can be clearly seen from a toy model.

Suppose $``$ is a second order differential operator, and $`y(x)`$ is a solution of the equation

$$y=0.$$
(3)
Brewin considers the function

$$\stackrel{~}{y}_\delta (x)=y(x)+\delta ^2f(x/\delta )$$
(4)
for some arbitrary scalar $`\delta `$ and an arbitrary, though bounded, function $`f(x)`$. Considering the difference between the solutions we find that

$$|y(x)\stackrel{~}{y}_\delta (x)|=𝒪(\delta ^2),$$
(5)
indicating that the solutions differ only by “second order” terms. Noting that $`\stackrel{~}{y}_\delta (x)`$ is a solution of some other related equation $`\stackrel{~}{}_\delta `$,

$$\stackrel{~}{}_\delta \stackrel{~}{y}_\delta =0,$$
(6)
where $`\stackrel{~}{}_\delta `$ is also a second order differential operator (though different from $``$), we find that the “residual” of $`\stackrel{~}{y}_\delta (x)`$ with respect to the original operator is

$$|\stackrel{~}{y}_\delta |=𝒪(1).$$
(7)
This toy model embodies precisely the properties observed by Brewin – second order convergence of the solutions , with no corresponding convergence observed in the Regge equations when they are evaluated on exact solutions interpolated from the continuum . The discrepancy between the two solutions in this toy model is seen to be a wave-like disturbance with frequency proportional to $`1/\delta `$.

In the continuum limit it is reasonable to expect that the discrete Regge equations (or a weighted average over them) approach a system of differential equations. Moreover, consideration of the Einstein equations leads us to expect that the limiting form of the Regge equations is a set of second order non-linear equations. Furthermore, previous numerical experiments suggest that Regge calculus is a second order method; we expect $`p=2`$. The toy model considered above leads us to expect that if the lattice solutions differ from the continuum solutions by terms with frequencies proportional to $`1/\delta `$, it is not unreasonable that the residuals of the Regge equations remain roughly constant as the resolution is improved.

This is precisely the behaviour observed by Brewin ; the solutions converge even though the residual of the equations do not. This explanation relies on the existence of high frequency, low amplitude waves in the simplicial solutions; a possibility not ruled out by any of the largely low-resolution applications of Regge calculus to date. In fact, one study hints at the existence of precisely this type of wave-like structure in the Regge solutions . In the next section we construct high resolution solutions of the Regge equations in order to gain insight into the fine-scale behaviour of the lattice solutions.

## 4 Numerical solution of the Regge equations

In this section we solve the Regge equations for the vacuum Kasner cosmology using a $`(3+1)`$-dimensional formulation of Regge calculus described elsewhere . The initial value problem is solved and the lattice is evolved subject to simplicial lapse and shift conditions.

The Kasner cosmology is a homogeneous anisotropic exact solution of the vacuum Einstein equations, and is the prototype for generic velocity dominated cosmological singularities. It also provides a model for epochs in the mixmaster cosmology between bounces. The metric may be written in the form

$$ds^2=dt^2+t^{p_1}dx^2+t^{p_2}dy^2+t^{p_3}dz^2,$$
(8)
where the field equations reduce to a pair of algebraic constraints,

$$p_1+p_2+p_3=p_1^2+p_2^2+p_3^2=1.$$
(9)
The Kasner exponents $`(p_1,p_2,p_3)`$ are a one-parameter set, and the solutions include the flat-spacetime case $`(0,0,1)`$.

The previous application of Regge calculus to the Kasner cosmology was viewed as a benchmark of the numerical technique, and as such the convergence of the solutions was an important issue. It was found that the lattice solution converged to the continuum as the second power of the typical lattice spacing . This earlier study used a fixed lattice resolution, containing 128 vertices arranged in two offset cubic grids each consisting of $`4\times 4\times 4`$ vertices. Reducing the overall scale of the lattice, combined with the $`T^3`$ topology of the Kasner solution, allowed the convergence rate to be estimated.

In this paper we directly subdivide the lattice whilst fixing the size of the spatial region. This is a more general and demanding convergence test, as any long wavelength modes excluded in the previous study are now resolvable. Since the major issue is convergence we refine the lattice in only one dimension, but reduce the scale of the remaining spatial axes to avoid the creation of long skinny triangles. Such long and skinny triangular elements are known to cause instabilities in finite element calculations.

This effectively “one-dimensional” model allows convergence analysis to be conducted whilst containing the computational scale of the problem. All calculations are performed on a lattice containing two offset cubic grids which are fully subdivided into triangles, tetrahedra and four-simplices. Each cubic grid contains $`2^n\times 4\times 4`$ vertices ($`N2^n=4,8,\mathrm{},1024`$), with the $`x`$-axis having a fixed length of $`X=10`$, and the remaining axes scaled to compensate; $`Y=Z=40/2^n`$. This ensures that the spatial edges surrounding a vertex are all of the same order of magnitude. The typical spatial resolution in each model is thus $`\delta 10/2^n`$.

The full four-dimensional two-slice initial value problem is solved at $`t=1`$, after which the lattice is evolved to $`t=2`$ using the $`(3+1)`$-dimensional Sorkin evolution scheme and “geodesic slicing” conditions; that is, with unit lapse and zero shift. The time step is chosen for each resolution to give a Courant factor of $`dt/l0.2`$; for $`N=512`$ we chose $`dt=3.9\times 10^3`$. However, the main results of this paper were found to be largely insensitive to the choice of timestep, provided it satisfied the Courant condition. For full details of the evolution and initial value algorithms, see Gentle and W. Miller . The results of the Regge evolution are compared with the continuum solution by performing a least-squares fit of the function

$$L(t)=L_0t^{p_r}$$
(10)
to the time evolution of each class of the axis-aligned spatial edges, $`L_x,L_y`$ and $`L_z`$. In agreement with previous studies, the lattice solution is found to remain spatially homogeneous throughout the evolution to within roundoff error. The typical standard deviation of the edges also remained below the limit of numerical accuracy.

Figure 1 displays the fractional difference between the early evolution of the Regge cosmology and the exact solution for the axisymmetric case with $`p_1=p_2=2/3,p_3=1/3`$. Although the average magnitude of the discrepancy is small, there is a clear oscillation superimposed upon a general trend. The figure displays data for three lattices consisting of $`128\times 4\times 4`$, $`256\times 4\times 4`$ and $`512\times 4\times 4`$ vertices.

An estimate of the rate at which the Regge solution converges to the continuum is given in figure 2. After performing a least squares fit of equation (10) to the simplicial solution in the region $`t=[1,2]`$, we plot the fractional difference between the fitted power $`p_r`$ and the analytic value $`p_c`$ for each spatial axis. It is clear that this power law fit approaches the continuum solution as the second power of the lattice spacing. Figure 3 shows the scaled error in the least squares fit, giving an indication of how well the power law (10) describes the Regge solution. Again, we find second order convergence.

The data represented in figure 3 contains additional information. It is clear from figure 1 that there is a strong wave-like oscillation superimposed upon the general error curve. Fitting equation (10) to such data will draw out the averaged solution, while the error in the fit will reflect the magnitude of any discrepancies from that average behaviour. Thus we see that both the oscillations and the underlying trend in figure 1 converge to zero as at least the second power of the grid spacing.

In light of the discussion in the previous section, the frequency of the oscillations evident in figure 1 is also of interest. It is clear from the figure that the frequency increases with the lattice resolution. Figure 4 shows this in a more quantitative manner, displaying the Fourier transform of a small segment of the fractional error. It is apparent that the frequency of the waves varies linearly with the number of vertices, or alternatively, that the wavelength of the oscillations is proportional to $`1/N`$.

Together, these $`(3+1)`$-dimensional convergence results show that the numerical Regge solution is a second order accurate approximation to the corresponding solution of the Einstein equations. Moreover, the leading order difference between the two solutions contains high frequency, bounded oscillations. We have shown that the frequency of these waves is proportional to $`N`$ (or $`1/\delta `$), while their magnitude reduces as $`1/N^2`$ ($`\delta ^2`$). This is precisely the variation postulated by Brewin , and strongly supports his explanation of the apparent non-convergence of the Regge equations.

Finally, we note that the Kasner solution displays only temporal oscillations, with the spatial three-geometries remaining homogeneous to high accuracy throughout the evolution. This is likely a result of the high degree of symmetry of the Kasner spacetime, in which the constant time spatial hypersurfaces are flat. In a more general setting we would expect to see similar high frequency, bounded oscillations in both space and time.

## 5 Conclusion

We have carried out a rigorous convergence study of a particular solution to the Regge equations, and shown that the solutions do indeed converge to the corresponding Einstein solutions in the limit of very fine discretisations.

Following a suggestion of Brewin we examined the behaviour of the error terms, and found convincing evidence for the existence of high frequency, low amplitude oscillations superimposed upon the continuum solution. Brewin previously suggested that waves of precisely this form in the simplicial solution could explain the lack of convergence observed in the residual of the Regge equations when they are evaluated on the interpolated continuum solution.

Together these results suggest that solutions of the Regge equations are generally second order accurate approximations to the corresponding Einstein spacetimes, with the discrepancy between the two solutions consisting of high frequency, low amplitude waves. Although these waves prevent the residual of the Regge equations converging to zero when evaluated on interpolated Einstein solutions , they do not affect the overall second order accuracy of the simplicial solutions.

We expect that all generic Regge simulations in vacuum will produce similar results. That is, all simulations will contain high frequency oscillations, linked to the inverse of the discretisation scale, which reduce in magnitude as the second power of that scale. These oscillations do not appear to induce instabilities in the numerical evolution of simplicial lattices. Indeed, previous studies have indicated that the oscillations gradually decay as the evolution proceeds . However, M. Miller  has reported instabilities arising in linearised Regge calculus on asymmetrical grids.

## Acknowledgements

We are grateful to Todd Lane for assistance with the preparation of the figures, and the Australian Research Council for a Small Grant to support this work.

## References


text:
# Consequences of Interstellar Absorption for Models of Anomalous X-ray Pulsars

## 1 Introduction

In the past few years, $`X`$-ray observations have revealed a new class of $`X`$-ray sources, the so-called anomalous $`X`$-ray pulsars (AXPs) (Mereghetti & Stella 1995; van Paradijs, Taam & van den Heuvel 1995). AXPs are sources of pulsed $`X`$-ray emission, with persistent luminosities $`L_x10^{35}10^{36}`$ erg $`\mathrm{s}^1`$ and soft spectra. Their periods lie in a very narrow range, between 6 and 12 seconds, and their characteristic ages are of order $`10^310^5`$ years. Their distances are usually very large, as inferred from the typically large column densities to them ($`N_\mathrm{H}10^{22}`$ cm-2).

Due to the long period of these objects, it is clear that their $`X`$-ray luminosity cannot be rotation-powered. Two broad classes of models have emerged to explain their $`X`$-ray emission. In one, the radiation is assumed to be powered by accretion from binary companions of very low mass (Mereghetti & Stella 1995), from the debris of a disrupted high mass companion (van Paradijs, Taam & van den Heuvel 1995; Ghosh, Angelini & White 1997), from the interstellar medium (Wang 1997), or from the debris falling back after the supernova explosion (Chatterjee, Hernquist & Narayan 2000; Alpar 1999, 2000; Marsden et al. 1999). In the second class of models, AXPs are hypothesized to be isolated, ultramagnetized neutron stars, with field strengths in the range $`10^{14}10^{15}`$ G, a.k.a. “magnetars” (Duncan & Thompson 1992). The $`X`$-ray luminosity could then be powered by either magnetic field decay (Thompson & Duncan 1996; Heyl & Kulkarni 1998) or by residual thermal energy (Heyl & Hernquist 1997a,b).

There are several observational constraints that can be used to discriminate among various models. The presence of an optical counterpart and observable Doppler shifts in the frequency of the $`X`$-ray pulses would suggest the presence of an optical counterpart, but this has not been observed so far. Optical and longer wavelength emission (Hulleman et al 2000; Perna, Hernquist & Narayan 2000) would indicate the presence of an accreting fossil disk. Any viable model, however, must also be consistent with the observed high pulse fractions and small emitting areas inferred for the $`X`$-ray radiation from these objects.

In this Letter we examine whether the observed pulse fractions and inferred emitting areas are consistent with a model for AXPs in which their $`X`$-ray emission is powered by thermal radiation from highly magnetized neutron stars. We show that the amplitudes of the pulsations for the total flux observed in a given band are strongly sensitive to the amount of absorption to the source, and, more specifically, increase with it. To understand this point, one first needs to notice that, for typical local blackbody (BB) emission, the energy-dependent pulse fractions increase with energy, as a consequence of the fact that the hardness of the blackbody spectrum becomes larger at higher energies. This point was well elucidated by Page (1995) and Page & Sarmiento (1996): the ratio of the BB fluxes emitted at energy $`E`$ by two regions of areas $`A_1`$ and $`A_2`$ with respective temperatures $`T_1`$ and $`T_2`$ is

$$\frac{F_{\mathrm{BB}}(E,T_1)}{F_{\mathrm{BB}}(E,T_2)}=\frac{A_1}{A_2}\frac{\mathrm{exp}(E/kT_2)1}{\mathrm{exp}(E/kT_1)1},$$
(1)
which is an increasing function of $`E`$ if $`T_1>T_2`$. This property remains valid for atmospheric models with beaming $`(\mathrm{cos}\theta )^n`$, and, furthermore, the increase of the pulse fractions with energy becomes stronger with increasing $`n`$.

Absorption by the intervening interstellar medium (ISM) does not affect the energy dependent pulse fractions (though this is rigorously true only for a perfect detector), as photons emitted at a given energy are reduced by the same amount independently of the phase of the rotation. However, pulse fractions measured from flux which is integrated over a given band do depend on absorption. Absorption, in fact, reduces the low-energy part of the flux with respect to the high-energy one. As a result, the main contribution to the pulse fractions is weighed more and more by the high energy component as the amount of absorption is increased.

We compute the expected pulse fractions and inferred areas for typical parameters of AXPs and a thermal cooling model for their emission, which includes processing through an atmosphere. We show that the inferred properties are consistent with observations.

## 2 Model

Our model is that of a neutron star (NS) cooling through an accreted envelope and with an intense, dipolar magnetic field. Heyl & Hernquist (1998b) showed that, for $`B_p10^{14}`$ G, the flux transmitted through the envelope is well described by $`F\mathrm{cos}^2\psi `$, where $`\psi `$ is the angle between the radial direction and the magnetic field. For a dipolar field, $`\mathrm{cos}^2\psi =4\mathrm{cos}^2\theta _p/(3\mathrm{cos}^2\theta _p+1)`$ (Greenstein & Hartke 1983), with $`\theta _p`$ being the magnetic colatitude (angle between radial direction and magnetic pole). It is given by

$$\mathrm{cos}\theta _p=\mathrm{cos}\theta \mathrm{cos}\alpha +\mathrm{sin}\theta \mathrm{sin}\alpha \mathrm{cos}\varphi ,$$
(2)
where $`\alpha `$ is the angle that the magnetic pole makes with the line of sight, and $`(\theta ,\varphi )`$ are the usual colatitude and longitude in spherical coordinates. In our calculation, we assume $`\alpha =90^{}`$ (orthogonal rotator). Therefore our results on the pulsation amplitudes must be interpretated as upper limits. We assume that the local emission, $`n(E,T)`$, is a blackbody spectrum modified by the presence of an atmosphere, for which we adopt the semianalytical model of Heyl & Hernquist (1998a), but with the inclusion of limb darkening 111The model is that of an accreted atmosphere with a magnetized iron envelope, but note that the composition of the envelope does not really affect any of the points made here.. This leads to a beaming of the radiation $`\mathrm{cos}\delta `$, with $`\delta `$ being the angle between the normal to the NS surface and the direction of the photon trajectory. We will also consider cases with a more intense beaming $`(\mathrm{cos}\delta )^n`$, but not extreme (i.e. we restrict our attention to $`n2`$). Static general relativistic effects (i.e. light deflection and gravitational redshift) are fully taken into account in our calculations.

Let us define $`e^{\mathrm{\Lambda }_s}=\sqrt{1R/R_s}`$, where $`R`$ is the radius of the NS star, and $`R_s=2GM/c^2`$ its Schwarzschild radius. Here we assume $`M=1.4M_{}`$, and consider the range of radii $`2(R/R_s)4`$, all compatible with the currently available models for the NS equation of state. Let $`D`$ be the distance from the star to the observer, and $`N_\mathrm{H}`$ the intervening column density. The flux measured by an observer at infinity is then given by

$$f(E;N_\mathrm{H})=\frac{\pi R_{\mathrm{}}^2\sigma T_{p,\mathrm{}}^4}{4\pi D^2}\frac{1}{kT_{p,\mathrm{}}}e^{\sigma (E)N_\mathrm{H}}_0^12x𝑑x_0^{2\pi }\frac{d\varphi }{2\pi }I_0(\theta ,\varphi )n[Ee^{\mathrm{\Lambda }_s};T_s(\theta ,\varphi )],$$
(3)
in units of phot cm-2 s-1 keV-1. Here $`x=\mathrm{sin}\delta `$, $`R_{\mathrm{}}Re^{\mathrm{\Lambda }_s}`$, and $`T_{p,\mathrm{}}T_pe^{\mathrm{\Lambda }_s}`$, where $`T_p`$ is the temperature at the pole. The general relativistic effects of light deflection are taken into account through the ray-tracing function $`\theta =\theta (x)`$, as described in Page (1995). The function $`I_0(\theta ,\varphi )`$ sets the normalization of the flux at each point $`(\theta ,\varphi )`$ of the NS surface

$$I_0(\theta ,\varphi )=\frac{4\mathrm{cos}^2\theta _p}{3\mathrm{cos}^2\theta _p+1}(\mathrm{cos}\theta _p)^{0.4},$$
(4)
where, following Heyl & Hernquist (1998b), we have assumed a further dependence of the flux on $`B^{0.4}`$. The local temperature on the star is given by $`T_s(\theta ,\varphi )=T_p[I_0(\theta ,\varphi )]^{1/4}`$, while for the photoabsorption cross section $`\sigma (E)`$ we use the analytical fit provided by Morrison & McCammon (1988).

Let now $`E_{\mathrm{min}}`$ and $`E_{\mathrm{max}}`$ define the boundaries of the band of observation, and let us define $`F(N_\mathrm{H})=_{E_{\mathrm{min}}}^{E_{\mathrm{max}}}𝑑Ef(E;N_\mathrm{H})`$. Then, if the column density to the source is $`N_\mathrm{H}`$, the pulse fraction in that given band is

$$pf(N_\mathrm{H})=\frac{F^{\mathrm{max}}(N_\mathrm{H})F^{\mathrm{min}}(N_\mathrm{H})}{F^{\mathrm{max}}(N_\mathrm{H})+F^{\mathrm{min}}(N_\mathrm{H})}.$$
(5)
In Figure 1, we plot this function for fixed values of $`T_pT_{eff}(\psi =0)`$ and $`n`$, while exploring the dependence on $`R/R_s`$. Our results agree with those of DeDeo, Psaltis & Narayan (2000) in the limit of $`N_\mathrm{H}0`$. However, for large column densities, our plots show that the pulsation amplitudes can significantly increase. Therefore, neglecting this effect can lead to misleading interpretations.

In Figure 2, we consider a typical range of column densities inferred for AXPs and show, for an intermediate value of $`R/R_s`$, the dependence of pulse fractions on temperature and degree of beaming. We have considered a range of temperatures that is typical for AXPs, if one accounts for the fact that processing through the atmosphere typically results in an overestimate of the effective temperature by a factor of $`23`$. Note that the pulsation amplitudes in a given energy band increase as $`T_p`$ decreases. This is due to the fact that, as $`T_p`$ decreases, the spectrum is shifted towards lower energies. This in turn leads, for the same amount of absorption, to a larger suppression of softer photons, and hence to an increase of the pulse fractions. Figure 2 clearly shows that a thermal model for highly magnetized NSs can produce pulsation amplitudes which are consistent with those observed so far (see Chakrabarty et al. 2000 for a summary).

Another property that must be accounted for by any viable model for AXPs is the inferred size of the emitting region. To estimate the effective emitting areas that an observer would deduce for an highly magnetized, cooling neutron star as in our model, we generated a spectrum using Equation (3). We then used the software XSPEC to convolve it with the instrument response (we used CHANDRA) and generate a set of simulated data. These were then fitted using a blackbody function. For a star of radius $`R=\{2,3,4\}R_s`$ (corresponding to $`R=\{8.25,12.37,16.4\}`$ km), we found, respectively222These values are not very dependent on $`T_p`$ and $`n`$, within the range of values considered for the AXPs., $`R_{\mathrm{}}`$ {1.3,2.3,3.5} km, significantly smaller than the stellar radius, but comparable to the parameters inferred for AXPs.

## 3 Conclusions

We have examined the properties of a model in which the $`X`$-ray emission of AXPs is thermal radiation from highly magnetized cooling neutron stars. We have shown that the amplitudes of the pulsations in a given energy band are an increasing function of the column density to the source. If the typically high column densities inferred for AXPs are taken into account, then even a moderate beaming can lead to pulse fractions comparable to those observed so far.

We have used the software XSPEC to generate simulated data from our model and then used these data to find the best-fit blackbody spectrum. We find that the best fit implies an effective radius of the star which is typically $`56`$ times smaller than the actual radius. This is again in agreement with the radii of the emitting regions inferred from the observations of AXPs.

Our results have clear implications for efforts to discern the nature of a variety of X-ray point sources. AXPs have variously been interpreted either as accreting sources, or as isolated, strongly magnetized neutron stars. If interstellar absorption is not taken into account, then the large observed pulse fractions appear to be inconsistent with those predicted by thermal cooling models, as found in recent work by DeDeo et al. (2000). However, as we have demonstrated, interstellar absorption is able to boost the pulsed fraction measured by a distant observer relative to that which would be seen otherwise to the point that a thermal cooling model can easily accommodate the pulsed fractions associated with observed AXPs, if the high inferred column densities to these sources are properly taken into account. Consequently, we do not believe that pulsed fractions discriminate between thermal models for AXP emission from those on which the X-ray luminosity is accretion-powered.

Likewise, interstellar absorption complicates efforts to constrain the sizes of thermally emitting areas on the surfaces of neutron stars. The spectrum of the X-ray point source discovered at the center of the Cas A supernova remnant (Tannanbaum 1999) has now been discussed separately by a number of groups (Pavlov et al. 1999; Umeda et al. 1999; Chakrabarty et al. 2000). Blackbody fits to the observed spectra yield emitting areas that are only a small fraction of a neutron star surface. These results cannot be explained by other models which interpret the spectra using hydrogen atmospheres, but which ignore both the non-uniform temperature distribution expected for a strongly magnetized neutron star and the consequences of interstellar absorption. As we have indicated, when all these various effects are accounted for, there does not appear to be a discrepancy between the effective emitting area of a strongly magnetized neutron star and the values which have been derived for objects like Cas A and the AXPs.


text:
# Monte Carlo Simulations of Globular Cluster Dynamics

## 1. Introduction

The first Monte Carlo methods for calculating the dynamical evolution of star clusters in the Fokker-Planck approximation were developed more than 30 years ago. They were first used to study the development of the gravothermal instability (Hénon 1971a,b; Spitzer & Hart 1971a,b). More recent implementations have established the Monte Carlo method as an important alternative to direct $`N`$-body integrations (see Spitzer 1987 for an overview, and the article by Giersz in this volume). The main motivation for our recent work at MIT was our realization a few years ago that the latest generation of parallel supercomputers now make it possible to perform Monte Carlo simulations for a number of objects equal to the actual number of stars in a globular cluster (in contrast, earlier work was limited to using a small number of representative “superstars,” and was often plagued by high levels of numerical noise). Therefore, the Monte Carlo method allows us to do right now what remains an elusive goal for $`N`$-body simulations (see, e.g., Aarseth 1999, and the article by Makino in this volume): perform realistic, star-by-star computer simulations of globular cluster evolution. Using the correct number of stars in a dynamical simulation ensures that the relative rates of different dynamical processes (which all scale differently with the number of stars) are correct. This is particularly crucial if many different dynamical processes are to be incorporated, as must be done in realistic simulations (cf. the article by Heggie in this volume).

Our implementation of the Monte Carlo method is described in detail in the papers by Joshi, Rasio, & Portegies Zwart (2000), Joshi, Nave, & Rasio (2000), and Joshi & Rasio (2000). We adopt the usual assumptions of spherical symmetry (with a 2D phase space distribution function $`f(E,J)`$, i.e., we do not assume isotropy) and standard two-body relaxation in the weak scattering limit (Fokker-Planck approximation). In its simplest version, our code computes the dynamical evolution of a self-gravitating spherical cluster of $`N`$ point masses whose orbits in the cluster are specified by an energy $`E`$ and angular momentum $`J`$, with perturbations $`\mathrm{\Delta }E`$ and $`\mathrm{\Delta }J`$ evaluated on a timestep that is a fraction of the local two-body relaxation time. The cluster is assumed to remain always very close to dynamical equilibrium (i.e., the relaxation time must remain much longer than the dynamical time). We have performed a large number of test calculations and comparisons with direct $`N`$-body integrations, as well as direct integrations of the Fokker-Planck equation in phase space, to establish the accuracy of our basic treatment of two-body relaxation (Joshi et al. 2000a). Fig. 1 shows the results from a typical comparison between Monte Carlo and $`N`$-body simulations. Our main improvements over Hénon’s original method are the parallelization of the basic algorithm and the development of a more sophisticated method for determining the timesteps and for computing the two-body relaxation from representative encounters between neighboring stars. Our new method allows the timesteps to be made much smaller in order to resolve the dynamics in the cluster core more accurately.

## 2. Summary of Recent Results

Our recent work has focused on the addition of more realistic stellar and binary processes to the basic Monte Carlo code, as well as a simple but accurate implementation of a static tidal boundary in the Galactic field (Joshi et al. 2000b). As a first application, we have studied the dependence on initial conditions of globular cluster lifetimes in the Galactic environment. As in previous Fokker-Planck studies (Chernoff & Weinberg 1990; Takahashi & Portegies Zwart 1998), we include the effects of a power-law initial mass function (IMF), mass loss through a tidal boundary, and single star evolution, and we consider initial King models with varying central concentrations. We find that the disruption and core-collapse times of our models are significantly longer than those obtained with previous 1D (isotropic) Fokker-Planck calculations, but agree well with more recent results from direct $`N`$-body simulations and 2D Fokker-Planck integrations (see also the article by Takahashi in this volume). In agreement with previous studies, our results show that the direct mass loss due to stellar evolution causes most clusters with a low initial central concentration to disrupt quickly in the Galactic tidal field. The disruption is particularly rapid for clusters with a relatively flat IMF. Only clusters born with high central concentrations or with very steep IMFs are likely to survive to the present and undergo core collapse.

In another recent study, we have used our Monte Carlo code to examine the development of the Spitzer “mass stratification instability” in simple two-component clusters (Watters, Joshi, & Rasio 2000). We have performed a large number of dynamical simulations for star clusters containing two stellar populations with individual masses $`m_1`$ and $`m_2>m_1`$, and total masses $`M_1`$ and $`M_2<M_1`$. We use both King and Plummer model initial conditions and we perform simulations for a wide range of individual and total mass ratios, $`m_2/m_1`$ and $`M_2/M_1`$, in order to determine the precise location of the stability boundary in this 2D parameter space. As predicted originally by Spitzer (1969) using simple analytic arguments, we find that unstable systems never reach energy equipartition, and are driven to rapid core collapse by the heavier component. These results have important implications for the dynamical evolution of any population of primordial black holes or neutron stars in globular clusters. In particular, primordial black holes with $`m_2/m_110`$ are expected to undergo very rapid core collapse independent of the background cluster, and to be ejected from the cluster through dynamical interactions between single and binary black holes (see Portegies Zwart & McMillan 2000 and references therein). We have also used Monte Carlo simulations of simple two-component systems to study the evaporation (or retention) of low-mass objects in globular clusters, motivated by the surprising recent observations of planets and brown dwarfs in several clusters (Fregeau et al. 2000).

Much of our current work concerns the treatment of dynamical interactions with primordial binaries. We are in the process of completing a first study of globular cluster evolution with primordial binaries (Joshi & Rasio 2000), based on the same set of approximate cross sections and recipes for dynamical interactions used in the Fokker-Planck simulations of Gao et al. (1991). Typical results are illustrated in Fig. 2. The heating of the cluster core generated by a small population of primordial binaries can support the cluster against core collapse for very long times.

The addition of binary stellar evolution processes will allow us to study in detail the dynamical formation mechanisms for many exotic objects, such as X-ray binaries, millisecond radio pulsars, and cataclysmic variables, which have been detected in large numbers in globular clusters. For example, exchange interactions between neutron stars and primordial binaries can lead to common-envelope systems and the formation of short-period neutron-star / white-dwarf binaries that can become visible both as ultracompact X-ray binaries and binary millisecond pulsars with low-mass companions (see, e.g., Camilo et al. 2000, on observations of 20 such millisecond radio pulsars in 47 Tuc; Rasio, Pfahl, & Rappaport 2000 present a preliminary study of this dynamical formation scenario, based on simplified Monte Carlo simulations).

We are also currently working on incorporating a more realistic treatment of tidal interactions, and, in particular, tidal shocking through the Galactic disk (based on Gnedin, Lee, & Ostriker 1999). Tidal shocks can accelerate significantly both core collapse and the evaporation of globular clusters, reducing their lifetimes in the Galaxy (Gnedin & Ostriker 1997).

Future work will include a fully dynamical treatment of all strong binary-single and binary-binary interactions (exploiting the parallelism of the code to perform separate numerical 3- or 4-body integrations for all dynamical interactions) as well as a fully dynamical treatment of tidal shocking (performing short, $`N`$-body integrations for each passage of the cluster through the Galactic disk or bulge).

#### Acknowledgments.

This work is supported in part by NSF Grant AST-9618116, NASA ATP Grant NAG5-8460, and by an Alfred P. Sloan Research Fellowship. Our computations are performed on the Cray/SGI Origin2000 supercomputer at Boston University under NCSA Grant AST970022N.

## References

Aarseth, S.J. 1999, PASP, 111, 1333

Camilo, F., Lorimer, D.R., Freire, P., Lyne, A.G., & Manchester, R.N. 2000, ApJ, in press \[astro-ph/9911234\]

Chernoff, D.F. & Weinberg, M.D. 1990, ApJ, 351, 121

Fregeau, J.M., Joshi, K.J., Portegies Zwart, S., & Rasio, F.A. 2000, AAS Meeting 196, #41.11

Gao, B., Goodman, J., Cohn, H., & Murphy, B. 1991, ApJ, 370, 567

Gnedin, O.Y., & Ostriker, J.P. 1997, ApJ, 474, 223

Gnedin, O.Y., Lee, H.M., & Ostriker, J.P. 1999, ApJ, 522, 935

Hénon, M. 1971a, Ap. Space Sci., 13, 284

Hénon, M. 1971b, Ap. Space Sci., 14, 151

Joshi, K.J., Rasio, F.A., & Portegies Zwart, S. 2000a, ApJ, in press (Aug 20) \[astro-ph/9909115\]

Joshi, K.J., Nave, C., & Rasio, F.A. 2000b, ApJ, submitted \[astro-ph/9912155\]

Joshi, K.J., & Rasio, F.A. 2000, in preparation

Portegies Zwart, S.F., & McMillan, S.L.W. 2000, ApJ, 528, L17

Rasio, F.A., Pfahl, E.D., & Rappaport, S.A. 2000, ApJ, 532, L47

Spitzer, L., Jr. 1969, ApJ, 158, L139

Spitzer, L., Jr. 1987, Dynamical Evolution of Globular Clusters (Princeton: Princeton University Press)

Spitzer, L., Jr. & Hart, M.H. 1971a, ApJ, 164, 399

Spitzer, L., Jr. & Hart, M.H. 1971b, ApJ, 166, 483

Takahashi, K., & Portegies Zwart, S.F. 1998, ApJ, 503, L49

Watters, W.A., Joshi, K.J., & Rasio, F.A. 2000, ApJ, in press (Aug 10) \[astro-ph/9912457\]

### Ivan King:

You showed that binaries stabilize a cluster for much longer than a Hubble time. Yet we see 15–20% of clusters in a state of core collapse. How do you reconcile this?

### Fred Rasio:

Indeed, if the type of evolution shown in Fig. 2 applied to all globular clusters, there would be no “core-collapsed” clusters in the Galaxy. However, the timescale on which real clusters will exhaust their primordial binary supply and undergo (deep) core collapse depends on a number of factors not considered here: the initial primordial binary fraction (some clusters may have much fewer binaries than the 10% assumed in Fig. 2), the orbit of the cluster in the Galaxy (the simulation of Fig. 2 is for an isolated cluster, but mass loss and tidal shocking can accelerate the evolution dramatically), the stellar IMF (the cluster shown in Fig. 2 contains all equal-mass stars and binaries; a more realistic mass spectrum will also accelerate the evolution), etc.

However, the simple picture that emerges from Fig. 2 may well, to first approximation, describe the dynamical state of most Galactic globular clusters observed today. Note that, for a cluster in the stable “binary burning” phase of Fig. 2, the ratio of half-mass radius to core radius $`r_\mathrm{h}/r_\mathrm{c}210`$ (for $`t10`$Gyr and $`t_{\mathrm{rh}}0.110`$Gyr), which is precisely the range of values observed for the $`80\%`$ of clusters that have a well-resolved core and are well-fitted by King models. Some of these clusters may have gone in the past through a brief episode of “moderate core collapse” (as shown around $`t60t_{\mathrm{rh}}`$ in Fig. 2). Yet, I do not believe that they should be called “core-collapsed” or “post-core-collapse” (nor would they be classified as such by observers). Unfortunately some theorists will even call “core-collapsed” clusters that have just reached the initial phase of binary burning ($`t1050t_{\mathrm{rh}}`$ in Fig. 2). Since the core has just barely contracted by a factor $`23`$ by the time it reaches this phase, it seems hardly justified to speak of a “collapsed” state.