text:
# A Unified Description of the Timing Features of Accreting X-ray Binaries

## 1. INTRODUCTION

Since the first detailed timing studies of accreting X-ray binaries made with the EXOSAT and Ginga satellites, the aperiodic- and quasi-periodic variability and spectral properties of low-magnetic-field ($`B10^{10}`$ G) neutron stars (hereafter simply neutron stars) and of accreting black-hole candidates (hereafter BHC) have been considered in a rather different fashion, despite their striking similarities. In the case of neutron stars, the sources were soon classified into two separate classes (the “Z” and the “atoll” sources) based on their spectral and timing behavior (see van der Klis 1995a and references therein). This classification applied, with few exceptions, to all bright systems; the fainter neutron stars, although not exhibiting the full complement of phenomenology, seemed to fit in as well. In the case of the BHCs, things appeared to be more complex and from the beginning a classification as precise as in the previous case appeared to be problematic (see Tanaka & Lewin 1995; van der Klis 1995a). Despite this complexity, however, some common ground was found between the phenomenologies of neutron stars and BHCs (van der Klis 1994a, 1994b). In particular, strong similarities were pointed out in the properties and dependence on X-ray luminosity and spectral shape of the powerful broad band noise seen in both neutron stars and BHCs, as well as between some of the QPOs. These similarities were interpreted in terms of particular source states occurring in both neutron stars and BHCs depending, presumably, on mass accretion rate.

With the advent of the Rossi X-ray Timing Explorer (RXTE), our knowledge of the properties of the aperiodic variability of accreting X-ray binaries took a substantial step forward. The power-spectral densities of neutron stars revealed two previously unknown QPOs with frequencies in the kHz range which, in addition to broad continuum components, led to power spectra containing up to five simultaneous quasi-periodic oscillations (QPOs). A number of characteristic frequencies could thus be identified for each power spectrum. The kHz oscillations were soon found to behave in a rather regular way and their phenomenological properties were fitted into the tight scheme organizing the neutron-star systems (see van der Klis 2000 for a review). BHCs showed fewer QPOs, superimposed on strong continuum components, which provided a more limited number of characteristic frequencies. Even though some new QPOs were discovered, such as the “hecto-hertz” QPOs (Morgan, Remillard & Greiner 1997; Remillard et al. 1999a, 1999b; Homan et al. 2001) and some $`<`$100 Hz QPOs turned out to be clearly similar to those in neutron stars, BHC phenomenology proved once again more difficult to classify. There are even indications that the similarities between neutron stars and BHCs disappear in the newly accessible frequency range above $``$100 Hz (see, e.g., Sunyaev & Revnivtsev 2000).

Recently, two studies have revisited the above question and identified a number of power-spectral components with frequencies that follow the same correlations in both neutron stars and BHCs. Wijnands & van der Klis (1999a, hereafter WK99) considered the break frequencies of the broad-band noise and the centroid frequencies of the low-frequency QPOs in the 0.02–70 Hz range and found these frequencies to follow a similar correlation in neutron stars and BHCs covering 2.5 decades in frequency. Psaltis, Belloni, & van der Klis (1999, hereafter PBK99) systematically studied the numerous QPOs and peaked noise components in the 0.1–1200 Hz range and identified two features whose characteristic frequencies follow a correlation covering three orders of magnitude in frequency. For bright neutron star sources, these two features are the “horizontal branch oscillation” (HBO) and the lower kHz QPO, respectively. For BHCs and low-luminosity neutron stars, they are the the low-frequency QPO and a broad peaked component at frequencies between 1 and 10 Hz.

These studies suggest that the same types of variability occur in both neutron-star and black-hole sources, which can severely constrain theoretical models for their variability properties. Two major classes of models exist to explain the variable-frequency QPOs in neutron stars: the beat frequency models (e.g., Alpar & Shaham 1985; Strohmayer et al. 1996; Miller, Lamb, & Psaltis 1998), in which the spin-frequency of the neutron star provides one of the characteristic frequencies, and the relativistic precession models (e.g., Stella, Vietri, & Morsink 1999; Psaltis & Norman 2001), in which all QPO frequencies arise from general relativistic frequencies in the accretion disk. The latter class of models do not depend explicitly on the properties of the compact object and are therefore also applicable to the case of the variable-frequency QPOs observed in BHCs. At the same time, discoseismic models offer a promising explanation for the “hecto-Hertz” QPOs in these sources (e.g., Nowak et al. 1997; Perez et al. 1997; see Wagoner 1999 for a review).

Similarly to the various models for the quasi-periodic variability properties of accreting compact objects, various attempts have been made for understanding the physical origin of the broadband aperiodic variability in both types of sources, which in fact constitutes most of the signal in the power spectra of BHCs. This noise continuum can often roughly be described by an $`1/f`$ noise with a characteristic low-frequency break and may, therefore, be produced by a large variety of physical mechanisms. Different models have considered the contribution of time-variable inhomogeneities (“blobs”) in the accretion disks (e.g., Bao and Ostgaard 1995; Takeuchi, Mineshige, & Negoro 1995), short-lived magnetic flares in the accretion-disk coronae (Poutanen & Fabian 1999), erratic injections of soft photons in a Comptonizing medium (Kazanas & Hua 1999), MHD instabilities in the accretion disks (Hawley & Krolik 2001), etc.

From the early days of X-ray astronomy it was clear that the power spectrum of the BHC Cyg X-1 is not a pure power law (Terrell 1972; Nolan et al. 1981). With EXOSAT it became evident that the power spectra of BHCs (Belloni & Hasinger 1990a) as well as atoll sources (Hasinger & van der Klis 1989) are not smooth, but show “bumps and wiggles”. These have been variously described by allowing breaks or exponential cutoffs in power-law fits, by adding broad Lorentzian components, and combinations of these approaches (see references above and, e.g., Berger and van der Klis 1998). For the low-luminosity neutron stars, recent work by Olive et. al (1998) and Barret et al. (2000) revealed power spectra that are quite similar to those of BHCs and atoll sources in the low hard state. All these power spectra are flat-topped, become steeper towards higher frequency and show bumps and wiggles. An approach to the description of these power spectra of neutron stars and BHCs adopted by a number of authors (e.g., Miyamoto et al. 1991; Olive et al. 1998; Nowak 2000) includes no power-law or broken power-law components in the fits (except perhaps a single steep power law at the very lowest frequencies), but uses only a sum of Lorentzian components, some of which are broad (with coherence factor Q$``$1). These authors demonstrated that the power spectra can be well modeled without the need for power-law components; the bumps and wiggles correspond to “shoulders” fitted with individual Lorentzian components.

A limitation of the traditional approaches is that peaked components (QPO) and broad-band components (noise) are treated a priori as different, and therefore are fitted with intrinsically different models. This applies to both BHC and NS systems. In this paper, we expand and elaborate on the description of the power spectra of accreting compact objects as the sum of a small number of distinct Lorentzian components, proposing it as a model to fit all timing features, independent of their Q values. Here, we apply it to power density spectra of a sample of sources in specific states, i.e. the Low/hard state of BHCs and the low-luminosity state of X-ray bursters, but we plan to investigate its applicability on a wider range of sources and states. We seek to offer a unified paradigm for the fitting and interpretation of the power-density spectra of both accreting neutron stars and black holes. This will help to unambiguously identify similar variability components in neutron-star and black-hole sources and thereby to constrain models for quasi-periodic and aperiodic variability of accreting compact objects. We apply this method to six additional sources: the low-luminosity X-ray bursters 1E 1724$``$3045, SLX 1735-269, and GS 1826-24 (some of these data have already been presented by Olive et al. 1998 and Barret et al. 2000), the recently discovered high-latitude X-ray transient XTE J1118$`+`$480, which is most probably a BHC (McClintock et al. 2001; Wagner et al. 2001; for the X-ray data see also Revnivtsev, Sunyaev & Borozdin 2000), the bright and peculiar NS system Cir X-1 (see Shirey et al. 1996), and the high-luminosity Z-source GX 17+2 (see Homan et al. 2002). Finally, we discuss the current status of the correlations between timing features of accreting X-ray binaries in the light of our new approach.

