text:
# Understanding the long-term spectral variability of Cygnus X-1 with BATSE and ASM observations

## 1 Introduction

Cygnus X-1 is an archetypical black-hole X-ray binary discovered in 1964 (Bowyer et al. 1965). Due to its brightness and the persistency of emission, its observational and theoretical studies have been of great importance for understanding of the process of accretion onto black holes in general.

A specific important tool for testing theoretical models is the study of spectral variability. Pointed observations of Cygnus X-1 with X-ray and $`\gamma `$-ray telescopes have provided detailed information about its broad-band spectra (Gierliński et al. 1997, hereafter G97, Gierliński et al. 1999, hereafter G99, Di Salvo et al. 2001, hereafter D01, Frontera et al. 2001, hereafter F01, McConnell et al. 2002, hereafter M02), and its short time scale variability (e.g., Revnivtsev, Gilfanov, & Churazov 2000; Maccarone, Coppi, & Poutanen 2000). However, many black hole sources, including Cygnus X-1, also show dramatic spectral changes such as spectral transitions which occur on the time scales of weeks. Understanding the causes of these variations may require long uninterupted observations which are not available. The only observationally possible alternative is to monitor a source by all-sky monitors such as the ASM detectors aboard RXTE and the BATSE detectors aboard CGRO.

Cygnus X-1 is an excellent target for such studies. It has been monitored by the ASM since 1996, and by the BATSE from 1991 to 2000. During almost 5 years (1996–2000) the monitoring was simultaneously by both detectors. Spectral variability in the ASM data alone have been studied by Wen, Cui, & Bradt (2001), Reig, Papadakis, & Kylafis (2002) and Smith, Heindl, & Swank (2002), and some aspects of spectral variability in the BATSE data have been studied by Crary et al. (1996). Some partial analyses of joint ASM-BATSE variability have been done by Zhang et al. (1997, hereafter Z97) for the 1996 state transition, and by Brocksopp et al. (1999a). Here we present a comprehensive joint analysis of spectral variability in the full currently available ASM data and of all of the BATSE data.

## 2 Data analysis

To ensure both a sufficient overlap between the observations by ASM and by BATSE and sufficient statistics, we have chosen to use data binned on one-day timescale. The ASM data (average daily count rates and their standard errors) have been taken from the public database at xte.mit.edu/ASM\_lc.html. The data analyzed here cover the time interval from 1996 January 5 to 2002 June 6 (MJD 50087–52431), yielding 2127 days with usable measurements. In order to convert counts to energy units, we have constructed our own matrix based on comparison with pointed observations taken on the days with the ASM coverage, see Appendix A. In this way, we obtain energy fluxes in the photon energy, $`E`$, ranges of 1.5–3 keV, 3–5 keV, and 5–12 keV.

For the BATSE data, we generally followed the standard Earth occultation analysis methodology (Harmon et al. 2002). Counting rate spectra from individual occultation steps were averaged over one-day intervals. Energy fluxes in the 20–100 keV and 100–300 keV ranges were derived for each day by fitting a power-law spectrum (with variable photon spectral index, $`\mathrm{\Gamma }`$, defined by the photon flux $`E^\mathrm{\Gamma }`$) separately to the hardware channel ranges that corresponded most closely to the chosen energy intervals. To minimize interference from Cygnus X-3, we excluded data from days on which the angle between Cygnus X-3 and the closest point on the Earth’s limb at the time of Cygnus X-1 occultation was $`<2.5\mathrm{°}`$. We also excluded days on which fewer than four occultation steps were measured, since these are more likely to be affected by systematic errors.

The above procedure yields 2729 days with usable data from 1991 April 25 to 2000 May 22 (MJD 48371–51686). In Appendix A, we compare the BATSE data with those from the same pointed measurements as used for the ASM, and find them to be in a good general agreement.

The data sets above include 1187 days of usable data with simultaneous ASM-BATSE coverage. Using these data, we study here correlations between the ASM and BATSE fluxes and spectral slopes.

Often, hardness of spectra is described by the ratio of count rates or fluxes, $`R`$, in two adjacent channels. However, even if energy fluxes are used, this value depends on the width of the channels. To describe the hardness in a more objective way, we define here the effective photon power-law spectral index, $`\mathrm{\Gamma }`$, by the condition of the ratio of the energy fluxes in the power law spectrum, $`R`$, being equal to that in the observed spectrum. This index is analogous to colors used in stellar astrophysics. When two adjacent channels have boundary energies of $`E_1<E_2<E_3`$, the relationship between $`\mathrm{\Gamma }`$ and $`R`$ is given by,

