text:
# Time Lags in Low Mass X-ray Binaries

## 1. Introduction

TLs have been reported so far from both accreting black holes (BHs, Cyg X-1, GX339-4) and NSs (4U0614+09, 4U1705-44, Ford et al. 1999). The amplitude of these TLs is typically 20 ms at 1 Hz for NSs and BHs (Ford et al. 1999). The origin of these TLs is currently under debate. Uniform comptonization models predicting constant TLs are however ruled out by the data. Comptonization in a non-uniform medium might account for the observed TLs (Kazanas et al. 1997), but in this case, their magnitude implies that the Comptonizing cloud extends to very large radii, which in turn poses an energy problem. Recently, Poutanen and Fabian (1999) have proposed a “magnetic flare” model that could reproduce the observed PDS and TL magnitude without requiring a large size region.

Using proprietary and archival RXTE/PCA data, we have initiated a systematic study of TLs from NSs, with the aim of determining whether TLs were associated with a singular spectral and/or timing state. We characterize the spectral state of a source using color-color diagrams (CCDs). For the timing state, we compute the PDS over a broad X-ray energy band (2-40 keV). To start with, we have selected two NSs that are always in a low state (1E1724-3045, GS1826-238, Barret et al. 2000), and two that are more variable; undergoing occasionally low states (4U1705-44, 4U1728-34). For the latter sources, we have sampled both their low and high states. Their low state PDS are typical of “island” state PDS, being characterized by a flat top below $`\nu _{\mathrm{Break}}`$ and a broad QPO-like feature in the declining part of the PDS above $`\nu _{\mathrm{Break}}`$ (at $`\nu _{\mathrm{QPO}}`$). TLs were computed using the techniques described in Nowak et al. (1999a) between the 7-40 keV hard and 2-7 keV soft photons.

## 2. Time Lags versus spectral/timing state

Figure 1 shows the CCD of the four sources. Fig. 2 shows the corresponding PDS (top) and TL spectra (bottom). TLs were detected only from 4U1728-34 and 4U1705-44, and in their lowest/hardest intensity states (namely their “island” state). This is the first report of TLs from 4U1728-34. No TLs were detected in their high states, with upper limits of 0.1-0.01 seconds between 1 and 10 Hz; i.e. a factor of 10 larger than the values detected during their low states. No TLs were detected from the two steady low state sources, and the upper limits we derived are lower than the observed values for 4U1728-34 and 4U1705-44, indicating that the non detection of similar TLs is not due to a lack of sensitivity.

Looking at Fig. 1 and 2, a few points can be drawn. First, TLs are not associated with a singular spectral state; 4U1728-34 and 1E1724-3045 occupy the same region of the CCD, and only the former shows TLs. Second, although the overall shape of their PDS is broadly similar, there is one noticeable difference that shows up very clearly in the $`\nu F\nu `$ representation of the PDSs; that TLs are associated with a timing state in which the whole PDS is shifted towards high frequencies ($`\nu _{\mathrm{Break}}`$ and $`\nu _{\mathrm{QPO}}`$ are a factor of 10 larger for 4U1728-34 and 4U1705-44 than for 1E1724-3045 and GS1826-238). Third, when $`\nu _{\mathrm{Break}}`$ and $`\nu _{\mathrm{QPO}}`$ are high, TLs are significantly detected at frequencies between $`\nu _{\mathrm{Break}}`$ and up to or slightly above $`\nu _{\mathrm{QPO}}`$. Fourth, although of lower significance than the effect observed in the two BHs Cyg X-1 and GX339-4 (Nowak et al. 1999a,b), there is an indication that the TL decreases with frequency, especially for 4U1728-34. Finally, TLs do not depend upon the intensity of the aperiodic variability, as the four PDSs of Fig 2. have comparable integrated RMS (ranging from 17 to 25% in the 2-40 keV band).

## 3. Conclusions

We have discovered TLs in the low state of 4U1728-34, and confirmed the previous detection of TLs from 4U1705-44 (Ford et al. 1999). In our attempt to associate the presence of TLs with a singular spectral state, we have shown that for the same spectral state (same position in the CCD) some sources do have TLs whereas some do not. The presence of TLs must thus be tracked somewhere else.

Despite having a limited sample, we have found that TLs are detected when the characteristic low state PDS shows excess power at high frequencies (with a $`\nu _{\mathrm{Break}}`$ $`1`$ Hz instead of $`0.10.2`$ Hz). Recently, Pottschmidt et al. (2000) showed that the shot relaxation time in the hard state of Cyg X-1 (scaling as $`\nu _{\mathrm{Break}}`$<sup>-1</sup>) anticorrelates with the TL amplitude. We cannot test the presence of this effect within the data set used here. However, if the same anticorrelation applies to NSs, it could provide an explanation for the non detection of smaller TLs from the sources with the lowest $`\nu _{\mathrm{Break}}`$. If $`\nu _{\mathrm{Break}}`$ and $`\nu _{\mathrm{QPO}}`$ are somehow related to the position of the inner disk radius within the corona, as possibly suggested by recent observations (Revnivtsev et al. 2000), the closer the disk gets to the central object (i.e the deeper inside the corona), and the larger the TLs are. In this picture, it is interesting to note that this subtle change in the accretion geometry is not reflected in the CCD, implying similar parameters for the Comptonizing cloud and input seed photon energy for all sources.

To conclude, the above study, although not addressing the origin of the TL itself has provided some hints about the conditions under which they are produced. Further work is needed to determine whether TLs are indeed always associated with the peculiar timing state described in this paper. Confirming this result would be extremely valuable for constraining any theoretical models attempting to reproduce the TL properties of these systems.

We thank E. Ford and A. Zdziarski for helpful comments on this paper.

## REFERENCES

Barret, D. et al., 2000, ApJ, 533, 329
Ford, E. C. et al., 1999, ApJL, 512, L31
Kazanas, D., Hua, X. & Titarchuk, L., 1997, ApJ, 480, 280
Nowak, M. et al., 1999a, ApJ, 510, 874
Nowak, M. et al., 1999b, ApJ, 517, 355
Poutanen, J. & Fabian, A., 1999, MNRAS, 306, L31
Pottschmidt K. et al., 2000, A&A, 357, L17
Revnivtsev, M., Gilfanov, M., Churazov, E., 2000, A&A submitted, (astro-ph/9910423)


format: text

text:
# Field-induced long-range order in the 𝑆=1 antiferromagnetic chain

\[

## Abstract

The quasi-one dimensional $`S=1`$ antiferromagnet in magnetic field $`H`$ is investigated with the exact diagonalization of finite chains and the mean field approximation for the interchain interaction. In the presence of the single-ion anisotropy $`D`$, the full phase diagram in the $`HT`$ plane is presented for $`HD`$ and $`HD`$. The shape of the field-induced long-range ordered phase is revealed to be quite different between the two cases, as observed in the recent experiment of NDMAP. The estimated ratio of the interchain and intrachain couplings of NDMAP ($`J^{}/J10^3`$) is consistent with the neutron scattering measurement.

PACS Numbers: 75.10.Jm, 75.30.Kz, 75.40.Cx, 75.50.Ee

\]

Since Haldane conjectured the presence of the spin excitation gap in the one-dimensional integer-$`S`$ Heisenberg antiferromagnets, the $`S=1`$ antiferromagnetic chain has attracted a lot of interest. Many theoretical and numerical works indicated the disordered ground state of the system. The mean field approximation for the interchain interaction combined with the exact diagonalization of finite chains revealed that such a disordered ground state should be realized even in the quasi-one-dimensional $`S=1`$ antiferromagnet if the ratio of the interchain and intrachain couplings holds $`zJ^{}/J0.05`$ where $`z`$ is the number of adjacent chains. The criterion was satisfied for the material Ni(C2H8N2)2NO2(ClO4), abbreviated NENP. In fact the Haldane gap had already been observed in the temperature dependence of the magnetic susceptibility and the high-field magnetization measurements of NENP.

The effective field theory and the size scaling analysis with the finite chain calculation indicated that the field-induced transition occurs from the disordered ground state to the gapless Tomonaga-Luttinger liquid phase at some critical external field $`H_{c1}`$ where the Haldane gap vanishes. With the interchain interaction, the gapless phase would become the ordered phase. Thus the field-induced antiferromagnetic long-range order is expected to appear at $`H_{c1}`$ in the quasi-one-dimensional $`S=1`$ system. NENP, however, didn’t exhibit such a field-induced spontaneous symmetry breaking, because the staggered moment is forced to appear even by a smaller external field through the alternation of the $`g`$-tensors. Recently the field-induced transition to the long-range order was discovered in the high-field magnetization measurement of novel quasi-one-dimensional $`S=1`$ systems Ni(C5H14N2)2N3(ClO4) and Ni(C5H14N2)2N3(PF6), abbreviated NDMAZ and NDMAP, respectively. The magnetic susceptibility measurement indicated the presence of the easy-plane single-ion anisotropy $`D`$ for both systems and the value was reported as $`D0.3J`$ for NDMAP. The experimental phase diagram of both systems in the $`HT`$ plane exhibited a quite different shape of the ordered phase between the two cases when the external field $`H`$ is (a)parallel and (b)perpendicular to the principal axis of $`D`$. This is probably because the universality of the symmetry breaking is different between the two cases. Thus it would be interesting to present the theoretical $`HT`$ phase diagram. In this paper, using the exact diagonalization of finite chains and the mean field approximation for the interchain interaction, the phase diagrams for $`HD`$ and $`HD`$ are presented, to compare with the experimental ones and give some suggestions for future experiments.

In order to investigate the magnetization process of the quasi-one-dimensional $`S=1`$ antiferromagnet with the planar single-ion anisotropy, we consider the Hamiltonian,

$``$ $`=`$ $`J{\displaystyle \underset{j}{}}{\displaystyle \underset{i}{}}𝐒_{j,i}𝐒_{j,i+1}+J^{}{\displaystyle \underset{j,j^{}}{}}{\displaystyle \underset{i}{}}𝐒_{j,i}𝐒_{j^{},i}`$ (1)
$`+`$ $`D{\displaystyle \underset{j}{}}{\displaystyle \underset{i}{}}(S_{j,i}^z)^2H{\displaystyle \underset{j}{}}{\displaystyle \underset{i}{}}S_{j,i}^\alpha ,`$ (2)
where $`J`$ and $`J^{}`$ are the intrachain and interchain couplings, respectively, $`j`$ is the index of a chain and $`_{j,j^{}}`$ means the sum over all the pairs of nearest-neighbor chains. We set $`J=1`$ and use the units such as $`g\mu _B=1`$ throughout the paper. Only the $`XY`$-like anisotropy $`D>0`$ is investigated, to compare with the experimental results of NDMAP. To consider the two cases (a)$`HD`$ and (b)$`HD`$, we set $`\alpha =z`$ and $`y`$, respectively. When the antiferromagnetic order is possibly induced by the external field $`H`$, $`xy`$-plane is an easy plane for (a) while $`x`$-axis is an easy axis for (b). Thus the universality class of the symmetry breaking should be different between the two cases; the former corresponds to the $`XY`$ model ($`U(1)`$) while the latter the Ising one ($`Z_2`$).

In the magnetization process at low temperatures two phase transitions are expected to occur at critical fields $`H_{c1}(T)`$ and $`H_{c2}(T)`$. At $`T=0`$ the Haldane gap vanishes at $`H_{c1}(0)`$ and the magnetization is saturated at $`H_{c2}(0)`$. The antiferromagnetic long-range order in the vertical direction to $`H`$ should appear between $`H_{c1}(T)`$ and $`H_{c2}(T)`$. Based on the mean field approximation for the interchain interaction, the critical field $`H_{c1}(T)`$, as well as $`H_{c2}(T)`$ is determined by the solution $`H`$ of the equation