## 2. A LORENTZIAN DECOMPOSITION OF THE OBSERVED POWER SPECTRA

Following Nowak (2000), we describe the power spectra of both neutron-star and black-hole sources by a sum of Lorentzian components, given individually by

$$P(\nu )=\frac{r^2\mathrm{\Delta }}{\pi }\frac{1}{\mathrm{\Delta }^2+(\nu \nu _0)^2},$$
(1)
where $`r`$ is the integrated fractional rms (over $`\mathrm{}`$ to $`+\mathrm{}`$) of each Lorentzian and $`\mathrm{\Delta }`$ its HWHM. We also display them as $`\nu P_\nu `$ vs. $`\nu `$ (see Belloni et al. 1997; Nowak 2000), where $`P_\nu `$ is the power density in rms normalization (e.g., van der Klis 1995b) at frequency $`\nu `$ and quote as their “characteristic frequency” the frequency $`\nu _{max}`$, at which they

attain their maximum in $`\nu P_\nu `$, i.e.,

$$\nu _{max}=\sqrt{\nu _0^2+\mathrm{\Delta }^2}=\nu _0\sqrt{1+1/(4Q^2)}$$
(2)
(see Belloni et al. 1997), where $`Q\nu _0/2\mathrm{\Delta }`$. This is the frequency around which the component contributes most of its power per logarithmic frequency interval and is a measure for the highest frequency covered by the Lorentzian. Note that the maximum in $`\nu P_\nu `$ does not correspond to the centroid frequency $`\nu _0`$ of the Lorentzian but $`\nu _{max}\nu _0`$. The difference is small for narrow QPOs, but becomes large in the case of broad features.

In Figure 1 we show a Lorentzian with $`\nu _0`$=0.1 Hz and different values of the coherence $`Q`$. The left panel shows the usual $`P_\nu `$ vs. $`\nu `$ representation whereas the right panel shows $`\nu P_\nu `$ vs. $`\nu `$. The shift in $`\nu _{max}`$ as Q decreases is evident. As the half-width $`\mathrm{\Delta }`$ increases at constant $`\nu _0`$ in order to produce a broader feature, the value of $`\mathrm{\Delta }`$ progressively takes over from $`\nu _0`$ in determining the characteristic frequency $`\nu _{max}`$. The logarithmic FWHM in $`\nu P_\nu `$ of the feature asymptotically approaches an upper limit of $`\mathrm{log}\frac{2+\sqrt{3}}{2\sqrt{3}}`$. For very low values of $`Q`$, the function approaches a zero-centered Lorentzian ($`Q`$=0), which in $`\nu P_\nu `$ peaks at $`\mathrm{\Delta }`$. For comparison, the dashed line in Fig. 1 corresponds to a zero-centered Lorentzian with $`\mathrm{\Delta }`$=0.1 Hz.

For fitting very broad features, even those with a centroid frequency that is formally larger than zero, a zero-centered Lorentzian model is often adequate. As a matter of fact, this is what a fitting procedure does: if a power spectral component is narrow ($`Q`$1), such a routine converges to a value of $`\nu _0`$ corresponding to the centroid of the observed peak, setting $`\mathrm{\Delta }`$ equal to its half-width. If, on the other hand, the component is wide ($`Q<`$1), the routine chooses a value of $`\nu _0`$ approaching 0 and set $`\mathrm{\Delta }`$ to fit the characteristic frequency of the feature (see Fig. 1). Of course, if the component is too broad, one never obtains a good fit (this would be the case if the power spectra were pure power laws). We discuss below how to deal with the intermediate cases.

Traditionally, the term QPO has been restricted to features with $`Q>`$2 (see van der Klis 1995a), while less coherent components have been classified as ‘noise’. Broad peaked components are sometimes called ‘peaked noise’. As noted by van der Klis (1995b) this classification is arbitrary; there is a continuous spectrum of intermediate cases, and indeed, there is evidence that a component can sometimes change from being a clear band-limited noise component into a clear QPO or vice versa (Di Salvo et al. 2000). Noise components can represent characteristic variability frequencies just as QPOs do. However, there is an ambiguity in how characteristic frequencies are usually reported for narrow (i.e., QPOs) and broad (i.e., noise) variability components. For QPOs, two parameters are usually quoted, i.e., the centroid frequency $`\nu _0`$ and the FWHM in $`P_\nu `$, but it is $`\nu _0`$ which is usually called the “characteristic frequency” of the component. For noise components, only one parameter is usually quoted, for example one describing the width of the component, such as the HWHM $`\mathrm{\Delta }`$, or some power-law break frequency $`\nu _{break}`$. (Simulations show that fitting zero-centered Lorentzians with broken power laws and vice-versa yields similar values for the $`\mathrm{\Delta }`$ of the Lorentzian and the $`\nu _{break}`$ of the broken power-law, but with systematic differences, see Section 4.2.) So, in intermediate cases, the definition of the characteristic frequency is ambiguous; it depends on whether the component is considered to be a QPO or noise.

Based on the discussion above, we adopt the following algorithm for describing the power spectra. (1) We fit a model that is the sum of $`n`$ Lorentzians and involves no power-law or broken power-law components. (2) As the characteristic frequency of each component we adopt its peak frequency in the $`\nu P_\nu `$ spectrum, i.e., $`\nu _{max}\sqrt{\nu _0^2+\mathrm{\Delta }^2}`$. As noted, if a component is narrow, its characteristic frequency $`\nu _{\mathrm{max}}`$ corresponds to the central frequency $`\nu _0`$ of the Lorentzian (see eq. ), whereas, if it is broad, the characteristic frequency approaches $`\mathrm{\Delta }`$. For very low-Q components it makes little difference to either the derived characteristic frequency or the fit to the data to set $`\nu _00`$.

Of course, in the standard “lifetime broadening” interpretation of a Lorentzian as the power spectrum of an exponentially damped harmonic oscillator, the Lorentzian centroid and width correspond to two different physical time scales (respectively, the oscillation period and damping time), whereas in the corresponding interpretation of a zero-centered Lorentzian the width is a measure for the exponential decay time of an exponential shot. Our description of both these types of components in terms of a single characteristic frequency $`\nu _{max}`$ is motivated primarily by empirical considerations and has no immediate interpretation within the damped oscillator picture. We note, however, that in interpretations where broad power spectral components arise by the superposition of variations covering a range of frequencies rather than by lifetime broadening, $`\nu _{max}`$ is just a measure for the highest frequency represented in that range, or the one which contributes the highest power. Identifying, in an accretion disk model, frequencies with radii in the disk would lead to a link between $`\nu _{max}`$ and the inner radius of the disk annulus contributing to the component. More insight on the connection to a physical model might come from the comparison of the fractional rms values of the different components in addition to their characteristic frequencies. As discussed in Sect. 5, two important observational facts are that the relative contribution of the broad components in low-frequency systems is such as to produce an approximate power law distribution, and that the highest frequency component in BH systems is considerably suppressed with respect to NS systems.

Our choice of components and parameters poses the problem of comparison with previous results. In this paper, we apply conversions from $`\nu _0`$ and $`\mathrm{\Delta }`$ into $`\nu _{max}`$ only for the data from Nowak (2000). All other data discussed when dealing with global correlations will have to be examined for complete consistency, but since most of them involve kHz QPOs, which have a rather high Q, it follows that $`\nu _{max}\nu _0`$, so the effects are small. The Lorentzian components of the two power spectra in Nowak (2000) with reasonably constrained parameters have $`Q`$ values ranging from 0.6 to 2.0, corresponding to correction for their characteristic frequencies of 1.9 to 1.03 respectively (see eq. ). As will be discussed below, we believe that fitting all components with Lorentzians will eventually facilitate comparison across different types of sources, as it allows us to treat broad noise components and QPOs in formally the same way without having to decide in advance how to interpret each component. This makes it easier to compare components in different sources and allows us to make the connection with the fits to power spectra of neutron-star sources with much higher characteristic frequencies.

## 3. OBSERVATIONS

