text:
# Delocalized Quasiparticles in the Vortex State of an Overdoped High-𝑇_𝑐 Superconductor Probed by 63Cu NMR

\[

## Abstract

We report the spin Knight shift ($`K_s`$) and the nuclear spin-lattice relaxation rate ($`1/T_1`$) in the vortex state as a function of magnetic field ($`H`$) up to 28 T in the high-$`T_c`$ superconductor TlSr2CaCu2O6.8 ($`T_c`$=68 K). At low temperatures well below $`T_c`$, both $`K_s`$ and $`1/T_1`$ measured around the middle point between two nearest vortices (saddle point) increase substantially with increasing field, which indicate that the quasiparticle states with an ungapped spectrum are extended outside the vortex cores in a d-wave superconductor. The density of states (DOS) around the saddle point is found to be $`\kappa N_0\sqrt{H/H_{c2}}`$, with $`\kappa `$=0.5$``$0.7 and $`N_0`$ being the normal-state DOS.

\]

The vortex state in unconventional superconductors, such as the high transition-temperature (high-$`T_c`$) superconducting copper oxides which have a $`d_{x^2y^2}`$-wave gap symmetry , is of great interest. Among various issues, the quasiparticle state associated with the vortices still remains far from being well understood. In an s-wave superconductor with isotropic gap, the low-lying quasiparticle states are localized within the vortex core whose radius is $`\xi =v_F/\pi \mathrm{\Delta }_0`$, with energies $`E=\mu \mathrm{\Delta }_0^2/E_F`$ ( $`\mu `$=1/2,3/2….) , where $`\mathrm{\Delta }_0`$ is the energy gap, $`v_F`$ and $`E_F`$ are the velocity and the energy level at the Fermi surface, respectively. This prediction was verified experimentally by the scanning tunneling microscopy (STM) . In contrast, the gap for d-wave superconductivity vanishes in the nodes which leads to a divergent $`\xi `$ along the nodal directions. The corresponding quasiparticle state is, therefore, expected to have very different properties. Recent theories suggest that the quasiparticle states may be extended far outside the cores along the nodal directions . However, this prediction has not yet been unambiguously confirmed by experiments. Heat capacity measurement at low fields in the vortex state of YBa2Cu3O7 (YBCO) has revealed a $`\sqrt{H}`$ dependence of the specific heat coefficient , which was taken as evidence for supporting the theoretical prediction . But heat capacity also contains contributions from the vortex cores which may have a $`\sqrt{H}`$-like variation in the low field regime as observed in s-wave superconductors and pointed out by Sonier et al . Our earlier NMR measurement on an underdoped cuprate  and thermal transport measurement  in the low $`H/H_{c2}`$ regimes are consistent with the existence of finite density of quasiparticle states outside the cores, but no quantitative estimation of the field-induced states has been available. The most serious question about the quasiparticle picture in a d-wave superconductor is raised recently by STM experiments that find no quasiparticle states outside the cores . Clearly, more spatially-resolved and bulk-sensitive studies, and also at higher fields, are needed to resolve this issue.

We report in this Letter measurements of the spin susceptibility and the nuclear spin-lattice relaxation rate ($`1/T_1`$) outside the vortex cores as a function of $`H`$ up to $`H`$=28 T (or $`H/H_{c2}`$=0.65) in a d-wave high-$`T_c`$ cuprate. We find compelling evidence for the quasiparticles in the vortex state being delocalized. Unlike specific heat measurement that picks up contributions from quasiparticles both inside and outside the vortex cores, NMR technique can measure separately the quasiparticle states in- and out-side the vortex core. In the mixed state, those nuclei located in between two nearest vortex cores have the largest distribution probability that gives rise to a singularity in the field distribution function, while those in the cores contribute to a ”knee” that is far away from the singularity . In practice, however, both the singularity and the knee are usually smeared out owing to, e.g.  lattice disorders, resulting in a more symmetric NMR line shape. Nevertheless, some information about the position selectiveness is retained; the peak of the NMR line that results from smeared singularity corresponds to the nuclei located outside the vortex cores and in the superconducting state.

The sample chosen for this study is a slightly overdoped cuprate TlSr2CaCu2O6.8 ($`T_c`$=68 K). This compound has several advantages over the most-extensively-studied YBCO family with regards to our purpose. First, the Knight shift arising from the spin susceptibility when field is applied parallel to c-axis is large , which allows an accurate measurement of its change in the superconducting state in the presence of magnetic fields. This is to be compared to the YBCO family where the corresponding shift is negligibly small because of the accidental cancellation of the hyperfine field, so that measurement along this direction is not possible . Second, the spin susceptibility in the normal state above $`T_c`$ is almost temperature independent, which helps to evaluate the density of states (DOS) in the normal state. Third, the upper critical field $`H_{c2}`$ 43 T is relatively small , so that the field-induced DOS at a given $`H`$, if any, should be larger. Finally, there is only one crystallographic site of Cu, which greatly simplifies NMR spectroscopy and enhances data accuracy.

The polycrystals were c-axis aligned under a field of 16 T and fixed with a polymer. For all measurements, the external field is applied along the c-axis. Figure 1 shows a typical NMR line shape at $`T`$=4.2 K (gray curve). The dotted curve is the corresponding theoretical field distribution for 75 vortex lattice (VL) obtained as described below. The position-dependent field $`H(r)`$ in the sample was calculated by using the London model ,

$`H(r)`$ $`=`$ $`\overline{H}{\displaystyle \underset{l,m}{}}{\displaystyle \frac{\mathrm{exp}(G_{lm}^2\xi ^2/2)\mathrm{exp}(i\stackrel{}{G_{lm}}\stackrel{}{r})}{1+G_{lm}^2\lambda ^2}}`$ (1)
where $`\overline{H}`$ is the averaged (externally applied) field, $`\xi `$ and $`\lambda `$ are the coherence length and penetration depth, respectively. The summation runs over all reciprocal VLs $`\stackrel{}{G_{lm}}=2\pi \sqrt{\frac{\overline{H}Sin\beta }{\varphi _0}}(m\widehat{x}+\frac{1}{Sin\beta }(lmCos\beta )\widehat{y})`$, where $`\beta `$ is the angle between two primitive VL vectors and $`\widehat{x}`$ and $`\widehat{y}`$ are the unit vectors of the reciprocal VL. We adopt $`\beta =75^{}=5\pi /12`$, which is the value found in neutron scattering  and STM  experiments. The NMR spectrum is the spatial average of a Dirac delta function, $`f(H^{})=\delta (H^{}H(r))_r`$. By use of the relation of $`H_{c2}=\varphi _0/2\pi \xi ^2`$, eq. (1) is left with a single parameter, $`\lambda `$, that is to be determined by experiment. In the theoretical curve of Fig. 1 (dotted curve) calculated from eq. (1), the singularity $`S`$ corresponds to the middle point between two nearest vortices (saddle point), $`M`$ and $`m`$ correspond to the vortex core and the middle point between four vortices, respectively. The solid thick curve is a calculation by convoluting eq. (1) with a Lorentzian broadening function $`p(H^{})=\frac{\sigma ^2}{\sigma ^2+4\pi ^2H^2}`$, with a fixed value of $`\sigma /\pi `$=170 Oe that is the full width at half maximum (FWHM) of the NMR spectrum at $`T`$=70 K (above $`T_c`$). In performing the fitting, both the experimental spectrum and the theoretical field distribution curve are normalized to unity. From the peak position, we determined the Knight shift, $`K_c`$, around the saddle point. We also obtained $`\lambda `$=950 $`\pm `$ 100$`\AA `$ at $`T`$=4.2 K, and hence the diamagnetic field $`H_{dia}=H_{saddle}\overline{H}`$ that depends on the magnitude of externally applied field. Similar value of $`\lambda `$ was obtained from the analysis of 203Tl lineshape which has a narrower width, with FWHM = 50 Oe at $`T`$=70 K and 140 Oe at $`T`$=4.2 K, respectively.

The value of $`\lambda `$=950$`\pm `$100 $`\AA `$ is slightly smaller than that estimated for the YBCO family (1200$``$1400 $`\AA `$) . This is consistent with the current compound being overdoped. Also, we have confirmed that the spatially dependent information is indeed retained, even though the singularity and the knee are smeared out. This was done by measuring $`1/T_1`$ in different positions. $`1/T_1`$ is larger by a factor of two in the higher end of the spectrum shown in Fig. 1 around the knee, than that measured at the peak.

Now we turn to the main result, the $`H`$ dependence of the Knight shift and $`1/T_1`$ measured at the spectrum peak. Figure 2 shows the typical data sets for the $`T`$ and $`H`$ dependencies of $`K_c`$. In the normal state $`K_c`$ shows very weak $`T`$-dependence and does not depend on $`H`$ above $`T`$=85 K. In particular, $`K_c`$ is constant (=1.42$`\%`$ ) in the temperature range of 85K $`T`$ 150 K. In contrast, at low temperatures below $`T_c`$, the shift strongly depends on $`H`$. This is clearly shown in Fig. 3 where data at $`T`$=4.2 K are plotted. The solid circles show the data corrected by the diamagnetic shift $`H_{dia}/\overline{H}`$ obtained as described earlier. It can be seen clearly that $`K_cK_{dia}`$ increases with increasing field. This corrected shift is composed of the spin part ($`K_s`$) and the part due to orbital susceptibility $`K_{orb}`$, respectively, the latter being $`H`$-independent in the range of field studied. Note that $`K_s`$ is related to the spin susceptibility and hence the DOS, $`N(E)`$ as $`K_sN(E)(\frac{f(E)}{E})𝑑E`$, where $`f(E)`$ is the Fermi distribution function. In particular, $`K_s`$ at low $`T`$ is proportional to the DOS at the Fermi level, $`N(0)`$. Thus, our result indicates that magnetic field induces a finite $`N(0)`$, in the superconducting state.

