text:
# On the X–ray time lags in the black hole candidates

## 1 Introduction

Strong variability of the X–ray flux is a generic property of the black hole binaries in the so called hard state (see e.g. van der Klis 1994). While it is usually assumed that X–ray emission is produced via Comptonization of the low frequency photons by an optically thin cloud of hot electrons (e.g. Sunyaev and Truemper 1979) the nature of the strong variations is less well understood. One of the possible clues to the origin of the variability may be connected with the presence of the time lags between the light curves in different energy bands. The time lags were first found in the early observations of Cygnus X-1 (Priedhorsky et al. 1979, Nolan et al. 1981) and confirmed in the numerous subsequent observations of Cygnus X-1 and other black hole candidates. In the vast majority of cases hard photons lag with respect to the soft ones. From observations it was found that (i) the value of the time lag depends on frequency, (ii) time lag approximately logarithmically depends on the energy separation between two bands (e.g. Miyamoto & Kitamoto 1989), (iii) cross-correlation function between two bands peaks at zero lags within $``$1 ms (Maccarone et al. 2000 and references therein). Comprehensive description of the observed Cygnus X-1 variability and time lags are given by Nowak et al., 1999a,b and Pottschmidt et al. 2000. A recent reviews of the observed lags and theoretical models are given by Poutanen 2001a,b.

We discuss below the time lags using RXTE data on Cygnus X-1 as the example. In Section 2 we show that observed energy dependence of the time lags excludes significant contribution of the standard reflected component to the time lags in Cygnus X-1. In Section 3 we introduce simple phenomenological model of the time lags, which naturally explains the logarithmic energy dependence of the time lags. In Section 4 we discuss the possibility of explaining the lags as due to diffusive propagation of the perturbations in the accretion flow. In Section 5 we discuss implications of the results. Section 6 summarizes our finding.

## 2 Time lags due to light crossing time of the reflector

For hard X–ray state of Cygnus X-1 one of the popular models assumes that an optically thick and geometrically thin disk is truncated at some radius from the compact object and most of the observed X–ray emission is coming from the more compact optically thin region located closer to the black hole (e.g. Esin, McClintock and Narayan 1997, Zdziarski et al. 1999). A fraction of the hard X–ray photons emitted by the central region can be intercepted by the optically thick disk and reemitted again. This re-emitted/reprocessed component, also called “reflected” component, has a very distinct spectral shape (e.g. Basko, Sunyaev & Titarchuk 1974, George & Fabian 1991). If matter in the optically thick disk is neutral or weakly ionized than the spectrum of the reflected component is harder (in the standard X–ray band) than the illuminating spectrum and it also contains prominent iron fluorescent $`K_\alpha `$ line. Because of the delay of photons which traveled from the emission region to the reflector and then to the observer hard lags may naturally appear (e.g. Poutanen 2001b).

In the presence of the direct and reflected components the observed spectrum $`S(E,t)`$ at time $`t`$ can be written as:

$`S(E,t)=D(E,t)+R(E,t)`$ (1)
where $`D(E,t)`$ and $`R(E,t)`$ are the time dependent spectra of the direct and reflected components respectively. Assuming that the direct component vary in intensity only and not in shape:

$`D(E,t)=A(t)D(E)`$ (2)
and neglecting the dependence of the reflected component spectral shape on the incident and emission angles the reflected component can be written as:

$`R(E,t)={\displaystyle \frac{\mathrm{\Omega }}{2\pi }}R(E){\displaystyle T(\tau )A(t\tau )𝑑\tau }`$ (3)
where $`R(E)`$ are the energy spectra of the reflected component calculated for the reflector subtending solid angle of $`2\pi `$, $`A(t)`$ is the normalization of the direct component and $`T(t)`$ is the reflected component transfer function, which characterizes the time dependence of the intensity of the reflected component in response for a infinitely short flare of the direct component. In the above expression $`T(t)`$ is normalized such that $`_0^{\mathrm{}}T(t)𝑑t=1`$ and the contribution of the reflected component to the time averaged spectrum is expressed through the standard factor $`\frac{\mathrm{\Omega }}{2\pi }`$, where $`\mathrm{\Omega }`$ is the total solid angle subtended by the reflector.
The expression (1) can be rewritten as:

$`S(E,t)=A(t)D(E)+{\displaystyle \frac{\mathrm{\Omega }}{2\pi }}R(E){\displaystyle T(\tau )A(t\tau )𝑑\tau }`$ (4)
Further assuming that contribution of the reflected component to the total spectrum is small one can write the expression for a frequency resolved spectrum, introduced in Revnivtsev et al. 1999:

$`S(E,f)=|\widehat{S}(E,f)|`$
$`|\widehat{A}(f)|D(E)\times \sqrt{1+2{\displaystyle \frac{\mathrm{\Omega }}{2\pi }}{\displaystyle \frac{R(E)}{D(E)}}Re[\widehat{T}(f)]}`$
$`|\widehat{A}(f)|\times (D(E)+{\displaystyle \frac{\mathrm{\Omega }}{2\pi }}R(E)Re[\widehat{T}(f)])`$ (5)
The equivalent width of the iron line in the frequency resolved spectrum is then obviously:

$`EW(f)EW_t\times Re[\widehat{T}(f)],`$ (6)
where $`EW_t`$ is the equivalent width in the total source spectrum. From the above expression it is clear that the equivalent width in the frequency resolved spectrum should decrease if the transfer function is smooth at a given frequency. E.g. when the light crossing time of the reflector is large the equivalent width should be low at the corresponding time scales. Using these arguments Revnivtsev et al. 1999, Gilfanov et al. 2000 set an upper limit on the effective size of the reflector in Cygnus X-1 of $``$50 $`R_g`$. We note here that the above expression assumes linear relation between variabilities in the reflected and direct components in the form of the convolution (3). If on the contrary the variability of the direct and reflected components are completely incoherent/uncorrelated then the contribution of the reflected component to the frequency resolved spectrum is proportional to square of the solid angle subtended by the reflector $`(\frac{\mathrm{\Omega }}{2\pi })^2`$.

One can further test the hypothesis of the large effective size of the reflector using the time lags. Assuming again that the light curve in any energy band can be expressed through (4) the expected phase lag $`\varphi (E,f)`$ between given energy $`E`$ and much lower energy (where contribution of the reflected component is negligible and the spectrum is solely due to direct component) can be written as:

$`tg(\varphi (E,f))={\displaystyle \frac{Im[\widehat{S}^{}(E,f)\widehat{D}(E,f)]}{Re[\widehat{S}^{}(E,f)\widehat{D}(E,f)]}}`$
$`{\displaystyle \frac{\frac{\mathrm{\Omega }}{2\pi }R(E)Im[\widehat{T}]}{D(E)+\frac{\mathrm{\Omega }}{2\pi }R(E)Re[\widehat{T}]}}{\displaystyle \frac{\mathrm{\Omega }}{2\pi }}{\displaystyle \frac{R(E)}{D(E)}}Im[\widehat{T}(f)]`$ (7)
Therefore in this approximation the frequency dependence of the phase/time lags is due to imaginary part of the Fourier transform of the transfer function. Energy dependence of the lags is simply the ratio of the reflected and direct spectral components. Since $`D(E)`$ is assumed to be smooth function of energy the time/phase lag as the function of energy should contain prominent iron line with the equivalent width comparable with the equivalent width of the line in the reflected component, i.e. much larger equivalent width than in the total spectrum.

In order to assess possible contribution of the reflected component in the observed time lags in Cygnus X-1 we calculated lags in the narrow energy channels. For our analysis we use publicly available data of Rossi X-ray Timing Explorer observations P10238 performed between Mar. 26, 1996 and Mar. 31, 1996 with total exposure time $``$70 ksec (we used only observations with all 5 PCU turned on). We used PCA data in the “Generic Binned” mode, having $`\frac{1}{64}`$sec$``$16 msec time resolution in 64 energy channels covering the whole PCA energy band (B\_16ms\_64M\_0\_249). The lags were calculated for every channel above 3.7 keV, with respect to the count rate in the 2.8-3.7 keV energy band. The resulting energy dependence of time lag at the frequencies of the order of 2.5 Hz is shown in Fig.1111Time lags at other frequencies in the range from $``$0.6 to $``$10 Hz have approximately similar energy dependence (see Fig.3). For comparison we plot in the same figure an expected energy dependence of the lags due to the reflected component, according to equation (7). Normalization of the expected curve is arbitrary. Here we used as $`D(E)`$ the averaged spectrum of the source during this observation. Reflected component $`R(E)`$ was calculated using XSPEC V11 (Arnaud 1996) model pexrav of Magdziarz & Zdziarski 1995 plus a Gaussian at the energy of 6.4 keV. The resulting spectrum has been convolved with the PCA energy response matrix. Note that expected lags were calculated using equation (7), i.e. with respect to the energy range where contribution of the reflected component was assumed to be zero. Therefore expected lags do not go to zero at the energy $`<`$4 keV. Two curves shown in Fig.1 differ in the normalization of the Gaussian added to the reflected continuum. For upper curve we used relative normalization of the Gaussian and reflected continuum obtained from the fit to the total Cygnus X-1 spectrum during observation 10238-01-08-00 (Gilfanov et al., 1999). For the lower curve the normalization of the Gaussian was chosen such that the equivalent width of the line with respect to the reflected continuum is $``$1 keV. The latter value is characteristic for an angle averaged reflection from a neutral matter with a normal abundance of heavy elements.

One can see from Fig.1 that the lags due to the reflected component should contain a prominent feature at the energy $``$6.4 keV. We stress again that this feature should be much more prominent in the lags than in the averaged spectrum, where iron line provides at most 10% relative deviation of the observed spectrum from a power law fit. Clearly no such feature is present in the observed lags. We therefore conclude that the observed lags are predominantly not due to the reflected component.

One can use the same argument in order to constrain geometry of the reflecting region. Indeed lack of the iron line feature in the observed energy dependence of the time lags imply certain constrains on the combination of the reflected component strength and the duration of the time delay due to the finite size of the reflector. In the Fig.3 we show the observed energy dependence of the time lags at three representing frequencies: 0.6, 2.5 and 10 Hz. For comparison we show in the same plot expected time lags for different geometries of the reflector and different strength of the reflected component. Four basic models have been considered:

* A: A point isotropic source above the flat disk with the hole in the middle. The height $`h`$ of the source above the disk plane is set to 10$`R_g`$, the size of the hole $`R_{in}`$ in the disk is set to 50$`R_g`$. Here $`R_g`$ is the gravitational radius for a 10$`M_{}`$ black hole. The factor $`\mathrm{\Omega }/2\pi `$, calculated for this geometry, is $``$ 0.2.
* B: The same geometry as for A, but the strength of the reflected component was enhanced, so that effective parameter $`\mathrm{\Omega }/2\pi `$ is 0.54. This is possible if e.g. the source of the primary emission is anisotropic and more emission is coming towards the disk than to the observer. Particular value of the normalization of the reflected component in this model was taken from the fit to the total Cygnus X-1 spectrum during observation 10238-01-08-00 (Gilfanov et al., 1999).
* C: The same geometry as for A, but for $`R_{in}=10R_g`$.The factor $`\mathrm{\Omega }/2\pi `$, calculated for this geometry, is $``$ 0.7

The transfer functions in the form $`\frac{\mathrm{\Omega }}{2\pi }T(t)`$ for models A and C are shown in Fig.2 for two values of the disk inclination $`\theta `$ ($`0^{}`$ and $`50^{}`$). Model B obviously has the same shape of the transfer function as model A and is different only in normalization. For comparison we also show the transfer function for a concave disk with the dependence of disk height over radius adopted from Shakura and Sunyaev (1973). For the frequency range considered here such concave disks do not differ from the flat disks.

The lags expected in these models (for inclination $`\theta =50^{}`$) are shown in Fig.3. The model curves were calculated for the equivalent width of the Fe $`K_\alpha `$ line with respect to the reflected continuum obtained from the fit of the time averaged source spectrum. As was mentioned above lack of prominent iron feature in the energy dependence of the time lags suggests that lags are largely not due to the reflection. Moreover lags expected in models B and C at 10 Hz must produce clear feature at 6.4 keV in contradiction to the data. Even for model A expected lags can be only marginally “hidden” in the observed lags. There are however many uncertainties and simplifications in the models discussed above. In particular for simplicity we neglected the dependence of the reflected component spectrum on the incident angle. Also the surface of the reflector (upper few Thomson optical depths) may not be smooth and flat, but wrinkled. Furthermore ionization of the reflector surface will modify the spectrum of the reflected component. Finally all relativistic effects associated with reflection from the inner regions of the accretion disk (see e.g. Campana and Stella 1995) have been ignored. Account for all these effects may modify the normalization and shape of the energy dependence of the time lags by a factor of few.

Perhaps the most dubious assumption is the “linear” response of the reflected component to the variations of the illuminating flux. Recent analysis (e.g. Nayakshin and Kallman 2000) shows that shape of the reflected component can vary in a complicated fashion in response to increase of the primary flux. Therefore conservative conclusions of this section are that i) “linear” reflected component, containing prominent fluorescent line at 6.4 keV, is not responsible for the observed lags and ii) some of the source/reflector geometrical models can be definitely excluded.

## 3 Time lags due to variations of the power law index

The spectrum of Cygnus X-1 in the 2–20 keV energy band can be reasonably well approximated by a power law. Deviations, in particular due to the reflected component, are present, but their amplitude relative to the power law component is usually less than 10% in this energy range. We assume below that at any moment of time the spectrum can be represented as a power law:

$`S(E,t)=A(t)E^{\alpha (t)}`$ (8)
where $`A(t)`$ is the normalization of the power law as a function of time and $`\alpha (t)`$ is the photon index, which is also a function time. For simplicity we set $`\alpha (t)=\alpha _0\beta (t)`$, where $`\alpha _0`$ is a photon index of the time averaged spectrum and $`\beta (t)`$ is the time variable part of the photon index. Assuming that $`\beta (t)\alpha _0`$ the spectrum of the source can be rewritten as:

$`S(E,t)=A(t)E^{\alpha _0}(1+\beta (t)lnE)`$ (9)
The phase lag $`\varphi `$ at a given frequency $`f`$ can then be written through the Fourier transforms $`\widehat{S}^{}(E,f)`$ of the source light curves at two different energies $`E_1`$ and $`E_2`$ as:

$`sin\varphi ={\displaystyle \frac{Im\left[\widehat{S}^{}(E_1,f)\widehat{S}(E_2,f)\right]}{|S(E_1,f)||S(E_2,f)|}}`$
$`{\displaystyle \frac{Im\left[(\widehat{A}^{}+\widehat{(\beta A)}^{}lnE_1)(\widehat{A}+\widehat{(\beta A)}lnE_2)\right]}{\widehat{A}^{}\widehat{A}}}`$
$`{\displaystyle \frac{Im\left[\widehat{A}^{}\widehat{(\beta A)}\right]}{\widehat{A}^{}\widehat{A}}}\times ln{\displaystyle \frac{E_2}{E_1}}`$ (10)
In the above expression we neglected terms of the order of $`O(\beta ^2)`$. Thus the assumption that the spectrum has a power law shape at any moment of time and the variations of slope are small automatically implies that phase lags (and consequently time lags) scales as $`ln\frac{E_2}{E_1}`$.

### 3.1 Universal slope flux correlation?

It is clear that if variations of slope $`\beta (t)`$ and normalization $`A(t)`$ are random and incoherent then the time lag is zero. If these two quantities are correlated then nonzero lag may appear. The simplest assumption would be a linear correlation between slope and normalization (i.e. between slope and flux):

$`\beta (t)=\gamma \times (A(t)A_0)`$ (11)
where $`A_0`$ is the normalization of the time averaged spectrum and $`\gamma `$ is the coefficient of linear proportionality between variations of normalization and slope.

The expression for phase lags is then further reduces to:

$`sin\varphi ={\displaystyle \frac{Im\left[\widehat{A}^{}\widehat{(A^2)}\right]}{\widehat{A}^{}\widehat{A}}}\times \gamma ln{\displaystyle \frac{E_2}{E_1}}`$ (12)
In this approximation the sign and the amplitude of the lags would simply reflect properties of the flux variations on different time scales. Unfortunately the function (12), calculated for the observed light curves, changes the sign several times within the frequency range where hard lags are observed. This means that simple linear correlation between the flux and the slope of the spectrum can be excluded or that sign of the correlations between flux and slope also has to change sign in order to keep the sign of (12) the same at all frequencies. We note here that change of the slope/flux correlations sign at different time scales was indeed reported for Cygnus X-1 (Li, Feng & Chen 1999). Anyway it is clear that simple linear correlation between the slope and the flux fails to produce observed hard lags and more complicated models are to be invoked.

## 4 Time lags due to accretion

In this section we construct simple model which attributes the lags to the propagation of perturbations in the accretion flow from the outer to inner region. As was discussed by Lyubarskii (1997) and Churazov, Gilfanov and Revnivtsev (2001) strong variability of the X–ray flux over very broad dynamic range of time scales suggests that different time scales are introduced to the accretion flow at different distances from the compact source. One might think e.g. in terms of the effect on the mass accretion rate of the magneto-hydrodynamic turbulence, which serves as a source of the viscosity in the accretion flow through the fluctuating magnetic stresses (e.g. Balbus and Hawley 1991, Hawley, Gammie and Balbus 1995, Brandenburg, Nordlund, Stein and Torkelsson 1995). The perturbations introduced at large distance from the compact object are then propagated down to the region of main energy release and cause observed variations of the X–ray flux. If the propagation of the perturbations in radial direction is due to viscous diffusion then the very fact that we observe strong variations of the X–ray flux means that propagation/diffusion time from a given radius $`r_0`$ is comparable to or shorter than the time scale of perturbations, introduced to the flow at this radius $`r_0`$. Otherwise these perturbations would be completely washed out before reaching the innermost region. This in turn means that lags could naturally appear in such situation if locally emitted spectrum is a function of radius $`r`$. Of course energy budget of the accretion flow far from the innermost region is small, but the observed lags are also small, of the order of 0.1 radian. We then tested feasibility of this model using following assumptions:

* Statistically independent (incoherent) perturbations are introduced to the accretion flow at different distances from the compact source. The characteristic time scale of these perturbations is a function of radius.
* We further assume that the shape of the locally emitted spectra depends on radius and it is getting progressively softer with increase of the distance from the compact object.

The light curve at a given energy E is:

$`L(E,t)=`$
$`{\displaystyle D(r_0,t_0)G(r,r_0,tt_0)ϵ(r)S(E,r)𝑑r𝑑r_0𝑑t_0}`$ (13)
where $`D(r_0,t_0)`$ is an initial perturbation introduced to the flow at radius $`r_0`$ and time $`t_0`$, $`G(r,r_0,tt_0)`$ is the Green function, which describes propagation of the perturbations from radius $`r_0`$ to $`r`$, $`S(E,r)`$ is the shape of the locally emitted spectrum at the radius $`r`$ and $`ϵ(r)`$ is the total luminosity of the accretion flow emitted at a given radius. Assumption of statistical independence of perturbations introduced at different radii implies that quantities like power density or cross spectra can be calculated independently for each initial radius and then averaged over range of initial radii. The Fourier transform $`\widehat{L}(E,f)`$ of the X–ray light curve is then:

$`\widehat{L}(E,f)={\displaystyle \widehat{D}(r_0,f)\widehat{G}(r,r_0,f)ϵ(r)S(E,r)𝑑r𝑑r_0}`$ (14)
where $`\widehat{D}(r_0,f)`$ is the Fourier transform of the perturbations introduced to the flow at the radius $`r_0`$, $`\widehat{G}(r,r_0,f)`$ is the Fourier transform of the Green function, $`f`$ is the frequency. From the above equation it is clear that in this approximation time lags are solely due to the properties of the Green function and the dependence of the spectral shape on radius. For our toy model we choose the simplest possible version with $`ϵ(r)1/r^2`$ per unit $`dr`$. In order to calculate phase and time lags we created two artificial light curves in the “soft” and “hard” bands assuming simple laws of emissivity in these bands:

$`S(r)=ϵ(r)S(E_{soft},r)={\displaystyle \frac{1}{r^2}}`$
$`H(r)=ϵ(r)S(E_{hard},r)={\displaystyle \frac{1}{r^2}}\times h(r)`$ (15)
where $`h(r)`$ is a decreasing function of the radius. For simplicity we assume that propagation of the perturbations can be characterized by the Green function of Lynden-Bell and Pringle (1974) describing diffusion in the geometrically thin disk approximation (see also Lyubarskii 1997). I.e.