In order to test our approach on power spectra containing complex broad-band components, we analyzed RXTE/PCA observations of one black-hole candidate and five neutron-star systems, covering a wide range of characteristic frequencies. The black-hole candidate is the recently discovered X-ray transient XTE J1118$`+`$480 (Remillard et al. 2000, Revnivtsev, Sunyaev & Borozdin 2000). This source was in a low hard state during these observations The neutron-star systems are three low-luminosity X-ray bursters (see Barret et al. 2000): 1E 1724$``$3045 (see also Olive et al. 1998), SLX 1735$``$269, and GS 1826$``$24, the peculiar system Cir X-1 (see Shirey et al. 1996), and the bright Z-source GX 17+2 (see Homan et al. 2002). The three bursters are usually considered to be atoll sources which are always in the low (island) state. The observation log is reported in Table 1.

For each observation, usually consisting of several continuous intervals corresponding to separate RXTE orbits, we produced a series of power spectra from data segments lasting 64 or 128 seconds (depending on the total length of the observation) and a time binning of 1/1024 seconds, normalized after Leahy et al. (1983). These spectra were then averaged together and logarithmically rebinned. The contribution of the Poissonian statistics was removed (see van der Klis 1995a) and the spectra were converted to squared fractional rms (Belloni & Hasinger 1990b; van der Klis 1995b). The average power spectra from different observations were compared and, if they were compatible, they were further averaged together. The total resulting numbers of power spectra averaged can be found in Table 1. Notice that our dataset B for XTE 1118$`+`$480 is the same as that analyzed by Revnivtsev, Sunyaev, & Borozdin (2000), and that our dataset C for 1E 1724$``$3045 corresponds to the data presented by Olive et al. (1998). Also, the power spectra from datasets M (SLX 1735$``$269) and N (GS 1826$``$24) have been presented in Barret et al. (2000). To each of our final average power spectra we fitted a combination of Lorentzians as described in the previous section. For the fits we used XSPEC V11.0.19. The only exception to the procedure outlined above is the

dataset for GX 17+2, for which we used one of the power density spectra from Homan et al. (2002) for the fits.

In practice, it turned out that all power spectra (with the exception of that for GX 17+2, see below) could be fitted with the following components (see Fig. 2, left panel):

* A zero-centered low-frequency Lorentzian L<sub>b</sub> fitting the low-frequency end (the flat top) of the band-limited noise visible in all power spectra, with characteristic (“break”) frequency $`\nu _b`$.
* Two zero-centered Lorentzians L and L<sub>u</sub> fitting the high-frequency end of the band-limited noise, with (“lower” and “upper”) characteristic frequencies $`\nu _{\mathrm{}}`$ and $`\nu _u`$.
* One or two Lorentzians fitting the region around the low-frequency QPO. The profile of this QPO is sometimes more complex than a simple Lorentzian, consisting of a relatively narrow core L<sub>LF</sub> (the “LF QPO”) at characteristic frequency $`\nu _{LF}`$, and a broader “hump” component L<sub>h</sub> with characteristic frequency $`\nu _h`$.

It is evident that, in this description, the band-limited noise is fitted with three broad (and in our case always zero-centered) Lorentzians. For the lowest-frequency Lorentzian, a flat top is actually seen in the power spectra, while the true behavior at low frequencies of the other two components is masked by components that are more powerful there. The lowest-frequency Lorentzian, L<sub>b</sub>, produces the break at $`\nu _b`$ in the band-limited noise, the highest-frequency one, L<sub>u</sub>, produces the second break at $`\nu _u`$, and the middle one, L, essentially “fills the hole” between the two, so that between $`\nu _b`$ and $`\nu _u`$ a relatively smooth, approximately $`1/f`$, spectrum is created, with some “bumps and wiggles”, as described by, e.g., Hasinger and van der Klis (1989) and Belloni and Hasinger (1990a), and on which the LF QPO is superimposed. We found that always $`\nu _b<\nu _{LF}\mathrm{}<\nu _h<\nu _{\mathrm{}}<\nu _u`$; when both the hump L<sub>h</sub> and the LF QPO L<sub>LF</sub> are present, their centroid frequencies $`\nu _h^0`$ and $`\nu _{LF}^0`$ typically coincide, but as the hump is broader, usually $`\nu _{LF}<\nu _h`$. The three components together describing the noise, L<sub>b</sub>, L and L<sub>u</sub>, were always required by the data, although $`\nu _u`$ was not in all cases well-constrained. Of the two components used to describe the low-frequency QPO region, L<sub>LF</sub> and L<sub>h</sub>, at least one was always required. Notice that the frequencies for L<sub>LF</sub> and L<sub>h</sub> presented in Tables 2 and 3 are centroid frequencies and HWHM, which are the values obtained by fitting the data: for these two components the characteristic frequencies discussed in the text have subsequently been computed by applying Equation 2.

For the case of GX 17+2, which has much higher characteristic frequencies, the main differences with respect to the scheme described above are: (a) the L<sub>b</sub> component is very broad but not zero-centered; (b) the L<sub>LF</sub> component shows two additional peaks corresponding to a sub-harmonic and a first overtone (see Homan et al. 2002); (c) the two higher-frequency components L and L<sub>u</sub> are in this case narrow kHz QPO peaks, as expected; (d) two more broad components are needed at frequencies below $`\nu _u`$.

### 3.1. XTE J1118$`+`$480

The power spectra of XTE J1118+480 was fitted with four components (see Fig. 3), the LF QPO already reported by Revnivtsev, Sunyaev & Borozdin (2000), and the three zero-centered Lorentzians L<sub>b</sub>, L, and L<sub>u</sub> for the band-limited noise. The best-fit parameters are reported in Table 2. All components are significant. The values of $`\chi ^2`$ are high, but an examination of the residuals shows no systematic trend. The $`Q`$ value for the LF QPO is around 7.

### 3.2. 1E 1724$``$3045

To fit the ten power spectra of 1E 1724$``$3045, five Lorentzian components are necessary. In the first three power spectra (C,D,E), the component in the LF QPO region is rather broad, with a $`Q`$ of 0.7–0.8 (the same component can be seen in Olive et al. 1998). In the remaining power spectra, an additional narrow peak appears, with $`Q>`$1, sometimes as high as $`Q15`$. By simple comparison of the spectra, we identify the narrow component as the LF QPO L<sub>LF</sub>, and the broad component as the “hump” L<sub>h</sub>. In Figure 4 two representative power spectra are shown: in the power spectrum of observation H the narrow LF QPO component L<sub>LF</sub> is clearly visible while in D it is absent. Note that in the $`\nu P_\nu `$ representation the hump peaks at a higher frequency than the corresponding

narrow LF QPO. The centroid frequencies of these two components are the same, so in a $`P_\nu `$ plot their peaks coincide. For observation C the narrow LF QPO is marginally present, but not required by the fit, while for observation E, like D, it is not detected (see Table 3).

### 3.3. SLX 1735$``$269

For SLX 1735-269, we produced only one power spectrum, using the same data as Barret et al. (2000). There are additional datasets in the RXTE public archive, but their total exposure time is not sufficient for the generation of good power spectra. The power spectrum is shown in Figure 5. A four-Lorentzian model was used to fit the data. The best fit parameters are reported in Table 3. As noted by Barret et al. (2000), they are remarkably similar to those of some observations of 1E 1724$``$3045. Since there is no narrow component, we identify the second Lorentzian as the “hump” component L<sub>h</sub>.

### 3.4. GS 1826$``$24

We obtained five power spectra of GS 1826$``$24. The first one (observation N) corresponds to that presented by Barret et al. (2000). The two X-ray bursts present in the data of observation N (Barret et al. 2000) as well as four more X-ray bursts present in the remaining datasets were excluded from the analysis. The “hump” component was present in all spectra; in four of the power spectra a narrow LF QPO was also observed. Notice that in this source there are significant differences between the centroid frequencies of L<sub>h</sub> and L<sub>LF</sub>. In observations Q and R only

a lower limit could be obtained on the characteristic “upper” frequency $`\nu _u`$. Figure 6 shows the two power spectra of GS 1826$``$24 with the highest and lowest characteristic frequencies. Notice the narrow LF QPO component present in the power spectrum of observation N and absent from observation Q. The best fit parameters can be found in Table 3.