This conclusion is corroborated by the $`H`$-induced increase of $`1/T_1`$ measured at the spectrum peak. Figure 4 shows the typical data sets of $`1/T_1`$ under various fields from $`H`$=0 to 28 T. The measurement was done by using a single saturation pulse  and at the NQR spectrum peak of 22.5 MHz for $`H`$=0. As seen in the figure, $`1/T_1`$ becomes proportional to $`T`$ below $`T`$ 10 K, and its magnitude increases with increasing $`H`$. This is seen more clearly in Fig. 5 where $`1/T_1T`$ at $`T`$=4.2 K is displayed. In terms of DOS, $`T_1`$ is expressed as $`1/T_1=\alpha (T)N(E)^2f(E)(1f(E))𝑑E`$, where $`\alpha (T)`$ is the enhancement factor due to antiferromagnetic (AF) correlation. For a d-wave superconductor , $`\alpha (T)`$ was shown to have a very weak $`T`$ variation below $`T_c`$ and $`\alpha (0)0.9\alpha (T_c)`$. Therefore, $`1/T_1TN(0)^2`$ at low $`T`$; the $`T_1`$ result also indicates that $`H`$ induces a finite DOS outside the vortex cores.

Next we evaluate more quantitatively the DOS around the saddle point, $`\delta N_H`$. First, one notices from Fig. 5 that at high fields of $`H`$5 T, $`1/T_1`$ is linear in $`H`$ which suggests that the induced DOS is proportional to $`\sqrt{H}`$. In fact, in this field regime, a functional of

$`\delta N_H`$ $`=`$ $`\kappa N_0\sqrt{H/H_{c2}}`$ (2)
can fit both the Knight shift and $`1/T_1`$ quite well, where $`N_0`$ is the DOS in the normal state. In Fig. 3, the solid curve is a fit to $`K_cK_{dia}=K_0+K_H\sqrt{H/H_{c2}}`$, with resultant fitting parameters $`K_0=1.00\pm 0.03\%`$ and $`K_H=0.32\pm 0.05\%`$. The solid line in Fig. 4 is $`\frac{1}{T_1T}=0.018+3.175\frac{H}{H_{c2}}`$ \[Sec-1K-1\]. On the other hand, $`1/T_1T`$ remains finite at $`H`$=0 and rises more slowly with increasing $`H`$ at low fields below 5 T. We ascribe the zero-field residual $`1/T_1T`$ to that coming from the residual DOS, $`\delta N_{imp}`$, due to scattering by impurity or crystal imperfection, which is estimated as $`\frac{\delta N_{imp}}{N_0}=\sqrt{\frac{(T_1T)_{Tc}}{0.9(T_1T)_{4.2K}}}`$=0.14. Then the slower rise of $`1/T_1T`$ below 5 T can be naturally understood as vortex flux preferentially occupying the impurity sites where a finite DOS has already been produced . In the impurity-dominated regime, a recent calculation  shows that the field-induced DOS has an asymptotic form of $`\delta N_H\frac{H}{H_{c2}}ln(\frac{H_0}{H})`$ where $`H_0`$ is a constant ($`\frac{E_f\mathrm{\Delta }_0}{900}`$K-1). Indeed, the data of $`H`$5 T can be fitted by such functional as seen in Fig. 5, where the resulting dotted curve is $`\frac{1}{T_1T}=0.229+1.650(\frac{H}{H_{c2}}ln\frac{270}{H})^2`$ in unit of Sec-1K-1. The same $`\delta N_H`$ functional fitting to the shift of $`H`$5 T with $`H_0`$=270 T predicts, as shown by the dotted curve in Fig. 3, the shift to remain at 1.03% when $`H`$=0. Since this value is a sum of the impurity-induced spin shift and $`K_{orb}`$, namely, $`1.03\%=\frac{\delta N_{imp}}{N_0}(1.42\%K_{orb})+K_{orb}`$, we obtain $`K_{orb}`$=0.96% and hence $`K_s`$=0.46%, with uncertainties of $`\pm `$0.03% (absolute value) . Then, $`\kappa =K_H/K_s`$ is estimated to be 0.70$`\pm 0.06\%`$. On the other hand, from the $`T_1`$ result and following the same procedure of estimating $`\delta N_{imp}`$, the prefactor $`\kappa `$ is estimated to be 0.50$`\pm 0.03\%`$. This difference of $`\kappa `$ by 30% might arise from the choosing of the AF enhancement factor $`\alpha (0)0.9\alpha (T_c)`$ that underestimates $`\delta N_H`$ from the $`T_1`$ result. The analysis of the Knight shift data could also overestimate $`\kappa `$; if a triangular model is assumed for the VLs, $`\kappa `$ would be reduced by $``$0.1.

To summarize the above discussions, it is important to stress that both the $`H`$-enhanced Knight shift and $`1/T_1`$ observed around the saddle point can be consistently accounted for by a DOS of $`\delta N_H=\kappa N_0\sqrt{H/H_{c2}}`$, with a conservative estimate of $`\kappa =0.5`$ 0.7. We therefore conclude that the quasiparticle states are extended outside the vortex cores and that they have a continuous excitation spectrum without a gap. This feature is in sharp contrast with that in an s-wave superconductor. Volovik first pointed out that the spatially averaged DOS associated with a single vortex in a d-wave superconductor is $`\delta N_H=\kappa N_0\sqrt{H/H_{c2}}`$, where $`\kappa `$ is on the order of unity . Although no calculation of the DOS outside the vortex cores is available to allow a direct comparison to our results, the substantially smaller $`\kappa `$ than unity found in our experiment could be due to the effects of vortex lattice, its disorder , or due to induction of second component of the order parameter .

In conclusion, from $`1/T_1`$ and the spin Knight shift measured around the saddle point in a slightly overdoped high-$`T_c`$ cuprate, we find that the quasiparticle states are extended outside the vortex cores with ungapped spectrum in a d-wave superconductor, which is in contrast to the localized states within the vortex cores in s-wave superconductors. The density of such quasiparticle states near the saddle point, in the field range of 5 T$`<H<`$28 T, follows $`\kappa N_0\sqrt{H/H_{c2}}`$ where $`N_0`$ is the normal-state DOS. The prefactor $`\kappa `$ is estimated to be in the range of 0.5$``$0.7. Although such states have not been found in the STM experiments, one may be able to reconcile the STM and NMR results by considering the momentum-dependent tunneling matrix in the STM experiments which likely vanishes in the nodal directions .

We thank W.G. Clark, Y. Shimakawa and Y. Kubo for collaborations in the early stage of this work. G.-q.Z thanks in particular A.V. Balatsky, Y. Matsuda, M. Ogata and S.H. Pan for several useful conversations, and partial support by a Grant-in-Aid for Scientific Research No. 11640350 from MEXT of Japan and by the Visiting Scientist Program of National High Magnetic Field Laboratory which is supported by NSF (DMR-9527035) and by the State of Florida. Additional support from COE research grant No. 10CE2004 of MEXT is acknowledged.


format: text

text:
# Pulsations and orbital modulation of the intermediate polar 1WGA J1958.2+3232

## 1 Introduction

Intermediate polars (IPs) are semi-detached interacting binaries in which a magnetic white dwarf accretes material from a Roche-lobe filling, usually late-type, main sequence companion star. The accretion flow from the secondary proceeds towards the white dwarf either through an accretion disc, an accretion stream, or some combination of both (known as disc overflow accretion), until it reaches the magnetospheric radius. Here the material attaches to the magnetic field lines and follows them towards the magnetic poles of the white dwarf. The infalling material that originates from an accretion disc takes the form of arc-shaped accretion curtains, standing above the white dwarf surface. At some distance from this surface, the accretion flow undergoes a strong shock, below which material settles onto the white dwarf, releasing X-ray to optical emission. Since the magnetic axis is offset from the spin axis of the white dwarf, this gives rise to the defining characteristic of the class, namely X-ray (and usually optical) emission pulsed at the white dwarf spin period. Additionally, as the X-ray ‘beam’ sweeps around the system, there is the possibility that some fraction of the emission will be reprocessed from structures such as the companion star or a bulge at the edge of an accretion disc. This will give rise to further optical emission pulsed at the lower orbital sideband of the spin frequency, namely the spin frequency of the white dwarf in the reference frame of the binary. Several IPs show a dominant optical pulsation at this orbital sideband frequency. X-ray pulsation at the orbital sideband frequency is due to an intrinsic modulation arising from pole-switching in the case of stream-fed accretion. Several IPs show this type of modulation (at least at some epochs) too. About twenty confirmed intermediate polars are now recognised with a similar number of candidate systems having been proposed. Comprehensive reviews of various aspects of their behaviour are given by Patterson (Patt (1994)), Warner (War95 (1995)), Hellier (Hell95 (1995); Hell96 (1996)) and Norton (Nor95 (1995)).

## 2 Observational history

1WGA J1958.2+3232 was first identified as a $`12`$ min X-ray pulsator by Israel et al (Is98 (1998)) from archival ROSAT PSPC light curves. Israel et al (Is99 (1999)) then identified the optical counterpart to the system as a 16th magnitude star, and claimed it was of spectral type B0Ve, indicating that the system was a Be/X-ray binary. Negueruela et al (Neg00 (2000)) subsequently showed that this optical counterpart had a spectrum consistent with that of a magnetic cataclysmic variable – an intermediate polar. Uslenghi et al (Us00 (2000)) later obtained time resolved white-light optical photometry of the star and discovered a $`12`$ min optical modulation with the same period as the X-ray modulation. They also saw some evidence that the pulse period may be twice this value. Next, Zharikov et al (Zha01 (2001)) obtained further time resolved R-band photometry and spectroscopy from which they detected a possible orbital period of 4.36 h as well as confirming the pulsation period as 733 sec. Most recently, Uslenghi et al (Us01 (2000)) detected circular polarization from the source in the R and I bands, with evidence for possible modulation of the polarization at twice the previously observed pulsation period.

## 3 Observations and data reduction

### 3.1 Optical photometry

Time resolved photometric data were obtained at the 1 m Jacobus Kapteyn Telescope (JKT) between 9–15 July 2000, using the SITe2 CCD and UBVRI filters. For the majority of the observations, the CCD was windowed to readout a 900 $`\times `$ 400 pixel region around the target, corresponding to a field-of-view of $`4.5^{}\times 2.0^{}`$. Exposures in different filters were interleaved each night, with a dead-time between exposures (for readout and filter change) of 30–40 s. Landolt standard fields were also observed through the relevant filters on each night. The observing log is shown in Table 1.

Bias subtraction and flat-fielding were carried out using routines in iraf on all target and standard frames. Aperture photometry using further iraf routines was then carried out on the target as well as a number of near-by comparison stars, indicated in Fig. 1. Lightcurves containing differential photometry of 1WGA J1958.2+3232 with respect to each of the comparison stars were constructed. Extinction coefficients were then derived for each filter used on each night by examining the variation in observed magnitude of the comparison stars as a function of airmass. By combining these with the observations of Landolt standards, this in turn allowed the derivation of zero points in each filter for each night. Finally, these corrections were applied to the lightcurves of the target, resulting in the calibrated lightcurves shown in Fig. 2.