$`\chi _{\mathrm{st}(1\mathrm{D})}^\rho (T,H)={\displaystyle \frac{1}{zJ^{}}},`$ (3)
where $`\chi _{\mathrm{st}(1\mathrm{D})}^\rho (T,H)`$ is the staggered susceptibility of the one-dimensional system ($`J^{}=0`$) and $`\rho `$ (=$`x`$, $`y`$ or $`z`$) specifies the easy component of the antiferromagnetic order. To consider the two cases (a) and (b) we set $`(\alpha ,\rho )=(z,x)`$ and $`(y,x)`$, respectively. Using the numerical diagonalization analysis, we obtain the all the eigenvalues and eigenvectors of finite chains and calculate $`\chi _{\mathrm{st}(1\mathrm{D})}^\rho (T,H)`$. Since the rotational symmetry doesn’t remain in any directions for $`HD`$, the available chain length is at most $`L=8`$ under the periodic boundary condition because of limitation of the memory size. At finite temperature, however, the correlation length is much smaller than that at $`T=0`$ ($`\xi `$5 to 6 lattice constants for $`H=D=0`$). Thus even the result for $`L=8`$ is expected to indicate some qualitative features of the bulk system. In fact the small oscillation in the $`H`$ dependence of $`\chi _{\mathrm{st}(1\mathrm{D})}^\rho (T,H)`$ due to the finite size effect only appears at $`T0.1`$ for $`L=8`$. Based on obtained $`H_{c1}(T)`$ and $`H_{c2}(T)`$ from the calculation for $`L=8`$, we discuss on the properties of the phase diagram in the $`HT`$-plane in the following.

To investigate the finite size effect, we show in Fig. 1 $`H_{c1}(T)`$ and $`H_{c2}(T)`$ of the isotropic system ($`D=0`$) with $`zJ^{}=0.05`$ determined by the equation (3) for $`L=`$8, 6 and 4 in the $`HT`$-plane. The lower and upper boundaries of the antiferromagnetic ordered phase are given by $`H_{c1}(T)`$ and $`H_{c2}(T)`$, respectively. Those lines also give the $`H`$ dependence of the Néel temperature $`T_\mathrm{N}`$. Fig. 1 indicates that estimated $`T_\mathrm{N}`$ is higher for larger $`L`$, which suggests that the finite-chain analysis tends to underestimate the ordered phase. On the other hand, the mean field approximation generally overestimates it. Thus on the present analysis the finite size effect and the fluctuation around the mean field are expected to suppress the errors due to each other.

In the following analysis we put $`D=0.3`$, realistic for NDMAP. The phase diagrams in the $`HT`$-plane obtained from $`\chi _{\mathrm{st}(1\mathrm{D})}^\rho (T,H)`$ for $`L=8`$ are shown in Figs. 2(a) and (b), for $`HD`$ and $`HD`$, respectively. The phase boundaries for $`zJ^{}=`$0.05, 0.04 and 0.03 are given in each case. Although $`H_{c2}(T)`$ exhibits an oscillation due to finite size effect, Fig. 2(a) clarifies an interesting feature; $`H_{c1}(T)`$ and $`H_{c2}(T)`$ are quite asymmetric for $`HD`$. In contrast, $`H`$ dependence of $`T_\mathrm{N}`$ is close to a symmetric bell shape for $`HD`$. According to our calculation for other values of $`D`$, the difference in the shape of $`T_\mathrm{N}`$ becomes clearer for larger $`D`$. With increasing $`D`$, $`T_\mathrm{N}`$ becomes more asymmetric for $`HD`$, while closer to a completely symmetric bell shape for $`HD`$. The asymmetric behavior of $`T_\mathrm{N}`$ for $`HD`$ is consistent with the $`H`$ dependence of the critical exponent $`\eta `$ in the transverse spin correlation function $`S_0^xS_r^x(1)^rr^\eta `$ of the isotropic $`S=1`$ chain at $`T=0`$. In the ground-state magnetization process of the system the antiferromagnetic spin correlation measured by $`\eta `$ rapidly grows from $`H_{c1}`$ to the maximum while slowly decays toward $`H_{c2}`$. Thus such an asymmetric behavior is expected to be characteristic of the system with $`U(1)`$ symmetry. On the other hand, the case (b) with only $`Z_2`$ symmetry is supposed to be close to a classical spin system where the spin wave excitation spectrum around $`m1`$ is similar to that around $`m0`$.

The phase diagrams for $`HD`$ and $`HD`$ ($`zJ^{}=0.03`$) are shown together in Fig. 3. They exhibit several characteristic behaviors observed in the experiments of NDMAP and NDMAZ; (i)At lower temperatures $`H_{c1}(T)`$ is smaller for $`HD`$. (ii)At intermediate $`H`$ $`T_\mathrm{N}`$ is higher for $`HD`$. (iii)The two $`T_\mathrm{N}`$ curves have two intersections, although the upper one has not been observed because such a strong magnetic field cannot be realized. (i) and (ii) imply that the ordered state resulting from the $`U(1)`$ ($`Z_2`$) symmetry breaking is less (more) sensitive to the external magnetic field, while less (more) stable against the thermal fluctuation. Our calculation for various values of $`D`$ indicates that the intersection of two $`T_\mathrm{N}`$ curves appears only for $`D>0`$. Thus the observed intersection is also an evidence of the positive $`D`$ for NDMAP and NDMAZ. It is consistent with the sign of $`D`$ determined by the susceptibility measurements.

The qualitative consistency of the present results with the experiments of NDMAP encourages us to estimate the strength of the interchain coupling $`J^{}`$ of NDMAP by fitting the calculated phase boundaries to the measured one. The calculated phase boundaries ($`H_{c1}(T)`$ and $`H_{c2}(T)`$) for $`zJ^{}=`$0.010, 0.012 and 0.014 are shown together with the experimental results of NDMAP in Figs. 4 (a)$`HD`$ and (b)$`HD`$, where dashed curves are guides for the eyes. We set $`D=0.3`$. The temperature is normalized as $`J`$=30K and the magnetic field is normalized such that $`H_{c1}(0)`$ coincides with the experimental results. In Fig. 4 (a) it is difficult to fit the calculated $`T_\mathrm{N}`$ curve to the experimental one, because the finite size effect is still large. Thus we determine $`zJ^{}`$ by the agreement of the maximum value of $`T_\mathrm{N}`$ as $`zJ^{}=0.012\pm 0.002`$. On the other hand, the fit of the calculated $`T_\mathrm{N}`$ to the measured one is much better for $`HD`$ in Fig. 4(b). It leads to the same conclusion $`zJ^{}=0.012\pm 0.002`$. However, the error of the estimation might be much larger, because we neglect the finite size correction and the fluctuation around the mean field of the interchain couplings. Since the number of the adjacent chains is usually more than 2, the result $`zJ^{}/J0.012`$ leads to $`J^{}/J10^3`$. It is consistent with the most dominant interchain coupling estimated from the recent neutron scattering experiment. Using the same mean field approximation with the transfer matrix calculation of the chain, a little smaller value of $`zJ^{}/J`$ was estimated. It is a reasonable difference, because the size correction of the present analysis would shift the result to a smaller value. However, it is difficult to determine which result is better, because some analyses beyond the mean field approximation should modify it into a larger value. Fig. 4 (b) suggests that the maximum of $`T_\mathrm{N}`$ of NDMAP could be detected with a little larger external field even for $`HD`$. It would also enable us to perform a more precise estimation of $`J^{}`$ and confirm the difference in the shape of $`T_\mathrm{N}`$ between $`HD`$ and $`HD`$.

Recently a sign of the field-induced long-range order was also observed in the NMR measurement of another quasi-one-dimensional $`S=1`$ system (CH3)4NNi(NO2)3, abbreviated TMNIN. Since $`J`$ is smaller ($`J12`$K), the full magnetization curve has already been obtained. Thus the full phase diagram in the $`HT`$-plane will possibly be determined in the near future. It would clarify more detailed features of the field-induced long-range order.

In summary, based on the mean field approximation for the interchain interaction and numerical diagonalization of finite chains, we obtained the phase diagram in the $`HT`$-plane of the quasi-one-dimensional $`S=1`$ antiferromagnet. The results revealed that the shape of the field-induced ordered phase is quite different between $`HD`$ and $`HD`$. Fitting the calculated phase boundary to the experimental result lead to the estimation of the interchain coupling of NDMAP as $`zJ^{}/J0.012`$, which well agrees with the estimation by the recent neutron scattering measurement ($`J^{}/J10^3`$).

We thank T. Goto, Z. Honda, K. Katsumata and I. Harada for fruitful discussions. The computation in this work has been done using the facilities of the Supercomputer Center, Institute for Solid State Physics, University of Tokyo. This research was supported in part by Grant-in-Aid for the Scientific Research Fund from the Ministry of Education, Science, Sports and Culture (11440103).


text:
# Are passive protostellar disks stable to self-shadowing?

## 1 Introduction

It has become widely accepted that the infrared excess observed in the spectral energy distribution of many T-Tauri stars is caused by thermal emission from a circumstellar dusty disk (Adams & Shu, Adams and Shu (1986); Adams, Lada & Shu, Adams et al. (1987)). However, the spectra predicted by theoretical models of both accretion-powered disks (Lynden-Bell & Pringle, Lynden-Bell and Pringle (1974)) and flat reprocessing disks (Friedjung, Friedjung (1985)) have an infrared spectrum that goes as $`\lambda F_\lambda \lambda ^{1.33}`$. This is steeper than the $`\lambda F_\lambda \lambda ^{0.5\text{to}1.0}`$ typically observed from T-Tauri stars (Rucinski, Rucinski (1985); Rydgren & Zak, Rydgren and Zak (1987)). A natural way to explain the relatively warm dust at large radii (which is needed to increase the slope of the model spectrum), is to invoke a flaring disk geometry, in which $`H/RR^\gamma `$, with $`\gamma >0`$ (Kenyon & Hartmann, Kenyon and Hartmann (1987)). Because of the angle between the disk surface and the star, the disk intercepts considerable stellar flux even at large radii, and therefore acquires a shallower temperature profile than the $`TR^{3/4}`$ found from the accreting and/or flat disk models. The flaring index $`\gamma `$ can be determined self-consistently by equating the intercepted flux with the emitted blackbody flux, under the assumption that the disk is vertically isothermal. These models yield $`\gamma =2/7`$ (Chiang & Goldreich, Chiang and Goldreich (1997), henceforth CG97; D’Alessio et al. D’Alessio et al. (1999), henceforth DCHLC99). Apart from producing the flatter SEDs, it was recognized that the optically thin surface layers of the disk will have much higher temperatures than the interior (Calvet et al.Calvet et al. (1991); Malbet & Bertout Malbet and Bertout (1991); CG97). In addition to producing an extra component in the SED, this optically thin layer is also able to explain certain dust features in emission which would otherwise be expected in absorption. A more direct piece of evidence for the flaring nature of protostellar disks was provided by the HST images of HH30 (Burrows et al. Burrows et al. (1996)).

Despite the successes of the flaring disk model in explaining observations of T-Tauri stars, there are still some theoretical issues that remain to be addressed. In this paper we address two questions: are the flaring disk models stable against perturbations in vertical height and temperature, and are they unique solutions to which a time-dependent evolution of the passive disk settles down? The issue of stability has been addressed before by DCHLC99. They considered perturbations in the temperature of the disk, and found that they damp out as they propagate inward. However, their analysis is based on the assumption that the disk vertical height quickly responds to changes in the temperature, or in other words that the vertical sound crossing time $`t_{\text{dyn}}`$ is much smaller than the cooling time $`t_{\text{therm}}`$. This is true for the inner regions of the disk, but breaks down at large radii. In this Letter we study the other extreme case: the case of $`t_{\text{dyn}}t_{\text{therm}}`$, which may occur in the outer regions of the disk. We start from the same equations as DCHLC99 and CG97, but replace the static equation for the vertical disk height with a simplified one-zone dynamic model. We derive the pertubation equations and dispersion relations from them, and find that under these circumstances the disk becomes unstable.

