format: text

text:
# Quasi-Periodic X-ray Oscillations in the Source Sco X-1

##

The RXTE observations of Scorpius X-1 in 1996–1999 are presented. The properties of its quasi-periodic X-ray oscillations are studied in detail. The results obtained are used for analysis in terms of the transition-layer model (TLM) and the relativistic-precession model (RPM) for a slowly rotating neutron star. Theoretical predictions of the two models are compared and their self-consistency is verified. The tilt of the magnetosphere to the accretion-disk plane, the neutron-star mass, and its angular momentum are determined in the model frameworks.

*Key words:* accretion, neutron stars, Scorpius X-1, X-ray sources, low-mass binary X-ray sources, quasi-periodic oscillations.

## 1 Introduction

Scorpius X-1 is a low-mass X-ray binary with an accreting neutron star and one of the brightest X-ray sources. Its spectral properties suggest that Scorpius X-1 belongs to the class of Z-type sources (Hasinger and van der Klis 1989), whose characteristic feature is a Z-shaped track in the color–color diagram. In this interpretation, the spectral properties are presented in the hard–soft color indices, each of which is defined as the harder-to-softer flux ratio in the corresponding energy band. The Z-shaped track is commonly divided into three parts called branches: the horizontal (HB, the upper part of the diagram), normal (NB, the intermediate part), and flaring (FB, the lower part) branches. The position on the Z track is generally believed to be associated with the rate of accretion in the direction from HB to FB. Six sources are currently known to exhibit Z tracks in the color–color diagram: Scorpius X-1, Cygnus X-2, GX 17+2, GX 5-1, GX 340+0, and GX 349+2.

The power-density spectra (Fourier transforms of the flux) of Z-type sources exhibit low-frequency (5–100 Hz) X-ray quasi-periodic oscillation (QPO) peaks. The names of the QPOs correspond to the branch with which their origin is identified: horizontal- (HBO), normal- (NBO), and flaring-branch (FBO) oscillations. HBOs (15–100 Hz) can also be detected in the NB spectral state. However, as one recedes from HB, the statistical significance of the QPO peaks decreases and they become undetectable. When moving along the Z track (from HB to FB) in its NB–FB segment, a QPO peak detectable in the range 5–20 HZ (NBO/FBO) emerges in the power-density spectra. For all the currently known Z-type sources, QPOs were also found in the range 200–1100 kHz (van der Klis 2000). Two kHz QPO peaks ($`\nu _1`$ and $`\nu _2`$ are the lower and upper peaks, respectively) can be simultaneously observed with a frequency difference of $`200`$–400 Hz. For a series of successive observations (i.e., on short time scales from several hours to several days), an increase in the flux is accompanied by an increase in the frequencies of the two peaks, but their difference ($`\mathrm{\Delta }\nu =\nu _2\nu _1`$) is not conserved (van der Klis et al. 1997).

The observations of Scorpius X-1 revealed all the QPO types characteristic of the low-frequency ($`<100`$ Hz) range in the power-density spectrum. At high frequencies, two kHz QPO peaks ($`\nu _1and\nu _2`$) were found in the RXTE/PCA data \[first detected by van der Klis et al. (1996)\]. Here, we study the QPO properties in detail. Our results are used for analysis in terms of the transition-layer model (TLM) and the relativistic-precession model (RPM). We compare theoretical predictions of the two models and verify their self-consistency. The X-ray flux variability of the source at low frequencies (0.1–128 Hz) is investigated by taking into account the power-law behavior of the power-density spectrum at frequencies above and below the break frequency.

## 2 DATA AND OBSERVATIONS

For our time analysis, we used data from the PCA (Proportional Counter Array) instrument (Jahoda et al. 1996) onboard the RXTE observatory (Bradt et al. 1993) retrieved from the Goddard Space Flight Center Electronic Archive.

The source Scorpius X-1 was observed by the RXTE observatory during eleven series of pointing observations (10056, 10057, 10059, 10061, 20053, 20426, 30035, 30036, 30406, 40020, 40706): in February and May 1996, in March, April, and August 1997, in January, February, and from May until July 1998, and in January and July 1999. The observations of Scorpius X-1 over this period correspond to three different observational epochs of RXTE/PCA (1, 3, and 4 in the adopted classification), for which the boundaries of the PCA energy channels were changed.

To construct the power-density spectra, we used observational data with a resolution of $``$122 or $`244`$ $`\mu `$m ($`2^{13}`$, $`2^{12}`$ s) from the zeroth to 87th or 249th PCA energy channels. These ranges (0–87 and 0–249) correspond to a detectable photon flux $`1.523.5`$, $`1.932.3`$, and $`2.038`$ keV for epochs 1, 3, and 4, respectively, or extend to $`60`$ keV if the entire accessible PCA energy range (channels 0–249) was used.

Of all the observations, we used only those during which the angle between the source direction and the Earth’s horizon was more than 10. During the observations of Scorpius X-1, all five proportional counters were not always switched on to record events. If the operating condition of one of the counters changed during a continuous observation (whose duration did not exceed the duration of one orbit and was, on average, $`33.5\times 10^3`$ s), then the time interval during which the total count rate changed abruptly was excluded from our analysis.

To analyze the variability of Scorpius X-1, we constructed its power-density spectra (van der Klis 1989) in the frequency range 0.03125–2048 Hz (see Fig. 1). The properties of the low- and high-frequency flux variability were investigated in the bands 0.1–256 and 256–2048 Hz, respectively. No correction was made for the background radiation and dead time.

We analyzed individual observations with a duration up to $`3.5\times 10^3`$ s. Fitting the power-density spectra by a constant and by a power law at frequencies below and above the break frequency did not yield acceptable results (according to the $`\chi ^2`$ test). The main reason was the absence of a sharp break and the resulting uncertainty in its measurement. The model in which at frequencies much higher ($`\nu /\nu _{break}1`$) and much lower ($`\nu /\nu _{break}1`$) than the break, each part of the spectrum could be fitted by its own power law and the transition between them was not jumplike (previously, this model was successfully used for a time analysis of the Z-type source Cygnus X-2; see, e.g., Kuznetsov 2001), proved to be more suitable:

$$P(\nu )=A\nu ^\alpha [1+(\nu /\nu _{\text{break}})^\beta ]^1.$$
(1)
The power-density spectra (see Fig. 1) were fitted by this model in the $`0.1`$$`128`$-Hz band with the additional introduction of one or two Lorenz lines to allow for the QPO peaks and their harmonics. To take into account the PCA dead-time effect, which causes the overall level to be shifted to the negative region \[because of this effect, the Poissonian noise level subtracted from all spectra differed from $`2.0`$, in units of the Leahy normalization; see van der Klis (1989) and Vikhlinin et al. (1994) for more details\], we added a constant to the general model.

In searching for the kHz (in the range $`5001200`$ Hz) QPO peaks $`\nu _1`$ and $`\nu _2`$, we analyzed the power-density spectra at high frequencies. A constant with the addition of one or two Lorenz lines was used as the model. Individual observations in which the detection confidence level of kHz QPO peaks exceeded $`3\sigma `$ were used for the subsequent analysis.

Our analysis of the low-frequency flux variability in Scorpius X-1 revealed QPOs at a high confidence level in all the individual observations that corresponded to the HB or NB spectral states. Note that the HBO frequency is one of the main frequencies in most models, including the TLM and RPM. The two models (see below) unequivocally establish the correspondence between the three QPO peaks: the HBOs, lower, and upper kHz QPOs.

To obtain the most accurate model parameters (the angle $`\delta `$ for the TLM and the neutron-star mass $`M_{NS}`$ and angular momentum $`a`$ for the RPM), which are invariants in each of the theories under consideration, it was necessary to calculate the three simultaneously observed main QPO frequencies (i.e., $`\nu _{HBO},\nu _1,\nu _2`$; see Fig. 1) as accurately as possible. Therefore, of all the preselected data in which both kHz QPO peaks ($`\nu _1`$ and $`\nu _2`$) were found, we used only those in which the main HBO peak ($`\nu _{HBO}`$) and its second harmonic ($`\nu _{2HBO}`$) were detected at a confidence level above $`4\sigma `$. In such observations, the two kHz QPO peaks were also reliably detected.

## 3 LOW-FREQUENCY X-RAY FLUX VARIABILITY

### 3.1 The Power of the HBO peaks

In Fig. 2, rms is plotted against HBO frequency. The flux variability (i.e., $`rms^2`$) at the peak corresponds to the integral of the Lorenz line that fits best the QPO peak in the power-density spectrum. For the Z-type sources, to which Scorpius X-1 belongs, an increase in HBO frequency is accompanied by a decrease in rms (see, e.g., van der Klis 2000). As a result, the detection confidence level of HBO peaks is lower at high frequencies ($`5060`$ Hz). Figure 2 shows the range of HBO frequencies found for Scorpius X-1. During all observations, the deviation from the mean $`\nu _{HBO}43.5`$ Hz did not exceed $`5`$ Hz (i.e., it changed by no more than $`10\%`$ during all observations).

### 3.2 Horizontal-Branch Oscillations and Their Harmonics

In Scorpius X-1, quasi-periodic HBO peaks were found at low frequencies ($`<100`$ Hz). Together with the oscillations at the main frequency $`\nu _{HBO}`$, we also found QPO peaks at frequencies close to the multiple ones. Note that the detection confidence level of HBO peaks decreases with frequency (see Fig. 2). That is why no multiple peaks were found (at a confidence level $`>4\sigma `$) for which $`\nu _{HBO}>48`$ Hz (the corresponding oscillation frequency of the second harmonic is $`\nu _{2HBO}>96`$ Hz).

The frequency ratio $`r_{2/1}=\nu _{HBO}/\nu _{2HBO}`$ of the multiple HBO harmonics is shown in Fig. 3. We clearly see that $`r_{2/1}`$ differs from 2.0: only a third of the 39 data points are above 2.0 (the dotted line in Fig. 3). Fitting by a constant yields $`r_{2/1}=1.965\pm 0.004`$, $`\chi _{red}^2=3.45`$ (in what follows, $`\chi _{red}^2=\chi ^2/d.o.f`$ is the reduced $`\chi ^2`$ value per degree of freedom), while fitting by a constant with a break yields $`r_{2/1}=1.976\pm 0.004`$, $`\chi _{red}^2=1.05`$ (and $`r_{2/1}=0.05\times (\nu _{HBO}45.31)`$ in the range of a linear dependence). When the scatter of data points exceeds the stochastic scatter and when the $`\chi _{red}^2`$ value is much larger than unity (as in the case of a constant fit), it should be assumed that the uncertainty could be produced by stochastic processes rather than by the statistics alone. We can then ignore the measurement error and determine the mean $`\overline{r}_{2/1}=\frac{1}{N}r_{2/1}`$ and the rms deviation from it $`\sigma =\sqrt{\frac{{\scriptscriptstyle (r_{2/1}\overline{r}_{2/1})^2}}{(N1)}}`$. In this case, we obtained the following ratio of the second and first harmonics: $`\overline{r}_{2/1}=1.973`$, $`\sigma =0.050`$. Under such an assumption, the harmonic ratio is compatible with 2.0, within the error limits.

Let us consider in more detail the case where the scatter of $`r_{2/1}`$ values about the mean is determined by the statistics alone and is incompatible with 2.0, within the error limits. Resonances in nonlinear oscillatory systems can emerge at frequencies at which an external periodic force acts on the oscillator. This requires that its frequency $`\gamma `$ satisfy the condition

$$\gamma =\omega _0\frac{p}{q},$$
(2)
where $`p`$ and $`q`$ are integers and $`\omega _0`$ is the cyclic frequency of free system oscillations.