$`G(r,r_0,t){\displaystyle \frac{r^{1/4}}{t}}exp\left[{\displaystyle \frac{r_0^{1/2l}+r^{1/2l}}{4t}}\right]I_l\left[{\displaystyle \frac{r^{1/4l}r_0^{1/4l}}{2t}}\right]`$ (16)
where $`I_\nu (x)`$ is the Bessel function of imaginary argument. For the parameter $`l`$ we use the value of $`1/3`$ appropriate for a standard $`\alpha `$ disk with constant ratio $`\frac{H}{r}`$, where $`H`$ is a half thickness of the disk. In this equation and below we measure radius $`r`$ in units of $`r_{min}`$, which is supposed to be the characteristic radius of the inner region of the accretion flow where most of the energy is released. Time is measured in units of the characteristic diffusion time $`k^2\times r_{min}^{1/2l}`$, where $`k`$ is the effective diffusion coefficient.

Two simple versions of the function $`h(r)`$ (hereafter “hardness” function) are shown in Fig.5. The dotted line shows the case when slope of the locally emitted spectrum is the same at all distances from the compact object, but it abruptly hardens within innermost region. The solid line shows the case when locally emitted spectrum gradually hardens as the distance from the compact object decreases. In Fig.5 light curves in the soft and hard energy bands are shown, assuming $`h(r)`$ function with the abrupt jump. Here initial radius $`r_0`$ was set to 30 and initial perturbation was assumed to be $`\delta `$ function in time and radius. Initially locally emitted spectrum is soft and it hardens when perturbation spreads all the way down to the minimal radius. After that moment hardness of spectrum remains practically unchanged.

We then calculated phase lags fixing the initial radius $`r_0=30`$ and varying the characteristic frequency of the initial perturbation. The dependence of the phase lag on frequency is shown in the Fig.7. Solid line in this figure shows the phase lag (in radians) as a function of frequency expressed in units of inverse characteristic diffusion time for a given initial radius $`r_0^{3/2}`$. These lags were calculated for $`h(r)`$ function with a jump (see Fig.5), but for the smoother $`h(r)`$ the resulting lags are practically the same. Dotted line in the same figure shows the value of the transfer function power $`|\widehat{G}(r,r_0,f)ϵ(r)𝑑r|^2`$. It characterizes to what extent an amplitude of a perturbation at a given frequency $`f`$ is suppressed in the observed light curve due to diffusive spreading. Values of the order of $`10^2`$ imply that “local” perturbations at $`r_0`$ with an amplitude of the order of unity would result in $``$10% variations in the innermost region where bulk of energy is released.

In order to calculate phase lags associated with the perturbations generated at different radii we made further simplifying assumption that at every initial radius $`r_0`$ instabilities generate single frequency which scales as a characteristic diffusion time at this radius222Note that for the $`\alpha `$ disk with constant ratio of the disk height to the radius the characteristic diffusion time and Keplerian time scales similarly. I.e. $`f(r_0)=Kr_0^{3/2}`$, where $`K`$ is a parameter. From Fig.7 one would expect the lags to increase with the increase of $`K`$ unless $`K`$ is too large. The resulting phase lags for two versions of the “hardness” function are shown in the Fig.8. The resulting phase lags weakly grows with frequency and reach $``$ 0.1 radian before falling sharply. The amplitude of the phase lags and the dependence of the lags on frequency are thus regulated by two major parameters: i) the strength of the spectral changes as a function of radius and ii) the relation between the time scale of perturbations, locally generated at radius $`r_0`$, and diffusion time scales from $`r_0`$ down to much smaller radii. We stress that under assumption that at each radius $`r_0`$ perturbations with single characteristic frequency are generated the phase lags do not depend on the shape of the power density spectra and are solely determined by the Green function and the “hardness” function.

In the model, considered above, one would expect that power density spectra in different energy bands should have somewhat different shape due to the influence of the “hardness” function. In the Fig.9 we show the ratio of the power density spectra in the “hard” and “soft” bands as a function of frequency. This ratio is again independent on the shape of the power density spectrum and is governed by the the Green function and the “hardness” function. One can see that, as expected, the ratio of the power density spectra increases with frequency.

For completeness we also calculated time lags using the same choice of the “hardness” function, but instead of diffusion Green function we substitute Green function appropriate for a free fall from given radius $`r_0`$. This is equivalent to the assumption that distinct blobs of emitting matter are falling towards the compact source with the velocity $`vr^{1/2}`$. The resulting lags are shown in Fig.10. Since in this case no viscous spreading is present, the absolute value of the lags is higher for the same choice of the “hardness” function and parameter $`K`$.

Finally we stress that if the locally emitted spectrum at any radii is a power law then observed lags will show logarithmic dependence on energy even if the total observed spectrum (integrated over all radii) is not well represented by a power law as assumed in the previous section. This follows from our assumption that perturbations generated at different radii are statistically independent.

### 4.1 “Antilags” due to reflected component

As we see in Section 2 observed energy dependence of the time lags contradicts the assumption that hard lags are predominantly produced by the time delay of the reflected component with respect to the primary continuum. In particular (see Fig.4) instead of a prominent hump at the energy of the iron fluorescent line the observed lags seem to be slightly suppressed here. This behavior might hint that reflected component somehow works in opposite direction and reduces the lags. We now discuss how such suppression could be introduced within the frame of the model considered above.

An assumption that locally emitted spectrum has a power law shape with an index slowly varying with the distance (so that hardest spectra are near the center) from the compact object will naturally lead to a hard lags with a logarithmic dependence of lags on energy. Let us now assume that locally emitted spectra also contain a reflected component. Assuming that overall geometry of the accretion flow resembles a truncated optically thick disk followed by an optically thin flow one would expect that the reflection fraction is higher for spectra emitted at larger distances from the compact object. I.e. softer spectra (emitted at larger distance from the compact object) have larger reflection fraction than harder spectra. In our model this means that reflected component is leading the harder spectrum. In this situation one would expect slight modification of the logarithmic dependence of lags on energy due to reflected component. The expected shape of the lags dependence on energy in this case is shown in Fig.4 with the solid line. For simplicity we assumed that lags are due to two power law spectra with the power law indices of 2.4 and 1.8 and the reflection fraction of $``$1 and $``$0.5 respectively. I.e. the light curve is:

$`L(E,t)=S_1(E)\delta (t)+S_2(E)\delta (t\mathrm{\Delta }t)`$ (17)
Reflected components in both spectra were smeared with a $``$0.8 keV Gaussian. In the limit of small frequencies one can easily write an expression for the time lag as a function of energy for the light curve described by eq. (17). Detailed shape of the energy dependence is affected by strength of the reflection features in both spectra, amount of smearing applied, relative normalization of both spectra. One illustrative example, calculated for one particular choice of these parameters, is shown in Fig.4 with a solid curve. This curve is not supposed to closely reproduce the data, but just to demonstrate that the trend is in right direction.

## 5 Discussion

The properties of the hard lags, observed in Cyg X-1 can be summaries as following:

* In the hard state of Cygnus X-1 lags are observed very significantly in the frequency range of few $`10^2`$–few 10 Hz. Phase lag is a slightly growing function of frequency with a maximum phase lag of $`0.1`$ rad between 3 and 10 keV (see recent data of Nowak et al. 1999a, Pottschmidt et al. 2000, also Fig.11).
* An energy dependence of the lags is approximately logarithmic (Miyamoto et al. 1988, Nowak et al. 1999a, see also Fig.4) with the slight, but statistically significant, deviations from this dependence.
* The cross correlation function peaks at zero lag (Maccarone et al. 2000).
* The coherence function is close to unity (Nowak et al. 1999a), at least in the middle of the frequency range, where lags are observed with high significance

Several different branches of models have been suggested as an explanation for the observed hard lags (see Poutanen 2001a,b for recent reviews). Our model belongs to the class of “propagation” models as defined by Nowak et al. 1999b. Other models involving “propagation” are those of Kato 1989, Nowak et al. 1999b, Böttcher and Liang 1999, Misra 2000, (see also Manmoto et al., 1996). In fact our version of the model is close to the original suggestion of Miyamoto et al. 1988, Miyamoto and Kitamoto 1989 that a clump of matter … drifts from the outer (soft-X-ray-emitting) region to the inner (hard-X-ray emitting) region.

The model is based on two major assumptions: a) different time scales are introduced to the accretion flow at different distances from the compact object and propagation time scales are comparable to the time scales of perturbations and b) locally emitted spectra are power laws with the photon indices smaller in the innermost region.

The logarithmic energy dependence of the lags follows directly from the assumption (b). Zero shift of the cross correlation function and coherence function close to unity naturally appear because most of the emission is coming from the innermost region and therefore most of the flux at all energies is released at the same time (see Fig.6). The dependence of the lags on frequency and maximal values of the lags are tunable parameters of the model (through the dependence of the photon index on the radius – the “hardness” function $`h(r)`$ and through the relation of the perturbation and propagation time scales). Examples shown above (see e.g. Fig.8) show that with reasonable assumptions one can approximately reproduce both the dependence of lags on frequency and the maximal value of the phase lag.

In addition this model provides natural explanation of the ratio of the power density spectra in different energy bands as a function of frequency (see Fig.9 and 11, upper panel) and offers a qualitative explanation for the deviations of the energy dependence of the lags from a simple logarithmic law (Fig.4). Broadly the model is consistent with the geometrical model in which the inside truncated optically thick disk is followed by the optically thin flow. At the distances larger than the disk truncation radius an optically thin flow seems to be also present and perhaps has a form of a corona sandwiching the disk (see arguments in Churazov et al. 2001). The latter assumption allows one to extend this model to the soft state of Cyg X-1 without any modifications.

Our assumption that each frequency in the observed flux variability is associated with some particular radius is of course a gross oversimplification. However due to suppression of the high frequency variations in the course of diffusion such situation may in practice materialize in the accretion flow if the times scales of locally induced perturbations at given radius are high, e.g. higher than the diffusion frequency (Churazov et al. 2001). From Fig.11 it is clear that there are several “humps” in the power density spectrum of Cyg X-1. Assuming that this is an indication that perturbations arising at some particular radius of the accretion flow are dominating at a range of frequencies one would then expect the phase lags to rise faster with frequency (compare Fig.8 and 7). It seems that such behavior is indeed observed as a distinct regions of fast phase lag rise (or time lag shoulders – see Nowak et al. 1999a,b) in the Cyg X-1 data (Fig.11).

Particular emission mechanism is not very important as long as only the origin of lags is considered. E.g. emission could be due to magnetic flares as in the picture of Galeev, Rosner, Viana 1979 or Poutanen and Fabian 1999. In our model flares spectra need not to evolve (from soft to hard) themselves. They could be short and could occur at different distances from the compact object. The observed variability could then be not due to individual flares but rather due to change of the number of active flares at a given moment of time. The model requires however that the spectra of a more distant flares to be softer than those of the flares in the innermost region of the accretion flow.