We will start the analysis in Section 2 with a description of the family of solutions of a static irradiated optically thick disk. The flared disk model mentioned above is just one of these solutions. The other solutions either diverge at, or are conical ($`HR`$) beyond some radius. In Section 3 we will discuss the stability under the two extreme conditions ($`t_{\text{therm}}t_{\text{dyn}}`$ and $`t_{\text{therm}}t_{\text{dyn}}`$). In Section 4 we discuss the consequences of this instability, and when it is expected to show up.

## 2 Static equations for an irradiated disk

A disk with flaring angle $`\alpha `$ intercepts a flux of $`F_{\text{irr}}=\alpha (R_{}/R)^2\sigma _\text{R}T_{}^4\text{erg}\text{cm}^2\text{s}^1`$. The flaring angle $`\alpha `$ is defined as the angle between the grazing incident stellar light and the surface of the flaring disk, and is assumed to be small. If we assume the disk to be vertically isothermal, the emitted flux is $`F_{\text{emit}}=\sigma _\text{R}T_\text{e}^4`$, where $`T_\text{e}`$ is the temperature of the disk. Equating $`F_{\text{irr}}`$ with $`F_{\text{emit}}`$ yields

$$T_\text{e}=\alpha ^{1/4}\left(\frac{R_{}}{R}\right)^{1/2}T_{}.$$
(1)
When $`RR_{}`$ the flaring angle is given by

$$\alpha =R\frac{d}{dR}\left(\frac{H_\text{s}}{R}\right),$$
(2)
where $`H_\text{s}`$ is the height of the surface of the disk above the midplane. For a static solution the vertical pressure balance equation for an isothermal disk yields a Gaussian density profile of the form

$$\rho =\frac{\mathrm{\Sigma }}{H_\text{p}\sqrt{2\pi }}\mathrm{exp}\left(\frac{z^2}{2H_\text{p}^2}\right),$$
(3)
where $`\mathrm{\Sigma }`$ is the surface density, and $`H_\text{p}`$ is the pressure scale height

$$H_\text{p}=\sqrt{\frac{kT_\text{e}R^3}{\mu _\text{g}m_\text{u}GM_{}}},$$
(4)
where $`\mu _\text{g}`$ is the mean molecular weight and $`m_\text{u}`$ is the unit atomic mass. The surface height $`H_\text{s}`$ and the pressure scale height $`H_\text{p}`$ are related by a number of the order of a few, dependent on the opacity. We assume this to be approximately a constant $`\chi H_\text{s}/H_\text{p}`$. CG97 take $`\chi 4`$ for their analytic calculation of the disk structure. If we eliminate $`T_\text{e}`$ from the above set of equations we arrive at

$$\frac{1}{\chi }\left(\frac{\mu _\text{g}m_\text{u}GM_{}}{kT_{}R_{}}\right)^4\left(\frac{R_{}}{R}\right)^2\left(\frac{H_\text{p}}{R}\right)^8=R\frac{d}{dR}\left(\frac{H_\text{p}}{R}\right).$$
(5)
Defining the dimensionless quantities $`rR/R_{}`$ and $`hH_\text{p}/R`$, one arrives at

$$r\frac{dh(r)}{dr}=\frac{C}{r^2}h(r)^8,$$
(6)
where $`C`$ is

$$C=\frac{1}{\chi }\left(\frac{\mu _\text{g}m_\text{u}GM_{}}{kT_{}R_{}}\right)^4.$$
(7)
This ordinary differential equation (ODE) has a 1-parameter family of solutions:

$$h(r)=\left(\frac{7C}{2}\frac{1}{r^2}+K\right)^{1/7},$$
(8)
where $`K`$ is an arbitrary real constant. For $`K<0`$ the solutions diverge at $`r=\sqrt{7C/2K}`$, for $`K>0`$ the solutions become constant beyond $`r=\sqrt{7C/2K}`$, and for $`K=0`$ one obtains a powerlaw solution. These three classes of solutions are shown in Fig.(1).

The power law solution is given by

$$h(r)=\left(\frac{2}{7C}\right)^{1/7}r^{2/7},$$
(9)
which is the solution given by CG97. All the other solutions asymptotically tend to this power law solution, Eq.(9), as $`r0`$. But as $`r`$ goes to infinity the solutions deviate from this power law at some radius, and either diverge to $`h\mathrm{}`$ at a particular radius, or asymptotically approach a conical geometry $`h(r)\text{const}`$, i.e. $`H_\text{s}R`$.

The direction in which the solutions deviate from each other is interesting. Since the inner regions of the disk (small radii) cast their shadow over the outer regions (large radii), the “causality” of the system is pointing outwards. A change in the structure of the disk at small radii has its repercussions on the disk at large radii, but not vice versa. In principle this means that the (numerical) integration of the ODE (Eq.5) must be done from inside-to-outside (i.e. in causal direction). However, in this direction the solutions tend to strongly deviate from each other. A slight change in integration constant at the start of the integration causes a very large difference at large radii. For stable integration of the ODE (Eq.5) one should place the boundary condition at the outside and integrate inwards. Though mathematically correct, it is physically meaningless, and may hint to an intrinsic instability of the time-dependent equations.

## 3 Time-dependent equations for irradiated disk

The stability of the flaring (power law) solution can be studied analytically using time-dependent linear perturbation analysis. We study two regimes, one of which has been studied before by DCHLC99.

### 3.1 Temperature perturbations at hydrostatic equilibrium

In the analysis of DCHLC99 the cooling time scale of the disk was assumed to be much larger than the orbital time scale, which is true for the inner regions of the disk. We may then assume that the disk will always be in hydrostatic equilibrium, and we follow perturbations of the temperature. Following DCHLC99 we write

$$\mathrm{\Gamma }\frac{dT_\text{e}}{dt}=F_{\text{irr}}F_{\text{emit}},$$
(10)
where $`\mathrm{\Gamma }`$ is the thermal heat capacity. After substituting Eqs.(2,4,7), the definitions of $`F_{\text{irr}}`$ and $`F_{\text{emit}}`$, and defining the perturbation $`T_\text{e}(r,t)=T_\text{e}^0(r)(1+\stackrel{~}{\psi }(r,t)r^2)`$, where $`T_\text{e}^0(r)`$ is the temperature belonging to the static flaring disk solution Eq.(9), one finds the following perturbation equation (see DCHLC99 for details):

$$\frac{\stackrel{~}{\psi }}{\overline{t}}=\frac{7}{4}\left(\frac{T_\text{e}^0}{T_{}}\right)^3r\frac{\stackrel{~}{\psi }}{r},$$
(11)
where $`\overline{t}(\sigma _\text{R}T_{}^3/\mathrm{\Gamma })t`$. Eq.(11) is an advection equation. Apparently a perturbation moves inwards, and its amplitude damps, since $`\stackrel{~}{\psi }`$ is conserved in amplitude along the inward moving characteristic, and therefore the comoving amplitude of $`\stackrel{~}{\psi }r^2`$ becomes smaller (DCHLC99).

### 3.2 Hydrodynamic perturbations at thermal equilibrium

When the cooling time of the disk is much shorter than the orbital time, one can assume that the disk is always in thermal equilibrium, and the disk height will react to changes in the irradiation flux in a hydrodynamic way, instead of being in vertical hydrostatic balance. In principle one should solve the time-dependent vertical hydrodynamics of the disk (1-D hydrodynamics). But for the present analysis we simplify the picture, and assume that the density profile is always given by Eq.(3), but with $`H_\text{p}`$ a function of time. This assumption can be justified, since in the early (linear) stages of a possible instability one has $`dv_z/dtz`$, which preserves the Gaussian shape of the density profile. By neglecting quadratic terms in the vertical velocity, one can then derive

$$\frac{d^2H_\text{p}}{dt^2}=\frac{kT_\text{e}}{\mu _\text{g}m_\text{u}}\frac{1}{H_\text{p}}\frac{GM_{}}{R^3}H_\text{p}.$$
(12)
When we define the dimensionless time $`\stackrel{~}{t}\sqrt{GM_{}/R_{}^3}t`$, a dimensionless form of the above equation is

$$\frac{^2h(r,\stackrel{~}{t})}{\stackrel{~}{t}^2}=\frac{C^{1/4}}{r^{5/2}h(r,\stackrel{~}{t})}\left(r\frac{h(r,\stackrel{~}{t})}{r}\right)^{1/4}\frac{h(r,\stackrel{~}{t})}{r^3},$$
(13)
with $`C`$ as defined in Eq.(7)

We now analyze the time-dependent behavior of linear perturbations of the power law solution (Eq.9). We take

$$h(r,\stackrel{~}{t})=\left(\frac{2}{7C}\right)^{1/7}r^{2/7}[1+\psi (r,\stackrel{~}{t})].$$
(14)
Substituting this in Eq.(13) yields (exactly):

$$\frac{^2\psi }{\stackrel{~}{t}^2}=\frac{1}{r^3}\left\{\frac{1}{(1+\psi )}\left(1+\psi +\frac{7}{2}r\frac{\psi }{r}\right)^{1/4}(1+\psi )\right\}.$$
(15)
For $`\psi 1`$ this becomes

$$\frac{^2\psi }{\stackrel{~}{t}^2}=\frac{7}{8r^3}\left\{2\psi r\frac{\psi }{r}\right\}.$$
(16)
This equation has solutions of the form

$$\psi =Ar^2e^{\sigma \stackrel{~}{t}}\mathrm{sin}\left(\omega \stackrel{~}{t}kr^3+\varphi _0\right),$$
(17)
with the following dispersion relations:

$`\omega ^2`$ $`=`$ $`\sigma ^2,`$ (18)
$`k`$ $`=`$ $`(16/21)\omega \sigma .`$ (19)
For inwards propagating waves (i.e. $`k/\omega <0`$) one has $`\sigma >0`$, and hence an exponentially growing mode111After submission the author became aware of work by E. Chiang (to be published) which reports on similar results.. The instability growth rate is inversely proportional to the wavelength of the perturbation, which indicates that the instability will be dominated by the shortest wavelengths. It should be noted that the validity of the equations breaks down at the scale of the disk vertical height. For modes with larger $`k`$ one should take radial radiative diffusion into account. This is, however, beyond the scope of the present set of equations.

## 4 Conditions for instability

The instability only sets in under certain physical conditions. The most obvious conditions for this instability to occur are:

1. Radiative cooling/heating time scale short compared to vertical sound crossing time, i.e. $`t_{\text{dyn}}t_{\text{therm}}`$.
2. Thermal dust-gas coupling is strong enough to keep the gas cooling time shorter than the vertical sound crossing time.
3. The equatorial temperature is dominated by irradiation, i.e. viscous dissipation by accretion is small in comparison.
4. The disk is optically thick to starlight in radial (grazing) direction.
5. The irradiation occurs through flaring; i.e. the reflected stellar flux from a scattering corona is weak in comparison to the direct interception of stellar light through flaring.

We have not investigated whether disk models that obey all the above points, but are vertically optically thin in the near/mid infrared, are also unstable.

Condition 1 is fulfilled for $`RR_{\text{inst}}`$, where $`R_{\text{inst}}`$ can be derived by equating the vertical sound crossing time with the thermal time: $`t_{\text{dyn}}=t_{\text{therm}}`$. We have