The power of the main peak must be maximal among all of the harmonics, because the intensity of resonance phenomena rapidly decreases with increasing $`p`$ and $`q`$ (Landau and Lifshitz 1965). In addition, if the frequency of the periodic external force $`\gamma `$ differs from $`\omega _0`$ by a small value $`\epsilon `$ (in the simplest case, $`p=1`$ and $`q=1`$), then in nonlinear oscillations, this can lead to a displacement of the rms maximum from the oscillator resonance frequency (Landau and Lifshitz 1965). In Scorpius X-1, the largest rms among all the HBO harmonics corresponds to the frequency $`\nu _{HBO}`$ (see Fig. 2). We assume from the outset that the latter (and this assumption proved to be valid) was the fundamental harmonic and the deviation of the frequency ratio of the second and first harmonics from 2.0 is admitted in the forced oscillations of a nonlinear oscillator under the action of a periodic force.

Note that in this case, the frequency of the observed low-frequency oscillation peak is $`\nu _{HBO}\omega _0/2\pi `$. This is particularly important in comparing the self-consistency of the theoretical models under consideration, which interpret the emergence of QPOs (to be more precise, the following three QPOs: $`\nu _{HBO}`$, $`\nu _1`$, and $`\nu _2`$) in binary systems with neutron stars.

If the deviation of the ratio of the multiple HBO harmonics $`\nu _{2/1}`$ from 2.0 is assumed to be attributable to forced nonlinear oscillations, then the observed $`\nu _{HBO}`$ must be close but not exactly equal to $`\omega _0/2\pi `$, the assumed frequencies of free oscillations of the oscillator. Note that $`\omega _0/2\pi `$ rather than $`\nu _{HBO}`$ corresponds to the low-frequency QPOs in the models considered below.

For this reason, for our analysis in terms of each of the models considered below, we use two values for the low-frequency QPOs: $`\nu _{HBO}`$ and $`\overline{\nu }_{HBO}`$. These values are the fundamental HBO harmonic and the mean between half the frequency of the second harmonic $`\nu _{2HBO}`$ and the fundamental harmonic $`\nu _{HBO}`$. For the latter case, from all the observations, we chose those in which the confidence level of each HBO peak was higher than $`4\sigma `$.

## 4 THEORETICAL MODELS

For our analysis in terms of the two models in question, we used the average spectra and the main frequencies: $`\nu _{HBO}`$ (as well as $`\overline{\nu }_{HBO}=1/2(\nu _{HBO}+1/2\nu _{2HBO})`$, the mean between the fundamental harmonic and half the second harmonic), $`\nu _1,\nu _2`$. The following approximate asymptotic relations between the frequencies are known from the observations of low-mass X-ray binaries (Stella et al. 1999): (1) $`\mathrm{\Delta }\nu =\nu _2\nu _1\nu _2^2`$; and (2) $`\nu _{HBO}\nu _1`$. These observed features are consistent with each of the theories under consideration.

### 4.1 The Transition-Layer Model (TLM)

We consider the motion of a clump of matter, the QPO source on the accretion-disk surface, in a Keplerian orbit around a neutron star (NS) in terms of this model. The magnetospheric axis is assumed to be not aligned with the normal to the disk surface but makes an angle $`\delta `$. After multiple passages through a slightly tilted magnetosphere, the clump comes under the effect of Coriolis forces. These forces cause the main Keplerian oscillation frequency $`\nu _K`$ to split up into two oscillation modes: radial ($`\nu _h`$) and perpendicular to the disk plane ($`\nu _L`$). The two modes are the solution of the equation for the rotation of a body in a noninertial frame of reference \[see Osherovich and Titarchuk (1999) for more details on its derivation and solution\]. In this model, the lower and upper kHz QPO peaks correspond to the Keplerian ($`\nu _K\nu _1`$) and hybrid ($`\nu _h\nu _2`$) frequencies, respectively. The relation between the kHz QPO peaks is given by

$$\nu _h=[\nu _K^2+(\mathrm{\Omega }/\pi )^2]^{1/2},$$
(3)
where $`\mathrm{\Omega }`$ is the angular velocity of the magnetosphere. The oscillation mode perpendicular to the disk surface is defined as

$$\nu _L=(\mathrm{\Omega }/\pi )(\nu _K/\nu _h)\mathrm{sin}\delta .$$
(4)
Since all three characteristic frequencies $`\nu _L`$, $`\nu _K`$, and $`\nu _h`$ are known from observations, the sought-for angle $`\delta `$ can be calculated by using (4):

$$\delta =\mathrm{arcsin}[(\nu _h^2\nu _K^2)^{1/2}(\nu _L\nu _h/\nu _K)].$$
(5)
To a first approximation, the angular velocity of the magnetosphere $`\mathrm{\Omega }`$ is constant. A more accurate equation that describes the dependence of $`\mathrm{\Omega }`$ on radius can be derived in the multipole magnetic-field approximation (see Osherovich et al. 1984). Assuming the contribution of the quadrupole component to the magnetic-field strength in the equatorial plane to be negligible and taking into account the dipole and octupole components, the final equation for the angular velocity can be represented as (Osherivich and Titarchuk 1999)

$$\mathrm{\Omega }/2\pi =C_0+C_1\nu _K^{4/3}+C_2\nu _K^{8/3}+C_3\nu _K^4,$$
(6)
with $`C_2=2(C_1C_3)^{1/2}`$.

To reconstruct the magnetospheric profile for Scorpius X-1, according to Eq. (3), it will suffice to have data only for the high-frequency range of the flux variability, i.e., $`\nu _1`$ and $`\nu _2`$. Therefore, we used all the data in which kHZ QPO peaks were found by removing the additional condition that the main HBO peak $`\nu _{HBO}`$ and its second harmonic $`\nu _{2HBO}`$ should be detected in the power-density spectrum at low frequencies at a confidence level above $`4\sigma `$.

The angular velocity is plotted against Keplerian frequency in Fig. 4. Using (6), we managed to obtain the magnetospheric profile suggested by the model under consideration with the following parameters: $`C_0=354`$ Hz, $`C_1=3.54\times 10^2`$ Hz-1/3, $`C_2=9.99\times 10^6`$ Hz-5/3, and $`C_3=7.21\times 10^{10}`$ Hz-3. In fitting the data, we assumed $`C_2`$ to be a free parameter, whereas in the model, $`C_2=`$2$`(C_1C_3)^{1/2}=`$1.01$`\times `$10-5. The similar values of $`C_2`$ determined by two independent methods are a weighty argument for the approximation used for the magnetospheric profile.

Using (5), we obtained $`\delta `$ for each triplet of frequencies: $`\nu _{HBO}`$, $`\nu _1`$, and $`\nu _2`$. In Fig. 5 the HBO frequency $`\nu _{HBO}`$ is plotted against the angle $`\delta `$ between the normal to the disk surface and the magnetospheric axis. Among all three frequencies used in calculating $`\delta `$, the strongest correlation is observed between $`\delta `$ and $`\nu _{HBO}`$. For this reason, the dependence of $`\delta `$ on low-frequency QPOs is analyzed in the model under consideration (below, we show that this is also true for the relativistic precession model), because the model invariants must be conserved, irrespective of the measured QPO frequency range.

Table 1 gives data on the constant fit to $`\delta `$ in two cases: when $`\nu _L=\nu _{HBO}`$ (see Fig. 5) and $`\nu _L=\overline{\nu }_{HBO}=1/2(\nu _{HBO}+1/2\nu _{2HBO})`$ (see Fig. 6), i.e., with the inclusion of the second HBO harmonic. The $`\chi _{red}^2`$ value is given in each case. The mean angle $`\overline{\delta }`$ and its rms deviation $`\sigma _D`$ are also given in Table. 1.

To reduce the possible stochastic scatter of data points (which probably dominates over the statistical scatter), we performed averaging depending on the HBO frequency. The upper panels in Figs. 5 and 6 show the data obtained in this way. Whereas the data in Fig. 5 reveal a clear correlation between $`\delta `$ and $`\nu _{HBO}`$, the data points in Fig. 6 virtually lie on the same line. If a larger number of points in each bin were averaged (in our case, we combined them by five), then $`\delta `$ (see Fig. 6, where $`\delta `$ is plotted against $`\overline{\nu }_{HBO}`$) would be even more compatible with a constant, according to the $`\chi ^2`$ test.

Let us compare the rms deviation between the HBO harmonic ratio $`r_{2/1}`$ determined above and $`\delta `$ (see Table 1). It turns out that $`\sigma _{r_{2/1}}2.5\%`$, while $`\sigma _\delta 2.0\%`$. Since, according to (5), $`\delta \nu _{HBO}`$ at small angles and since $`\nu _{HBO}`$ increasingly deviates from $`\omega /2\pi `$ to larger values with decreasing $`r_{2/1}`$, the correlation for $`\delta `$ (clearly seen only in Fig. 5) can result from the emergence of nonlinear oscillations considered in the previous sections.

### 4.2 The Relativistic Precession Model (RPM)

The high-velocity motion of matter in a strong gravitational field can generate oscillations attributable to general-relativity effects. The RPM (see, e.g., Morsink et al. 1999) considers the motion of a point mass around a gravitating center. If a particle moves in an orbit that does not lie exactly in the equatorial plane of a compact object but is inclined at an infinitesimal angle, then the particle orbit will precess.

For a nonrotating gravitating center (i.e., in the Schwarzschild approximation), the expression for the particle angular velocity matches the classical formula for Keplerian motion. However, for a NS or a black hole with an intrinsic angular momentum, the azimuthal frequency can be expressed (Bardeen et al. 1972; Stella et al. 1999) in a system of units with $`G=1`$ and $`c=1`$ as

$$\nu _\varphi =\sqrt{\frac{M}{r^3}}\left[2\pi \left(1+a\sqrt{\frac{M}{r^3}}\right)\right]^1.$$
(7)
Below, we consider only the case where the compact object and the test particle corotate: $`\nu _\varphi >0,a>0`$ ($`a`$ is the relative angular momentum). The epicyclic frequency $`\nu _r`$ together with the azimuthal frequency $`\nu _\varphi `$ determine the orbital periastron rotation, $`\nu _{per}\nu _\varphi \nu _r`$, while the frequency in a direction perpendicular to the disk plane $`\nu _\theta `$ determines the nodal precession, $`\nu _{nod}|\nu _\varphi \nu _\theta |`$. The corresponding equations were derived by Okazaki et al. (1987) and Kato (1990):

$$\nu _r=\nu _\varphi \sqrt{16\frac{M}{r}+8a\sqrt{\frac{M}{r^3}}3\frac{a^2}{r^2}},$$
(8)
$$\nu _\theta =\nu _\varphi \sqrt{14a\sqrt{\frac{M}{r^3}}+3\frac{a^2}{r^2}}.$$
(9)
In contrast to the TLM, the Keplerian rotation frequency in the RPM corresponds to the upper QPO peak, $`\nu _2=\nu _\varphi `$, while the observed periastron precession frequency corresponds to the lower QPO peak, $`\nu _1=\nu _{per}\nu _\varphi \nu _r`$. In the Kerr approximation, $`\nu _\varphi \nu _\theta `$ (for the Schwarzschild case with $`a=0`$, the frequencies are equal: $`\nu _\varphi =\nu _\theta `$), the HBOs are identified with nodal precession: $`\nu _{nod}|\nu _\varphi \nu _\theta |`$. Some authors (e.g., Stella et al. 1999) believe that $`\nu _{HBO}`$ is an even harmonic of $`\nu _{nod}`$. In contrast to $`\nu _{HBO}=\nu _{nod}`$, this relation causes the mass of the compact source in Eqs. (7)–(9) required by the RPM to decrease.

In the RPM approximation under consideration, the NS mass and relative angular momentum are invariants and do not depend on Keplerian frequency. Thus, we can determine $`M_{NS}`$ and $`a`$ from Eqs. (7)–(9) using the three observed main frequencies: $`\nu _{HBO}`$, $`\nu _1`$, and $`\nu _2`$.

As in the section devoted to the TLM, here, the invariant (more specifically, the NS mass $`M_{NS}`$ suggested by the model) was also determined in two cases: $`\nu _{nod}=\nu _{HBO}`$ and $`\nu _{nod}=\overline{\nu }_{HBO}=1/2(\nu _{HBO}+1/2\nu _{2HBO})`$ (i.e., with the inclusion of the second HBO harmonic). Since the derived $`M_{NS}`$ proved to be most dependent (as in the TLM) on low-frequency HBOs, Figs. 7 and 8 show $`M_{NS}`$ as a function of $`\nu _{HBO}`$.