## 6 Conclusions

We show that some properties of the black hole candidates spectral variability in the X–ray band (in particular time lags) can be explained by a simple phenomenological model. This model assumes that a) different time scales are introduced to the accretion flow at different distances from the compact object and propagation time scales are comparable to the time scales of perturbations and b) locally emitted spectra are power laws with the photon indices smaller in the innermost region.

We also show that energy dependence of the time lags excludes the possibility that lags (in the 0.1 – 10 Hz range) are predominantly due to the delay of the reflected component with respect to the illuminating continuum. This conclusion is derived under assumption of a “standard” neutral reflected component which linearly responds to the variations of the illuminating flux.

## 7 Acknowledgements

This research has made use of data obtained through the High Energy Astrophysics Science Archive Research Center Online Service, provided by the NASA/Goddard Space Flight Center. O.Kotov acknowedges partial support by RFBR grant 00-15-96649.


format: text

text:
# Hubble’s Nobel Prize

## 1 Introduction

The astronomer Edwin Powell Hubble never won the Nobel Prize because he died unexpectedly, at the age of almost 64, on September 28, 1953, due to a cerebral thrombosis. His long time family physician assured (!) his wife, Grace Hubble, that the death had been “instantaneous and without pain” (Christianson 1995 — hereafter CHR).

Is there a Nobel Prize for Astronomy? No, there is not! Hubble would have been the first one to brake the old tradition and to change the statutes of the award, as I show below. In the light of this possibility, one might also wonder what Hubble would win the Prize for, given his many scientific achievements.

But before that, in the next section, I describe another first rank award series, also under the auspices of the Royal Swedish Academy of Sciences, namely, the Crafoord Prize. This Prize includes astronomy in the list of awards. In the final section, I comment on the story behind the Nobel award to Hubble and speculate on the choice that would be made by the Nobel Committee, from among his many fundamental investigations in astronomy, as the Prize statement.

## 2 The Crafoord Prize

The Anna-Greta Holger Crafoord Fund was established in 1980 to promote basic scientific research through yearly donations to the Royal Swedish Academy of Sciences. Holger Crafoord (1908-1982) was very active in Swedish industry. From 1964 and on, he developed and manufactured the artificial kidney, a sort of biological dialyzer that would become of vital importance in the world. His company also developed a series of medical instruments that contributed to earn him an enormous fortune.

In 1976 he became an honorary doctor of medicine at the University of Lund, followed by his wife, Anna-Greta Crafoord (1914-1994), in 1987.

Specifically, the purpose of the Fund is to promote and award research in the fields not covered by the Nobel Prizes in natural sciences, namely, mathematics, geosciences, biosciences (with special emphasis in ecology and rheumatoid arthritis), and astronomy. The awards follow a closed cycle based on the annual sequence:

1. mathematics,
2. geosciences,
3. biosciences,
4. astronomy,
5. geosciences,
6. biosciences,
7. mathematics.

The first Crafoord Prize was awarded in 1982 to V.I. Arnold, from Moscow State University, for his contribution to the theory of non-linear differential equations. The Prize amounts to $500,000 US, a gold medal, and a diploma.

The Crafford Prize is thus every six years assigned to astronomy. The first recipient was Lyman Spitzer, Jr. (1985), then Allan R. Sandage (1991) and in 1997 there were two winners, Fred Hoyle and Edwin E. Salpeter.

More information on the Crafoord and Nobel Prizes are found in the electronic pages of the Royal Swedish Academy of Sciences at the address www.kva.se/eng/pg/prizes/index.asp.

## 3 Hubble’s Nobel

As early as the 1930s (CHR), Fred Hoyle, a frequent guest of the Hubbles at Pasadena, informed Hubble that there was a move, known of in England, by the Nobel Prize Committee in the direction of a legal amendment to the award statutes to make it possible for Hubble to be honoured with this major distinction in natural sciences. The very same rumor was also heard from Nobel laureate Robert Millikan, Caltech’s celebrated physicist.

Of course, Hubble had already been awarded many distinctive prizes, the highest being the Barnard medal, a charge of the National Academy of Sciences, which was granted to him in 1935 at Columbia University. The medal was established in 1895 and is awarded once every five years. All of Hubble’s predecessors were Nobel laureates, among them Roentgen, Rutherford, Einstein, Bohr and Heisenberg (CHR). In this case Hubble was also distinguished by being both the first American and the first astronomer to win the medal. The citation in his award was for his “important studies of nebulae, particularly of the extragalactic nebulae which provide the greatest contribution that has been made in recent years to our observational knowledge of the large-scale behavior of the Universe”.

By 1949 nothing had yet happened, but with the support of the 5-metre Palomar telescope, already in operation, Hubble’s work had gained much publicity, which could have triggered a decision by the Nobel Committee (CHR). The final word was soon given when Enrico Fermi and Subrahmanyan Chandrasekhar joined their colleagues in the Committee unanimously voting Hubble the 1953 Prize in physics (CHR). But it was too late, Hubble’s death came first. Incidentally, one should recall that the Nobel prize is not awarded posthumously.

These are the facts. Let us now speculate on an alternate universe, one in which Hubble did survive that unfortunate September afternoon.

The immediate question is what Hubble would be cited for in the Nobel award. The conservative approach usually adopted by the Nobel Committee is well known. The classical example is Einstein’s. Awarded the 1921 prize, he was specially cited for his discovery of the law of the photoelectric effect. It is needless to say that both the Special and General Relativity Theories, which had already been put forward, were not explicitly mentioned.

Allan Sandage, a Crafford laureate and the greatest of Hubble’s followers, in a paper celebrating the centennial of the birth of Hubble (Sandage 1989), enumerates Hubble’s four central accomplishments undertaken from 1922 to 1936. Sandage adds that any one of them guarantees Hubble a place in the history of modern science. They are:

(a) the morphological sequence of galaxy types,

(b) the discovery of Cepheids in NGC6822, with parallel work in M31 and M33, settling decisively the question of the nature of galaxies,

(c) the determination of the homogeneity of the distribution of galaxies, averaged over many solid angles, and

(d) the linear velocity-distance relation.

All of the four but one may be blurred with controversies, not only in modern times but also, and certainly, in Hubble’s time. Nowadays it is recognized that the so-called Hubble “tuning-fork” scheme in (a) applies mainly to close and bright galaxies, that the homogeneity in (c) is broken by the presence of enormous structures like “walls”, “streams” and “voids” of galaxies, and, finally, even Hubble himself never clearly advocated the idea of a velocity-distance relation in (d); rather he usually put it as a relation between spectral shift and distance, as is evident throughout the classical book The Realm of the Nebulae (Hubble 1936). In short, these three otherwise fundamental advances in modern astronomy seem to collide with the spirit of Nobel citations, i.e., they are not solid statements about nature.

This shortcoming is not the case, by any means, with item (b) above. Hubble’s work here is the end point of the great debate about the nature of the nebulae. With the unambiguous determination of the distances to the “spiral nebulae” to be much larger than the dimensions of our stellar system, Hubble proved definitely that they are extragalactic, and as a result laid down the foundations of a new branch of research, namely, extragalactic astronomy. The other three points mentioned by Sandage are immediate consequences of this major realization, which were soon recognized as such by the genius of Hubble.

Thus, the Nobel prize to Edwin Powell Hubble goes for his contribution to the definitive understanding of the nature of the nebulae and for the creation of a new era of scientific investigation.

## 4 References

Christianson G E 1995 Edwin Hubble Mariner of the Nebulae, The University of Chicago Press, Chicago

Hubble E P 1936 The Realm of the Nebulae, Yale University Press, New Haven

Sandage A R 1989 The Journal of the Royal Astronomical Society of Canada 83, 351

Acknowledgments — I wish to thank Dr. P.S.S. Guimarães for reading the manuscript and for useful suggestions.


text:
# Strong gravitational lensing by multiple galaxies

## 1 Introduction

Gravitational lensing by galaxies is usually modelled using a single lens. Observational results from the CLASS survey (Myers et al. 1995) suggest a probability $`p_\mathrm{s}`$ for gravitational lensing of a source at $`z>1`$ of about $`10^3`$, comparing well with the theoretical prediction by Pei (1995). The probability that two galaxies both lie close enough to the line of sight to strongly lens a background galaxy might be expected to be roughly $`p_\mathrm{s}^210^6`$ and so one might conclude that it is very unlikely for a source to be lensed strongly by two galaxies at different redshifts. However, even though the number of known lens systems to date is still relatively small – about 50 strong lens systems are known (Falco et al. 1999) – there is strong evidence that single lens models cannot explain the image geometries and magnification ratios in all cases. First, the observed image positions and magnifications in some systems, like the ‘Cloverleaf’, require either an unreasonable large mass for the lens or a substantial component of external shear (Kneib et al. 1998; Kneib, Cohen & Hjorth 2000; Soucail et al. 2000) in order to explain the details of the observed image geometry. In some cases, the lensing galaxy has been found to be part of a compact group of galaxies, and including the potential of the group greatly improves the lens model (Keeton & Kochanek 1997; Keeton, Kochanek & Seljak 1997; Kundić et al. 1997; Möller & Natarajan 2000). Second, it has been shown recently for two lens systems detected in the CLASS survey that the strong lensing potential is due to two galaxies at different redshifts, suggesting that double lens systems are not that uncommon (Koopmans et al. 1999; Koopmans & Fassnacht 1999).
In the next years, it is expected that using a new generation of instruments, lens surveys will increase the number of known lens systems by at least a factor of about 10 (Blain 1996; 2000). The large sample of strong lens systems should then contain a significant number of cases in which a background source is lensed by two foreground galaxies.
Gravitational lensing by multiple galaxies has been investigated before by Kochanek & Apostolakis (1988, KA98 hereafter) and Seitz & Schneider (1992). In KA98 the authors studied lensing by two lenses at the same and at different redshifts using a very similar ray-tracing routine to that used in this paper. However, due to the severe limits imposed by the computing power available at the time, their work was necessarily restricted to the study of only a very small fraction of the parameter space and did not include lens evolution. Multiple lensing has also been investigated in the context of microlensing, where the lenses are point masses (Lewis et al. 1993; Schneider, Ehlers & Falco 1992; Wambsganss, Witt & Schneider 1992). In a cosmological context, multiple lensing is naturally incorporated in simulations which combine ray-tracing methods with N-body simulations (Jain, Seljak & White 2000; Wambsganss, Cen & Ostriker 1998). However, due to the resolution limit of the N-body simulations, such work is concerned with the lensing effect of large-scale structure, on cluster and super-cluster scales, and not with the strong lensing effect due to individual galaxies.
In this paper, we investigate the expected fraction of double galaxy lens systems and their statistical and characteristic properties. In section 2 we outline the ray-tracing method and the model of the evolving galaxy population used. In section 3 we show how double lens systems differ qualitatively from single lens systems in terms of both lensing cross sections and image configurations. In section 4 we estimate the probability of double lensing as a function of source redshift. In section 5 we look at the statistical image geometries, such as the ratio of quadruple:double images and discuss the expected distribution of the masses and redshifts of the lenses. We discuss some more elaborate multiple lens models qualitatively in section 6 and look at some observational issues in Section 7.

