text:
# Variability of the X-ray P Cygni Line Profiles from Circinus X-1 Near Zero Phase

## 1 Introduction

The discovery of broad P Cygni line profiles from the luminous X-ray binary Circinus X-1 (hereafter Cir X-1), which we reported in Brandt & Schulz (2000; hereafter Paper I), was the first time such lines had been observed clearly in the X-ray spectrum of a cosmic object. The P Cygni lines demonstrated the presence of a high-velocity outflow, in line with other suggestions of outflow in this system (e.g., Johnston, Fender, & Wu 1999). However, the nature of Cir X-1 in general is still poorly understood and, despite advances in recent years, there remains uncertainty about even the most basic properties of this system. Since its discovery (Margon et al. 1971), it has appeared bright and variable in X-rays exhibiting a period of 16.6 days (Kaluzienski et al. 1976), and it has thus been a frequent target of most X-ray observatories. The compact object in the Cir X-1 system is thought to be a neutron star (Tennant, Fabian, & Shafer 1986) that can radiate at super-Eddington luminosities during times of strong mass transfer. Its heavily reddened optical counterpart (e.g., Moneti 1992) shows strong, asymmetric H$`\alpha `$ emission; this emission is variable and appears to arise from several sites in the system including the accretion disk (e.g., Whelan et al. 1977; Mignani, Caraveo, & Bignami 1997; Johnston et al. 2001). The system shows two arcminute-scale radio jets (Stewart et al. 1993), and an arcsecond-scale asymmetric jet (Fender et al. 1998) suggests the presence of relativistic outflow from the source. Cir X-1 is often included among the “Galactic microquasar” X-ray binaries (Mirabel 2001).

The nature of the companion star is an important issue for the interpretation of the observed X-ray P Cygni lines. Specifically, one must address if the P Cygni lines could be made by absorption of X-rays from the neutron star in a high-velocity wind from the companion; such a scenario is likely for the high-mass X-ray binary Cygnus X-3, which also shows X-ray P Cygni lines (Liedahl et al. 2000) and has a Wolf-Rayet star companion. The large velocities observed for the X-ray P Cygni lines of Cir X-1 (up to $`\pm 1900`$ km s-1) would require a high-mass companion, probably of spectral type O.

Several attempts have been made to determine the spectral type of the companion star. While a low-mass companion is favored by most of the recent studies (e.g., Stewart et al. 1991; Glass 1994), the constraints remain weak due largely to the heavy interstellar reddening. Additional evidence for a low-mass companion comes from observations of the correlated X-ray spectral and timing properties, observations of type I X-ray bursts, and theoretical considerations. For example, Shirey, Bradt, & Levine (1999a) and Shirey, Levine, & Bradt (1999b) have demonstrated in extensive studies of RXTE data that Cir X-1 exhibits spectral branches in the hardness-intensity diagram that can be identified with the horizontal, normal, and flaring branches of “Z” type low-mass X-ray binaries (LMXBs; e.g., Hasinger & van der Klis 1989; Schulz, Hasinger, & Trümper 1989). Qu, Yu, & Li (2001) have also noted some similarities to the LMXBs GX $`51`$ and Cyg X-2. Furthermore, Tauris et al. (1999) considered the (somewhat uncertain) kinematic properties of the Cir X-1 system and argued that the companion is a low-mass ($`2.0`$ $`M_{}`$), unevolved star.

Given that the stellar wind from the likely low-mass companion cannot explain the observed properties of the X-ray P Cygni profiles from Cir X-1, we interpreted these profiles in the context of an equatorial wind driven from the accretion disk by a combination of Compton heating and radiation pressure (see Paper I). The fact that Cir X-1 can radiate with high luminosity relative to its Eddington luminosity ($`L/L_{\mathrm{Edd}}`$) makes it a natural system in which to expect observable outflows, because of the larger amount of photon pressure available per unit gravitational mass. Accretion-disk winds in X-ray binaries have been discussed both theoretically and observationally (e.g., Begelman, McKee, & Shields 1983; Raymond 1993; Chiang 2001; Proga & Kallman 2001), and the P Cygni lines seen from Cir X-1 match the lines predicted by Raymond (1993) reasonably well. We also note that Iaria et al. (2001a, 2001b) have recently found evidence for a large column density of ionized gas along the line of sight based on observations of ionized iron K edges; the same gas might well produce both the X-ray P Cygni lines and the ionized iron K edges.

P Cygni lines are a common property of the ultraviolet spectra of cataclysmic variables (CVs) possessing an accretion disk. In these systems, the line profiles are associated with high-velocity polar outflows (e.g., Córdova & Howarth 1987), in contrast to the equatorial geometry proposed for the accretion-disk wind in Cir X-1. The proposed equatorial geometry was motivated by evidence that the Cir X-1 system is viewed in a relatively edge-on manner. Although eclipses are not observed, the spectral variability caused by the observed X-ray absorption is best explained by a model with a relatively edge-on accretion disk (Brandt et al. 1996; Shirey et al. 1999b). In this model, X-rays created near the neutron star reach the observer via two light paths: (1) a direct, but often absorbed, light path that can intersect the outer bulge of the accretion disk, and (2) an indirect, electron-scattered light path that avoids the absorption associated with the direct light path. This model can naturally explain the large observed changes in the column density of the absorbing gas without corresponding changes in its apparent covering fraction (see Brandt et al. 1996 for details). We note that radio imaging of SS433, which shares several similarities with Cir X-1, provides strong evidence for an equatorial wind-like outflow from this system (Blundell et al. 2001).

In this second paper, we present results from a variability study of the X-ray P Cygni lines from Cir X-1. As in Paper I the observations were performed with the High-Energy Transmission Grating Spectrometer (HETGS) onboard Chandra at “zero phase.” Zero phase in Cir X-1 is thought to be associated with the periastron passage of the neutron star, and near this phase Cir X-1 is highly variable in X-ray intensity. We present analyses of two observations, which include a re-assessment of the data presented in Paper I representing a high-flux state of Cir X-1 and a second observation during which Cir X-1 was in a low-flux state with strong variability. In Schulz & Brandt (2001; hereafter SB2001) we already presented a broad-band continuum fit to some of the data as well as showed temporal changes of the Si XIV line. The analysis here focuses on a comparison of the bright lines in the two zero-phase flux states and their variability with intensity and time.

Throughout this paper, we adopt a distance to Cir X-1 of 6 kpc (Stewart et al. 1993; Case & Bhattacharya 1998) and an interstellar column density of $`2\times 10^{22}`$ cm-2 (e.g., Predehl & Schmitt 1995).

## 2 Chandra Observations

Cir X-1 was observed with the Chandra HETGS (C.R. Canizares et al., in preparation) on 2000 February 29 (starting at 22:09:50 UT) for 30 ks (observation I) and on 2000 June 8 (starting at 04:13:01 UT) for another 15 ks (observation II).111For details on Chandra and the HETGS see the Chandra Proposers’ Observatory Guide at http://asc.harvard.edu/udocs/docs/. Both observations were performed near zero phase of the binary orbit. They were not precisely at the same phase; observation I started at phase 0.9917, and observation II started at phase 0.9923 (adopting the ephemeris of Glass 1994). In Paper I we showed that observation I occurred during the onset of a luminous flare. Figure 1 puts observation II similarily into perspective with the long-term RXTE ASM light curve. Here it appears that we observed the source during a prominent intensity dip. Observation I is also different from observation II because it suffered from a misalignment of the applied subarray during the observation which led to the loss of half of the data. Otherwise observation II was performed with a similar instrumental setup; the zeroth order was blocked from being telemetered, and the same subarray with CCDs S0 and S5 inactive was applied leading to a frame time of 1.7 s.

### 2.1 Raw-count spectra

Observation I was the basis for the discovery of the P Cygni profiles, and most of its basic properties were described in Paper I. We reprocessed these data using several improved calibration data products. The first one concerns a small correction to the wavelength scale, which arose from a re-assement of the detector pixel size and which led to a scale factor of 500 ppm in the wavelength scale. This leads to a maximal difference of 0.0075 Å at 15 Å compared to the previous definition. For shorter wavelengths this change is of the order of the previous accuracy of the wavelength determination. A second change came from the application of the CCD quantum efficiency, which has now been provided with a finer energy grid. The benefits here are that instrumental edges, for our interest especially the CCD Si-edge, have been accounted for more accurately and are less likely to produce spectral artifacts. A third improvement came from correction of a slight mismatch of the energy grids in the ancillary data files and the grating response. All data were processed the same way with the exception that in observation II we additionally added positive and negative orders together. In this respect we re-gained essentially the same exposure as in observation I. The raw-count spectra in Figure 2 (top) indicate a very high number of counts during observation I. The MEG 1st order spectra have $`1.4\times 10^6`$ events, the HEG 1st order spectra have $`1.5\times 10^6`$ events, and the MEG 3rd order spectra have $`3.9\times 10^5`$ events.

In observation II (middle) the MEG 1st order spectra have $`3.3\times 10^5`$ events, the HEG 1st order spectra have $`2.1\times 10^5`$ events, and the MEG 3rd order spectra have $`4.0\times 10^4`$ events. The difference in intensity is well illustrated by the raw-count spectra shown in Figure 2. Throughout most of the spectral bandpass, the number of counts in observation II is significantly lower in each bin compared to observation I. This had the favorable effects that we did not suffer from any frame dropouts during data transmission and were much less affected by pile-up. The number of counts per bin in observation II drops below $`1`$ above 12 Å, which basically reflects the noise level. In Figure 2 we therefore terminate the spectra at 12 Å for better visibility of the details at shorter wavelengths, although the MEG spectrum in observation I continues up to 15.5 Å (see Paper I).

### 2.2 Pile-up issues

The amount of pile-up in the spectra makes detailed modeling of the X-ray continuum difficult. Pile-up affects the MEG and HEG spectra over a considerable bandpass where the number of incident counts per frame in a given pixel is very high. In order to get a lower limit of the number of detected counts per frame and pixel a rough estimate of its dependence on wavelength, one can divide the count spectra in Figure 2 by the number of CCD frames obtained during the observation period and the number of pixels included in each bin. Values significantly above 0.01 indicate that the spectra are likely affected by pile-up. In general to obtain more precise information there is no straightforward procedure, but there are several other leads that allow one to identify pile-up affected regions in the spectra. High photon densities are only reached in the 1st orders where the effective area of the instrument is highest. In the MEG 1st order spectrum this happens below $``$ 11 Å, while in the HEG 1st order spectrum it happens below $``$ 8 Å. All higher orders are not affected by pile-up. However, they still suffer from contamination of piled-up photons from the 1st orders. Pile-up in the 1st orders creates a depletion of photons at a specific wavelength. Since it is a CCD effect, multiple photon coincidences appear in higher CCD channels since they are recorded as an event that has the sum of the incident energies. In the HETGS these channels also coincide with the locations of the higher order grating spectra. In this respect three first-order photons, for example, appear as one photon in 3rd order at 1/3 of the original wavelength. This means that the MEG 3rd order below $`2`$ Å is contaminated by pile-up from the MEG 1st order. Since no viable model yet exists to correct for this effect, we will avoid all areas in the spectra that suffer substantially from pile-up. In the MEG 1st order the range below 2 Å is also affected by some contamination from zero-order scattering and is excluded as well. For observation I, for example, this leaves the following wavelength bands for spectral modeling: 2.4–8.5 Å in the MEG 3rd order, 8.2–10.1 Å in the HEG 1st order, and above 11 Å in the MEG 1st order. The wavelength range below 2 Å in the HEG 1st order is mildly affected by pile-up, and here we have to add corrections (see §3.5). The available energy bands are thus 0.8–5.2 keV and 6.2–8 keV. Observation II shows considerably lower count rates, and its spectra are much less affected by pile-up and zero-order scattering; here we can use almost the entire 1.0–8.0 keV energy band.