Table 2 gives a constant fit to the data (shown in Figs. 7 and. 8). We clearly see a large $`\chi _{red}^2`$ value for each $`\nu _{nod}`$: compared to Table. 1, which gives a fit to $`\delta `$ (in the TLM), the corresponding $`\chi _{red}^2`$ values are larger by a factor of $`36`$. Nevertheless, the rms deviation $`\sigma _{M_{NS}}`$ is also $`2.0\%`$, as is $`\sigma _\delta `$ in the TLM. In contrast to Fig. 4, in which $`\delta `$ is plotted against $`\overline{\nu }_{HBO}`$ and is essentially compatible with a constant, we clearly see a linear correlation of $`M_{NS}`$ with $`\overline{\nu }_{HBO}`$ in Fig. 8: as $`\overline{\nu }_{HBO}`$ increases from $`39.5`$ Hz to $`46.5`$ Hz, the required $`M_{NS}`$ decreases from $`2.8M_{}`$ to $`2.6M_{}`$.

The relativistic precession model makes it possible to estimate the compact-source mass and to determine the Keplerian orbital radius. For comparison, Fig. 9 shows the mass–radius relations for various equations of state for neutron stars taken from Miller et al. (1998).

### 4.3 Comparison between the Assumed Invariants of the Models under Consideration

To check whether the two models are self-consistent and to determine whether the derived invariants $`\delta `$ and $`M_{NS}`$ correlate with the HBO frequency (see Fig. 5–9), we used Spearman’s nonparametric correlation test. The coefficient $`r_s`$ equal to unity in magnitude points to a total correlation between the two quantities. In each case under consideration, we give the probability (prob) that the correlation found is the result of statistical fluctuations. This test is also convenient in the case where it is necessary to establish which of the several correlations under consideration is strongest. Our estimates are given in Tables 3 and 4. Note that the data in Table 4 correspond to the case where the statistical scatter was slightly reduced through averaging. Thus, the RPM model (see Table 4) is inconsistent ($`r_s=1`$; the total negative correlation between $`M_{NS}`$ and $`\overline{\nu }_{HBO}`$) and the assumed NS mass is not an invariant in this model.

Allowance for the second harmonic in determining the low-frequency QPO peak (more specifically, $`\overline{\nu }_{HBO}=1/2(\nu _{HBO}+1/2\nu _{2HBO})`$) results in a reduction of the correlation coefficient $`r_s`$ for $`\delta `$ in the TLM, while the reverse is true for $`M_{NS}`$. In particular, we see from Fig. 8 that the data points are well fitted by a straight line: $`\chi _{red}^2=4.31`$ (see Table 2, where $`\chi _{red}^2=16.3`$ for the same case but for a constant fit to the data). Figures 6 and 8 (upper panels) show that averaging can reduce the stochastic scatter when determining the model invariants. However, for the TLM (see Fig. 6), $`\delta `$ is virtually constant, $`\overline{\delta }=5.55`$, with the rms deviation $`\sigma _{\overline{\delta }}=0.03`$ (which is $`0.5\%`$), while for the RPM (see Fig. 8), $`M_{NS}`$ strongly correlates with $`\overline{\nu }_{HBO}`$ and is incompatible with a constant (the rms deviation is $`\sigma _{\overline{M}_{NS}}1.7\%`$, which roughly corresponds to $`\sigma _{M_{NS}}`$; see Table 2).

## 5 DISCUSSION AND RESULTS

Both theoretical models suggest the existence of invariants. These invariants are the tilt of the magnetosphere to the disk normal $`\delta `$ for the TLM and the NS mass $`M_{NS}`$ and its relative angular momentum $`a`$ for the RPM in the Kerr approximation. Since the RPM did not consider the quadrupole moment of the compact source and disregarded its oblateness (which can give a contribution of $`1015\%`$ to the nodal precession frequency), we managed to avoid the inclusion of various theoretical equations of state for neutron stars in our analysis. However, it was of interest to compare our results with the constraints imposed on the equations of state. Figure 9 shows some of these constraints taken from Miller et al. (1998). The curves correspond to the mass–radius relation for nonrotating neutron stars. Although the angular velocity for Scorpius X-1 is nonzero, the correction for spin (here, $`\nu _s200400`$ Hz) can be disregarded.

The assumed NS mass $`M_{NS}`$ and the Keplerian orbital radius $`R`$ of the particle rotating around it (see Fig. 9) derived from Eqs. (7)–(9) must be consistent with the constraints imposed by the equations of state. More specifically, $`M_{NS}<M_{EOS}`$, $`R>R_{EOS}`$ (here, $`M_{EOS}`$ and $`R_{EOS}`$ are the NS mass and radius for various equations of state). Only the equation of state “L” suggest that a massive NS (in our case, $`2.7M_{}`$) could exist and it agrees with the derived $`M_{NS}`$ in the two cases ($`\nu _{nod}=\nu _{HBO}`$ and 2$`\nu _{nod}=\nu _{HBO}`$). At the same time, the other two equations of state (“A” and “UU”) are incompatible with the RPM results in the Kerr (i.e., with NS rotation) approximation. Moreover, we obtained $`M_{NS}`$ for Cygnus X-2 that was compatible, within the error limits ($`M_{NS}=2.69\pm 0.10`$ for $`\nu _{nod}=\nu _{HBO}`$ and $`M_{NS}=2.29\pm 0.09`$ for $`2\nu _{nod}=\nu _{HBO}`$; Kuznetsov 2002), with $`M_{NS}`$ for Scorpius X-1 (see Table 2). Such a coincidence is not an argument for the relativistic precession model either.

In addition, we showed that a nonparametric test for correlation between the assumed NS mass in the RPM and the HBO frequency yields a positive result. Thus, the RPM is not a self-consistent model in which the invariants are conserved and the derived mass of the compact object is large enough for the existing equations of state for neutron stars. At the same time, for the transition-layer model, the derived value is compatible with a constant under certain conditions (e.g., using the theory of nonlinear oscillations; for more details, see Kuznetsov and Titarchuk 2002). If the uncertainty in the calculated angle $`\delta `$ does not exceed $`2\%`$ (which is lower than the rms deviation for $`\delta `$ in Scorpius X-1; see Table 1), then its value may be assumed to be constant.

The observed inconsistency of the RPM may stem from the fact that the orbital inclination of the test particle in the model simplification is assumed to be infinitesimal. An exact solution for an arbitrary inclination is given in Sibgatullin (2001). In such an analysis, the frequency $`\nu _{HBO}`$ for fixed mass and angular momentum can change by a factor of $`3`$ as the inclination changes from 0 to $`\pi /2`$ (for the marginally stable orbit at $`M_{NS}=2M_{}`$ and $`\nu _\varphi =1200`$ Hz). This dispenses with the need to take the nodal precession frequency $`\nu _{nod}`$ to be equal to 1/2 of $`\nu _{HBO}`$ to obtain more acceptable NS-star masses.

In addition, it is worth noting that the angle $`\delta `$ in the RPM may not be strictly an invariant. Higher energy release in the disk can affect the curvature of its inner region. The tilt angle $`\delta `$ between the magnetospheric axis and the normal to the disk surface can then differ for different source fluxes. As a result, the angles $`\delta `$ can differ for the same $`\nu _{HBO}`$ (see Figs. 5 and 6).

## ACKNOWLEDGMENTS

This study was supported in part by the Program of the Russian Academy of Sciences “Astronomy: Nonstationary Astronomical Phenomena”, the Russian Foundation for Basic Research (project no. 02-02-17347), and grant no. 00-15-99297 of the President of Russia. I used the RXTE data retrieved from the Goddard Space Flight Center Electronic Archive. I wish to thank L. Titarchuk, B. Stone, N. White, J. Swank, Ph. Newman, and J. Repaci for the opportunity to work with the RXTE archival data on CR-ROMs. I am grateful to the referee for helpful remarks and comments.

Translated by V. Astakhov


text:
# The largest black holes and the most luminous galaxies

## 1 Introduction

Recent progress in reverberation mapping of active galacic nuclei (AGNs) allowed the first meaningful correlation between the broad line region (BLR) size ($`R_{\mathrm{BLR}}`$) and the black hole (BH) mass in more than 30 objects. This provided a simple way to calculate BH masses for a large number of sources and resulted in a flood of papers on this topic. Some papers (e.g. Vestergaard 2002, hereafter V02; McLure and Jarvis 2002) investigated, in great detail, the wavelength dependence of the $`R_{\mathrm{BLR}}LM`$ relationship and provided useful ways for adopting the method to other wavelength bands. This opens the way for the study of BH masses in large samples of high luminosity high-z quasars.

All the new BH mass estimates are based on a single relationship obtained for a single sample of 34 AGNs for which BLR sizes are available from decade long reverberation mapping campaigns. More than half the sample was observed at the Wise observatory over a period of about 12 years (Kaspi et al. 2000, hereafter K00). Other objects have been monitored in other observatories and in several “AGN watch” campaigns (Netzer & Peterson 1997; Peterson 2001). The main findings are a significant $`R_{\mathrm{BLR}}\lambda L_\lambda (5100)`$ relationship ($`L_\lambda (5100)`$ is the monochromatic luminosity at 5100Å) and the confirmation that the BLR gas is in virial motion (e.g. Peterson and Wandel 2000). These, plus the (model dependent) conversion of the observed full-width at half maximum (FWHM) of various emission lines into 3-D gas velocities, are sufficient to derive the mass of the central BH.

This letter discusses the mass of the largest BH in the universe; those found in the centers of the most luminous quasars. It follows the works of Laor (1998; 2001), McLure & Jarvis (2002), Woo & Urry (2002) and others who used such methods for obtaining BH masses beyond the original K00 sample. The paper addresses also the Shields et al. (2002) new results and extends the mass estimates to much larger quasar samples. Section 2 presents new mass calculations for a large number of sources and §3 illustrates the new correlations found. Section 4 discusses the new results in light of the available information on the largest, most luminous elliptical galaxies and the epochs of quasars and galaxy formation.

## 2 The largest BH

### 2.1 BH Mass measurements

New mass estimates have been obtained for a large number of AGNs using the $`R_{\mathrm{BLR}}L`$ relationship obtained from the K00 sample; the only sample available for such calibration. This relationship is given, schematically, by

$$R_{\mathrm{BLR}}=c_1L_\lambda ^\gamma $$
(1)
which results in the following mass estimate:

$$M_{\mathrm{BH}}=c_2L_\lambda ^\gamma [FWHM]^2.$$
(2)
Here $`c_1`$ and $`c_2`$ are constants that include the flux normalization and various assumptions about the velocity field in the BLR. The slope $`\gamma `$ is derived from the reverberation campaigns results and is in the range 0.5-0.7 (see below). The expression in eqn. 2 can be used to derived “single epoch” masses that combine the constants $`\gamma `$ and $`c_2`$ with observed FWHM of certain emission lines in individual objects. The method has been described in various papers including K00, V02, and McLure & Jarvis (2002). Its more useful applications are based on measured $`\lambda L_\lambda `$(5100) and FWHM(H$`\beta `$) for low redshift sources (the quantities used by K00) and the combination of $`\lambda L_\lambda `$(1350) and FWHM(C iv$`\lambda 1549`$) for high redshift objects. V02 has looked into the inter-calibration of the two and supplied the expressions that are used in this work except for a small correction in the value of $`c_2`$ that was introduced to adjust her constants to the cosmology assumed here: $`H_0=70`$ km s-1 Mpc-1, $`\mathrm{\Omega }_m=0.3`$ and $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$. McLure & Jarvis (2002) provided similar expressions for $`MgII\lambda 2798`$ which are not used in this work.

### 2.2 The sample