## 2 Method

### 2.1 Ray tracing

The method used in this paper is based upon the ray-tracing routines developed and described in Möller (1997) and Möller & Blain (1998, MB98 hereafter). The statistical work in this paper requires a large number of multi-plane ray-tracing calculations which made it necessary to modify the routines in two ways:

1. the deflection angle is calculated in two planes at different redshifts using the full multi-lens equation (Schneider, Ehlers & Falco 1992). For two lenses, the lens equation becomes:

$$\stackrel{}{\beta }=\stackrel{}{\theta }\stackrel{}{\alpha _\mathrm{A}}\frac{D_{\mathrm{AS}}}{D_{\mathrm{OS}}}\stackrel{}{\alpha _\mathrm{B}}\frac{D_{\mathrm{BS}}}{D_{\mathrm{OS}}},$$
(1)
where $`\stackrel{}{\beta }`$ is the source position, $`\stackrel{}{\alpha _\mathrm{A}}`$ is the deflection due to the near lens A, $`\stackrel{}{\alpha _\mathrm{B}}`$ is the deflection due to the far lens B and $`D_{\mathrm{AS}}`$, $`D_{\mathrm{BS}}`$ and $`D_{\mathrm{OS}}`$ are the angular diameter distances between source plane S and lens A, lens B and the observer respectively.
2. the code uses an adaptive grid to find all the images and determine the magnifications and shears on both the image and source planes. In order to generate the adaptive grid on the image plane, a coarse grid of triangles is first lensed and a coarse magnification map is obtained, as described in MB98. Each element on the coarse grid is then divided again into $`N\times N`$ sub-elements, where $`N`$ is chosen to be proportional to the magnification: $`N=2+(\mu 1)/33`$, for $`\mu <100`$ and $`N=5`$ otherwise. The resulting adaptive grid is then divided into triangles and lensed. The position of the resulting images are then found and the magnifications calculated in the usual way. If necessary the procedure can be repeated, but we found that a single iteration was sufficient for this application. This improves the efficiency and achievable resolution of the simulations by more than an order of magnitude. A different choice of $`N(\mu )`$ is possible, but we did not find a scaling that provides a significantly better performance than the one we used.

In this paper we mainly consider spherical pseudo-isothermal mass distributions (PIMD) as deflectors, but the routine can deal with any parametric spherical or elliptical lens profile.

### 2.2 The lens profiles

Since we are interested mainly in a qualitative investigation of double lensing involving two galaxy lenses we consider here only simple spherical mass distributions. Galaxy lenses are often modelled using a singular isothermal sphere (SIS). Such a model is, however, unphysical due to the infinite central surface mass density and total mass. Real galaxies are more realistically modelled using a pseudo-isothermal profile with a projected surface mass density of the form (Kneib et al. 1996, Natarajan & Kneib 1997):

$$\mathrm{\Sigma }(R)=\frac{\mathrm{\Sigma }_0r_0}{1r_0/r_\mathrm{c}}\left(\frac{1}{\sqrt{r_0^2+R^2}}\frac{1}{\sqrt{r_\mathrm{c}^2+R^2}}\right),$$
(2)
where we choose a small core radius $`r_0=0.1\mathrm{kpc}`$ and a cut off radius of $`r_\mathrm{c}=100\mathrm{kpc}`$. The total mass enclosed is then,

$$M_{\mathrm{tot}}=\underset{r\mathrm{}}{lim}M(r)=2\pi \mathrm{\Sigma }_0r_0r_\mathrm{c},$$
(3)
and the deflection angle at impact parameter $`R`$ is

$$\left|\stackrel{}{\alpha }_{\mathrm{A},\mathrm{B}}\right|=\frac{4GM_{\mathrm{tot}}}{(r_\mathrm{c}r_0)Rc^2}\left[\sqrt{r_0^2+r^2}\sqrt{r_\mathrm{c}^2+r^2}+(r_\mathrm{c}r_0)\right].$$
(4)
For realistic, small core radii, $`r_00.1\mathrm{kpc}`$, the lensing cross section does not differ significantly from that of a truncated singular isothermal profile. However, since a singularity in the mass profile reduces the number of images by one, the number of images and image geometries are different in a singular and non-singular model.

### 2.3 The lens population

In order to investigate the statistical properties of double lensing, we created a list of 5000 lens pairs along the line of sight to a given source at redshift $`z_\mathrm{s}`$. The sample is determined randomly, using a simple Monte–Carlo sampling method (Press et al. 1988) from a Press–Schechter distribution function (Press & Schechter 1974),

$$dN(M,z)=\frac{\rho }{\sqrt{\pi }}\frac{\gamma }{M^2}\left(\frac{M}{M^{}}\right)^{\gamma /2}\mathrm{exp}\left[\left(\frac{M}{M^{}}\right)^\gamma \right]dM,$$
(5)
where $`\rho `$ is the mean smoothed density of the universe, $`\gamma =1+n/3`$ relates to the initial power spectrum index $`n`$, where we choose $`n=1`$, corresponding to a scale-invariant spectrum, and

$$M^{}(z)=M_0^{}(1+z)^{2/\gamma }$$
(6)
is a characteristic bound mass at redshift $`z`$. We convert the mass from the Press–Schechter function into a velocity dispersion for the PIMD model of the lens assuming a cut-off radius for the lens mass distribution of $`r_\mathrm{c}=100\mathrm{kpc}`$. Thus, the conversion from mass $`M`$ to velocity dispersion $`\sigma _\mathrm{v}\sqrt{GM/2r_\mathrm{c}}`$. With this conversion an $`M^{}`$ galaxy has an approximate maximum velocity dispersion of $`245\mathrm{km}\mathrm{s}^1`$ if $`M^{}=3.6\times \mathrm{\hspace{0.17em}10}^{12}\mathrm{M}_{}`$ (Blain, Möller & Maller 1999). As shown in Fig.1 this cut-off radius is chosen so that the largest fraction of lenses will have a separation of $`1^{\prime \prime }`$ as is found in the CLASS survey (Helbig et al. 1999). For SIS lenses, and hence also, approximately, for PIMD lenses with small cores, the Einstein radius $`R_\mathrm{E}`$ increases with the velocity dispersion $`\sigma _\mathrm{v}`$ as $`R_\mathrm{E}\sigma _\mathrm{v}^2`$. Note that there is a degeneracy between the choices of $`r_\mathrm{c}`$ and $`M^{}`$; neither value is constrained significantly by observations. The particular values are chosen so that the value of $`M^{}`$ agrees with Blain et al. (1999).

1. the mass of both lenses has to lie in a range $`M_{\mathrm{min}}<M<M_{\mathrm{max}}`$, where the values of $`M_{\mathrm{min}}`$ and $`M_{\mathrm{max}}`$ are chosen so that less than 2 per cent of the single lensing cross section is due to objects lying outside this range. As shown in Fig. 2, the fraction of lenses with masses less than $`10^{13}\mathrm{M}_{}`$ for a value of $`M^{}=3.6\times 10^{12}\mathrm{M}_{}`$ and $`z_\mathrm{s}=2`$ is greater than 99 per cent.
2. All lenses with an Einstein radius $`R_E<0.1^{\prime \prime }`$ are excluded from the sample. As shown in Fig. 1, lenses with smaller Einstein radii contribute less than 2 per cent to the total lensing cross section.

Overall, the error on the calculated lensing cross section for double lenses that is introduced by these selections is less than 5 per cent.

### 2.4 Placement of lenses

We choose random pairs of galaxies from the Press–Schechter distribution function which are placed at uniformly random positions on the lens plane inside a $`4^{\prime \prime }\times 4^{\prime \prime }`$ field. The magnifications and image configurations for each case are then obtained using the ray-tracing method. The parameters for the simulation are summarized in Table 1.

## 3 Properties of Individual double lens systems

### 3.1 Caustics and magnification of point sources

Two deflectors that are in close proximity to each other will modify their respective caustic structures in a way that depends on their redshifts, mass profiles and separation. As a first step toward an understanding of double lensing, we used our routine to study the qualitative dependence of the magnification maps on the parameters of such systems.
Panel (b) of Figs. 3, 4, 5 & 6 show magnification maps and critical lines for two PIMD lenses at different redshifts for four different lens separations. The parameters for the figures are summarised in Table 2.

These maps demonstrate four generic regimes that occur for double lenses. It is immediately apparent that the two lenses have a strong effect on each other – the magnification maps differ significantly from those of isolated spherical lenses. In particular, there are large high-magnification regions along the caustics. In Fig. 3 the two lenses are very close to each other and produce a joint caustic that is extended in a direction perpendicular to their separation. In Fig. 4 the lenses are further apart, but their individual Einstein radii still overlap. The caustic encloses a large area of high magnification in the source plane. In Fig. 5 the lenses are separated by more than the sum of their individual Einstein radii. In this geometry the caustic is extended primarily along the direction of the lens separation. In Fig. 6 the separation is sufficiently large that the Einstein radii do not overlap. However, there is still a ‘trail’ of high magnification between the two lenses in the source plane and there is sufficient shear to distort their individual caustics into astroids.
In Fig. 7 we show the cross section for magnification $`\mu `$ of point sources above a threshold value $`A`$ for the four configurations of Figs. 3-6.

Between magnifications of $`\mu 20`$ and $`\mu 50`$ there is an increase of up to a factor of ten as compared to the sum of the cross sections of individual lenses. The range of magnification for which there is a significant increase in the cross section depends on the degree of overlap between the Einstein radii of the lenses.
The results show qualitatively that double lens systems are much more likely to produce high magnifications of $`\mu >20`$ for point sources as compared with isolated lenses. This discrepancy between the magnification distribution for double lens systems and the sum of the curves for the individual lenses demonstrates that magnification bias for double lens systems is expected to be large; hence the rough estimate of the double lensing probability as $`10^6`$ given in the introduction is too small.