$`t_{\text{dyn}}`$ $``$ $`\sqrt{{\displaystyle \frac{R^3}{GM_{}}}},`$ (20)
$`t_{\text{therm}}`$ $``$ $`{\displaystyle \frac{k\mathrm{\Sigma }}{2m_\text{u}\sigma _\text{R}T_\text{e}^3}}.`$ (21)
The surface density $`\mathrm{\Sigma }(R)`$ (in units of $`\text{g}/\text{cm}^2`$) is a free function of the model, as long as the disk remains optically thick to the incident stellar radiation. Equating $`t_{\text{dyn}}=t_{\text{therm}}`$ yields the radius $`R_{\text{inst}}`$ beyond which the instability may set in:

$$\frac{R_{\text{inst}}}{\text{AU}}=3\times 10^5\left(\frac{M_{}}{M_{}}\right)^{13/3}\left(\frac{L_{}}{L_{}}\right)^4\left(\frac{H_\text{s}}{H_\text{p}}\right)^4\mathrm{\Sigma }^{14/3},$$
(22)
where $`L_{}`$ is the luminosity of the star.

The instability found here is an instability of the highly simplified analytic flaring disk models. Including more realistic physics could perhaps stabilize the disk. Some elements of realistic physics that should be considered in future work before one can claim that the instability is real are:

1. 2-D axisymmetric hydrodynamics.
2. Radial radiative transfer of stellar radiation instead of the inclination angle formula.
3. 1-D vertical radiative transfer, because the top layers have a higher temperature (see CG97).
4. Radial diffusion of radiation, or even fully 2-D radiative transfer, because radial exchange of energy might counter-act the shadowing effect which causes the instability.

## 5 Discussion

This Letter presents an analysis of the structure and stability of irradiated non-accreting protostellar disks based on a highly simplified set of equations. In the context of these equations it is concluded that the flaring non-accreting reprocessing disk model is unstable when the vertical sound crossing time is larger than the heating/cooling time scale. It should be noted, however, that detailed multidimensional numerical modeling is required to verify whether this instability is real or merely an artifact of the over-simplified equations.

The instability has a simple physical interpretation. Consider a flaring disk solution, and perturb it with an enlargement of the scale height at some point. The ‘sunny side’ of the hill will get overheated and expand vertically in a run-away fashion. The ‘shadow side’ will receive insufficient radiation and collapse. As the perturbations grow into the non-linear regime they start to cast a ‘real’ shadow ($`\alpha <0`$) over the disk at larger radii, thus completely depriving this part of the disk of irradiation. At this point the validity of the simplified equations breaks down, and what happens in this non-linear regime is unclear. Disk self irradiation (e.g. Bell Bell (1999)) might become important at this stage.

The reason why we find an instability, while DCHLC99 found only damping modes is because they considered only perturbations in the disk where the thermal time scale exceeds the vertical sound crossing time scale. While for actively accreting disks this condition is always satisfied, for passive disks the reverse may be true. This was recognized by DCHLC99, but an analysis of modes on a dynamic time scale in the outer regions of the protostellar disks was not carried out, hence their different conclusion.

###### Acknowledgements.

I wish to thank H. Spruit, C. Dominik, P. Armitage, S. Doty, J. Papaloizou, E. v. Dishoeck, G-J. v. Zadelhoff and the referee P. D’Alessio for useful suggestions and comments.


text:
# Young Stellar Groups, Runaway Stars, and Pulsars

## 1. Introduction

Most O and B stars in the Galactic disk are found in young star clusters or OB associations, where they often reside in binaries. About 30% of the O stars and 5–10% of the B stars are found in the field. These objects almost invariably have large space velocities (30–200 km s-1). Most of these ‘runaway’ stars appear to be single. Many have high He abundance, and rotate fast.

Two scenarios for the formation of runaway stars are considered viable: (i) a supernova explosion in a massive binary (Zwicky 1957; Blaauw 1961), in which the primary becomes a compact object, and the secondary moves away with a velocity comparable to its pre-explosion orbital velocity, and (ii) dynamical ejection from a dense stellar cluster (Poveda et al. 1967; Gies & Bolton 1986). Both scenarios can explain the properties of the ensemble of runaway stars observed to date. As a result, the relative importance of the two scenarios has been debated extensively in the past decades.

In order to find specific rather than statistical evidence for either scenario, we need to establish a common origin for either a runaway star and a pulsar or a runaway star and a dense stellar cluster. The availability of Hipparcos milli-arcsec (mas) global astrometry for the nearby stars, and of pulsar astrometry provided by VLBI and timing measurements, therefore stimulated us to carry out a systematic study of the past orbits of the nearby runaways and pulsars. The full study will be reported elsewhere (Hoogerwerf et al. 2000); here we restrict ourselves to a brief summary of some of the main results.

## 2. A binary supernova in Upper Scorpius

Blaauw (1952) proposed that the O9.5V runaway star $`\zeta `$ Oph originated in the nearby association Sco OB2. It is a single star, appears to be a blue straggler, and has He abundance $`Y0.40`$ and $`v_{\mathrm{rot}}\mathrm{sin}i=348`$ km s-1. The space motion of $`\zeta `$ Oph is $``$40 km s-1 relative to Sco OB2, and the star traversed the 5 Myr old subgroup Upper Scorpius $``$1 Myr ago. The kinematics of the expanding HI shell which surrounds Upper Scorpius is consistent with a supernova explosion around the same time, and the present-day mass function of the subgroup indicates that about 40 $`M_{}`$ is missing (de Geus 1992). It is therefore natural to speculate that the supernova occurred in a massive binary, and that $`\zeta `$ Oph is the former secondary. Mass transfer during the binary evolution prior to the explosion would then be the cause of the increased He abundance, the large rotation velocity, and the blue straggler nature. Can we find the compact object?

The catalog of pulsars with distances and proper motions maintained by Taylor, Manchester & Lyne (1993) contains seven objects within 1 kpc with measurement accuracy better than 10%. The radial velocities $`v_{\mathrm{rad}}`$ of these pulsars are unknown, so we retraced their paths in the gravitational field of the Galaxy by numerical orbit integration for a range of $`\pm `$ 500 km s-1 in $`v_{\mathrm{rad}}`$. Taking into account the characteristic age $`P/2\dot{P}`$ of the pulsars leaves only one likely candidate: PSR J1932+1059. If its $`v_{\mathrm{rad}}=200\pm 50`$ km s-1, then it passed through Upper Scorpius about 1 Myr ago, on its way towards the Galactic plane (which it recently crossed). If the pulsar did not originate in Upper Scorpius, then it must have been launched from the Galactic plane some $``$ 50 Myr ago. This is unlikely given a characteristic age of about 3 Myr.

Having identified the most likely pulsar, we computed the absolute minimum separation $`D_{\mathrm{min}}`$ between the orbits of $`\zeta `$ Oph and PSR J1932+1059, while sampling the error distributions in position and velocity. The main uncertainties are the errors in the parallax of $`\zeta `$ Oph ($`\pi =7.1\pm 0.7`$ mas) and the pulsar ($`\pi =5\pm 1.5`$ mas, Campbell 1995), and the remaining range in $`v_{\mathrm{rad}}`$. We found that the distribution of $`D_{\mathrm{min}}`$ is consistent with $`\zeta `$ Oph and the pulsar being in the same location $`1.0\pm 0.1`$ Myr ago in Upper Scorpius (Figures 1 & 2). This is compelling dynamical evidence for the binary supernova scenario.

The value of $`P/2\dot{P}`$ is an uncertain age indicator, and the $`3`$ Myr for PSR J1932+1059 is consistent with the dynamical age of 1 Myr derived here. While the current period is 0.22 seconds, the implied period at birth is 0.18 seconds.

Pulsars are expected to receive a kick velocity $`\stackrel{}{v}_{\mathrm{kick}}`$ at birth. Observations of the ensemble of pulsars suggest that $`\stackrel{}{v}_{\mathrm{kick}}`$ is a few hundred km s-1 (Hartman 1997; Hansen & Phinney 1997). Assuming that $`\zeta `$ Oph and PSR J1932+1059 originated in the same binary allows a determination of $`\stackrel{}{v}_{\mathrm{kick}}`$. Based on the angle between the orbits of the two objects, and the magnitude of the space velocities, we find $`|\stackrel{}{v}_{\mathrm{kick}}|=350\pm 50`$ km s-1. New VLBI observations of the proper motion and parallax of PSR J1932+1059 which are being obtained by Campbell (priv. comm.) will make it possible to improve this estimate further.

## 3. A stellar encounter in Orion

Blaauw & Morgan (1954) drew attention to the isolated stars AE Aur (O9.5 V) and $`\mu `$ Col (O9.5V/B0V), which run away in opposite directions from the association Ori OB1 with space velocities of about 100 km s-1 each (Figure 3). Blaauw & Morgan speculated that both stars originated in the same event in Ori OB1, some 2.5 Myr ago. Gies & Bolton (1986) suggested that this event was an encounter between two hard binaries, which also produced the eccentric massive binary $`\iota `$ Ori (O9III+B1III).

We used the globally accurate Hipparcos proper motions and parallaxes, and published radial velocities to retrace the orbits of $`\iota `$ Ori and the two runaways. We carried out $`2.5\times 10^6`$ experiments to sample the error distribution in the measurements (the largest uncertainties are in the parallax, as the objects are at distances of $`\mathrm{}>450`$ pc). We found that the distribution of the radius $`D_{\mathrm{min}}`$ of the minimum volume enclosing the three objects is consistent with AE Aur, $`\mu `$ Col and $`\iota `$ Ori being simultaneously in one location $`2.5\pm 0.2`$ Myr ago. This is strong evidence for the dynamical ejection scenario.

It is of interest to ask, where did the binary-binary encounter occur? We assumed that the center-of-mass motion of the four stars involved in the encounter was identical to the mean motion of the parent cluster. By using conservation of linear momentum, we then computed the current location and mean motion of the cluster. This is somewhat sensitive to the masses of AE Aur and $`\mu `$ Col, but all data are consistent with the binary-binary encounter taking place in the nascent Trapezium cluster (Figure 3). This identification is supported by the age of $``$2 Myr (Palla & Stahler 1999), the high stellar density, the short dynamical time as evidenced by mass segregration (Hillenbrand & Hartmann 1998), and the high binary fraction (Weigelt et al. 1999) of the Trapezium cluster.

## 4. Conclusions and prospects

The results from the previous two sections demonstrate that runaway stars are created both by binary supernova explosions and by dynamical ejection, and that it is now possible to identify individual events.

Enlarging the sample would provide a probe of the recent history of massive star formation in the Galactic disk, and would constrain the fraction of binaries in young clusters (Portegies Zwart 1999). We therefore took all 3622 O–B5 stars in the Hipparcos Catalog with (i) measured radial velocity $`v_{\mathrm{rad}}`$ (1118), (ii) parallax $`\pi >2\sigma _\pi `$ (which translates to a distance less than about 700 pc), (iii) accurate proper motion $`\sigma _\mu /\mu <0.1`$, and (iv) space velocity $`\stackrel{}{v}>30`$ km s-1. This produced 56 runaway candidates. Similarly, we took all the pulsars in the Taylor et al. (1993) catalog with (i) measured distance and proper motions (94), a distance estimate of $`D<1`$ kpc, and (iii) accurate proper motion $`\sigma _\mu /\mu <0.1`$. This leaves only seven pulsars. We then retraced the orbits of all these objects, to identify the sites of origin. This work is in progress.