Figure 2 (bottom) shows the exposure-corrected spectra for observations I and II. The top set of spectra with the higher fluxes are from observation I, and the bottom set are from observation II. The spectra are integrated over the entire exposures. Without pile-up the blue, green and red spectra of both sets should overlap. This qualitatively gives an impression of the degree of pile-up, which is indicated by how much the blue and green spectra deviate from the red spectra.

For the analysis of the lines we mostly rely on the HEG 1st order spectra for statistical reasons. This requires some corrections in areas where pile-up appears strong, which is between 2 and 7 Å. From the difference of the exposure-corrected spectra and the pile-up free 3rd order spectra we can determine an energy dependent flux correction factor, which we later use to correct the equivalent widths of the lines. The model of the HEG 1st order spectrum of observation II in Figure 3 has these factors applied, and the figure shows that we get good agreement. The additional systematic uncertainty in the determination of equivalent widths can be up to 20$`\%`$ is the worst case.

### 2.3 Determining fluxes and luminosities

In general we follow the procedure to determine the luminosity outlined in SB2001 and choose our model spectrum from previous fits to ASCA data (Brandt et al. 1996); the model consisted of two blackbody components with both intrinsic partial-covering absorption and interstellar absorption of $`2.0\times 10^{22}`$ cm-2 using cross-sections from Morrison & McCammon (1983). The two observations can be compared to the low-count-rate state (observation II) and high-count-rate state (observation I) in Brandt et al. (1996). We emphasize, however, that the current HETGS spectra do not qualify for a detailed physical discussion of the continuum, and we use the fits mainly as a parameterization of the continuum emission in order to perform the line analysis. We also draw some basic information from the change of the fits from observation I to observation II.

For observation I we find a soft blackbody temperature of $`kT=0.66\pm 0.17`$ keV \[partial covering $`N_\mathrm{H}=(8.2\pm 0.4)\times 10^{23}`$ cm-2\] and a hard blackbody temperature of $`kT=1.22\pm 0.28`$ keV \[partial covering $`N_\mathrm{H}=(2.3\pm 0.2)\times 10^{22}`$ cm-2\]; the partial covering fraction was 45$`\%`$. We find a total absorbed X-ray flux of $`2.0\times 10^8`$ erg cm-2 s-1 with the soft component contributing 61$`\%`$ between 2–8 keV. This corresponds to a luminosity of $`1.8\times 10^{38}`$ erg s-1. During observation II we could fix the the partial covering columns to the values found in observation I; we find a covering fraction of 89$`\%`$. The temperature for the soft blackbody component was $`kT=0.16\pm 0.09`$ keV, and for the hard blackbody component it was $`kT=1.55\pm 0.21`$ keV. The total absorbed X-ray flux dropped to $`4.8\times 10^9`$ erg cm-2 s-1 with a contribution of the soft component of $`<1`$%. The corresponding luminosity is $`2.8\times 10^{37}`$ erg s-1 . We emphasize again that the model parameters we have empirically derived to parameterize the continuum are not optimal for comparison with previous results (Brandt et al. 1996; Shirey et al. 1999). A full pile-up model will be required for detailed HETGS continuum fitting, and even then the limited high-energy sensitivity will remain an obstacle. As we discuss in §3.5, the correction for absorption during observation II is generally uncertain, and hence the true luminosity of Cir X-1 is difficult to constrain during this observation. It is plausible that the true luminosity during observation II is comparable to that during observation I. Such behavior would be consistent with that observed by ASCA (Brandt et al. 1996) and RXTE (Shirey et al. 1999b) during dips. Evidence for heavy and complex absorption in the Chandra data is found in the depths of the neutral Fe K (see §3.5) and maybe the Si K edges. Both observations show a neutral Si K edge with a depth of $`\tau =0.15`$. The large depth of this edge required us to add an additional edge to our fits. However, although current uncertainties in the improved instrument efficiencies cannot account for such a large effect, it was also observed recently in HETGS spectra of other bright sources (Miller et al. 2002).

## 3 P Cygni Line Analysis

### 3.1 Basic line measurements

For the line analysis we binned all spectra into 0.005 Å bins (which oversamples the HEG resolution by a factor of two). All the analysis in this section has been performed on the raw data. The spectra are integrated over the entire exposure. We use lines detected in the HEG 1st and MEG 3rd orders, except in observation I where we use the MEG 1st order above 10 Å. We detect a number of resonance line features in both observations, of which the brightest clearly show P Cygni profiles (Paper I). Figure 3 shows the spectrum of observation II from 1.7–9 Å with all identified features marked. Table 1 includes all line features detected with a signal-to-noise ratio larger than three for both observations. Most of the lines are from H-like and He-like species, although we do detect some Li-like and B-like species from elements with high atomic numbers. For some of the lines the P Cygni shape is less pronounced, and we find them either mostly in emission or absorption (see below). Notable in observation II is the detection of the forbidden line of He-like Mg XI, which cannot appear as a P Cygni line since it is not a dipole transition.

Table 1 lists each ion in column 1, the transition name in column 2 (“H-like” denotes ions with only one electron left, “He-like” with two, and so on), and the predicted wavelength from the SPEX (Mewe 1994) line list in column 3. In column 4 we list the measured wavelength, which represents the intersection between the emission and the absorption with the continuum. The intersection with the continuum is generally a good estimate of the rest-wavelength position. This is underlined by the fact that on average the profiles appear symmetric. In order to quantify the properties of the P Cygni lines, we normalized each profile by dividing the data by the underlying local continuum (where we simply applied a polynomial fit) and then fitted the absorption and emission components with Gaussians. Within the statistical uncertainties, simple Gaussian functions proved to be sufficient for the fits. The velocities listed in column 5 represent the difference between the measured wavelength and the peak position of the Gaussian fitting the absorption part of the line (note, of course, that the maximum velocity to which absorption is observed is higher than this); lines without clear P Cygni shapes do not have velocity values listed. Most emission peaks appear slightly above the expected wavelength in column 3; in a few cases they match well, and here we have probably slightly underestimated the velocity. We find velocities of 200–1900 km s-1.

We determined the equivalent widths of the lines from the fits. They are listed in column 6 of Table 1. The values are measured from HEG 1st order data unless noted otherwise. The equivalent width of a P Cygni line, $`W_{\mathrm{pc}}`$, is defined as the sum of the equivalent widths of the emission and the absorption (taking both to be positive quantities). For the analysis below we also list the contribution of absorption as a fraction of $`W_{\mathrm{pc}}`$ ($`f_\mathrm{a}`$ in column 7 of Table 1). The Fe-line area below 2 Å is treated separately below.

### 3.2 High-flux versus low-flux states

In this section we analyze the variations of the P Cygni lines between observation I and observation II. The measured wavelengths (as well as the ones for the peak and valley positions) of the same transitions in both observations agree to within the stated uncertainties. The most dramatic change is observed in the equivalent widths of the lines. In Figure 4 we compare the P Cygni profiles for the lines with the highest equivalent widths in both observations, which are the Lyman $`\alpha `$ lines of S XVI (left), Si XIV (middle), and Mg XII (right). Pile-up effects for S XVI and Si XIV cannot account for the strong variability of the lines. The equivalent widths in observation I are generally lower than in observation II throughout the whole bandpass. The equivalent widths of the three Lyman $`\alpha `$ lines in Figure 4 differ by a factor of $`4`$, whereas the other lines typically differ by a factor of $`2`$.

### 3.3 Emission lines

The fluxes of the emission lines are shown in column 9 of Table 1. The lines in observation I appear bright with fluxes from $`1.5\times 10^4`$ ph cm-2 s-1 at 10.63 Å to $`2.2\times 10^3`$ ph cm-2 s-1 at 1.85 Å. The general trend of declining line flux with increasing wavelength is largely due to absorption by the interstellar column density. The brightest lines appear to be the Lyman $`\alpha `$ lines of Fe XXVI, S XVI, and Si XIV, as well as He-like Fe XXV. In observation II the lines’ fluxes dropped significantly to a range from $`1.5\times 10^5`$ ph cm-2 s-1 to $`6.2\times 10^4`$ ph cm-2 s-1. There is a trend for the largest decline in line flux to happen toward longer wavelengths; above 6.5 Å the decline is $``$ 12–20 while below 5 Å the decline is $``$ 3–5.

Of particular interest is the He-like Mg XI triplet, because it is the only one we can clearly resolve and appears free of contaminating features; for comparison, the Si XIII triplet is confused with the Mg XII Lyman $`\gamma `$ line as well as the Si K edge. He-like triplets are sensitive to temperature and density and thus can be used as plasma diagnostics; we use the calculations by Porquet & Dubau (2000). Calculation of the $`G=(i+f)/r`$ ratio ($`i`$ is the flux of the intercombination line, $`f`$ the flux of the forbidden line, and $`r`$ the flux of the recombination line) and $`R=i/f`$ ratio allows estimation of the temperature and density, respectively (Gabriel & Jordan 1969). Our best constraints come from observation II since during this observation we detect all the lines in the triplet as shown in Figure 5. We compute a $`G`$ ratio of $`2.9\pm 1.0`$ and an $`R`$ ratio of $`1.1\pm 0.34`$, implying a temperature of $`(3.8\pm 1.1)\times 10^6`$ K and a density of $`(1.9\pm 0.9)\times 10^{13}`$ cm-3. A $`G`$ ratio of $`3`$ indicates a photoionized plasma dominated by recombination (Bautista & Kallman 2000) and justifies the use of the Porquet & Dubau (2000) calculations, which are specifically designed for a recombining plasma. During observation I we only detect the intercombination line at 9.23 Å and use the $`1\sigma `$ uncertainty in the continuum as an upper limit for the detection of the forbidden and resonance lines. We estimate an $`R`$ ratio of $`>0.5`$, implying a density of $`<4\times 10^{13}`$ cm-3. The $`G`$ ratio is not well constrained during observation I.