### 3.2 Image geometries

The images of point sources and one small, extended source are shown in panels (b) of Figs. 3-6 for the different source positions marked in the panel (a) to the left. In these panels, *all* images are shown, irrespective of their magnification. Those images which lie in the central region of the lenses will be strongly de-magnified and will, in most cases, not be observable. Discounting those images the panels show that, due to the extended shape of the caustics, four–image geometries are common. Point sources that lie in the high magnification region between the two lenses produce characteristic aligned triple images. Extended sources in the same region produce straight arcs with counterimages that lie on the opposite side of one of the lens centres.
In small–separation double lens systems with different redshifts the more distant lens can lie within the Einstein radius of the nearer lens, and so the more distant lens can itself be multiply imaged, leading to a total of six observed lensed images. The position of the images of the more distant lens is marked by a cross in panels (b) of Figs. 3-6. Also, for such small–separation systems, highly magnified sources are likely to have high image multiplicities, with three or five magnified images. For systems of intermediate separations, as shown in Fig. 4, three magnified images may appear in a peculiar triangular configuration, with large and nearly equal separations between the three images. Systems with intermediate separations are also the most likely to produce triple aligned images or straight arcs, as in Fig. 5. Large separation double lens systems, as in Fig. 6, produce image configurations very similar to those of elliptical lens galaxies.

## 4 Statistics of double lensed systems

### 4.1 Definition of double lenses

In order to estimate the number of lens systems in which a background source is lensed by two foreground objects, it is first necessary to adopt a clear definition of what is meant by a “double lens”. The number of systems in which a second galaxy introduces only external shear is certainly larger than the number of systems for which the background object would be multiply imaged by both lenses individually. Here we will define two regimes. A “weakly coupled double lens” is a lens system for which the caustics for each of the two lenses do not merge, but for which there is still a significant effect on the individual caustic structures. A “strongly coupled double lens” is a lens system in which the caustic structures merge. In other words, strongly coupled systems have a single connected multiple imaging region in the source plane, whereas weakly coupled systems do not. This distinction has the advantage that it is easy to classify objects in each regime from the topology of the inner caustic lines. A line connecting the two lenses will cross the high magnification caustic line less than twice only in the strongly coupled case. For example, the lens in Fig. 6 is a weakly coupled double lens, whereas that in Fig. 4 is a strongly coupled double lens.

### 4.2 Double lensing probability

The relative probability of double lensing is given by the cross section ratio of double lensing to lensing by the individual galaxies. The cross section due to double lensing is given by

$$_\stackrel{}{\mathrm{p}}_0^{z_s}\sigma (\stackrel{}{p},z_s)n(\stackrel{}{p})\frac{dV_{\mathrm{co}}}{dz}𝑑z𝑑\stackrel{}{p},$$
(7)
where $`\sigma (\stackrel{}{p})`$ is the cross section for multiple images by a double system with parameters $`\stackrel{}{p}`$, $`n(\stackrel{}{p})`$ is the comoving number density of such systems and $`V_{\mathrm{co}}`$ is the comoving volume at redshift $`z`$. In our model the parameters of the system are the two redshifts, $`z_1`$ and $`z_2`$, the masses, $`m_1`$ and $`m_2`$, and the separation of the two lenses $`\mathrm{\Delta }\theta `$. The comoving number density of objects

$$n(\stackrel{}{p})=n(m_1,z_1)f(m_2,z_2,\mathrm{\Delta }\theta ),$$
(8)
where $`f(m_2,z_2,\mathrm{\Delta }\theta )`$ is the probability of finding another lens with mass $`m_2`$ and redshift $`z_2`$ at a distance $`\mathrm{\Delta }\theta `$. We assume no spatial correlation and so $`f(m_2,z_2,\mathrm{\Delta }\theta )=f(m_2,z_2)\times f(\mathrm{\Delta }\theta )`$. To proceed further it is necessary to compute the lensing cross section $`\sigma (\stackrel{}{p},z_s)`$. Even though this is possible analytically in the case of simple lenses (KA98), the necessary formalism is cumbersome and cannot be extended to more elaborate lens models. Since the ray-tracing code described above is both extremely fast and accurate, we use it to obtain the form of $`\sigma (\stackrel{}{p},z_s)`$ numerically. Instead of sampling the function $`\sigma (\stackrel{}{p},z_\mathrm{s})`$ at regular intervals, to obtain an approximate functional form, we solve the complete integral in equation 5 numerically in a Monte–Carlo fashion. This is done by sampling lenses randomly from a Press–Schechter distribution as described in Section 2 and obtain the value of $`\sigma (\stackrel{}{p},z_s)`$ for each lens system using ray tracing.

### 4.3 Numerical results

Using the lens population and ray tracing lensing routines described in Section 2, we obtain the expected number of double lens systems numerically for ten source redshifts in the range $`0.5<z_s<10`$. We calculate the total cross section for strong lensing by two lenses by summing the cross sections for multiple imaging calculated for each lens pair in our sample of 10,000 lenses for each source redshift, and then multiply this by the probability that two lenses with the given properties are found inside the $`4\mathrm{}\times 4\mathrm{}`$ field. To obtain a conservative double lensing probability, we assume that there is no spatial correlation.
To test our procedure we also calculated the total multiple lensing cross sections due to the individual lenses. Fig. 8. shows the theoretical single lensing probability and the results from our simulations for the different source redshifts.

As can be seen, the results agree well with the theoretical predictions. The errors are statistical and due to the finite number of lenses used – the numerical error introduced by the ray-tracing routine is several orders of magnitudes smaller.
Also shown in Fig. 8 is the probability of double lensing as a function of source redshift, for both weak and strongly coupled double lens systems. The dashed and dotted curves shown are obtained by squaring the probability for single lensing and normalising to the double lensing probability at $`z_\mathrm{s}=10`$. Not surprisingly, the double lensing probability is consistently lower than that for single lenses. However, at redshifts larger than about 2 it is above the simple estimate from the introduction. This shows that the change in caustic structure due to the double lens potential cannot be neglected. Most noteworthy is that for source redshifts of about 5, about one in twenty lenses is expected to be a strongly coupled double lens. Also, the number of double lens systems at redshifts $`z_\mathrm{s}12`$ is in good qualitative agreement with the number of double lenses discovered so far in a radio–selected sample (two in a sample of 50).
We have not assumed any correlation in our lens sample. As the galaxy two-point correlation function is positive on these scales (Peebles 1993), our result is therefore likely to be an underestimate of the number of double lens systems in the case that the two lenses are at the same redshift. For a correlation of the form

$$\xi (r)=\left(\frac{r}{r_0}\right)^\alpha ,$$
(9)
with $`\alpha 1.8`$ and $`r_0=5h^1\mathrm{Mpc}`$ at $`z=0`$ we estimate that the two–point correlation increases the relative number of systems with separations below 1” by a factor of about 4. This means, that there will be an increase in the double lensing probability for correlated lenses at the same redshift relative to that for uncorrelated lenses of roughly the same factor. As shown in Fig. 7 double lenses are much more likely to produce high magnifications than individual lenses. Since flux limited lens surveys are biased towards the discovery of highly magnified images, the ratio of double lens systems to single lens systems in these surveys is likely to be higher than our prediction by a factor of a few.

## 5 Statistical properties

### 5.1 Lens properties

We have shown in the previous section that a large sample of gravitational lens systems is likely to contain a significant fraction of double lens systems. It is instructive to look at the properties of the double lens pairs in the sample.
Fig. 9 shows the contribution to the total lensing cross section as a function of the ratio of the two individual lens masses in panel (a) and as a function of the ratio of the two lens redshifts in panel (b). It can be seen that most of the double lens pairs in this sample have lenses lying at similar redshifts and with similar masses. This is not surprising, as the cross section for lensing of sources at a given redshift peaks at an optimal value of the redshift and mass of the lens galaxy. Positive spatial correlation will increase the number of double lens systems that are likely to lie at similar redshifts further. Note, however, that for uncorrelated pairs the efficency of lensing as a function of redshift is far less peaked for the second lens than that for a single lens (see, for example, Fig. 5 in Fukugita et al. 1992). In fact, the likelihood that the ratio of redshifts of uncorrelated double lenses is 2 or larger is $`50`$ per cent.

### 5.2 Statistical image geometries

Image geometries and magnification bias are strongly dependent on the mass profile of the lens, and so a two lens system is expected to lead to unusual and complex image geometries. In Section 2, we showed some individual image geometries for extended sources. To determine the statistical properties of the images, we generated images for small sources in the double lenses in our sample. The histogram of the cross section for the number of images is shown in Fig. 10.

We show the results for two different sets of observational selection criteria. In panel (a) the statistics include all images, in panel (b) only images with minimum separation of $`0.05\mathrm{}`$ and a maximum magnification ratio of 100 are included and in panel (c) only images with minimum separation of $`0.05\mathrm{}`$ and maximum magnification ratio of 20 are included. The histograms show clearly the increased cross section towards three ($`15\%`$) and four ($`10\%`$) images. Individual spherical PIMD lenses can only produce two magnified and one demagnified image. Strongly coupled double lens systems can produce three or more magnified images. Strongly and weakly coupled double lenses both lead to quadruple systems and to a small fraction of five- and six-image systems. However, since the fraction of double lens systems will be small, and the vast majority of all lens galaxies are expected to have some effective elliptical profile, which increases the cross section for the formation of four or more images in a similar way, the effect of double lenses on the overall image statistics in a large sample of lenses will be small. Note that the image configuration shown in Fig. 5(b) is characteristic of double lens systems and hard to reproduce in most single lens models. Also, three strongly magnified images, without an additional counter-image, is a clear sign of a possible double lens system. Magnification bias in any flux limited lens sample will increase the relative number of lens systems with high image magnifications and hence high image multiplicities. This effect will increase the fraction of systems with 3 or more images by a few per cent.

## 6 Extending double lens models

### 6.1 Double lens plus external shear

Our above analysis does not assume any correlation between the lens positions, however, the spatial positions of galaxies on the sky are in reality correlated, and most galaxies occur in small groups or clusters. To show qualitatively how other galaxies in the environment affect double lensing, we model the potential perturbation due to a nearby group or cluster as an external shear that acts on the double lens system. In Fig. 11 we show the magnification maps for the double lens system shown in Fig. 5, now modified by an external shear, which is assumed to be perpendicular to the alignment of the lenses in Fig. 11(a) and parallel to the alignment of the lenses in Fig. 11(b).