Three AGN samples have been used in this work: 1) the LBQS sample (Forster et al. 2000 and references therein), 2) A sample of 104 high redshift high luminosity quasars with ground-based spectrophotometry ($`L_\lambda (1350)`$) and good FWHM(C iv$`\lambda 1549`$) measurements, and 3) the new $`L_\lambda (5100)`$ and FWHM(H$`\beta `$) listed by Shields et al (2002). Many of the sources in the second sample are UM quasars and the raw data can be found in MacAlpine and Feldman 1982, Baldwin, Wampler & Gaskell (1987) and Baldwin (1977).

Forster et al. (2000) supplied monochromatic luminosities and FWHMs for many emission lines in about 1000 LBQS quasars. Since many of the sources have been observed through relatively small aperture, and under poor weather conditions, it was decided to use the Bj magnitudes that are much more accurate. This follows Green et al. (2001) who studied the Baldwin relationship in this sample and obtained monochromatic luminosities using the same method. All fluxes have been corrected for galactic reddening using the Green et al. procedure. A major assumption here, and in Green et al. (2001), is that the observed continuum can be described by a single $`L_\nu \nu ^\alpha `$ power-law with $`\alpha =0.5`$. This approximation neglects the possible dependence of $`\alpha `$ on source luminosity which may affect the $`LM`$ relationship (see §4). Forster et al. (2001) provided several different measurements of FWHM(C iv$`\lambda 1549`$) with and without the narrow line component. The “single” component fit was used and the “broad only” fit was checked to verify that the results are not sensitive to this choice. A handful of sources with FWHM(H$`\beta `$)$`<1,000`$ km s-1 or with FWHM(C iv$`\lambda 1549`$)$`>20,000`$ km s-1 were removed from the sample since those were considered unreliable or affected too much by the narrow emission line. As for the second C iv$`\lambda 1549`$ sample, no galactic reddening was applied and the same assumption about $`L_\nu `$ was used. In this case there is no significant dependence on $`\alpha `$ since the original papers quote the observed flux at around rest wavelength of 1450Å.

The above samples are optically selected and suffer from various selection effects. This is of no real consequence to the main goal of the paper which is to derive the mass of the largest known BHs. It may affect, however, the derived $`ML`$ correlations (§4).

## 3 The $`LM`$ relationship for high luminosity AGNs

BH masses have been calculated using equation 2 and the normalizations derived by K00 and V02 adjusted to the cosmology chosen here. The determination of the slope $`\gamma `$ is crucial for the present work and will be discussed prior to presentation of the new results.

We start from the original K00 sample to which we apply two statistical methods for finding $`\gamma `$: the Akritas & Bershady (1996) BCES estimator (for which we only consider the BCES bisector) and the fitexy method described in Press et al. (1992). The merits of the different methods have been discussed, extensively, in several papers and will not be repeated here. Our experience shows that that the differences between the slopes obtained by the different methods are larger than the formal uncertainties on the slopes of each method. The K00 sample adjusted to the new cosmology gives $`\gamma =0.58\pm 0.12`$ for the BCES bisector estimator and $`\gamma =0.68\pm 0.03`$ for the fitexy method. The two are formally consistent with each other and $`\gamma (\mathrm{BCES})`$ was adopted here.

Since the purpose of this work is to extrapolate to very large $`L`$, we also experimented with removing the lowest luminosity objects from the sample. Removing the three objects with $`\lambda L_\lambda (5100)<10^{43}`$ ergs s-1 resulted in $`\gamma (\mathrm{BCES})=0.71\pm 0.21`$ and $`\gamma (fitexy)=0.69\pm 0.03`$. Removing the seven objects with $`\lambda L_\lambda (5100)<10^{43.7}`$ ergs s-1 resulted in $`\gamma (\mathrm{BCES})=0.58\pm 0.19`$ and $`\gamma (fitexy)=0.74\pm 0.04`$. All these results suggest that the two methods are consistent with each other and the slope cannot be determined to an accuracy better than about 0.15. The value adopted for illustrating the results of this work is the smaller one found for the entire K00, $`\gamma =0.58`$. The implications for the case of larger or smaller $`\gamma `$ are discussed in §4.

Shields et al. (2002) suggested the use of the “physically motivated” value of $`\gamma =0.5`$. The strongest argument for using this value is the suggestion by Netzer & Laor (1993) that the outer boundary of the BLR is determined by the dust sublimation radius which is similar to the measure $`R_{\mathrm{BLR}}`$ to within a factor $`2`$. There are several problems in applying this idea to the present mass determination. First, the “reverberation radius” is determined by the responsivity of H$`\beta `$ to changes in the ionizing luminosity, $`L_{\mathrm{ion}}`$, which is smaller than the bolometric luminosity that determines the dust sublimation radius. In addition, $`\gamma =0.5`$ means the same BLR ionization parameter for low luminosity Seyferts and the highest luminosity quasars. This has never been shown to be the case in large QSO samples. Thus, more work is required to justify this theoretical value of $`\gamma `$.

The masses computed with the $`\gamma =0.58`$ slope are presented in Fig. 1. The diagram contains data for 505 QSOs with C iv$`\lambda 1549`$ measurements and 219 source with H$`\beta `$ measuremets. The luminosity range is roughly $`\lambda L_\lambda (1350)=10^{4447.5}`$ ergs s-1. Also shown is the best regression line (see below) and the mass range of $`\pm \sigma _M`$ around the median calculated in luminosity bins of 0.3 dex. The largest BHs are found in sources with $`z>2`$ with $`M_{\mathrm{BH}}10^{10.2}`$ $`M_{}`$ (15 with mass exceeding $`10^{10}`$ $`M_{}`$). Using $`\gamma =0.68`$ (the slope found with the $`fitexy`$ method), raise this number to about $`10^{10.4}`$ $`M_{}`$ (62 with mass exceeding $`10^{10}`$ $`M_{}`$).

Fig. 2 shows $`M_{\mathrm{BH}}`$ as a function of redshift for the same sample under the same assumptions.