Both samples are severely incomplete. The Hipparcos Catalog is complete only for $`V<7.39`$ (depending on latitude), and less than a third of the O and B stars in it have a measured radial velocity. Because of beamed radio emission, we cannot observe all pulsars, and not all of those that do radiate in our direction have been found; of these, only a few have an accurately measured proper motion. Even so, to date we have found the parent stellar group for two more pulsars (in addition to PSR J1932+1059) and for $``$20 runaway stars (10 of which are new). The two pulsars both originate in the $`50`$ Myr old association Per OB3. This is not surprising, as the peak of pulsar production occurs for the lower-mass B stars, i.e., in old groups. Twenty of the 56 runaway stars can be linked to the nearby associations. Other runaways in the sample presumably originate in open clusters, including a newly-discovered pair of stars running away in opposite directions from the region of the $`\lambda `$ Ori cluster. Further progress requires accurate space motions for the young open clusters in the Solar neighbourhood. Some of the remaining runaways (and pulsars) must have originated beyond 700 pc, where our knowledge of the parent groups is poor.

Much can be learned by measuring radial velocities for the nearby O and B field stars, and by increasing the sample of pulsars with measured astrometry. At the current accuracy of about 1 mas in parallax and 1 mas/yr in proper motion, a full census is possible of the nearby 500 pc. Extension to a larger volume requires $`\mu `$as astrometry, which will be provided by future space observatories such as DIVA, FAME, and GAIA.

#### Acknowledgments.

It is a pleasure to thank A. Blaauw, R. Campbell and M.A.C. Perryman for stimulating conversations.

## References

Blaauw, A. 1952, BAN, 11, 414

Blaauw, A., & Morgan, W. W. 1954, ApJ, 119, 625

Campbell, R. M. 1995, PhD Thesis, Harvard University

de Geus, E. J. 1992, A&A, 262, 258

Gies, D. R., & Bolton, C. T. 1986, ApJS, 61, 419

Hansen, B. M. S., & Phinney, E. S. 1997, MNRAS, 291, 569

Hartman, J. W. 1997. A&A, 322, 127

Hillenbrand, L. A., & Hartmann, L. W. 1998, ApJ, 492, 540

Hoogerwerf, R., de Bruijne, J. H. J., & de Zeeuw, P. T. 2000, A&A, submitted

Lyne, A. G., Anderson, B., & Salter, M. J. 1982, MNRAS, 201, 503

Palla, F., & Stahler, S. W. 1999, ApJ, 525, 772

Portegies Zwart, S. 2000, preprint \[astro-ph/0005245\]

Poveda, A., Ruiz, J., & Allen, C. 1967, Bol Obs Tonantzintla, 28, 86

Taylor, J. H., Manchester, R.N., & Lyne, A.G. 1993, ApJS, 88, 529

Weigelt, G., Balega, Y., Preibisch, T., Schertl, D., Schöller, M., & Zinnecker, H. 1999, A&A, 347, L15

de Zeeuw, P. T., Hoogerwerf, R., de Bruijne, J. H. J., Brown, A. G. A., & Blaauw, A. 1999, AJ, 117, 354

Zwicky, F. 1957, Morphological Astronomy, Springer–Verlag, Berlin

### Pavel Kroupa:

1) Dynamical ejection from the Trapezium cluster $``$2.5 Myr gives an important alternative age estimate that is larger than the typically quoted age. 2) My N-body calculations show that up to 50% of all stars with $`M>8M_{}`$ are found at $`R>2R_{\mathrm{tidal}}`$ after a few Myr. 3) Ejected low-mass binaries will have highly eccentric orbits. This allows for a test of tidal circularisation theory of fully convective stars.


text:
# The Influence of Central Black Holes on Gravitational Lenses

## 1 INTRODUCTION

Recent observations suggest that many if not all nearby galaxies host massive central black holes (e.g. Kormendy & Richstone 1995; Magorrian et al. 1998). The existence of such central black holes can be accommodated by hierarchical structure formation theories (e.g. Silk & Rees 1998; Kauffmann & Haehnelt 2000). The effects of central black holes on gravitational lensing have not been studied in detail, although it is commonly believed that the central singularity can suppress the central image. This is important because lensing theories predict that any non-singular lens should have an odd number of images while the observed lenses always show an even number of images (Burke 1981; see also Schneider et al. 1992 for a general review on gravitational lensing); the only possible exception is APM 08279+5255 (Ibata et al. 2000). The lack of central images could indicate either that the potential is singular or that the central surface density in lenses is so high that the central image is highly demagnified and thus unobservable in current surveys (Narayan et al. 1984). The purpose of this paper is to clarify, through simple examples, the effects of central black holes on lensing properties. We show that the presence of black holes can not only suppress but also create additional images. These additional images may be observable in some cases. Central black holes also introduce somewhat unusual critical curves and caustic structures. The outline of this paper is as follows. In section 2, we first outline the lensing basics and derive some analytical results for the critical curves, caustics and image properties for an isothermal sphere plus a black hole. We then generalize to the case of an isothermal ellipsoid plus a black hole. In section 3, we discuss the implications of black holes on gravitational lenses, including the constraints on the core radius.

## 2 EFFECTS OF CENTRAL BLACK HOLES ON LENSING

We model a lensing galaxy as an isothermal ellipsoid plus a central black hole. The isothermal ellipsoid model is not only analytically tractable but also consistent with models of individual lenses, lens statistics, stellar dynamics, and X-ray galaxies (e.g. Fabbiano 1989; Maoz & Rix 1993; Kochanek 1995, 1996b; Grogin & Narayan 1996; Rix et al. 1997). The surface density distribution of an isothermal ellipsoid is given by

$$\frac{\mathrm{\Sigma }}{\mathrm{\Sigma }_{\mathrm{cr}}}=\frac{1}{2q}\frac{1}{\sqrt{x^2+y^2/q^2+r_\mathrm{c}^2}},$$
(1)
where $`r_\mathrm{c}`$ is the core radius, $`q`$ is the axial ratio, and $`\mathrm{\Sigma }_{\mathrm{cr}}=c^2D_\mathrm{s}/(4\pi GD_\mathrm{l}D_{\mathrm{ls}})`$ is the critical surface density, $`D_\mathrm{l}`$, $`D_\mathrm{s}`$ are angular diameter distances from the observer to the lens and source, respectively, and $`D_{\mathrm{ls}}`$ is the angular diameter distance from the lens to the source. Notice that all the lengths ($`x,y,r_\mathrm{c}`$) are expressed in units of the critical radius, $`R_{\mathrm{cr}}`$,

$$R_{\mathrm{cr}}=D_\mathrm{l}b,b4\pi \left(\frac{\sigma }{c}\right)^2\frac{D_{\mathrm{ls}}}{D_\mathrm{s}},$$
(2)
where the critical angle $`b`$ is the angle extended by the critical radius on the sky ($`b0.23`$ arcsecond for typical lens galaxies; e.g. Jackson et al. 1998), and the velocity dispersion $`\sigma `$ in eq. (2) is related to, but not identical to, the observable line of sight velocity dispersion; we shall ignore this minor complication in our analysis and simply treat it as a parameter (see Keeton, Kochanek & Seljak 1997 for further discussions).

The lensing properties of the isothermal ellipsoid have been given by Kassiola & Kovner (1993), Kormann, Schneider & Bartelmann (1994), and Keeton & Kochanek (1998). The lens equation of an isothermal ellipsoid plus a black hole is

$`u`$ $`=`$ $`x{\displaystyle \frac{1}{\sqrt{1q^2}}}\mathrm{tan}^1\left({\displaystyle \frac{\sqrt{1q^2}x}{\phi +r_\mathrm{c}}}\right)m{\displaystyle \frac{x}{r^2}},`$
$`v`$ $`=`$ $`y{\displaystyle \frac{1}{\sqrt{1q^2}}}\mathrm{tanh}^1\left({\displaystyle \frac{\sqrt{1q^2}y}{\phi +q^2r_\mathrm{c}}}\right)m{\displaystyle \frac{y}{r^2}},`$ (3)
where $`\phi ^2=q^2(r_\mathrm{c}^2+x^2)+y^2,r^2=x^2+y^2`$, and the dimensionless black hole mass is given by

$$m=\frac{M_{\mathrm{bh}}}{M_0},M_0\frac{\pi \sigma ^2R_{\mathrm{cr}}}{G},$$
(4)
with $`M_{\mathrm{bh}}`$ being the mass of the black hole. Physically, $`M_0`$ is the mass of the galaxy contained within a cylinder with radius $`R_{\mathrm{cr}}`$ and hence $`m_0`$ is essentially the ratio of the black hole mass to the mass of the galaxy in the inner parts (within $`R_{\mathrm{cr}}`$). The magnification ($`\mu `$) of a given image is given by

$$\mu ^1=\frac{u}{x}\frac{v}{y}\frac{v}{x}\frac{u}{y}.$$
(5)
For definitiveness, we shall adopt a lens redshift of 0.5 and a source redshift of 2. In the Einstein-de Sitter universe ($`\mathrm{\Omega }_0=1`$ with no cosmological constant), the angular and length units are given by