Figure 5 shows that the Mg XI triplet probably does not suffer from blending, which would mainly stem from possible lower-shell absorption events. These are unlikely given the high ionization states observed. Fe L emission lines are unlikely to appear as blends, because Kallman et al. (1996) predicted that recombination emission is dominated by K-shell ions and not L-shell ions. We cannot exclude the possibility that resonance scattering affects the line ratios; this would tend to increase the flux of the resonance line and thus lower the $`G`$ ratio. The consequence would be that the temperature is higher or, in the extreme case, that the Porquet & Dubau (2000) calculations are not applicable because the plasma is not predominately recombining. However, the resonance line does not appear very strong in the spectrum, and a correction of $`G`$ due to resonance scattering seems unlikely to be highly significant.

### 3.4 Absorption column densities

The measurements of the absorption-line equivalent widths allow estimation of the ionic column densities in the absorbing medium. For an absorption line on the linear part of the curve of growth, the relation between equivalent width and column density is

$$\frac{W_\lambda }{\lambda }=\frac{\pi e^2}{m_\mathrm{e}c^2}N_\mathrm{j}\lambda f_{\mathrm{ij}}=8.85\times 10^{13}N_\mathrm{j}\lambda f_{\mathrm{ij}}$$
(1)
where $`N_\mathrm{j}`$ is the column density of ion species $`j`$, $`f_{\mathrm{ij}}`$ is the oscillator strength of the transition, and $`W_\lambda `$ is the equivalent width of the absorption line (e.g., Spitzer 1978; see also Lee et al. 2001). In Table 1 we list the column density, $`N_{\mathrm{abs}}`$, for each ion species derived using this formula.

Several effects may lead to systematic underestimates of column densities when the formula above is used. Given the observed P Cygni profiles, for example, it is likely that line emission fills in some of the absorption; it is difficult to correct for the resulting column-density underestimates precisely. Furthermore, line saturation can also lead to an underestimate of the column density. The observed lines could be saturated even though they do not appear “black” because of multiple light paths to the observer (e.g., the electron scattering of some X-rays around the absorbing material; see §1). The observed lines might also have unresolved substructure composed of narrower, saturated lines. In fact, since we observe lines from high-order transitions in the Lyman series, the Lyman $`\alpha `$ lines are probably saturated; the oscillator strengths decrease by a factor of $`30`$ from the 1s-2p (Lyman $`\alpha `$) to the 1s-4p (Lyman $`\delta `$) transitions. We therefore see the column densities in Table 1 more as lower limits, and we note that the column densities derived from the high-order Lyman transitions probably have the least systematic error. For the determination of column densities we thus rely only on the highest-order Lyman transition detected.

In order to estimate equivalent hydrogen column densities, we must correct for the ionization fraction of a given ion. From the analysis of the Mg XI triplet, it is likely that we observe a photoionized plasma. In the following we use the calculations by Kallman & Bautista (2001) for such a plasma which provide a representative ionization structure. For a temperature of $`3.8\times 10^6`$ K (see §3.3), we find an ionization parameter of $`\xi =L/nR^220,000`$ erg cm s-1, which is consistent with the fact that we observe mostly H-like ions. We find a similar value of $`\xi `$ when we compare our abundance ratios \[AR=$`N_{\mathrm{abs}}`$(H-like)/$`N_{\mathrm{abs}}`$(He-like)\] with their calculations. Table 2 shows the ionization fractions for Mg, Si, S, and Fe and the derived equivalent hydrogen column densities as determined from the observed abundance ratios; the column densities are given for both a solar-abundance distribution from Morrison & McCammon (1983) and an interstellar-abundance distribution from Wilms, Allen, & McCray (2000).

The equivalent hydrogen column densities in Table 2 appear to increase significantly between observation I and observation II. However, as discussed in §3.3, we only have limited constraints on the ionization structure of the gas during observation I because we do not detect all the lines in the Mg XI triplet. We cannot rule out, for example, a reduction in the ionization fraction during observation II that explains all the observed changes. Some of the derived column densities should be taken with caution. It is unusual, for example, that we observe a strong Mg XII 1s-3p transition, whereas the 1s-2p transition is missing. In fact if we discard the result for this line, the derived column density using Mg would be similar to those computed using the other elements.

TABLE 2
IONIZATION FRACTIONS AND EQUIVALENT HYDROGEN COLUMN DENSITIES        $`N_\mathrm{H}^\mathrm{c}`$   $`N_\mathrm{H}^\mathrm{c}`$  Ion ARa $`f_{\mathrm{ion}}^\mathrm{b}`$ ($`10^{22}`$ cm-2) ARa $`f_{\mathrm{ion}}^\mathrm{b}`$ ($`10^{22}`$ cm-2)   Observation I Observation II  Fe XXVI 0.68 0.327 1.5 (1.9) 1.7 0.491 4.6 (5.6)  Fe XXV $`\mathrm{}`$ 0.514 $`\mathrm{}`$ $`\mathrm{}`$ 0.277 $`\mathrm{}`$  S XVI 2.4 0.480 0.7 (1.1) 8.7 0.290 6.4 (9.9)  S XV $`\mathrm{}`$ 0.199 $`\mathrm{}`$ $`\mathrm{}`$ 0.032 $`\mathrm{}`$  Si XIV 5.1d 0.310 0.7 (1.4) 8.2 0.244 7.9 (15.8)  Si XIII $`\mathrm{}`$ 0.035 $`\mathrm{}`$ $`\mathrm{}`$ 0.019 $`\mathrm{}`$  Mg XII 14.5d 0.237 3.0 (4.7) 22.4d 0.162 7.2 (11.5)  Mg XI $`\mathrm{}`$ 0.016 $`\mathrm{}`$ $`\mathrm{}`$ 0.007 $`\mathrm{}`$     a$`N_{\mathrm{abs}}`$(H-like)/$`N_{\mathrm{abs}}`$(He-like)
bIonization fractions from Kallman & Bautista (2001).
cEquivalent hydrogen column density for a solar-abundance distribution (Morrison & McCammon 1983) and an interstellar-abundance distribution (Wilms, Allen, & McCray 2000; in parentheses).
dEstimated lower limits.

### 3.5 The Fe K line region

The Fe K line region is especially interesting for spectroscopy because Fe has quite high fluorescent yields (Auger destruction is not as dominant as in lower $`Z`$ elements) as well as ions that can survive at high ionization parameters. Figure 6 (top) shows the Fe K region from 1.7–2.0 Å as observed with the HEG during observation I. The model (red curve) includes the continuum level as fitted over the entire available bandpass (see §2.3). In addition we applied an empirical correction to account for pile-up (see §2.2) which we computed by matching featureless continuum stretches in the 1.97–3.00 Å range. Below 1.97 Å additional Gaussians were fitted locally and then added to the model. At the short wavelength side of the Fe edge we additionally adjusted the model to match the continuum level. We observe an Fe K line centered at 1.94 Å, which is consistent with the charge states Fe I–X (e.g., House 1969). The line has an equivalent width of 1.9 mÅ. The spectrum also shows P Cygni lines of Fe XXV and Fe XXVI with the properties given in Table 1. We resolve a neutral Fe edge at 1.75 Å with an optical depth of $`\tau =0.035\pm 0.017`$, implying an equivalent hydrogen column density of $`3.4\times 10^{22}`$ cm-2. After subtraction of the interstellar column density, the column densities of neutral and ionized gas appear similar.

Figure 6 (bottom) shows the same spectral region for observation II. The model was calculated in a similar way to that for observation I; however, no correction for pile-up was needed. We detect an Fe K line with a similar position and equivalent width as that observed during observation I. However, we do not see obvious P Cygni lines from Fe XXV and Fe XXVI but mostly resonance absorption with significantly increased equivalent width (see Table 1). The resonance-line positions are the same as observed during observation I indicating little change in velocity. The fact that we do not see much emission from ionized Fe in observation II may at first seem surprising, since we see P Cygni lines in lower $`Z`$ ion species. However, if one compares the fluxes of the Lyman $`\alpha `$ lines and He-like lines of S, Si, and Mg in Table 1, then the decrease in flux is consistent with that observed in these other ions.

A substantial increase in the depth of the neutral Fe K edge is observed between observation I and observation II. During observation II the depth is $`\tau =0.108\pm 0.028`$. To derive this value we had to take into account an unidentified feature at 1.73 Å; because of this feature we had to force the edge to match the high-energy continuum. It is difficult to derive a column density from the edge depth. Due to the partial-covering nature of the absorption and the fact that the neutral absorbing column density is expected to be $`>10^{24}`$ cm-2 during observation II (see §4.4 of Shirey et al. 1999b), the true column density is likely to be substantially larger than expected from simple consideration of the apparent edge depth. A column density of $`>10^{24}`$ cm-2 will absorb much of the direct continuum at wavelengths longward of the edge, and due to this “saturation” effect the true depth of the edge (and hence the column density) can be severely underestimated.

## 4 Short-term Variability of the Lines

We have also studied short-term variability of the P Cygni lines on time scales within the 30 ks and 15 ks exposure times. To do this, we apply a “sliding-window” technique, where we define a fixed exposure-time window that is sufficiently long to receive a significant signal from the bright lines. At a set of equally spaced times during the observation (separated by a constant “time step”), we then create spectra using the data within the exposure-time window (the exposure-time window is centered on each of the equally spaced times). In the case of observation I, we selected an exposure-time window length of 5 ks and a time step of 1 ks (note that the resulting spectra have some overlap). Within a total exposure time of 30 ks, this allows us to create 24 spectra. Using these spectra, we then created three-dimensional “spectral maps” where the $`x`$-axis is the wavelength of each spectrum in 0.015 Å steps, the $`y`$-axis is the time in 1.0 ks steps, and the $`z`$-axis is the photon flux in each pixel.

Figure 7 shows spectral maps for three line regions: the Fe line region (left), the Si XIV Lyman $`\alpha `$ line region (middle), and the Mg XII Lyman $`\alpha `$ line region (right). The colors represent the flux levels from low (green to black) to high (orange to yellow) on a logarithmic scale. The scale ranges from $``$ 0.1–0.8 ph cm-2 s-1. The color changes in these maps in the vertical direction show the variability of the spectra. This variability appears somewhat smoothed since the windows overlap.

The spectral continuum is clearly changing; this is indicated by the large-scale color distribution in each map. At short wavelengths (left map), the continuum stays high throughout most of the observation and then drops rapidly in the last 5 ks. At intermediate wavelengths (middle map), the continuum remains low throughout the first one-third of the observation, increases rapidly around the middle of the observation, and remains high throughout most of the rest of the observation. At longer wavelengths (right map), this trend is continued; the continuum stays low throughout the first half of the observation. This qualitative assessment of the continuum variability is consistent with the RXTE results of Shirey et al. (1999a) who find the spectra moving along spectral branches in the X-ray color-color diagram.

The three maps in Figure 7 show that the shapes of the P Cygni lines are highly variable on time scales of 2–3 ks. In the time maps at the top of Figure 7 we observe a correlation of the emission and an anti-correlation of the resonance absorption with the local continuum flux. When the local continuum is low, the emission-line parts of the P Cygni profiles are weak; they grow stronger with increasing continuum flux. The absorption-line parts show the opposite behavior; their amplitudes are larger at lower continuum fluxes, and absorption practically vanishes at high flux levels. This behavior on short time scales reflects the same general pattern we observed by comparing the high-flux state in observation I with the low-flux state in observation II. This also explains why for Fe XXV the emission line reaches its peak flux about 7 ks before we observe the peak in the Si XIV emission line. The line flux is directly connected with the local continuum flux, and in this case the spectral change pattern creates this time difference.