### 3.5. Cir X-1

A detailed description of the power density spectra of Cir X-1 can be found in Shirey et al. (1996, 1998). The L<sub>LF</sub> component in this source was seen to vary from 1 to $``$12 Hz (see also PBK99). We selected two observations among those reported by Shirey et al. (1996), one at relatively low ($``$1.3 Hz) and one at relatively high ($``$8.5 Hz) characteristic frequencies. In order to avoid effects due to the variability of the source, for observation U we selected only the last RXTE orbit, corresponding to the lowest characteristic frequencies. Both power density specta can be fitted with three Lorentzian components (see Table 4 and Fig. 7)): a zero-centered one (L<sub>b</sub>), a narrow one (L<sub>LF</sub>), and a broad one at high frequencies (L).

### 3.6. GX 17+2

In order to include one example of a NS source with high characteristic frequencies (and kHz QPOs), we selected one observation of the Z-source GX 17+2 from the sample analyzed by Homan et al. (2002), corresponding to $`S_z`$ between 0.5 and 0.6 (middle of the horizontal branch, see Homan et al. 2002). Notice that, since the two kHz peaks in this observation are quite well separated from the other components, our fits below 200 Hz do not affect the kHz QPOs, so that for these high-frequency peaks we could adopt the values reported in Homan et al. (2002). We obtained the best fit with a model consisting of six Lorentzian components (see Fig. 8): a broad ($`Q`$=0.6) component which we identify with L<sub>b</sub>, a narrow one with a sub-harmonic and a first harmonic (see Fig. 8) identified with L<sub>LF</sub>, and two additional broad ($`Q<1`$) components with characteristic frequencies below $`\nu _b`$, which we indicate as L$`{}_{}{}^{}{}_{b}{}^{}`$ and L$`{}_{}{}^{\prime \prime }{}_{b}{}^{}`$. Notice that the L$`{}_{}{}^{}{}_{b}{}^{}`$ component is not very significant (about 3$`\sigma `$). The best fit parameters can

be found in Table 4.

### 3.7. Summary of the fits and comparison with power-law models

The three zero-centered Lorentzians L<sub>b</sub>, L and L<sub>u</sub> are required in many spectra to fit the band-limited noise. A narrow low-frequency QPO L<sub>LF</sub> is detected in the black hole candidate, and in some of the neutron star spectra (preferentially those with the lower characteristic frequencies – note that the frequencies are low in the BHC as well). The broad hump component is required in most neutron-star spectra, but not in the BHC. In 1E 1724$``$3045, when L<sub>LF</sub> and L<sub>h</sub> are simultaneously present, their centroid

frequencies are identical, but in GS 1826$``$24 these frequencies are usually different.

In order to compare the quality of our current fits with that of other models for broad components, we also fitted the power spectra first with a combination of a broken power law and one Lorentzian for the narrow LF QPO peak. This model has 7 free parameters (compared with the 9 of our model). We applied this model to two representative power spectra (observations C of 1E 1724$``$3045 and A of XTE J1118$`+`$480). The best fit $`\chi ^2`$ values were 1220 (/75 dof) and 478 (/75 dof) respectively, to be compared with our values 166(/73 dof) and 156(/73 dof). Clearly, the multi-Lorentzian model fits much better. A model consisting of a twice-broken power law (Belloni and Hasinger 1990a) plus a LF QPO, with 9 free parameters, gives similar $`\chi ^2`$ values as our multi-Lorentzian fit. As noted in §2, the advantage of our current approach is that it facilitates comparison across source types.

## 4. GLOBAL CORRELATIONS

As PBK99 and Nowak (2000) noted, the correlations observed between the frequencies of the LF QPO and both kHz QPOs in luminous neutron stars can be extended to lower frequencies by relations between the LF QPO and two broad Lorentzians respectively taking over the role of the lower kHz QPO (seen in several low-luminosity neutron stars and in BHCs) and the upper kHz QPO (seen in Cyg X-1 and GX 339$``$4). Our results may fit in with this. The frequencies $`\nu _b`$, $`\nu _h`$ or $`\nu _{LF}`$, $`\nu _{\mathrm{}}`$, and $`\nu _u`$ define a model that is formally similar to that for luminous neutron stars. In Figure 2 the left panel, the four frequencies are those of our fit to the BHC XTE 1118$`+`$480. In the right panel, they represent a power spectrum of a luminous neutron-star Z source (e.g., van der Klis et al. 1997); here the high-frequency components have been made narrower. The band-limited noise L<sub>b</sub> and the low-frequency QPO $`L_{LF}`$ are both common also at high luminosities (van der Klis 1995a). A broad peaked noise component (“bump” or “wiggle”) similar to our L component is seen in several higher luminosity sources as well (e.g., Shirey et al. 1998; PBK99) and may be a low-frequency manifestation of the lower kHz QPO (PBK99). Finally, L<sub>u</sub> was first found by Nowak (2000) in Cyg X-1 and GX 339$``$4 and tentatively identified as a low-frequency version of the upper kHz QPO. The presence of such a component was already noted previously (but not interpreted in the same way) by Olive et al. (1998) for 1E 1724$``$3045 (our observation C; see also Yoshida et al. 1993).

The power spectra can be more complex than the sum of four Lorentzians (see, e.g., Ford and van der Klis 1998; Di Salvo et al. 2001; van Straaten et al. 2000 and our GX 17+2 power spectrum). In particular, the low-frequency QPO often shows multiple harmonics, the second sometimes being the strongest one (e.g., Belloni et al. 1997). In the following we use the frequencies of the strongest QPO peaks in the power spectrum (but see Sect. 4.2). We try to identify the four components discussed here in all source classes, and investigate if their frequencies follow or extend the correlations of WK99 and PBK99. If these components at very different frequencies are thus phenomenologically related, perhaps they are physically related. In addition to our results, we include additional new points on Cyg X-1 and GX 339$``$4 (Nowak 2000), XTE J1550$``$564 (Homan et al. 2001), 4U 1915$``$05 (Boirin et al. 2000), and 4U 1728$``$34 (Di Salvo et al. 2001).

### 4.1. Correlation between the flat-top power level and the break frequency

In Cyg X-1 the break frequency $`\nu _b`$ and the power level of the flat part of the band-limited noise P<sub>flat</sub> are anti-correlated (Belloni & Hasinger 1990). From their figure, it is evident that P$`{}_{\mathrm{flat}}{}^{}1/\nu _b`$. The high-frequency end of the power spectrum remains approximately unchanged. This is seen in other BHCs (Yoshida et al. 1993; van der Klis 1994a; Méndez & van der Klis 1997), in the neutron star 4U 0614+09 (van Straaten et al. 2000), and also clearly in our two sources showing significant changes in their characteristic frequencies (Fig. 9 and Fig. 10).

The Lorentzian description of the band-limited noise suggests a clue to this effect. From Equation (1), for a zero-centered Lorentzian, $`P_{\mathrm{flat}}=r^2/\pi \mathrm{\Delta }`$, i.e., when the fractional rms $`r`$ of the Lorentzian is constant, $`P_{\mathrm{flat}}`$ is proportional to $`1/\mathrm{\Delta }`$, and hence, as $`\nu _0=0`$, to $`1/\nu _{max}`$. Hence, the observed effect can be reproduced by allowing the characteristic frequency of L<sub>b</sub> to vary, while keeping its total fractional rms constant. The dotted lines in Fig. 10 represent the predicted $`P_{\mathrm{flat}}`$ vs. $`\nu _b`$ relation for constant $`r`$. The definition of $`r`$ used here (Eq. 1) involves integration from

$`\mathrm{}`$ to $`+\mathrm{}`$, so the true values of the rms are a factor $`\sqrt{2}`$ less. Most points approximately follow such correlations, with a difference between sources in normalization (i.e., different $`r`$).

So, the total rms of L<sub>b</sub> is rather constant in each of our sources but different between sources. The neutron stars all cluster near an rms of $``$20%, whereas the two BHCs have values of $``$30% and $``$40%, respectively, consistent with what would be expected if rms correlates with a basic compact-object property such as mass. The lower rms in neutron-star sources could also result from the contribution of a non-varying component from the stellar surface.