We model the source of the external shear as a point mass of $`7\times 10^{12}\mathrm{M}_{}`$ at a distance of 14 arcsec. In Fig. 11(a) the external shear “bends” the caustic structures in its direction, but the effect on the caustic area is small. It is only when the external perturber lies close to and is nearly perfectly aligned with the lens pair, as in Fig. 11(c), that the size of the caustic is affected significantly. Such a situation is likely to be rare, however, and so we expect that, on average, the effect on lensing statistics due to external shear from perturbers is likely to be small. However image positions and magnifications can be affected greatly by even moderate external shear, and so external perturbers have to be taken into account in accurate lens modelling of individual double lens systems. A more detailed treatment of the effect of external perturbers like groups of galaxies on the properties of galaxy lenses can be found in Keeton, Kochanek & Seljak (1997) and Möller, Natarajan & Kneib (in preparation).

### 6.2 Spiral galaxies

In the above analysis we kept the number of parameters as small as possible and used the PIMD as a simple but reasonable lens model. Including spiral discs into the lens model significantly changes the lensing behaviour of individual systems, as shown by Maller, Flores & Primack (1997) and MB98. The statistical lensing properties of spiral lenses has been investigated by Bartelmann & Loeb (1998), Keeton & Kochanek (1998) and Blain, Möller & Maller (1999). These studies showed that there is a significant effect on the properties of lensing by a single galaxy. We show qualitatively how the presence of spiral discs will affect double lensing in Fig. 12.

Both lenses are spiral galaxies similar to the Milky Way, which consists of an PIMD halo and an inclined thin exponential disc, of central surface mass density $`\mathrm{\Sigma }_0`$ and scale length $`r_\mathrm{s}`$, containing about 10 per cent of the halo mass. The effect of the discs on the lens properties depends strongly on the relative alignment of the discs. If they are perpendicular to one another, then the shear along the caustics is greatly reduced and the area enclosed by the caustics shrinks drastically. If both discs are aligned with each other and the angle between their major axes is less than about 60 deg, then the caustic lines can enclose a very large high–magnification area as seen in Fig. 12(a). It is difficult to assess the effect of spirals on double lens systems in a statistical sense due to the large parameter space. The effect of discs on double lens statistics is expected to be small if the alignment of the discs is uncorrelated, as the asymmetry that is introduced in the potential by the disc will *on average* counteract the ellipticity introduced due to a second lens. Individual double spiral galaxy lenses will in general have different properties than double spherical or elliptical lenses. Spiral galaxies are at least twice as abundant as elliptical galaxies and therefore, since the inclination effect due to the discs enhances the high magnification cross section, double lens systems containing one or two spiral galaxies are likely. Observations of such lens systems could provide strong constraints on the mass profile of one or both lenses.

## 7 Observations of double lenses

The observed properties of double lenses are likely to depend more strongly on the wavelength of observation than those of single lenses. This is due to the fact that light from the images and the more distant lens (or its images if it is multiply lensed by the nearer lens) is likely to be superimposed on the light distribution of the near lens galaxy. In addition, dust in the interstellar medium of the lens galaxies could lead to significant extinction of one or more of the images. Double lens systems will only be easily observable if both lenses are relatively faint at the wavelength of observation. This would be the case, for example, in the radio or sub-mm waveband if the two lenses are elliptical galaxies and the background source is a small young star forming galaxy. It would not be the case for optical observations of two massive spiral galaxies lensing a distant galaxy. In this case images are likely to be too faint relative to the lenses to be observable at most wavelengths. As quasar radio surveys are not affected by extinction, a complete lens sample from a radio survey, like the CLASS survey, would therefore be especially suited to observe double lens systems. The Planck mission will discover many thousands distant sources serendipitously at sub-mm wavelengths of which a fraction of order 10 per cent could be lensed (Blain 1998a). We predict that about 5 per cent of these lens systems will be double lenses. The Atacama Large Millimetre Array (ALMA) which will observe at sub-mm wavelengths would be well suited to detect the images of distant sources (Blain 1998b;2000). The subarcsecond resolution of ALMA would resolve the individual images and allow their direct study, which would greatly improve the accuracy of lens mass models.

## 8 Conclusions

In this paper we have investigated gravitational lensing of high-redshift background sources by more than one galaxy along the line of sight. Using a Press–Schechter halo distribution and a ray-tracing code we have estimated the number of double lens systems that are to be expected in a large lens sample. We have discussed the properties of such double lens systems and investigated more complicated double lens models qualitatively. In summary, our main results are:

1. In a cosmology with $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$ and $`\mathrm{\Omega }_\mathrm{M}=0.3`$, about 2-5 per cent of all multiply imaged sources at $`z2`$ or higher are expected to be lensed by more than one lens along the line of sight.
2. The second lens induces a strong asymmetry in the effective lensing potential. This leads to a significant change in the caustic structure. The cross section for high magnification of point sources increases significantly due to this effect.
3. Double lenses lead to a significant fraction of lens systems with three ($`15\%`$) and four ($`10\%`$) images, and can lead to five- and six-image configurations.
4. The two lenses in a double lens system are likely to be of similar mass and redshift.
5. Additional external shear acting on a lens pair can modify the caustic structure and the image geometries in individual systems significantly, and needs to be included in lens modelling just as for single lens systems.
6. Double spiral galaxy lenses can have a large high-magnification cross section, if their postion angles are aligned with each other, and the inclination of both of the discs towards the line of sight is higher than about 65 deg.
7. Future lens surveys, especially in the submillimetre wavebands will contain a significant number of double lens systems.

## Acknowledgements

OM acknowledges the Lensnet TMR network and AWB the Raymond and Beverly Sackler Foundation for support. We thank Priya Natarajan, Daniel Mortlock and the referee, Jean-Paul Kneib for useful comments on the manuscript.


text:
# A Feynman graph selection tool in GRACE system

## 1 Introduction

Using an automatic Feynman graph calculation package, we can generate the information of all Feynman graphs for given processes. Sometimes it is necessary to select graphs from the set of all graphs by some conditions. However, it is not so easy to select them correctly by hand when a huge number of graphs is involved, such as higher order corrections or SUSY processes.

A program grcsel selects out a subset of Feynman graphs from the set of graphs, generated by GRACE, according to given selection conditions. The output information of selected graphs is written in the same format as that of the original set. This enables us to generate Feynman amplitudes within GRACE in the same procedure as for all graphs. So we can perform cross section calculation, gauge invariance check, event generation and so on for the selected ones. grcsel helps us:

1. to find decay-graphs and evaluate signal/background ratio,
2. to check the accuracy of approximated calculation,
3. to confirm precision of the calculation,
4. to reduce the calculation time,
5. to develop kinematics routine.

## 2 Overview of grcsel

grcsel consists of three parts: a steering-part defines basic functions of graph selections and reads input files, an interpreter-part parses and evaluates commands, and a utility-part handles subsets of particles, vertices and graphs.

Once a physics process and the order of calculation are fixed, Feynman graphs are generated by grc program with specified Feynman rules described in physics model file. The information on graphs generated are stored in a file named out.grf (we call the format of this file .grf format).

grcsel reads the physics model file and out.grf and selects graphs according to a kind of propagator, characteristics of graph topology, a type of vertex or a graph number. grcsel outputs those selected graphs in the same format as out.grf. Successively this output file can be used as the input to source code generation for Monte Carlo integration or event generation. We can also use grcsel again reading output of previous execution of grcsel. The schematic view of how grcsel works in GRACE system is shown in Fig. 1.

### 2.1 Running grcsel

Graph selection starts by the program grcsel:

grcsel

This program requires out.grf file by default. The graph selection commands are read through standard input, which may be given interactively or by a script file. With a script file where grcsel commands are prepared, we can redirect that file:

grcsel \< command.in

To use another input .grf format file, e.g. out1.grf, instead of out.grf, we can add the filename after grcsel command as:

grcsel out1.grf \< command.in

### 2.2 grcsel command

In a script file there are a series of grcsel commands such as declaration of variables and basic functions to specify the selection conditions or operators. grcsel has 14 basic functions in total. Three of them return a subset of graphs in accordance with specified selection condition and two functions output set of graphs. They are summarized in Table 1.

In grcsel, graphs, particles or vertices are treated as elements of a set of type gset, pset or vset, which are defined as a set of graphs, particles or vertices, respectively. Set variables have to be declared at first with their types. Operations on sets are available and are shown in Table 2.

## 3 Example

In the following example, graphs with $`\nu _e`$ propagator connected to the initial electron and final $`W^{}`$ are selected among graphs of $`e^+e^{}W^+W^{}\gamma `$ process. Selected graphs are output into a file named out1.grf.

```
% e+ e- --> W+ W- Photon
%
% out1.grf : with neutrino propagator
% at the vertex of initial electron and
% final W-.
%

gset gs0, gs1;

gs0 = ~[];   % all graphs

gs1 = cutprop(gs0, [‘‘nu-e’’], [0,3]);

outgset(‘‘out1.grf’’, gs1);

quit;

```

In Fig. 2, selected graphs are shown.

## 4 Remarks

grcsel has been developed in the framework of GRACE 2.1.7.4. It can handle tree and 1-loop graphs and it supports standard and MSSM physics model. grcsel is included in a distribution kit of GRACE 2.1.7.4.

Acknowledgments

We wish to thank the members of MINAMI-TATEYA collaboration for continuous discussions and many kinds of support. We are also grateful to express our sincere gratitude to Prof. Y.Shimizu for the valuable suggestions and continuous encouragements. Authors appreciate Prof. Y.Watase for the encouragements. This work was supported in part by the Grant-in Aid (No. 12680363, 10640285, 10680366 and 11440083) of Monbu-sho, Japan.


text:
# Mass Density Profiles of LSB Galaxies

## 1 Introduction

Low Surface Brightness (LSB) galaxies are dark-matter dominated galaxies where the stellar populations only make a small contribution to the observed rotation curves. It is therefore straight-forward to compare the observed rotation curves of these galaxies with those derived from numerical cosmological simulations, where the dark matter is the dominant component.