Figure 7 also shows the spectra used to make the spectral maps. The time sequence of the Si XIV line at 6.18 Å shows that the absorption part of the profile has vanished by about 12 ks into observation I, which would imply a column density change of Si XIV from about twice the amount in Table 1 to nearly zero on a time scale of $`3`$ hr. The fact that the flux of the local continuum is anti-correlated with the absorption amplitude suggests that changes in the ionization fraction, rather than changes in column density, are primarily responsible for the changes in absorption strength (although we cannot strictly rule out the latter).

The variability of the Si XIV line during observation II is shown in Figure 8. Although during observation II the overall X-ray flux is substantially lower than during observation I, the ionization parameter appears to be varying over approximately the same range. The emission-line fluxes during observation II, on the other hand, are lower by significant factors (see §3.3).

Figure 8 also illustrates the correlation of spectral hardness and absorption strength for observation II. We computed a hardness ratio by dividing the flux in the 1.7–3 Å band by that in the 3–8 Å band; this ratio is similar to the soft-band hardness ratio computed by Shirey et al. (1999b). The flux is computed in the 1.7–8 Å band. We used the same sliding-window technique as used for the spectral maps, and thus in Figure 8 the sequence of data points in the hardness-ratio diagram can be directly matched to the sequence of line profiles. The data points in Figure 8 are connected by dotted lines indicating the time sequence. In the inset of Figure 8 we also plot the softest (spectrum 1) and hardest (spectrum 2) spectra during observation II, which occurred near the beginning and the end of the observation, respectively. Figure 8 shows that while the flux is generally decreasing throughout the observation, the hardness ratio increases until it reaches a plateau at $`2.3`$. This spectral behavior is consistent with that studied during dip periods with RXTE (Shirey et al. 1999b). The more the continuum hardens, the more the absorption in the P Cygni profiles decreases.

During observation I we see similar variations of the P Cygni line on similar time scales. However, here the spectrum is dominated by the soft component, and spectral variations of the hard component cannot be directly measured.

## 5 Discussion

### 5.1 Refined constraints on the accretion-disk wind model

In Paper I we argued that the P Cygni profiles are most likely produced by a high-velocity outflow from an accretion disk viewed in a relatively edge-on manner. The central X-ray source illuminates the disk and produces a wind driven by both thermal and radiation pressure (e.g., Begelman et al. 1983). The intermediate-temperature ($`5\times 10^6`$ K) part of this wind produces the H-like and He-like X-ray lines (e.g., Raymond 1993; Ko & Kallman 1994) with P Cygni profiles. Our current data remain generally consistent with this scenario, although we recognize that the strong wind variability seen from Cir X-1 is not usually accounted for in the predominantly time-independent theoretical studies; we use these studies as general guides to interpretation rather than for strict comparison purposes.

Our new analyses, combined with other recent work on accretion-disk winds and photoionized plasmas, have allowed us to extend some of the findings in Paper I. For example, one of the potential problems with the disk-wind model noted in §3 of Paper I can now be ameliorated. There we noted that likely values for the wind’s launching radius ($`r_{\mathrm{launch}}10^5`$ km), launching density ($`n_{\mathrm{launch}}>10^{15}`$ cm-3), and ionization parameter ($`\xi 1000`$ erg cm s-1) led to a wind that was optically thick to electron scattering; most line photons attempting to traverse the wind would then be Compton-scattered out of the line. However, our analysis in §3.4, which uses the calculations of Kallman & Bautista (2001) rather than those of Kallman & McCray (1982), now suggests that a significantly higher ($`\xi 20,000`$ erg cm s-1) ionization parameter is likely for much of the gas creating the P Cygni lines. The required column density along the line of sight is thus correspondingly reduced. Hydrogen column densities of $`10^{23}`$ cm-2 can now be accomodated with only moderate shielding of the wind from the full X-ray continuum or clumping of the wind.

The density measurement in the line-emission region of $`2\times 10^{13}`$ cm-3 using the Mg XI triplet (see §3.3) provides a further independent constraint on the disk-wind model. Following the arguments in §3 of Paper I with the new ionization parameter of $`\xi 20,000`$ erg cm s-1, the expected launching density is $`n_{\mathrm{launch}}>5\times 10^{13}`$ cm-3. The measured density in the line-emission region is somewhat, but not greatly, lower than the launching density. Note that since the launching density depends upon the assumed launching radius, the general consistency found above supports our adopted $`r_{\mathrm{launch}}10^5`$ km (chosen to make the observed terminal velocity of the wind comparable to the escape velocity at $`r_{\mathrm{launch}}`$; see §3 of Paper I).

Gratings observations of other neutron star LMXBs thus far have not generally revealed P Cygni lines as prominent as those observed from Cir X-1 (e.g., Cottam et al. 2001; Marshall et al. 2001; Paerels et al. 2001; Schulz et al. 2001; Sidoli et al. 2001). It is difficult to account for this difference quantitatively, but it may be due to a combination of mass accretion rate (relative to the Eddington rate) and inclination. The high, and perhaps super-Eddington, mass accretion rate of Cir X-1 allows an unusually powerful outflow to be driven from the accretion disk, and a relatively edge-on geometry is optimal for the viewing of an equatorial outflow.

### 5.2 Line and continuum variability

Given the results from ASCA and RXTE (Brandt et al. 1996; Shirey et al. 1999b), the large difference in overall X-ray flux between observation I and observation II is likely to be due mostly to a change in the amount of absorbing material along the line of sight. As discussed by several authors (see §1), this material is probably associated with a thickening of the accretion disk induced by the strong mass transfer occurring at zero phase. The change in Fe K edge depth observed by Chandra is consistent with a large change in the amount of absorbing material, although it is difficult to constrain the exact column density during observation II. The fact that we observe the same basic P Cygni lines during observation I and observation II supports the idea that the underlying ionizing continuum source has not changed dramatically in strength. Furthermore, the fact that the line equivalent widths are significantly larger during observation II than observation I (see Figure 4) is generally consistent with the idea that blockage of direct continuum emission is occurring during observation II. However, it is worth noting that the emission-line fluxes do drop significantly between observation I and observation II (see §3.3). This can be understood without invoking continuum changes if the absorbing material screens ionizing photons from reaching a significant fraction of the line-emitting gas (this screening could be stronger at longer wavelengths, explaining the trend with wavelength noted in §3.3).

However, while it seems clear that changes in absorption play a critical role in causing much of the variability of Cir X-1, some changes in underlying continuum strength and shape are probable as well. These are particularly likely during observation I since it was made while Cir X-1 was undergoing an X-ray flare above its “base” level (see Figure 1 of Paper I). Our data are not ideal for probing the nature of the continuum variability in detail, but Shirey et al. (1999a) have demonstrated that Cir X-1 shows spectral variations that classify it as a “Z” type LMXB. The Chandra data strongly suggest that much of the short-term P Cygni line variability is caused by spectral changes in the local continuum which lead to changes in the ionization level of the wind (see §4), although it is difficult to rule out entirely changes in the wind geometry as well. The observed short-term changes of the lines can be dramatic, causing them to vary from being almost completely in absorption to almost completely in emission. Modeling this variability in detail will be a challenge.

### 5.3 Comparisons with the spectral features seen from superluminal jet sources

The Fe K line region in the low-flux state spectrum of Cir X-1 shows strong, blueshifted resonance absorption features from Fe XXV and Fe XXVI. Significant line emission during the low-flux state is only seen from neutral or nearly neutral Fe. Such spectral behavior is reminiscent of that seen from Galactic superluminal jet sources containing a black hole like GRS 1915+105 (Kotani et al. 2000; Lee et al. 2001) and GRO J1655–40 (Ueda et al. 1998; Yamaoka et al. 2001). Similar Fe absorption has also now been observed in GX 13+1 (Ueda et al. 2001), a neutron star LMXB that has been classified as an atypical atoll source (Hasinger & van der Klis 1989), and MXB 1659–298 (Sidoli et al. 2001), also probably an atoll source; these results show that these spectral features can appear independent of the nature of the compact object and the existence of jets. In most of these cases it has been argued that the absorption lines come from a highly ionized plasma layer extending to large heights above the equatorial plane (Yamaoka et al. 2001). Specifically, a geometrically thick accretion flow at $`10^4`$ Schwarzschild radii is assumed, which in the case of GRO J1655–40 amounts to $`10^5`$ km. This is of the same magnitude as our launching radius for the wind in Cir X-1, and it is worth considering the possibility that the material creating the Fe X-ray absorption lines in superluminal sources might generally be outflowing rather than inflowing. An equatorial wind scenario appears consistent with the gratings observations of GRS 1915+105 (Lee et al. 2001); if the outflow velocity for GRS 1915+105 is assumed to be the same as that for Cir X-1, the outflow would need to be less inclined relative to the line of sight to be consistent with the 770 km s-1 upper limit on the blueshift. The strong Fe XXVI absorption feature from GRS 1915+105 shows variations on similar time scales to those observed from Cir X-1.

Superluminal jet sources containing black holes have not shown as many spectral features as we observe from Cir X-1; typically only features from highly ionized iron are seen. This could plausibly be explained if their winds are more highly ionized due to their larger luminosities and different continuum shapes, so that iron is the only abundant element not fully stripped of its electrons.

We thank all the members of the Chandra  team for their enormous efforts. We gratefully acknowledge the financial support of CXC grant GO0-1041X (WNB, NSS) and Smithsonian Astrophysical Observatory contract SV1-61010 for the CXC (NSS).


format: text

text:
# 1 Introduction

## 1 Introduction

The problem of intelligent information processing in the ever-expanding universe has recently been investigated by several authors (Krauss and Starkman 1999; Ćirković and Bostrom 2000). These considerations have entered a new phase, building upon advances in observational cosmology and the strong foundations of the pioneering studies in cosmological prediction by Rees (1969), Dyson (1979), Frautschi (1982), Tipler (1986) and Barrow and Tipler (1986). In the same time, interest in the fundamental physical limits of computation has grown, largely motivated by tremendous advances in computer science (e.g. Lloyd 2000). Continuation of such trend of technological progress will lead humanity, sooner or later, into a stage of highly advanced galactic civilization, in recent years conventionally dubbed the posthuman era. Apart from the epochal importance of any such possible development for biological and social sciences, it has a strong bearing on the nascent astrophysical discipline of physical eschatology, for the reasons suggested boldly by Dyson (1979), when he wrote:

> It is impossible to calculate in detail the long-range future of the universe without including the effects of life and intelligence. It is impossible to calculate the capabilities of life and intelligence without touching, at least peripherally, philosophical questions.

