A Statistical Method to Search for Recoiling Supermassive Black Holes in Active Galactic Nuclei
aa r X i v : . [ a s t r o - ph . GA ] O c t Mon. Not. R. Astron. Soc. , 1–9 (2015) Printed 9 July 2018 (MN L A TEX style file v2.2)
A statistical method to search for recoiling supermassiveblack holes in active galactic nuclei
P. Raffai , ⋆ , Z. Haiman , and Z. Frei , Institute of Physics, E¨otv¨os University, 1117 Budapest, Hungary MTA-ELTE EIRSA “Lend¨ulet” Astrophysics Research Group, 1117 Budapest, Hungary Department of Astronomy, Columbia University, New York, NY 10027, USA
Released 2015 Xxxxx XX
ABSTRACT
We propose an observational test for gravitationally recoiling supermassive black holes(BHs) in active galactic nuclei, based on a correlation between the velocities of BHsrelative to their host galaxies, | ∆ v | , and their obscuring dust column densities, Σ dust (both measured along the line of sight). We use toy models for the distribution ofrecoil velocities, BH trajectories, and the geometry of obscuring dust tori in galacticcentres, to simulate 2 . × random observations of recoiling quasars. BHs with recoilvelocities comparable to the escape velocity from the galactic centre remain bound tothe nucleus, and do not fully settle back to the centre of the torus due to dynamicalfriction in a typical quasar lifetime. We find that | ∆ v | and Σ dust for these BHs arepositively correlated. For obscured (Σ dust >
0) and for partially obscured (0 < Σ dust ∼ < . / m ) quasars with | ∆ v | >
45 km / s, the sample correlation coefficient betweenlog ( | ∆ v | ) and Σ dust is r = 0 . ± .
02 and r = 0 . ± .
02, respectively. Allowingfor random ± / s errors in | ∆ v | unrelated to the recoil dilutes the correlationfor the partially obscured quasars to r = 0 . ± .
004 measured between | ∆ v | andΣ dust . A random sample of ∼ > ,
500 obscured quasars with | ∆ v | >
45 km / s wouldallow rejection of the no-correlation hypothesis with 3 σ significance 95 per cent ofthe time. Finally, we find that the fraction of obscured quasars, F obs ( | ∆ v | ), decreaseswith | ∆ v | from F obs ( <
10 km / s) & . F obs (cid:0) > km / s (cid:1) . .
4. This predictedtrend can be compared to the observed fraction of type II quasars, and can furthertest combinations of recoil, trajectory, and dust torus models.
Key words: black hole physics — methods: observational — galaxies: active —galaxies: nuclei.
Both theoretical models (Begelman, Blandford & Rees1980; Volonteri, Haardt & Madau 2003) and observations(e.g. Comerford et al. 2009) suggest that it is common forsupermassive black holes (SMBHs) in galactic centres toform binaries that gradually lose their energy through radi-ating gravitational waves (e.g. Haehnelt 1994; Sesana et al.2005). Numerical simulations of black hole (BH) binarymergers (see e.g. Healy, Lousto & Zlochower 2014 and ref-erences therein) suggest that the merger remnant SMBHsreceive a recoil velocity of typically several hundred (andin some spin and mass-ratio configurations up to thou-sands) of km/s, due to highly anisotropic gravitational-wave emission in the final merger phase. Taking into ac-count the spatial distribution of mass in the galactic centre ⋆ E-mail: praff[email protected] (PR); region, and corresponding dynamical friction, this meansthat the SMBH can engage in a damped oscillating mo-tion (Madau & Quataert 2004; Komossa & Merritt 2008;Tanaka & Haiman 2009) with an amplitude comparable to,or exceeding the O (10-100 pc) size of the optically thickdusty molecular torus (‘dust torus’) believed to be surround-ing galactic centres (for an overview, see H¨onig 2008). Fora large recoil, the SMBH can escape from the galactic cen-tre and can remain wandering within the dark matter halo(Guedes et al. 2009, 2011). As accreting material can staygravitationally bound to the moving SMBH, and thus theSMBH can remain active for 10 − yr after the mergerevent (e.g. Loeb 2007), kinematic and spatial signatures ofthe recoil could be found in the spectra of quasars (QSOs; seee.g. Bonning, Shields & Salviander 2007; Blecha et al. 2011;Komossa 2012; Blecha et al. 2015). These signatures couldconfirm the gravitational recoil of merger remnant SMBHs,and could probe the SMBH binary parameters, recoil tra- c (cid:13) Raffai, Haiman, & Frei jectory models, as well as the spatial distribution and com-position of matter in the galactic centre region.In this paper, we propose that in the presence ofdust tori obscuring galactic nuclei, the gravitational re-coil of active merger remnant SMBHs should introduce acorrelation between dust column mass densities along theline of sight (Σ dust ) and magnitudes of line-of-sight pe-culiar velocities of SMBHs relative to their host galax-ies ( | ∆ v | ). Proxies to estimate both of these quantitiescan be measured from observable features of QSO spectra(Bonning, Shields & Salviander 2007; Ledoux et al. 2015).As pointed out earlier by Komossa & Merritt (2008), re-coiling BHs can spend a significant fraction of their timeoff-nucleus, possibly reducing the fraction of fully obscured(‘type II’) QSOs; we also follow up on this suggestion.Using a selected combination of models of gravita-tional recoil, SMBH trajectories, and obscuring dust tori,we demonstrate the feasibility of detecting Σ dust – | ∆ v | cor-relations by simulating a set of 2 . × random observationsof recoiled QSOs with a Monte Carlo method. We charac-terize the strength of correlation between Σ dust and | ∆ v | forobscured (i.e. Σ dust >
0) QSOs in two different | ∆ v | inter-vals, and estimate the number of obscured QSOs that couldbe used to reject the hypothesis of no correlation betweenΣ dust and | ∆ v | with 3 σ significance. As the strength of thecorrelation between Σ dust and | ∆ v | , as well as the underly-ing Σ dust ( | ∆ v | ) relation depends on the presumed combina-tion of models, observational studies on QSO spectra couldprovide an opportunity for testing chosen combinations ofthese models. We also calculate the fraction of QSOs ob-scured by their dust tori, F obs ( | ∆ v | ), which, compared tothe observed fraction of type II (i.e. obscured) QSOs couldprovide an independent test for a chosen combination of re-coil, SMBH trajectory, and dust tori models. We propose touse the SDSS-DR10 Quasar Catalog (Pˆaris et al. 2014) toperform these tests in the near future, and will report onthis observational search in a separate publication.The paper is organized as follows. In §
2, we describe thegravitational recoil and SMBH trajectory models we used inour Monte Carlo simulations. In §
3, we describe our imple-mentation of a smooth dust torus model. We present ourresults in § §
5. Finally, we offer our conclusions and summarize the im-plications of this work in § Anisotropic emission of gravitational waves emitted by co-alescing SMBHs carry away linear momentum, resulting ina recoil of the merger remnant SMBH in the opposite di-rection. Numerical simulations of this process have beencarried out for merging SMBHs with equal and unequalmasses, zero and non-zero spins aligned or counter-alignedwith the orbital angular momentum, and with spins point-ing in random directions with equal probability (see e.g.Healy, Lousto & Zlochower 2014 and references therein).To keep our Monte Carlo simulation as general andrealistic, but at the same time, as computationally cheapas possible, following Tanaka & Haiman (2009), we adoptedthe analytical formulae given in equation (4) of Baker et al. (2008) to construct the distribution of recoil velocity mag-nitudes, v recoil (note that Baker et al. 2008 gives resultsvery similar to more recent, slightly modified formulae inLousto et al. 2012 and Healy, Lousto & Zlochower 2014).The direction of the recoil velocity depends on the ori-entation of the spin and orbital angular momentum vectors.It is plausible that prior to merger, the BHs accrete a signifi-cant amount of gas coherently from a circumbinary accretiondisc, whose inner regions are aligned with the binary’s or-bital plane (e.g Ivanov, Papaloizou & Polnarev (1999)). Onethen expects that the spin angular momentum vectors at thetime of the merger may be aligned with the orbital angularmomentum of the binary (Bogdanovi´c, Reynolds & Miller2007; although there could still be significant misalign-ment at merger for fast-spinning and unequal-mass bina-ries; Gerosa et al. 2015). Thus, the kick direction may notbe random, and may lie preferentially in the plane of thecircumbinary disc (and also in the binary orbital plane).However, we emphasize that here we are interested only inthe kick direction relative to the orientation of the largerscale nuclear torus. It is much less clear whether the bi-nary orbital plane and/or accretion disc is aligned with thesymmetry plane of this torus. This depends on the trans-fer of angular momentum between large (100 pc) and small(sub-parsec) scales, which is sensitive to turbulence, star for-mation, and feedback processes (Dubois et al. 2014). Obser-vationally, parsec-scale maser discs (which could be takenas proxies for the accretion discs that determine the kickdirection) appear to be oriented randomly with respect tothe plane of their host galaxies (Kormendy & Ho 2013). Wetherefore simply assume that the kick direction is random,i.e. not aligned with the symmetry plane of the torus. Thus,we chose directions of recoil velocities randomly from a uni-form distribution covering a whole sphere.Using the fitting formulae of Baker et al. (2008), wecalculated v recoil for a given pair of masses ( m , ) and di-mensionless spin vectors of the merging SMBHs ( ~α , ≡ c~S , /Gm , , where ~S , are the spins of the SMBHs, c isspeed of light, and G is Newton’s constant). Observationsof active galactic nuclei (AGNs) indicate that a significantnumber of SMBHs have | ~α , | > .
9, although a second pop-ulation of SMBHs with 0 . < | ~α , | < . m , > × M ⊙ SMBHs (see Reynolds 2013). In our simu-lations, for simplicity, we first assigned each merging SMBHto one of the two | ~α , | populations with equal probability,and then drew a random | ~α , | value from the corresponding(i.e. | ~α , | ∈ [0 . ,
1] or | ~α , | ∈ [0 . , . m and m . The fact that we chose the massof the merger remnant as M = m + m means that in ourMonte Carlo simulations M values are under-represented atthe lowest ( ∼ M ⊙ ) and at the highest ( ∼ M ⊙ ) val-ues, and overrepresented at values in between, compared tothe SMBH mass function given in Aller & Richstone (2002).Fig. 1 shows the histograms M and v recoil values we obtainedby randomizing 2 . × pairs of merging SMBHs.Using the above set of v recoil ’s with uniformly dis-tributed random directions, we simulated the resulting c (cid:13) , 1–9 method to search for recoiling SMBHs in AGNs (M [M Sun ]) SMBH mass functionSimulated M values (cid:1) (cid:0) (v recoil [km/s]) Figure 1.
Left: histogram of 2 . × simulated merger rem-nant SMBH masses M = m + m (red curves), where the valuesof 2 m and 2 m were randomized from the SMBH mass func-tion given in Aller & Richstone (2002) and represented by darkblue curves. As a result of the randomization process, M val-ues are under-represented at the lowest ( ∼ M ⊙ ) and at thehighest ( ∼ M ⊙ ) values, and overrepresented at values in be-tween, compared to the SMBH mass function. Right: histogramof the corresponding recoil velocity magnitudes, v recoil , of themerger remnant SMBHs, obtained using the framework presentedin Baker et al. (2008), and assuming that merging SMBHs havedimensionless spin magnitudes, | ~α , | , from two distinct popu-lations (i.e. | ~α , | ∈ [0 . ,
1] and | ~α , | ∈ [0 . , . | ~α , | values are uniformly dis-tributed within each of the two | ~α , | intervals, and that ~α , have uniformly distributed random directions. SMBH trajectories based on the model presented inMadau & Quataert (2004). The initial position of the SMBHand the origin of our coordinate system were chosen to co-incide with the galactic centre. This model assumes thatthe SMBH is embedded in a spherical stellar bulge, andis decelerated by dynamical friction. The one-dimensionalvelocity dispersion of stars in the bulge were calculatedusing the empirical M – σ relation (Tremaine et al. 2002)as σ D = (cid:0) / √ (cid:1) × (cid:0) M/ . × M ⊙ (cid:1) /
200 km / s. Ac-cording to observations by Barth, Greene & Ho (2005) andGreene & Ho (2006), this relation extends to active SMBHswith masses as low as M ∼ M ⊙ . We simulated the trajec-tories with a ∆ t = 10 yr time resolution up to a randomlypicked final time between T ∈ [∆ t, × ∆ t ]. The maximumduration of T max ≡ × ∆ t = 3 × yr was chosen tobe comparable with the maximum known duration of QSOactivity (see Martini 2004 for a review). The simulation wasterminated before T if the orbit of the SMBH has decayedand the SMBH settled back at the galactic centre (i.e. if theradial distance, d , and radial velocity, v r , of the SMBH con-verged to d < − pc and v r < − km / s, respectively).The position, velocity, and ∆ v of the SMBH in its final statewere calculated and stored, as if they were results of a QSOobservation made at a random time during the QSO activ-ity after the recoil. Histograms of d and v r values obtainedfor the 2 . × QSOs are shown in Fig. 2. Note that weassumed that recoiled SMBHs remain active along their en-tire trajectory, and also that they are active at the time themock observation is made.Guedes et al. (2009) followed the motions of mergerremnant SMBHs with fixed masses in various configurationsof dark matter haloes using numerical simulations. Theyconcluded that due to the asymmetric and inhomogeneousmass distribution in the halo, the SMBH trajectories suffer − − (VAL) VAL = d [pc]VAL = v r [km/s] Figure 2.
Histogram of distance d (in pc; blue curve), and ra-dial velocity magnitude, v r (in km/s; red curve), of the 2 . × SMBHs, both measured from the galactic centre at the randomlypicked time of a mock observation. Both histograms show bi-modality. The peaks at lower d and v r values (log (VAL) < − (VAL) > d and v r values correspond to SMBHs whose initial recoil veloc-ity was high enough to make them gravitationally unbound andescape from the galactic centre region. large deviations from the Madau & Quataert (2004) trajec-tories. SMBHs might not even return to the galactic centre,after the SMBHs reach the first turning point in their os-cillating motion. To check how much this could affect theresults of our Monte Carlo simulations, we simulated an in-dependent set of 10 QSOs, and examined the fraction ofSMBHs that reach the first turning point before the mockobservation is made. We found that for the entire QSO sam-ple, this fraction is ≃
31 per cent, while it is reduced to ≃
11 and ≃ | ∆ v | > / s and | ∆ v | >
45 km / s, respectively. The reason for the full pop-ulation having a higher fraction is that once the SMBHsmake their first U–turn, their orbits decay rapidly due todynamical friction. As a result, many of these end up with | ∆ v | < / s and | ∆ v | <
45 km / s at the mock observationtime (as shown in Fig. 2). As we will show in § | ∆ v | > / s show a Σ dust – | ∆ v | correlation.Furthermore, below we propose to study only the subset ofQSOs with | ∆ v | >
45 km / s, due to the limitations imposedby inevitable random velocity errors. We conclude that inthis sub-sample of QSOs, the effects of asphericity and in-homogeneities will reduce the predicted correlations only by ≈ ∼
10 timeslonger for the oscillations to be damped. Here, we neglectthis effect, leaving its evaluation to future work. However, wenote that it would likely result in a larger fraction of SMBHsremaining displaced from the nucleus, and possibly a largerfraction of QSOs exhibiting the correlations we propose. c (cid:13) , 1–9 Raffai, Haiman, & Frei
According to the unified scheme of AGNs, galactic nuclei areobscured by optically and geometrically thick dusty molec-ular tori, with the amount of obscuration depending on theviewing angle (Antonucci 1993). Even though dust typicallyonly constitutes ∼ ρ d (see their equation 8). The effective potential createdby the central SMBH and the angular momentum distribu-tion in the stellar core makes ρ d axisymmetric around thegalactic centre.To cover the range of masses M ∗ of the nuclear starclusters in the examples in Schartmann et al. (2005) and theobservations presented in Leigh, B¨oker & Knigge (2012), wechose M ∗ = 10 β M ⊙ , where β was drawn from a uniformdistribution in the range β ∈ [0 , M dust = (5 . × M ⊙ ) × [ M ∗ / (2 × M ⊙ )] in order to reproduce parameter valuesof the example torus given in table 1 in Schartmann et al.(2005).The radius of the stellar core, R c , was chosen to be R c = 2 × ( a log [ M ∗ /M ⊙ ] − b ) pc , (1)where the dimensionless constants a = 0 . b =1 . r eff = 0 . R c , for late- and early-type galaxies in Figure 13in Georgiev & B¨oker (2014). Note that we assumed a con-stant stellar mass-to-light ratio, M ∗ /L V , where L V is thetotal luminosity of the nuclear star cluster in the V band.The radius R T of the torus was chosen such that thesizes of the torus and of the nuclear star cluster are compa-rable, with R T = 5 pc for M ∗ = 2 × M ⊙ (see table 1 inSchartmann et al. 2005): R T = 3 . × − R c + 2 .
33 pc . (2)For all other parameters in our torus simulations, wefollowed the methods in Schartmann et al. (2005). Specifi-cally, we set the outer radius of the torus R out = 3 R c , theexponent of the angular momentum distribution in the stel-lar core γ = − .
5, and the turbulent velocity of the cloudsbuilding up the torus v t ≈ σ ∗ , where we used equation (9) in Figure 3.
An example for the adopted geometry and interiordensity structure along a cross-section of a torus surrounding anAGN, based on Schartmann et al. (2005). The mass of the nuclearstar cluster containing the torus and producing its dust contentwas chosen to be M ∗ = 10 M ⊙ , while the mass of the centralSMBH was set to M = 10 M ⊙ . Schartmann et al. (2005) to calculate the velocity dispersionof the stars in the nuclear star cluster, σ ∗ .The density distribution of a torus in this model is fullydetermined by two parameters: the mass of the nuclear starcluster, M ∗ , and the mass of the central SMBH, M . Usingthe random values of M ∗ and M , for each torus we calculatedˆ ρ d ( R, z ) ≡ ρ d ( R, z ) /ρ (see Eq. 8 in Schartmann et al. 2005)in a R out × R out rectangular cross-section along the R − z plane with a linear resolution of R out / ρ d ( R, z ) = 0wherever ˆ ρ d ( R, z ) < ˆ ρ d ( R out , ρ such thatthe total mass of the torus with the resulting ρ d ( R, z ) equals M dust . In our simulation, only the dust column mass density,Σ dust , was recorded as the final output, which was calculatedby numerically integrating ρ d ( R, z ) along the line of sightfrom each SMBH position to the observer for a randomlyoriented torus.For illustration, Fig. 3 shows a visualization of the ge-ometry and interior density structure ( ρ d ) along a cross-section of a torus with M ∗ = 10 M ⊙ and M = 10 M ⊙ .Additionally, in Fig. 4, we show the histogram of the outerradii, R out , of the 2 × dust tori. This shows that theouter radii follow a ∝ /R out distribution, ranging from R out , min ≃
10 pc to R out , max ≃
180 pc.
We are now ready to present the results of the Monte Carlosimulation of 2 . × random observations of recoilingQSOs. In Fig. 5, we show the histogram of dust columnmass densities (Σ dust ) obtained by integrating the mass den-sity of the randomly oriented dust tori along the line-of-sightto the recoiling SMBH (dark blue curve). For reference, wealso show a Σ dust histogram for a second set of simulationsof 2 . × QSOs with recoil velocities set to v recoil = 0(light red curve). We have found that the total number ofunobscured QSOs (i.e. with Σ dust = 0) is ∼ ,
700 ( ≃ ∼ ,
200 ( ≃
18 per cent) for the v recoil > c (cid:13) , 1–9 method to search for recoiling SMBHs in AGNs
10 100 2001001000 R out [pc]
Figure 4.
Histogram of outer radii, R out , of the 2 × simulateddust tori. R out values follow a ∝ /R out distribution, rangingfrom R out , min ≃
10 pc to R out , max ≃
180 pc. (cid:2) (cid:3) (cid:4) Σ dust [(cid:5)(cid:6)(cid:7) ] v recoil > (cid:8) (cid:9)(cid:10)(cid:11)(cid:12) v recoil = 0 km/s Figure 5.
Histogram of dust column mass densities (Σ dust ) for2 . × random observations of QSOs, obtained by integratingthe mass density of the randomly oriented dust tori along theline of sight to the recoiling, off-centre SMBHs (dark blue curve).For reference, we show the same histogram for a similar sample ofQSO without any recoil ( v recoil = 0; light red curve). The recoilingQSOs are much more likely to be seen completely unobscured(there are ∼ ≃
61 per cent) and ∼ ≃
18 per cent)QSOs with Σ dust = 0 in the samples with and without recoil,respectively). The recoiling QSOs also show a tail of high Σ dust values up to twice the maximum value for non-recoiling QSOs, asexpected (see text). and v recoil = 0 samples, respectively. QSOs in the v recoil = 0sample do not have Σ dust values above Σ dust ≃ / m be-cause such QSOs are never obscured by more than a halfcross-section of their dust torus, and a half cross-section ofany dust tori in the Schartmann et al. (2005) model have amaximum possible dust column density of Σ torus ≃ / m .Recoiling QSOs ( v recoil > / s) can obtain up to twicehigher values of Σ dust (i.e. up to Σ dust ≃
12 g / m ), whichcorresponds to the maximal obscuration of a QSO locatedin the ‘equatorial plane’, behind the dust torus (and for thegeometrically largest tori). This configuration is rare, withonly 409 ( ≃ . v recoil > dust > / m .Fig. 6 shows Σ dust versus | ∆ v | for the 2 . × recoilingQSOs. We have calculated the sample correlation coefficient, r , between Σ dust and log ( | ∆ v | ), defined by r ≡ h Σ dust log ( | ∆ v | ) i − h Σ dust ih log ( | ∆ v | ) i p h Σ i − h Σ dust i p h log ( | ∆ v | ) i − h log ( | ∆ v | ) i , (3)restricted to various ranges of | ∆ v | min | ∆ v | | ∆ v | max .Here h ... i refers to averaging over the sample of 2 . × QSOs, or its various subsets. Fig. 7 shows r as a function of | ∆ v | min (left-hand pan-els) and | ∆ v | max (right-hand panels) for sub-samples of ob-scured QSOs with | ∆ v | > | ∆ v | min and | ∆ v | < | ∆ v | max ,respectively. As can be seen in the figure, obscured QSOswith | ∆ v | < / s show no correlation between their Σ dust and | ∆ v | values ( r < ≃
0, and the corresponding p –valuefor the hypothesis of no correlation is p ≃ . | ∆ v | > / s show a significant correlation,with r = 0 . ± .
01 and p ≃
0. The correlation coeffi-cient increases with | ∆ v | min until it reaches its maximum at | ∆ v | min = 45 km / s with r = 0 . ± . | ∆ v | min and | ∆ v | max to vary simultaneouslydoes not change the result that r is the highest for the 3 , | ∆ v | >
45 km / s. Depending on the sig-nal to noise, the measurement error on | ∆ v | of individualQSOs is a few tens of km/s (for example from spectral linefitting, or by using the centroids of individual broad lines;see e.g., Ju et al. 2013; Shen et al. 2013). We therefore re-strict our further investigations to the 3824 obscured QSOswith | ∆ v | >
45 km / s. For this subset, we have found the ap-proximate relation Σ dust (log {| ∆ v |} ) = (1 . ± .
1) g / m × (log {| ∆ v | / [1 km / s] } ) − (1 . ± .
2) g / m , although ofcourse this relation depends strongly on the underlying setof models. We also note that in addition to the randommeasurement errors on | ∆ v | , broad line centroids can shiftby O (100 km/s) between pairs of observations only weeksapart, for physical reasons unrelated to recoil (Ju et al. 2013;Shen et al. 2013). These shifts, as well as the measurementerrors, are not expected to correlate with Σ dust . However,they introduce random errors on the Σ dust (log ( | ∆ v | )) re-lation, and can therefore increase the number of QSOs re-quired to detect the underlying Σ dust –log ( | ∆ v | ) correla-tion (see discussion below).We next estimate the minimum number of QSOs, N min ,required for a significant detection of the correlation. Todo this, for each trial N ∈ [10 , N QSOs, se-lected from among the 3824 obscured QSOs recoiling with | ∆ v | >
45 km / s, (ii) computed the p -value of the hypothesisof no correlation between Σ dust and | ∆ v | for each of these1000 samples, and (iii) measured the fraction of samples forwhich p < × − (corresponding to the rejection of theno-correlation hypothesis at 3 σ confidence). This fraction, F σ , is shown in Fig. 8 (blue crosses, right-hand panel), to-gether with the average h r i over the 1000 random samples,as a function of N (blue crosses, left-hand panel).As a test of these statistics, we have repeated the cal-culations of both h r i and F σ for two additional sequencesof 1000 random samples of N ∈ [10 , | ∆ v | < / s; theresults are shown by the red empty circles in Fig. 8). In thesecond set, we computed random Σ dust values for stationarySMBHs ( v recoil = 0), and associated a | ∆ v | with each case,drawn randomly from the | ∆ v | distribution for the recoiling( v recoil >
0) SMBHs. The results for these ‘scrambled’ data”are shown in Fig. 8 by the filled green circles.As can be seen in Fig. 8, neither the | ∆ v | < / s,nor the ‘scrambled data’ sets show correlations, i.e. h r i ≃ F σ ≃ N . On the other hand, for obscuredQSOs with | ∆ v | >
45 km / s, we find h r i ≃ .
28 for all N s,and F σ increases monotonically with N from F σ ≃ F σ ≃ N ∈ [10 , c (cid:13) , 1–9 Raffai, Haiman, & Frei
Figure 6.
The line-of-sight dust column density, Σ dust , versusthe (absolute value) of line-of-sight velocity offset | ∆ v | relativeto the host galaxy, for each of the 2 . × recoiling QSO. Forcompleteness, we show the entire range of velocities we simulated,down to unmeasurably low values. To illustrate how far each re-coiling QSO has travelled from the galactic centre by the timeof the mock observation, we sorted the data into three groups,based on their radial offset d from the galactic centre: QSOs withthe smallest offset ( d < d < d > dust values: Σ dust = 0 (correspond-ing to the ejected QSOs with the largest offsets, viewed in anunobscured direction), Σ dust ≃ / m (the maximum possibleobscuration in the Schartmann et al. (2005) model by half of thetorus, for QSOs whose orbits have fully decayed back to the centreby dynamical friction and are viewed along the equator), and atΣ dust ≃
12 g / m (the maximum possible obscuration, for QSOslocated behind the full torus). r (cid:21)(cid:22) (cid:25)(cid:26) (cid:27)(cid:28) (cid:29)(cid:30) (cid:31) | ∆ v min [km/s] p !" r % & ’ ( ∆ )* max [km/s] p Figure 7.
Sample correlation coefficient r between Σ dust andlog ( | ∆ v | ), with its 1 σ uncertainty (upper panels, blue curves),and the corresponding p values for the hypothesis of no correlation(lower panels, red curves). The left-hand (right-hand) panels showboth r and p as a function of | ∆ v | min ( | ∆ v | max ), for sub-samplesof obscured QSOs with | ∆ v | > | ∆ v | min ( | ∆ v | < | ∆ v | max ). QSOswith | ∆ v | ∼ > / s show a positive correlation, with r as high as r = 0 . ± .
02 for | ∆ v | min ≈
45 km / s. that the hypothesis of no correlation between Σ dust andlog ( | ∆ v | ) could be rejected with a 3 σ significance 95 percent of the time by observing a random sample of ∼ | ∆ v | >
45 km / s. Note that in our sim-ulations, 3824 QSOs ( ≃ . | ∆ v | >
45 km / s, thus gathering a sample of ∼
260 suchQSOs may require observing as many as ∼ +,- F σ QSOs in sample | ∆ v| ≥
45 km/s| ∆ v| < 5 km/sScrambled data < r > | ∆ v| ≥
45 km/s| ∆ v| < 5 km/sScrambled data Figure 8.
Left-hand panel: the expectation value h r i of the cor-relation coefficient between Σ dust and log ( | ∆ v | ), computed in1000 random samples of N QSOs. The results are shown as a func-tion of N in the range N ∈ [10 , | ∆ v | >
45 km / s (blue crosses), obscuredQSOs with | ∆ v | < / s (red circles), and QSOs without recoil,but with | ∆ v | values assigned randomly from the | ∆ v | distribu-tion for recoiling SMBHs (‘scrambled data’; filled green circles).While h r i ≃ N in the last two data sets, the first data setshows a stable correlation with h r i ≃ .
28 for all N . Right-handpanel: the fraction F σ of the 1,000 random samples in which the‘no correlation’ null hypothesis can be rejected at > σ signifi-cance. While F σ ≃ | ∆ v | < / s and ‘scrambled’ datasets, it increases monotonically with N from F σ ≃ | ∆ v | >
45 km / s data set, rising above F σ > .
95 for N > for the fact that not all QSOs experienced a recoil in thepast).To examine the impact of the sign of the line-of-sightvelocity offset, in Fig. 9 we show Σ dust versus | ∆ v | in thesubsample of 17157 obscured QSOs with | ∆ v | > / s,and indicate whether they are moving away from (∆ v > v <
0, shown as redfull circles) the observer. First, we note that receding QSOsare slightly more common, with an excess of 1405 and 1162in the | ∆ v | > / s, and | ∆ v | >
45 km / s samples, respec-tively, over the approaching cases. As Fig. 9 shows, Σ dust val-ues of the receding QSOs also tend to be higher, especiallywhen the velocity offset is large ( | ∆ v | ∼ >
100 km / s). This isbecause receding QSOs are more likely to be found behind,rather than in front of the tori. In order for a SMBH to beapproaching with a large line-of-sight velocity while beinglocated behind the torus, it would have had to reach its firstU-turn. As mentioned above, the majority of SMBHs withlarge recoil velocities do not make their first U-turn withina QSO lifetime (or within our mock observation times).In Fig. 10, we compare the numerically constructedprobability density functions of various parameter distribu-tions of QSOs within the whole sample, and within the sub-sample of obscured QSOs with | ∆ v | >
45 km / s. Based onthe comparison, we conclude that AGNs where the centralSMBH has a mass around M ≃ (0 . − − × M ⊙ , apositional offset from the centre of the host galaxy between d offset ≃ − − pc, and having a nuclear star cluster withmass M ∗ & M ⊙ , have the highest probability of being inthe subset of obscured QSOs with | ∆ v | >
45 km / s. The factthat high-mass nuclear star clusters are over-represented inthe subset of QSOs showing Σ dust – | ∆ v | correlation (see theupper-right panel in Fig. 10) suggests that such QSOs aremore likely to be found in active galaxies with higher masses.We attribute this to the fact that in more massive galaxies, c (cid:13) , 1–9 method to search for recoiling SMBHs in AGNs . / | ∆
27 89:;<= Σ d?@ABCDE ] ∆ v > 0 km/s ∆ v < 0 km/s Figure 9. Σ dust versus | ∆ v | for the 17157 obscured QSOs with | ∆ v | > / s. Receding (i.e. ∆ v >
0) QSOs are shown byblue crosses, while approaching QSOs (i.e. ∆ v <
0) are shownby red filled circles. We also show the Σ dust (log ( | ∆ v | )) =(1 . ± .
1) g / m × (log ( | ∆ v | / [1 km / s])) − (1 . ± .
2) g / m curve (green solid line) we obtained by fitting Σ dust versuslog ( | ∆ v | ) values for obscured QSOs with | ∆ v | >
45 km / s. Re-ceding SMBHs are slightly more common among the obscured(Σ dust >
0) QSOs, their velocity offsets extend to higher values,and they tend to be more heavily obscured (larger Σ dust ), espe-cially at large velocity offsets (∆ v ∼ >
100 km / s). See text for theorigin of these differences. F G H
IJKL MN OPSQR ]) TUV W log XYZ \^_‘a ]) b c log (v recoil [km/s]) efg h log (d oijklm npqrs Figure 10.
The probability distribution functions of different pa-rameters of the QSOs within the entire sample (blue curves), andwithin the subsample of obscured QSOs with significant correla-tions (i.e. with | ∆ v | >
45 km / s; red curves). The four differentpanels, clockwise from top left, show the SMBH masses ( M ), themasses of the nuclear star clusters ( M ∗ ), recoil velocities ( v recoil ),and the offsets of QSOs from the centres of their host galax-ies ( d offset ). The difference between the red and blue curves re-veal that QSOs with M ≃ (0 . − − × M ⊙ , M ∗ & M ⊙ , v recoil ≃
100 km / s, and d offset ≃ − − pc are the mostover-represented in the sample showing strong correlations be-tween Σ dust and | ∆ v | . the recoiling BHs spend longer times at off-centre distancescomparable to the size of the obscuring torus. This host-sizedependence of the proposed correlation is another potentialtestable prediction of our toy model.Finally, we have studied the fraction of obscured QSOswithin the whole sample of QSOs, F obs ( | ∆ v | ). We havefound that this fraction decreases monotonically with | ∆ v | from F obs ( <
10 km / s) & . F obs (cid:0) > km / s (cid:1) . . t uwx yz{ } ∆ ~(cid:127) (cid:128)(cid:129)(cid:130)(cid:131)(cid:132)(cid:133)(cid:134) (cid:135)(cid:136)(cid:137) Figure 11.
The fraction of obscured (Σ dust >
0) sourcesamong the entire sample of 2 . × QSOs. The trend of F obs ( | ∆ v | ) decreasing with | ∆ v | from F obs ( <
10 km / s) & . F obs (cid:0) > km / s (cid:1) . . observed fraction of obscured (e.g. type II) QSOs, which canpotentially provide an independent test for a chosen combi-nation of recoil, trajectory, and dust tori models. The results in the previous section show that a positive cor-relation may be present between Σ dust and | ∆ v | . We em-phasize that this conclusion is based on a specific set ofidealized toy models – more realistic models will undoubt-edly lead to different predictions. Nevertheless, our resultsshould motivate an observational search; more realistic andflexible models will likely be needed to interpret any correla-tion discovered in actual data. We plan to present the resultsfrom a search in the SDSS DR10 data base in a forthcomingpublication.Here we address two natural caveats: while proxies canbe used to estimate both Σ dust and | ∆ v | , these estimates willhave uncertainties we have so far ignored. Although we haveno a priori reason to expect that these uncertainties will becorrelated, they will represent noise, which will reduce thestatistical significance of any correlation, and thus increasethe number of QSOs required for a detection. Furthermore,we have so far assumed that a clean QSO sub-sample withΣ dust > | ∆ v | >
45 km / s can be constructed. Morerealistically, observational errors will pollute any such selec-tion (e.g. QSOs with | ∆ v | >
45 km / s can be missed, andQSOs with | ∆ v | <
45 km / s can be mistakenly included).The line-of-sight velocity offset ( | ∆ v | ) can be estimatedusing the centroids or peaks of broad lines, relative to thoseof the narrow lines. The level of obscuration (Σ dust ) canbe estimated using QSO colours: a larger obscuring columnshould result in a redder continuum (see e.g. Ledoux et al.(2015) and references therein). An alternative proxy for theobscuring column could be the equivalent width of the broadlines themselves. Although the geometry of the broad lineregion is not fully understood, type II quasars have weakor no broad emission lines. In the standard unified modelof AGN (Antonucci 1993), this is attributed to the obscura-tion of the high-velocity broad-line region by the dusty torus(note that the continuum is still visible and can therefore not c (cid:13) , 1–9 Raffai, Haiman, & Frei be spatially colocated with the broad-line emitting region).This implies that the equivalent width of these broad linesmust vary with the level of the obscuration.Here we estimate the impact of having random errors onboth | ∆ v | and Σ dust . As mentioned above, we expect | ∆ v | ,when determined using line centroids, to suffer physical vari-ations of order σ ∆ v ≈
100 km / s. We therefore repeated ouranalysis above, except we added a random additional ∆ v component to the velocity offset of each simulated QSO,drawn from a Gaussian distribution with a standard devia-tion of σ ∆ v . For any given σ ∆ v , we then selected the subsetof the 2 . × QSOs with Σ dust > | ∆ v | >
45 km / s,as before. Finally, we measured the expectation value of thecorrelation coefficient h r i between Σ dust > | ∆ v | in1000 randomly drawn sets of N QSOs from this subset, anddetermined the minimum number N min of QSOs required toinfer r > σ significance 95 per cent of the time.The result of this random-velocity-error exercise isshown in Fig. 12 for both Σ dust –log ( | ∆ v | ) and Σ dust – | ∆ v | correlations. The right-hand panel in this figure shows N min as a function of the random error σ ∆ v . As expected, thenumber of QSOs required to detect a correlation increasemonotonically with σ ∆ v , rising to N min = 3560 for theΣ dust – | ∆ v | correlation and for σ ∆ v = 100 km / s; this is morethan an order of magnitude increase over the N min = 260in the idealized case for the Σ dust –log ( | ∆ v | ) correlationwithout any velocity-offset error. The left-hand panel inFig. 12 shows the corresponding values of the correlationcoefficients; h r i decreases monotonically with σ ∆ v , reaching h r i = 0 . ± .
016 for the Σ dust – | ∆ v | correlation and for σ ∆ v = 100 km / s. This is a decrease from the h r i = 0 . ± . dust –log ( | ∆ v | ) correlation bya factor of three. Additionally, we have calculated h r i forpartially obscured QSOs only (i.e. QSOs with 0 < Σ dust ∼ < . / m ), again assuming σ ∆ v = 100 km / s, and we havefound it to be h r i = 0 . ± . | ∆ v | and Σ dust .Similarly to the errors in | ∆ v | , we have studied theimpact of errors in Σ dust . Here, there are several poten-tial concerns. First, in the presence of uncertainties, asharp selection Σ dust > dust comparable to the maxi-mum possible column density of a half-torus according tothe Schartmann et al. (2005) model (see Fig. 6), Σ torus ≃ / m (i.e. for configurations corresponding to true typeII QSOs), the broad lines may be undetectable, or mayhave too low S/N for a reliable centroid measurement. We here repeated our fiducial analysis, except that we re-placed the selection criterion Σ dust > min < Σ dust < Σ max . We chose Σ min = 0 . / m ≈ . × Σ torus andΣ max = 2 . / m ≈ . × Σ torus . The number of QSOs withinthis range (and also satisfying ∆ v >
45 km/s as before) is N = 1436, or ∼ . . × QSOs. With this cut, we find that the correlations remainsignificant: h r i = 0 . ± .
026 for the Σ dust –log ( | ∆ v | )correlation, with a p -value of p = 2 × − .We also evaluated the impact of random Gaussian errors Note that the broad lines and the optical continuum arethought to arise from different spatial locations, and can sufferdifferent levels of obscuration; we ignore this complication here. σ ∆ v [km/s] < r > ( F σ = . ) log (| ∆ v|)| ∆ v| 0 50 100 15010 σ ∆ v [km/s] N m i n ( F σ = . ) log (| ∆ v|)| ∆ v| Figure 12.
The figure shows the impact of random errors in thevelocity offset; the error is assumed to be Gaussian with a stan-dard deviation of σ ∆ v . The right-hand panel shows the increasein the minimum number of QSOs ( N min ) required to detect cor-relations between Σ dust and log ( | ∆ v | ) (blue solid curves), orbetween Σ dust and | ∆ v | (red dashed curves) in a subset of QSOswith Σ dust > | ∆ v | >
45 km / s, as σ ∆ v increases. Curves inthe left-hand panel, with ± σ errors, show the corresponding de-crease in the correlation coefficient h r i (see text for more details).Note that measurement errors below σ ∆ v
45 km / s results withhigher h r i and lower N min values for the Σ dust –log ( | ∆ v | ) cor-relation, while the Σ dust – | ∆ v | correlation appears to be strongerwhen σ ∆ v >
45 km / s. in Σ dust > σ Σ . The results aresimilar to those shown in Fig. 12: we find that h r i decreases,and the minimum number of QSOs required to detect thecorrelation between Σ dust and log ( | ∆ v | ) ( N min ) increaseswith σ Σ . For the case of σ Σ = 0 . / m = 0 . × Σ torus , wefind h r i = 0 . ± .
05, with a p -value of p ≈ − , and N min = 290. Thus, random Gaussian errors in Σ dust > dust –log ( | ∆ v | ) correlation.Based on the above, we conclude that statistical errorson estimating the velocity offsets and obscuring column den-sities can diminish the correlations predicted in the idealizedtoy model. This is mostly due to the impact of these errorson the purity of selection of a sub-sample, intended to in-clude only the most correlated QSOs. In the examples above,we find that this increases the number of QSOs required fora detection by an order of magnitude, but there is no in-dication that these errors will prevent a detection of thecorrelations in a large sample (tens of thousands) of QSOs.Finally, as the distribution of SMBH masses in QSOsobservable in different survey projects may vary signifi-cantly, we have considered how using higher SMBH masses(see Greene & Ho 2007, 2009 and Shen & Kelly 2012) af-fect our results. In these tests, we have also taken into ac-count the possibility of a correlation between the SMBHmass and the mass of the stellar core (Ferrarese et al. 2006;Graham et al. 2011; Kormendy & Ho 2013). We have foundthat the increased SMBH masses lead fewer SMBHs escap-ing the galactic nucleus (an effect that should increase thenumber of SMBHs in the population showing Σ dust – | ∆ v | correlation), however they also lead to significantly smallertorus sizes, resulting with a net decrease in the correlationcoefficient of a factor of ∼
2. This demonstrates that al-though the strength of the Σ dust – | ∆ v | is sensitive to the un-derlying assumptions of the applied models, the Σ dust – | ∆ v | correlation itself is a robust result. c (cid:13) , 1–9 method to search for recoiling SMBHs in AGNs A large fraction of luminous QSOs are believed to be trig-gered by merger events, and if QSO activity often followsthe coalescence of SMBHs in the merging nuclei, then a largefraction of all QSOs may have experienced recent gravita-tional recoil. While a handful of QSOs with large spatialor velocity offsets have been identified as recoil candidates(Komossa 2012), the ubiquity of this phenomenon remainspoorly known.Here, we proposed a statistical technique to search fora population of recoiling SMBHs with smaller velocity off-sets among luminous QSOs, based on a positive correlationbetween these velocity offsets and the column densities ofobscuring dust tori. The correlations are introduced by thedamped oscillating motions of SMBHs with typical recoilspeeds, which are comparable to the escape velocities fromgalactic nuclei. These SMBHs spend a significant time off-nucleus, and at distances comparable to the size of obscuringtori in located in the nuclei.We have demonstrated, using simple models of gravi-tational recoil, SMBH trajectories, and the geometry of ob-scuring dust tori, that observing a random sample of a fewthousand partially obscured QSOs, with line-of-sight veloc-ity offsets | ∆ v | >
45 km / s, could allow a significant detec-tion of the correlation between the line of sight dust columndensity Σ dust and | ∆ v | (or log ( | ∆ v | )).The existence of this correlation is found to be robustwithin the simple models we choose, and detectable even inthe face of random errors on Σ dust and | ∆ v | . However, weexpect that the tightness of the correlation, the underlyingΣ dust − | ∆ v | relation, and the properties of the QSO subsetshowing correlations can depend strongly on the model as-sumptions. Nevertheless, our results should motivate search-ing for correlations in real data; a positive detection wouldallow testing various combinations of recoil, BH trajectory,and dust tori models. A further test could be possible bycomparing the predicted fraction of obscured QSOs to theobserved fraction of type II QSOs. We propose to carry outthese tests in the future using catalogues of observationaldata for QSOs, such as the SDSS-DR10 QSO catalogue. ACKNOWLEDGEMENTS
The authors would like to thank Bence B´ecsy, GergelyD´alya, and ´Akos Sz¨olgy´en for useful discussions and for valu-able comments on the manuscript. P´eter Raffai is grateful forthe support of the Hungarian Academy of Sciences throughthe ’Bolyai J´anos’ Research Scholarship programme. We alsoacknowledge financial support from the NASA AstrophysicsTheory Program under grant no. NNX11AE05G (to ZH)and from OTKA under grant no. 101666 (to ZF).
REFERENCES
Aller M. C., & Richstone D., 2002, AJ, 124, 3035Antonucci R., 1993, Ann. Rev. A. & A., 31 473Baker J. G. et al., 2008, ApJ, 682, L29Barth A. J., Greene J. E., & Ho L. C., 2005, ApJ, 619,L151 Begelman M. C., Blandford R. D., & Rees M. J., 1980,Nature, 287, 307Blecha L. et al., 2011, MNRAS, 412, 2154Blecha L. et al., 2015, preprint (arXiv:1508.01524)Bogdanovi´c T., Reynolds C. S., & Miller M. C., 2007, ApJ,661, L147Bonning E. W., Shields G. A., & Salviander S., 2007, ApJ,666, L13Comerford J. M. et al., 2009, ApJ, 698, 956Dubois Y. et al., 2014, MNRAS, 440, 2333Ferrarese L. et al., 2006, ApJ, 644, L21Georgiev I. Y., & B¨oker T., 2014, MNRAS, 441, 3570Gerosa D. et al., 2015, MNRAS, 451, 3941Graham A. W. et al., 2011, MNRAS, 412, 2211Greene J. E., & Ho L. C., 2006, ApJ, 641, L21Greene J. E., & Ho L. C., 2007, ApJ, 667, 131Greene J. E., & Ho L. C., 2009, ApJ, 704, 1743Gualandris A., & Merritt D., 2008, ApJ, 678, 780Guedes J. et al., 2009, ApJ, 702, 890Guedes J. et al., 2011, ApJ, 729, 125Haehnelt M. G., 1994, MNRAS, 269, 199Healy J., Lousto C. O., & Zlochower Y., 2014, Phys. Rev.D, 90, 104004H¨onig S. F., 2008, Dr. rer. nat. dissertation, RheinischeFriedrich-Wilhelms-Universit¨atIvanov P. B., Papaloizou J. C. B., & Polnarev A. G., 1999,MNRAS, 307, 79Ju W. et al., 2013, ApJ, 777, 44Komossa S., 2012, Adv. Astron., 2012, 364973Komossa S., & Merritt, D., 2008, ApJ, 689, L89Kormendy J., & Ho L. C., 2013, ARA&A, 51, 511Ledoux C. et al., 2015, A&A, 580, A8Leigh N., B¨oker T., & Knigge C., 2012, MNRAS, 424, 2130Loeb A., 2007, Phys. Rev. Lett., 99, 041103Lousto C. O. et al., 2012, Phys. Rev. D, 85, 084015Madau P., & Quataert E., 2004, ApJ, 606, L17Martini P., 2004, in Ho L. C., ed., Coevolution of BlackHoles and Galaxies, Cambridge Univ. Press, Cambridge,p. 170Pˆaris I. et al., 2014, A&A, 563, A54Reynolds C. S., 2013, Class. Quantum Gravity, 30, 244004Schartmann M. et al., 2005, A&A, 437, 861Sesana A., Haardt F., Madau P., & Volonteri M. 2005, ApJ,623, 23Shen Y., & Kelly B. C., 2012, ApJ, 746, 169Shen Y., Liu X., Loeb A. & Tremaine S., 2013, ApJ, 775,49Tanaka T. & Haiman Z., 2009, ApJ, 696, 1798Tremaine S. et al., 2002, ApJ, 574, 740Volonteri M., Haardt F., & Madau P., 2003, ApJ, 582, 559 c (cid:13)000