$$b=1^{^{\prime \prime }}\left(\frac{\sigma }{250}\right)^2,R_{\mathrm{cr}}=3.6h^1\mathrm{kpc}\left(\frac{\sigma }{250}\right)^2,$$
(6)
where $`h`$ is the Hubble constant in units of $`100\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$ and $`\sigma `$ is in units of $`\mathrm{km}\mathrm{s}^1`$. The dimensionless black hole mass for the given lens and source redshifts is

$$m=2.4\times 10^3h\left(\frac{\sigma }{250}\right)^{3.27},$$
(7)
where we have used the tight correlation of the black hole mass with velocity dispersion found by Ferrarese & Merritt (2000), $`M_{\mathrm{bh}}4\times 10^8M_{}(\sigma /250)^{5.27}`$ (see Gebhardt et al. 2000 for a different scaling.)111We neglect any possible evolution of the black hole mass with time.

The lens equation (eq. 3) can be readily solved numerically. In general there are several curves in the image plane along which the magnification is infinite ($`\mu ^1=0`$). These are called “critical curves”, and they map to “caustics” in the source plane. Caustics mark discontinuities in the number of images, so in order to determine the number of images produced by a lens model, it is sufficient to examine the caustics. We shall consistently employ this idea in this paper. In order to gain some analytical insights into the influence of central black holes, we start with the spherical case (i.e. $`q=1`$) in the next subsection.

### 2.1 Spherical Isothermal Models with Black Holes

For the spherical case ($`q=1`$), the surface density is given by

$$\frac{\mathrm{\Sigma }}{\mathrm{\Sigma }_{\mathrm{cr}}}=\frac{1}{2\sqrt{r^2+r_\mathrm{c}^2}}.$$
(8)
This profile has been studied quite extensively by Hinshaw & Krauss (1987).

The lens equation for this profile in the presence of a black hole is quite simple,

$$r_\mathrm{s}=r\alpha (r),\alpha (r)=\frac{\sqrt{r^2+r_\mathrm{c}^2}r_\mathrm{c}}{r}+\frac{m}{r},$$
(9)
where $`r_\mathrm{s},r_\mathrm{c}`$ are the source position and core radius in units of the critical radius and $`m`$ the dimensionless mass of the central black hole given in eq. (4).

#### 2.1.1 Critical Curves and Caustics

The critical curves are given by the equation

$$\mu ^1=\frac{r_\mathrm{s}}{r}\frac{dr_\mathrm{s}}{dr}=0.$$
(10)
The solution satisfies either $`r_\mathrm{s}/r=0`$ or $`dr_\mathrm{s}/dr=0`$. One can verify that the first condition ($`r_\mathrm{s}/r=0`$) always yields one single critical curve,

$$r^2=\frac{1}{2}\left[\mathrm{\hspace{0.17em}1}2r_\mathrm{c}+2m+\sqrt{(12r_\mathrm{c})^2+4m}\right],$$
(11)
which maps to a degenerate caustic point at the origin. The second condition ($`dr_\mathrm{s}/dr=0`$) can be manipulated into a polynomial

$`r^6+(2m2r_\mathrm{c}+r_\mathrm{c}^2)r^4`$ $`+`$
$`(m^22mr_\mathrm{c}+r_\mathrm{c}^2+2mr_\mathrm{c}^22r_\mathrm{c}^3)r^2`$ $`+`$
$`m^2r_\mathrm{c}^22mr_\mathrm{c}^3`$ $`=`$ $`0.`$ (12)
This equation is cubic in $`r^2`$, and so can be solved analytically. One finds that beyond a critical black hole mass, there is no physical solution (i.e. positive solution in $`r^2`$), while below the critical mass, there are always two physical solutions. The critical mass can be found by examining eq. (2.1.1); some algebra yields the critical mass

$$m_{\mathrm{cr}}=r_\mathrm{c}(1+r_\mathrm{c})3\left(\frac{r_\mathrm{c}^2}{2}\right)^{2/3}.$$
(13)
For a given core radius, the central images are always suppressed by the black hole if its mass exceeds the critical value, $`m_{\mathrm{cr}}r_\mathrm{c}`$ when $`r_\mathrm{c}1`$. This can be understood intuitively as follows: when the mass is comparable to $`r_\mathrm{c}`$, the mass of the black hole is comparable to the galaxy mass enclosed by the core radius, so the black hole makes the deflection angle nearly flat at $`rr_\mathrm{c}`$, similar to the singular isothermal case, for which we have only two images (cf. Fig. 2). Notice, however, that even when the black hole mass is below the critical value, there is a region in the source plane where the central image is suppressed (see below).

We can obtain approximate but more illuminating solutions for the positions of critical curves and caustics when $`r_\mathrm{c}1`$ and $`mm_{\mathrm{cr}}r_\mathrm{c}`$. The first assumption ($`r_\mathrm{c}1`$) is supported by the lack of central images in lenses (e.g. Wallington & Narayan 1993; Kochanek 1996a; see also §2.1.2), while the second condition is valid for most galaxies as can be seen from eq. (7), and if $`r_\mathrm{c}0.002`$ ($`10`$ pc for $`\sigma =250\mathrm{km}\mathrm{s}^1`$). Under these two assumptions, the inner most critical curve satisfies $`rr_\mathrm{c}`$ (justified below), the lens equation (eq. 9) simplifies to

$$r_\mathrm{s}=(1\kappa _\mathrm{c})r\frac{m}{r},$$
(14)
where $`\kappa _\mathrm{c}1/(2r_\mathrm{c})`$ is the central surface density of the galaxy, in the absence of the black hole. From this simplified lens equation, one can readily show that the innermost critical curve forms where $`dr_\mathrm{s}/dr=0`$, i.e.

$$r_{\mathrm{crit},1}=\sqrt{\frac{m}{\kappa _\mathrm{c}1}}.$$
(15)
The second assumption (i.e. $`mm_{\mathrm{cr}}`$) ensures that our condition, $`rr_\mathrm{c}`$, is satisfied. The corresponding caustic in the source plane is

$$r_{\mathrm{caust},1}=2\sqrt{m(\kappa _\mathrm{c}1)}=2mr_{\mathrm{crit},1}^1.$$
(16)
The middle critical curve in general does not satisfy the condition $`rr_\mathrm{c}`$. To derive its approximate position, we simplify eq. (2.1.1) by setting $`m=0`$, after which we are left with a quadratic equation in $`r^2`$. One finds that, to the leading order of $`r_\mathrm{c}`$, the positions of the critical curves and caustics are

$$r_{\mathrm{crit},2}\sqrt{r_\mathrm{c}},r_{\mathrm{caust},2}12\sqrt{r_\mathrm{c}}=12r_{\mathrm{crit},2}.$$
(17)
From eq. (11) one finds that the outer most critical curve and its corresponding caustic are at

$$r_{\mathrm{crit},3}1r_\mathrm{c},r_{\mathrm{caust},3}=0.$$
(18)
For a source located outside $`r_{\mathrm{caust},2}`$ there are only two images. When a source is between $`r_{\mathrm{caust},1}`$ and $`r_{\mathrm{caust},2}`$, there are four images; compared with the case without a black hole, one additional image has been created very close to the black hole. For a source inside $`r_{\mathrm{caust},1}`$, there are only two images, the central image that would be present without a black hole has been suppressed. The radius $`r_{\mathrm{caust},1}`$ can therefore be regarded as a measure of the zone of influence of the black hole: for any source $`r_\mathrm{s}<r_{\mathrm{caust},1}(2m/r_\mathrm{c})^{1/2}`$, the central image is destroyed. Notice that this zone of influence scales as $`m^{1/2}`$, identical to the scaling of the Einstein radius of an isolated black hole with mass.

Fig. 1 illustrates the evolution of critical curves and caustics with the mass of the black hole increasing from zero to 0.03. The core radius is taken to be 0.05, a somewhat arbitrary value which satisfies the upper limit inferred from the lack of central images (e.g. Wallington & Narayan 1993; Kochanek 1996a). When there is no black hole at the center (Fig. 1a) the critical curves are two circles, the outer critical curve maps to a degenerate caustic point at the origin, while the inner critical curve maps to the solid circle. Fig. 1b shows the configuration when we include a black hole mass with $`m=0.002`$, a value expected for a galaxy with $`\sigma 250\mathrm{km}\mathrm{s}^1`$ (cf. eq. 7). In addition to the two critical curves found when $`m=0`$, an additional tiny critical curve appears around the origin (too small to be seen), which maps to the inner solid circle. The positions of the three critical curves and their corresponding caustics are well approximated by eqs. (15-18). For a source moving from infinity to the origin, it would first have two images, then four images as the source moves across the first caustic. Compared with the case without a black hole, one additional image has been created very close to the black hole. The number of images decreases to 2 again as the source moves across the second caustic. Notice that even in this case where the black hole mass is substantially below the critical mass $`m=m_{\mathrm{cr}}0.0175`$, the zone of influence determined by $`r_{\mathrm{caust},1}0.28`$ is substantial. Within the zone of influence, the central image is completely suppressed. When the black hole mass is further increased to $`m=0.01`$ (Fig. 1c), the inner-most critical curve becomes larger while the middle critical curve shrinks, their corresponding caustics also approach each other. When the black hole mass is increased to the critical mass $`m_{\mathrm{cr}}`$ (eq. 13), the two critical curves merge into a single curve, so do their corresponding caustics. When the black hole mass is increased above the critical value (an example is shown for $`m=0.03`$ in Fig. 1d), one is left with a single critical curve and a degenerate caustic point at the center. In this case, there are only two images no matter where the source is. Notice, however, that the black hole mass for the panels 1c and 1d are somewhat unrealistic (cf. eq. 7.)

The change of image numbers may first appear puzzling, but can be easily understood in the conventional diagram of the deflection angle, $`\alpha (r)`$, vs. radius, shown in Fig. 2. Three examples are shown for $`r_\mathrm{c}=0.05`$ corresponding to three black hole masses, $`m=0,0.01`$, and 0.03. For each source position, the image positions are given by the the intercepts of the deflection angle curve with the straight line $`f(r)=rr_\mathrm{s}`$. For example, it is easy to see why for $`m=0.01`$, there are three critical curves. It is clear from this figure that the central black hole does not perturb the two outer images very much, as expected, it only strongly perturbs the central images.

#### 2.1.2 Image Positions and Magnifications

In this subsection, we will obtain approximate image positions and magnifications when there are four images, i.e., when the source position is between $`r_{\mathrm{caust},1}`$ and $`r_{\mathrm{caust},2}`$. We are particularly interested in the case when the source is close to $`r_{\mathrm{caust},1}`$. For such cases, there are usually two bright outer images and two faint central images; one of the fainter central image is created by the black hole. We would like to address the question whether such images are observable.

Using the resultant technique as discussed in Witt & Mao (1995), one can verify that the four image magnifications must satisfy an exact relation

$$\mu _1+\mu _2+\mu _3+\mu _4=2.$$
(19)
Note that all the magnifications include their parities (signs).

The brightest image has positive parity and is located in regions where $`rr_\mathrm{c}`$. For this image, the lens equation can be simplified into $`r^2r_\mathrm{s}r+r_\mathrm{c}r=0`$, where we have only retained terms to the order of $`(r_\mathrm{c}/r)^2`$. From this equation, we can obtain the position and magnification of the brightest image

$$r_11+r_\mathrm{s}r_\mathrm{c},\text{and}\mu _1^1\frac{r_\mathrm{s}}{r_1^3}(r_1^2+mr_\mathrm{c})\frac{r_\mathrm{s}}{r_1}.$$
(20)
Similarly, we can obtain the position for the bright outer image that has negative parity

$$r_21+r_\mathrm{s}+r_\mathrm{c},\text{and}\mu _2^1\frac{r_\mathrm{s}}{r_2}.$$
(21)
For singular isothermal spheres ($`r_\mathrm{c}=0`$), the magnifications and positions are exact for these two images, and they satisfy $`\mu _1+\mu _2=2`$.

For the two central images, we can estimate their positions and magnifications if they satisfy $`rr_\mathrm{c}`$. In this case, the lens equation can be approximated by eq. (14) and the image magnification is given by

$$\mu =\frac{1}{(1\kappa _\mathrm{c})^2m^2/r^4}.$$
(22)
For very small $`m`$ we can obtain for the positions and magnifications for the third and fourth image

$$r_3\frac{r_\mathrm{s}}{\kappa _\mathrm{c}1}\text{and}\mu _3\mu _\mathrm{c},\mu _\mathrm{c}\frac{1}{(\kappa _\mathrm{c}1)^2}.$$
(23)
$$r_4\frac{m}{r_\mathrm{s}}\text{and}\mu _4^1\frac{r_\mathrm{s}}{r_4^3}(r_\mathrm{s}r_4+2m)=\frac{r_\mathrm{s}^4}{m^2}.$$
(24)
Note that the magnification for the third image $`\mu _3\kappa _\mathrm{c}^2=4r_\mathrm{c}^2`$, is independent of the source position when the black hole mass is small. Previous studies neglected the role of central black holes and therefore the lack of central images immediately implies an upper limit on the core radius. If the brightest image is magnified by a factor of few, and no central image is seen with a dynamic range (DR) of 100, then this would imply $`r_\mathrm{c}0.1`$, which is roughly the published upper limits on the core radius from lenses (e.g. Wallington & Narayan 1993; Kochanek 1996a). Hence the lack of central images in lenses may mean either that the core radius is small (see §2.1.2) or that the black hole mass is quite massive; we return to this point in the discussion.

The fourth central image is usually the closest to the black hole and is usually very faint. However, when the source lies very close to $`r_{\mathrm{caust},1}`$, the second term in the denominator of eq. (22), which is the square of the shear ($`\gamma _{\mathrm{BH}}`$) from the black hole, can become comparable to $`(1\kappa _\mathrm{c})`$, and then the magnification is of order unity or higher. The question is: how large is this region? To estimate this, we combine eqs. (15) and (22) and find that the position of the image is related to its magnification by

$$\left(\frac{r}{r_{\mathrm{crit},1}}\right)^4=\left(1\frac{\mu _\mathrm{c}}{\mu }\right)^1,$$
(25)
where $`\mu _\mathrm{c}`$ is the magnification of an image formed well within the core radius in the absence of the black hole and is given in eq. (23). A small displacement $`\mathrm{\Delta }r`$ from the critical curve is related to the magnification through

$$\mathrm{\Delta }r=\frac{\mu _\mathrm{c}}{4\mu }r_{\mathrm{crit},1}.$$
(26)
This displacement is related to a corresponding displacement of the source in the source plane, through eq. (14). By definition the first order term of a Taylor expansion of the lens equation vanishes near the caustic. The second term, after some algebra, gives

$$\mathrm{\Delta }r_\mathrm{s}=\frac{m}{r_{\mathrm{crit},1}^3}\times \mathrm{\Delta }r^2=\frac{\sqrt{m}\mu _\mathrm{c}^{7/4}}{16\mu ^2}.$$
(27)
The probability of the source in a multiply-imaged system, lying inside this region (i.e. having a magnification $`>\mu `$), is therefore

$$P_{\mathrm{BH}}(>\mu )=\frac{2r_{\mathrm{caust},1}\mathrm{\Delta }r_\mathrm{s}}{r_{\mathrm{caust},2}^2}=\frac{m\mu _\mathrm{c}^{3/2}}{4\mu ^2(12\sqrt{r_\mathrm{c}})^2}.$$
(28)
For any realistic situation, this probability is exceedingly small ($`10^6`$); the inclusion of magnification bias and extended source structure does not substantially enhance the probability of this faint central image to have magnifications above unity.

In case of high dynamic-range (DR) maps, however, central images with $`\mu 1`$ can still be detected. In such cases eq. (28) is no longer valid and we have to find a new estimate for the probability. To do this, we insert eq. (25) into the simplified lens equation, eq. (14), and arrive at

$$r_\mathrm{s}=\frac{1}{2}r_{\mathrm{caust},1}\left[\left(1\frac{\mu _\mathrm{c}}{\mu }\right)^{1/4}+\left(1\frac{\mu _\mathrm{c}}{\mu }\right)^{+1/4}\right].$$
(29)
For $`\mu \mathrm{}`$, the source position goes to $`r_{\mathrm{caust},1}`$, as expected. In case $`\mu =\mu _4=10^3`$ (DR$`10^3`$; the magnification is negative because the parity of the central image created by the BH is $``$1), the source position $`r_\mathrm{s}1.15r_{\mathrm{caust},1}`$, assuming $`r_\mathrm{c}=0.05`$ and $`m=0.002`$. The probability of observing this image is then $`2\pi r_{\mathrm{caust},1}(r_\mathrm{s}r_{\mathrm{caust},1})/(\pi r_{\mathrm{caust},2}^2)`$, which is $``$8% for the above given values of $`r_\mathrm{c}`$ and $`m`$. If $`r_\mathrm{c}0.05`$, the currently known sample of 20–30 radio graviational lens systems could therefore contain several systems with central images that are observable, if the dynamic range in radio observations can be improved by a factor of 5–10. To separate the two central images, which are only a few milli-arcseconds apart, one obviously requires high-resolution VLBI observations.

### 2.2 Elliptical Isothermal Density Distribution With Black Holes

The spherical case we studied in the previous subsection is idealized, since galaxies nearly always show some ellipticity. In this section, we will study the critical curves and caustics for a more realistic elliptical density distribution plus a black hole. The critical curves can be found by solving eq. (5) numerically and the caustics are found using the lens mapping (eq. 3). Fig. 3 shows the critical curves and caustics structure as a function of black hole mass for an elliptical density mass distribution with $`q=0.7`$ and again $`r_\mathrm{c}=0.05`$. When there is no black hole at the center (Fig. 3a), one sees the well known critical curves and caustics. When a source is inside the diamond caustic, there are five images. When a source is outside the diamond but inside the elliptical caustic, there are three images. When the source moves further outside, there is only one single image. The image configurations for three source positions are illustrated in the same panel. Usually the central images are quite faint and difficult to observe. Now as we include a black hole mass with $`m=0.002`$ (cf. eq. 7), an additional small critical curve appears very close to center (too small to be seen), and this maps to the middle elliptical caustic. A source outside the diamond but inside this middle elliptical caustic has two images. This is analogous to Fig. 1b. When a source moves from infinity to the origin, the image number changes according to $`2424`$, with each caustic crossing either increases or decreases the number of images by two. When the black hole mass is further increased to $`m=0.01`$ (Fig. 3c) the inner-most critical curve becomes larger while the middle critical curve shrinks, their two corresponding caustics also approach each other. When the black hole mass is further increased to or above a critical value, $`0.02`$ (a value slightly larger than that in the spherical case), the two inner critical curves and their corresponding caustics cancel each other, and one is left with an elliptical critical curve and a diamond caustic at the center (an example is shown for $`m=0.03`$); this caustic separates the outer two image regions from the inner four image regions. In either region, the central images have been suppressed by the black hole. The behavior in the elliptical density distribution case is qualitatively similar to the spherical case, the main difference is that the ellipticity breaks the spherical symmetry and changes the central degenerate caustic point in the spherical case into the central diamond (compare Figs. 1 and 3.); the caustic change due to ellipticity is similar to the case without central black holes (see e.g. Schneider et al. 1992).

## 3 Discussion

As clearly demonstrated in Figs. 1–3, central black holes in gravitational lenses introduces qualitatively new features in the critical curves and caustics. For realistic black hole masses, the black hole introduces a new region in the source plane that has four images. However, this additional four-image region does not help to resolve the problem of the apparent excess of quadruple lenses relative to double lenses (King & Browne 1996; Kochanek 1996b). The reason is that these four image configurations have two very faint images at the center and two bright images outside, which are very different from the observed quadruple lenses where the images are roughly on a circle from the lens center. The image closest to the black hole is the faintest image; for the case shown in Fig. 3b, this image is about a factor of $`10^5`$ times fainter than the brightest one. As shown in section 2.1.1, there is only a negligible probability that the central images can reach magnifications of order unity. However, we find that there is a non-negligible probability ($``$8% for $`r_\mathrm{c}=0.05`$ and $`m=0.002`$) that the central image can reach magnification $`10^3`$, with the two distinct central images separated by about $`10^3b`$ ($``$ milli-arcseconds). Such images can in principle be detected with future high dynamic range and high resolution VLBI observations. Note that this probability obviously depends on the values of $`m`$ and $`r_\mathrm{c}`$, and is uncertain since for the latter only upper limits are available at present.

The massive black holes also complicate the interpretation of core radius inferred due to the lack of central images. The central image will be completely suppressed for all source positions if the black hole mass exceeds the critical value as given in eq. (13). This critical mass is reached when the ratio of the mass of the black hole to the galaxy mass in the inner parts (more precisely, $`M_0`$ in eq. 4) is approximately equal to the core radius in units of the critical angle $`b`$. For a galaxy with velocity dispersion of $`250\mathrm{km}\mathrm{s}^1`$, the mass of the central black hole is $`4\times 10^8M_{}`$, and $`M_01.6\times 10^{11}h^1M_{}`$. The black hole is super-critical when $`r_\mathrm{c}0.002`$; in physical units, this implies an upper limit on the core radius of about $`10`$ pc for typical lens configuration with $`R_{\mathrm{cr}}\mathrm{few}\mathrm{kpc}`$. This upper limit is usually not reached except for lens systems with favorable image geometry such as B1030+074 (Norbury et al. 2000, in preparation.) However, the presence of central black holes can also suppress the central image even when it is much below the critical value in some regions (see Figs. 1b and 3b). The zone of influence of the black hole is characterized by the $`r_{\mathrm{caust},1}`$ (eq. 16) in the spherical case: when the source is inside the zone ($`r_\mathrm{s}<r_{\mathrm{caust},1}(2m/r_\mathrm{c})^{1/2}`$), the central image is destroyed. For a given source position $`r_\mathrm{s}`$, the central image disappears when the black hole mass satisfies

$$m\frac{r_\mathrm{c}r_\mathrm{s}^2}{2}\frac{r_\mathrm{c}}{2(\mu _11)^2},$$
(30)
where in the last step we have replaced the source position on the right hand side by the magnification of the brightest image $`\mu _1`$ using eq. (20). For a galaxy with $`\sigma =250\mathrm{km}\mathrm{s}^1`$, $`m0.002`$, and $`\mu _14`$, the central image will be destroyed when $`r_\mathrm{c}<0.04`$ ($`100\mathrm{pc}`$ for a typical critical radius of a few kpc). More generally, the influence of black holes on the core radius constraint has to be obtained from detailed modelling. Eq. (30) also highlights that the best lenses to search for central images are systems with low magnifications. In such systems, the source is closer to outer caustic ($`r_{\mathrm{caust},2}`$) and hence more likely to be outside the zone of influence. As a result, the central images are less likely to be suppressed by the black hole. Unfortunately, these systems also have large flux ratios between the two bright out images (cf. eqs. 20-21). Since lens surveys favor systems with large magnifications due to magnification bias and systems with small flux ratio, it is not surprising that very few lenses discovered so far show central images, perhaps as a result of both small core radii and central black holes in gravitational lenses.

## Acknowledgments

We thank Ian Browne for helpful discussions and Chuck Keeton for insightful criticisms on a draft of the paper.


text:
# Non-LTE effects in neutron star atmospheres

## Abstract

We present the first non-LTE model atmospheres for neutron stars (NS). We study their structure and NLTE effects on the emergent thermal radiation of old isolated NS.

Institut für Astronomie und Astrophysik, Univ. Tübingen, Germany

Our models assume plane-parallel geometry, hydrostatic and radiative equilibrium. The radiation transfer equations are solved simultaneously with the NLTE rate equations by using an Accelerated Lambda Iteration technique (Werner & Dreizler 1999). Blanketing by millions of iron lines is considered with an Opacity Sampling method. The models neglect magnetic fields, which is a good approximation as long as B$`<10^8`$$`10^{10}`$G, depending on $`T_{\mathrm{eff}}`$ and chemical composition (Zavlin et al. 1996). Iron opacities are calculated from Kurucz line lists, and bound-free cross-sections from Opacity Project data.

We have computed a set of three models with different chemical composition: pure H, He, and Fe models. We chose $`T_{\mathrm{eff}}`$=200 000 K, $`\mathrm{log}g`$ =14.39, corresponding to a NS with M=1.4 M and R=10 km. Fig. 1 shows their emergent fluxes. They display distinctly harder high energy tails than a blackbody. Fig. 2 shows details of the model structures. Due to NLTE effects, temperatures show a slight increase in the surface layers. Departures from LTE are most pronounced in the iron model. They affect the line formation regions and hence fluxes in the strongest Fe line cores, typically of the order 10%. The general spectral shape is unaffected. NLTE effects might be even more important at higher $`T_{\mathrm{eff}}`$.

Werner K., Dreizler S. 1999, Journal of Computational and Applied Math., in press
Zavlin V.E., Pavlov G.G., Shibanov Y.A. 1996, A&A 315, 141


text:
# 𝑛(𝑘⃗) Structure and Matrix Elements in Photoemission: Results from Insulating Ca2CuO2Cl2

\[

## Abstract

The photon energy dependence of angle-resolved photoemission data on Ca2CuO2Cl2 is presented, and the role of the matrix element is investigated. It is shown that despite strong variations in the matrix element, the loss in intensity as one crosses the antiferromagnetic zone boundary is a robust feature of the Mott insulator. As expected for a two dimensional system, the dispersion remains independent of the photon energy even when the matrix element does not.

\]

The question of how a metal evolves into an insulator is one of the most fundamental in solid state physics. For non-interacting electrons, the Fermi surface shrinks and eventually disappears as a band is filled. When electron-electron interactions dominate, the situation is less clear. By considering the strong Coulomb interaction, Mott qualitatively described how a material predicted by band theory to be a metal would in fact be an insulator. However, it remains unclear as to how the details of the electronic structure evolves from a half-filled metal to a Mott insulator. One of the most characteristic features of any metal is its Fermi surface which can be defined by the surface of steepest descent in the electron momentum distribution function, $`n(\stackrel{}{k})`$. Thus an equivalent question to the one above is how does the Fermi surface of a metal vanish as strong electron correlations drive the system into an insulator?

In the context of specific many-body models such as the Hubbard model, it has been shown that a structure in $`n(\stackrel{}{k})`$ survives even when the on site Coulomb U drives the system insulating, albeit the discontinuity in $`n(\stackrel{}{k})`$ which existed in the metal has been washed out. This effect is linked to the fact that $`n(\stackrel{}{k})`$ reflects the underlying Fermi statistics of the electronic system. For the specific case of a two dimensional square lattice that resembles the CuO2 planes of the cuprates, there is a drop in $`n(\stackrel{}{k})`$ across a line that is close to the antiferromagnetic zone boundary.

It is believed that the information under the sudden approximation of the momentum distribution function, $`n(\stackrel{}{k})`$, can be extracted from angle-resolved photoemission spectroscopy (ARPES) data via the relation $`n(\stackrel{}{k})`$ = $`A(\stackrel{}{k},\omega )f(\omega )𝑑\omega `$.  Here $`A(\stackrel{}{k},\omega )`$ is the spectral function, and $`f(\omega )`$ is the Fermi function. In a real experiment, the $`n(\stackrel{}{k})`$ structure would be further modified by the photoionization cross-section. In the metallic state of optimally doped Bi2Sr2CaCu2O8+δ(BSCCO) the steepest descent of $`n(\stackrel{}{k})`$ gives a Fermi surface consistent with traditional ARPES analysis methods, despite the complication of matrix elements. The intriguing new result is that the $`n(\stackrel{}{k})`$ pattern of the insulator, Ca2CuO2Cl2(CCOC), is strikingly similar to the $`n(\stackrel{}{k})`$ pattern seen in BSCCO. This realization, coupled with many-body theoretical calculations on various forms of the Hubbard model, suggests that the insulator pattern contains information that is related to a Fermi surface which has been destroyed by strong electron-electron interactions thus giving a qualitative concept of a remnant Fermi surface as the surface of steepest descent in $`n(\stackrel{}{k})`$. Although the detailed shape of the remnant Fermi surface is not crucial, it acts to emphasize a robust feature which we observe in the insulator. While this may not be a rigorous definition, as the Fermi surface is only defined for a metal, this idea allows a practical connection from the pseudogap seen in underdoped cuprates to the properties of the insulator.

A photon energy dependence study of this issue is important to extract the intrinsic $`n(\stackrel{}{k})`$ structure from the raw data which can be affected by the matrix element. Because the single particle spectral function obtained from the half filled CCOC insulator will not suffer from possible complications arising due to intrinsic electronic inhomogeneity, bilayer splitting, or superstructure as in other cuprates, this study may be able to provide insight on some of the recent controversy regarding the Fermi surface seen in BSCCO.

In this report we show that despite the aberrations caused by the matrix element effect, the drop in $`n(\stackrel{}{k})`$ as one crosses the antiferromagnetic zone boundary is a robust feature of CCOC. This result differs from the conclusion of a recent study based on data from only a single cut across the magnetic zone boundary. We also show that the dispersion remains independent of photon energy even when the matrix element does not. This is an expected result from an ideal two dimensional system.

The measurements are conducted at Beamline V-3 of the Stanford Synchrotron Radiation Laboratory (SSRL). The angular resolution is $`\pm 1^{}`$. The CCOC samples are oriented prior to the experiment by the Laue method and are cleaved in situ. The light is incident at $`45^{}`$ to the surface normal and has an in-plane component of the polarization along $`k_x`$, which is also parallel to the Cu-O bonds. The chamber pressure is better than $`5\times 10^{11}`$ torr. Data is taken in the $`\stackrel{}{k}`$-space octant $`(0,0)(\pi ,0)(\pi ,\pi )(0,0)`$ and at temperatures ranging from 100 to 185K, well below the Néel temperature of 250K. The insulating samples are checked for charging by varying the photon flux. All data are normalized by the photon flux which is measured by a Au mesh upstream from the sample.

We begin with a short discussion of matrix element effects in photoemission. Extracting the single particle spectral function from ARPES measurements is complicated by the fact that the measured photoemission intensity is a product of the spectral function and the matrix element. In interacting electron systems, it is impossible to calculate the matrix element exactly, thus further complicating the ARPES analysis. In general it is a function of the experimental geometry, photon energy, and the electronic wave function. However, symmetry arguments can be very powerful in understanding some of the general properties of the matrix element. Since its details are not well understood, the objective in a given photoemission study must be to focus on only those features of the data which are robust against variations in the matrix element.

For the $`n(\stackrel{}{k})`$ plots on CCOC, the spectra are integrated from -0.5 to 0.3 eV binding energy relative to the valence band maximum. Spectra are taken at the crosses, and the data is symmetrized about $`k_x`$=$`k_y`$ to better illustrate the remnant Fermi surface. Figure 1 shows the original comparison between metallic BSCCO and insulating CCOC in panels (A) and (B). In both cases the geometry of the experimental setup is identical. Assuming a $`d_{x^2y^2}`$ orbital symmetry we would expect a vanishing cross section as one approaches the line $`k_x`$=0. This is observed in both data sets. In BSCCO the only drop in intensity which is not naturally explained by matrix elements is where a Fermi surface crossing has occurred. As seen in the figure the intensity drop matches the traditional method for determining a Fermi surface crossing by following the dispersion by eye is indicated by the dots. A similar drop in intensity is observed in the insulator. Although the feature is less well defined here, the striking resemblance it bears to the metal suggests that the origin is similar, and hence it is qualitatively described as a remnant Fermi surface. This allows us to make a natural link between the dispersion in the insulator and the pseudogap in the underdoped regime of the high-Tc cuprates.

As mentioned above, we must still show the robustness of this feature to demonstrate the validity of this concept. Panels B through E of Figure 1 show $`n(\stackrel{}{k})`$ patterns of CCOC taken at photon energies of 25.2, 29, 32.3, and 41eV. At first glance they may appear different. One clearly sees the position of maximum intensity shift. Thus the position of this maximum along the antiferromagnetic zone boundary can be attributed to matrix element effects. However, the expected suppression for a state with $`d_{x^2y^2}`$ symmetry as one approaches the zone center is observed in all the data. In fact, a careful look at the data shows that they are very similar. The exact shape of the remnant Fermi surface may change slightly, but at all photon energies used there is a loss of spectral weight as one crosses the approximate AF zone boundary from $`(0,\pi )`$ to $`(\pi ,0)`$. Panel F plots the remnant Fermi surface as determined by figures 1B through 1E. The shaded region in F gives a rough estimate for the shape and uncertainty of the remnant Fermi surface. Thin black lines in panels B through E form the border of the shaded region shown in F. It may appear that the remnant Fermi surface is more hole like or electron like depending on the photon energy chosen, but the broadness and the minor variability due to matrix elements prevent one from clearly identifying the remnant Fermi surface as either case, as illustrated by the hashed region. However, the point which is clear from the data is that globally, there is always an unexpected loss in intensity as one crosses the region in the Brillouin zone indicated by the shaded area in figure 1F, which represents the remnant Fermi surface. Indeed, the remnant Fermi surface is a robust feature of the insulator.

Although the $`n(\stackrel{}{k})`$ image plots can provide a wealth of information and are extremely good for summarizing data, it is important to look at the raw data to fully appreciate what information is being given by the image plots. This is done in Figure 2. Panels a through d (e through f) plot the EDC’s from (0,0) to $`(\pi ,\pi )`$ ($`(\pi ,0)`$ to $`(\pi /2,\pi /2)`$) taken at 25.2, 29, 32.3, and 41eV photon energy respectively. We find the spectra are qualitatively similar. This is true even at 41eV where the peak is poorly defined throughout the zone. To examine them more closely, Figure 3 plots both the dispersion of the peak position and $`n(\stackrel{}{k})`$ together for the (0,0) to $`(\pi ,\pi )`$ cut. Only slight differences exist between the 4 photon energies. It is apparent that the matrix element can be responsible for shifting the maximum in intensity relative to the maximum in dispersion. Such an effect is seen between the 29 and 41eV data, which may be rigorously compared. A similar shift of spectral weight has been reported by Haffner et al. in Sr2CuO2Cl2. These authors conclude that there is no $`n(\stackrel{}{k})`$ structure or remnant Fermi surface. We believe the intensity shift is a small modulation due to the a matrix element compared with the robust underlying $`n(\stackrel{}{k})`$ structure seen throughout the zone at all photon energies. To illustrate this point we also plot the spectra, dispersion, and $`n(\stackrel{}{k})`$ along $`(\pi ,0)`$ to $`(0,\pi )`$. From the raw spectra in figure 2 and the peak positions shown in figure 3 it is clear that the dispersion does not change. Meanwhile the intensity of the feature varies randomly along this cut. In the extreme case between 25.2eV and 29eV the intensity is increasing as one approaches $`(\pi /2,\pi /2)`$ for the former, and decreasing for the latter. This completely contrasts the situation along (0,0) to $`(\pi ,\pi )`$, where there are only slight shifts of the intensity profile, but the global structure remains the same.

The cuprates are generally believed to be two dimensional electronic systems. In such an ideal case, the electronic states probed by ARPES are independent of photon energy, and hence $`A(\stackrel{}{k},\omega )`$ extracted from the spectra will remain the same. Thus it is natural to question this assumption when photon energy dependent shifts in the physics are reported. This may seem particularly true for the half filled cuprates which possess three dimensional long range magnetic order. However, in figures 2 and 3 we have shown that the dispersion of the insulator is essentially independent of photon energy. Three dimensionality may still play a role in causing a small shift in dispersion with photon energy, but it is a safe assumption to treat the electronic structure of the half filled insulator as basically two dimensional.

Recently, the role of the matrix element effects is also under examination regarding the Fermi surface of BSCCO. Several groups have reported the existence of an electron-like Fermi surface centered at (0,0). These groups report a subtle change of dispersion near the flat band region around $`(\pi ,0)`$, causing a change in the shape of the Fermi surface. Largely due to the difficulty in understanding the effects caused by photoelectron matrix elements, there is no concensus in attempting to resolve this issue. The features in BSCCO are better defined than the relatively broad features in CCOC, and one might suspect that, as a result, they may be more susceptible to photon energy dependence. However, the idea that the dispersion should be independent of photon energy for an ideal two dimensional system is the same, regardless of the resolution. A change in the dispersion can not be a trivial matrix element effect. Possible explanations for a dispersion which changes with photon energy include: dispersion in $`k_z`$, bi-layer splitting, and phase separation. With the exception of Pb doped samples BSCCO also has the additional complication of superstructure in the Bi-O planes. In BSCCO, it would be surprising for the dispersion in $`k_z`$ to have a large effect considering the two dimensional nature of the cuprates. In reality though it is likely that all of these effects play roles which result in a potentially very complex picture.

Interestingly, models which have multiple component electronic structures where the wavefunctions and dispersions differ could have a simple explanation for the observed controversy. In such a case, variations in the matrix element could now lead to different interpretations of the observed ARPES spectra at different photon energies. This is possible because there may now be a different matrix element for each component, and each matrix element may vary differently with photon energy. Thus causing one component to be more dominant at one photon energy while hidden at another. This would naturally lead to different interpretations as a function of photon energy.

We have shown that the loss in intensity as one crosses the antiferromagnetic zone boundary is a robust feature of the insulator, Ca2CuO2Cl2, which can not be explained solely by matrix element effects. However, the photon energy dependence does underscore the qualitative rather than quantitative nature of the remnant Fermi surface concept. We argue that much physics can be learned in spite of the effect which matrix elements can have in ARPES, as long as care is taken to properly sort out the intrinsic versus the extrinsic physics. In particular, the resulting connection between the d-wave like dispersion and the pseudogap in the underdoped regime is robust.

We thank C. Dürr, M.S. Golden, and J. Fink for very open and stimulating discussions. This work was performed at Stanford Synchrotron Radiation Laboratory which is supported by the DOE Office of Basic Energy Science, Division of Chemical Sciences and Material Sciences. It is also supported by the ONR grant N00014-98-1-0195 and NSF grant DMR-9705210.