### 3.2 Optical photo-polarimetry

Time resolved B- and R-band circular photopolarimetry of 1WGA J1958.2+3232 was obtained at the 2.56m Nordic Optical Telescope (NOT) between the nights of 28 February and 2 March, 2001. The Andalucia Faint Object Spectrograph and Camera (ALFOSC) was used with a 1/4-waveplate and a calcite block. The o and e images of the star were separated by 15 arc sec such that the images of 1WGA J1958.2+3232 do not overlap with other stars in the field.

The R-band data of 1WGA J1958.2+3232 were obtained on 1 March 2001 between UT 5.94 and UT 6.81, covering 2.14 spin cycles whilst the B-band data were obtained on 3 March 2001 between UT 5.86 and UT 6.78, covering 2.25 spin cycles. The seeing was between 1 and 1.5 arc seconds on both nights and exposure times were 58 sec for the B-band and 45 sec for the R-band images. On each night all the images were taken as one long sequence of exposures, giving a time-resolution of less than one minute.

The CCD-frames were bias-subtracted and flat-fielded with standard iraf-routines and aperture photometry was carried out with daophot (Stetson 1992). Measurements were also made for some other bright stars in the field of 1WGA J1958.2+3232 to calibrate the zero-point of the circular polarization. The sign of the circular polarization was checked by observing polarization standards $`\mathrm{Grw}+70^{}8247`$ and LP 790–29. The instrumental polarization was found to be slightly dependent on the position of the star in the CCD-image. Corrections (1.2% in the B-band, and 0.9% in the R-band) for instrumental polarization were made based on measurements of unpolarized field stars near the position of the object in the CCD-frame.

### 3.3 Optical spectroscopy

Two spectra of 1WGA J1958.2+3232, each spanning the wavelength range $`36006700`$Å, were obtained at the 2.5m Isaac Newton Telescope (INT) on 21 July 2000 using the Intermediate Dispersion Spectrograph, 500mm camera, R150V grating and TEK5 CCD. This combination yields a dispersion of 3.17Å pixel-1. The exposure time for each spectrum was 1800s and the two exposures began at 02:02UT and 02:34UT respectively. The resulting frames were bias subtracted then flat-fielded using tungsten lamp flats; the spectra were optimally extracted and calibrated using CuAr arcs. Fig. 3 shows the average spectrum obtained, with the continuum normalised to unity, and emission lines identified.

The average INT spectrum shows a similar range of emission lines to that reported by Negueruela et al (Neg00 (2000)). The Balmer series lines up to at least H$`\zeta `$ are all in emission, as are lines due to HeI (4472Å, 4922Å, 5876Å) and HeII (4686Å, 5412Å). Unlike the spectra reported by Negueruela et al (Neg00 (2000)) however, the lines are not double peaked.

### 3.4 X-ray observations

We have extracted the ASCA (Tanaka, Inoue & Holt 1994) data of 1WGA J1958.2+3232 from the HEASARC. The observation was carried out on 1998 May 15 over an elapsed time of about 15.5 hrs. Data were taken with all 4 co-aligned instruments, two GIS’s and two SIS’s, all in standard configuration. We have applied the standard screening criteria (excluding Earth occultation, passages through the South Atlantic Anomaly, and times of high background), resulting in $``$30 ksec of good data in each instrument. We have extracted light curves in 16s bins from each instrument in low (0.4–2 keV for SIS, 0.7–2 keV for GIS) and high (2–10 keV) energy bands, both from the source and from nearby source-free background region. We have scaled the latter using the detector areas, and subtracted from the former, to create net light curves; we have then summed the net light curves in the low and the high bands, and also created a total (0.4/0.7–10 keV) light curve.

## 4 Results

The average magnitudes and colours of 1WGA J1958.2+3232 derived from our seven nights of JKT data are listed in Table 2. formal uncertainties are dominated by the intrinsic variability of the source. We note that these magnitudes are consistent with the UBV values reported by Uslenghi et al (Us00 (2000)) and obtained about 1 year before our data.

### 4.1 Optical power spectra

The JKT lightcurves in magnitudes were converted into flux units prior to time series analysis. Fig. 4 shows the power spectrum of the longest dataset, i.e. the V-band lightcurve, in two panels containing the low and high frequency regimes. Frequencies corresponding to the peaks were measured by fitting Gaussian profiles to the relevant portions of the power spectra.

The power spectrum is dominated by three sets of signals: one set centred around a frequency of $`117.8`$ d-1; one set around $`55.5`$ d-1; and a final set around $`5.5`$ d-1. Each set of peaks consists of a series of sidebands due to daily aliases either side of a central strong peak. Clearly these three sets of frequencies are related to the pulse period, roughly twice that period, and the orbital period of the system as previously detected. There is also the possibility that one of the short periods is related to a beat between the white dwarf spin and system orbital periods. However, it is important to note that the strongest peaks in each of the three sets are not harmonically related to each other. A solution to the identification of the frequencies is as follows.

* We identify the strongest peak in the highest set of frequencies with twice the spin frequency of the white dwarf, i.e. $`2\omega =117.819\pm 0.005`$ d-1 or $`P_{\mathrm{spin}}/2=733.33\pm 0.03`$ s. This is identical with the pulse period previously detected by Israel et al (Is98 (1998)) and Uslenghi et al (Us00 (2000)).
* We identify one of the peaks in the middle set of frequencies as the beat frequency between the white dwarf spin and system orbital frequencies. However, exactly which of the peaks is to be identified with the beat frequency is not immediately clear.
* We identify one of the peaks in the lowest set of frequencies with the orbital frequency of the system. However, as above, exactly which of the peaks is to be identified with the orbital frequency is not immediately clear.

The strongest peak in lowest set of frequencies corresponds roughly with the orbital period claimed by Zharikov et al (Zha01 (2001)), i.e. $`\mathrm{\Omega }5.5`$ d-1 or $`P_{\mathrm{orb}}4.4`$ h. However, if this is chosen, it implies that the beat frequency is the –2 day alias of the strongest peak in the middle set of frequencies, and this peak clearly has a very low significance since there are many other peaks of comparable strength in the spectrum. Conversely, if the strongest peak of the middle set is the beat frequency, it implies that the orbital frequency is the –2 day alias of the strongest peak in the lowest set. This peak too has a relatively low significance. Given the relative weakness of the –2 day aliases of the strongest peaks in each of these two sets of frequencies, the most likely combination of peaks to represent the system frequencies is as follows.

* The –1 day alias of the strongest peak in the lowest set of frequencies is the orbital frequency of the white dwarf, $`\mathrm{\Omega }=4.455\pm 0.005`$ d-1 or $`P_{\mathrm{orb}}=5.387\pm 0.006`$ h. This is the –1 day alias of the orbital period previously detected by Zharikov et al (Zha01 (2001)).
* The –1 day alias of the strongest peak in the middle set of frequencies is the beat frequency between the white dwarf spin and the orbit, $`(\omega \mathrm{\Omega })=54.434\pm 0.005`$ d-1 or $`P_{\mathrm{beat}}=1587.24\pm 0.14`$ s.

More extreme aliases combinations are probably unlikely, since the power at these aliases is low, although such combinations are possible in principle.

Fig. 5 shows the JKT power spectra for the low frequency regime in each of the U, B, V, R and I bands. The orbital frequency detected in the V-band (indicated by the dotted line) is also seen clearly in the R and B-band power spectra. Broad peaks in the I and U-band power spectra are consistent with the same period, possibly shifted to higher frequencies, but this apparent shift is probably an artifact of the small amount of data in these bands. We therefore detect the orbital period of the system in all five wave-bands.

Fig. 6 shows the same JKT power spectra, but this time in the higher frequency regime, for all five wave-bands. Both the B and I-band lightcurves contain signals consistent with twice the spin frequency (i.e. $`2\omega `$) seen in the V band, indicated by the right-most dotted line. In the U band, broad peaks appear either side of this location, whilst in the R band, no signal is apparent. The two left-most dotted lines indicate the proposed spin frequency ($`\omega `$) and the beat frequency ($`\omega \mathrm{\Omega }`$). The beat frequency is not clearly detected in any waveband other than V, whilst the spin frequency itself may be detected in the R band, and possibly in the B band.

### 4.2 Ephemerides

From the phasing of the peaks in the Fourier transform of the JKT V-band lightcurve, we are able to determine spin and orbital ephemerides, based on our observations of:

Half the spin period:

$$\mathrm{HJD}=2451739.318445(2)+0.0084876(4)\mathrm{E}$$
(1)
Orbital period:

$$\mathrm{HJD}=2451739.36435(2)+0.2245(3)\mathrm{E}$$
(2)
where the zero phase is at minimum light in each case.

We can also combine our determination of the pulse period with earlier measurements to estimate the spin down rate of the white dwarf. The most accurate previous determination of the pulse period is that by Uslenghi et al (Us00 (2000)) who measured it as $`733.24\pm 0.02`$s in June 1999. Combined with our determination of $`733.33\pm 0.03`$s from July 2000, this yields a rate of change for the period of:

$$\dot{P}=(2.6\pm 1.2)\times 10^9$$
(3)
Of the other confirmed intermediate polars, PQ Gem, V1223 Sgr and AE Aqr are also spinning down, but their $`\dot{P}`$ values are at least two orders of magnitude smaller than this (e.g. Patterson Patt (1994); see also http://lheawww.gsfc.nasa.gov/users/mukai/iphome /iphome.html). Similarly, the systems which are spinning up (BG CMi, DQ Her, EX Hya, AO Psc, GK Per) have $`\dot{P}`$ values which are also of the order $`10^{11}`$. FO Aqr has been seen to change from spin down to spin up over the last several years, but again the rate of period change is much smaller than the apparent change in 1WGA J1958.2+3232. Indeed the inferred $`\dot{P}/P`$ for 1WGA J1958.2+3232 is of the order of $`10^4`$ yr-1, which is what is expected for a neutron star accretor rather than a white dwarf (e.g. King & Williams KW83 (1983)). It is likely therefore that the measured pulse periods for 1WGA J1958.2+3232 are not as accurate as claimed.

### 4.3 Optical modulation profiles

The JKT U, B, V, R and I lightcurves were folded, using the ephemerides above, at half the proposed spin period, at the proposed beat period, and at the proposed orbital period of the system. After folding and binning, fluxes were converted back to magnitudes, yielding the profiles shown in Fig. 7. The orbital modulation is clearest in the V-band, but a similar profile with similar phasing is also apparent in the B and R-bands. The orbital profiles in the U and I-bands are more sparse but also consistent with the V-band profile. The modulation at half the proposed white dwarf spin period is clearest in the B and V-bands. The other bands are too noisy to draw any conclusions. Although an apparent modulation is seen in the I-band, it is out of phase with those in the B and V-bands. Of the folds at the beat period, only those in the B and V bands are convincing, although the profile in the R band is also reasonable. Much of the noise evident in all the U and I band folded profiles is due to the relatively small amount of data in these bands.

### 4.4 Optical polarimetry

The extent of the polarimetric lightcurves are not sufficient to allow any period searching. Instead, the reduced B- and R-band circular polarization data and light curves were averaged into 12 phase bins over one spin cycle ($`P_{\mathrm{spin}}=1466`$ sec), according to the ephemeris in Section 4.2. The light curves shown in Figures 8a and 9a represent differential magnitudes with respect to the bright field star H (see Figure 1). We note that these B and R band pulse profiles appear to show a single peak when folded on the inferred spin period, rather than the double peaked pulse profile that would be expected from the photometry reported previously and elsewhere in this paper. However, since the data span only just over two spin cycles, the profiles may not be as representative as other more extensive data sets.

The phase binned B-band curves (Figure 8) show that the circular polarization is about +0.6% between spin phases 0.3 and 0.8, contemporaneous with a 0.1 magnitude fading, whilst the rest of the time the polarization is close to zero. The mean circular polarization is $`+0.32\%\pm 0.10\%`$, a clear $`3\sigma `$ detection. In the R-band, the circular polarization is negative (Figure 9b) with a mean level of $`1.99\%\pm 0.11\%`$. At phase 0.0 the level of polarization is around –1% and the level increases to –5% (near spin phase 0.4). Just after phase 0.5 the circular polarization drops quickly to the –1% level.

We also note that the circular polarization light curves of the flux comparison star H are consistent with zero in both wavebands, varying typically between $`\pm 0.1\%`$ in both B and R.

### 4.5 X-ray pulsation

Figure 10 shows the power spectrum of the ASCA X-ray lightcurve. A strong signal is clearly detected at a frequency of $`2\omega `$, but no other harmonics or frequencies related to either the orbital or beat frequencies of 1WGA J1958.2+332 are seen. The ASCA data in two energy bands folded at both the pulse period and twice the pulse period (i.e. the spin period of the white dwarf) are shown in Figures 11 and 12. There is only a slight difference between the two pulses apparent in Figure 12, as might be expected from the lack of a sub-harmonic in the power spectrum. 1WGA J1958.2+3232 appears to be a classic ’double-peaked’ pulse profile intermediate polar. The modulation depth is close to 100% in the soft energy band, and is less in the hard band, indicating that photoelectric absorption is a major contributor to the overall profile. There are indications of some structure to the pulse profiles, beyond a simple sinusoidal modulation. For instance, the ‘notches’ apparent near the pulse peaks are reminiscent of those seen in AO Psc and V1223 Sgr (Norton 1995). However, in this case they are only marginally significant given the error bars shown in the profiles.

## 5 Discussion and Summary

Optical photometry from the JKT has allowed us to confirm the detection of an optical pulsation in 1WGA J1958.2+3232 at a period of 733.33 s and define a refined ephemeris. ASCA observations also confirm the presence of an X-ray pulsation at the same period as seen in the optical. Our detection of a second optical pulsation period of 1587.24 s, identified with the beat period of the system, implies that the spin period of the white dwarf is twice the short pulsation period, i.e. $`P_{\mathrm{spin}}=1466.66`$ s and that the orbital period is the –1 day alias of the strongest low frequency peak in our power spectrum, i.e. $`P_{\mathrm{orb}}=5.387`$ h. This is also the –1 day alias of the previously reported orbital period. We detect an orbital modulation in each of the U, B, V, R and I bands for the first time.

We also confirm the presence of circularly polarized emission from this source, which is thus only the fifth IP to exhibit such behaviour. The circular polarization is negative in the R-band and positive in the B-band, and in each case shows evidence for variation across the spin cycle. The level of polarization in the R-band, at $`1\%`$ to $`5\%`$, is greater than that seen by Uslenghi et al (2001) who saw a mean level of $`0.5\%`$, whilst our detection of positive circular polarization in the B-band is a first detection in that band.

Since the spin period of the white dwarf is twice the pulsation period observed in both optical and X-ray flux, the pulse profile is actually double peaked. Following Norton et al (Nor99 (1999)), a double peaked X-ray pulse profile indicates that two-pole accretion in 1WGA J1958.2+3232 conspires to produce two peaks per revolution of the white dwarf. In this model, the white dwarf has a relatively low magnetic field strength, and consequently an accretion disc which is truncated only relatively close to the white dwarf. This leads to large footprints of the accretion curtains and an optical depth to X-rays across the accretion regions which is largest in a direction parallel to the white dwarf surface, and smallest perpendicular to the surface. (This is the reverse of the standard accretion curtain model which leads to a single peaked X-ray pulse profile.) Confirmation of this model would require pulse phase resolved optical spectroscopy in order to measure the direction of flow of the accreting material corresponding to the phases of pulse maximum and minimum. Since this model also predicts that the accretion flow is essentially via a disc, an X-ray beat modulation would not be expected to be present, in accord with all observations so far, including the ASCA data presented here.

The different polarities of the circular polarization observed in the B- and R-band data, together with light curve variations, also suggest two pole accretion. We suggest that between spin phases 0.0 and 0.5 the negative polarization accretion pole is seen, and after phase 0.5 emission from this pole reduces. Between spin phase 0.3 and 0.8, the positive polarization accretion pole is seen in the B-band. The different colours of the two accretion poles may indicate that the poles have significantly different magnetic field strengths, possibly as a result of an offset dipole field structure. We note however, that the detection of circular polarization from an IP is generally taken to be an indication of a relatively strong magnetic field, in contradiction to the implication of the double-peaked pulse profile model. There is though no evidence for a strong magnetic field in this system. The figure of $`B=8`$ MG stated by Uslenghi et al (2001) is merely a value chosen for illustrative purposes and is in no way ‘derived’ from their polarization measurements. We suggest that the magnetic field strength of 1WGA J1958.2+3232 is low enough to allow a small disc truncation radius, and consequently produce a double peaked X-ray pulse profile, but high enough that geometrical effects allow a detection of polarized emission. For instance, in a system with a more symmetrical field pattern than is implied here, the (positive) polarization from one pole may be effectively cancelled out by that (negative) from the other pole, so yielding a net polarization close to zero and hence undetectable. It may be that the only reason we detect polarization in 1WGA J1958.2+3232 is that the two poles have sufficiently different emission properties that a net polarization remains. In this interpretation, the question of whether each IP exhibits polarized emission is not only dependent upon the strength of the white dwarf magnetic field, but also on the geometry of the system and that of its magnetic field.

###### Acknowledgements.

The data analysis reported in this paper was carried out using facilities provided by PPARC, Starlink and the Open University Research Committee. Astronomical computing at the OU is supported under PPARC grant PPA/G/O/2000/00037; HQ was also supported by PPARC research grant GR/L64621 during this work. SK wishes to thank the Finnish Academy of Sciences and Letters (Academia Scientiarum Fennica) for support. HL wishes to acknowledge support from grants 71355 and 44011 by the Finnish Academy. The authors thank Ian Howarth for acquiring and reducing the optical spectra reported herein. The Jakobus Kapteyn Telescope and Isaac Newton Telescope are part of the Isaac Newton Group of Telescopes on the island of La Palma. They are operated on behalf of the UK Particle Physics and Astronomy Research Council and the Netherlands Organisation for Scientific Research, in the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofísica de Canarias. The Nordic Optical Telescope is operated on the island of La Palma jointly by Denmark, Finland, Iceland, Norway and Sweden, in the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofísica de Canarias. The polarization data presented here have been taken using ALFOSC, which is owned by the Instituto de Astrofísica de Andalucia (IAA) and operated at the Nordic Optical Telescope under agreement between IAA and the NBIfAFG of the Astronomical Observatory of Copenhagen. This research has made use of data obtained from the High Energy Astrophysics Science Archive Research Center (HEASARC), provided by NASA’s Goddard Space Flight Center. iraf is written and supported by the iraf programming group at the National Optical Astronomy Observatories (NOAO) in Tucson, Arizona. NOAO is operated by the Association of Universities for Research in Astronomy (AURA), Inc., under cooperative agreement with the National Science Foundation.


text:
# RUNAWAY MERGING OF BLACK HOLES: ANALYTICAL CONSTRAINT ON THE TIMESCALE

## 1 INTRODUCTION

While black holes (BHs) arising from stellar evolution have masses of orders $`10^0`$$`10^1M_{\mathrm{}}`$, those found in galaxy nuclei have masses of orders 106–10$`{}_{}{}^{9}M_{\mathrm{}}^{}`$ (Kormendy & Richstone 1995; Kormendy & Gebhardt 2001). The possible missing link, i.e., an intermediate-mass BH, was recently discovered in a starburst galaxy M82 as a source of compact X-ray emission (Kaaret et al. 2001; Matsumoto et al. 2001; see also earlier references therein). The observed strong variability implies that the source is an accreting BH. The observed luminosity of $`10^{41}`$ ergs s-1 implies that the BH mass is greater than $`10^3M_{\mathrm{}}`$ if the emission is isotropic and its luminosity is below the Eddington limit. The BH lies at the 2 µm secondary peak (Matsushita et al. 2000), which is an unusually active site of ongoing star formation. Hence the BH was formed during the ongoing starburst and is less massive than a so-called super star cluster, i.e., a building block of a starburst region. The typical mass of a super star cluster is of order $`10^6M_{\mathrm{}}`$.

The formation of an intermediate-mass BH is attributable to successive merging of stellar-mass BHs within a super star cluster (Taniguchi et al. 2000; Matsushita et al. 2000; Ebisuzaki et al. 2001; Miller & Hamilton 2001). The cluster sinks into the galaxy center due to dynamical friction and upon evaporation releases the intermediate-mass BH. If there are more than one such BHs, they merge with each other. Thus an intermediate-mass BH could eventually evolve to a supermassive BH and is crucial to studying a possible evolutionary connection between a starburst and an active galactic nuclei (Mouri & Taniguchi 2002 and references therein). The latter harbors a supermassive BH as the central engine.

Since the 2 µm secondary peak of M82 has a starburst age of order $`10^7\mathrm{yr}`$ (Satyapal et al. 1997), the formation of the intermediate-mass BH should have occurred within this short duration. There is a possibility of runaway merging, which is known for collisions of BHs and those of normal stars (Quinlan & Shapiro 1989; Lee 1993; Portegies Zwart et al. 1999). Massive objects segregate into the cluster core due to dynamical friction, increasing the merging probability. The probability itself is higher for more massive objects. Their merging products are even more massive. Thus the growth rate is higher than that expected from the merging rate at the initial time.

However, runaway merging is a nonlinear process and hence has been studied only for specific cases using numerical methods, e.g., $`N`$-body simulation. Theoretical understanding of the underlying physics is insufficient. Also, the existing numerical results are not applicable to formation of an intermediate-mass BH. A more general analytical approach, albeit less realistic, is desirable.

We make an analytical analysis of the runaway merging. In our idealized model, BHs of a single mass $`m_0`$ exist initially with a number density $`n_0`$ and a velocity dispersion $`v_0`$. These BHs are assumed to merge with each other via two-body interactions, i.e., energy loss due to gravitational radiation (§2). The BH mass distribution is studied using the Smoluchowski equation, a master equation describing the evolution of number densities of various-mass particles that merge with each other via two-body interactions. We obtain an analytical constraint on the timescale of runaway merging (§3), and apply the result to BHs in a super star cluster with $`m_0=30M_{\mathrm{}}`$, $`n_0=10^6\mathrm{pc}^3`$, and $`v_0=1\text{km\hspace{0.17em}s}^1`$ (§4).

## 2 MERGING VIA GRAVITATIONAL RADIATION

When two BHs become close, energy loss due to gravitational radiation could exceed the orbital kinetic energy, and a binary could form. The binary immediately merges through subsequent gravitational radiation. For this process, we estimate the cross section and then the rate coefficient.

We consider an interaction of two BHs where the BH masses are $`m_i`$ and $`m_j`$, the initial relative velocity is $`v_{\mathrm{rel}}`$, and the impact parameter is $`b`$. If the two BHs pass away narrowly, the orbit is close to parabolic. We thereby assume that the eccentricity $`e`$ is slightly greater than 1, and consider the leading term for $`e1`$. This is equivalent to an assumption of significant gravitational focusing. The distance of closest approach $`r_{\mathrm{min}}`$ is

$$r_{\mathrm{min}}=\frac{b^2v_{\mathrm{rel}}^2}{2G(m_i+m_j)}.$$
(1)
For an unperturbed parabolic ($`e=1`$) or hyperbolic ($`e>1`$) orbit, the quadrupole formalism gives the total energy of the gravitational radiation $`\delta E`$ at Newtonian order as

$$\delta E=\frac{8}{15}\frac{G^{7/2}}{c^5}\frac{(m_i+m_j)^{1/2}m_i^2m_j^2}{r_{\mathrm{min}}^{7/2}}g(e),$$
(2)
with $`g(e)=425\pi /(32\sqrt{2})`$ at $`e=1`$ (Turner 1977). The condition for the two BHs to pass away is

$$\delta E<\frac{1}{2}\frac{m_im_j}{m_i+m_j}v_{\mathrm{rel}}^2.$$
(3)
Equations (1)–(3) yield the minimum impact parameter $`b_{\mathrm{min}}`$ for the BHs to pass away. Then we obtain the merging cross section $`\sigma _{\mathrm{mer}}=\pi b_{\mathrm{min}}^2`$ as

$$\sigma _{\mathrm{mer}}=2\pi \left(\frac{85\pi }{6\sqrt{2}}\right)^{2/7}\frac{G^2(m_i+m_j)^{10/7}m_i^{2/7}m_j^{2/7}}{c^{10/7}v_{\mathrm{rel}}^{18/7}}.$$
(4)
Quinlan & Shapiro (1989) derived the same formula with a somewhat different reasoning. The cross section $`\sigma _{\mathrm{mer}}`$ is larger than that corresponding to the Schwarzschild radius $`\sigma _{\mathrm{Sch}}\pi (2Gm/c^2)^2(v_{\mathrm{rel}}/c)^{18/7}\sigma _{\mathrm{mer}}`$. Thus our Newtonian approximation is sufficiently accurate.

Since the merging cross section $`\sigma _{\mathrm{mer}}`$ is smaller than that for gravitational scattering, $`\sigma _{\mathrm{sca}}\pi (Gm/v_{\mathrm{rel}}^2)^2(c/v_{\mathrm{rel}})^{10/7}\sigma _{\mathrm{mer}}`$, the BH motions are assumed to be in thermal equilibrium as a zeroth-order approximation. The equipartition of the kinetic energy is assumed to be achieved as

$$\frac{1}{2}m_iv_i^2=\frac{1}{2}m_0v_0^2\text{for any }m_i.$$
(5)
Here $`v_i`$ and $`v_0`$ are three-dimensional velocity dispersions for BH masses $`m_i`$ and $`m_0`$, respectively. The probability distribution of the initial relative velocity $`v_{\mathrm{rel}}`$ is given by

$$P(x)=\left(\frac{x}{2\pi }\right)^{1/2}\mathrm{exp}\left(\frac{x}{2}\right)\text{with}x=\frac{3m_im_jv_{\mathrm{rel}}^2}{(m_i+m_j)m_0v_0^2}.$$
(6)
Equations (4) and (6) yield the average $`v_{\mathrm{rel}}\sigma _{\mathrm{mer}}`$ as the rate coefficient $`R_{i,j}`$ for merging of BHs of masses $`m_i`$ and $`m_j`$:

$$R_{i,j}=\frac{AG^2m_0^2}{c^3}\left(\frac{v_0}{c}\right)^{11/7}\left(\frac{m_im_j}{m_0^2}\right)^{15/14}\left(\frac{m_i+m_j}{m_0}\right)^{9/14},$$
(7)
with $`A=\mathrm{\Gamma }(5/7)(2\pi )^{11/14}\mathrm{\hspace{0.17em}3}^{1/2}\mathrm{\hspace{0.17em}85}^{2/7}33.33`$. However, in this estimate, mass segregation has not been taken into account. Massive BHs sink into the core of the star cluster and have a high number density even if their total number is small. The scale height of BHs of a mass $`m_i`$ is proportional to $`m_i^{1/2}`$ (multimass King’s model; Sigurdsson & Phinney 1995). Thus mass segregation enhances the merging rate by a factor $`(m_im_j/m_0^2)^{3/2}`$. The rate coefficient is then proportional to $`(m_im_j/m_0^2)^{18/7}`$ instead of $`(m_im_j/m_0^2)^{15/14}`$.

The expression (7) is too complicated to be studied analytically using the Smoluchowski equation. We use the inequality $`m_i+m_j2(m_im_j)^{1/2}`$ and obtain a simplified expression for the merging rate coefficient:

$$R_{i,j}\frac{BG^2m_0^2}{c^3}\left(\frac{v_0}{c}\right)^{11/7}\left(\frac{m_im_j}{m_0^2}\right)^{\lambda /2},$$
(8)
with $`B=\mathrm{\Gamma }(5/7)\pi ^{11/14}\mathrm{\hspace{0.17em}2}^{10/7}\mathrm{\hspace{0.17em}3}^{1/2}\mathrm{\hspace{0.17em}85}^{2/7}52.04`$. The exponent $`\lambda `$ is $`39/14`$ or $`81/14`$, respectively, if we ignore or incorporate the mass segregation.

We consider two-body interactions alone and thus ignore three-body interactions. They result in BH binaries. The BH binary undergoes subsequent interactions with other BHs and increases its binding energy. At each of the interactions, a fraction of the change in binding energy is converted into the translational kinetic energy. The BH binary eventually merges by emitting gravitational radiation or escapes from the star cluster (Sigurdsson & Hernquist 1993; see also Portegies Zwart & McMillan 2000). Although a numerical calculation based on a reliable BH mass function is required to draw a definitive conclusion, we expect that the former occurs much more often than the latter in a massive cluster with a large escape velocity if the BHs are massive enough ($`10M_{\mathrm{}}`$; Miller & Hamilton 2001). The merging timescale is short for a binary of massive BHs. The escape of a massive BH binary is due to a collision with a massive BH. Before such a rare collision, the BH binary would merge. Our present model would accordingly underestimate the merging rate. This possible underestimation is not serious because we are to obtain an upper limit on the timescale of the runaway merging.

## 3 SMOLUCHOWSKI EQUATION

Suppose that particles of a single mass $`m_0`$ exist with a number density $`n_0`$ at the initial time $`t=0`$. They merge with each other via two-body interactions. The merging of particles of masses $`m_i`$ and $`m_j`$ has a rate coefficient $`R_{i,j}`$ and results in a particle of a mass $`m_i+m_j`$. We consequently have particles of masses $`m_i=m_0i`$ ($`i=1,2,3,\mathrm{}`$) and number densities $`n_i`$. If nondimensional terms are defined as

$$\stackrel{~}{n}_i=\frac{n_i}{n_0},\stackrel{~}{R}_{i,j}=\frac{R_{i,j}}{R_{0,0}},\text{and}\stackrel{~}{t}=n_0R_{0,0}t,$$
(9)
the evolution of the number densities is described by the Smoluchowski equation:

$$\frac{d\stackrel{~}{n}_i}{d\stackrel{~}{t}}=\frac{1}{2}\underset{j+k=i}{}\stackrel{~}{R}_{j,k}\stackrel{~}{n}_j\stackrel{~}{n}_k\stackrel{~}{n}_i\underset{j=1}{\overset{\mathrm{}}{}}\stackrel{~}{R}_{i,j}\stackrel{~}{n}_j.$$
(10)
Here the first sum accounts for the increase of $`\stackrel{~}{n}_i`$ due to merging of particles satisfying the condition $`m_j+m_k=m_i`$ while the second sum accounts for the decrease of $`\stackrel{~}{n}_i`$ due to merging of particles of a mass $`m_i`$ with those of any mass. The initial condition is $`\stackrel{~}{n}_1=1`$ and $`\stackrel{~}{n}_i=0`$ for $`i2`$ at $`\stackrel{~}{t}=0`$.

The Smoluchowski equation is based on a mean-field approximation and thereby ignores spatial fluctuations, e.g., mass segregation. It is nevertheless possible to incorporate consequences of the fluctuation into the rate coefficient $`R_{i,j}`$. The number densities $`n_i`$ are then regarded as those averaged over the fluctuation scale for particles of a mass $`m_0`$.

The analytical solution of equation (10) is known only for $`\stackrel{~}{R}_{i,j}=1`$, $`(i+j)/2`$, and $`ij`$ (see Hayashi & Nakagawa 1975). We study the solution for $`\stackrel{~}{R}_{i,j}=(ij)^{\lambda /2}`$ with $`\lambda 2`$ and relate it to the solution for $`\stackrel{~}{R}_{i,j}(ij)^{\lambda /2}`$, which corresponds to our rate coefficient (8), using the $`\mathrm{}`$-th moment:

$$\stackrel{~}{M}_{\mathrm{}}=\underset{i=1}{\overset{\mathrm{}}{}}i^{\mathrm{}}\stackrel{~}{n}_i.$$
(11)
The total mass density and average mass of the particles are given respectively by $`\stackrel{~}{M}_1`$ and $`\stackrel{~}{M}_2/\stackrel{~}{M}_1`$ (Hayakawa & Hayakawa 1988). From equation (10), the moment $`\stackrel{~}{M}_{\mathrm{}}`$ evolves as

$$\frac{d\stackrel{~}{M}_{\mathrm{}}}{d\stackrel{~}{t}}=\frac{1}{2}\underset{i,j}{}\stackrel{~}{n}_i\stackrel{~}{n}_j\stackrel{~}{R}_{i,j}\left[(i+j)^{\mathrm{}}i^{\mathrm{}}j^{\mathrm{}}\right].$$
(12)
We have $`d\stackrel{~}{M}_1/d\stackrel{~}{t}=0`$, which means the conservation of the total mass, $`\stackrel{~}{M}_1=1`$. Multiplying both sides of equation (12) for $`\stackrel{~}{R}_{i,j}=(ij)^{\lambda /2}`$ and $`\mathrm{}=2`$ by $`\stackrel{~}{M}_2^\lambda `$ gives

$$\frac{d\stackrel{~}{M}_2^{(\lambda 1)}}{d\stackrel{~}{t}}(\lambda 1)\text{for }\lambda 2.$$
(13)
Here we have used the mass conservation $`\stackrel{~}{M}_1=1`$ and also a relation originating in the Schwarz inequality:

$$\frac{_{i=1}^{\mathrm{}}i^{\lambda /2+1}\stackrel{~}{n}_i}{_{i=1}^{\mathrm{}}i\stackrel{~}{n}_i}\left(\frac{_{i=1}^{\mathrm{}}i^2\stackrel{~}{n}_i}{_{i=1}^{\mathrm{}}i\stackrel{~}{n}_i}\right)^{\lambda /2}\text{for }\lambda 2.$$
(14)
The equality holds at $`\lambda =2`$. Equation (13) has the solution:

$$\stackrel{~}{M}_2\left[\frac{1}{1(\lambda 1)\stackrel{~}{t}}\right]^{\frac{1}{\lambda 1}}\text{for }\lambda 2.$$
(15)
Thus the average mass $`\stackrel{~}{M}_2/\stackrel{~}{M}_1=\stackrel{~}{M}_2`$ becomes infinity at a finite time $`\stackrel{~}{t}_{\mathrm{rm}}(\lambda 1)^1`$, owing to emergence of infinite-mass particles. We regard $`\stackrel{~}{t}_{\mathrm{rm}}`$ as the timescale of runaway merging, which is shorter than the timescale $`\stackrel{~}{t}=1`$ corresponding to the initial merging rate $`n_0R_{0,0}`$ (eq. ). Even if $`\stackrel{~}{R}_{i,j}=(ij)^{\lambda /2}`$ is replaced with $`\stackrel{~}{R}_{i,j}(ij)^{\lambda /2}`$, the solution (15) remains the same.

The total number density of the particles is given by the zeroth moment $`\stackrel{~}{M}_0`$. Hence it would appear natural to adopt $`\stackrel{~}{M}_1/\stackrel{~}{M}_0=\stackrel{~}{M}_0^1`$ as the average mass. If both sides of equation (12) for $`\stackrel{~}{R}_{i,j}(ij)^{\lambda /2}`$ with $`\lambda 2`$ and $`\mathrm{}=0`$ are multiplied by $`\stackrel{~}{M}_0^{\lambda 2}`$, it follows that $`\stackrel{~}{M}_0`$ becomes zero at a finite time $`\stackrel{~}{t}_{\mathrm{rm}}^{}2(\lambda 1)^1`$. However, before the time $`\stackrel{~}{t}_{\mathrm{rm}}^{}`$, the moment equation $`\stackrel{~}{M}_1=1`$ becomes invalid (Hendriks, Ernst, & Ziff 1983). This is the reason why we have adopted $`\stackrel{~}{M}_2/\stackrel{~}{M}_1`$ as the average mass.

The emergence of infinite-mass particles occurs in more general cases. If the rate coefficient $`\stackrel{~}{R}_{i,j}`$ scales as $`\stackrel{~}{R}_{ai,aj}=a^\lambda \stackrel{~}{R}_{i,j}=a^\lambda \stackrel{~}{R}_{j,i}`$, infinite-mass particles emerge within a finite timescale for $`\lambda >1`$ (Hayashi & Nakagawa 1975). Lee (1993) used this property to interpret runaway merging of BHs found in $`N`$-body simulations. Since the timescale of the runaway merging had not been estimated analytically, we have obtained the analytical constraint (see Hendriks et al. 1983 for similar discussion in terms of sol-gel phase transition). If the rate coefficient behaves as $`\stackrel{~}{R}_{i,j}i^\mu j^\nu `$ with $`\nu >1`$ for $`j\mathrm{}`$, infinite-mass particles emerge and high-order moments become infinity just after the initial time $`\stackrel{~}{t}=0`$ (van Dongen 1987). This is the case in our rate coefficient (8). We have ignored moments $`\stackrel{~}{M}_{\mathrm{}}`$ with $`\mathrm{}>2`$ because the Smoluchowski equation (10) is merely used as an idealized model in the present work.

## 4 Discussion

The runaway merging of BHs is caused by two mechanisms, which are equally important. First, the merging rate $`R_{i,j}`$ is intrinsically large for massive BHs (eq. ). This is because the energy loss due to gravitational radiation $`\delta E`$ is large. Massive BHs have low velocities (eq. ). If the initial relative velocity $`v_{\mathrm{rel}}`$ is low or the BH masses $`m_i`$ and $`m_j`$ are large, the gravitational focusing is significant and the distance of closest approach $`r_{\mathrm{min}}`$ is small (eq. ). If the distance $`r_{\mathrm{min}}`$ is small or the BH masses $`m_i`$ and $`m_j`$ are large, the energy loss $`\delta E`$ is large (eq. ). Second, low-velocity massive BHs sink into the cluster core and have an enhanced number density. The first and second mechanisms correspond to $`\lambda =39/14`$ and 3, respectively. They are comparable and individually exceed the critical value for the onset of runaway merging, $`\lambda =1`$ (§3).

Equations (8), (9), and (15) yield a constraint on the timescale $`t_{\mathrm{rm}}`$ of the runaway merging as

$`t_{\mathrm{rm}}`$ $``$ $`{\displaystyle \frac{1}{\lambda 1}}{\displaystyle \frac{c^3}{BG^2m_0^2n_0}}\left({\displaystyle \frac{v_0}{c}}\right)^{11/7}`$ (16)
$`=`$ $`\left({\displaystyle \frac{m_0}{30M_{\mathrm{}}}}\right)^2\left({\displaystyle \frac{n_0}{10^6\mathrm{pc}^3}}\right)^1\left({\displaystyle \frac{v_0}{1\text{km\hspace{0.17em}s}^1}}\right)^{11/7}\times \{\begin{array}{cc}4\times 10^7\mathrm{yr}\text{no mass segregation}\hfill & \\ 2\times 10^7\mathrm{yr}\text{mass segregation.}\hfill & \end{array}`$ (19)
The mass $`m_0=30M_{\mathrm{}}`$, the number density $`n_0=10^6\mathrm{pc}^3`$, and the velocity dispersion $`v_0=1\text{km\hspace{0.17em}s}^1`$ provide the most favorable conditions for the runaway merging among those expected for BHs in actual super star clusters (see below). We accordingly obtain an upper limit of order $`10^7\mathrm{yr}`$ on the timescale $`t_{\mathrm{rm}}`$. This upper limit is comparable to the starburst age observed for the 2 µm secondary peak of the starburst galaxy M82 (§1). Since a starburst duration is of order $`10^8\mathrm{yr}`$ (see de Grijs, O’Connell, & Gallagher 2001 for the case of a fossil starburst M82 B), it is generally possible that an intermediate-mass BH is formed in an ongoing starburst region.

The BH mass $`m_0=30M_{\mathrm{}}`$ has been adopted because stars heavier than 30–$`40M_{\mathrm{}}`$ evolve directly to BHs with no supernova explosion (Brown & Bethe 1994; Fryer 1999). The BH mass is equal to that of the progenitor. Although the less massive stars could evolve to BHs, they explode as supernovae and blow off most of their mass. A low-mass BH is easily ejected from the star cluster by three-body interactions and does not contribute to the intermediate-mass BH (§2).

The number density $`n_0=10^6\mathrm{pc}^3`$ and the velocity dispersion $`v_0=1\text{km\hspace{0.17em}s}^1`$ of the BHs have been adopted in the light of a nearby super star cluster, R136 in the Large Magellanic Cloud (see Portegies Zwart et al. 1999). The central mass density of R136 is $`10^6M_{\mathrm{}}\mathrm{pc}^3`$. We assume the Salpeter initial mass function $`dn/dmm^{2.35}`$ with the cutoffs at 0.1 and $`100M_{\mathrm{}}`$, which yields the number fraction 0.04% for stars heavier than $`30M_{\mathrm{}}`$ and the average mass $`0.4M_{\mathrm{}}`$. We also assume the multimass King’s model (§2). Then the above mass density corresponds to the number density $`n_010^6\mathrm{pc}^3`$ for BHs heavier than $`30M_{\mathrm{}}`$. Here wind mass loss from stars has been ignored as a zeroth-order approximation. The mass loss rate observed for a $`30M_{\mathrm{}}`$ star is $`10^6M_{\mathrm{}}\mathrm{yr}^1`$ (Becker 1999). Since the lifetime of a $`30M_{\mathrm{}}`$ star is $`10^6`$ yr, the total mass loss is at most a few times $`M_{\mathrm{}}`$. On the other hand, the stellar velocity dispersion of R136 is estimated to be about $`10\text{km\hspace{0.17em}s}^1`$, which corresponds to the velocity dispersion $`v_01\text{km\hspace{0.17em}s}^1`$ for BHs of a mass $`30M_{\mathrm{}}`$. If runaway merging occurs among stars before they evolve to BHs (§1), the number density $`n_0`$ could be higher and hence the upper limit on the timescale $`t_{\mathrm{rm}}`$ could be lower.

The timescale $`t_{\mathrm{ms}}`$ of mass segregation is less than that of the runaway merging. If the average stellar mass $`m`$ is $`0.4M_{\mathrm{}}`$, the mass density $`\rho `$ is $`10^6M_{\mathrm{}}\mathrm{pc}^3`$, and the velocity dispersion $`v`$ is $`10\text{km\hspace{0.17em}s}^1`$, the star cluster has the relaxation timescale $`t_{\mathrm{rlx}}v^3/(\pi G^2m\rho )=4\times 10^7\mathrm{yr}`$, which yields $`t_{\mathrm{ms}}(m/m_0)t_{\mathrm{rlx}}=5\times 10^5\mathrm{yr}`$ for BHs of a mass $`m_0=30M_{\mathrm{}}`$ (Spitzer 1969).

Super star clusters have ranges of structure parameters. Existing few observations of super star clusters other than R136 tend to imply $`v_0>1\text{km\hspace{0.17em}s}^1`$ (Ho & Filippenko 1996; Smith & Gallagher 2001). Observations of globular clusters, i.e., likely descendants of super star clusters, tend to imply $`n_0<10^6\mathrm{pc}^3`$ (Pryor & Meylan 1993). In these cases, the timescale $`t_{\mathrm{rm}}`$ of the runaway merging is long (eq. ). The formation of an intermediate-mass BH in an ongoing starburst region occurs only under rare favorable conditions, i.e., high number density $`n_0`$ and low velocity dispersion $`v_0`$ of stellar-mass BHs. Also, the star cluster as a whole has to be sufficiently massive in order to provide a sufficient number of the stellar-mass BHs. In fact, only one intermediate-mass BH heavier than $`10^3M_{\mathrm{}}`$ has been found among more than 100 super star clusters of the starburst galaxy M82, although there might exist some intermediate-mass BHs that are not accreting the gas and hence are not observable.


text:
# Formation of Turbulent Cones in Accretion Disk Outflows and Application to Broad Line Regions of Active Galactic Nuclei

## 1 INTRODUCTION

Non-spherical outflows are a ubiquitous feature of AGN. In this Letter we discuss the stability of such an outflow launched normally from the disk with a decaying power law velocity profile. The motivation for this work is to understand the conditions under which the flow may be unstable to cloud formation, and to infer the structure of the region in which such clouds may reside. Elvis (2000) has argued that data require AGN BLR clouds to reside in a narrow biconical structure. Our work herein supports this possibilty.

A variety of launching mechanisms for AGN accretion disk outflows have been studied: radiatively accelerated outflows (Arav et al., 1994), hydrodynamic line-driven winds (Proga et al., 1999, 2000), hydromagnetic disk winds (Königl & Kartje, 1994), (Pelletier & Pudritz, 1992), thermal wind-type outflows in low luminosity AGNs (Pietrini & Torricelli-Ciamponi, 2000), and others. In this work we parameterize the disk wind independent of the launching mechanism.

A disk wind outflow can become unstable from linear and nonlinear perturbations in velocity, density, temperature, ionization state, etc. These can lead to cloud formation. The cloud model is a leading, though not unanimously accepted, paradigm for the structure of the AGN BLR regions: models exist with (Elvis, 2000; Urry & Padovani, 1995) and without (Murray & Chiang, 1995; Murray et al., 1995) clouds. The main reason for the lack of agreement is the uncertainty about cloud survival in AGN enviroments (Mathews & Ferland, 1987). They may be unstable to evaporation (e.g., see (Pier & Voit, 1995)), or be insufficiently supported by pressure and be dynamically unstable to shredding (Mathews & Ferland, 1987; Poludnenko et al., 2002). Among the suggested solutions has been magnetic confinement of clouds (Rees, 1987; Bottorff & Ferland, 2000).

However a key point is that a typical cloud needs only to exist long enough to reprocess radiation (Rees, 1987; Celotti et al., 1992; Kuncic et al., 1996, 1997). If clouds are destroyed thereafter and new clouds are formed, then the system can exist in dynamical equilibrium even though each cloud loses its identity swiftly. This may apply to the broad line region. Short cloud survival and regeneration times can be associated with turbulence (Bottorff & Ferland, 2001; Bottorff et al., 2000). The balance between formation and destruction maintains a constant BLR cloud number density.

Our analysis suggests that a nonlinear wind instability resulting from the combination of launching profile + radiation pressure can lead to the formation of a turbulent biconical zone that might manifest itself observationally as the BLR. This biconical structure is the region where flow streamlines cross. Such a bicone is remarkably similar to that inferred empirically by Elvis (2000) and from numerical simulations by Proga et al. (2000). Growing observational data corroborates the existance of such a structure (e.g. NGC 1068 (Arribas et al., 1996), Mrk 3 (Ruiz et al., 2001), etc.)

We formulate the problem in section 2.1 and describe the flow field in section 2.2. We analyze the linear and nonlinear stability of the flow and a possible scenario for thermal instability resulting in formation of the BLR clouds in section 2.3. In section 3 we discuss the results and their implications for the AGN structure and dynamics.

## 2 DISK WIND MODEL

### 2.1 Formulation of the problem

All of our calculations are performed in cylindrically symmetric geometry, however we do not include axial rotation. At the origin we assume a black hole of mass $`M_{BH}`$. Fig. 1 shows the setup.

We take a Keplerian, geometrically thin, optically thick accretion disk to be located along the abscissae with the distance $`\lambda `$ measured along the disk. The central engine is taken to be a point-like source at the origin coincident with the black hole. We exclude disk surface illumination by the central source including all disk properties into the disk-wind launch velocity profile.

Our wind originates at the disk surface with launch velocity profile

$$v_{z0}(\lambda )=v^{}\left(\frac{\lambda }{\lambda ^{}}\right)^n,$$
(1)
where $`v^{}`$ is the wind velocity at the innermost launch point of the disk with coordinate $`\lambda ^{}`$ and $`n`$ is an arbitrary non-positive number111We discard positive values of $`n`$ as unphysical.. We assume that the wind propagates into a low density medium.

After the wind material leaves the disk, it is subjected to gravity from the black hole and radiation from the central engine. Proga et al. (1999, 2000) discuss one possible approach to line driven winds in the context of numerical modeling. To capture the most important qualitative features of the flow and to simplify the analysis we assume a continuum driven wind and use only Thomson opacities. This assumption somewhat underestimates the driving force since it excludes the contribution of line-driving. However, on the other hand, it is an overestimate since typically a significant fraction of the central source spectrum lies in the Compton scattering regime. To compensate for that we introduce a factor $`f_T<1`$, intended to decrease the total luminosity $`L`$. The choice of a specific value of $`f_T`$ depends on a particular situation being considered.

In our model, a fluid parcel of density $`\rho `$ leaves the disk with the velocity of equation (1), determined by the coordinate of its launch point, and finds itself in a centrally symmetric, conservative potential field. For Thomson scattering the radiative force per unit volume is

$$\stackrel{}{F}^{rad}(r)=\frac{1}{c}\frac{\stackrel{}{r}}{|\stackrel{}{r}|}\alpha _\nu F_\nu 𝑑\nu =\frac{n\sigma _{T,e}}{c}F\frac{\stackrel{}{r}}{|\stackrel{}{r}|},$$
(2)
and the effective potential field is

$$U(r)=U_{rad}+U_{grav}=\frac{\rho GM_{BH}}{r}\left(\frac{f_TL}{L_{Edd}}1\right)=\rho \frac{\alpha }{r},$$
(3)
where $`L_{Edd}=1.2510^{38}erg/s(M_{BH}/\mathrm{M}_{})`$ is the Eddington luminosity. Note, that here we give an expression for the potential energy per unit volume.

### 2.2 Description of the flow field

Finding a trajectory and velocity of a fluid parcel in the outflow reduces to the problem of the motion of a particle in a centrally symmetric repulsive potential, described by (3). The trajectory is a hyperbola (Landau & Lifshitz, 1976), and in polar coordinates is

$$r(\varphi )=\frac{p}{e\mathrm{cos}\varphi 1}.$$
(4)
Here $`\varphi [0,\pi /2]`$. Eccentricity of the orbit $`e`$ and parameter of the orbit $`p`$ are defined as

$$e=1+\frac{p}{r_0},p=2K\left(\frac{r_0}{r^{}}\right)^{2n+1}r_0,$$
(5)
where $`r_0=\lambda _0`$ is the coordinate of the fluid element launch point. The dimensionless quantity

$$K=\lambda ^{}(v^{})^2/2\alpha ,$$
(6)
is the ratio of the initial kinetic energy of a fluid particle launched at $`\lambda ^{}`$ to its potential energy and measures the relative importance of the disk wind vs. central radiation source. ($`K`$ and $`n`$ are the key parameters in our analysis). Equation (4) can be considered as an implicit solution for $`\lambda _0`$ as a function of $`z`$ and $`\lambda `$. For any point $`(\lambda ,z)`$ we can uniquely determine the launch point of the streamline, $`\lambda _0`$. Fig. 1 shows a sample streamline pattern.

The velocity distribution along a streamline is

$$\{\begin{array}{ccc}v_\lambda (\lambda _0,\lambda )\hfill & =\hfill & \sqrt{\frac{\alpha }{\lambda _0}(e+1)}\frac{\mathrm{sinh}\xi }{e\mathrm{cosh}\xi +1},\hfill \\ v_z(\lambda _0,\lambda )\hfill & =\hfill & \sqrt{\frac{\alpha }{\lambda _0}(e^21)(e+1)}\frac{\mathrm{cosh}\xi }{e\mathrm{cosh}\xi +1}.\hfill \end{array}$$
(7)
Parameter $`\xi [0,+\mathrm{})`$ is defined by $`\mathrm{cosh}\xi =(\lambda /\lambda _0)(e+1)e`$.

### 2.3 Stabilty analysis

We base our linear stability analysis for the inviscid case on the *Fjørtoft theorem*, which is derived from the inviscid limit of the Orr-Sommerfeld equation, i.e. Rayleigh equation (Yih, 1969; Maslowe, 1985). We have studied the outflows with the values of $`2.5n0.0`$ and $`0.05K2.5`$. In all the regimes studied, the disk wind is linearly stable to infinitesimal velocity perturbations.

However a study of equation (4) shows that a non-linear instability can be present in a subregion of the flow. Depending on $`n`$, a streamline inclination angle at any given point $`(r,\lambda )`$ is either a monotonically increasing or monotonically decreasing function of $`\lambda _0`$. The condition for the change in the character of monotonicity can be obtained: using equation (4) we can determine the inclination angle of a given streamline $`\mathrm{tan}\mathrm{\Theta }(\lambda _0,\lambda )=dz/d\lambda .`$ Using

$$\frac{d}{d\lambda _0}\left(\underset{\lambda \mathrm{}}{lim}\mathrm{tan}\mathrm{\Theta }(\lambda _0,\lambda )\right)=0,$$
(8)
we find that the monotonicity changes for $`n=1/2`$, i.e. the Keplerian disk is the borderline case. This means that for $`n>1/2`$, each successive streamline is less inclined towards the abscissae axis than the preceding one and eventually intersects all of the preceding streamlines. For $`n1/2`$ the inclination angle is a monotonically non-increasing function of $`\lambda _0`$ and then the streamlines never intersect.

The intersection of two streamlines can be approximated as an intersection of two jets of vanishing thickness whose interaction leads to a resultant turbulent jet with an inclination angle lying between that of the two original jets (e.g., see (Landau & Lifshitz, 1959)). Therefore, streamline intersection causes the biconical zone of compressible turbulence in the overall disk-wind outflow. The extent of such a zone can be estimated: since the streamline inclination angle is a monotonically increasing function of $`\lambda _0`$, the streamline originating at the innermost radius of the disk is the most inclined one. Therefore, it delineates the lower boundary of the turbulent biconical zone. Its inclination angle is

$$\mathrm{tan}\mathrm{\Theta }_{min}=\underset{\lambda \mathrm{}}{lim}\frac{dz}{d\lambda }(\lambda ^{},\lambda )=2\sqrt{K+K^2}.$$
(9)
An exact shape of the upper boundary of the bicone can be found by solving $`dz/d\lambda _0=0`$ for $`\lambda _0`$ as a function of $`\lambda `$ and then substituting it into the equation (4) to find the exact dependence $`z(\lambda )`$ for the upper boundary. Instead, we give an approximate expression for the inclination angle of the outer boundary, found by analogy to equation (9):

$$\begin{array}{ccc}\mathrm{tan}\mathrm{\Theta }_{max}\hfill & =\hfill & lim_\lambda \mathrm{}\frac{dz}{d\lambda }(\lambda _0,\lambda )\hfill \\ & =\hfill & 2\sqrt{K\left(\frac{\lambda _0}{\lambda ^{}}\right)^{2n+1}+K^2\left(\frac{\lambda _0}{\lambda ^{}}\right)^{2(2n+1)}}.\hfill \end{array}$$
(10)
This defines the basic geometry of the bicone: the *inclination angle* $`\mathrm{\Theta }=(\mathrm{\Theta }_{max}+\mathrm{\Theta }_{min})/2`$ and the *divergence angle* $`\mathrm{\Delta }\mathrm{\Theta }=\mathrm{\Theta }_{max}\mathrm{\Theta }_{min}`$. See Figure 1.

Typically, the bicone has a divergence angle $`\mathrm{\Delta }\mathrm{\Theta }`$ that slowly increases with radius. Therefore equation (10), which gives the limiting inclination angle of the streamline originating from the point $`\lambda _0`$, depends on the choice of such a particular streamline. However, the increase in the divergence angle is typically small. Moreover, a particular streamline used to determine $`\mathrm{\Theta }_{max}`$ can be determined independently from the properties of the accretion disk or the estimated extent of the outflow.

Three other effects will enter in a more detailed hydrodynamic treatment that influence the divergence angle of the bicone. An effect which tends to make the divergence angle smaller than that determined using equations (9) and (10) arises because the resultant of two intersecting jets lies at an intermediate angle between the original two. Each successive interaction of resultant jets would then tend to narrow the overall conical structure. In addition, the change in the optical depth of the outflow after the formation of the turbulent bicone would shade the outer parts of the outflow. This will make the bicone narrower by increasing the inclination angle of the outer streamlines due to the reduced radiation pressure. In competition with both of these above effects is the turbulence itself, which could broaden the bicone.

The initial density fluctuations in the turbulent biconical zone cannot themselves represent the BELR clouds; the density and temperature contrasts are too small compared to the ambient warm highly ionized gas (WHIM in the nomenclature of Elvis (2000)). However, thermal instability in such a medium can transform initial turbulent overdensities into the BELR clouds (see (Burkert & Lin, 2000; Hennebelle & Pérault, 1999) for the discussion of two-phase medium formation via linear thermal instability in optically thin regions). The formation of the compressible turbulent bicone from intersecting flows can produce inhomogeneities with density contrast up to 4 in the adiabatic case and higher if one allows for the possibility of radiative cooling. Such high density contrasts (compared to the linear density fluctuations discussed by Burkert & Lin (2000)) can serve as a seed for thermal instability. The instabilty proceeds in the nonlinear regime from the onset. Rapid cooling would then produce clumps. The maximum clump density reached at a given scale is determined by how much cooling can occur during an eddy turnover time on that scale, i.e. during an average clump survival time.

This mechanism can lead to a dynamical equilibrium described by the balance between in situ formation, destruction, and re-formation of clumps comoving with the ambient flow. As long as any given clump survives long enough to reprocess radiation, problems of clump confinement against thermal conduction and destruction via high velocity differences between the clumps and the ambient medium are eliminated. The dynamical equilibrium of the turbulent flow, ensures a constant number density of BELR clouds in the turbulent biconical zone over observation durations.

Finally, there is a site where a linear instability of the Kelvin-Helmholtz type might operate, namely the shear layer between the disk wind outflow and the infalling (or stationary) circumpolar material, that fills the cone of the outflow. Simulations of Proga et al. (2000) show the formation of clumpy structure in that boundary shear layer. They attribute it to the Kelvin-Helmholtz instability. However, in regimes for which the outflow develops the turbulent bicone, such linear instability will be completely overwhelmed by the nonlinear instability induced by the intersection of streamlines. In cases when the bicone is not formed, the outflow density decreases from the maximum values near the plane of the disk to nearly the values of the circumpolar infalling material in the upper parts of the outflow with the density transition occurring over the large range of angles. BLR clouds are not expected in this case. 222Note, that even in the cases when the turbulent bicone is absent the upper boundary of the outflow is still roughly defined by the streamline launched at $`\lambda ^{}`$..

## 3 DISCUSSION

We considered the linear and nonlinear stability of a disk outflow and found that it is linearly stable to infinitesimal velocity perturbations. However, when $`n>1/2`$ and a supereddington radiation source is present, nonlinear instability leads to the formation of a biconical zone of compressible turbulence. High density contrasts in this zone may trigger thermal instability and lead to the further condensation of clumps into BELR clouds. Turbulence plays a key role in establishing the dynamical equilibrium which maintains a steady cloud number density, even though each cloud is short lived.

In the regimes when the outflow can develop the turbulent bicone, the geometry is remarkably similar to that described by Elvis (2000). For a disk outflow with $`K0.05`$ and an initial launch velocity profile slightly flatter than Keplerian, namely $`n=0.45`$, the turbulent bicone inclination angle is $`\mathrm{\Theta }28^o`$ with the divergence angle $`\mathrm{\Delta }\mathrm{\Theta }7^o`$. This is in quantitative and qualitative agreement with Elvis (2000): (1) no absorbers are along the line of sight passing above the outflow; (2) along the line of sight passing through the turbulent bicone broad absorption and emission line features will be observed; (3) narrow absorption lines (NAL) will be seen along sight lines that fall inside of the bicone. Figure 1 shows the geomtry of the outflow and the emerging bicone with the angle values found here.

The outflow geometry in our picture depends only on dimensionless quantities $`K`$ and $`n`$. However, in addition to $`\mathrm{\Theta }`$ and $`\mathrm{\Delta }\mathrm{\Theta }`$ determined above, we can also infer the inner scale by matching observed cloud velocities. We assume the same value of $`n`$ and $`K`$ as in the previous paragraph. We take $`M_{BH}=10^9\mathrm{M}_{}`$ and broad absorption line velocity $`v_{BAL}=10^4`$ km $`\mathrm{s}^1`$(the velocity observed when the line of sight passes directly along the bicone towards the central source). In our model $`v_{BAL}`$ is the terminal velocity at infinity of material launched with the largest initial velocity and potential, i.e. at the point $`\lambda _0=\lambda ^{}`$. Setting $`\lambda _0=\lambda ^{}`$ in (7) and using the expressions for $`v_\lambda (\lambda ^{},\lambda )`$ and $`v_z(\lambda ^{},\lambda )`$ in the limit $`\lambda \mathrm{}`$, we find

$$\lambda ^{}=\frac{2\alpha (1+K)}{v_{BAL}^2}.$$
(11)
Assuming $`L=1.5L_{Edd}`$ and $`f_T=0.7`$ the corresponding value of $`\lambda ^{}=\lambda _01000a.u.`$, which matches Elvis (2000). The corresponding maximum launch velocity is $`2000`$ km $`\mathrm{s}^1`$. Therefore, the observed velocity in the NAL region is $`1000`$ km $`\mathrm{s}^1`$.

In principle, the paradigm presented here applies not only to luminous AGN but to any source with a centrally symmetric potential, a sufficiently luminous central radiation source, and a disk wind, e.g. disk winds in young stellar objects.

We have not discussed the role of magnetic fields, disk rotation, metallicity, wind density fall off, or the physics of line driving. These should be considered in future work.

The authors thank V. Pariev for extremely important comments and suggestions. AYP thanks N. Murray and M. Elvis for valuable discussions. AYP and AF acknowledge support from NSF grant AST-9702484 and NASA grant NAG5-8428. EGB acknowledges support from DOE grant DE-FG02-00ER5460. The authors acknowledge DOE support from the Laboratory for Laser Energetics.