The data in Fig. 1 suggest a simple linear dependence of the form $`M_{\mathrm{BH}}L^\beta `$. This has been tested by performing a linear regression analysis using the same two methods described earlier. The procedure used for calculating the errors is the following: For $`L_{1350}`$, the assumption is of a constant error of 0.15 dex representing the measurement uncertainty, the extrapolation in wavelength and the typical range in luminosity due to continuum variability. This number does not affect the resulting slope $`\beta `$ in any significant way. As for the mass, this was done using standard error propagation combining all errors due to the uncertainties in $`L`$ and in FWHM (line width uncertainties are given in Forster et al. 2000). The combined error for this case is typically 0.15-0.25 dex. No uncertainties are listed for FWHM(C iv$`\lambda 1549`$) in the second quasar sample and for FWHM(H$`\beta `$) in Shields et al. (2002). A uniform error of 0.2 dex in $`M_{\mathrm{BH}}`$ was assumed in those cases. The errors are relatively large and are expressed in logarithmic form (i.e. 0.5(log(x+dx)$``$log(x$``$dx), see Lyons, 1991).

Table 1 lists several slopes obtained by the two methods for our standard case of $`\gamma =0.58`$ and for $`\gamma =0.68`$, the slope obtained by the $`fitexy`$ method. Given the various biases and unknowns, it is reasonable to assume that the real uncertainty in $`\beta `$ is at least as large as $`\pm 0.15`$. With this uncertainty, the slopes of the C iv$`\lambda 1549`$ sample and the entire sample are barely consistent with each other and the slopes of the H$`\beta `$ sample and the entire sample are indistinguishable. The scatter in slope is probably due to the very different luminosity range of the C iv$`\lambda 1549`$ and the H$`\beta `$ samples. A second approach that was tries assumed a uniform uncertainty in $`M_{\mathrm{BH}}`$ of 0.3 dex for all objects. This gave very similar results. The overall conclusion is that for luminous AGNs, $`M_{\mathrm{BH}}L^{0.9\pm 0.15}`$.

The $`M_{\mathrm{BH}}L`$ correlation found here is very different from the one found in K00. The reason is probably the incompleteness of the small K00 sample which resulted in a biased sampling of FWHM(H$`\beta `$) vs. $`\lambda L_\lambda (5100)`$ not representing the parent population. Indeed, K00 found FWHM(H$`\beta `$)$`L^{0.27}`$ while in the samples under study the correlation is much flatter. The FWHM–luminosity dependence in various samples will be addressed in a separate paper (Corbett et al. 2003).

The tight $`M_{\mathrm{BH}}L`$ relationship enables the study of the Eddington ratio, $`L/L_{\mathrm{Edd}}`$, in these samples. The observed $`ML`$ relationship suggests a very weak, if any dependence of $`L/L_{\mathrm{Edd}}`$ on $`L`$ or on $`M_{\mathrm{BH}}`$. This impression is confirmed by a formal statistical analysis. Since the results are marginal, they will not be presented here. Another important issue is the mean $`L/L_{\mathrm{Edd}}`$. This depends on the distribution in this property as well as on the exact conversion from $`\lambda L_\lambda `$ to bolometric luminosity and the value of $`\gamma `$. Assuming first $`\gamma =0.58`$ and $`L_{\mathrm{Edd}}=9\lambda L_\lambda (5100)`$, as in K00, gives a median $`L/L_{\mathrm{Edd}}`$ of 0.53. The composite spectra published recently by Telfer et al. (2002) suggest a different conversion with $`L_{\mathrm{Edd}}5\lambda L_\lambda (5100)`$. This translates to a median of 0.28. The above values are transformed to 0.33 and 0.18 for the case of $`\gamma =0.68`$. In both cases the distribution wis wide, covering about a factor 10 in $`L_{\mathrm{Edd}}`$. Thus, the choice of $`\gamma =0.58`$ results in a large mean $`L/L_{\mathrm{Edd}}`$ and a large number of sources with super Eddington luminosities. As explained by Woo and Urry (2002), the implications to the derived $`M_{\mathrm{BH}}L`$ relation are very importance (see §4).

## 4 Discussion: the largest BHs and the most luminous galaxies

The new results presented here suggest that the largest BHs are situated in the most luminous quasars that are, typically, the highest redshift sources. At the extreme end of the distribution we find BH masses of order $`5\times 10^{10}`$ $`M_{}`$ if $`\gamma =0.7`$, and $`1.1\times 10^{10}`$ $`M_{}`$ if $`\gamma =0.5`$. This is greater than obtained so far in large samples. The recent work by Shields et al. (2002) aimed at the calibration of the the \[O iii\] $`\lambda 5007`$ line width as a bulge mass estimator. The method is based on the close agreement between FWHM(\[O iii\] $`\lambda 5007`$) and the stellar velocity dispersion $`\sigma _{}`$ at low luminosity and the K00 mass estimates at higher luminosity. Using this method and $`\gamma =0.5`$ (their Table 2) they find one object with $`M_{\mathrm{BH}}`$ exceeding $`10^{10}`$ $`M_{}`$ and several others approaching this mass. As shown in Fig. 1, the C iv$`\lambda 1549`$ samples includes many more sources with such large masses.

Before addressing the cosmological consequences we note the various factors influencing the $`ML`$ relationship and likely reasons for overestimating $`M_{\mathrm{BH}}`$.

1. The K00 sample covers a limited luminosity range and all mass estimates corresponding to $`\lambda L_\lambda (1350)>10^{46}`$ are necessarily obtained by extrapolation. Since this is the only sample available so far, there is no independent way to verify the largest masses until successful reverberation mappings are obtained for higher luminosity AGN. Moreover, as explained in §2, the slope of the $`R_{\mathrm{BLR}}`$-$`L`$ relationship is uncertain. The slope chosen here ($`\gamma =0.58`$) is close to the middle of the range. Its increase to 0.7 will increase the mass at the high luminosity end by a factor of about 2.5.
2. The largest new mass estimates are based on the measured $`\lambda L_\lambda (1350)`$ which is scaled to the K00 luminosity assuming the same spectral energy distribution (SED) for high and low luminosity AGNs. This assumption has never been tested in large quasar samples. The data for such test are already available (Telfer et al. 2002) but the results are not yet known. Intrinsic reddening, in the quasar host galaxy, is another potential complication related to the inter-calibration of optical and UV luminosities.
3. The FWHM(C iv$`\lambda 1549`$) may not reflect the virial motion of the BLR gas in high luminosity quasars.
4. The samples used here suffer from various selection effects. This influence only slightly the largest derived masses but can affect much more the $`LM_{\mathrm{BH}}`$ correlation. For example, magnitude limited samples may not include the less luminous quasars, those with the smallest $`L/L_{\mathrm{Edd}}`$. This results in a false impression of a very strong $`ML`$ correlation. Woo and Urry (2002) investigated this idea in great detail and concluded that all strong $`ML`$ correlations obtained so far suffer from such a selection effect.

The main conclusion of this work is that the largest BH masses are found in the highest luminosity quasars. The masses of such BHs can reach the extreme values of $`10^{10.310.6}`$ $`M_{}`$, depending on the value of $`\gamma `$. Using recent conversions to host galaxy properties one finds $`M_{\mathrm{bulge}}10^{13.113.4}`$ $`M_{}`$ (Kormendy & Gebhardt 2001), $`M_{\mathrm{B},\mathrm{bulge}}25`$ mag (Kormendy 2001) and $`\sigma _{}`$ exceeding 800 km/sec (Tremaine et al. 2002). Such galaxies have never been observed and are not predicted to exist by standard galaxy formation theories.

In principle, this is still consistent with the observations since the sources with the largest $`M_{\mathrm{BH}}`$ are the most luminous ones and will completely out-shine any host galaxy. Thus, there is no direct way to rule out the existence of such galaxies. However, the theoretical implication are in conflict with recent ideas that the largest galaxies attain their mass through a series of mergers, a process that operates continuously to redshift 2 or smaller. A similar difficulty is found for the BH growth since those same theories (e.g. Haehnelt and Kauffman 2000; Yu and Tremaine 2002) assume that galactic nuclei BHs increase their mass up to redshifts smaller than 2 by the same series of mergers (or, perhaps, only through large mergers). Thus, the largest BHs are predicted to be associated with the most massive galaxies at $`z<2`$, in conflict with the data in Fig. 2. It is clear that active BH with $`M_{\mathrm{BH}}>10^{10}`$ $`M_{}`$ are not found in the local universe. It is also clear that dormant BHs of this mass, or the galaxies with extreme properties that are supposed to host such BHs, have never been found. The whereabout of the huge BH formed at $`z3`$ is thus unknown.

A more plausible suggestion is that some or all the conversion factors used to obtain the galactic mass, magnitude and $`\sigma _{}`$ from the BH mass, that are based on measurements in the local universe, cannot be extrapolated to high luminosity high redshift objects. Perhaps they are only valid at $`z<2`$, after galaxies and nuclear BHs have accumulated most of their mass. If correct, this would mean that some “normal looking” galaxies contain extremely massive BHs. A similar suggestion by Laor (2001) involves a dependence of $`M_{\mathrm{BH}}/M_{\mathrm{bulge}}`$ on the BH mass or the absolute magnitude of the host galaxy.

To conclude, either the measurements of BH masses presented here for the most luminous quasars are grossly overestimated, because of the reasons described above, or else the relationships between BH masses and various properties of their host galaxies at high $`z`$ are very different from those measured in the local universe. A second conclusion, which is less certain because of various selection effects, is that for AGNs, $`M_{\mathrm{BH}}L^{0.9\pm 0.15}`$.

The work described in this paper is based primarily on a decade-long AGN monitoring project at the Wise observatory. I am grateful to many of my colleagues and students that helped in making this into a very successful project. Special thanks go to Dan Maoz, Shai Kaspi and Ohad Shemmer who led various parts of the project and without whom it would have been impossible to bring it to completion. Useful discussions with Ari Laor are gratefully acknowledged. This work is supported by the Israel Science Foundation grant 545/00.


text:
# Asymmetric binaries and hierarchy of states in self-gravitating systems

## Abstract

We complete the analysis of the equilibrium shapes of rotating self-gravitating gases initiated in Phys. Rev. Lett. 89 031101 (2002) by studying the formation and thermodynamics of asymmetric binaries and by introducing a family of order parameters to classify the possible equilibrium states.

Gravity-dominated systems such as classical self-gravitating gases are known to possess highly non-standard thermodynamic properties binney&tremaine87 ; padmanabhan90 ; padmata . The long-range nature of the Newtonian potential prevents thermodynamic functionals from being extensive gallavotti99 , causing such features as the inequivalence of statistical ensembles, the occurrence of negative specific heats and of collapse transitions, and the existence inhomogeneous ground states lyndenbell67 ; lyndenbell68 ; thirring70 ; devega96 ; chavanis01 . These attributes are not peculiar to gravitating gases alone. Similar traits mark the behaviour of systems with long-range, non-integrable potentials barre01 , of finite systems or of systems of a size comparable to the interaction range, such as atomic clusters or highly excited (“hot”) nuclei gross82 ; gross140 , are even present in certain short-ranged models compagner  and in general at phase separations gross174 . Besides the relevance to astrophysics and cosmology strictly speaking, understanding the equilibrium state of self-gravitating systems is hence important also from a broader, fundamental viewpoint gross174 .

In a recent work prl ; epjb , to be later referred to as I, we investigated the possible equilibrium shapes of three-dimensional rotating self-gravitating gases using a microcanonical mean-field theory in which the system is enclosed in a finite spherical volume $`V`$ with the total energy $`E`$ and the total angular momentum $`𝑳`$ conserved. The singularity of the Newtonian potential at short distances was cured by the use of Lynden-Bell statistics lyndenbell67 . It turned out that, along with the known “gaseous” and “collapsed” (a single, dense cluster) states, the system can form symmetric “binaries” (two identical clusters) if $`E`$ is sufficiently low, that is if gravity is strong enough, and rotation is sufficiently fast. Moreover, at intermediate $`E`$ and $`𝑳`$, a large phase-coexistence region with negative specific heat appears, where the system is found in dense structures (either a single cluster or a binary) embedded in a vapor.

The purpose of this paper is to complete the study of I by describing the formation and thermodynamics of asymmetric binaries (which are far more common in the universe than symmetric ones). Furthermore, a family of order parameters will be introduced to discern the possible equilibria and allow for a classification of the equilibrium states. Some “exotic” equilibria (with a very low entropy) that can be obtained using these order parameters are presented as an example.

We begin by recalling the key steps in the theory of I. The Hamiltonian is ($`i,j=1,\mathrm{},N`$)

$$H=\frac{1}{2m}\underset{i}{}p_i^2Gm^2\underset{i<j}{}|𝒓_i𝒓_j|^1$$
(1)
with $`𝒓_iV^3`$ denoting the position and $`𝒑_i^3`$ the momentum of the $`i`$-th particle, $`m`$ being its mass. $`G`$ is the gravitational constant. The task is to find the particle density profiles that maximize the Boltzmann entropy $`S=k\mathrm{log}W`$, where the state sum $`W`$ is given by

$$W=\frac{\alpha }{N!}\delta (HE)\delta (𝑳\underset{i}{}𝒓_i\times 𝒑_i)d^{3N}𝒓d^{3N}𝒑$$
(2)
with $`\alpha `$ a constant. $`W`$ is a function of $`E`$ and $`𝑳`$, as are $`S`$ and the desired density profiles. The radius $`R`$ of the box is taken to be fixed. Adhering to the notation of I, we shift to the reduced variable $`𝒙=𝒓/R`$, which we shall indicate as $`𝒙=(x_1,x_2,x_3)`$ (Cartesian representation) or $`𝒙=(x,\theta ,\varphi )`$ (spherical representation). Particles are assumed to have “hard cores”, so that the density profile $`c`$ is normalized as $`c(𝒙)𝑑𝒙=\mathrm{\Theta }`$, with $`0c1`$. $`\mathrm{\Theta }`$ is a fixed parameter ($`0\mathrm{\Theta }4\pi /3`$) which throughout this study, as in I, is taken to be $`\mathrm{\Theta }=0.02`$. A brief discussion about the dependence of the results on $`\mathrm{\Theta }`$ is given at the end of the paper. Finally, we take $`m=1`$ and measure energy and angular momentum in units of $`GN^2/R`$ and $`\sqrt{RGN^3}`$, respectively.

Assuming Lynden-Bell statistics to mimic the hard cores lyndenbell67 , at stationary points of the entropy surface the density profile $`c`$ satisfies the condition (see I for details)

$$\mathrm{log}\frac{c(𝒙)}{1c(𝒙)}=\frac{\beta }{\mathrm{\Theta }}\frac{c(𝒙^{\mathbf{}})}{|𝒙𝒙^{\mathbf{}}|}𝑑𝒙^{\mathbf{}}+\frac{1}{2}\beta (𝝎\times 𝒙)^2\mu $$
(3)
where

$$\beta =\frac{3/2}{[E\frac{1}{2}𝑳^T𝕀^1𝑳\mathrm{\Phi }[c]]},$$
(4)
$`𝝎=𝕀^1𝑳`$ is the angular velocity, and $`\mu `$ is a Lagrange multiplier implementing the constraint on $`\mathrm{\Theta }`$. $`𝕀`$ stands for the inertia tensor, with elements ($`a,b=1,2,3`$)

$$I_{ab}[c]=\frac{1}{\mathrm{\Theta }}c(𝒙)(x^2\delta _{ab}x_ax_b)𝑑𝒙$$
(5)
(in units of $`NR^2`$), while $`\mathrm{\Phi }`$ is the gravitational potential. In general, the solutions of Eq. (3), whose number and form depend on $`E`$ and $`𝑳`$, can be written in spherical coordinates as

$$c(𝒙)=\underset{l=0}{\overset{\mathrm{}}{}}\underset{m=l}{\overset{l}{}}b_{lm}(x)Y_{lm}(\theta ,\varphi )$$
(6)
where the “weights” $`b_{lm}(x)`$ have to be determined by solving (numerically) the system of integral equations

$$b_{lm}(x)=g(x,\theta ,\varphi )Y_{lm}(\theta ,\varphi )d\mathrm{cos}\theta d\varphi $$
(7)
$$g(x,\theta ,\varphi )=\left[1+e^{\frac{\beta }{\mathrm{\Theta }}_{l,m}u_{lm}(x)Y_{lm}(\theta ,\varphi )\frac{1}{2}\beta \omega ^2x^2\mathrm{sin}^2\theta +\mu }\right]^1$$

Some details of this task are discussed in I.

If $`𝑳`$ is taken to lie along the $`3`$-axis and if only even harmonics ($`l=0,2,4,\mathrm{}`$) are included in the series (6), entropy-maximizing solutions of (3) can be divided in two main classes (see I): (a) axial-rotationally symmetric ones, e.g. gas-cloud (high $`E`$, low $`L`$), distorted gas clouds or disks (high $`E`$, high $`L`$), and single clusters (low $`E`$, low $`L`$); (b) axial-rotationally asymmetric ones, like identical double clusters or symmetric binaries (SBs, low $`E`$, high $`L`$). The former are independent of the $`\varphi `$ angle, and thus for them the series (6) involves only the *zonal* harmonics, with $`m=0`$. The latter, instead, must depend on $`\varphi `$ and the lowest-order term in (6) that can provide this dependence is the *sectoral* harmonics $`Y_{2,2}`$. In fact, on physical grounds, SBs can be characterized by a non-zero value of the order parameter

$$D_{2,2}=\left|x^4b_{2,2}(x)𝑑x\right|$$
(8)
which is proportional to the difference between the diagonal components of the inertia tensor in the $`1`$ and $`2`$ directions. Using (5) and (6) it is indeed easy to see that

$$D_{2,2}=\mathrm{\Theta }\sqrt{\frac{15}{16\pi }}|I_{11}I_{22}|$$
(9)
The inclusion of odd harmonics leads to new solutions, namely asymmetric binaries (ABs). A sample of asymmetric binaries obtained in this way (with the addition of fixing the center of mass at $`𝒙=\mathrm{𝟎}`$) is shown in Fig. 1.

As it was to be expected, such solutions appear at low energies at intermediate angular momenta ($`L0.46`$) and exist in a small range of values of $`L`$ up to $`L0.5`$. This whole region lies in the mixed phase with negative specific heat of I, hence asymmetric binaries at equilibrium are embedded in a vapor and they don’t exist as a “pure” phase. Moreover, in this range they are significantly more probable than symmetric binaries, as shown in Fig. 2, where the entropy of the different stable solutions at $`L=0.46`$ is shown.

A simple physical argument suffices to explain why at high angular momenta ABs must have a much lower entropy than SBs, if they exist at all. SBs have a larger moment of inertia than ABs, hence a smaller rotational energy. As a consequence, for SBs the “random” kinetic contribution $`EL^2/(2I)`$ is larger than for ABs. Thus they have a larger entropy. So ABs are expected to be strongly suppressed or significantly less probable at high angular momenta. Indeed, we were unable to find stable asymmetric binaries for $`L>0.5`$. The fact that ABs tend to form the smaller cluster in proximity of the boundary of the box is due to the fact that ABs start to form by evaporating small masses from a large, dense cluster. The conservation of the center of mass relegates them to the border. In well formed binaries, the peak of density in the smaller cluster is separated from the boundary.

In Fig. 3 the caloric curve $`\beta =S/E`$ versus $`E`$ is shown at $`L=0.46`$.

The jump occurring at $`E=0.2`$, where ABs become more favorable than SBs, is an artifact of the asymptotics employed in the mean-field theory of I. (Being $`S`$ a smooth and multiply-differentiable function, such a feature is clearly spurious.)

Again on physical grounds, the natural order parameter for AB solutions is $`x_1^3:=x_1^3c(𝒙)𝑑𝒙`$, because in spatially asymmetric structures it must be non-zero because of evident symmetry reasons (the two clusters are aligned along the $`x_1`$ axis), at odds with what must happen in symmetric structures. Now using (6) one can see by working out explicitly the integrals involving spherical harmonics morse  that

$$\begin{array}{c}x_1^3=\sqrt{\frac{12\pi }{25}}x^5b_{1,1}(x)𝑑x\sqrt{\frac{6\pi }{175}}x^5b_{3,1}(x)𝑑x+\hfill \\ \hfill +\sqrt{\frac{2\pi }{35}}x^5b_{3,3}(x)𝑑x\end{array}$$
(10)
In principle, all of the above terms should be considered. However, in any spatially symmetric structure (gas, single cluster, SBs) all of them are zero, because odd harmonics do not contribute, as discussed above. Each of them is non-zero only in presence of spatial asymmetries, as in ABs, hence any of them can be used as an order parameter. In particular, we find it convenient to choose the quantity

$$D_{3,1}=\left|x^5b_{3,1}(x)𝑑x\right|$$
(11)
whose behaviour at $`L=0.46`$ is shown in Fig. 4. AB solutions appear at $`E=0.18`$, where the order parameter jumps discontinuously from zero to a finite value, and becomes more favorable than SBs at $`E=0.2`$.

Eq.s (8) and (11) suggest that, in general, one can classify the different equilibrium states by introducing the family of order parameters defined by

$$D_{l,m}=\left|x^{l+2}b_{l,m}(x)𝑑x\right|$$
(12)
with integers $`l=0,1,2,\mathrm{}`$ and $`m=l,l+1,\mathrm{},l`$. For small $`l`$, each of these functions has a simple physical interpretation.

For $`l=0`$, $`D_{0,0}`$ is simply a constant proportional to $`\mathrm{\Theta }`$, since, recalling that $`c(𝒙)𝑑𝒙1=\mathrm{\Theta }`$,

$$c(𝒙)𝑑𝒙=2\sqrt{\pi }x^2b_{0,0}(x)𝑑x$$
(13)
For $`l=1`$, one has that $`D_{1,1}`$, $`D_{1,0}`$ and $`D_{1,1}`$ are nothing but the coordinates of the center of mass (modulo a proportionality constant):

$$x_1^{CM}x_1=\sqrt{\frac{4\pi }{3}}x^3b_{1,1}(x)𝑑x$$

$$x_2^{CM}x_2=\sqrt{\frac{4\pi }{3}}x^3b_{1,1}(x)𝑑x$$
(14)
$$x_3^{CM}x_3=\sqrt{\frac{4\pi }{3}}x^3b_{1,0}(x)𝑑x$$

For $`l=2`$, the situation is slightly more involved. $`D_{2,0}`$ can be easily seen to be linked to the quantity

$$x^23x_3^2=\sqrt{\frac{16\pi }{5}}x^4b_{2,0}(x)𝑑x$$
(15)
which is a convenient order parameter for disk (deformed gas) solutions. In fact, for an hypothetical configuration being a perfect sphere (say a homogeneous gas cloud filling the whole $`V`$), one would have $`D_{2,0}=0`$ evidently. When rotation deforms the cloud by shrinking it along the $`3`$-axis, one would have $`D_{2,0}0`$. As discussed in I, $`D_{2,2}`$ is instead the order parameter discerning axial-rotationally symmetric structures from axial-rotationally asymmetric ones (it is non zero for both SBs and ABs). It is particularly convenient for SBs because of its connection to the inertia tensor components, see (9). As for the remaining $`D_{2,m}`$’s, they can be easily seen to be related to the off-diagonal components of the inertia tensor (5).

It is clear that higher-order $`D`$’s are all connected to homogeneous polynomials in $`x_1`$, $`x_2`$, $`x_3`$, which are in turn connected to the Cartesian representation of spherical harmonics morse . For example, for $`l=3`$, $`D_{3,0}`$ is connected to

$$2x_3^33x_1^2x_33x_2^2x_3=\sqrt{\frac{16\pi }{7}}x^5b_{3,0}(x)𝑑x$$
(16)
Their physical interpretation for larger $`l`$ is however less straightforward.

Using these order parameters, solutions can be classified hierarchically, a procedure that is particularly effective at high $`L`$. In our setting, with the center of mass fixed at $`𝒙=\mathrm{𝟎}`$ one has the following picture. Homogeneous gas clouds have all $`D_{lm}`$’s equal to zero except $`D_{0,0}`$. For deformed clouds (disks) only $`D_{l,0}`$ with $`l`$ even are non zero. For symmetric binaries (which are connected to “gas”-like solutions continuously, see I), also $`D_{2,2}0`$ (and in general all $`D_{l,m}`$ with $`l`$ and $`m`$ even), all other $`D`$’s being zero, while for asymmetric binaries $`D_{3,1}`$ is also not zero. This classification can be continued further, for the order parameters (4) allow to construct all possible solutions of (3), by adding a proper field that selects the solution for which the chosen order parameter is not zero. As an example, we show in Fig. 5 two of several equilibria that can be obtained in this way at very high angular momentum. Of course, such “exotic” solutions of (3), if stable, have much lower entropy than those discussed in I and here.

To summarize, we have carried further the equilibrium study of rotating self-gravitating gases of classical non-overlapping particles started in I by describing the formation and thermodynamics of asymmetric binaries. Moreover, we introduced a family of order parameters that allow to characterize all the different shapes that can be encountered (homogeneous gas, disk, symmetric binary, asymmetric binary, along with more exotic states), and discussed their physical meaning in simple cases. The present study does not substantially alter the global phase diagram presented in I, since asymmetric binaries do not exist as a pure thermodynamic state but only in the phase coexistence region. Among the problems that remain open we would like to mention the following two. First, the dependence on $`\mathrm{\Theta }`$. Our results have been obtained at fixed $`\mathrm{\Theta }=0.02`$. Different $`\mathrm{\Theta }`$’s may lead to different results. By letting $`\mathrm{\Theta }0`$ one obtains the case of overlapping particles (hard cores are neglected and Lynden-Bell statistics is replaced by Boltzmann statistics in space), which is a classical problem in cosmology dating back to antonov62 . The rotating case is dealt with in npb . For Boltzmann particles we haven’t found any evidence of asymmetric binaries. On the other hand, when $`\mathrm{\Theta }`$ grows ($`\mathrm{\Theta }4\pi /3`$) the system may become too closely packed and the interesting features one observes for small $`\mathrm{\Theta }`$ might disappear at some point. The second problem is the dependence on $`R`$. Throughout this paper we have assumed that $`R`$ is fixed but it would be interesting to analyze this problem in a context in which $`R`$ varies, for instance to mimic an expanding-universe scenario.


text:
# Signatures of AGN model for UHECR

## Abstract

We demonstrate that the energy spectra of Ultra High Energy Cosmic rays (UHECR) as observed by AGASA, Fly’s Eye, HiRes and Yakutsk detectors, have the imprints of UHE proton interaction with the CMB radiation in the form of the dip at $`E1\times 10^{19}`$ eV, of the beginning of the GZK cutoff, and of very good agreement with calculated spectrum shape. We argue that these data, combined with small-angle clustering and correlation with AGN (BL Lacs), point to the AGN model of UHECR origin at energies $`E<1\times 10^{20}`$ eV. Our consideration includes also the case when correlation with BL Lacs is excluded from the analysis. The excess of the events at $`E>1\times 10^{20}`$ eV , which is observed by AGASA (but absent in HiRes data) can be explained by another component of UHECR, e.g. by UHECR from superheavy dark matter.

The nature of signal carriers and of the sources of UHECR are not yet established (for recent reviews see ). The most natural primary particles are extragalactic protons. Due to interaction with the CMB radiation the Ultra High Energy (UHE) protons from extragalactic sources are predicted to have a sharp steepening of energy spectrum, so called GZK cutoff . For uniformly distributed sources, the GZK cutoff is characterized by energy $`E_{1/2}`$ where the integral spectrum calculated with energy losses taken into account becomes twice lower than the power-law extrapolation from low energies, $`E_{1/2}=5.3\times 10^{19}`$ eV .

The particles with energies higher than $`1\times 10^{20}`$ eV are undoubtly observed. There are at least two “golden” events at energies $`23\times 10^{20}`$ eV  with very reliable energy determination (see also discussion in Ref.). However, the real contradiction with the existence of the GZK cutoff is observed only by AGASA (see Fig.1). Data of HiRes  are in good agreement with presence of the GZK cutoff (Fig. 1). The data of other detectors are not as conclusive, though a few events with energy higher than $`1\times 10^{20}`$ eV are observed there (see ) for a review and  for recent discussion).