Early observation of dwarf and LSB galaxies showed that their rotation curves rose less steeply than predicted by numerical simulations based on the Cold Dark Matter (CDM) hypothesis (Moore, 1994; Flores & Primack, 1994; de Blok & McGaugh, 1997; McGaugh & de Blok, 1998). In the CDM model, halos are characterised by a steep central cuspy power-law mass density distribution $`\rho (r)r^\alpha `$. Initial simulations indicated that $`\alpha =1`$ (Navarro, Frenk & White, 1996). More recent results indicate a more extreme value $`\alpha =1.5`$ (e.g. Moore et al., 1998, 1999; Bullock et al., 1999). Rotation curves of dwarf and LSB galaxies, however, show a more solid-body-like rise consistent with a mass distribution dominated by a central constant-density core ($`\alpha 0`$), and hence are inconsistent with the CDM predictions. Similar conclusions have been reached by Salucci (2001) and Salucci & Borriello (2000) for high surface brightness disk galaxies.

The conclusions regarding LSB galaxies were based on Hi observations with limited spatial resolution, and the data were in part affected by beam smearing. Even though McGaugh & de Blok (1998) showed that the steep signatures of the rotation curves implied by CDM could not be hidden by any reasonable amount of beam smearing, there were later suggestions that the observed data *were* consistent with the CDM predictions (van den Bosch et al., 2000; van den Bosch & Swaters, 2000) if proper beam smearing corrections were applied. Swaters, Madore & Trewhella (2000) published optical rotation curves of five LSB galaxies from the de Blok, McGaugh & van der Hulst (1996) sample and found that in several cases the inner rotation curve slopes were somewhat steeper than found from the Hi curves. It is thus conceivable that the data could be reconciled with CDM models once beam smearing corrections are properly taken into account. This is not borne out by improved data, as we show in this Letter.

## 2 The Data

McGaugh, Rubin & de Blok (2001) (MRdB), de Blok, McGaugh & Rubin (2001) (dBMR) and de Blok & Bosma (2001) present high-resolution hybrid H$`\alpha `$/Hi rotation curves of LSB galaxies, and show that a large fraction of them are characterised by a slow, solid-body-like rise. They compare (cuspy) CDM NFW halos (Navarro, Frenk & White, 1996) and pseudo-isothermal (core-dominated) halos (e.g. Begeman, Broeils & Sanders, 1991) with the data and find that the pseudo-isothermal model is statistically a much better fit to the data. Most NFW model fits suffer from systematic effects resulting from the fitting program trying to reconcile $`v(r)r`$ data with a steep $`v(r)r^{1/2}`$ prediction. The discrepancies between Hi and optical data are found to be less severe than initially suggested in Swaters, Madore & Trewhella (2000): there is reasonable agreement between new optical data and the old Hi curves in the majority of the cases (MRdB). The main conclusion drawn from the new data is that LSB galaxies cannot be well fitted by CDM rotation curves. Here we derive the mass density profiles that are required to give rise to the observed rotation curves. These can be compared directly to those predicted by theoretical models.

We use the sample of LSB and dwarf galaxy rotation curves described in dBMR, MRdB and de Blok & Bosma (2001). This sample includes the five LSB galaxies originally presented in Swaters, Madore & Trewhella (2000) that were re-derived in dBMR. To this we have added the 4 LSB galaxies from the Ursa Major sample of Verheijen (1997) with $`\mu _0(B)>22.0`$ mag arcsec-2 with reliable Hi rotation curves.

## 3 Results

In principle, one can invert an observed rotation curve to determine the parent mass distribution. In practice, this procedure can be unstable for thin disks (Binney & Tremaine, 1987) (but see Sackett 1997). For LSB galaxies, the disk component is negligible since the potential is dominated by the dark matter halo. We therefore invert the observed rotation curves assuming a spherical mass distribution, which is a straightforward and robust procedure. From $`^2\mathrm{\Phi }=4\pi G\rho `$ and $`\mathrm{\Phi }=GM/r`$ one can derive the mass density $`\rho (r)`$:

$$4\pi G\rho (r)=2\frac{v}{r}\frac{v}{r}+\left(\frac{v}{r}\right)^2,$$
(1)
where $`v`$ is the rotation velocity and $`r`$ is the radius.

Implicit in this procedures are assumptions: that these galaxies are all dark matter dominated, that the gas motions are circular in a planar disk and that the spectrum samples the nucleus and major axis. We also assume that the galaxies are symmetric. (The latter is a good assumption: MRdB, dBMR and de Blok & Bosma (2001) omitted asymmetric galaxies from the samples, whereas the effects of any mild, residual asymmetries were incorporated in the errorbars.)

Here we ignore mass contributions of the stellar and gas components. As long as these do not dominate the dynamics, as is the case in LSB galaxies, this “minimum disk” assumption is a good one. A minimum disk also produces an upper limit on the steepness of the halo profile as inclusion of gas and stars will tend to flatten the derived slopes. See also dBMR and de Blok & McGaugh (1997). In a forthcoming paper we show the results for non-minimum disk hypotheses, which are not substantially different from the ones derived here.

Figure Mass Density Profiles of LSB Galaxies shows the derived mass density profiles. Over-plotted are the profiles of the best fitting pseudo-isothermal and NFW halo as listed in dBMR and de Blok & Bosma (2001). The errors on the data points are derived by rigorously carrying through the errors in the rotation curve data points. The arrows indicate the size of the seeing disk for each galaxy. (For the 4 Verheijen 1997 galaxies we indicate the radio beam size.)

The shape of the mass density profiles can generally be characterised by two components: an outer one with an isothermal slope of $`2`$ and a more shallow one in the inner parts. After determining the “break-radius” where the slope changes most rapidly, we determine the slope of the inner component using a weighted least-squares fit. Table 1 collects the values of the inner slope. The range over which the power-law is fitted is indicated by the range over which the fit is drawn in Fig. Mass Density Profiles of LSB Galaxies. The uncertainty $`\mathrm{\Delta }\alpha `$ is determined by re-measuring the slope twice, once by including the first data point outside the break-radius, and once by omitting the data point at the break-radius. The maximum difference between these two values and the original slope is adopted as the uncertainty. The isothermal mass model generally follows the derived mass profile very well. The NFW models usually fail dramatically in the inner parts. For the galaxies, the signature of the shallow inner slope is usually already present well outside the seeing disk. Optical “beam smearing” can therefore not cause these shallow inner slopes.

Minimum disk does provide relevant limits on the inner slopes as is particularly well illustrated by UGC 6614. This is a bulge-dominated giant LSB galaxy, and we expect the bulge to contribute significantly to the dynamics at small radius for any plausible stellar mass-to-light ratio. However, even in the minimum disk case we find a shallow slope $`\alpha =0.3`$. This means that for any non-minimum disk situation the dark matter must be depressed even further away from the NFW case (i.e. the true slope must be even flatter).

Values for the inner slopes (Fig. Mass Density Profiles of LSB Galaxies) are asymmetric. In Fig. Mass Density Profiles of LSB Galaxies we have distinguished between galaxies where the profile is well-resolved and those where the turn-over in the profile occurs at or within the seeing radius. The most unresolved galaxies tend to have the most negative values of $`\alpha `$. The resolved galaxies define a well-determined peak at $`\alpha =0.2\pm 0.2`$ that is inconsistent with CDM predictions that $`\alpha =1.5`$.

One possible point of concern is the wing of steeper slopes extending to $`\alpha =1.8`$, where the extreme values originate from the Verheijen (1997) UMa LSB profiles. Does this mean that despite all of the above, there are LSB galaxies that *can* be well fitted with CDM models? And does the tendency of the UMa LSB galaxies to have steep slopes indicate systematic effects in the new LSB data? The answer to both questions is negative, as the following analysis shows.

Table 1 lists the radius in kpc, $`r_{\mathrm{in}}`$, of the innermost data point of each profile. For the LSB sample generally $`r_{\mathrm{in}}<1`$ kpc. For the UMa galaxies we find larger $`r_{\mathrm{in}}`$: three of the four have $`r_{\mathrm{in}}=1.5`$ kpc. In Fig. Mass Density Profiles of LSB Galaxies, we plot the values of $`r_{\mathrm{in}}`$ versus the inner slope $`\alpha `$. Also drawn are the logarithmic slopes as a function of radius for pseudo-isothermal halos with core-radii $`R_C=0.5,1,2`$ kpc, as well as a NFW model and a CDM $`r^{1.5}`$ model (Moore et al. 1999), both of the latter converging to a slope $`\alpha =3`$ in the (far) outer parts. These two models are chosen to have parameters $`c=8`$ and $`V_{200}=100`$ km s-1 to approximately match the four UMa galaxies. However, this choice is not critical.

Galaxies with small values of $`r_{\mathrm{in}}`$ ($`0.15`$ kpc) show clear evidence of a core ($`\alpha 0`$), whereas galaxies with larger values of $`r_{\mathrm{in}}`$ exhibit steeper slopes. Fig. Mass Density Profiles of LSB Galaxies shows that distribution is consistent with an isothermal halo with a core radius of a few kpc, whereas the NFW and CDM models do not match the data at all.

Hence, only galaxies with small values of $`r_{\mathrm{in}}`$ measure the core. Larger values sample a transition zone where the slope is changing from $`\alpha =0`$ (center) to $`\alpha =2`$ (outer isothermal regions). In the zone between $`1`$ and $`10`$ kpc the slopes of the pseudo-isothermal and CDM models are approximately equal, so large values of $`r_{\mathrm{in}}`$ might erroneously lead to the conclusion that measured slopes are consistent with CDM. The four UMa galaxies (and some LSB galaxies) have $`r_{\mathrm{in}}`$ in this transition zone and thus show steep slopes. We have modelled the beam smearing or seeing effects potentially present in the optical data (discussed in a forthcoming paper) and find that we can strongly exclude the possibility that these affect the results down to resolutions of $`0.1^{\prime \prime }`$. We thus predict that higher spatial resolution data (with smaller values of $`r_{\mathrm{in}}`$) will also detect cores in the less well-resolved galaxies.

Similar arguments apply to the beam smearing corrected Hi curves in van den Bosch et al. (1999) and van den Bosch & Swaters (2000). With values $`r_{\mathrm{in}}1`$ kpc these data trace *not* the inner slope but instead the steep slope at the turnover of the constant-density core.

## 4 Conclusions

Mass density profiles of LSB galaxies exhibit inner slopes that are best described by a power-law $`\rho (r)r^\alpha `$ with $`\alpha =0.2\pm 0.2`$. This result implies that halos of LSB galaxies are dominated by cores. This result is inconsistent with the value $`\alpha =1.5`$ predicted by CDM models. The steep slopes found for some LSB galaxies arise when the innermost data point is sampling the transition region between core and outer $`\alpha =2`$ isothermal region, not the core itself. Our data are not consistent with CDM predictions, but suggest that LSB galaxies have halos which contain cores of radii of order 1 kpc.