### 4.2. Correlation between the low-frequency QPO and break frequency

This correlation (WK99) contains data from both BHCs and neutron stars, to which we add our data except SLX 1735$``$269 (which does not have a LF QPO, see Table 3), the two accurate points from Table 1 in Nowak (2000), GX 339-4 (full dataset) and Cyg X-1 (4-15 keV data), corrected for the low-Q shift (see Sect. 2). and Q-corrected GX 339-4 points from Nowak, Wilms & Dove (2002), for which we plot their Q2 component as our l<sub>h</sub>. For Cyg X-1 we used the first “harmonic” of the peak, which is by

far the strongest (see Nowak 2000 and Sect. 4). We also included data for 4U 1728-34 (Di Salvo et al. 2001), but not for 4U 1915-05 (Boirin et al. 2000) or XTE J1550-564 (Homan et al. 2001), as these authors report parameters for the QPOs only.

The result is shown in Figure 11. The two frequencies follow the same trend for over three orders of magnitude, but the individual correlations in different sources are significantly different. Various systematic effects affect the data. The characteristic frequency of the break depends on energy and hence on source spectrum and energy band (Belloni et al. 1997). WK99 used broken power laws, while we used zero-centered Lorentzians. By simulations, we find that this causes over- or underestimates of the break frequency by $`15`$%, unlikely to explain the largest of the offsets we observe. Some discrepancies are consistent with factors of 2 or 4, suggesting that that some of the QPO frequencies are not fundamental. Indeed, the sub-harmonic in GX 17+2 lies on the WK99 correlation.

Some sources show an additional broad component, L<sub>h</sub>, at a slightly higher characteristic frequency than L<sub>LF</sub> (see Sect. 3), which might represent the non-Lorentzian wings of the narrow LF QPO. Its characteristic frequency $`\nu _h`$ is well correlated with $`\nu _{LF}`$: the corresponding points are added in Figure 11. This correspondence between a narrow QPO and a broader component at a similar centroid frequency may be related to the cases in which a narrow QPO appears near the break of a band-limited noise (e.g., Morgan, Remillard, & Greiner 1997 in GRS 1915+105; Di Salvo et al. 2001 in 4U 1728$``$34; Wijnands & van der Klis

1999b and Revnivtsev, Borozdin & Emelyanov 1999 in XTE J1806–246; Homan et al. 2001 in XTE J1550–564). Notice that this L<sub>h</sub> component is clearly present in the GX 339-4 data from Nowak, Wilms, & Dove (2002): when only one component is present, it is L<sub>h</sub>, when two (or more) are present, they are L<sub>LF</sub> and L<sub>h</sub> (see Fig. 11).

### 4.3. Correlation between low-frequency QPO and L

Figure 12 shows the correlation of PBK99, together with the Nowak (2000), Boirin et al. (2000), Homan et al. (2001) Di Salvo et al. (2001) data. and Nowak, Wilms, & Dove (2002) data. Most points agree quite well with the correlation, except two GS 1826$``$24 points and the Homan et al. (2001) points (see below).

The Nowak (2000) points would be further above the correlation without Q-correction. Most of the GX 339-4 points by Nowak, Wilms & Dove (2002) lie a factor of 2 above the main correlation, again an indication that the frequency adopted for L<sub>LF</sub> might not be that of the fundamental. For the QPO frequencies of XTE J1550-564 (Homan et al. 2001) we used the strongest of the low-frequency harmonic peaks for L<sub>LF</sub> and the 185$``$285 Hz QPOs for L. These points do not follow the main correlation of PBK99, but extend its “second branch”, but points corresponding to the beginning of the outburst (Cui et al. 1999) in the original PBK99 correlation, follow the main branch. The points of 4U 1728$``$34 by Di Salvo et al. (2001) and of 4U 1915$``$05 by Boirin et al. (2000) also follow partly the main branch and partly the second branch.

Di Salvo et al. (2001) concluded that the feature plotted on this second branch vs. $`\nu _{\mathrm{}}`$ is in fact, in our current terminology, $`L_b`$, and hence the main and second branch in PBK99 together represent the WK99 relation, in

accordance with their observation that L<sub>b</sub> changes from band-limited noise into a QPO when its frequency increases. Indeed, if we plot $`\nu _b`$ vs. $`\nu _{\mathrm{}}`$, the points cluster around the extension of the second branch to lower frequencies (not shown in Fig. 12).

### 4.4. Correlation between the frequencies of components L and L<sub>u</sub>

To this correlation (Fig. 13) between upper and lower kHz QPO, we added new data for 4U 1728$``$34 (Di Salvo et al. 2001) and 4U 1915-05 (Boirin et al. 2000), as well as our $`L_u`$ components and those by Nowak (2000). Nowak’s BHC points (with Q-correction) and SLX 1735$``$269 follow the low-frequency extension of the correlation, but the BHC XTE J1118$`+`$480 and the neutron stars 1E 1724$``$3045 and GS 1826$``$24 are above it. Perhaps this discrepancy is related to the presence of other components in the power spectrum in the few 100 Hz range. A broad Lorentzian with an aprozimately constant $`100200`$ Hz frequency occurs in low-luminosity neutron stars (see, e.g., Ford & van der Klis 1998; van Straaten et al. 2000; Di Salvo et al. 2001) and might also occur in these sources (van Straaten et al. 2001).

## 5. DISCUSSION

Applying the model described in Sect. 2, we have shown that the power spectra of low-luminosity burst sources, of a BHC in the low state, of a Z source, and of Cir X-1 can be decomposed into the sum of a small number ($`45`$) of Lorentzian components. We obtain, for each power spectrum, four characteristic frequencies: $`\nu _b`$, $`\nu _{LF}`$ or $`\nu _h`$, $`\nu _{\mathrm{}}`$, and $`\nu _u`$, which include the broad noise components as well as the narrow QPOs identified in previous studies. In this manner, the identification of similar power spectral components between different sources (as done in PBK99) becomes easier and free of potentially subjective interpretation (see, however, Nowak 2000 and §2 for a discussion

of some complications related to the wings and harmonic structure of the low-frequency QPO).

The presence in some power spectra of the additional L<sub>h</sub> component needs to be discussed in some more detail. This component seems to be closely related to L<sub>LF</sub>. Sometimes only one of these two components is observed; sometimes both are present. When they appear simultaneously, their centroid frequencies are approximately the same in 1E 1724$``$3045 but not in GS 1826$``$24. As already mentioned, there are cases of power spectra in the literature (see, e.g., Morgan, Remillard, & Greiner 1997), for which the characteristic frequency of a band-limited noise component is rather close to that of a low-frequency QPO. Moreover, Di Salvo et al. (2001) actually observed the transformation of one into the other, which also involved ’transitory’ power spectra where both a narrow and a broad component seem to be present simultaneously, with the QPO sitting on the break of the noise spectrum. It seems likely that the connection between $`\nu _h`$ and $`\nu _{LF}`$ is a similar case. Based on their frequency coincidence, the two components may be related physically as well. However, unlike Di Salvo et al. (2001), we observe that, in the case discussed here, there is a preference for the presence of the narrow LF QPO when the characteristic frequencies are low, not when they are high. That the narrow LF QPO component tends to disappear in the spectra with the highest characteristic frequencies suggests that caution is necessary when identifying this LF QPO with those seen at higher luminosities. It is possible that, after one narrow Lorentzian disappears, another one comes into prominence, as observed in 4U 1728-34 (di Salvo et al. 2001). A more detailed discussion of the characteristics of this component in a number of BHCs will be presented in a forthcoming paper.

We have also studied the correlations between the properties of the power-spectral components and extended the correlations found earlier for a wide range of other accreting sources. Our results are consistent with those of earlier studies, but also indicate some discrepancies, which may arise from the fact that different sources may follow similar but not identical correlations (see also Psaltis et al. 1998 where this was demonstrated quantitatively for the case of kHz QPOs in bright neutron-star sources). Before discussing in detail the nature of these discrepancies, a more homogeneous analysis using identical models would be necessary for all sources.