Certainly, the most important property of intelligence is its capacity for information processing. As it has been recognized for a long time, this form of information processing is still entirely within limits of physical, specifically thermodynamical, laws. This conclusion does not necessarily entail reductionism, since consciousness may still contain uncomputable and therefore irreducible elements in addition to the conventionally established ones (e.g. Penrose 1989). And thermodynamics in general is, in the limit of very long timescales, determined by properties of the universe as a whole, i.e. by astrophysics and cosmology. Our goal in this note is to consider very simple cases of maximal informational resources available to advanced civilizations, either extraterrestrial or posthuman. We do not enter into any possible distinction between the two, since the ages of humanity, Earth and our Galaxy, as well as the pace of chemical evolution, suggest that there may be intelligent communities much older than ours even in our cosmological neighbourhood (Livio 1999).

## 2 An estimate of the maximal quantity of processed information

Let us consider various energy fields as potential sources of energy for information processing, using the assumption (”Cosmic Sum rule”) that $`\mathrm{\Omega }=\mathrm{\Omega }_\mathrm{\Lambda }+\mathrm{\Omega }_b+\mathrm{\Omega }_{\mathrm{CDM}}=1`$, where the first term corresponds to the vacuum energy, second to the baryonic matter, and the third to the (non-baryonic) cold dark matter (for the general cosmological considerations, see Peebles 1993). The best present estimates approach values $`\mathrm{\Omega }_\mathrm{\Lambda }0.7`$, $`\mathrm{\Omega }_{\mathrm{CDM}}0.25`$ and $`\mathrm{\Omega }_b0.05`$ (e.g. Ćirković and Bostrom 2000, and references therein). One cannot do anything useful with the vacuum energy.111On the contrary, it plausibly harms the energy budget by preventing, through the so-called cosmological ”no-hair” theorem, the development and subsequent exploitation of the cosmological shear during the Hubble expansion (Gibbons and Hawking 1977). As far as CDM is concerned, it could conceivably be used as an energy source, since the annihilation of these cosmions and anticosmions (present in approximately equal numbers according to the standard theory) would produce potentially usable energy (for some consequences of annihilation, see Kaplinghat, Knox and Turner 2000). However, depending on the mass spectrum of cosmions, their galactic density is rather small, and since their interactions are by definition very weak, their gathering and separation for annihilation will pose huge engineering problems. If we consider the model of ”posthumanity soon”, than only usable matter field is the baryonic matter, in the local environment concentrated in the form of planetary systems.

We may use the Brillouen’s (1962) equality for the upper limit of processable quantity of information (in bits):

$$I_{\mathrm{max}}=\frac{\mathrm{\Delta }E}{k_BT\mathrm{ln}2}=1.05\times 10^{16}\frac{\mathrm{\Delta }E}{T}.$$
(1)
Here, $`k_B`$ is the Boltzmann constant. Both $`\mathrm{\Delta }E`$ and $`T`$ are functions of time, as well as of cosmological parameters $`\mathrm{\Omega }`$, $`\mathrm{\Lambda }`$ and $`H_0`$. There are several important consequences of the Brillouen’s inequality for cosmology. The amount of information available for processing in any causally connected region of finite proper size with trivial topology is necessarily finite.

As for the working temperature of posthuman computers, we may wish to consider two somewhat extreme models: (I) the first with $`T=T_{\mathrm{CMB}}=2.730\pm 0.014`$ K (Staggs et al. 1996)—indicating quick reaching of posthuman stage; additional cooling would require diverting the precious energy resources from computing itself—and (II) the second with

$$T=T_{\mathrm{vac}}=\frac{\mathrm{}c}{k_B}\sqrt{\frac{\mathrm{\Lambda }}{12\pi ^2}}=3.2985\times 10^{30}h\left(\frac{\mathrm{\Omega }_\mathrm{\Lambda }}{0.7}\right)^{\frac{1}{2}}\mathrm{K}$$
(2)
indicating posthumanity in the asymptotic limit of physical eschatology. In this equation, $`\mathrm{\Lambda }`$ is the vacuum energy-density corresponding to the dimensionless cosmological density fraction $`\mathrm{\Omega }_\mathrm{\Lambda }`$, and $`h`$ is the dimensionless Hubble constant ($`H_0100h`$ km s-1 Mpc-1).

Now we can write for the expendable energy:

$$\mathrm{\Delta }E=\underset{V}{}q(\mathrm{\Omega }_b\rho _{\mathrm{crit}})c^2𝑑V,$$
(3)
where $`V`$ is the proper volume available to posthuman civilization and $`q`$ is the energy extraction efficiency (between 0 and 1, and optimistically over 0.5). The usage of $`\mathrm{\Omega }_b`$ is particularly useful when we consider possibility of posthuman civilizations of truly intergalactic size. Currently there is some confusion about the exact value of $`\mathrm{\Omega }_b`$, since the primordial nucleosynthesis inferences apparently conflicts with conclusions drawn from the microwave background anisotropy observations (e.g. Kaplinghat and Turner 2001). It is to be expected that this problem will be solved soon, but in the meantime we note that the exact value is inessential for our purposes. If we restrict ourselves to the model (I), that is, ”posthumanity soon”, we may neglect spatial variation of $`q`$ and $`\mathrm{\Omega }_b`$ (which would be present in other cases, say $`q`$ is certainly different in intergalactic space from the one within galaxies), and using the definitional relation for $`\mathrm{\Omega }_b`$ we get

$$\mathrm{\Delta }Eqc^2\underset{V}{}\rho _b𝑑V.$$
(4)
Now, plugging this into (1), we obtain

$$I_{\mathrm{max}}=\frac{qc^2}{k_B\mathrm{ln}2}\frac{\underset{V}{}\rho _b𝑑V}{T}=3.5\times 10^{36}q\underset{V}{}\rho _b𝑑V,$$
(5)
in bits. If we wish to consider truly short-term posthuman civilization, we may state that the value of the integral is equal to

$$\underset{V}{}\rho _b𝑑V=2\times 10^{33}\frac{\overline{M}}{M_{}}n\mathrm{grams},$$
(6)
where $`n=1,2\mathrm{}`$ is the number of planetary systems controlled by our prototype posthuman civilization, and $`\overline{M}`$ is the average mass of a planetary system in the Galaxy. The fraction in (6) is likely to be less or about 1.5, when mass of the hidden matter (like comets in the Oort cloud whose total mass is still uncertain; see Weissman 1983) in our planetary system is taken into account, as well as the possibility of harvesting some interstellar matter between systems. Taking all this into account, we reach the estimate of

$$I_{\mathrm{max}}=7\times 10^{69}qn\frac{\overline{M}}{M_{}}\mathrm{bits}.$$
(7)
An analogous estimate can be obtained for the case (II) of posthumans in the far future. This case will certainly be prone to much larger uncertainties of not only quantitative, but also qualitative nature. Therefore, we shall here give just a sketch and postpone the detailed analysis to a forthcoming study. To this end, we may generalize the expression (5) taking into account the changes in both resources and methods available to advanced civilizations during long cosmological scales

$$I_{\mathrm{max}}(t)=\frac{q(t)c^2}{k_B\mathrm{ln}2}\underset{t_0}{\overset{t}{}}\frac{dt}{T(t)}\underset{V(t)}{}\rho _b𝑑V.$$
(8)
Here, $`V(t)`$ denotes the physical (proper) volume of space available to the advanced civilization under consideration, and $`T(t)`$ is the evolution of the cosmic temperature ($`T(t_0)=T_{\mathrm{CMB}}=2.730\pm 0.014`$ K). In this case, it is important to notice that the evolution of the equilibrium temperature of the universe in the context of physical eschatology is a problem not exactly solved to this day (approximate results for Einstein-de Sitter model are presented in Adams and Laughlin 1997), mainly because at very late epochs controversial sources (like the proton decay and Hawking radiation of an uncertain number of decaying black holes) come into play. Fortunately, we may be in position to know the asymptotic limit of the process of cooling of the universe, since it is determined by the temperature associated with the cosmological event horizon (Gibbons and Hawking 1977) in eq. (2). Since it can be shown that—if current estimates of the cosmological constant are correct—we have already entered the exponentially expanding (quasi)de Sitter phase (Ćirković and Bostrom 2000), we may use the simplest approximation of temperature decreasing as $`T(t)=T(t_0)\mathrm{exp}[H(tt_0)]`$, leading to the formal equalisation with (2) in $`t=t_0+(1/H_0)\mathrm{ln}[T(t_0)/T_{\mathrm{vac}}]t_0+69.58/H_0`$, i.e. in about 1140 Gyr (for $`h=0.6`$). As shown by Adams and Laughlin (1997; although only for the Einstein-de Sitter case), this is too pessimistic, since other photon sources will replace CMB as the determinants of the temperature of the universal heat bath. It will take probably many orders of magnitude larger time for the universe to reach the de Sitter temperature, but what is important is that there may be no further cooling. Thus, Dyson’s (1979) idea about using a special form of hybernation for the expression in eq. (1) to diverge is unfeasible. That said, we leave considerations of the detailed solutions of eq. (8) to a subsequent work, which will benefit from our increased understanding of the future thermal history of the universe.

## 3 Discussion

We have estimated (within an order of magnitude) the information processing power of advanced extraterrestrial or future posthuman communities. Our estimate is conservative in the sense that we have ignored the complicated and not yet fully understood issue of dissipationless computation (see, for instance, Porod et al. 1984) and assumed classical limit of $`kT\mathrm{log}2`$ dissipation per logical step. In addition, we have neglected the important issues of information transmission, noise and error correction (for a preliminary treatment of these, see the pioneering study of Sandberg 2000). All these issues should be covered in future, more realistic treatments.

It is worth noticing that numerical estimates reached above are quite conservative in comparison to the proposed fundamental information bounds, such as the holographic bound (e.g. Sussking 1995)

$$I_{\mathrm{max}}^h=\frac{A}{4L_{Pl}^2\mathrm{ln}2},$$
(9)
(where $`L_{Pl}\sqrt{G\mathrm{}/c^3}2\times 10^{33}`$ cm is the Planck length) or the Bekenstein (1973, 1981) bound

$$I_{\mathrm{max}}^B=\frac{2\pi RE}{c\mathrm{}\mathrm{ln}2},$$
(10)
where $`R`$ is the radius and $`E`$ the total mass-energy of the information cache. For instance, the equivalent of a holographic bound for the case of entire planetary systems from eq. (7) gives $`I_{\mathrm{max}}^h9.8\times 10^{76}(n\overline{M}/M_{})^2`$ bits, and the application of the Bekenstein bound to the same case yields $`I_{\mathrm{max}}^B5.13\times 10^{71}q(n\overline{M}/M_{})(R/1\mathrm{cm})`$ bits (where one should keep in mind that in order for this bound to be meaningful, the size of the memory $`R`$ has to be larger from its gravitational radius, which is about $`3\times 10^5`$ cm per Solar mass). Parenthetically, we note that most of the treatments of cosmological limitations on computation to be found in the literature (e.g. Tipler 1986; Sandberg 2000) and on the Internet use the Bekenstein bound, which is certainly more realistic and practical, but does not look entirely sound from the conceptual point of view. The reason for such a disadvantage of the Bekenstein bound compared to the holographic bound is that the Newtonian gravitational constant $`G`$ figures explicitly in the latter and not in the former. Since we believe that quantum gravitational degrees of freedom may contain a huge amount of information (and could be potentially exploited by advanced communities), a liberal approach would favor the usage of the holographic bound.