In this Letter we shall demonstrate that the observed spectra have the imprints of UHE proton interaction with the CMB radiation in the form of the dip at the energy $`E1\times 10^{19}`$ eV, produced by $`p+\gamma _{\mathrm{CMB}}p+e^++e^{}`$ interaction, in the form of the beginning of the GZK cutoff, and in the form of good agreement between predicted and observed spectra. We argue that at least at energies up to $`1\times 10^{20}`$ eV the data (spectrum, small-scale clustering , and probably the correlation with BL Lac sources  can be explained in the model with AGN as the sources of UHE protons.

Calculating the spectra, we shall use the cosmological parameters as follows from recent observations : flat universe with $`\mathrm{\Omega }_{\mathrm{tot}}=1`$ and $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$. At small redshifts $`z`$, neutrinos, baryons and CDM behave as non-relativistic matter with $`\mathrm{\Omega }_m=0.3`$. We shall use the Hubble constant $`H_0=70`$ km/s Mpc. The relation between time and redshift is given by

$$dt=\frac{dz}{H_0(1+z)\sqrt{\mathrm{\Omega }_m(1+z)^3+\mathrm{\Omega }_\mathrm{\Lambda }}},$$
(1)
The spectrum of UHE protons in the model with uniform distribution of the sources and with the power-law generation spectrum can be calculated using the formalism of Ref., with the continuous energy losses from Refs. (note the difference in formulae due to cosmology with $`\mathrm{\Lambda }`$ term):

$`J_p(E)`$ $`=`$ $`(\gamma _g2){\displaystyle \frac{c}{4\pi }}{\displaystyle \frac{_0}{H_0}}{\displaystyle _0^{z_m}}{\displaystyle \frac{dz_g(1+z_g)^{m1}}{\sqrt{\mathrm{\Omega }_m(1+z_g)^3+\mathrm{\Omega }_\mathrm{\Lambda }}}}`$ (2)
$`\times `$ $`\left[E_g(E,z_g)\right]^{\gamma _g}{\displaystyle \frac{dE_g(z_g)}{dE}},`$ (3)
where $`z_m`$ is a maximum redshift in the evolution of the sources, $`z_g`$ is a redshift at generation and $`E_g(z_g)`$ is energy of a proton at generation, for present ($`z=0`$) energy $`E`$; $`_0=n_0L_p`$ is cosmic ray (CR) emissivity at $`z=0`$ ($`n_0`$ and $`L_p`$ are space density of the sources and their CR luminosity, respectively). As the general case we assume cosmological evolution of the sources given by $`(z)=_0(1+z)^m`$, where the absence of evolution corresponds to $`m=0`$. All energies in Eq.(3) are given in GeV and luminosities in GeV/s. Dilation of energy interval is given by , modified by Eq.(1) as

$`{\displaystyle \frac{dE_g(z_g)}{dE}}=(1+z_g)`$ (4)
$`\times \mathrm{exp}\left[{\displaystyle \frac{1}{H_0}}{\displaystyle _0^{z_g}}{\displaystyle \frac{dz(1+z)^2}{\sqrt{\mathrm{\Omega }_m(1+z)^3+\mathrm{\Omega }_\mathrm{\Lambda }}}}\left({\displaystyle \frac{db_0(E^{})}{dE^{}}}\right)\right]`$ (5)
where $`b_0(E)=dE/dt`$ is the energy loss due to interaction with CMB photons at $`z=0`$. Derivative $`db_0(E^{})/dE^{}`$ at $`z=0`$ ( given in Ref.) is taken at energy $`E^{}=(1+z)E_g(E,z)`$.

For particles with energies $`E1\times 10^{17}`$ eV, the maximum redshift for evolution of CR sources $`z_m`$ is not important if it is larger than 4. In Fig.1 the calculated spectra are compared with the data of AGASA, HiRes, Fly’s Eye and Yakutsk detectors. Note that the theoretical spectra are the same in all four panels. We assume the generation spectrum $`E^2`$ at $`EE_c`$ and $`E^{\gamma _g}`$ at $`EE_c`$ with a spectrum cutoff at $`E_{\mathrm{max}}`$. The calculations for this case are easily generalized from Eq.(3), as it is done by Eq.(10) from Ref.. As parameters we have chosen $`E_{\mathrm{max}}=1\times 10^{21}`$ eV, $`E_c=1\times 10^{18}`$ eV, $`\gamma _g=2.7`$ for non-evolutionary ($`m=0`$) model and $`\gamma _g=2.5`$ for evolutionary ($`m=3`$) model. The required CR emissivity for both models is $`_0(2.53.5)\times 10^{46}`$ ergs/Mpc3yr. The choice of $`E_c`$ is motivated by the observed chemical composition. We assume that transition to extragalactic component occurs at the observed second knee -, at $`E_24\times 10^{17}`$ eV. The observed rigidity cutoff for protons, $`E/Z=(45)\times 10^{15}`$ eV according to KASCADE data , corresponds to the cutoff energy for iron nuclei $`E1\times 10^{17}`$ eV. Remaining gap $`(14)\times 10^{17}`$ eV can be filled by ultraheavy nuclei with charge up to $`Z=92`$ . The transition to the lighter chemical composition with a large fraction of protons has been observed at $`E>3\times 10^{17}`$ eV by AGASA , FE , Yakutsk  and HiRes  detectors. The data of HiRes  show that at $`E>6\times 10^{17}`$ eV the protons can be the dominant component.

Fig.1 shows that that the signatures of interaction of UHE protons with the CMB radiation, the dip and the beginning of the GZK cutoff, are seen in the data.

The most natural sources of the observed UHE protons are AGN. The required CR emissivity meet well the local emissivity of AGN, e.g. that of Seyfert galaxies is of order $`_{\mathrm{Sy}}n_{\mathrm{Sy}}L_{\mathrm{Sy}}1\times 10^{48}`$ ergs Mpc3yr. The protons can be accelerated in AGN up to energies of order $`10^{21}`$ eV . An interesting possibility is acceleration in the jets by unipolar induction . The correlation of UHECR with BL Lacs , i.e. with AGN whose jets are directed towards us, strongly supports this mechanism.

Two sets of observational data favor rectilinear propagation of UHE signal carriers in the universe.

The first set is the above-mentioned correlation  with BL Lacs. The second one is the small-angle clustering . Its most natural interpretation is given  in terms of statistically occasional arrival of two or three particles from a compact source. Such an interpretation needs rectilinear propagation of the primaries and large number of sources .

Since the proton origin of UHECR is almost proved, the correlation with BL Lacs directly implies the rectilinear propagation of UHE protons. Thus, these correlations become supersensitive tools to measure extragalactic magnetic fields. Below we shall discuss what is the scale of this field and whether this field can be already excluded.

Magnetic field must not produce the angular deflection larger than angular resolution of sources in the detectors, which is typically $`\vartheta _{\mathrm{res}}2.5^{}`$. The correlation is found in the energy range $`(48)\times 10^{19}`$ eV, for which the largest attenuation length is $`l_{\mathrm{att}}1000`$ Mpc. The required upper limit for the magnetic field, which is homogeneous on this scale, is $`B_l2\times 10^{12}l_{1000}^1`$ G, where $`l_{1000}`$ is attenuation length for $`4\times 10^{19}`$ eV protons in units of 1000 Mpc. For a magnetic field with small homogeneity length $`l_{\mathrm{hom}}`$ the required upper limit is

$$B\frac{E\vartheta _{\mathrm{res}}}{e\sqrt{l_{\mathrm{att}}l_{\mathrm{hom}}}}6\times 10^{10}\text{G},$$
(6)
where the numerical value is given for $`l_{\mathrm{att}}1000`$ Mpc and $`l_{\mathrm{hom}}10`$ kpc.

We argue that these fields are not excluded.

The observed Faraday rotations give only the upper limits on large scale extragalactic magnetic field . All known mechanisms of generation of the large scale cosmological magnetic field results in extremely weak magnetic field $`10^{17}`$ G or less (for a review see ). The strong magnetic field can be generated in compact sources, presumably by dynamo mechanism, and spread away by the flow of the gas. These objects thus are surrounded by magnetic halos, where magnetic field can be estimated or measured. The strong magnetic fields of order of $`1\mu `$G are indeed observed in galaxies and their halos, in clusters of galaxies and in radiolobes of radiogalaxies. As an example one can consider our local surroundings. Milky Way belongs to the Local Group (LG) entering the Local Supercluster (LS). LG with a size $`1`$Mpc contains 40 dwarf galaxies, two giant spirals (M31 and Milky Way) and two intermediate size galaxies. The galactic winds cannot provide the appreciable magnetic field inside this structure. LS with a size of 10 – 30 Mpc is a young system where dynamo mechanism cannot strengthen the primordial magnetic field. In fact LS is filled by galactic clouds submerged in the voids. The vast majority of the luminous galaxies reside in a small number of clouds: 98 % of all galaxies in 11 clouds . Thus, accepting the hypothesis of generation of magnetic fields in compact sources, one arrives at the perforated picture of the universe, with strong magnetic fields in the compact objects and their halos (magnetic bubbles produced by galactic winds) and with extremely weak magnetic fields outside. However, even in this picture there is a scattering of UHE protons off the magnetic bubbles and the scattering length is $`l_{\mathrm{sc}}1/\pi R^2n`$, where $`R`$ is the radius ofa magnetic bubble and $`n`$ is their space density. Among different structures, the largest contribution is given by galaxy clusters which can provide $`l_{\mathrm{sc}}(12)\times 10^3`$ Mpc.

Leaving the correlation of UHECR with AGN (BL Lacs) as an open problem for future observations we turn now to the alternative of excluding these correlations from analysis. The small-angle clustering then can be probably explained in the other extreme case of very strong magnetic field due to lensing effect . Rectilinear propagation of UHE protons is not needed anymore.

Till now we discussed mostly UHECR at $`E<1\times 10^{20}`$ eV. At higher energies, as AGASA data show, there might be an excess of events. They should be interpreted as the new component. One possibility is given by Superheavy Dark Matter (SDMP) . The spectrum of UHE particles, produced at the SHDM decay is now reliably calculated using the different methods . All calculations give very similar spectrum with $`\gamma _g1.92.0`$. This spectrum is shown in Fig.1. It describes well the observed AGASA excess. The primary particles are predicted to be UHE photons, which at energy $`E>1\times 10^{20}`$ eV cannot be excluded by available experimental data, in particular by the inclined Haverah Park showers

In conclusion, the observed energy spectra reveal the signatures of interaction of UHE protons with the CMB in the form of the dip, beginning of the GZK cutoff and the good agreement with the predicted spectrum. Combined with small-angle clustering  and correlation with BL Lacs , these data require the rectilinear propagation of UHE protons and AGN as their sources. The correlation with AGN (BL Lacs) becomes thus the most sensitive method of measuring very weak extragalactic magnetic fields. In case this correlation is not confirmed by future observational data, and the small-angle clustering is explained by some more sophisticated phenomena (e.g. by magnetic lensing), AGN remain the favorable UHECR sources, but as one of the possible candidates.

This work was supported in part by INTAS (grant No. 99-01065). We are grateful to Pierre Sokolsky for active discussion and for providing us with HiRes data. Many thanks to Masahiro Teshima for the data of AGASA, and to M.I.Pravdin for the data of Yakutsk. We gratefully acknowledge the useful discussions with Michael Kachelriess and Igor Tkachev.


text:
# |𝑉_{𝑐𝑏}| and |𝑉_{𝑢𝑏}| Determinations from Inclusive Semileptonic 𝐵 Decays

## 1 INTRODUCTION

$`|V_{cb}|`$ and $`|V_{ub}|`$ are two of the fundamental constants of the Standard Model. They need to be determined by experimental measurements. These can be done with leptonic, semileptonic, or nonleptonic decays of $`B`$ mesons. In this talk I will discuss the determinations of $`|V_{cb}|`$ and $`|V_{ub}|`$ from inclusive semileptonic B decays. An observable in these processes is schematically related to $`|V_{cb}|`$ or $`|V_{ub}|`$ in the form

$$\mathrm{observable}=|V_{c(u)b}|^2T,$$

where $`T`$ is a quantity derived from theory. Theory may also be involved in the experimental analysis to obtain the observable. With the theoretical relationship between the observable and $`|V_{cb}|`$ or $`|V_{ub}|`$, the measured quantity can be converted into a value of $`|V_{cb}|`$ or $`|V_{ub}|`$. The main obstacle in the theoretical description of the weak decay processes is long-distance strong interaction effects. There are two classes of observables: (1) routine observables, such as branching fractions and lifetimes, and (2) theory-motivated observables, such as the differential decay rate $`d\mathrm{\Gamma }/dw`$ at zero recoil ($`w=1`$) in the exclusive semileptonic decay $`BD^{}l\nu `$ for determining $`|V_{cb}|`$. The latter is motivated from theory because its relationship with $`|V_{cb}|`$ or $`|V_{ub}|`$ is theoretically clean and model-independent in some sense.

## 2 THEORY OF INCLUSIVE SEMILEPTONIC $`B`$ DECAYS

Strong interaction effects on inclusive semileptonic $`B`$ decays are contained in the hadronic tensor

$`W_{\mu \nu }`$ $`=`$ $`{\displaystyle \frac{1}{2\pi }}{\displaystyle d^4ye^{iqy}B\left|[j_\mu (y),j_\nu ^{}(0)]\right|B}`$ (1)
$`=`$ $`g_{\mu \nu }W_1+{\displaystyle \frac{P_\mu P_\nu }{M_B^2}}W_2i\epsilon _{\mu \nu \alpha \beta }{\displaystyle \frac{P^\alpha q^\beta }{M_B^2}}W_3`$
$`+{\displaystyle \frac{q_\mu q_\nu }{M_B^2}}W_4+{\displaystyle \frac{P_\mu q_\nu +q_\mu P_\nu }{M_B^2}}W_5.`$
where the charged current $`j_\mu (y)=\overline{q}(y)\gamma _\mu (1\gamma _5)b(y)`$ with $`q=c,u`$. $`P`$ and $`q`$ denote the momenta of the $`B`$ meson and the virtual $`W`$ boson, respectively. The central issue in the theory of inclusive semileptonic $`B`$ decays is to compute the structure functions $`W_{15}`$, which are a priori independent. In the last several years, theoretical tools have been developed to address the issue.

### 2.1 Light-cone (LC) approach

The starting point of the LC approach  to inclusive B decays is the light-cone expansion. Since the $`B`$ meson is heavy, $`M_B\mathrm{\Lambda }_{\mathrm{QCD}}`$, light-cone separations between the currents dominate the decay dynamics. In analog to the analysis of deep inelastic scattering, the light-cone expansion of the matrix element in the neighborhood of the light-cone, $`y^2=0`$, provides a systematic way of calculating nonperturbative QCD effects:

$$B|[j^\mu (y),j^\nu (0)]|B=\underset{n}{}C_n^{\mu \nu }(y,P)(y^2)^n.$$
(2)
At leading twist, the structure functions are related to a single universal distribution function:

$`W_1`$ $`=`$ $`2[f(\xi _+)+f(\xi _{})],`$ (3)
$`W_2`$ $`=`$ $`{\displaystyle \frac{8}{\xi _+\xi _{}}}[\xi _+f(\xi _+)\xi _{}f(\xi _{})],`$ (4)
$`W_3`$ $`=`$ $`{\displaystyle \frac{4}{\xi _+\xi _{}}}[f(\xi _+)f(\xi _{})],`$ (5)
$`W_4`$ $`=`$ $`0,`$ (6)
$`W_5`$ $`=`$ $`W_3,`$ (7)
where $`\xi _\pm =(q^0\pm \sqrt{|\stackrel{}{q}|^2+m_q^2})/M_B`$. The distribution function is defined as

$`f(\xi )`$ $`=`$ $`{\displaystyle \frac{1}{4\pi }}{\displaystyle \frac{d(yP)}{yP}e^{i\xi yP}}`$ (8)
$`\times B|\overline{b}(0)y\gamma U(0,y)b(y)|B|_{y^2=0}.`$
It gives the probability of finding a $`b`$-quark with momentum $`\xi P`$ inside the $`B`$ meson with momentum $`P`$ .

The leading nonperturbative QCD effect is encoded in the distribution function. What is known about the distribution function? First, the normalization of it is exactly known to be $`_0^1𝑑\xi f(\xi )=1`$, because of the $`b`$-quark number conservation in strong interactions. Second, the gross shape of it is determined. In the free quark limit, $`f(\xi )=\delta (\xi m_b/M_B)`$. The mean value and variance of the distribution function have been calculated using the heavy quark effective theory. Consequently, the distribution function is known to be sharply peaked around $`m_b/M_B`$ with a narrow width of order $`\mathrm{\Lambda }_{\mathrm{QCD}}/M_B`$.

### 2.2 Comparison with the heavy quark expansion (HQE) approach

Another approach to inclusive $`B`$ decays is the heavy quark expansion approach . The starting point of the HQE approach is the local operator product expansion. It exploits the large b-quark mass to perform a local operator product expansion ($`y0`$) of the time-ordered product of currents:

$$T[j^\mu (y),j^\nu (0)]=\underset{d}{}C_d^{\mu \nu }(y)O_d(0).$$
(9)
The starting point of the HQE approach is different from that of the LC approach. The HQE approach is based on a short-distance ($`y0`$) expansion in local operators of increasing dimension, while the LC approach is based on a non-local light-cone ($`y^20`$) expansion in matrix elements of increasing twist. The light-cone expansion includes not only local contributions ($`y=0`$), but also non-local contributions ($`y0`$). Thus, the LC approach contains more dynamic effects.

In addition, the HQE approach assumes quark-hadron duality and has to use the quark phase space, while the LC approach does not rely on quark-hadron duality and uses the physical hadron phase space. The singularities appear at the endpoints of the lepton energy spectra in both $`BX_ul\nu `$ and $`BX_cl\nu `$ calculated in the HQE approach, implying the breakdown of the operator product expansion near the boundaries of phase space. There are no endpoint singularities in the lepton energy spectra calculated in the LC approach.

The partial resummation of the heavy quark expansion introduces a different distribution function (“shape function”) . It should be noted that the correction to the leading contribution given in terms of the shape function in the HQE approach is of order $`\mathrm{\Lambda }_{\mathrm{QCD}}/m_b`$, while the correction to the leading contribution given in terms of the distribution function in the LC approach is of order $`(\mathrm{\Lambda }_{\mathrm{QCD}}/M_B)^2`$.

### 2.3 A manifestation of quark-hadron duality violation

The physical phase space at the hadron level is larger than the phase space at the quark level. Therefore, calculations assuming quark-hadron duality cannot account for the rate due to the extension of phase space from the quark level to the hadron level. This rate missing is a clear manifestation of the violation of quark-hadron duality. This is the case for the HQE approach.

Since the LC approach does not rely on quark-hadron duality and uses the hadron phase space, a comparison of the decay rates calculated in the LC and HQE approaches can quantitatively demonstrate how large duality is violated. It has been found that the total rate for the inclusive semileptonic decay $`BX_cl\nu `$ ($`BX_ul\nu `$) calculated in the LC approach is about $`14\%`$ ($`12\%`$) larger than the HQE approach . The nonperturbative QCD correction changes sign: It increases the total rate in the LC approach, while it decreases the total rate in the HQE approach. Therefore, there is significant duality violation in the HQE approach. In particular, it cannot include the phase space effect. From the measured inclusive semileptonic branching fractions and the lifetime of the $`B`$ meson, the LC approach would give rise to smaller values of $`|V_{cb}|`$ and $`|V_{ub}|`$ than the HQE approach, as I will discuss in more detail below.

## 3 CONVENTIONAL METHODS

Conventionally, the routine observables, the branching fractions and the lifetime of the $`B`$ meson, have been used to determine $`|V_{cb}|`$ and $`|V_{ub}|`$ from inclusive semileptonic $`B`$ decays. These observables are related to $`|V_{cb}|`$ or $`|V_{ub}|`$ in the form

$`{\displaystyle \frac{(BX_cl\nu )}{\tau _B}}`$ $`=`$ $`|V_{cb}|^2\gamma _c,`$ (10)
$`{\displaystyle \frac{(BX_ul\nu )}{\tau _B}}`$ $`=`$ $`|V_{ub}|^2\gamma _u.`$ (11)
Theory is needed to compute $`\gamma _c`$ and $`\gamma _u`$. Besides parametric uncertainties, additional uncertainties in the theoretical calculations of $`\gamma _c`$ and $`\gamma _u`$ stem from the assumption of quark-hadron duality in the HQE approach and the detailed shape of the distribution function, which is unknown at present, in the LC approach.

The uncertainty due to the assumption of quark-hadron duality is irreducible. The errors usually quoted in the determinations of $`|V_{cb}|`$ and $`|V_{ub}|`$ from inclusive semileptonic $`B`$ decays using the HQE calculations should be interpreted with caution. The significant uncertainty from the assumption of quark-hadron duality has not been included in the errors. In contrast, the uncertainty due to the shape of the distribution function in the LC approach can be reduced. The distribution function is universal: It incorporates bound state effects in inclusive $`B`$ decays, including semileptonic , radiative , and nonleptonic  inclusive $`B`$ decays. Measurements of decay distributions and moments in these processes will impose constraints on the distribution function. Especially, the $`\xi _u`$ spectrum \[$`\xi _u=(q^0+|\stackrel{}{q}|)/M_B`$\] in $`BX_ul\nu `$ and the photon energy spectrum in $`BX_s\gamma `$ are directly proportional to the distribution function. These spectra are most sensitive to the shape of the distribution function. Measurements of these spectra may thus provide the most stringent constraint on the distribution function. Moreover, other nonperturbative QCD methods, such as lattice QCD, may help to determine the shape of the distribution function, too.

Including the uncertainties from the input parameters and the shape of the distribution function, the calculations in the LC approach yield

$$\gamma _c=49\pm 9\mathrm{ps}^1,\gamma _u=76\pm 16\mathrm{ps}^1.$$

These results allow reliable determinations of $`|V_{cb}|`$ and $`|V_{ub}|`$ without the assumption of quark-hadron duality. The accuracies could be further improved, as discussed above.

A particular problem exists in the $`|V_{ub}|`$ determination, namely a very large $`BX_cl\nu `$ background. Various kinematic cuts can be used to suppress the background. The fractions of events with the different cuts are listed in Table 1.

Imposing kinematic cuts causes additional theoretical uncertainties, since the extrapolation in phase space to obtain the desired observable involves theory. The hadronic invariant mass ($`M_X`$) cut is most efficient. It retains a vast majority of phase space, so that only a small extrapolation to the entire phase space is needed. The other kinematic cuts are less efficient and would cause larger theoretical uncertainties because of a considerable dependence on dynamic subtlety.

## 4 A NEW METHOD FOR THE PRECISE DETERMINATION OF $`|V_{ub}|`$

In the light-cone limit, $`y^20`$, the $`b`$-quark number conservation leads to the sum rule

$$S_0^1𝑑\xi _u\frac{1}{\xi _u^5}\frac{d\mathrm{\Gamma }}{d\xi _u}=|V_{ub}|^2\frac{G_F^2M_B^5}{192\pi ^3}$$
(12)
for the charmless inclusive semileptonic decay $`BX_ul\nu `$. The sum rule establishes a theoretically clean relationship between the observable $`S`$ and $`|V_{ub}|`$. The sum rule has the following advantages: (1) it is model-independent (in particular, independent of the distribution function), (2) it has no reliance on quark-hadron duality, (3) no free parameters (such as $`m_b`$) except for $`|V_{ub}|`$ enter the sum rule, and (4) it receives no perturbative QCD corrections. Therefore, by using the sum rule, a more reliable and precise value of $`|V_{ub}|`$ can be determined than by using the conventional method.

Of course, to obtain the theory-motivated observable $`S`$ in $`BX_ul\nu `$ one must deal with the large $`BX_cl\nu `$ background. The strategy appears to be: (1) apply the kinematic cut $`M_X<M_D`$ (or other realistic cuts), (2) measure the weighted spectrum $`\xi _u^5d\mathrm{\Gamma }/d\xi _u`$, (3) extrapolate it to entire phase space to obtain the integral $`S`$, and (4) determine $`|V_{ub}|`$ from the observable $`S`$ by using the theoretically clean relationship (the sum rule) between them.

This method for determining $`|V_{ub}|`$ is analogous to the determination of $`|V_{cb}|`$ from the observable, $`d\mathrm{\Gamma }/dw`$ at zero recoil, in the exclusive semileptonic decay $`BD^{}l\nu `$. Both use a theoretically clean relationship, which is derived from QCD symmetries, between the observable and $`|V_{ub}|`$ or $`|V_{cb}|`$. Both require an extrapolation in phase space to obtain the theory-motivated observable.

There are two sources of theoretical uncertainties in this method for the $`|V_{ub}|`$ determination, one from higher twist corrections to the sum rule and another from the extrapolation to obtain $`S`$. In principle, the light-cone expansion systematically takes into account higher twist corrections. The quantitative estimate using the heavy quark effective theory has shown that the error on $`|V_{ub}|`$ from higher twist corrections to the sum rule amounts to around $`1\%`$ . The error on $`|V_{ub}|`$ due to the shape of the distribution function used for the extrapolation has been assessed at the $`6\%`$ level . With improved knowledge of the shape of the distribution function, as discussed before, one could reduce this error.

## ACKNOWLEDGMENTS

This work was supported by the Australian Research Council.