Our description of the power spectra in terms of just Lorentzian components characterized by a single characteristic frequency $`\nu _{max}`$ puts noise and QPO phenomena on equal footing. This approach is in accordance with recent findings from the analysis of RXTE data of the atoll source 4U 1728$``$34 by Di Salvo et al. (2001), who showed that in some cases a broad flat-top noise component turns into a narrow, peaked component, therefore becoming a “standard” QPO. The fact that $`\nu _{max}`$ of components of widely differing $`Q`$ seems to maintain (at least on some occasions) a simple relation to other parameters suggests that perhaps the power distribution of these components is related to the extent of disk annuli whose inner radius, setting $`\nu _{max}`$, varies in a simpler fashion than their width (§2). It might also be related to the presence of a lifetime broadening mechanism that operates at a characteristic timescale related to the dynamical frequencies in the disk, such as viscous dissipation (see, e.g., Psaltis & Norman 2001). Finally, in a shot-noise model of the broad-band variability, large changes in the $`Q`$ value of a component can be produced by a change in the correlation time between different shots (see, e.g., Vikhlinin, Churazov & Gilfanov 1994), accounting thus for such changes. However, as noted in Sect. 2, life-time broadening mechanisms do not naturally involve $`\nu _{max}`$ as a basic parameter; rather, centroid and width reflect aspects of the physical process which are a priori independent.

Independent of the particular physical mechanism that is responsible for the various power-spectral components, our results strongly suggest that the low-frequency Lorentzian (or flat-topped power-law component) is not different from the other components. Indeed, not only can narrow QPOs turn into flat-topped components, as discussed above, but also the frequency $`\nu _\mathrm{b}`$ is correlated to the frequencies of all narrow QPOs. The latter property, together with the fact that in BHCs the low-frequency “noise” accounts for almost the entire power in X rays, suggests that this component is generated at the same region in the accretion flow as the other QPOs (see, e.g., Psaltis & Norman 2001).

Our analysis, with the help of the improved statistics in the RXTE/PCA data as compared to previous, less sensitive missions suggests that no $`1/f`$ noise component is needed for the interpretation of power spectra in BHCs and low-luminosity bursters. However, the superposition of Lorentzian components does produce an overall power spectrum that in neutron stars resembles an $`1/f`$ distribution over a range of frequencies, and in the BHC XTE J1118+480 is clearly steeper, but also roughly follows a power law. The relative total rms strengths of the various Lorentzians are such that a relatively smooth power spectrum (in the $`P_\nu `$ vs. $`\nu `$ representation) is produced. This property needs to be addressed by theoretical models as much as the values for the characteristic frequencies do.

Within the context of our results, we can also discuss the recent suggestion by Sunyaev & Revnivtsev (2000) of a possible way of distinguishing between systems hosting a black hole and a weakly-magnetic neutron star based on their power spectra. According to these authors, neutron-star systems can show significant power above $`300`$ Hz, while BHCs do not. The relevant effect can be seen comparing Figure 3 to Figure 4: the power spectrum of 1E 1724-3045 contains much more power at high frequencies than that of XTE 1118+480, as was already noted by Sunyaev & Revnivtsev (2000) and Revnivtsev, Sunyaev & Borozdin (2000) for these two sources respectively.

Our power-spectral decomposition suggests that the difference in the power spectra of BHCs and neutron stars is caused by two effects. The first is that all characteristic frequencies are lower in the case of the BHC, probably reflecting the mass dependence of the dynamical time scale in the accretion flow, as was also noted by Sunyaev & Revnivtsev (2000). The second is more unexpected: in the power spectrum of XTE 1118$`+`$480 the fractional rms of the Lorentzian with the highest frequency is considerably lower with respect to the other components. This effect is not easy to understand. However, the fact that the component is significantly detected in this BHC, albeit weaker than in neutron stars, does seem to suggest that it is not a feature that requires a neutron-star surface for its generation; apparently a component with these characteristics can also be produced by the accretion disk around a black hole.

In conclusion, we suggest that the underlying physical mechanisms producing similar components in all these sources may be the same and that the scatter in their different correlations is due to a combination of differences in the analysis method and differences in the fundamental system parameters. Although not all timing features in black-hole and neutron-star low-mass X-ray binaries are included in this representation (e.g., the VLFN, the normal branch oscillations in Z and perhaps some atoll sources—but see PBK—, the $`1`$ Hz QPO in dipping neutron stars, the $`100`$ Hz QPOs in atoll sources—but see §4.4—, and the HS power-law noise as well as the very-low frequency QPOs in some BHCs), the ones included here are important and their understanding would be a major step forward in the study of variability in accreting X-ray sources.

We thank S. Campana, L. Stella, and C. Dullemond for useful discussions, D. Barret for help with importing power spectra into XSPEC, and J. Homan for having provided the power density spectrum of GX 17+2. T. B. acknowledges the hospitality of Harvard-Smithsonian CfA and of MIT, and thanks the Cariplo Foundation for financial support. D. P. has been supported by a post-doctoral fellowship of the Smithsonian Institution and by the NASA Long Term Space Astrophysics program under grant NAG 5-9184. This work was supported in part by the Netherlands Organization for Scientific Research (NWO).


format: text

text:
# Slices of matrices — a scenario for spectral theory

## The basic definition

Let $`S`$ be a real $`n\times n`$ symmetric matrix with simple spectrum $`\sigma (S)=\{\lambda _1>\mathrm{}>\lambda _n\}.`$ When is a matrix simultaneously an orthogonal and an upper triangular conjugation of $`S`$? More precisely, we consider the *slice* $`_S`$ through $`S`$, defined to be the intersection of the sets $`\{Q^TSQ`$, for arbitrary real orthogonal matrices $`Q`$ with $`detQ=1\}`$ and $`\{RSR^1,`$ for arbitrary real upper triangular matrices $`R`$ with $`\text{diag}R>0\}.`$

How large is a slice? Clearly, if $`Q^TSQ=RSR^1`$ then $`SQR=QRS`$, and hence, since $`S`$ has simple spectrum, $`QR=f(S),`$ for some (real) polynomial $`f`$. The matrices $`Q`$ and $`R`$ are uniquely determined from $`f(S)`$ if $`f(S)`$ is invertible: this is the standard *QR-factorization* of a matrix (). For convenience, we write $`f(S)=QR=\left[f(S)\right]_Q\left[f(S)\right]_R`$.

Slices have been appearing in disguised form in the literature of numerical analysis and integrable systems, and questions about the geometry of slices have come up intermittently. By putting these three aspects side by side, we expect to convince the reader that the concept is indeed a natural one.

## Slices and the computation of eigenvalues

Francis, in his fundamental work on the $`QR`$ algorithm () considered the following map between matrices. Take $`S`$ invertible symmetric, factor $`S=QR`$ and define the *$`QR`$ step* $`S^{}=RQ`$. It is clear that $`S^{}`$ is symmetric with the same spectrum of $`S`$, since $`S^{}=Q^TSQ`$. But more is true: $`S^{}=RSR^1`$. In particular, as Francis had already pointed out, $`S`$ and $`S^{}`$ have the same bandwidth (the reader should compare this argument with the usual one ). Also, $`S^{}`$ belongs to the slice through $`S`$: in the notation for the elements of $`_S`$ presented in the introduction, $`S^{}`$ is associated to the function $`f(x)=x`$.

Numerical analysts also know well that $`S^{(k)}=\left[S^k\right]_Q^TS\left[S^k\right]_Q`$ equals the matrix obtained by applying $`k`$ times the $`QR`$ step starting from $`S`$: in other words, $`S^{(k)}`$ is the $`k`$-th step of the $`QR`$ iteration starting from $`S=S^{(0)}`$. Again, $`S^{(k)}`$ belongs to $`_S`$, and is associated to $`f(x)=x^k`$.

Instead, by taking $`f(x)=x^{1/k}`$, the resulting matrix $`S^{(1/k)}`$ is, in a precise sense, the $`1/k`$-th $`QR`$ step. Taking the limit ()

$$\underset{k\mathrm{}}{lim}(S^{(1/k)}S^{(0)})k=\underset{ϵ0}{lim}\frac{S^{(ϵ)}S}{ϵ}$$