There are several reasons for pursuing this problem further. One is related to the confidence we may have in predictions of nonlinear dynamical systems by future computing devices. It is a well-known fact that (apparent) complexity of nonlinear systems increase at an exponential rate with time. No matter how advanced computing devices are employed in order to predict the behaviour of such systems, the prediction will break down at a particular timescale. One of the main factors determining this timescale is certainly the maximal amount of information which can be processed whatsoever during the computations necessary for prediction. But probably the most interesting issue to be solved by calculations such as these is the question of information cost of running large-scale detailed simulations of human environment. There are reasons to believe that advanced posthuman civilizations will run such simulations, which are sometimes aptly called ”ancestor simulations” (Prof. Nick Bostrom, manuscript in preparation). It is clear, however, that the depth (or complexity) of such a simulation would be very sensitive to the computing power and information resources of a posthuman ”director”. In order to assess viability of this scenario, one should attempt to reasonably predict those quantities. On the other hand, it remains for computer and cognitive scientists, as well as sociologists, to answer the deep question of minimal informational cost of any realistic simulation of human consciousness and society.

It is a pleasure to express gratitude to Prof. Nick Bostrom, whose deep insights into all sorts of futurological and philosophical questions has largely motivated the study of this topic. Srdjan Samurović is acknowledged for invaluable technical help.

References

Adams, F. C. and Laughlin, G.: 1997, Rev. Mod. Phys. 69, 337.

Barrow, J. D. and Tipler, F. J.: 1986, The Anthropic Cosmological Principle (New York: Oxford University Press).

Bekenstein, J.: 1973, Phys. Rev. D 7, 2333.

Bekenstein, J.: 1981, Phys. Rev. Lett. 46, 623.

Brillouin, L.: 1962, Science and Information Theory (New York: Academic Press).

Ćirković, M. M. and Bostrom, N.: 2000, Astrophys. Space Sci. 274, 675.

Dyson, F.: 1979, Rev. Mod. Phys. 51, 447.

Frautschi, S.: 1982, Science 217, 593.

Gibbons, G. W. and Hawking, S. W.: 1977, Phys. Rev. D 15, 2738.

Kaplinghat, M., Knox, L., and Turner, M. S.: 2000, Phys. Rev. Lett. 85, 3335.

Kaplinghat, M. and Turner, M. S.: 2001, Phys. Rev. Lett. 86 385.

Krauss, L. M. and Starkman, G.: 2000, Astrophys. J. 531, 22.

Livio, M.: 1999, Astrophys. J. 511, 429.

Lloyd, S.: 2000, Nature 406, 1047.

Peebles, P. J. E.: 1993, Principles of Physical Cosmology (Princeton: Princeton University Press).

Penrose, R.: 1989, The Emperor’s New Mind (Oxford: Oxford University Press).

Porod, W., Grondin, R. O., Ferry, D. K. and Porod, G.: 1984, Phys. Rev. Lett. 52, 232.

Rees, M. J.: 1969, Observatory 89, 193.