$$R=\{\begin{array}{cc}\frac{E_3^{2\mathrm{\Gamma }}E_2^{2\mathrm{\Gamma }}}{E_2^{2\mathrm{\Gamma }}E_1^{2\mathrm{\Gamma }}},\hfill & \mathrm{\Gamma }2\text{;}\hfill \\ \frac{\mathrm{ln}E_3/E_2}{\mathrm{ln}E_2/E_1},\hfill & \mathrm{\Gamma }=2\text{.}\hfill \end{array}$$
(1)
The uncertainty of $`\mathrm{\Gamma }`$ is given by the propagation of errors, $`\mathrm{\Delta }\mathrm{\Gamma }=\mathrm{\Delta }R/|\mathrm{d}R/\mathrm{d}\mathrm{\Gamma }|`$ where $`\mathrm{\Delta }R`$ is the uncertainty of $`R`$. This allows us to obtain two effective indices for the ASM data, and one for BATSE. The BATSE effective index has a substantially better accuracy than either of the fitted BATSE indices.

## 3 Lightcurves

Figure 1a presents the 1.5–12 keV ASM and 20–300 keV BATSE energy flux light curves. We have rebinned the original ASM and BATSE one-day average data to assure that the error on each flux point is $`<10\%`$ and $`<20\%`$, respectively. This reduces the number of data points by only 5 and 7, respectively, and the relative accuracy of most of the values in Figure 1a is much better than the above ones. Figure 1b shows the corresponding hardness in the 3–12 keV and 20–300 keV bands, as given by the effective index of equation (1). We have chosen here the 3–12 keV index in order to avoid the effect of variable absorption by the line-of-sight matter, which affects mostly the spectrum below 3 keV.

In Figure 1, we see a soft-state transition in 1994, observed by both BATSE and OSSE (e.g., G99), very similar in form to the famous state transition of the year 1996 (Z97), as illustrated in more detail in Figure 2. Unfortunately, Cygnus X-1 was not observed at lower energies during this transition.

During the simultaneous operation of ASM and BATSE, we see the well-known overall anticorrelation of the fluxes from the two instruments. Interestingly, the 3–12 keV ASM index was rather similar to the BATSE one during the 1996 soft state. We also in Figure 1a a number of flares in the ASM flux which were not accompanied by any significant drops in the BATSE flux.

In addition to the 1996 soft state, we also see three shorter episodes when the ASM flux became high and its spectrum very soft, with $`\mathrm{\Gamma }(312\mathrm{keV})>2.4`$, around MJD 51850, 51900, and 51940, reported by Cui, Feng, & Ertmer (2002). More recently, another soft state lasting for about two months occured around MJD $`52150`$–52210. After this, Cygnus X-1 went to an intermediate \[i.e., still soft, with $`\mathrm{\Gamma }(312\mathrm{keV})2`$\] state for about a month. This was followed by yet another occurence of the soft state, which started around MJD 52230 and has lasted almost continuously (with some short intervals of an intermediate state around MJD 52350) for $`200`$ days so far, i.e., much longer than other state transitions of Cygnus X-1 observed so far. These recent state transitions are shown in detail in Figure 3. It is worth noticing that the spectrum during the last soft state was substantially softer than during the 1996 soft state, in the range of $`\mathrm{\Gamma }(312\mathrm{keV})3`$–3.5. Very interestingly, the spectrum since MJD 52160 until now (MJD 52431), i.e., for $`>270`$ days, has shown $`\mathrm{\Gamma }(312\mathrm{keV})2`$, i.e., much softer than the usual $`\mathrm{\Gamma }(312\mathrm{keV})1.7`$. Unfortunately, no monitoring at higher energies was available after MJD 51686, including those highly interesting state transitions.

We note that there have been recent reports of short (with duration $`1`$ day), very strong, flares from the direction of Cygnus X-1 measured at photon energies $`20`$ keV, at which range the flux exceeded the typical, hard-state, flux by a factor of $`10`$ (Stern, Beloborodov, & Poutanen 2001; Golenetskii et al. 2002). These flares occured on MJD 49727, 49801, 51289 and 52329. We have looked into characteristics of the lightcurves around those days, and have found, unfortunately, no common characteristic that would allow some insight into the nature of the flares. For example, the last event occured in a soft state whereas the three previous ones, in hard states.

## 4 Spectral Correlations

### 4.1 The ASM data

Based on the variability from the ASM, we distinguish here three main spectral states. We define the soft and hard states by the conditions (see Fig. 1) of $`\mathrm{\Gamma }(312\mathrm{keV})>2.4`$, $`<2.1`$, respectively, and an intermediate state by $`\mathrm{\Gamma }(312\mathrm{keV})=2.1`$–2.4. Our definition of the soft state corresponds, e.g., to MJD between approximately 50230 and 50307 (1996 May 27–Aug. 12). These definitions may differ somewhat from other ones used in the literature.

Wen et al. (1999) found strong correlations of the ASM flux and hardness with the orbital phase within the hard state, with a peak of the hardness and a dip in the flux around the phase zero (i.e. the companion star in front of the X-ray source), see also Brocksopp et al. (1999a). This effect appears to be due to absorption by a partially ionized wind from the companion. On the other hand, no phase modulation has been found in the soft state. In order to distinguish between this effect and an intrinsic spectral variability, we separately consider the hard state data in the orbital phase between 0.85 and 0.15 and outside this phase interval. We use the ephemeris of Brocksopp et al. (1999b).

We find the logarithms of the fluxes in the three ASM channels strongly correlate with each other, linearly within each state. However, the coefficients for each of the linear dependence are different in each of the three states, as well as the ones for the 1996 soft state are different than those for the 2000–02 soft states. This type of variability is most conveniently shown on flux-index diagrams (analogous to color-magnitude diagrams in stellar astrophysics), Figures 4a, b.

In Figures 4a, b, we see two distinct branches: one formed by the hard state, and the other one corresponding to the soft states, connected by the intermediate state. The hard state is slightly extended towards low fluxes and hard spectra due to local absorption close to the phase zero. Since the bound-free opacity of cosmic medium generally decreases towards increasing energies, the phase effect on the index and flux in the 3–12 keV range is significantly weaker than on those in the 1.5–5 keV one. However, it is clear that absorption can by no means explain the very strong hardness-flux anticorrelation in the hard state. A similar hardness-flux anticorrelation is commonly seen in Seyferts (e.g., Nandra et al. 2000), including some radio-loud ones (e.g., Zdziarski & Grandi 2001).

We note that the negative correlation in the 1.5–5 keV range is stronger than that in the 3–12 keV range, but the latter correlation is still relatively strong. This behavior is consistent with the spectrum pivoting at an energy above the ASM range. As we find by analyzing the BATSE data, the pivot point is at an energy between 20 and 100 keV. On the other hand, the horizontal spread in Figures 4a, b corresponds to a second variability pattern, in which the X-ray spectrum simply moves up and down without changing its shape. These two variability patterns are explained theoretically in §7.1.

On the other hand, the hardness-flux correlation is positive in the soft states, as found before by Wen et al. (2001). A simple interpretation of this type of variability is a two-component spectrum, with a constant component at low energies (e.g., a blackbody or a disk blackbody) and a variable high-energy tail. Then, the stronger the tail, the higher the flux and the less of the blackbody cutoff is observed, making the spectrum harder. A theoretical interpretation of this behavior is given in §7.2.

We note that the existence of two main branches in the index-flux dependence is independent of the conversion from counts to energy units. If we plot the analogous figures for the raw ASM counts, the obtained shapes are similar to those in Figure 4, see Reig et al. (2002).

Given the above dependencies, there is a strong linear correlation between the spectral indices in the 3–12 keV and 1.5–5 keV energy ranges. Figure 5 shows a comparison of the two quantities, analogous to color-color diagrams often shown for X-ray binaries. However, unlike neutron-star binaries showing changes in the sign of this correlation (e.g., Schulz, Hasinger, & Trümper 1989), here the two hardnesses (or $`\mathrm{\Gamma }`$) are positively correlated in all states, with only the linear coefficients being state-dependent. Thus, the X-ray variability of Cygnus X-1 is only on the diagonal branch of its color-color diagram.

Interestingly, the 2000–02 soft states show distinctly different dependencies between the two indices and the 1.5–5 keV index and flux from those for the 1996 soft state. In Figures 4a and 5, we see that $`\mathrm{\Gamma }(1.55\mathrm{keV})`$ is generally lower (i.e., the spectrum is harder) and the flux is lower at a given hardness during the 2000–02 soft states. On the other hand, Figure 4b shows that the 3–12 keV fluxes and slopes during the two periods of the soft state are rather similar, but extending to lower fluxes and softer spectra. Also, the flux-hardness correlations in the later soft states remain positive. We discuss theoretical explanation of the differences between the two soft states in §7.3 below.

By inspection of Figure 1, we see that Cygnus X-1 in its hard state can spend hundreds of days (e.g., MJD 50700–50900) with its X-ray flux changing only by a factor of two and $`\mathrm{\Gamma }(312\mathrm{keV})`$ changing only by $`\pm 0.1`$. If a random time interval of 100 days is chosen within the hard state, the corresponding flux-index points would, in most cases, occupy a relatively small region within the diagonal hard-state strips in Figures 4a, b. The relatively large spread of the indices and fluxes within the hard state seen in these figures is mostly due to the flares and dips in the long-term X-ray light curve. Thus, although the overall anticorrelation between the flux and hardness is extremely strong, it is achieved only on very long time scales, of at least a year. The uniformity of this anticorrelation also indicates that the minor flares/dips in the light curve still belong to the hard state, indicating some common underlying physics (see §7.1). On the other hand, flux changes by a factor of two (corresponding to the horizontal spread of the region covered by the hard-state data in Fig. 4) occur on a time scale of tens of days or less. This explains the lack of a noticeable hardness-flux correlation within time intervals of 5.6 day (i.e., the orbital period) found by Wen et al. (2001).

In the soft states (especially the 1996 one), the correlated flux-index changes occur on a much shorter time scale of a day or so, as clearly seen in Figure 1. Thus, although the overall soft-state flux-hardness correlation is much weaker than the hard-state anticorrelation in Figures 4a, b, the former can be measured on much shorter time intervals, in agreement with the very strong correlation within 5.6-day intervals found by Wen et al. (2001).

### 4.2 The BATSE data

Now, we consider the BATSE data alone. We have found that the spectral state can be best defined for these data using the 100–300 keV flux (see Fig. 2). We then define the hard state by $`F(100300\mathrm{keV})>3.0`$ keV cm-2 s-1, and the soft state by $`F(100300\mathrm{keV})<2.2`$ keV cm-2 s-1, with the intermediate state in between. These definitions are in good agreement with those adopted for the ASM during the period of the simultaneous ASM-BATSE coverage (e.g., the 1996 soft state in the BATSE data corresponds to MJD 50232–50305, very similar to the range obtained from the ASM data in §4.1).

Figure 6a shows the effective 20–300 keV spectral index as a function of the 20-100 keV flux. In the hard state, we still see some anticorrelation between the hardness and flux, qualitatively similar to those in the 1.5–12 keV range. However, the anticorrelation is only weak, and, in particular, much weaker than the one for the 3–12 keV range (which, in turn, was weaker than that for the 1.5–5 keV range). We then see in Figure 6b that the 20–300 keV hardness becomes strongly positively correlated with the 100–300 keV flux. This behavior shows that the variable spectrum pivots at an energy $`100`$ keV.

The positive hardness-flux correlation in the 100–300 keV range in the hard state is also seen for the fitted spectral index in that range, as shown in Figure 7. Also, almost no correlation is seen for these quantities in the 20–100 keV range. This shows that the pivoting effect is not an artifact of our use of the effective 20–300 keV index. On the other hand, the fitted indices bear rather large statistical uncertainties. This prevents seeing any trends in the soft states (Fig. 7), where the BATSE fluxes are low and the number of days limited.

Similarly to the case of the ASM data, substantial variations in the range covered by the source on the $`\mathrm{\Gamma }`$-$`F`$ diagram during the hard state occur on long time scales. If, for example, we consider a typical 100-day interval during the hard state, the corresponding data points would cover a relatively small part of the cyan region in Figures 6a, b.

Figure 8 shows then the relationship between the 20–100 keV and 100–300 keV fluxes. We see that, in spite of the pivoting of the spectrum in the hard state in the BATSE energy range, the two fluxes are still positively correlated. This appears to be due to two effects. First, the pivot point is below 100 keV, and a part of the 20–100 keV flux will still positively be correlated with the 100–300 keV flux. Second, in addition to the pivoting pattern of variability, there is also variability with the spectrum simply moving up and down, as discussed in §4.1.

We also note that BATSE hardness ratio, $`F(100300\mathrm{keV})/F(20100\mathrm{keV})`$ \[related to the effective spectral index used here via eq. (1)\] shows significant variations by a factor of 3 in the hard state. This supports earlier findings of Stern et al. (2001).

In the soft state, $`\mathrm{\Gamma }(20300\mathrm{keV})`$ is only weakly changing and roughly independent of $`F(20100\mathrm{keV})`$, see Figure 6a. This shows that the 20–100 keV spectrum simply moves up and down with relatively small changes of the spectral shape, consistent with our interpretation of the variability in the ASM range being caused by the amplitude variability of the high-energy tail. In the 20–100 keV range, only the tail (moving up and down) is observed. This behavior is also consistent with the strong positive correlations of the two BATSE fluxes seen in Figure 8.

Finally, there appears to be also a weak positive correlation between $`\mathrm{\Gamma }(20300\mathrm{keV})`$ and $`F(100300\mathrm{keV})`$ in the 1996 soft state, see Figure 6b. This indicates that the spectrum above 100 keV softens at low flux levels (as indeed seen directly in the lightcurve, see Fig. 2). This points to the existence of a second variability pattern in the 1996 soft state, consisting of a softening of the slope of the tail when the tail amplitude is low.

### 4.3 ASM-BATSE correlations

We consider now the combined ASM-BATSE data. In general, virtually all quantities are strongly correlated with each other, in spite of the different specific times at which the two instruments measure the emission of Cygnus X-1 within a given day. This shows that the intraday variability is relatively weak compared to the variability on longer time scales. We define here the states based on $`\mathrm{\Gamma }(312\mathrm{keV})`$ (as in §4.1).

In Figure 9b, we compare the fluxes in the lowest and the highest energy channels, i.e., 1.5–3 keV vs. 100–300 keV. We see that they are anticorrelated in the hard state. On the other hand, the 100–300 keV flux is completely independent of the soft X-ray flux in the soft state. In Figure 9a, we see that the 20–100 keV flux becomes weakly positively correlated with the 1.5–3 keV flux in the hard state, while it remains independent of it in the soft state.

These correlations in the hard state are consistent with the spectrum pivoting at an energy of $`100`$ keV. In the soft state, the lack of correlation of the BATSE fluxes with the 1.5–3 keV one shows that the existence of a constant component at low energies and a variable and independent component at high energies. These findings are consistent with those obtained above while considering the ASM and BATSE data separately.

Figures 10a, b present the ASM-BATSE flux-index diagrams analogous to those for Figures 4 for the ASM data alone. Figure 10a shows that, in the hard state, the dependence of the BATSE hardness on the ASM flux is similarly anticorrelated as that of the ASM index on the ASM flux (Fig. 4). In the soft state, $`\mathrm{\Gamma }(20300\mathrm{keV})`$ is uncorrelated with the ASM flux, as well as it shows a substantial spread of values.

On the other hand, the 1.5–5 hardness is now positively correlated in the hard state with $`F(100300\mathrm{keV})`$ (Fig. 10b), opposite to the case of the ASM data alone. This is again consistent with existence of a pivot point below 100 keV, reversing the correlation with respect to that with a flux at energies $`100`$ keV. On the other hand, there is a moderate hardening of the ASM spectrum with the increasing $`F(100300\mathrm{keV})`$ in the soft state. This is explained by changing of the amplitude of the high energy tail with a constant soft, blackbody component which dominates the ASM flux.

Figure 11 shows an ASM-BATSE index-index correlation. The 3–12 keV hardness is positively correlated with the 100–300 keV one in the hard state, showing that the spectrum changes its local slope in the same direction at all measured energies. The BATSE index is substantially softer than that of the ASM, showing that the spectrum at a given time softens with increasing energy. On the other hand, there is no correlation in the soft state. This means that the tail component varies independently of the blackbody one.

## 5 Fractional Variability

An independent way of quantifying the variability of Cygnus X-1 is to consider the fractional variability as a function of energy. The intrinsic (i.e., after relatively small corrections for the measurement errors) root mean-square (rms) variability based on one-day averages in the entire data sets is 89%, 48%, 26%, 21%, 30% in the 1.5–3, 3–5, 5–12, 20–100, and 100–300 keV energy ranges, respectively (see relatively similar values obtained by Z97 for a period including the 1996 state transtion). However, this variability includes the state transitions. As we have found in §4 above, the spectral variability indicates different physical mechanisms operating in the two main states. Thus, we have computed the rms variability separately in the hard state (defined as above and excluding the phase range of 0.85–0.15 for the ASM data) and in the 1994, 1996 and 2000–02 soft states. The results are shown in Figure 12.

The rms dependence in the hard state is consistent with our finding of a pivot point in the 20–100 keV range, where the relative variability has the minimum. The maximum variability is in the softest channel, consistent with the variability being driven by a variable input of soft seed photons (see §7). On the other hand, the dependence in the 1996 soft state is consistent with a constant soft component and a variable tail with almost constant shape. Also, the 1994 soft state (not observed by the ASM) has very similar values of the rms, confirming our previous finding of the two states being very similar to each other. On the other hand, the 2000–02 soft states show an offset of the rms values with respect to those of the 1996 soft state. We have checked that it is not due to varying properties of the relatively large number of soft episodes during 2000–02. In particular, the rms values for only the most recent soft state (MJD $`>52250`$) are even larger than those shown on Figure 12. A possible explanation for the offset is a decrease of the characteristic energy of the constant soft component (see §7.3), which then moves the rms pattern to lower energies. In addition, the rms of the high-energy tail in the recent soft states appears higher than that in the 1996 soft state.

In the hard state, the shown values of rms correspond to the frequency range of $`10^8`$$`10^5`$ Hz. It is highly interesting to compare the power in the variability on those long time scales with that corresponding to short time scales. The rms values corresponding to seven RXTE/PCA observations in the hard state are given by Lin et al. (2000). They obtain the rms for 0.002–10 Hz in the 2–50 keV energy range of $`20`$–40% (with a relatively flat energy dependence, consistent with the variability at a constant spectral shape), i.e., values very similar to those obtained by us. On the other other hand, the high-frequency hard-state power spectrum per log of frequency of Cygnus X-1 has an overall maximum (usually with two nearby peaks) in the $`0.1`$–10 Hz range and it declines fast both towards lower and higher frequencies (e.g., Gilfanov, Churazov, & Revnivtsev 1999). This implies that the broad-band hard-state spectrum of Cygnus X-1 has at least two maxima with comparable power, one at $`1`$ Hz, and one somewhere below $`10^5`$ Hz.

## 6 Pointed Observations

Further insight into the character of spectral variability of Cygnus X-1 can be gained from comparing spectra from pointed observations taken at different times. Here we discuss information that can be derived from spectral variability between those observations and luminosities in different states.

### 6.1 Spectra

Figure 13 shows three broad-band spectra in the hard state, two in the soft state, and one in the intermediate state. The triangles and circles in Figures 411 correspond to four of those pointed spectra (1996 May 30–31, June 22, September 12, 1998 May 3–4).

In the hard state, we see two patterns of spectral variability. First, the 1991 June (black) and 1998 May (blue) spectra are almost parallel to each other, and vary only in the normalization being different by a factor of $`1.5`$. (Note that this 1991 spectrum is thus at about highest daily flux observed by BATSE.) Such a pattern was also seen within the 1991 observation, when the four spectra taken within one day (with each observation lasting $`2`$ hours) varied only in the normalization within a factor of $`2`$ (G97). We can identify this type of variability with the spread in the $`F`$-$`\mathrm{\Gamma }`$ correlations in §4 at constant values of $`\mathrm{\Gamma }`$, e.g., clearly seen in both the ASM (Fig. 4) and BATSE (Fig. 6) data.

Second, the 1996 September hard-state spectrum (cyan) is substantially softer than the other hard-state spectra and its high-energy cutoff is at a visibly lower energy. In a thermal Comptonization fit (see §7 below), $`kT57`$ keV whereas $`kT86`$ keV and 77 keV in the 1991 June 6 and 1998 May 3–4 spectra, respectively. The cyan and blue spectra intersect at $`50`$ keV. We can identify this type variability with the hardness-flux anticorrelation seen in the ASM data (Fig. 4), and the anticorrelation between the soft and hard X-ray fluxes (Fig. 9b).

In the soft and intermediate states, we see the tail (in the red, magenta and green spectra) moves up and down with relatively small changes in the spectral shape. The soft X-ray peak is almost the same in the two soft-state spectra (red and magenta), but the form of the intermediate spectrum above 3 keV implies a decrease of that peak in the intermediate state. Also, when the tail moves up, the X-ray spectrum hardens. These effects are consistent with the behavior of the soft and intermediate state ASM spectra (Fig. 4).

### 6.2 Bolometric fluxes and the nature of the state transitions

One of the outstanding puzzles of Cygnus X-1 has been the small difference, claimed to be $`50`$–70% (Z97) between the bolometric fluxes in the hard and the soft states. This requires fine-tuning of the hard-state accretion rate as well as fine-tuning of models of the transition (Esin et al. 1998). However, we stress here that the estimate by Z97 was based on their estimates of the effects of absorption and the bolometric fluxes based solely on the ASM and BATSE data.

To make further progress, we have compiled in Table 1 the intrinsic bolometric fluxes of Cygnus X-1 based on a number of highly accurate broad-band pointed observations, including those described in §6.1. Those data provide much better estimates of the unabsorbed bolometric flux than those possible with the ASM-BATSE data. We see that whereas the typical flux in the hard state is $`4`$$`5\times 10^8`$ erg cm-2 s-1, it is $`1.5\times 10^7`$ erg cm-2 s-1 in the soft state, i.e., $`3`$–4 times higher. We also note that the two spectra with the highest hard-state fluxes correspond to either the peak of intraday variability or a relatively soft continuum. Also, the flux in the intermediate state is just between those in the two main states. The fact that the total observed flux during the 1996 state transition remains almost unchanged (see Fig. 1), as pointed out by Z97, can be reconciled with the results in Table 1 by noting that whereas the effect of absorption on the measured bolometric flux is negligible in the hard state (with the spectrum peaking at $`100`$ keV, Fig. 13), it becomes crucial in the soft state, and very difficult to estimate with the ASM data alone (with the peak of the $`EF_E`$ spectrum at $`1`$ keV).

This result makes models of the state transitions based on a change of the accretion rate, $`\dot{M}`$, much more viable than with the previous result of Z97. In particular, a relatively large fraction of the accretion flow in the hard state advected (e.g., Narayan & Yi 1995) or forming outflows (e.g., Blandford & Begelman 1999), i.e., not efficiently radiating, becomes now possible.

## 7 Theoretical Interpretation

We interpret here the spectral variability of Cygnus X-1 in terms of Comptonization of soft X-ray blackbody seed photons. We use a Comptonization code by Coppi (1992, 1999), with its present xspec (Arnaud 1996) version (eqpair) described in G99. This model was also used to fit both states of Cygnus X-1 by Poutanen & Coppi (1998), and some states of GRS 1915+105 by Zdziarski et al. (2001). In general, the electron distribution can be purely thermal or hybrid, i.e., Maxwellian at low energies and non-thermal at high energies in the presence of an acceleration process. The electron distribution, including the temperature, $`kT`$, is calculated self-consistently from the luminosities in Comptonized photons, $`L_{\mathrm{hard}}`$, and in the blackbody seed photons irradiating the cloud, $`L_{\mathrm{soft}}`$. The plasma optical depth, $`\tau `$, includes a contribution from $`e^\pm `$ pairs (which may be negligible or not). The importance of pair production depends on the ratio of the luminosity to the characteristic size, $`r`$, which is usually expressed in dimensionless form as the compactness parameter, $`\mathrm{}L\sigma _\mathrm{T}/(rm_\mathrm{e}c^3)`$ (where $`\sigma _\mathrm{T}`$ is the Thomson cross section and $`m_\mathrm{e}`$ is the electron mass). We also include Compton reflection (Magdziarz & Zdziarski 1995) and Fe K fluorescent lines. We refer the reader to G99 for details.

We fit this model to the spectra from pointed observations (§6.1). Using the obtained models, we calculate the energy fluxes in the bands corresponding to those of the ASM and BATSE. The fluxes from the pointed observations have uncertainties negligible compared to those of the monitoring instruments and are very weakly dependent on the assumed model, given the energy resolution of the pointed instruments being much less than the width of any of the ASM-BATSE energy bands. Then, in order to compare theoretical variability patterns with those observed, we change some parameters in some of the best-fit models, and recalculate the fluxes. For the fitted BATSE indices, we create simulated BATSE data with models of pointed observations and theoretical models, and fit them with a power law spectrum. We use two BATSE responses corresponding to the extreme orientations of the detector with respect to the source and average the resulting indices. However, the differences in the two values of $`\mathrm{\Gamma }`$ are $`0.03`$, i.e., negligibly small.

### 7.1 The hard state

In the hard state, the Comptonizing plasma is nearly thermal (G97; Poutanen 1998; Zdziarski 2000). The electron distribution may contain a weak high-energy tail, but its contribution to spectral formation below several hundred keV is negligible (M02). If the compactness is low enough for pair production being negligible, the main parameters of a model are $`\tau `$, $`L_{\mathrm{hard}}`$, and $`L_{\mathrm{soft}}`$.

We have found that we can explain the anticorrelation of the X-ray flux and hardness (§4.1) by a variable luminosity of seed blackbody photons, $`L_{\mathrm{soft}}`$, irradiating a cloud of thermal plasma with an approximately constant energy dissipation rate (given by $`L_{\mathrm{hard}}`$) and a constant $`\tau `$. Under these conditions, the plasma electron temperature, $`kT`$, will adjust itself to the variable seed flux as to satisfy energy balance. Then, the higher the seed flux, the lower $`kT`$, and (since a decrease of $`kT`$ at a constant $`\tau `$ corresponds to an increase of $`\mathrm{\Gamma }`$) the softer the X-ray spectrum. Since the Comptonized spectrum joins at low energies to the peak of the spectrum of the seed photons, this will also correspond to an increase of the soft X-ray flux. The pivoting of the spectrum will occur somewhere below 100 keV causing the decrease of the high energy flux above 100 keV.

We have applied this model to the two hard-state observations of Cygnus X-1 by BeppoSAX (1996 Sept. 12, F01; 1998 May 3–4, D01), see §6.1. As discussed in those papers, Cygnus X-1 shows a soft X-ray excess in addition to the main thermal Comptonization component. The soft excess is well modeled by an additional Comptonization region. The seed photons are blackbody at a temperature of $`0.15`$ keV. This model provides a good fit to either of the data sets, as shown in Figure 13. The position of these data in Figures 411 are shown by the filled triangles and circles, respectively. We note that there may be some offsets of the BeppoSAX points with respect to the corresponding ones from ASM and BATSE due to calibration uncertainties (see Appendix A) and intraday variability.

We then vary only $`L_{\mathrm{soft}}`$ for each of the Comptonization regions, increasing and decreasing it by a factor of 2, and keep the compactness at a low enough value for pair production to be negligible. We do not adjust the model normalization. The predicted spectra for the 1998 May observation are shown in Figure 14a. We see the pivot point at $`50`$–90 keV. The resulting predictions for the fluxes and the spectral hardness are shown in Figures 411 by the solid black lines originating at the filled symbols (corresponding to the actual observations). In the case of the 1996 September observation, we show only the effect of the decreasing $`L_{\mathrm{soft}}`$ (i.e., hardening of the X-ray spectrum) since the observed spectrum is already rather soft, see Figure 13. Given the simplicity of our assumption, the agreement with the observed spectral variability is rather good, yielding support to the picture of the X-ray variability driven by a variable flux of irradiating soft X-ray photons. A possible accretion geometry corresponding to this variability pattern is shown in Figure 15 (see, e.g. Poutanen, Krolik, & Ryde 1997; Esin et al. 1998, Zdziarski 1998), and its discussion in given in §7.4.

The above model is the same as the one used to explain the X-ray spectral variability of 3C 120 by Zdziarski & Grandi (2001). Note that the pivot obtained for 3C 120 is at $`7`$–10 keV, which difference is due to both the power-law component of 3C 120, $`\mathrm{\Gamma }1.9`$, softer than $`\mathrm{\Gamma }1.7`$ typical for the hard state of Cygnus X-1, and the lower energy of the seed photons in AGNs than in black-hole binaries. Specifically, UV seed blackbody photons with a temperature of 10 eV were assumed for 3C 120, compared to $`150`$ eV fitted for the hard state of Cygnus X-1.

We note that an increase/decrease of $`\mathrm{\Gamma }`$ in this model corresponds to a substantial decrease/increase of $`kT`$, in the $`60`$–100 keV range in the case shown in Figure 14a. Such an effect is seen when we compare the two hard-state BeppoSAX spectra. The softer and the harder one have $`kT57`$ keV and 77 keV, respectively. However, we do not know how universal this correlation is, and the BATSE data are of insufficient accuracy to test this effect.

On the other hand, the hot Comptonizing plasma may be $`e^\pm `$ pair dominated. If it is mostly thermal, no signature of annihilation will be visible in the spectra (Zdziarski 1986; Maciołek-Niedźwiecki, Zdziarski, & Coppi 1995). Thus, we cannot detect the presence of pairs just by the spectral shape. However, their presence will have an effect on the spectral variability. Namely, a variable $`L_{\mathrm{soft}}`$ will cause changes in the pair production rate causing changes in the total Thomson optical depth, $`\tau `$. Pairs act as a thermostat (e.g., Malzac, Beloborodov, & Poutanen 2001). We have fitted a pair-dominated model (achieved at a high enough compactness, $`\mathrm{}500`$ and assuming a nonthermal fraction of 5% consistent with the results of M02) to the 1998 May data, and then changed $`L_{\mathrm{soft}}`$ by a factor of 2 as in the case of $`e^{}p`$ plasma. Indeed, the equilibrium $`kT`$ stayed almost constant at 72 keV in all three cases. The resulting variability pattern is shown in Figure 14b. We see that the obtained spectra indeed show no annihilation features (while pair annihilation is included in the eqpair model). The predictions of this model are shown by the dotted lines connected to the filled circle in Figures 411. We see that the present data are not sensitive enough to determine which of the two cases is preferred.

As we discussed in §4, there has to be another variability pattern in the hard state when the total luminosity varies with a more or less constant spectral shape. If the X-rays in Cygnus X-1 are produced in a hot inner accretion flow (as in Fig. 15), a higher accretion rate (and higher luminosity) corresponds to a higher optical depth, $`\tau `$, of the hot disk. If the geometry of the flow does not change, i.e., the feedback from the cold disk is constant, the seed photon luminosity (from the outer cool disk) is linearly proportional to the total luminosity. The resulting spectrum is then almost independent of luminosity, except that the high energy cutoff decreases with an increase of $`\tau `$ since the electron temperature will be adjusted to keep the Comptonization parameter $`y4\tau kT/m_\mathrm{e}c^2`$ constant (see e.g. Poutanen 1998 for a review).

We simulate this variability pattern by increasing and decreasing the total luminosity by $`\sqrt{2}`$ starting from the same BeppoSAX spectra of Cygnus X-1 and, at the same time, adjusting $`\tau `$ (with no pairs) according to predictions of models of hot accretion disks. If $`y=`$ constant, $`\tau L^{2/7}`$ and $`L^{1/6}`$ in the advection and cooling dominated cases (Zdziarski 1998). Here we assumed the former. A decrease of $`kT`$ with increasing $`L`$ at a constant $`\mathrm{\Gamma }`$ (since $`y`$ is constant) is consistent with the four data sets of G97, as well as it has been observed in two bright states of the black hole binary GX 339–4 (Zdziarski 2000; Wardziński et al. 2002). As we noted above this variability pattern should have a smaller amplitude comparing with the first pattern. The resulting spectra are shown in Figure 14c. The observed spectra that show similar variability are the hard state spectra (the black and blue ones) in Figure 13. The resulting predictions for the fluxes and the spectral slopes are shown in Figures 411 by the dashed lines connected to the filled circle. We can see that the two patterns describe well the full range of observed correlations.

We note that these two variability patterns give opposite predictions for the dependence of $`kT`$ on the flux at $`200`$ keV. In the case of varying $`L_{\mathrm{soft}}`$, $`kT`$ either increases or remains constant when flux increases, and the opposite is true for the case of varying the total $`L`$.

### 7.2 The 1996 soft state

For the 1996 soft state, we use the ASCA-RXTE observation (G99) and the BeppoSAX-CGRO one (F01; M02). In those papers, the data were already fitted with the eqpair model, in which a hybrid, thermal-nonthermal plasma Comptonizes seed, blackbody photons coming from an accretion disk (with the maximum blackbody temperature of $`kT_{\mathrm{soft}}`$. Now a substantial part ($`0.7`$ in those two observations) of the dissipation in the Comptonizing cloud goes into particle acceleration, in addition to thermal heating (see G99 for details). The points corresponding to those observations are shown by the open triangles and circles, respectively, in Figures 411.

Then, we vary $`L_{\mathrm{hard}}`$ by a factor of 1.5 up and down while keeping $`L_{\mathrm{soft}}`$ constant. The predicted spectra are shown in Figures 16a, b. This qualitatively explains the spectral variability in the soft-state region in Figures 411 as shown by the solid lines extending from the open symbols (corresponding to the actual observations). This variability pattern corresponds to a stable optically-thick accretion disk (Shakura & Sunyaev 1973) and variable flares on its surface. A likely corresponding accretion geometry is shown in Figure 17, and more discussion in given in §7.4. The disk in Cygnus X-1 is thermally and viscously stable (G99), supporting this picture. The model of constant disk and a variable tail component has also been found by Churazov, Gilfanov, & Revnivtsev (2001) to fit RXTE/PCA data for Cygnus X-1.

However, the above variability pattern of variable $`L_{\mathrm{hard}}`$ cannot explain the correlation of the BATSE spectral slope with the 100–300 keV flux (Figs. 2 and 6b) as well as substantial variability of that slope (Fig. 10a). We note that a softening of the $`50`$–500 keV spectrum with decreasing flux in the same energy range has also been observed by OSSE between two observations in the soft states in 1994 and 1996 (G99).

We have tested a number of models to explain this secondary variability pattern, with variable either the power-law index of the accelerated electrons, the fraction of the energy which goes to the acceleration, and the variable maximum Lorentz factor, $`\gamma _{\mathrm{max}}`$, of the accelerated electrons. We have found the first pattern does not provide a good fit to the data at all, whereas the last one provides the best fit. Changes of $`\gamma _{\mathrm{max}}`$ have little effect on the slope or the flux at low energies (in the ASM range), but markedly affect the BATSE fluxes and the slope. Starting from the best-fit model for the 1996 June data, we decrease $`\gamma _{\mathrm{max}}`$ from $`10^3`$ to 10. The resulting spectrum is shown by the dotted curve in Figure 16c. The corresponding predictions for the correlation diagrams are shown for the 1996 June spectrum only by the dotted lines in Figures 411.

### 7.3 The 2000–02 soft states

As shown in §4.1 and 5, the 2000–02 soft states show different spectral and timing properties from those of the 1996 soft state. In particular, Figure 12 shows an offset between the two ASM rms dependencies. In the framework of the above soft-state model (§7.2) consisting of two components, a constant disk soft blackbody and a variable tail, the rms offset is stronly suggestive of a reduction of the disk color temperature. Then, the variable tail starts to dominate at correspondingly lower energies. From Figure 12, we estimate the reduction of the characteristic color temperature, $`kT_{\mathrm{soft}}`$, to be by a factor of $`1.5`$–2.

This inference is, in fact, supported by the ASM colors. As shown in Figures 4 and 5, the later soft states have harder 1.5–5 keV spectra, as well as of a lower flux at a given hardness. This is consistent with a decrease of $`kT_{\mathrm{soft}}`$, in which case the 1.5–5 keV range has now a larger contribution from the tail, which has a harder spectrum than the high-energy cutoff of the disk blackbody, dominant in 1996 up to $`3`$ keV (see Fig. 16b).

We have modeled this behaviour using the eqpair model. Starting from the best-fit model for the 1996 June data, we decreased $`kT_{\mathrm{soft}}`$ from 370 eV (the best fit value) to 200 eV. This, indeed, has led to a hardening of the 1.5–5 keV spectrum. However, we also found that the 3–12 keV hardens as well, contrary to the data. This implies some change of the parameters of the Comptonizing plasma in addition to a reduction of $`kT_{\mathrm{soft}}`$. Unfortunately, the ASM data alone hardly constrain the shape of the Comptonization spectrum. Here, we found we can account for the ASM data by a purely phenomenological assumption of setting the reflection fraction to zero. We stress, however, that the actual character of the change of the hot plasma can only be determined by broad-band (from soft to hard X-rays) pointed observations (which were not available to us). E.g, the power-law index of the accelerated electrons, the fraction of the energy which goes to the acceleration, or $`\gamma _{\mathrm{max}}`$ could have changed.

The resulting effect on the 1.5–5 and 3–12 keV slopes is shown by the dashed line originating at the open circle in Figures 45. We see that indeed we can qualitatively account for the changes of the color and fluxes from the 1996 soft state to the 2000–02 ones. The corresponding spectrum is shown by the dashed curve in Figure 16c. We show this spectrum only for energies $`20`$ keV, since we cannot make any predictions for higher energies based on the ASM data alone.

Note that we adjusted the normalization of the spectrum to match the data. We also note the decrease of the disk temperature is not accompanied by a correspondingly large decrease in the luminosity, as seen, e.g., in Figures 1 and 3. This is possible if the temperature decrease was not caused by an increasing accretion rate but by a change of some physical conditions of the disk causing the color correction (the ratio of the color temperature to the effective one) to decrease. Any more detailed explanation of the changes requires pointed observations in the soft X-ray range.

The spectral variability within the 2000–02 soft states appears to follow the primary pattern found for the 1996 soft state, i.e., changes of $`L_{\mathrm{hard}}`$ at a constant $`L_{\mathrm{soft}}`$. We do not show this pattern separately in Figures 45 for the sake of clarity.

### 7.4 Discussion

Table 2 summarizes our findings above. We stress that the primary variability patterns in the hard and soft states are modelled by opposite behaviors of the hard and soft luminosities. This argues agains spectral models of the state transition related to a quantitative change of one parameter while keeping the geometry constant, e.g., variations in the relative fraction of the power dissipated in a disk corona (e.g., Young et al. 2001). Although the overall spectrum in the hard and soft states is indeed hard and soft, respectively, their spectral variability patterns are completely different.

In the soft state, the spectral variability is indeed fully consistent with the presence of flares (forming a corona) on the surface of an accretion disk, as schematically shown in Figure 17. The luminosity in the flares is strongly variable while that in the underlying disk is approximately constant. On the other hand, it is likely that the power in the flares is derived from dissipation in the underlying disk. In that case, $`L_{\mathrm{hard}}+L_{\mathrm{soft}}`$ would be approximately constant rather than just $`L_{\mathrm{soft}}`$. We have modelled this variability pattern as well, and have found it completely consistent with the data but barely distinguishable from the one for constant $`L_{\mathrm{soft}}`$ shown in Figures 411. Thus, we do not show it for the sake of simplicity. The 1996 and 2000–02 soft states both follow the variability pattern driven by changes of $`L_{\mathrm{hard}}`$, but the latter show the disk blackbody temperature lower than that of the former (see §7.3).

In contrast, the primary variability in the hard state corresponds to a variable flux in the seed photons while keeping the power supplied to electrons in the hot plasma approximately constant. In our opinion, the most likely geometry appears to be a hot inner flow surrounded by an overlapping thin, optically-thick, accretion disk supplying the soft photons, as schematically shown in Figure 15. Changes in the inner radius of the disk cause then variable flux of seed photons. The seed photons are likely to be mostly due to reprocessing of the hard radiation. As noted in §4, this variability takes place on rather long time scales, which are apparently needed to change the outer disk. On the other hand, changes in the local $`\dot{M}`$ in the hot disk cause the overall spectrum to move up and down, but with a constant shape determined by the approximately constant geometry.

An increase in the global $`\dot{M}`$ then causes a decrease in the inner radius of the outer disk, which is compatible with some physical scenarios (e.g., Esin et al. 1998). The data do allow a modest increase of $`L_{\mathrm{hard}}`$ (instead of its constancy) with increasing $`\dot{M}`$, but the present data are not sensitive enough to test it. When the thin disk moves close to the minimum stable orbit, the hot accretion flow collapses and is replaced by the active regions above the thin disk.

On the other hand, it is also possible that the hot plasma forms active coronal regions above a disk also in the hard state. The variable $`L_{\mathrm{soft}}`$ can be then related to changes in the bulk velocity of the active regions/flares (Beloborodov 1999; Malzac et al. 2001) and a constant power dissipated in the flares causes then the constancy of $`L_{\mathrm{hard}}`$. However, the reversal of this variability pattern in the soft state needs then to be explained.

We note that, for the sake of simplicity, we have neglected the well-known correlation between the spectral index and Compton reflection (Zdziarski, Lubiński, & Smith 1999; Gilfanov et al. 1999, 2000). These changes of the reflected component would have a relatively minor effect on our results. We stress, however, that the above correlation is completely consistent with the flux-hardness anticorrelation in the hard state. Namely, the cold medium providing the soft seed photons is most likely the same one as the reflecting medium (see Fig. 15). Then an increase of the flux of soft photons due to reprocessing of the hard radiation would also be accompanied by an increase of the reflection strength. Thus, the presence of one of those correlations makes it likely that the other one is also present.

## 8 Conclusions

We have presented the history of the flux and hardness of Cygnus X-1 from the entire course of its monitoring by the BATSE, and up to 2002 June by the ASM. In particular, we show in detail the ASM data from the recent soft states (MJD 51840 till 52431).

We find then numerous correlations between various ASM-BATSE fluxes and indices. In the hard state, we find the X-ray flux at $`<12`$ keV and $`>100`$ keV to be negatively and positively, respectively, correlated with the hardness in the same energy range. We find this is compatible with the variable slope of the overall spectrum and a pivot at $`50`$ keV. This behavior is well modeled by thermal Comptonization in a hot plasma with an approximately constant power irradiated by a variable flux of soft seed photons.

In all of the soft states, the X-ray flux at $`<12`$ keV is positively correlated with the hardness (as earlier found by Wen et al. 2001). However, in the case of the 1996 soft state, the correlation disappears above 20 keV, where the spectrum undergoes flux changes at an approximately constant shape. This is well modeled by hybrid, thermal/nonthermal Comptonizationin a hot plasma with a variable power irradiated by soft seed photons with a constant flux, i.e., the behavior opposite to that in the hard state. In addition, we find the overall 1.5–300 keV spectrum is moving up and down in the hard state, and the high energy tail in the soft state softens at low flux levels. The four patterns are summarized in Table 2.

We also compute the fractional variability amplitude as a function of photon energy in the two states. In the hard state, it has a minimum at the 20–100 keV range, consistent with the presence of the spectral pivot at $`50`$ keV, and the maximum at the 1.5–3 keV range, consistent with the variability being driven by a variable soft flux. In the soft state, the fractional variability is lowest in the 1.5–3 keV range and increases strongly towards higher energies, consistent with the presence of a stable soft (disk) component and a variable high-energy tail.

Our preferred model of the accretion flows in the hard state is an inner hot flow overlapping with an outer thin accretion disk with a variable inner radius. That variability is responsible for varying flux of seed photons and changes of the spectral slope. However, a dynamic corona above a cold disk is also possible. In both scenarios, variations of the local accretion rate at the constant geometry are then responsible for the changes of the luminosity occuring at an approximately constant shape. In the soft state, the data are fully compatible with variable nonthermal flares occuring above a stable thin disk.

We also study in detail differences between the 1996 and 2000–02 soft states. We find the ASM spectra of the latter to be harder in the 1.5–5 keV range and to have the rms variability significantly higher. Both effects can be modeled by a decrease of the color disk temperature from the 1996 soft state to the 2000–02 ones. On the other hand, we find the 1994 soft state (monitored by BATSE only) to have the same spectral and timing properties at energies $`>20`$ keV as the 1996 soft state.

Furthermore, we present a compilation of detailed broad-band spectra of Cygnus X-1 from pointed observations in the hard, soft, and intermediate states. The spectral variability observed in this data set fully supports our conclusions based on the ASM-BATSE data. Furthermore, we find the bolometric luminosity in the spectra from the pointed observations increases by a factor of $`3`$–4 from the hard state to the soft state. This amplitude of the changes of $`L`$ strongly supports models of the state transitions based on a change of the global $`\dot{M}`$.

This research has been supported by grants from KBN (5P03D00821, 2P03C00619p1,2) and the Foundation for Polish Science (AAZ). JP and AAZ acknowledge support from the Royal Swedish Academy of Sciences and the Polish Academy of Sciences through exchange programs. We thank the ASM and BATSE teams for their efforts in carrying out the monitoring, C. Done and T. Di Salvo for providing us with their BeppoSAX data, M. Gierliński for making the original versions of Figures 15 and 17, and G. Wardziński and P. Lachowicz for help with some aspects of the data analysis.

## Appendix A Comparison with Pointed Observations

The ASM results are given in counts per second, $`R_i`$, in the three adjacent channels, $`i=a,b,c`$, which nominal boundaries are 1.5, 3, 5, 12 keV. In general, conversion of those rates into energy flux units is uncertain, and depending on the spectral shape and the position of the source on the sky. Here, we have adopted a simple conversion based on comparison of fluxes from six pointed observations with the ASM average daily fluxes taken on the same days. The pointed observations are listed in Table 1, and correspond to MJD 50226 (RXTE, G99), 50233–50234 (ASCA-RXTE, G99), 50256, 50338 (BeppoSAX, F01), 50936–50937 (BeppoSAX, D01), MJD 51091 (BeppoSAX, C. Done, private communication). The first five of them are also shown in Figure 13.

We have then assumed a linear dependence of the energy flux from each channel on the count rate in the same one and in a neighbouring one. This takes into account the dependence of the conversion on the hardness of the incident spectra. We have then least-square fitted the conversion coefficients. In each of the three fits, we have given the same relative weight to each of the measurements, in order to avoid the fit being dominated by the data points with the highest flux (i.e., in the soft state). In the case of the channel $`b`$, we have found that allowing the dependence of the conversion on both neighbouring channel did not improve the fit quality with respect to the dependence on the count rate in the $`a`$, $`b`$ channels only.

In this way, we have obtained,

$$\left(\begin{array}{c}F_a\\ F_b\\ F_c\end{array}\right)=\left(\begin{array}{ccc}0.2552& 0.0324& 0\\ 0.0357& 0.2507& 0\\ 0& 0.0533& 0.3591\end{array}\right)\left(\begin{array}{c}R_a\\ R_b\\ R_c\end{array}\right)$$
(A1)
where $`F_i`$ and $`R_i`$ are in units of keV cm-2 s-1 and cm-2 s-1, respectively. We see in Figures 18a, b that the resulting overall agreement is very good, given that some of the differences have to be systematic, caused by the pointed observations being not taken exactly on the same times as the ASM observations and by calibration differences between different instruments.

The BATSE fluxes in energy units have been obtained in the way described in §2. Figure 19 shows the comparison with the pointed observations listed above. We see a generally good agreement between the two data sets.


format: text

text:
# Qualitative Analysis of Correspondence for Experimental Algorithmics

## 1 Introduction

Correspondence is a ubiquitous concept in the interpretation of spatial datasets. Correspondence establishes analogy, indicating objects that play similar roles with respect to some context. For example, correspondence between template and image features supports object recognition, correspondence among isobars in a weather map aids identification of pressure troughs and ridges, and, as shown in this paper, correspondence and lack thereof among level curves in datasets like Fig. 1 supports characterization of matrix properties for scientific computing applications.

This paper develops a novel qualitative analysis mechanism for extracting and utilizing correspondence in spatial data, with particular focus on model selection in applications with sparse, noisy data. We develop our mechanism in the context of experimental algorithmics, an emerging field that designs methodologies for empirically evaluating algorithms on realistic test problems, and interpreting and generalizing the results to guide selection of appropriate mathematical software. It is the preferred method of analysis in applications where domain knowledge is imperfect and for which our understanding of the factors influencing algorithm applicability is incomplete. For instance, when solving linear systems associated with finite-difference discretization of elliptic partial differential equations (PDEs), there is little mathematical theory to guide a choice between, say, a direct solver and an iterative Krylov solver plus preconditioner. An experimental algorithmics approach is to parameterize a suitable family of problems, and mine a database of PDE “solves” to gain insight into the likely relative performance of these two approaches . Thus experimental algorithmics is data-driven and aims to empirically capture the relationship between problem characteristics and algorithm performance.

Many tasks in experimental algorithmics involve assessing the eigenstructure of a given matrix. Eigenstructure helps characterize the stability, sensitivity, and accuracy of numerical methods as well as the fundamental tractability of problems. Recently, the spectral portrait (e.g. Fig. 1) has emerged as a tool for graphically visualizing eigenstructure. A spectral portait characterizes how the eigenvalues of a matrix change as perturbations (e.g. due to numerical error) are introduced in computations involving the matrix. Level curves in a spectral portait correspond to perturbations, and the region enclosed by a level curve contains all possible eigenvalues that are equivalent with respect to perturbations of a given magnitude. Analysis of level curves with respect to a class of perturbations reveals information about the matrix (e.g. nonnormality and defective eigenvalues) and the effects of different algorithms and numerical approximations.

Experimental algorithmics thus determines high-level properties from analysis of data from low-level computational experiments; we focus here on the particular case where such properties are extracted from graphical representations and where it is necessary to minimize the computational experiments performed (owing to the cost and complexity of conducting them). We pose the extraction of high-level properties as a model selection problem and show that correspondence can be exploited to drive data collection in order to discriminate among possible models. Each step in our framework is parameterized by domain knowledge of properties such as locality, similarity, and correspondence. Therefore the framework is generic both with respect to a variety of problems in experimental algorithmics (two case studies are presented here), and to problems in other domains (e.g. connections to weather data analysis are discussed). Furthermore, the approach leads to efficient and explainable data collection motivated directly by the need to disambiguate among high-level models.

## 2 Qualitative Analysis of Spatial Data

Our mechanism for qualitative analysis of correspondence is based on the Spatial Aggregation Language (SAL)  and the ambiguity-directed sampling framework . SAL programs apply a set of uniform operators and data types (Fig. 2) in order to extract multi-layer geometric and topological representations of spatial data. These operators utilize domain knowledge of physical properties such as continuity and locality, specified as metrics, adjacency relations, and equivalence predicates, to uncover regions of uniformity in spatially distributed data. Ambiguity-directed sampling focuses data collection so as to clarify difficult choice points in an aggregation hierarchy. It seeks to maximize information content while minimizing the number and expense of data samples.

As an example of aggregation, consider the construction of level curves in the spectral portrait in Fig. 1. The key steps in such an analysis (after ) are:

* The input field maps sample locations (e.g. on a uniform grid) to perturbation levels, representing allowable imprecision before a point becomes indistinguishable from an eigenvalue (computational details are discussed later).
* Aggregate points in a neighborhood graph (e.g. Delaunay triangulation or regular grid), localizing computation to spatially proximate points.
* Interpolate values at new locations from values at neighboring samples, in this case determining locations of points with perturbation level belonging to a discrete set (e.g. $`10^1`$, $`10^2`$, $`\mathrm{}`$).
* Classify neighboring similar-enough objects into equivalence classes with an equivalence predicate, in this case testing equality of field value.111While an implementation such as marching squares might combine interpolation, aggregation, and classification, we view them as conceptually distinct operations.
* Redescribe each equivalence class of lower-level objects as a single higher-level object, in this case abstracting connected points into curves. The curve can be represented more compactly and abstractly (e.g. with a spline) than its set of sample points.

As a consequence of redescription, the next aggregation level can treat curves as first-class objects, and aggregate, classify, and redescribe them, for example to find curves nested around a single eigenvalue. This higher-level process uses the same operators, but with different parameters specifying locality, equivalence, and abstraction. Ambiguity arises when, for example, not enough sample points are available to be confident in a curve’s location, or in the separation of two curves. Ambiguity-directed sampling then optimizes selection of new locations (e.g. near the ambiguity) for data collection in order to clarify the decision-making.

As this example illustrates, SAL and ambiguity-directed sampling provide a suitable vocabulary (e.g. distance and similarity metrics) and mechanism (bottom-up aggregation and top-down sampling) to uncover multi-level structures in spatial data sets. Successful applications include decentralized control design  weather data analysis , analysis of diffusion-reaction morphogenesis , and identification of pockets in $`n`$-dimensional space and decomposition of a field based on control influences .

## 3 Qualitative Analysis of Correspondence

Our correspondence mechanism builds on the relationship between lower-level and higher-level objects in a SAL hierarchy (refer again to Fig. 2). The mechanism has two key steps: (1) establish analogy as a relation among lower-level constituents of higher-level objects; (2) establish correspondence between higher-level objects as an abstraction of the analogy between their constituents. For example, in object recognition, analogy might match image and template features, and correspondence might abstract the analogy as a rigid-body transformation. Similarly, in level curve analysis, analogy might match sample points on neighboring curves by location and local curvature, and correspondence might abstract the match as a parameterized deformation of spline representations of the curves. The analogy between constituents is well-defined only because of the context of the higher-level objects; higher-level correspondence then captures a more global view of the local matches.

Tab. 1 outlines our correspondence mechanism. Traditional SAL operators collect and abstract groups of lower-level objects into higher-level objects, and establish pairs of higher-level objects for which correspondence is to be considered. Two pieces of domain knowledge are then applied:

indicates pairs (and confidence in the pairs) of lower-level objects that are analogous with respect to the higher-level objects they comprise. Examples include testing feature values for each pair of lower-level objects, hashing indices in the higher-level objects’ local coordinate systems, or explicitly constructing a spatial relation such as a triangulation. The predicate can enforce bijective analogy if appropriate. The analogize operator in Tab. 1 applies a predicate $`a`$ for the constituent objects $`l_1h_1`$ and $`l_2h_2`$ of neighboring (by $`g_h`$) higher-level object pairs $`h_1`$ and $`h_2`$; it returns a graph labeling the confidence of each analogous object pair.

abstracts an analogy relation on lower-level objects into a description of higher-level object correspondence. Examples include redescribing point-wise matches as a rigid-body transformation or a parameterized deformation. The abstraction captures confidence in the overall correspondence, for example computing root-mean squared distance or a Hausdorff metric between locations or features (e.g. local curvature) of analogous objects. Correspondence optimization thus entails adjusting the underlying analogy to maximize the abstracted confidence. Finally, while nearby constituents of one object often match nearby constituents of another, discontinuity (e.g. when an outer curve envelops two inner ones in Fig. 1) is a noteworthy event detectable with the abstraction mechanism. The correspond operator in Tab. 1 performs correspondence abstraction from an analogy $`g_l`$, applying a function $`c`$ to the subgraph of $`g_l`$ on the constituent objects $`\mathrm{}(h_1)`$ and $`\mathrm{}(h_2)`$ of neighboring (by $`g_h`$) higher-level objects $`h_1`$ and $`h_2`$; it returns a graph labeling each such higher-level object pair with the resulting correspondence spatial object.

An aggregation/correspondence hierarchy establishes a distribution of possible high-level models for an input instance, thereby posing a model selection problem: choose the one that (e.g. in a maximum-likelihood sense) best matches the data. Our mechanism supports model selection in two key ways. (1) The operators estimate and optimize confidence in correspondence. Since correspondence implies mutual support among parts of a model, it can allow relatively high-confidence model selection even with sparse, noisy data. (2) The operators bridge the lower-/higher-level gap. This allows weaknesses and inconsistencies detected in higher-level correspondence to focus lower-level data collection to be maximally effective for model disambiguation.

## 4 Applications in Experimental Algorithmics

We present two case studies applying our analysis framework to experimental algorithmics problems. For each we describe the underlying numerical analysis problem, our particular solution approach, and results. To the best of our knowledge, these are the first systematic algorithms for performing complete imagistic analyses (as opposed to relying on human visual inspection ), and which focus data collection and evaluate models until a high-confidence model is obtained.

### 4.1 Matrix Spectral Portrait Analysis

Our first case study focuses on the previously introduced task of matrix spectral portrait analysis (Fig. 1). Formally, the spectral portrait of a matrix $`𝒜`$ is defined as:

$$𝒫(z)=\mathrm{log}𝒜_2(𝒜zI)^1_2,$$
(1)
where $`I`$ is the identity matrix. The singularities of this map are located at the eigenvalues of the matrix, and the analysis determines the sensitivity of computation to numerical imprecision by analyzing how the map decreases moving away from the eigenvalues. As discussed in the introduction, the region enclosed by a level curve of a certain magnitude contains all points that act as “equivalent” eigenvalues under perturbations of that magnitude. For example, Fig. 1 illustrates that, under large enough perturbation, eigenvalues at 3 and 4 are indistinguishable; after larger perturbations, their joint contour also surrounds the eigenvalue at 2, and then the eigenvalue at 1. This illustrates that sensitivity analysis is solved by identifying values at which level curves merge.

Tab. 2 describes qualitative correspondence analysis of spectral portraits. Data are collected by computing Eq. 1; the analysis determines perturbation equivalence of eigenvalues by detecting curve merges via level curve correspondence. The first aggregation level generates samples on a coarse regular grid around the known eigenvalue locations, and then interpolates and abstracts level curves as in the spatial aggregation example. The second aggregation level finds correspondence among these curves from a Delaunay triangulation analogy of their constituent points. It tracks correspondence outward from the eigenvalues to establish a model of merge events, establishing confidence in such a model by good analogy (most points match) and good curve separation (regularly-spaced samples separate the two curves, providing evidence that contours do not merge at a smaller perturbation). Ambiguity-directed sampling generates additional data in order to separate curves and to ensure that each eigenvalue pair is merged at some perturbation.

Fig. 3 demonstrates the application of the mechanism to the companion matrix of the polynomial $`(x1)^3(x2)^3(x3)^3(x4)`$  (see also Fig. 1). The initial grid has at least one sample between each eigenvalue (in the example, a resolution of 0.5) and extends one unit beyond the bounding box of the eigenvalues. While correspondence is found even with this coarse grid (left of Fig. 3), the confidence is not high (only a few samples separate the curves), and no curve surrounds all eigenvalues. Ambiguity-directed sampling computes additional points on a finer, larger grid (right of Fig. 3), yielding the correct results. We have applied this same mechanism to a variety of polynomial companion matrices with different numbers and spacings of roots; in each case, the correspondence mechanism correctly identifies the correct model with high confidence after 1-3 subsamples and 1-3 grid expansions.

### 4.2 Qualitative Computation of Jordan Forms

Our second case study focuses on analysis of the Jordan decomposition of a matrix. Matrix decomposition is an important technique, revealing pertinent features of a matrix and supporting algorithmic techniques in areas including data analysis, PDEs, and linear algebra. The Jordan decomposition reveals the eigenstructure of a matrix as follows. Consider a matrix $`𝒜`$ that has $`r`$ independent eigenvectors with eigenvalues $`\lambda _i`$ of multiplicity $`\rho _i`$. The Jordan decomposition of $`𝒜`$ contains $`r`$ upper triangular “blocks,” as revealed by the diagonalization:

$$^1𝒜=\left[\begin{array}{cccc}𝒥_1& & & \\ & 𝒥_2& & \\ & & & \\ & & & 𝒥_r\end{array}\right],𝒥_i=\left[\begin{array}{cccc}\lambda _i& 1& & \\ & \lambda _i& 1& \\ & & & 1\\ & & & \lambda _i\end{array}\right]$$

The typical approach to computing the Jordan form leads to a numerically unstable algorithm ; taking extra care usually requires more work than the original computation! Recently, however, an experimental algorithmics approach, inferring multiplicity from a geometric analysis of eigenvalue perturbations, has proved successful . It is well known that the computed eigenvalues corresponding to the actual value $`\lambda _i`$ are given by:

$$\lambda _i+|\delta |^{\frac{1}{\rho _i}}e^{\frac{i\varphi }{\rho _i}},$$
(2)
where $`\lambda _i`$ is of multiplicity $`\rho _i`$, and the phase $`\varphi `$ of the perturbation $`\delta `$ ranges over {$`2\pi ,4\pi ,\mathrm{},2\rho _i\pi `$} if $`\delta `$ is positive and over {$`3\pi ,5\pi ,\mathrm{},2(\rho _i+1)\pi `$} if $`\delta `$ is negative. The insight of  is to graphically superimpose numerous such perturbed calculations so that the aggregate picture reveals eigenvalue multiplicity. The phase variations imply that computed eigenvalues lie on the vertices of a regular polygon with $`2\rho _i`$ sides, centered on $`\lambda _i`$, and with diameter influenced by $`|\delta |`$. For example, Fig. 4 shows perturbations for the 8-by-8 Brunet matrix with Jordan structure $`(1)^1(2)^1(7)^3(7)^3`$ , for $`\delta [2^{50},2^{40}]`$. The six “sticks” around the eigenvalue at 7 clearly reveal that its Jordan block is of size 3.222The multiplicity of the second eigenvalue at 7 is revealed at a smaller perturbation level. The “noise” in Fig. 4 is a consequence of having two Jordan blocks with the same eigenvalue and size, and a “ring” phenomenon studied in ; we do not attempt to capture these effects in this paper.

Tab. 3 describes qualitative correspondence analysis of Jordan form. Data are collected by randomly perturbing at a specified magnitude $`\delta `$; the analysis determines multiplicity by detecting symmetry correspondence in the samples. The first aggregation level collects the samples for a given $`\delta `$ into triangles. The second aggregation level finds congruent triangles via geometric hashing , and uses congruence to establish analogy among triangle vertices. Correspondence abstracts the analogy into a rotation about a point (the eigenvalue), and evaluates whether each point rotates onto another and whether matches define regular polygons. Ambiguity-directed sampling collects additional random samples as necessary. A third level then compares rotations across different perturbations, re-visiting perturbations or choosing new perturbations in order to disambiguate. Fig. 5 demonstrates this mechanism on the Brunet matrix discussed above. The top part uses a small set of sample points, while the bottom two parts use a larger set and illustrate a good vs. bad correspondence.

We organized data collection into rounds of 6-8 samples each and experimented with three policies on where next to collect data after completing a round: (1) at the same perturbation level, (2) at a higher perturbation level, or (3) at the same perturbation level unless the number of posited models increased (model “hallucination,” as in bottom of Fig. 5). We tested 10 matrices across 4-10 perturbation levels each, as described in . We varied a tolerance parameter for triangle congruence from 0.1 to 0.5 (effectively increasing the number of models posited) and determined the number of rounds needed to determine the Jordan form. Policy 1 required an average of 1 round at a tolerance of 0.1, up to 2.7 rounds at 0.5. Even with a large number of models proposed, additional data quickly weeded out bad models. Policy 2 fared better only for cases where policy 1 was focused on lower perturbation levels, and policy 3 was preferable only for the Brunet-type matrices. In other words, there is no real advantage to moving across perturbation levels! In retrospect, this is not surprising since our Jordan form computation treats multiple perturbations (irresp. of level) as independent estimates of eigenstructure.

## 5 Discussion

Our general mechanism for uncovering and utilizing correspondence has proved successful on challenging problems from experimental algorithmics. The mechanism leverages properties such as locality, continuity, and decomposability, which are exhibited by the applications studied here as well as many physical systems. Decomposability and locality allow us to qualify correspondence in a manner that drives data collection. Contintuity allows correspondence to obtain confidence in a model even with sparse, noisy data: consistent matches among nearby constituents mutually reinforce each other, allowing correspondence abstraction to detect and filter out inconsistent interpretations.

Correspondence can be an appropriate analysis tool for a variety of reasons. In the spectral portrait application, level curves summarize groups of eigenvalue computations, so their high-level correspondence aids characterization of the underlying computation. In the Jordan form application, however, the higher-level entity has no significance as a geometric object in numerical analysis terms, but correspondence is applicable due to the semantics of superposition. This semantics also leads to phenomena such as hallucination (given enough samples, any pattern can be found), requiring a more careful treatment of decomposition.

Our work is similar in spirit to that of  for weather data interpretation, and can be seen as a significant generalization, formalization, and application of techniques studied there for finding correspondence in meteorological data. Similarly, our correspondence framework captures and generalizes the computation required in object recognition, allowing the body of research developed there to be applied to a broader class of applications, such as experimental algorithmics. Compared to traditional manual analyses of graphical representations in experimental algorithmics, the algorithmic nature of our approach yields advantages such as model evaluation and targeted sampling. As with compositional modeling , we advocate targeted use of domain knowledge, and as with qualitative/quantitative model selection (e.g. ), we seek to determine high level models for empirical data. Our focus is on problems requiring particular forms of domain knowledge to overcome sparsity and noise in spatial datasets. A possible direction of future work is to explore if the inclusion-exclusion methodology popular in grid algorithms  is also useful for tracking correspondence.

Our long-term goal is to study data collection policies and their relationships to qualitative model determination. The notion of estimating problem-solving performance by collecting data (and vice versa) is reminiscent of reinforcement learning  and active learning . The decomposable nature of SAL computations promises to (i) support the design of efficient, hierarchical algorithms for model estimation and (ii) provide a deeper understanding of the recurring roles that correspondence plays in spatial data analysis.


text:
# Comment on “The Effect of Vacuum Polarization … in Strongly Magnetized Plasmas” by Ozel (astro-ph/0203449)

## 1 Introduction

There are two polarization modes for photons (with energy $`E`$ much less than the electron cyclotron energy) propagating in a strongly magnetized plasma: The extraordinary mode (X-mode) and the ordinary mode (O-mode). These two modes have very different opacities in the atmosphere of a magnetized neutron star (NS). In a recent paper (Lai & Ho 2002; hereafter LH02), we discussed the effect of resonant mode conversion due to vacuum polarization: A photon propagating outward in the NS atmosphere can convert from one polarization mode into another as it traverses the resonant density, $`\rho _vY_e^1\eta ^2(B/10^{14}G)^2(E/1\mathrm{keV})^2`$ g cm-3, where $`Y_e`$ is the electron fraction, and $`\eta 1`$ is a slowly varying function of the magnetic field $`B`$. This mode conversion is analogous to the Mikheyev-Smirnov-Wolfenstein mechanism for neutrino oscillation and is effective for $`E>\mathrm{a}\mathrm{few}\mathrm{keV}`$. Because the two photon modes have vastly different opacities, the vacuum-induced mode conversion can affect radiative transport and surface emission from strongly magnetized NSs. In Ho & Lai (2002) (hereafter HL02), this effect was explored further in numerical models of magnetar atmospheres (see HL02 for references of earlier works related to strong-field vacuum polarization and its effects on radiative transport).

In a recent paper (Ozel 2002; hereafter Ozel02), Ozel has dismissed the mode conversion effect discussed our papers as “mistakes”. The purpose of this note is to clarify why our arguments/calculations of this effect are correct. Indeed, there are a number of problems in Ozel02, and many of its criticism and description of our work are incorrect or inaccurate. We restrict our comments to problems in Ozel02 that are directly related to our work111We have previously communicated to Ozel in private our comments described in this note, but the comments were rejected by her. We therefore feel it is necessary to put this note on astro-ph..

## 2 Problems with Ozel (2002)

1. Throughout the paper (especially in the appendix), Ozel02 dismisses the MSW-like mode conversion effect first discussed in our papers. The arguments in Ozel02 are incorrect. Here is how to understand the effect. Adiabatic evolution of a state (whether it is quantum mechanical wavefunction or polarization state of a classical EM wave) means that the evolution is continuous, with no sudden change of the parameters that characterize the state (wavefunction or eigenvalue). Thus the different ways of describing the modes are of real significance. The mode eigenfunctions are characterized by the polarization ellipticity $`K`$: Eq. (2.27) of HL02 gives $`K_j=\beta \left[1+(1)^j(1+r/\beta ^2)^{1/2}\right]`$ for the X-mode and O-mode ($`j=1,2`$), and eq. (2.43) gives $`K_\pm =\beta \pm (\beta ^2+r)^{1/2}`$ for the plus-mode and minus-mode. Obviously, $`K_1`$ and $`K_2`$ are discontinuous across the vacuum resonance (note that $`\beta =0`$ at vacuum resonance, and $`\beta `$ changes sign across the resonance), while $`K_+`$ and $`K_{}`$ are continuous (see Fig.1 of LH02). Similarly, the indices of refraction (the eigenvalues) $`n_1`$ and $`n_2`$ are discontinuous, while $`n_+`$ and $`n_{}`$ are continuous across the resonance (in quantum mechanics, this is called “avoid crossing”). Thus, under adiabatic condition (i.e., when the density gradient of the medium is sufficiently gentle), the polarization state will evolve along the $`K_+`$ or $`K_{}`$ curve, rather than along the $`K_1`$ or $`K_2`$ curve. Since $`K_+`$ manifests as the O-mode at high density and as X-mode at low density (see Fig.1 of LH02), we have adiabatic mode conversion across the resonance. This conversion is not a matter of semantics, but has true physical significance: For example, without this conversion the mode opacity exhibits a spike at the vacuum resonance (see the left panel of Fig. 3 of HL02) but with conversion the opacity has a plateau-like transition (the right panel of Fig.3). As we show both analytically (see sect. 4 of LH02, where we discuss photon decoupling depth in the case of no mode conversion and the case with mode conversion) and numerically (see HL02), this will lead to a genuine difference in the radiative transfer and the emergent spectra.

Note that the adiabatic condition should not be confused with the condition of “Faraday depolarization”. The adiabatic condition requires the background density to vary slowly (the precise condition is derived in LH02), while Faraday depolarization implies that a modal description of radiative transport is possible (as first discussed by Gnedin & Pavlov 1974 and repeated in sect. 2 of Ozel02; the issue of Faraday depolarization was briefly discussed in sect. 2.4 of HL02, as well as in sect. 2.1 of Ho & Lai 2001). Adiabatic mode conversion exists even in a lossless medium (with no absorption and scattering) for which Faraday depolarization is always satisfied. The MSW analogy is: Neutrinos are essentially lossless in the sun, yet they can convert from one flavor to another.

Also, the adiabatic mode conversion should not be confused with mode switching due to scattering; the latter was treated in Ozel02, in HL02, and in many previous papers (see footnote 1 of LH02 and sect. 2.4 of HL02).

2. In addition to the conceptual problems in Ozel02 as discussed above, there are also several (less important) issues we would like to clarify in response to Ozel02:

(1) When mode conversion is neglected, the photon opacity exhibits a sharp feature at vacuum resonance, and in HL02 we adopted an approximate “equal-grid” numerical scheme to handle this feature. Ozel02 has criticized our “equal-grid” scheme. While it is certainly fine to criticize this, it should be noted that the limitation (including the one discussed in Ozel02) and usefulness of the equal-grid method are fully discussed in HL02 (sect. 5.1; this is also mentioned at the end of sect. 6.3, after we comment on previous works). In fact, the subtlety in resolving the density-dependent, narrow vacuum resonance feature (and in particular the fact that the region very close to the resonance contributes most to the optical depth) was first discussed in our papers \[see sect. 4 of LH02 (in particular footnote 2) and sect. 3, 4, 5.1 of HL02\] — there was no hint that this subtlety was appreciated in previous papers (including Ozel 2001; this is why we discussed in detail these issues in our paper, including comparison with previous works in sect. 6.3).

(2) We have emphasized several times in our paper \[e.g., see HL02, sect. 1 after eq. (1.2), the end of sect. 2.4, the second to last paragraph of sect. 7\] that a more rigorous treatment of the problem of radiative transfer in strong magnetic fields with vacuum polarization requires one to go beyond the modal description (as used in all NS atmosphere works so far) and to solve the transport equations for the Stokes parameters (i.e., including the “off-diagonal” terms of the photon density matrix).

## 3 Conclusion

We stand by the main results of LH02 and HL02, although one should keep in mind the limitations (as discussed explicitly in sect. 7 of HL02). Strong-field vacuum polarization introduces several novel features in the radiative transfer whose full, rigorous solution remains an important problem. Much work remains to be done in the area of highly magnetized NS atmospheres/surfaces in order to confront current and future observational data.

We thank Lars Hernquist and Ira Wasserman for useful comments. This work is supported in part by NASA grant NAG 5-8484 and an Alfred P. Sloan fellowship to D.L.