yields a vector field $`\dot{S}=X(S)`$, whose solution $`S(t)`$ starting at $`S(0)=S`$ at time $`k`$ equals $`S^{(k)}`$! The computation of this limit appears in a number of arguments in the subject. The expression for $`S^{(ϵ)}`$ is

$$S^{(ϵ)}=\left[S^ϵ\right]_Q^TS\left[S^ϵ\right]_Q,\text{ where }S^ϵ=\left[S^ϵ\right]_Q\left[S^ϵ\right]_R.$$

Evaluating the derivatives in $`ϵ`$, we learn that

$$\frac{d}{dϵ}S^{(ϵ)}=(\frac{d}{dϵ}\left[S^ϵ\right]_Q^T)S\left[S^ϵ\right]_Q+\left[S^ϵ\right]_Q^TS(\frac{d}{dϵ}\left[S^ϵ\right]_Q),$$

$$(\mathrm{log}S)S^ϵ=(\frac{d}{dϵ}\left[S^ϵ\right]_Q)\left[S^ϵ\right]_R+\left[S^ϵ\right]_Q(\frac{d}{dϵ}\left[S^ϵ\right]_R).$$

From the last equation,

$$\left[S^ϵ\right]_Q^T(\mathrm{log}S)S^ϵ\left[S^ϵ\right]_R^1=\left[S^ϵ\right]_Q^T(\frac{d}{dϵ}\left[S^ϵ\right]_Q)+(\frac{d}{dϵ}\left[S^ϵ\right]_R)\left[S^ϵ\right]_R^1,$$

which obtains

$$\left[S^ϵ\right]_Q^T(\mathrm{log}S)\left[S^ϵ\right]_Q\left[S^ϵ\right]_R\left[S^ϵ\right]_R^1=\mathrm{log}S^{(ϵ)}=\left[S^ϵ\right]_Q^T(\frac{d}{dϵ}\left[S^ϵ\right]_Q)+(\frac{d}{dϵ}\left[S^ϵ\right]_R)\left[S^ϵ\right]_R^1.$$

Now evaluate the derivatives at $`ϵ=0`$: we must have $`I=S^0=I.I=\left[S^0\right]_Q\left[S^0\right]_R`$ and $`S^{(0)}=S`$. The equation above yields

$$\mathrm{log}S=\frac{d}{dϵ}\left[S^ϵ\right]_Q|_{ϵ=0}+\frac{d}{dϵ}\left[S^ϵ\right]_R|_{ϵ=0}.$$

The two terms in the right hand side of the last equation are special matrices: they are respectively skew symmetric and upper triangular. Consider the (unique, linear) decomposition of a matrix $`M=\mathrm{\Pi }_aM+\mathrm{\Pi }_uM`$ as a sum of a skew symmetric and an upper triangular matrix. Then, from the expression for the derivative of $`S^{(ϵ)}`$, $`\frac{d}{dϵ}\left[S^ϵ\right]_Q|_{ϵ=0}=\mathrm{\Pi }_a\mathrm{log}S`$. The vector field which interpolates the $`QR`$ iteration then is $`X(S)=[S,\mathrm{\Pi }_a\mathrm{log}S]`$.

The iterations and flows defined above lie in $`_S`$. In particular, both preserve the eigenvalues of the initial condition, its symmetry and its bandwidth.

For numerical analysts, the asymptotic behavior of the $`QR`$ iteration, well known to Francis, is of capital importance. Say, for example, that $`J`$ is an arbitrary Jacobi matrix (i.e., a real, tridiagonal matrix whose entries $`J_{k,k+1}=J_{k+1,k},k=1,\mathrm{},n1`$ are strictly positive). Starting with $`J`$, the $`QR`$ iteration converges to a diagonal matrix $`D`$. Not only $`D`$ and $`J`$ have the same spectrum but $`D`$ must lie in the closure of $`_J`$, the slice through $`J`$. Steps of $`QR`$ type related to different functional parameters $`f`$ give rise to iterations which, starting from $`J`$, always converge to diagonal matrices, with diagonal entries consisting of the (distinct) eigenvalues of $`J`$ in an arbitrary order. The closure of the slice $`_J`$ clearly ought to contain additional points. In a nutshell, they correspond to *reducible* Jacobi matrices, i.e., matrices for which some entry $`J_{k,k+1}=J_{k+1,k}`$ is zero. This will be explained in the sequel.

## Slices and the Toda flows

We now describe briefly the Toda lattice, which is an integrable system for which slices appear as natural phase spaces. Consider $`n`$ particles on the line with positions $`x_k`$ and velocities $`y_k`$ evolving under the Hamiltonian

$$H=\frac{1}{2}\underset{k=1}{\overset{n}{}}y_k^2+\underset{k=1}{\overset{n1}{}}e^{(x_kx_{k+1})}.$$

This Hamiltonian was introduced as a model of wave propagation in one dimensional crystals (). It was Flaschka’s remarkable discovery () that this dynamical system is equivalent to the matrix differential equation

$$\dot{J}=[J,\mathrm{\Pi }_aJ],$$

where $`J`$ is the Jacobi matrix with nonzero entries

$$J_{k,k}=y_k/2\text{ and }J_{k,k+1}=J_{k+1,k}=e^{(x_kx_{k+1})/2}/2.$$

The fact that the differential equation is in the so called *Lax pair form* () implies that the evolution $`J(t)`$ is actually an orthogonal conjugation of the initial condition $`J(0)`$. But again more is true: the evolution stays within the set of Jacobi matrices, as we should expect, since there are no other physically significant variables in the problem (velocities essentially are the diagonal entries, and distances between particles with neighboring indices correspond to the off-diagonal entries). Indeed, it is not hard to check that

$$J(t)=\left[e^{tJ(0)}\right]_Q^TJ(0)\left[e^{tJ(0)}\right]_Q=\left[e^{tJ(0)}\right]_RJ(0)\left[e^{tJ(0)}\right]_R^1,$$

the celebrated solution of the Toda lattice by factorization (). To check this formula, proceed as in the computation of the limit of the previous section. Thus, again, the orbit $`J(t)`$ lies within the slice through $`J(0)`$.

Moser also computed the asymptotic behavior of the Toda lattice (): for $`f(x)=x`$, the orbit $`J(t)`$ starting from a Jacobi matrix $`J(0)`$ converges to a diagonal matrix, with eigenvalues disposed in decreasing order. This is in accordance with the relationship between Toda and $`QR`$ and has a natural physical interpretation. Diagonal entries are velocities: in the long run, particles move apart and tend to undergo uniform motion, each with speed given by a different eigenvalue. Faster particles move ahead, explaining the orderd outcome of the eigenvalues along the diagonal found in applications of the $`QR`$ iteration. The remarkable fact that asymptotic speeds both at $`\mathrm{}`$ and $`+\mathrm{}`$ are the same will not be relevant to us: this is strong evidence of the integrability of the system — each orbit ’remembers’ this data. Other matrices belonging to the boundary of the slice through $`J(0)`$ correspond to asymptotic behavior in which the system of particles breaks into essentially disconnected components, the so called *clustering*.

For a symmetric matrix $`S`$ and an arbitrary function $`g`$, the matrix equation $`\dot{S}=[S,\mathrm{\Pi }_ag(S)]`$ admits a similar solution by factorization: conjugate the initial condition by $`\left[e^{tg(S(0))}\right]_Q`$. It is this last fact which historically was responsible for relating the $`QR`$ iteration and Toda flows (, ): when $`g(x)=\mathrm{log}(x)`$, the resulting differential equation gives rise to the orbits interpolating the $`QR`$ iteration. For applications of these differential equations to numerical analysis, the reader may consult .

## Geometric aspects of slices