Sandberg, A.: 2000, Journal of Transhumanism 5 (available at
http://transhumanist.com/volume5/Brains2.pdf).

Staggs, S. T., Jarosik, N. C., Meyer, S. S., and Wilkinson, D. T.: 1996, Astrophys. J. 473, L1.

Susskind, L.: 1995, J. Math. Phys. 36, 6377.

Tipler, F. J.: 1986, Int. J. Theor. Phys. 25, 617.

Weissman, P. R.: 1983, Astron. Astrophys. 118, 90.


text:
# Untitled Document

Non–parametric Inference in Astrophysics
Larry Wasserman, Chris Miller, Bob Nichol, Chris Genovese, Woncheol Jang, Andy Connolly, Andrew Moore, Jeff Schneider and the PICA group. 111See www.picagroup.org for latest software, papers and memberships of the PICA group.

> We discuss non–parametric density estimation and regression for astrophysics problems. In particular, we show how to compute non–parametric confidence intervals for the location and size of peaks of a function. We illustrate these ideas with recent data on the Cosmic Microwave Background. We also briefly discuss non–parametric Bayesian inference.

1. Nonparametric Inference

The explosion of data in astrophysics provides unique opportunities and challenges. The challenges are mainly in data storage and manipulation. The opportunities arise from the fact that large sample sizes make nonparametric statistical methods very effective. Nonparametric methods are statistical techniques that make as few assumptions as possible about the process that generated the data. Such methods are inherently more flexible than more traditional parametric methods that impose rigid and often unrealistic assumptions. With large sample sizes, nonparametric methods make it possible to find subtle effects which might otherwise be obscured by the assumptions built into parametric methods. We begin by discussing two prototypical astrostatistics problems.

Problem 1. Density Estimation. Let $`X_1,\mathrm{},X_n`$ denote the positions of $`n`$ galaxies in a galaxy survey. Let $`f(x)dx`$ denote the probability of finding a galaxy in a small volume around $`x`$. The function $`f`$ is a probability density function, satisfying $`f(x)0`$ and $`f(x)𝑑x=1`$. We regard $`X_1,\mathrm{},X_n`$ as $`n`$ random draws from $`f`$. Our goal is to estimate $`f(x)`$ from the data $`(X_1,\mathrm{},X_n)`$ while making as few assumptions about $`f`$ as possible. Figure 1 shows redshifts from a pencil beam from the Sloan Digital Sky Survey. The figure shows several nonparametric density estimates that will be described in more detail in Section 3. The structure in the data is evident only if we smooth the data by just the right amount (lower left plot).222 The data involve selection bias since we can only observe brighter objects for larger redshifts. However, the sampling is fairly complete out to about $`z=0.2`$.

Problem 2. Regression. Figures 2 and 3 show cosmic microwave background (CMB) data from BOOMERaNG (Netterfield et al. 2001), Maxima (Lee et al. 2001) and DASI (Halverson 2001). The data consist of $`n`$ pairs $`(X_1,Y_1),\mathrm{},(X_n,Y_n)`$. Here, $`X_i`$ is multipole moment and $`Y_i`$ is the estimated power spectrum of the temperature fluctuations. If $`f(x)`$ denotes the true power spectrum then

$$Y_i=f(X_i)+ϵ_i$$

where $`ϵ_i`$ is a random error with mean 0. This is the standard regression model. We call $`Y`$ the response variable and $`X`$ the covariate. Other commonly used names for $`X`$ include predictor and independent variable. The function $`f`$ is called the regression function. The goal in nonparametric regression is to estimate $`f`$ making only minimal smoothness assumptions about $`f`$.

The main messages of this paper are: (1) with large data sets one can estimate a function $`f`$ nonparametrically, that is, without assuming that $`f`$ follows some given functional form; (2) one can use the data to estimate the optimal amount of smoothing; (3) one can derive confidence sets for $`f`$ as well as confidence sets for interesting features of $`f`$. The latter point is very important and is an example of where rigorous statistical methods are a necessity; the usual confidence intervals of the form “estimate plus or minus error” will not suffice.

The outline of this paper is as follows. Section 2 discusses some conceptual issues. Section 3 discusses kernel density estimation. Section 4 discusses nonparametric regression. Section 5 explains something that might be less familiar to astrophysicists, namely, nonparametric estimation via shrinkage. Section 6 discusses nonparametric confidence intervals. In Section 7 we briefly discuss nonparametric Bayesian inference. We make some concluding remarks in Section 8.

Notation: We denote the mean of a random quantity $`X`$ by $`E(X)`$, often written as $`X`$ in physics. The variance of $`X`$ is denoted by $`\sigma ^2Var(X)=E(XE(X))^2`$. A random variable $`X`$ has a Normal (or Gaussian) distribution with mean $`\mu `$ and variance $`\sigma ^2`$, denoted by $`XN(\mu ,\sigma ^2)`$, if

$$Pr(a<X<b)=_a^b\frac{1}{\sigma \sqrt{2\pi }}\mathrm{exp}\left\{\frac{1}{2\sigma ^2}(x\mu )^2\right\}𝑑x.$$

We use $`\widehat{f}`$ to denote an estimate of a function $`f`$.

2. Some Conceptual Issues

2.1. The Bias-Variance Tradeoff. In any nonparametric problem, we need to find methods that produce estimates $`\widehat{f}`$ of the unknown function $`f`$. Obviously, we would like $`\widehat{f}`$ to be close to $`f`$. We will measure closeness with squared error:

$$L(f,\widehat{f})=(f(x)\widehat{f}(x))^2𝑑x.$$

The average value of the error is called the risk or mean squared error (MSE) and is denoted by:

$$R(f,\widehat{f})=E_f\left[L(f,\widehat{f})\right].$$

A simple calculation shows that

$$R(f,\widehat{f})=\mathrm{Bias}_x^2𝑑x+\mathrm{Var}_x𝑑x$$

where $`\mathrm{Bias}_x=E[\widehat{f}(x)]f(x)`$ is the bias of $`\widehat{f}(x)`$ and $`\mathrm{Var}_x=Var[\widehat{f}(x)]=E[(\widehat{f}(x)E[\widehat{f}(x)])^2]`$ is the variance of $`\widehat{f}(x)`$. In words:

$$\mathrm{RISK}=\mathrm{BIAS}^2+\mathrm{VARIANCE}.$$

Every nonparametric method involves some sort of data-smoothing. The difficult task in nonparametric inference is to determine how much smoothing to do. When the data are over-smoothed, the bias term is large and the variance is small. When the data are under-smoothed the opposite is true; see Figure 4. This is called the bias-variance tradeoff. Minimizing risk corresponds to balancing bias and variance.

2.2. Nonparametric Confidence Sets. Let $`f`$ be the function of interest, for example, the true power spectrum in the CMB example. Assume that $`f`$ where $``$ is some very large class of functions. A valid (large sample) $`1\alpha `$ confidence set $`C_n`$ is a set $`C_n`$ such that

$$\underset{n\mathrm{}}{lim\; inf}\underset{f}{inf}Pr(fC_n)1\alpha $$

where $`n`$ is sample size. In words, $`C_n`$ traps the true function $`f`$ with probability approximately $`1\alpha `$ (or greater). In parametric models, confidence intervals take the form $`\widehat{\theta }\pm 2\mathrm{se}`$ where $`\widehat{\theta }`$ is an estimate of a parameter $`\theta `$ and se is the standard error of the estimate $`\widehat{\theta }`$. Bayesian interval estimates take essentially the same form. Nonparametric confidence sets are derived in a different way as we shall explain later in the paper.

If prior information is available on $`f`$ then it can be included by restricting $`C_n`$. For example, if it is thought that $`f`$ has at most three peaks and two dips, we replace $`C_n`$ with $`C_n`$ where $``$ is the set of functions with no more than three peaks and two dips.

Having constructed the confidence set we are then in a position to give confidence intervals for features of interest. We express features as functions of $`f`$, written $`T(f)`$. For example, $`T(f)`$ might denote the location of the first peak in $`f`$. Then

$$(\underset{fC_n}{inf}T(f),\underset{fC_n}{sup}T(f))$$

is a $`1\alpha `$ confidence interval for the feature $`T(f)`$. In fact, we can construct valid, simultaneous confidence intervals for many features of interest this way, once we have $`C_n`$. In section 6, we report such intervals for the CMB data.

Let us dispel a common criticism about confidence intervals. An oft cited but useless interpretation of a 95 per cent confidence interval is: if we repeated the experiment many times, the interval would contain the true value 95 per cent of the time. This interpretation leads many researchers to find confidence sets to be irrelevant since the repetitions are hypothetical. The correct interpretation is: if the method for constructing $`C_n`$ is used on a stream of (unrelated) scientific problems, we will trap the true value 95 per cent of the time. The latter interpretation is correct and is more scientifically useful than the former.

2.3. Where is the Likelihood? The likelihood function, which is a familiar centerpiece of statistical inference in parametric problems, is notably absent in most nonparametric methods. It is possible to define a likelihood and even perform Bayesian inference in nonparametric problems. But for the most part, likelihood and Bayesian methods have serious drawbacks in nonparametric settings. See section 7 for more discussion on this point.

3. Kernel Density Estimation.

We now turn to problem 1, density estimation. Let us start this section with its conclusion: the choice of kernel (smoothing filter) is relatively unimportant; the choice of bandwidth (smoothing parameter) is crucial; the optimal bandwidth can be estimated from the data. Let us now explain what this means.

Let $`X_1,\mathrm{},X_n`$ denote the observed data, a sample from $`f`$. The most commonly used density estimator is the kernel density estimator defined by

$$\widehat{f}(x)=\frac{1}{n}\underset{i=1}{\overset{n}{}}\frac{1}{h}K\left(\frac{xX_i}{h}\right)$$

where $`K`$ is called the kernel and $`h`$ is called the bandwidth. This amounts to placing a smoothed out lump of mass of size $`1/n`$ over each data point $`X_i`$. Excellent references on kernel density estimation include Silverman (1986) and Scott (1992).

The kernel is usually assumed to be a smooth function satisfying $`K(x)0`$, $`xK(x)𝑑x=0`$ and $`\tau x^2K(x)𝑑x>0`$. A fact that is well known in statistics but appears to be less known in astrophysics is that the choice of kernel $`K`$ is not crucial. The optimal kernel that minimizes risk (for large samples) is called the Epanechnikov kernel $`K(x)=.75(1x^2/5)/\sqrt{5}`$ for $`|x|<\sqrt{5}`$. But the estimates using another other smooth kernel are usually numerically indistinguishable. This observation is confirmed by theoretical calculations which show that the risk is very insensitive to the choice of kernel. In this paper we use the Gaussian kernel $`K(x)=(2\pi )^{1/2}e^{x^2/2}`$.

What does matter is the choice of bandwidth $`h`$ which controls the amount of smoothing. Figure 1 shows the density estimate with four different bandwidths. Here we see how sensitive the estimate $`\widehat{f}`$ is to the choice of $`h`$. Small bandwidths give very rough estimates while larger bandwidths give smoother estimates. Statistical theory tells us that, in one dimensional problems,

$`R(f,\widehat{f})`$ $`=`$ $`\mathrm{BIAS}^2+\mathrm{VARIANCE}`$
$``$ $`{\displaystyle \frac{1}{4}}h^4c_1A(f)+{\displaystyle \frac{c_2}{nh}}`$
where $`c_1=x^2K(x)𝑑x`$, $`c_2=K(x)^2𝑑x`$ and $`A(f)=(f^{\prime \prime }(x))^2𝑑x`$. The risk is minimized by taking the bandwidth equal to

$$h_{}=c_1^{2/5}c_2^{1/5}A(f)^{1/5}n^{1/5}.$$

This is informative because it tells us that the best bandwidth decreases at rate $`n^{1/5}`$ and leads to risk of order $`O(n^{4/5})`$. Generally, one cannot find a nonparametric estimator that converges faster than $`O(n^{4/5})`$. This rate is close to the rate of parametric estimators, namely, $`O(n^1)`$. The difference between these rates is the price we pay for being nonparametric.

The expression for $`h_{}`$ depends on the unknown density $`f`$ which makes the result of little practical use. We need a data-based method for choosing $`h`$. The most common method for choosing a bandwidth $`h`$ from the data is cross-validation. The idea is as follows.

We would like to choose $`h`$ to minimize the squared error $`(f(x)\widehat{f}(x))^2𝑑z=\widehat{f}^2(x)𝑑z2\widehat{f}(x)f(x)𝑑x+f^2(x)𝑑x`$. Since $`f^2(x)𝑑x`$ does not depend on $`h`$, this corresponds to minimizing

$$J(h)=\widehat{f}^2(x)𝑑z2\widehat{f}(x)f(x)𝑑x.$$

It can be shown that

$$\widehat{J}(h)=\widehat{f}^2(x)𝑑z2\frac{1}{n}\underset{i=1}{\overset{n}{}}\widehat{f}_i(X_i).$$

is an unbiased estimate of $`E[J(h)]`$, where $`\widehat{f}_i`$ is the “leave-one-out” estimate obtained by omitting $`X_i`$. Some algebra shows that

$$\widehat{J}(h)\frac{1}{hn^2}\underset{i}{}\underset{j}{}K^{}\left(\frac{X_iX_j}{h}\right)+\frac{2}{nh}K(0)$$
(1)
where $`K^{}(x)=K^{(2)}(x)2K(x)`$ and $`K^{(2)}`$ is the convolution of $`K`$ with itself. Hence, it is not actually necessary to compute $`\widehat{f}_i`$. We choose the bandwidth $`\widehat{h}`$ that minimizes $`\widehat{J}(h)`$. The lower left panel of figure 1 was based on cross-validation. An important observation for large data bases is that (1) can be computed quickly using the fast Fourier transform; see Silverman (1986, p 61-66).

4. Nonparametric Kernel Regression

Returning to the regression problem, consider pairs of points $`(X_1,Y_1),\mathrm{},(X_n,Y_n)`$ related by

$$Y_i=f(X_i)+ϵ_i.$$

The kernel method for density estimation also works for regression. The estimate $`\widehat{f}`$ is a weighted average of the points near $`x`$: $`\widehat{f}(x)=_{i=1}^nw_iY_i`$ where the weights are given by $`w_iK\left(\frac{xX_i}{h}\right)`$. This estimator is called the Nadaraya-Watson estimator. Figure 2 shows that estimator for the CMB data. Note the extreme dependence on the bandwidth $`h`$.

Once again, we use cross-validation to choose the bandwidth $`h`$. The risk is estimated by

$$\widehat{J}(h)=\frac{1}{n}\underset{i=1}{\overset{n}{}}(Y_i\widehat{f}_i(X_i))^2.$$

The first three panels in Figure 2 show the regression data with different bandwidths. The second plot is based on the cross-validation bandwidth. The final plot shows the estimated risk $`\widehat{J}(h)`$ from cross validation. Figure 3 compares the nonparametric fit with the fit by Wang, Tegmark and Zaldarriaga (2001).

Given the small sample size and the fact that we have completely ignored the cosmological models (as well as differential error on each data point) the nonparametric fit does a remarkable job. It “confirms,” nonparametrically, the existence of three peaks, their approximate positions and approximate heights. Actually, the degree to which the fit confirms the three peaks requires confidence statements that we discuss in section 6.

5. Smoothing by Shrinking

There is another approach to nonparametric estimation based on expanding $`f`$ into an orthogonal series. The idea is to estimate the coefficients of the series and then “shrink” these estimates towards 0. The operation of shrinking is akin to smoothing. These methods have certain advantages over kernel smoothers. First, the problem of estimating the bandwidth is replaced with the problem of choosing the amount of shrinkage which is, arguably, supported by better statistical theory than the former. Second, it is easier to construct valid confidence sets for $`f`$ in this framework. Third, in some problems one can choose the basis in a well-informed way which will lead to improved estimators. For example, Donoho and Johnstone (1994, 1995) and Johnstone (this volume) show that wavelet bases can be used to great advantage in certain problems.

Suppose we observe $`Y_i=f(x_i)+ϵ_i`$ where, for simplicity, we assume that $`x_1=1/n,x_2=2/n,\mathrm{},x_n=1`$. Further suppose that $`ϵ_iN(0,\sigma ^2)`$. Let $`\varphi _1,\varphi _2,\mathrm{}`$ be an orthonormal basis for $`[0,1]`$:

$$_0^1\varphi _j^2(x)𝑑x=1\mathrm{and}_0^1\varphi _i(x)\varphi _j(x)𝑑x=0\mathrm{when}ij.$$

For illustration, we consider the cosine basis: $`\varphi _1(x)1`$, $`\varphi _2(x)=\sqrt{2}cos(\pi x)`$, $`\varphi _2(x)=\sqrt{2}cos(2\pi x),\mathrm{}`$. Expand $`f`$ in this basis: $`f(x)_{j=1}^{\mathrm{}}\beta _j\varphi _j(x)_{j=1}^n\beta _j\varphi _j(x).`$ Estimating $`f`$ then amounts to estimating the $`\beta _j`$’s. Let $`Z_j=n^{1/2}_{i=1}^nY_i\varphi _j(i/n)`$. It can be shown that $`Z_jN(\theta _j,\sigma ^2),j=1,\mathrm{},n`$ where $`\theta _j=\sqrt{n}\beta _j`$. Once we have estimates $`\widehat{\theta }_j`$, we set $`\widehat{\beta }_j=n^{1/2}\widehat{\theta }_j`$ and $`\widehat{f}(x)=_{j=1}^n\widehat{\beta }_j\varphi _j(x).`$

How do we estimate $`\theta =(\theta _1,\mathrm{},\theta _n)`$ from $`Z=(Z_1,\mathrm{},Z_n)`$? A crude estimate is $`\widehat{\theta }_j=Z_j`$, $`j=1,\mathrm{},n`$. This leads to a very noisy (unsmoothed) estimate of $`f`$. Better estimates can be found by using shrinkage estimators. The idea – which goes back to James and Stein (1961) and Stein (1981) – is to estimate $`\theta `$ by shrinking the vector $`Z`$ closer to the origin. A major discovery in mathematical statistics was that careful shrinkage leads to estimates with much smaller risk. Following Beran (2000) we consider shrinkage estimators of the form $`\widehat{\theta }=(\alpha _1Z_1,\alpha _2Z_2,\mathrm{},\alpha _nZ_n)`$ where $`1\alpha _1\alpha _2\mathrm{}\alpha _n0`$ which forces more shrinkage for higher frequency cosine terms.

Let $`\alpha =(\alpha _1,\mathrm{},\alpha _n)`$ and let $`R(\alpha )`$ denote the risk of $`\widehat{\theta }`$ using shrinkage vector $`\alpha `$. An estimate of $`R(\alpha )`$, called Stein’s unbiased risk estimate (SURE), is

$$\widehat{R}(\alpha )=\underset{j}{}\left[\widehat{\sigma }^2\alpha _j^2+(Z_j^2\widehat{\sigma }^2)(1\alpha _j)^2\right]$$

where $`\sigma ^2`$ has been estimated by $`\widehat{\sigma }^2=\frac{1}{k}_{i=nk+1}^nZ_i^2`$ with $`k<n`$. Using appropriate numerical techniques, we minimize $`\widehat{R}(\alpha )`$ subject to the monotonicity constraint. The minimizer is denoted by $`\widehat{\alpha }`$ and the final estimate is $`\widehat{\theta }=(\widehat{\alpha }_1Z_1,\widehat{\alpha }_2Z_2,\mathrm{},\widehat{\alpha }_nZ_n)`$. Beran (2000) shows that the estimator obtained this way has some important optimality properties. Beran calls this approach REACT (Risk Estimation, Adaptation, and Coordinate Transformation). The estimated function $`\widehat{f}`$ turns out to be similar to the kernel estimator; due to space limitations we omit the plot.

6. Confidence Sets

When estimating a scalar quantity $`\theta `$ with an estimator $`\widehat{\theta }`$, it is common to summarize the uncertainty for the estimate by reporting $`\widehat{\theta }\pm 2se`$ where $`\mathrm{se}\sqrt{Var(\widehat{\theta })}`$ is the standard error of the estimator. Under certain regularity conditions, this interval is a 95 per cent confidence interval, that is,

$$Pr\left(\widehat{\theta }2\mathrm{se}\theta \widehat{\theta }+2\mathrm{se}\right).95.$$

This follows because, under the conditions alluded to above, $`\widehat{\theta }N(\theta ,\mathrm{se}^2)`$.

But the “plus or minus 2 standard errors” rule fails in nonparametric inference. Consider estimating a density $`f(x)`$ at a single point $`x`$ with a kernel density estimator. It turns out that

$$\widehat{f}(x)N(f(x)+\mathrm{bias},\frac{c_2f(x)}{nh})$$
(2)
where

$$\mathrm{bias}=\frac{1}{2}h^2f^{\prime \prime }(x)c_1$$
(3)
is the bias, $`c_1=x^2K(x)𝑑x`$ and $`c_2=K^2(x)𝑑x`$. The estimated standard error is

$$\mathrm{se}=\left\{\frac{c_2\widehat{f}(x)}{nh}\right\}^{1/2}.$$
(4)
Observe from (2) that $`(\widehat{f}(x)f(x))/seN(\mathrm{bias}/\mathrm{se},1)`$. If use the “estimate plus/minus 2 se” rule then

$`Pr\left(\widehat{f}(x)2\mathrm{se}f(x)\widehat{f}(x)+2\mathrm{se}\right)`$ $`=`$ $`Pr\left(2{\displaystyle \frac{\widehat{f}(x)f(x)}{\mathrm{se}}}2\right)`$
$``$ $`Pr\left(2N({\displaystyle \frac{\mathrm{bias}}{\mathrm{se}}},1)2\right).`$
If $`\mathrm{bias}/\mathrm{se}0`$ then this becomes $`Pr(2<N(0,1)<2).95`$. As we explained in Section 2, the optimal bandwidth is of the form $`h=cn^{1/5}`$. If you plug $`h=cn^{1/5}`$ this into (3) and(4) you will see that $`\mathrm{bias}/\mathrm{se}`$ does not tend to 0. The confidence interval will have coverage less than .95. In summary, “estimate plus/minus 2 standard errors” is not appropriate in nonparametric inference. There are a variety of ways to deal with this problem. One is to use kernels with a suboptimal bandwidth. This undersmooths the estimate resulting in a reduction of bias.

Another approach is based on the REACT method (Beran and Dumbgen, 1998). We construct a confidence set $`C_n`$ for the vector of function values at the observed data, $`𝐟_n=(f(X_1),\mathrm{},f(X_n))`$. The confidence set $`C_n`$ satisfies: for any $`c>0`$,

$$\underset{n\mathrm{}}{lim\; sup}\underset{𝐟_nc}{sup}|Pr(𝐟_nC_n)(1\alpha )|0$$

where $`a=\sqrt{n^1_ia_i^2}`$. The supremum is important: it means that the accuracy of the coverage probability does not depend on the true (unknown) function.

The confidence set, expressed in terms of the coefficients $`\theta `$, is

$$C_n=\{\theta :n^1\underset{j}{}(\theta _j\widehat{\theta }_j)^2\widehat{R}_r+n^{1/2}\widehat{\tau }z_\alpha \}$$

where $`z_\alpha `$ is such that $`P(Z>z_\alpha )=\alpha `$ where $`ZN(0,1)`$ and $`\widehat{\tau }`$ is a quantity computed from the data whose formula we omit here. Finally, the confidence set for $`f`$ is

$$𝒟_n=\{f:f=\underset{j}{}\beta _j\varphi _j:\beta _j=n^{1/2}\theta _j,\theta 𝒞_n\}.$$

Let us return to the CMB example. We constructed a 95 per cent confidence set $`C_n`$, then we searched over $`C_n`$ and found the possible number, location and heights of the peaks. We restricted the search to functions with no more than three peaks and two dips as it was deemed unlikely that the true power spectrum would have more than three peaks within the range of scales presently covered by the balloon experiments. Curves with just one or two peaks cannot be ruled out at the 95 per cent level i.e. they are still viable descriptions of the data but at a low statistical significance than three peaked models. The confidence intervals, restricted to three peak models, are as follows.

| Peak | Location | Height |
| --- | --- | --- |
| 1 | (118,300) | (4361,8055) |
| 2 | (377,650) | (1822,4798) |
| 3 | (597,900) | (1839,4683) |

The 95 per cent confidence interval for the ratio of the height of the second peak divided by the height of the first peak is $`(.21,1.4)`$. The 95 per cent confidence interval for the ratio of the height of the third peak divided by the height of the second peak is $`(.22,2.82)`$. Not surprisingly, the intervals are broad because the data set is small. The reader is referred to Miller et al (2002), for a more complete discussion of this work and our final results e.g. improvements in measurement error that are needed to get more precise confidence sets.

7. Nonparametric Bayes

There seems to be great interest in Bayesian methods in astrophysics. The reader might wonder if it is possible to perform nonparametric Bayesian inference. The answer is, sort of.

Consider estimating a density $`f`$ assumed to belong to some large class of functions such as $`=\{f:(f^{\prime \prime }(x))^2𝑑xC\}`$. The “parameter” is the function $`f`$ and the likelihood function is $`_n(f)=_{i=1}^nf(X_i).`$ Maximizing the likelihood leads to the absurd density estimate that puts infinite spikes on each data point. It is possible to put a prior $`\pi `$ over $``$. The posterior distribution on $``$ is well defined and Bayes theorem still holds:

$$Pr(fCX_1,\mathrm{},X_n)=\frac{_C_n(f)𝑑\pi (f)}{_{}_n(f)𝑑\pi (f)}.$$

Lest this seem somewhat abstract, take note that much recent work in statistics lately has led to methods for computing this posterior.

However, there is a problem. The parameter space $``$ is infinite dimensional and, in such cases, the prior $`\pi `$ is extremely influential. The result is that the posterior may concentrate around the true function very slowly. Worse, the 95 per cent Bayesian credible sets will contain the true function with very low frequency. In many cases the frequency coverage probability of the Bayesian 95 per cent credible set is near 0! Since high dimensional parametric models behave like nonparametric models, these remarks should give us pause before casually applying Bayesian methods to parametric models with many parameters.

The results that make these comments precise are fairly technical. The interested reader is referred to Diaconis and Freedman (1986), Barron, Schervish and Wasserman (1999), Ghosal, Ghosh and van der Vaart (2000), Freedman (2000), Zhao (2000) and Shen and Wasserman (2001). The bottom line: in nonparametric problems Bayesian inference is an interesting research area but is not (yet?) a practical tool.

8. Conclusion

Nonparametric methods are at their best when the sample size is large. The amount and quality of astrophysics data have increased dramatically in the last few years. For this reason, we believe that nonparametric methods will play an increasingly important role in astrophysics. We have tried to illustrate some of the key ideas and methods here. But we have really only touched on a few main points. We hope through our continued interdisciplinary collaboration and through others like it elsewhere, that the development of nonparametric techniques in astrophysics will continue in the future.

References

Barron, A., Schervish, M. and Wasserman, L. (1999). The consistency of posterior distributions in nonparametric problems. The Annals of Statistics. 27, 536-561.

Beran, R. (2000). REACT Scatterplot Smoothers: Superefficiency Through Basis Economy. Journal of the American Statistical Association.

Beran, R. and Dumbgen, L. (1998). Modulation Estimators and Confidence Sets. The Annals of Statistics, 26, 1826-1856.

Diaconis, P. and Freedman, D. (1986). On the consistency of Bayes estimates (with discussion). Ann. Statist. 14, 1-67.

Donoho, D.L. and Johnstone, I.M. (1994). Ideal spatial adaptation via wavelet shrinkage. Biometrika, 81, 425-455.

Donoho, D.L. and Johnstone, I.M. (1995). Adapting to unknown smoothness via wavelet shrinkage. Journal of the American Statistical Association, 90, 1200-1224.

Freedman, D. A. (1999). On the Bernstein-Von Mises theorem with infinite dimensional parameters. The Annals of Statistics, 27, 1119-1140.

Ghosal, S., Ghosh, J.K. and van der Vaart, A. (2000). Convergence rates of posterior distributions. The Annals of Statistics, 28, 500-531.

Halverson, N.W. et al. submitted to ApJ (2001) astro-ph/0104489

James, W. and Stein, C. (1961). Estimation with quadratic loss. In: Proc. Fourth Berkeley Symp. Math. Statist. Prob., 1, (J. Neyman, ed.) University of California Press, Berkeley. pp. 361–380.

Lee, A.T., et al. (2001) astro-ph/0104459

Miller, C. J., Nichol, R. C., Genovese, C., Wasserman, L., 2002, ApJ Letters, submitted (see astro–ph/0112049).

Netterfield, C.B., et al. (2001) astro-ph/0104460.

Scott, D.W. (1992). Multivariate Density Estimation: Theory, Practice, and Visualization. John Wiley, New York.

Shen, X. and Wasserman, L. (2001). Rates of convergence of posterior distributions. To appear: The Annals of Statistics.

Silverman, B.W. (1986). Density Estimation for Statistics and Data Analysis. Chapman-Hall: New York.

Stein, C. (1981). Estimation of the mean of a multivariate normal distribution. The Annals of Statistics, 9, 1135-1151.

Wang, X., Tegmark, M., and Zaldarriaga, M. (2001). To appear in Physcial Review D. astro-ph/0105091.

Zhao, L. H. (2000). Bayesian aspects of some nonparametric priors. The Annals of Statistics, 28, 532-552.