Detailed study of variables for which the Toda flow becomes especially simple led Moser () to a parametrization of Jacobi matrices. It turns out that Jacobi matrices have simple spectrum and its eigenvectors always have nonzero first coordinates — in particular, they can be normalized so as to have (strictly) positive first coordinates. Moser proved that the map taking Jacobi matrices to $`n`$-uples of distinct real numbers (the eigenvalues) and a point in the first octant of the unit sphere in $`^n`$ (the first coordinates of the normalized eigenvectors) is a diffeomorphism. He then showed that an appropriate choice of Toda flow (i.e., of functional parameter $`g`$) gives rise to orbits joining any two Jacobi matrices with the same spectrum. Thus, given a Jacobi matrix $`J`$, $`_J`$ is the set of all Jacobi matrices with the same spectrum as $`J`$. (We remind the reader that Jacobi matrices always have simple spectrum.) From Moser’s parametrization, *Jacobi slices* are diffeomorphic to $`^{n1}`$.

More generally, getting back to the definition of slices, we have seen that a matrix $`S^{}_S`$ is of the form $`S^{}=Q^TSQ=RSR^1`$, where $`f(S)=QR`$ for some function $`f`$. Keeping account of the requests on $`Q`$ and $`R`$, one may obtain a simple coordinatization of a slice. A matrix $`S`$ is *irreducible* if $`S`$ has no invariant subspace generated by a subset of the canonical vectors $`e_1,\mathrm{},e_n`$. Indeed, there is a diffeomorphism between positive polynomials up to scalar multiplication and elements of $`_S`$, for an irreducible symmetric matrix ().

Which matrices are in the boundary $`_J`$ of the slice through a Jacobi matrix $`J`$? Numerical analysts knew that the diagonal matrix $`\lambda =(\lambda _1,\mathrm{},\lambda _n)`$ and its permuted counterparts $`\lambda _\pi =(\lambda _{\pi (1)},\mathrm{},\lambda _{\pi (n)})`$. As we shall see, the simplest possible way of combining such $`n!`$ points to form a reasonable boundary turns out to be the topological description of $`_J`$. Let $`𝒫_\mathrm{\Lambda }`$ be the *permutohedron*, obtained by taking the convex closure of the points of the form $`\lambda _\pi ^n`$. Notice that actually $`𝒫_\mathrm{\Lambda }`$ has an interior of dimension $`n1`$, since the sum of the coordinates of any point in $`𝒫_\mathrm{\Lambda }`$ equals the trace of $`J`$. The permutohedron $`𝒫_\mathrm{\Lambda }`$ is homeomorphic to the closure of $`_J`$, as was proved in  with combinatorial arguments.

The figure below describes the situation for a $`3\times 3`$ Jacobi matrix $`J`$ with spectrum $`\sigma (J)=\{4,2,1\}`$. On the left is a topological representation of the closure of the slice $`_J`$. Interior points of the hexagon are points in $`_J`$. The boundary consists of two kinds of points: the six diagonal matrices, corresponding to the vertices, and the matrices which form the edges, which have either entry $`(1,2)`$ or $`(2,3)`$ equal to zero. On the right, the permutohedron associated to $`\sigma (L)`$ is projected on the $`(x,y)`$ plane. Notice that a map taking vertices to vertices, say, taking $`\text{diag}(1,2,4)`$ to $`(1,2,4)`$ and $`\text{diag}(2,1,4)`$ to $`(2,1,4)`$ may not behave so naïvely on all vertices, for continuity reasons.

The permutohedron associated to a Jacobi $`4\times 4`$ matrices is also drawn below. Faces have a clear meaning: there must be four hexagons associated to matrices for which entry $`(1,2)`$ equals zero (and hence contained a Jacobi $`3\times 3`$ block of fixed spectrum), six quadrilaterals corresponding to matrices with entry $`(2,3)`$ equal to zero, and four more hexagons for matrices with entry $`(3,4)`$ equal to zero. There are also $`4!`$ vertices.

The permutohedron is familiar to spectral theorists, from the well known Schur-Horn theorem (). Let $`S`$ be a symmetric matrix with simple spectrum $`\sigma (S)=\{\lambda _1,\mathrm{},\lambda _n\}`$, and consider the set $`𝒪_S=\{Q^TSQ,Q\text{ orthogonal }\}`$. Then the map $`S\text{diag}S`$ is surjective from $`𝒪_S`$ to $`𝒫_\mathrm{\Lambda }`$. Such a map is highly not injective. Slices are in a sense minimal sets in $`𝒪_S`$ in bijection with permutohedra.

The presence of a convex polytope raised the possibility of relating slices to moment maps of symplectic toroidal actions (,). Indeed, following this route, Bloch Flaschka and Ratiu obtained an *explicit* diffeomorphism between the closure of $`_J`$ and $`𝒫_\mathrm{\Lambda }`$, which we now describe. Let $`S`$ be a symmetric matrix with simple spectrum, and consider the spectral decomposition $`S=Q^T\mathrm{\Lambda }Q`$, where, as usual, $`Q`$ is orthogonal and $`\mathrm{\Lambda }`$ is diagonal. Now, order the diagonal entries of $`\mathrm{\Lambda }`$ in descending order: $`Q`$ is then defined up to a choice of sign for each column of $`Q^T`$ (since they are normalized eigenvectors of $`S`$). Thus, $`S=Q^TD\mathrm{\Lambda }DQ`$, for some diagonal matrix $`D`$ of signs. The matrix $`DQ\mathrm{\Lambda }Q^TD`$ is dependent on $`D`$, but its diagonal is not! The map $`S\text{diag}(Q\mathrm{\Lambda }Q^T)`$ is the *BFR map*. Its restriction to the closure of a Jacobi slice $`_J`$ is the required diffeomorphism to $`𝒫_\mathrm{\Lambda }`$ ().

The proof of the result above makes use of sophisticated machinery. However, once we knew what had to be proved, a simpler argument appeared, yielding a more general result (). Let $`S`$ be an arbitrary real, symmetric, irreducible matrix with simple spectrum. Call the diagonal matrices in the closure of $`_S`$ *accessible*, and their images under the BFR map the *accessible vertices*. Now let $`𝒫_S`$, the *spectral polytope* of $`S`$, be the convex closure of the accessible vertices of $`S`$.

Theorem:  The BFR map is a diffeomorphism between the closure of the slice $`_S`$ and the spectral polytope $`𝒫_S`$.

Generically, the spectral polytope of a symmetric matrix with simple spectrum is indeed the permutohedron associated to the $`n!`$ diagonal arrangements of the eigenvalues. But this is not the case in general. A self-contained description of the spectral polytope of an irreducible matrix $`S`$ is as follows. Consider the spectral decomposition, $`S=Q^T\mathrm{\Lambda }Q`$. We denote by $`Q_{\{r_1,\mathrm{},r_k\},\{c_1,\mathrm{},c_k\}}`$ the minor of $`Q`$ consisting of the entries in the intersection of rows with indices $`r_1,\mathrm{},r_k`$ with columns with indices in $`c_1,\mathrm{},c_k`$. For a permutation $`pi`$ in $`n`$ symbols, let $`\mathrm{\Pi }`$ denote the permutation matrix with entries $`\mathrm{\Pi }_{i,j}=\delta _{i,\pi (j)}`$.

Theorem:  A diagonal matrix $`\mathrm{\Lambda }_\pi =\mathrm{\Pi }^T\mathrm{\Lambda }\mathrm{\Pi }`$ is an accessible vertex of the slice $`_S`$ if and only if the minors $`Q_{\{\pi (1)\},\{1\}},Q_{\{\pi (1),\pi (2)\},\{1,2\}},\mathrm{},Q_{\{\pi (1),\mathrm{},\pi (n)\},\{1,\mathrm{},n\}}`$ have nonzero determinant.

As an example, let $`S=Q^T\mathrm{\Lambda }Q`$ be a $`3\times 3`$ matrix with eigenvalues $`\mathrm{\Lambda }=\text{diag}(4,2,1)`$, and so that the entry $`Q_{1,1}`$ equals zero, but no other minor of $`Q`$ has determinat equal to zero. The slice $`_S`$ and the spectral polytope have only four vertices: diagonal matrices with entry $`(1,1)`$ equal to 4 do not belong to $`\overline{_S}`$, and the vertices of the spectral polytope are $`(2,1,4),(2,4,1),(1,2,4)`$ and $`(1,4,2)`$$`\overline{_S}`$ is a quadrilateral.