Lighthouse in the Dust: Infrared Echoes of Periodic Emission from Massive Black Hole Binaries

TL;DR

This paper models infrared echoes of periodic emission from massive black hole binaries to distinguish between isotropic and anisotropic variability, aiding in identifying and characterizing such systems.

Contribution

It introduces models for IR reverberation of periodic AGN emission, accounting for dust torus geometry and Doppler effects, to interpret IR variability signatures.

Findings

IR variability can occur without optical/UV signals in Doppler-boosted cases.

The dust torus is constrained to be thin with an inner radius of 1-5 parsecs.

Models are consistent with WISE IR data for PG 1302-102.

Abstract

The optical and UV emission from sub-parsec massive black hole binaries (MBHBs) in active galactic nuclei (AGN) is believed to vary periodically, on timescales comparable to the binary's orbital time. If driven by accretion rate fluctuations, the variability could be isotropic. If dominated by relativistic Doppler modulation, the variability should instead be anisotropic, resembling a rotating forward-beamed lighthouse. We consider the infrared (IR) reverberation of either type of periodic emission by pc-scale circumbinary dust tori. We predict the phase and amplitude of IR variability as a function of the ratio of dust light crossing time to the source variability period, and of the torus inclination and opening angle. We enumerate several differences between the isotropic and anisotropic cases. Interestingly, for a nearly face-on binary with an inclined dust torus, the Doppler boost can produce IR variability without any observable optical/UV variability. Such orphan-IR variability would have been missed in optical searches for periodic AGN. We apply our models to time-domain WISE IR data from the MBHB candidate PG 1302-102 and find consistency with dust reverberation by both isotropically emitting and Doppler-boosted sources in the shorter wavelength W1-W2 (2.8-5.3 micrometer) bands. We constrain the dust torus to be thin (aspect ratio ~0.1), with an inner radius at 1-5 pc. More generally, our dust echo models will aid in identifying new MBHB candidates, determining their nature, and constraining the physical properties of MBHBs and their dust tori.

Tables

Table 1. Table 1: Parameters of the model and their fiducial values used in dust echo calculations unless stated otherwise.

Parameter	Meaning	Fiducial Value	Notes
Source Parameters
$L^0$	Bolometric luminosity	$7\times {10}^{46}$ erg s-1	Value inferred for PG 1302-102
$BC$	Bolometric correction from V band	$9$	Typical value for quasars
$P$	Variability period	$R_d/c$	Match light crossing time of torus
— — Additional source parameters needed only in Doppler case: — —
${\alpha }_{\nu }$	“in-band” Spectral index	1	Approximate PG 1302-102 value in optical
${\alpha }_{\mathrm{bol}}$	“bolometric” Spectral index	-1	Effective spectral index when integrating
			over entire spectrum (Eq. 4)
$\beta$	Orbital velocity$/c$	$0.07$	Approximate PG 1302-102 value for fiducial ${\alpha }_{\nu }$
$I$	Orbital inclination	$0$	0=edge-on; $\pi /2$=face-on
Dust Parameters
$R_d$	Inner edge of dust	$1.0$ pc	Approximate sublimation radius for PG 1302-102
$k$	Absorption/emission efficiency exponent	1.6	Approximate value for graphites (Draine & Lee, 1984)
${\nu }_0$	Efficiency cutoff frequency	$c/\left(2\pi a_{\mathrm{eff}}\right)$	Needed for efficiency $Q_{\nu }=\min [{\left(\nu /{\nu }_0\right)}^k,1]$
$a_{\mathrm{eff}}$	Grain size	$0.16$ $\mu$m	Set by ${\nu }_0$
— — Additional dust parameters needed only for torus case: — —
$J$	Inclination	$\pi /2$	0=edge-on; $\pi /2$=face-on
${\theta }_T$	Opening angle	$\pi /4$	0=sphere ($\cos {\theta }_T$=covering fraction)

Table 2. Table 2: Best-fit sinusoid parameters for model 𝒜 ​ sin ⁡ Ω ​ ( t − t 0 ) + ℬ 𝒜 Ω 𝑡 subscript 𝑡 0 ℬ \mathcal{A}\sin{\Omega\left(t-t_{0}\right)}+\mathcal{B} , with Ω = 2 ​ π / 1884 Ω 2 𝜋 1884 \Omega=2\pi/1884 day -1 . The reduced χ 2 superscript 𝜒 2 \chi^{2} statistic for the optical fit is consistent with that found by (D’Orazio et al., 2015b ) , but without the inclusion of damped random walk noise. The reduced χ 2 superscript 𝜒 2 \chi^{2} statistic is smaller for the W1 and W2 data because we have binned each epoch into one data point as explained in Appendix B.1 .

	V	W1	W2
$𝒜$ [mag]	${0.1269}_{-0.0004}^{+0.0004}$	${0.0862}_{-0.0271}^{+0.0158}$	${0.0813}_{-0.0260}^{+0.0157}$
$t_0$ [days]	${298.1}_{-0.7}^{+0.7}$	${715}_{-32}^{+26}$	${768}_{-27}^{+28}$
$ℬ$ [mag]	${14.8329}_{-0.0003}^{+0.0003}$	${11.39}_{-0.02}^{+0.01}$	${10.31}_{-0.02}^{+0.01}$
${\chi }_{\mathrm{red}}^2$	40.5	2.0	0.4

Table 3. Table 3: Measured quantities used to determine the nature of the IR light curves of PG 1302-102. The table is divided into three sections. Section I displays the quantities measured directly from the light curves with no model assumptions: the ratios of IR to V-band amplitudes, phase lags, and period. Section II displays quantities found from fitting a blackbody spectrum to WISE band fluxes. Case (i) utilizes W1, W2, and W3 bands, case(ii) uses only W1 and W2 bands. Section III displays quantities found from incorporating the quantities measured in the first two sections of the table into our dust-echo models. Different results are found depending on the source model: an isotropically varying source (Iso) or a Doppler-boosted source (Dop), and the data used in the spectral fit: case (i) and (ii) above. Note that t d / P subscript 𝑡 𝑑 𝑃 t_{d}/P is constrained both by the Φ Φ \Phi measurements in the first section of the table and also by R d subscript 𝑅 𝑑 R_{d} in the second section. Because the Φ Φ \Phi measurements yield t d / P subscript 𝑡 𝑑 𝑃 t_{d}/P values modulo half of a variability cycle (Eqs. 17 and 15 ), the combined constraint on t d / P subscript 𝑡 𝑑 𝑃 t_{d}/P can consist of multiple disconnected regions. The A IR / A subscript 𝐴 IR 𝐴 A_{\rm{IR}}/A value in section III combines both relative amplitude measurements in the first section of the Table and the source type which sets A 𝐴 A .

I Model Independent Measures
$\frac{A_{\mathrm{W1}}}{A_V}$	$\frac{A_{\mathrm{W2}}}{A_V}$	${\Phi }_{\mathrm{W1}}$ [cycles]	${\Phi }_{\mathrm{W2}}$ [cycles]	$P$ [days]
${0.68}_{-0.21}^{+0.12}$	${0.64}_{-0.20}^{+0.12}$	$0.10\to 0.30$	$0.20\to 0.36$	${1877}_{-3}^{+22}$
		$\{\pm$ Half Integer $\}$
II Spectral Fit Measures
	$T_d$ [K]	$\cos {\theta }_T$	B.C.	$R_d$ [pc]
case (i)	$577\pm 3$	${0.08}_{-0.03}^{+0.04}$	$9\pm 3$	$3.9\pm 1.0$ ($t_d/P=3.1\pm 0.8$)
case (ii)	$881\pm 14$	${0.05}_{-0.02}^{+0.02}$	$9\pm 3$	$1.9\pm 0.4$ ($t_d/P=1.5\pm 0.3$)
III Dust-Echo Modeling Measures
	$\frac{A_{\mathrm{IR}}}{A}$	$t_d/P$	$J$ [rad]	$I$ [rad]
Iso, case (i)	$0.46\pm 0.2$	$\{2.3\to 2.36$	$J\sim >1.0$	Unconstrained
		$2.6\to 2.86$
		$3.1\to 3.36$
		$3.6\to 3.86\}$
Iso, case (ii)		$\{1.2\to 1.36$		Unconstrained
		$1.6\to 1.8\}$
Dop, case (i)	${0.31}_{-0.11}^{+0.07}$	$\{2.35\to 2.61$	Not allowed for $I=0$
		$2.85\to 3.11$		(Figure 1 D’Orazio et al. (2015b))
		$3.35\to 3.61$
		$3.85\to 3.9\}$
Dop, case (ii)		$\left\{1.35\to 1.61\right\}$

Equations

$$F_{\nu }^{\mathrm{Iso}}=F_{\nu }^0\left[1+A\sin \Omega t\right],$$

$$F_{\nu }^{\mathrm{Dop}}=\frac{F_{\nu }^0}{{\left[\gamma \left(1-\frac{v_{||}^{\mathrm{obs}}}{c}\right)\right]}^{3-{\alpha }_{\nu }}},$$

$$P=\frac{2\pi a^{3/2}}{\sqrt{GM}}.$$

$$F^{\mathrm{Dop}}=\frac{F^0}{{\left[\gamma \left(1-\frac{v_{||}^{\mathrm{obs}}}{c}\right)\right]}^{3-{\alpha }_{\mathrm{bol}}}},$$

$$\pi a_{\mathrm{eff}}^2{\bar{Q}}^{\mathrm{src}}{F}^{\mathrm{iso}}(t,R_d)=4\pi a_{\mathrm{eff}}^2{\int }_0^{\infty }{Q}_{\nu }\pi B_{\nu }\left[T_d(t)\right]d\nu .$$

$$F_{\nu }^{\mathrm{grain}}$$

$$F_{\nu }(t)$$

$$t_{\mathrm{em}}$$

$$F_W(t)={\int }_0^{\infty }W(\nu )F_{\nu }(t)d\nu \sim {\int }_{{\nu }_{\min }}^{{\nu }_{\max }}{F}_{\nu }(t)d\nu$$

$$\frac{v_{||}}{c}$$

$${\bar{Q}}^{\mathrm{src}}{\int }_0^{\infty }{F}_{\nu }^{\mathrm{Dop}}(t,R_d,\theta ,\varphi )d\nu =4{\int }_0^{\infty }{Q}_{\nu }\pi B_{\nu }\left[T_d(t,\theta ,\varphi )\right]d\nu .$$

$${\left[D(t,\theta ,\varphi )\right]}^4\frac{L^0}{4\pi R_d^2}=4{\int }_0^{\infty }{Q}_{\nu }\pi B_{\nu }\left[T_d(t,\theta ,\varphi )\right]d\nu ,$$

$$L_{\mathrm{IR}}^{\mathrm{Iso}}(t)={\Sigma }_{d}4\pi a_{\mathrm{eff}}^2{\int }_0^{2\pi }{\int }_0^{\pi }\frac{L\left(t_{\mathrm{em}}(t,\theta ,\varphi )\right)}{16\pi R_d^2}{R}_d^2\sin \theta d\theta d\varphi$$

$$={\Sigma }_d\pi a_{\mathrm{eff}}^2{L}^0\left[1+A\mathrm{sinc}\left(\mathrm{\Omega }{\mathrm{t}}_{\mathrm{d}}\right)\sin \left(\mathrm{\Omega }\left(\mathrm{t}-{\mathrm{t}}_{\mathrm{d}}\right)\right)\right].$$

$$L(t)=L^0\left[1+A\sin \left(\Omega t\right)\right].$$

$$\frac{A_{\mathrm{IR}}^{\mathrm{Iso}}}{A}=\frac{1}{2\pi }\frac{P}{t_d}\sin \left[2\pi \frac{t_d}{P}\right],$$

$${\Phi }_{\mathrm{Iso}}$$

$$L_{\mathrm{IR}}^{\mathrm{Dop}}$$

$${\Phi }_{\mathrm{Dop}}=\frac{t_d}{P}-\mathrm{sign}\left(\frac{{\mathrm{A}}_{\mathrm{IR}}^{\mathrm{Dop}}}{\mathrm{A}}\right)\frac{1}{4}\mathrm{cycles},$$

$$\frac{A_{\mathrm{IR}}^{\mathrm{Dop}}}{\beta \cos I}=\frac{\mathrm{sinc}\left(2\pi \frac{{\mathrm{t}}_{\mathrm{d}}}{\mathrm{P}}\right)}{2\pi \frac{t_d}{P}}-\frac{\cos \left(2\pi \frac{t_d}{P}\right)}{2\pi \frac{t_d}{P}},$$

$$2\pi \frac{t_d}{P}=\tan \left(2\pi \frac{t_d}{P}\right).$$

$${\frac{t_d}{P}|}_{\mathrm{zeros}}\approx \frac{2m+1}{4}-\frac{1}{{\pi }^2(2m+1)}m=1,2,\mathrm{3....}$$

$$\frac{2(2\pi t_d/P)}{(2\pi t_d/P{)}^2-2}=\tan \left(2\pi \frac{t_d}{P}\right),$$

$$L_{\mathrm{IR}}^{\mathrm{Iso}}$$

$$\displaystyle L^{\rm{Dop}}_{\rm{IR}}=\Sigma_{d}\pi a^{2}_{\rm{eff}}L^{0}\cos{\theta_{T}}\bigg{\{}1+\beta\cos{I}\bigg{.}$$

$$\times \left(\frac{\mathrm{sinc}\left[\mathrm{\Omega }{\mathrm{t}}_{\mathrm{d}}\cos {\theta }_{\mathrm{T}}\right]}{\Omega t_d}-\frac{\cos \left[\Omega t_d\cos {\theta }_T\right]}{\Omega t_d}\right)\cos \left(\Omega \left[t-t_d\right]\right)\},$$

$$\frac{\Delta L_{IR}^{\mathrm{Iso}}/(L^0\cos {\theta }_T)}{\Delta L/L^0}\equiv \frac{A_{\mathrm{IR}}^{\mathrm{Iso}}}{A}$$

$$\frac{\Delta L_{IR}^{\mathrm{Dop}}/(L^0\cos {\theta }_T)}{\Delta L/L^0}\equiv \frac{A_{\mathrm{IR}}^{\mathrm{Dop}}}{\beta \cos I}$$

$$=\frac{\mathrm{sinc}\left(2\pi \frac{{\mathrm{t}}_{\mathrm{d}}}{\mathrm{P}}\cos {\theta }_{\mathrm{T}}\right)}{2\pi \frac{t_d}{P}}-\frac{\cos \left(2\pi \frac{t_d}{P}\cos {\theta }_T\right)}{2\pi \frac{t_d}{P}}.$$

$$L_{\mathrm{IR}}^{\mathrm{Iso}}$$

$$L_{\mathrm{IR}}^{\mathrm{Dop}}$$

$$\frac{A_{W1}}{A_{\mathrm{V}}}\approx 0.68_{-0.21}^{+0.12},\frac{A_{W2}}{A_{\mathrm{V}}}\approx 0.64_{-0.20}^{+0.12}.$$

Full text

††pagerange: Lighthouse in the Dust: Infrared Echoes of Periodic Emission from Massive Black Hole Binaries††thanks: We dedicate this work to the memory of Arlin Crotts, an expert on supernova light echoes who passed away on November 19, 2015.–B.2

1 Introduction

Massive black holes (MBHs) exist at the centers of most, if not all, galaxies (Kormendy & Richstone, 1995; Kormendy & Ho, 2013). Galactic mergers can deliver MBHs, as well as an ample supply of gas (Barnes & Hernquist, 1992, 1996; Barnes, 2002; Mayer, 2013), to the centers of newly coalesced galaxies where the MBHs form a binary. The interaction of massive black hole binaries (MBHBs) with gas and surrounding stars can drive the pair to sub-pc separations where gravitational radiation reaction drives the binary to coalescence (Begelman et al., 1980). Characterization of the population of such sub-pc binaries, through present electromagnetic (EM), and future gravitational wave (GW) channels will provide a powerful tool for understanding the mutual build-up of galaxies and their central black holes (e.g. Kormendy & Ho, 2013), the dynamics of gas and stars in galactic nuclei (e.g. Merritt & Milosavljević, 2005), and the low-frequency gravitational wave sky (e.g. Kocsis & Sesana, 2011; Shannon et al., 2015; Arzoumanian et al., 2016; Ravi et al., 2015).

Electromagnetic signatures of MBHBs can arise from their interaction with gas. Hydrodynamical simulations of gas discs surrounding compact MBHBs show that accretion rates onto a binary are generally comparable to the accretion rates onto a single black hole of an equivalent mass (Shi et al., 2012; D’Orazio et al., 2013; Farris et al., 2014; Shi & Krolik, 2015; Muñoz & Lai, 2016). Binary accretion can also be identified by its unique signatures. Depending on the ratio of BH masses, $q\equiv M_{2}/M_{1}\leq 1$, the accretion induced emission can be periodically modulated at predictable frequencies (Hayasaki et al. 2007, MacFadyen & Milosavljević 2008, Cuadra et al. 2009, Roedig et al. 2012, Noble et al. 2012, Shi et al. 2012, D’Orazio et al. 2013, Farris et al. 2014, Gold et al. 2014b, Gold et al. 2014a, Dunhill et al. 2015, Shi & Krolik 2015, Farris et al. 2015, D’Orazio et al. 2016, Muñoz & Lai 2016).

Tell-tale luminosity variations caused by an accreting MBHB can also occur due to special relativity alone. For a binary orbiting at relativistic speeds, relativistic Doppler boosting causes luminosity variations at the binary orbital period (D’Orazio et al., 2015a, b). For binaries with disparate mass ratios, q\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip 1.0pt\sim}\hss}\raise 1.0pt\hbox{<}}0.05, accretion has been found to be steady and dominated by the smaller, secondary BH (D’Orazio et al., 2013; Farris et al., 2014). In this case, for a sufficiently compact binary, Doppler boosting of the faster-moving secondary BH’s emission is expected to be the dominant source of variability. Even near-equal mass binaries, for which high- amplitude accretion variability is expected, may emit steadily in their rest frame. This would be the case if the minidiscs surrounding each BH act as buffers between accretion rate variations at the edge of the minidisc and regions further inside which dominate the emission (Tanaka & Haiman, 2013; Roedig et al., 2014). In this case too, Doppler boosting may be the dominant source of variability from accreting, close MBHBs.111For equal-mass binaries the two BHs would be seen out of phase, canceling the net Doppler-boost, unless one BH outshines the other, e.g., for eccentric binaries (Muñoz & Lai, 2016).

Both accretion rate variability and the Doppler boost scenario for periodic emission from MBHBs have recently been developed to interpret the MBHB candidate PG 1302-102 (Graham et al., 2015b; D’Orazio et al., 2015a; Charisi et al., 2015; D’Orazio et al., 2015b; Kun et al., 2015), a bright $z=0.3$ quasar exhibiting nearly sinusoidal periodicity in the V-band continuum. Given the measured binary mass, period, and spectral slope of PG 1302-102, D’Orazio et al. (2015b) showed that the observed amplitude of variability in the V-band and at ultra-violet (UV) wavelengths is consistent with that expected for Doppler boosting of emission from an accretion disc around the secondary BH. Further confirmation of the Doppler-boosting model for PG 1302-102 requires continued, long-term monitoring of the system at optical and UV wavelengths.

Measurements at other wavelengths can provide additional clues to the nature of the central engine of PG 1302-102, as well as of other sub-pc MBHB candidates (Graham et al., 2015a; Charisi et al., 2016; Liu et al., 2015; Zheng et al., 2016; Li et al., 2016). Here we focus on time-domain observations of PG 1302-102 in the infrared (IR). Jun et al. (2015, hereafter J15) have recently analyzed data from the Wide-field Infrared Survey Explorer (WISE) satellite to report a periodicity in the IR continuum of PG 1302-102 which is consistent with the optical period, but with a phase lag of $335\pm 153$ days in the W1 band (3.4 $\mu$m) and $524\pm 148$ days in the W2 band (4.6 $\mu$m). J15 attribute this phase lag to reverberation of the optical/UV continuum of PG 1302-102 by a surrounding dusty torus at $\sim$pc distances from the illuminating source. Another interesting feature – not discussed by J15 – is the diminished amplitude ($\sim 8\%$) of the IR fluctuations relative to the optical/UV variability ($\sim 13\%$).

In this work, we develop a toy model to interpret the findings of J15 and, in general, reverberated IR emission from a periodically variable central optical/UV source. While dust echoes of variable AGN have been considered in the past (e.g. Barvainis, 1992; Hönig & Kishimoto, 2011; Koshida et al., 2014; Hönig et al., 2017), the expected periodic nature of the brightness fluctuations of MBHBs, on tractable timescales of years, makes these systems unique. We here examine reverberation of an isotropically “pulsating” source, mimicking expectations from fluctuating accretion onto a MBHB. We contrast this scenario with reverberation of an anisotropically varying source, as expected from Doppler-boosted emission from a binary along its orbit. Such a source illuminates the surrounding dust similarly to a forward-beamed lighthouse, sweeping around at the binary’s orbital frequency. To our knowledge, reverberation from such a “rotating light-house” has not been previously explored. While in this paper we focus on dust echoes in the context of binary AGN, our models can be applied to other reverberating systems with a periodic, Doppler-boosted central engine.

We employ simple geometric models of a dusty torus that is optically thick in the optical/UV bands, but optically thin to its own IR dust emission, and is centered on the emitting source. We take into account the relative light travel time to different parts of the torus, and identify the ratio of dust light crossing time to the source variability period as the key parameter, setting the amplitude and phase of the reverberated periodic IR light curve.

We enumerate several differences between the echoes of anisotropic sources from their isotropic analogues. We find, for example, that the phase lags of the IR light curves in the two cases differ by a quarter cycle. Most importantly, for the Doppler boosted emission models, we find that, depending on the relative inclination angles of the binary’s orbital plane and of the symmetry plane of the dusty torus, variability can be present in both optical and IR, or in either one band but not the other. This means that present searches for MBHB candidates in optical surveys may have missed some candidates, and motivates extending such a search to include IR bands.

As an example of the utility of our models, we apply them to interpret the IR emission from the MBHB candidate PG 1302-102. We find that the shorter wavelength (\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip 1.0pt\sim}\hss}\raise 1.0pt\hbox{<}}2.8\rightarrow 5.3\mum) IR emission from PG 1302-102, for which long term time series data exists, is consistent with reverberation of either a Doppler-boosted or an isotropically varying central source.

While there are not enough epochs in the longer wavelength W3 ($12\mu$m) and W4 ($22\mu$m) WISE bands to fit reverberated light curve models, we use them to fit an IR spectral energy distribution (SED), allowing an estimate of the total reverberated IR luminosity, which places additional constraints on the dust radius and opening angle. Including the W3 band in the SED fitting, we find that the Doppler case is disfavoured and only nearly face-on dust tori illuminated by an isotropically varying source are allowed. We find that the bright emission in the W4 band is inconsistent with dust reverberation by the central source, and must be produced by another mechanism, such as cooler dust farther from the nucleus.

In the cases consistent with periodic dust reverberation, the dusty torus has to be thin; the size of the torus differs depending on the SED fit but has to be between 1-5 pc.

If further monitoring of longer W3 and W4 band emission from PG 13012-102 can confirm that W3 emission is dominated by reverberation by the central source, then the Doppler-boosted source may be ruled out. Otherwise, the cases of an isotropic and Doppler-boosted source are indistinguishable with current data. We find two main differences between the predictions of each. First, the Doppler-boosted scenario excludes dust tori that are within $\sim 55$ degrees of face on, while the isotropic scenario can allow face-on dust tori. However, the angular size of the reverberated-IR emitting region in PG 1302-102 is much too small for direct IR imaging to test this prediction. Second, each scenario predicts different dust torus radii. Both continued monitoring of the IR and optical/UV light curves and measurement of the reverberated-IR spectral energy distribution will narrow the constraints on the dust torus size, allowing differentiation between the Doppler boost and isotropic models.

The rest of this paper is organized as follows. We proceed in §2 by describing the MBHB and dust system. In §3 we develop models for IR emission from a dust region heated by periodically-varying isotropic and Doppler-boosted MBHB continua. In §4 we present analytic models to explore the parameter dependencies and their consequences, differentiating between isotropic and Doppler-boost scenarios. In §5 we apply our models to the MBHB candidate PG 1302-102. In §6 we summarize our findings and the limitations and possible extensions of the simple model developed here, and end by offering our conclusions.

2 Model Setup

2.1 Dust Torus

Parsec-scale dust structures surrounding the central engines of active galactic nuclei (AGN) are ubiquitous. The unification of type I (unobscured) and type II (obscured) AGN posits that the difference in AGN types is the viewing angle relative to an obscuring dust torus (Antonucci, 1993; Krolik & Begelman, 1988). This dust, heated by the central AGN source, emits in the IR and its properties have been investigated by high-resolution IR imaging as well as by modeling the IR SED.

High resolution IR imaging puts an upper limit of a few parsecs on the size of the emitting dust region (see, e.g. Elitzur, 2006, and references therein), while a wide variety of models of optically thin and thick, spatially smooth and clumpy dust distributions of various dust geometries and compositions have been employed to model the IR SEDs of AGN (see, e.g. the review by Netzer, 2015, and references therein). Neither imaging nor SED fitting, however, uniquely determines the dust properties or geometry.

We do not attempt to reproduce the full SEDs of AGN dust tori. Instead, our goal is to understand the basic properties of IR dust echoes that results from heating by a periodically variable optical/UV continuum, using a model as simple as possible. To facilitate this goal, we assume that the dust structure is centered on the illuminating source and absorbs all incident optical/UV radiation in a thin inner shell. We assume that the dust is optically thin to its own IR emission and consists of a species with a single grain size.

2.2 Central MBHB Source

The central source of optical/UV continuum which heats the dust and causes it to emit in IR is accretion-induced emission from a compact MBHB. For the expected BH masses ranging from $10^{6}\rightarrow 10^{10}M_{\rm{\odot}}$, the thermal emission from a steady-state accretion disc (Shakura & Sunyaev, 1973) has a modified blackbody spectrum that peaks in the X-rays (lower-mass BHs) to the optical (higher-mass BHs). This radiation can be efficiently absorbed to heat $\mu$m size dust particles (see §3.1).

Here we are interested in the effect of periodically variable emission from the MBHB. The simplest case is that of an isotropic, sinusoidally pulsating source,

[TABLE]

where $F^{0}_{\nu}$ is the average specific flux, $A$ is the amplitude and $P=2\pi/\Omega$ the period of the modulation, and the initial phase is taken to be zero. We note that, due to limb-darkening and geometric foreshortening, emission from an accretion disk is likely not spatially isotropic (e.g. Namekata & Umemura, 2016, and references therein), however, we employ the isotropic case as a toy model for binary accretion- induced variability, and more importantly, this choice serves as a control with which to compare to the Doppler-boosted case, which exhibits a distinct, time-dependent emission anisotropy.

Relativistic Doppler boost modulates the observed flux as

[TABLE]

where we assume that the emission is intrinsically steady and emanates from a single BH (the smaller, faster BH in the pair; this assumption is justified for PG 1302-102 in D’Orazio et al. 2015b), $c$ is the speed of light, $v^{\rm{obs}}_{||}$ is the projection of the source velocity along the observer’s line of sight, $\gamma=\left[1-(v_{s}/c)^{2}\right]^{-1/2}$ is the Lorentz factor of the secondary which orbits at speed $v_{s}=(1+q)^{-1}\sqrt{GM/a}$ for binary with mass $M$ and separation $a$,222Not to be confused with dust grain radius $a_{\rm{eff}}$ introduced below. and we assume throughout, for simplicity, that the binary is on a circular orbit. The variability period is the binary orbital period

[TABLE]

Eq. (2) further assumes that the rest-frame specific flux is a power-law $F^{0}_{\nu}\propto\nu^{\alpha_{\nu}}$ with slope $\alpha_{\nu}\equiv d\rm{ln}F_{\nu}/d\rm{ln}\nu$ at the observed frequency.

If the source is observed over all frequencies, as is approximately the case for dust absorbing a spectrum of optical/UV emission, then we emphasize here that the form of the source spectrum is not relevant. Using the Lorentz invariant quantity $F_{\nu}/\nu^{3}$ to integrate over all frequencies, the total observed flux becomes

[TABLE]

where, as long as we integrate over the entire source spectrum, we must have $\alpha_{\rm{bol}}=-1$. However, we introduce the $\alpha_{\rm{bol}}$ notation because we will also use $\alpha_{\rm{bol}}=4$ for the purposes of analytic approximation below.

Throughout we use the $\alpha_{\nu}$ notation when considering “in-band” observations of Doppler-modulated optical/UV emission, and the $\alpha_{\rm{bol}}$ notation when considering the spectrum-integrated emission seen by the dust.

Figure 1 illustrates combinations of MBHB periods and masses for which the secondary’s orbital velocity causes a significant modulation in the observed light curve. The figure assumes $\alpha_{\nu}=1$, an edge-on view of the binary, and a mass ratio of $q=0.05$, and shows that MBHBs with orbital periods of a few years and total masses of $\geq 10^{8}M_{\rm{\odot}}$ can exhibit $\geq 0.1$ mag modulations due to Doppler boosting.

3 Calculations of the Dust Echo

3.1 Reverberation of isotropically variable emission

3.1.1 Spherical dust shell

We start with the simplest case of a hollow spherical shell of dust with the source located at its center. We adopt spherical coordinates ($r,\theta,\phi$) with the source located at the origin, the observer at $(d,\pi/2,0)$ and the shell at $(R_{\rm{d}},\theta,\phi)$. Assuming that the dust is in radiative equilibrium with the heating source, we find the dust temperature as a function of time by equating the power absorbed by a dust grain to that radiated,

[TABLE]

Here $a_{\rm eff}$ is the effective (spherical) grain radius which describes the dust cross section for absorption and also the surface area for emission, $\pi B_{\nu}$ is the blackbody flux from a uniformly emitting dust grain at temperature $T_{d}$, $Q_{\nu}$, is the efficiency of dust absorption/emission and $\bar{Q}^{\rm src}$ is its average over the source spectrum.

Radiation with wavelength \lambda\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip 1.0pt\sim}\hss}\raise 1.0pt\hbox{<}}2\pi a_{\rm{eff}} is absorbed efficiently by dust grains. For longer wavelength radiation, grains of the same size become transparent. Hence for the absorption/emission efficiency we choose $Q_{\nu}=1$ for frequencies above a cutoff $\nu_{0}\sim c(2\pi a_{\rm{eff}})^{-1}$ and a power law fall off in efficiency for lower frequency radiation (e.g. Barvainis, 1987), $Q_{\nu}\equiv\rm{min}\left[(\nu/\nu_{0})^{k},1\right]$ where $k\geq 0$. 333 For this form of $Q_{\nu}$, the right-hand side of Eq. (5) can be written in terms of polylogarithmic functions. Draine & Lee (1984) fit $k=1.6$ for graphite grains. Because the efficiency for absorption is unity for high frequency radiation, above $\sim$ $1\mu$m, we take $\bar{Q}^{\rm src}=1$ throughout.

The observed flux from a single dust grain at temperature T is

[TABLE]

where $\theta_{s}=\theta_{c}$ is the angle subtended on the sky by a grain with radius $a_{\rm{eff}}$ at a distance $d$ from the observer. From the grain number density and the time-dependent dust temperature everywhere in the shell (Eq. 5), we compute the total observed specific flux summing over all dust grains as

[TABLE]

Here $\Sigma_{d}$ is the surface number density of the dust shell; $\Sigma_{d}\rightarrow\pi^{-1}a^{-2}_{\rm{eff}}$ in the assumed limit that all optical/UV radiation is absorbed by the sphere. We have evaluated the temperature at the retarded time $t_{\rm{em}}$; light emitted at the front of the dust shell at time $t$ and from the location ($R_{\rm{d}}$, $\theta$, $\phi$) at the earlier time $t_{\rm{em}}$ reach the observer at the same time at $t+(d-R_{\rm{d}})/c$. The top panel of Figure 2 illustrates this by drawing cross sections of the ellipsoids of constant light travel time (described by the equation for $t_{\rm{em}}$).

The total flux observed by an instrument with bandpass function $W(\nu)$ is

[TABLE]

where we assume for simplicity that $W(\nu)$ is a top-hat function with frequency limits $\nu_{\rm min}$ and $\nu_{\rm max}$.

3.1.2 Torus shell

The spherical dust shell of the previous section is described by a single geometric parameter, its radius $R_{d}$. We expand upon the spherical model by cutting out portions of the sphere to make a torus (or disc) which emits only from its inner edge. This torus model introduces a second and a third dust geometry parameter: the opening angle of the torus $\theta_{T}$ and the inclination of the torus to the plane perpendicular to the line of sight $J$. (Note that the opening angle sets the geometrical dust covering fraction $\cos{\theta_{T}}$.) These angles and the binary’s orbital inclination angle $I$ are drawn schematically in the bottom panel of Figure 2. To obtain the torus model, we simply set the emission to zero in the empty cone $J\pm\theta_{T}$ (and $[J+\pi]\pm\theta_{T}$) as shown in this illustration.

3.2 The Lighthouse: anisotropic, Doppler-boosted emission

To compute the effects of Doppler boost on the IR light curve we need to change the form of the source (optical/UV) flux. We use Eq. (2) for this purpose, with one important distinction: here the line-of-sight velocity $v_{||}$ is the line-of-sight speed of the secondary BH as observed by a dust grain at unit position $\mathbf{\hat{r}}_{\rm{dust}}$ in the dust shell. Written in barycentric coordinates $(r,\theta,\phi)$,

[TABLE]

where for brevity of notation we have written the binary period in terms of the angular orbital frequency $\Omega=2\pi/P$, $\phi_{0}$ is the $\phi$ coordinate of the secondary at the reference time $t=0$, and we have parameterized the secondary’s orbital velocity as $\beta\equiv a(1+q)^{-1}\Omega/c$, which depends on the binary mass ratio $q$, total mass $M$, and period $P$ through the binary orbital frequency $\Omega$ and separation $a$.

Because the emission is anisotropic, the dust temperature varies across the surface of the shell, and the analogue of Eq. (5) becomes

[TABLE]

To evaluate the LHS, we emphasize the assumption $\bar{Q}^{\rm{src}}_{\nu}=1$, that the dust absorbs all of the light emitted by the central source, or practically, that the majority of source emission is in the optical and UV. Then we may use Eq. (4) to write

[TABLE]

where $L^{0}$ is the bolometric source luminosity and our result is independent of the shape of the source spectrum. Hence, the Doppler case introduces only two new source parameters $\beta,$ and $I$ (taking the place of $A$ in the isotropic case). Once the dust temperature is found from these implicit equations, it can be used to find $F_{\nu}$ in the same way as derived for an isotropic source in §3.1 (in either the spherical shell or the torus geometry).

3.3 Model Parameters

The parameters of our model and their fiducial values are summarized in Table 1 for the convenience of the reader. These fiducial values are used in calculating the dust echo light curves below, unless stated in the text otherwise. The parameters are divided into two categories, describing the central MBHB source and the dust torus, respectively.

4 Analysis

We now identify the effect of the model parameters (see Table 1) on the IR light curves. For comparisons of the optical/UV continuum light curves with their reverberated IR counterparts, we focus on the average brightness, phase, and variability amplitude of the periodic IR emission relative to the optical/UV continuum.

In the general case, we compute IR light curves by numerically solving for the dust temperature in Eqs. (5) or (11) and then evaluating Eq. (8) for the total in-band flux. We begin by building up intuition for reverberation of periodic sources by analytically evaluating a few simplified cases.

4.1 Spherical dust shell

In the case of an optically thick dust sphere, none of the optical/UV continuum reaches the observer, but we begin by looking at this case, as it is the simplest. This can be regarded mainly as an academic exercise. However, the spherical case can be of practical interest, describing systems for which only periodic IR emission is observed. It can also be applicable if the dust sphere is only marginally optically thick, and/or if the dust distribution is clumpy on small scales, making the dust sphere porous (or if the sphere has a rare narrow ‘hole’ toward the observer).

4.1.1 Isotropically variable source

We first consider isotropic, sinusoidal emission by the central source. The IR luminosity evaluated at the retarded time $t_{\rm{em}}$ (Eq. 3.1.1), integrated over all frequencies, becomes

[TABLE]

Here $\rm{sinc}(x)=\sin{x}/x$ is the cardinal sine function, and we assumed the optical/UV luminosity of Eq. (1),

[TABLE]

with average luminosity $L^{0}$ and modulation amplitude $A$.444This analytic result can be generalized to arbitrary periodic functions by replacing $L(t)$ with a Fourier series expansion.

This simple expression in Eq. (12) gives some insight into the reverberation of a periodic continuum. We find, as expected, that the average luminosity in the IR is the same as in the optical/UV. A novel result is that the amplitude of modulation is diminished by a factor of

[TABLE]

where $A_{\rm{IR}}$ and $A$ are the amplitudes of the IR and optical/UV modulations, respectively, and we have defined the light-crossing time across the dust torus as $t_{d}\equiv R_{d}/c$. Additionally, the IR is modulated at the same period as the optical/UV continuum, but with a phase lag given by

[TABLE]

written in fractions of a cycle. When the IR fluctuation amplitude (Eq. 14) is positive, we recover the expected phase lag given by the light travel time from the central source to the dust shell. However, if $A_{\rm{IR}}<0$, there is an additional half-cycle phase change. We stress that this half-cycle phase shift occurs only for optical/UV light-curves that are reflection symmetric about their average.

Eq. (14) shows that the IR modulation amplitude is determined by the ratio $t_{d}/P$. For $t_{d}/P\rightarrow 0$, the light travel time across the dust shell is insignificant and the IR light-curve tracks the optical/UV with the same amplitude. As $t_{d}/P$ increases, the ratio $A_{\rm IR}/A$ decreases, falling to zero for $t_{d}/P\gg 1$, and whenever the light-crossing time of the dust shell is an integer multiple of the variability period, where $t_{d}/P=\frac{m}{2}$. The analytic result from Eq. (14) is shown in Figure 3 (black curve), along with the numerical results based on the equations in §3 (black crosses).

4.1.2 Doppler-modulated source

Recall from §3.2, that the dust sees a source flux modulated by the Doppler factor to the $4^{\rm{th}}$ power, i.e., $\alpha_{\rm{bol}}=-1$. Unfortunately, for this case, we are unable to find an exact analytic expression for the IR echo. However, we do find an exact analytic solution when $\alpha_{\rm{bol}}=4$, and we find that this exact solution is a good approximation to the better justified $\alpha_{\rm{bol}}=-1$ case (the two cases converge in the limit that the binary and the dust torus are face on). Here we present the exact $\alpha_{\rm{bol}}=4$ solution, and compare it to the numerical solutions for $\alpha_{\rm{bol}}=-1$ – justifying the use of the analytic $\alpha_{\rm{bol}}=4$ results.

When $\alpha_{\rm{bol}}=4$, the exact solution for the Doppler IR echo is

[TABLE]

and hereafter we assume that the binary orbital velocity is only mildly relativistic and allow the Lorentz factor $\gamma\rightarrow 1$.

The reverberated IR variation follows $\cos\left[\Omega\left(t-t_{d}\right)\right]$, in contrast to the optical/UV continuum (Eq. 3.2) which follows $\sin\left(\Omega t\right)$. This means that the IR echo has a phase lag of

[TABLE]

i.e. it is a quarter cycle out of phase with the isotropic case. This can be understood in analogy with the isotropic case, where the IR echo is delayed by $\Delta t=R_{d}/c$, which is the average light travel time difference between the front and back of the dust sphere, or equivalently the travel time between the front and halfway to the back ($\phi=\pi/2$ and $t_{\rm{em}}=t-R_{d}/c$ in Figure 2). This is also the case for the Doppler source, however, the Doppler flux seen by dust at $\phi=\pi/2$ is additionally one quarter of a cycle out of phase from the Doppler flux seen by the observer at $\phi=0$.

The relative amplitudes also differ from the isotropic case,

[TABLE]

where we have replaced A of the isotropic case by a Doppler analogue $\beta\cos{I}$. It is interesting to note that in the $\alpha_{\rm{bol}}=4$ case, $A^{\rm{Dop}}_{\rm{IR}}$ is the negative derivative, with respect to $2\pi t_{d}/P$, of $A^{\rm{Iso}}_{\rm{IR}}$. We are unaware of this being anything more than a mathematical coincidence.

Figure 3 compares analytic (Eq. 18) and numerical (§3 with $\alpha_{\rm{bol}}=4$) relative IR-variability amplitudes for a Doppler-boosted source to those of the isotropic source. Differently from the isotropic case, the modulation amplitude falls to zero for $t_{d}/P\rightarrow 0$. This is because the Doppler-boosted emission is observer-dependent, and emanating from a steady rest-frame source; in the limit that light travel time across the dust is insignificant, conservation of energy requires that the total emission integrated over a sphere does not vary in time. Observed IR variability only arises because different light travel times from different parts of the sphere cause the observer to see different cross sections of the dust sphere, in the plane perpendicular to the observer, heated at different look-back times (see Fig. 2).

The IR echo of the Doppler-boosted source becomes steady ($A_{\rm IR}=0$) at discrete values of $t_{d}/P$ that are offset by approximately a quarter cycle compared to the isotropic case, and given by

[TABLE]

The non-trivial solutions can be approximated by

[TABLE]

which converge to the condition $t_{d}\approx\left(\frac{m}{2}+\frac{1}{4}\right)P$ for $t_{d}\gg P$. The largest IR amplitudes occur for

[TABLE]

with the first three solutions at $t_{d}/P\approx 0.33,0.95,1.46$ and subsequent solutions approaching $t_{d}/P=m/2$, $m=4,5,6...$, coinciding with the zeros of the isotropic case.555This follows because of the derivative relation between the Doppler and Isotropic amplitudes (see discussion below Eq. 18).

Finally, in Figure 4 we show that the properties of the simple analytic solution discussed above captures the salient features of the more relevant “bolometric” $\alpha_{\rm{bol}}=-1$ case. One difference between the $\alpha_{\rm{bol}}=-1$ and $\alpha_{\rm{bol}}=4$ solutions is that the former depends on the binary inclination. The figure shows numerical solutions (as black, red, and blue stars) for both $\alpha_{\rm{bol}}=-1$ and $\alpha_{\rm{bol}}=4$. The exact analytic solution for $\alpha_{\rm{bol}}=4$ is shown by the solid black curve. For $\alpha_{\rm{bol}}=-1$, we plot the numerical solutions for two extremes, an edge-on (red) and a face-on (blue) binary. In the case of the face-on binary ($I=\pi/2$) we find relative amplitudes identical to the exact $\alpha_{\rm{bol}}=4$ solution. For the edge-on binary ($I=0$) the solutions for $\alpha_{\rm{bol}}=4$ and $\alpha_{\rm{bol}}=-1$ predict nearly the same $t_{d}/P$ for which zero and maximum relative amplitude occur, but the simple analytic $\alpha_{\rm{bol}}=4$ solution tends to over-predict the maximum relative amplitudes somewhat. In summary, Figure 4 shows that we may use the exact $\alpha_{\rm{bol}}=4$ solutions as a very good estimate of the maximum relative IR amplitudes from a Doppler-boosted binary. Finally, Figure 4 also plots approximate analytic solutions for $\alpha_{\rm{bol}}=-1$ (as dashed blue/red curves; see the Appendix A for details). 666 While it does not affect the analysis of the relative IR amplitude, we note that the phases of IR echoes with $\alpha_{\nu}<3$ are one half of a cycle out of phase with those with $\alpha_{\nu}>3$.

4.1.3 Absorption/emission efficiency

In the limit of our analytic solutions for the IR luminosity, where we integrate the dust emission over all frequencies, the dust absorption/emission efficiency does not affect our result (see Eqs. 5 and 11). However, when integrating over a specific IR band, the variation in dust temperature shifts the dust spectrum blue-ward and red-ward over a variability cycle. This introduces an additional filter-dependent IR variability, which depends on the absorption/emission efficiency of the dust – specifically on the frequency at which the cutoff in efficiency occurs relative to the observing band.

Figure 5 shows the effects of the adopted dust absorption/emission efficiency and finite frequency band. This figure plots the analytic solution and compares it with numerical solutions that integrate over a narrow range of frequencies from $2.8\mu$m to $4.0\mu$m (the WISE W1 band), with different values of the cutoff frequency $\nu_{0}$ in the absorption/emission efficiency (see Table 1). The small temperature changes are not enough to shift the dust spectrum across the observing band and boost the IR variability. Instead, we see that the effect of narrowing the observing band is to decrease the relative IR variability amplitude. This is likely because hotter dust emits less in the chosen IR observing band, decreasing the peak IR emission corresponding to when the central source is brightest, buffering the IR variability amplitude.

4.1.4 Light Curves

Numerically evaluating the expressions in §3, we next plot the optical/UV (source) and IR (echo) light curves for the spherical case in Figure 6. The top (bottom) panel of this figure assumes an isotropic (Doppler-boosted) source. We include the dust absorption/emission efficiency and integrate over all frequencies. To compute the IR Doppler light curves below, we use the effective ’broad-band’ $\alpha_{\rm{bol}}=-1$, while for the optical/UV Doppler light curves, we use the fiducial $\alpha_{\nu}=1$. This means that to compare relative variability amplitudes from a plot like Figure 6 to the predicted values for the Doppler case, as in Eq. 18, one must divide the IR amplitude by the optical amplitude and multiply by a factor of $(3-1)/(3-(-1))=1/2$ – this factor is a “bolometric correction” for the Doppler modulation amplitude from the optical-only band (measured by Astronomers on Earth) to the full optical/UV spectrum (seen by the dust).

The model parameters and their fiducial values are given in Table 1.

The IR amplitudes of the light curves in Figure 6 are in agreement with the predictions from Figure 3. Each IR light curve is computed for a dust sphere with radius given in units of $R_{0}=1$pc, as labeled in the figure legend. For the choice of $P=R_{0}/c$ this gives $t_{d}/P=0.8$ for the orange curve, $t_{d}/P=1$ for green, and $t_{d}/P=4/3$ for purple.

For the isotropic case, we confirm that the IR light curves lag the optical/UV continuum by the fraction of a cycle given in Eq. (15). Recall that for $A_{\rm{IR}}/A$ of different signs, the corresponding light curves are half a cycle out of phase. This means that the orange curve, for which $A_{\rm{IR}}/A<0$, is $0.5+0.8=1.3$ cycles behind the optical/UV (shifted to the right in Figure 6), while the purple curve, for which $A_{\rm{IR}}/A>0$, is $4/3$ cycles behind.

Comparison of the top and bottom panels of Figure 6 shows the predicted $1/4$-cycle lag between the isotropic and Doppler IR light curves (accounting for the half-cycle phase shifts discussed above).

Figure 7 shows the importance of the binary inclination in the Doppler case. The optical/UV and IR light curves both exhibit zero modulation amplitude as the binary approaches a face-on configuration. For the optical/UV this is because $v_{||}\rightarrow 0$ for a face-on binary. For the IR emission this is because, once we have averaged over all azimuthal angles (directions perpendicular to the line of sight), there is no time variation in dust temperature along the line of sight. Note that there will be line-of-sight variations in this azimuthally-averaged dust temperature even for a face-on binary if the dust distribution is not symmetric about the plane perpendicular to the line of sight and containing the source (e.g. for a misaligned torus; as discussed in the next section).

4.2 Torus

We next extend the results of the previous section to a torus geometry, i.e. by excising regions of the spherical shell with $\theta\leq\theta_{T}$. When the observer is looking down the axis of the torus ($J=\pi/2$: see Figure 2), we find

[TABLE]

[TABLE]

and the relative amplitudes are

[TABLE]

[TABLE]

Because the dust distribution is centered around the source, the phase lags are identical to those in the spherical shell case (Eqs. 15 and 17).

4.2.1 Isotropically variable source

Figure 8 explores the isotropic, torus solutions. The leftmost panel shows the IR to optical/UV amplitude ratios for different opening angles, for face on tori (analytic solutions as solid curves and numerical evaluations shown as crosses). The middle and right panels of Figure 8 show numerical solutions for tori inclined by $J=\pi/4$ and for $J=0$ (edge-on).

The effect of increasing the torus opening angle $\theta_{T}$ is two-fold. First, it decreases the total IR luminosity – this is not relevant for the fractional variability amplitude, but is important for the absolute IR luminosities. Second, it moves the location of the zeros of the IR amplitude to larger $t_{d}/P$. This is because the IR variability is nullified, in the isotropic case, when an integer number of variability periods matches the light crossing time along the line-of-sight dust structure. Depending on the tilt $J$ of the torus, $\theta_{T}$ changes this line-of-sight extent, and hence the values of $t_{d}/P$ for which $A_{\rm IR}=0$.

Specifically, for a face-on torus ($J=\pi/2$), the $\theta_{T}$-dependent line- of-sight extent of the dust shell is $2R_{d}\cos{\theta_{T}}$ (this can be discerned from Eq. (3.1.1) and visualized in Figure 2). As the torus is tilted away from $J=\pi/2$, the relationship between the closest and furthest points of the sphere along the line of sight becomes less dependent on $\theta_{T}$. This can be observed by comparing the left and middle panels of Figure 8.

To bracket the dependence on $\theta_{T}$, we consider the extreme cases of thin dust rings ($\theta_{T}\rightarrow\pi/2$). These are shown as green curves in Figure 8. Because $\lim_{x\to 0}\rm{sinc}{(x)}=1$, Eq. (14) tells us that the limit of a face-on ring ($\theta_{T}\rightarrow\pi/2$, $J=\pi/2$), recovers the optical/UV amplitude at all $t_{d}/P$, but at lower IR luminosity set by the covering fraction of the thin ring. This is simply because light-travel time effects are no longer important for a face-on ring (left panel of Figure 8); as $\theta_{T}\rightarrow\pi/2$ the $A_{\rm{IR}}/A$ curves stretch out further to the right at the expense of lower IR luminosity. Hence for a face-on ($J=\pi/2$) torus, the limiting behaviour is set by the black curve for a dust sphere, and a line at $A_{\rm{IR}}/A=1$ for a thin (zero-luminosity) dust ring.

In the limit of a thin, edge-on ring ($\theta_{T}\rightarrow\pi/2$, $J=0$) the solution for reverberated isotropic emission becomes

[TABLE]

where $J_{0}(z)$ is the zeroth-order Bessel function of the first kind, and the solution is valid for $a_{\rm{eff}}\ll R_{d}$. This solution is plotted in the right $J=0$ panel of Figure 8 and shows that, for a $J=0$ torus, the possible amplitudes are bracketed again by the two analytic solutions for a sphere and an edge-on ring. We plot the numerical solution for a tilted dust ring with $J=\pi/4$ in the middle panel of Figure 8, and find again that the dust ring (which suffers the least amount of averaging/cancellations) exhibits the largest relative IR amplitudes.

To investigate further the dependence on the torus inclination angle $J$, the left panel of Figure 9 shows $A_{\rm{IR}}/A$ for different inclinations, for the limiting case of a thin dust ring. We find that as the inclination increases from edge-on ($J=0$) to face-on ($J=\pi/2$), nearly the full range $0<A_{\rm{IR}}/A<1$ becomes available for all t_{d}/P\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip 1.0pt\sim}\hss}\raise 1.0pt\hbox{>}}0.4. In practice, the viable parameter space will be set by the total IR luminosity which places a lower limit on the thickness (i.e. covering fraction) of the dust ring.

4.2.2 Doppler-modulated source

Figure 10 explores the dependence of IR variability amplitude on dust geometry for a Doppler-boosted source with an edge-on binary inclination. We find similar $\theta_{T}$ dependence as in the isotropic case; for the face-on torus (left panel) the location of $A_{\rm{IR}}/A$ zeros stretches out along the $t_{d}/P$ axis as $\theta_{T}$ increases. However, contrary to the isotropic case, the limit of a $\theta_{T}\rightarrow\pi/2$ is $A_{\rm{IR}}/A=0$, rather than $A_{\rm{IR}}/A=1$. In both cases, this is a reflection of the intrinsic source behaviour in the limit that light-travel time effects become negligible, as discussed above.

As for the isotropically varying source, Fig. 10 shows the limiting cases of thin dust rings in the Doppler case (green lines). For a Doppler-boosted source (again using $\alpha_{\rm{bol}}=4$ following the discussion in §4.1.2), the edge-on dust ring produces the IR echo

[TABLE]

where $J_{1}(z)$ is the first-order Bessel function of the first kind.777Note that $-A^{\rm{Dop}}_{\rm{IR}}$ is the derivate of $A^{\rm{Iso}}_{\rm{IR}}$ with respect to $t_{d}/P$. This solution is plotted in the right ($J=0$) panel of Figure 10. The IR echo for a dust ring with an inclination of $J=\pi/4$ is shown in the middle panel of the same figure. Appendix A derives approximate edge-on ring solutions for the $\alpha_{\rm{bol}}=-1$ case, finding that the $\alpha_{\rm{bol}}=4$ solution here is adequate for discerning the behavior of IR dust echoes.

We explore the dependence on dust inclination further for a Doppler-boosted source in the right panel of Fig. 9. As in the isotropic case, we find that the main effect of increasing $J$ (towards face-on) is to stretch out $A^{\rm{Dop}}_{\rm{IR}}(t_{d}/P)$ in $t_{d}/P$. Figure 9 shows that even the purple $J=7\pi/16$ curve eventually reaches the same maximum relative amplitude as the $J=\pi/2$ curve ($\sim 0.6$), but at $t_{d}/P\sim 1.5$ instead of $t_{d}/P\sim 0.3$. This is simply because the closer to face on, the further away the dust ring needs to be to retain the same difference in light-travel times. Specifically, for a dust ring tilted by the angle $\delta\phi$ away from face on, the difference in light-travel time from the bottom vs the top of the ring is $\delta t=R_{d}/c\sin{\Delta\phi}$. So $R_{d}$ must be made large enough for $t_{d}/P\sin{\Delta\phi}\sim 0.3$, the value where the Doppler-boosted IR amplitudes reach a maximum in the $J=0$ case.

In summary, the largest IR variability amplitudes are realized for the thinnest, ring-like dust distributions. If the dust ring can be oriented at any inclination to the line of sight, then for isotropic illumination, the reverberated IR amplitude can be any fraction ($\leq 1$) of the optical/UV modulation amplitude for t_{d}/P\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip 1.0pt\sim}\hss}\raise 1.0pt\hbox{>}}0.4, and is restricted to be between one and a non-zero minimum value for t_{d}/P\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip 1.0pt\sim}\hss}\raise 1.0pt\hbox{<}}0.4. For the same orientable dust ring illuminated by an edge-on Doppler-boosted binary source, the IR variability amplitude has a maximum value of $\approx 0.6$ for t_{d}/P\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip 1.0pt\sim}\hss}\raise 1.0pt\hbox{>}}0.3, and smaller maximum values below t_{d}/P\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip 1.0pt\sim}\hss}\raise 1.0pt\hbox{<}}0.3. For a face-on binary inclination however the relative IR amplitude approaches infinity because the optical/UV variability drops to zero while the IR variability does not (see next subsection). The overall IR flux, and/or any other independent constraint on the thickness, size, and orientation of the dust structure will restrict the allowed IR variability amplitudes from those presented here for the thin dust rings.

4.2.3 Light curves

In Figure 11 we show the IR echoes (solid lines) and optical/UV light-curves (dashed lines) for various torus opening and inclination angles, choosing a fixed value of $t_{d}/P=0.6$. We recover the amplitudes shown in Figure 8 and 10 and observe the expected dimming of the IR light curves for larger torus opening angles. In both the isotropic and Doppler cases, the phase lag of the IR echo is independent of the dust geometry parameters $J$ and $\theta_{T}$. In the isotropic cases (left panels in Fig. 11) the half-cycle phase shift (§4.1.1) between the $J=\pi/2$ and $J=0$ curves is consistent with the corresponding signs of $A_{\rm{IR}}/A$ in Figure 8: for $J=\pi/2$, $A_{\rm{IR}}/A<0$ for the chosen $\theta_{T}$, while for $J=0$, $A_{\rm{IR}}/A>0$.

Finally, Figure 12 explores the effects of binary inclination in the Doppler torus model. With the freedom to orient the binary plane relative to a non-spherically symmetric ($\theta_{T}\neq 0$) dust structure through parameters $I$ and $J$, the possibility of generating IR variability with no observed optical/UV variability arises. Figure 12 demonstrates that when the binary is face on, there is no observed optical/UV variability (e.g. Eq. 3.2), but the IR variability persists.

For the case of a Doppler-boosted source, surrounded by a spherical dust shell, a face-on binary generates no IR variability because there is no time-changing emission between the front and back hemispheres of the dust shell, the integrated dust emission is constant. This back-to-front symmetry is broken in the case of the torus, as long as the dust is not mirror-symmetric around a plane perpendicular to the line of sight and containing the source, i.e., $\theta_{T}\neq 0$, $J\neq 0$ and $J\neq\pi/2$.

5 Application to MBHB candidate PG 1302-102

We now apply our models to the periodically variable quasar and MBHB candidate PG 1302-102. We use data published by J15 from the W1 ($3.4\mu$m) and W2 ($4.6\mu$m) bands of the WISE instrument. Measurements in the WISE W3 ($12\mu$m) and W4 ($22\mu$m) bands only exist at two epochs so we do not use these to fit IR light curves (see further discussion below), however, these bands can still be used to constrain the spectral energy distribution (SED). Data for W1 and W2 exist at six-month intervals from 2010 to July 2015 with a three-year gap during the hibernation of the instrument between 2011 and 2014. As in J15, we supplement the WISE W1 and W2 band data with an earlier data point from the AKARI mission (Ohyama et al., 2007, see J15 for details). Each epoch of WISE measurements consists of $\sim 10$ data points taken within a few days of each other. To fit models to the WISE data, we use a single data point at each epoch with mean and variance given by the average mean and variance of each data point at that epoch (see Figure 13).

5.1 Strategy

Using our models, we aim to determine whether the IR data from PG 1302-102 is consistent with reverberation from a periodic optical/UV source, and whether such a source is varying isotropically, or via the Doppler-boost mechanism. We also aim to constrain the parameters of a putative dusty torus.

Given the optical/UV period $P$ of the source, the primary observables are the relative brightness, variability amplitude, and phase of the IR light curves. We have shown that the amplitude ratio $A_{\rm{IR}}/A$ is dependent on the timescale ratio $t_{d}/P$ and on the dust torus opening angle $\theta_{T}$ and inclination $J$; the IR phase lag $\Phi$ is dependent on the ratio $t_{d}/P$. For a Doppler-boosted source, the IR amplitude also depends on binary inclination $I$ and orbital velocity (given by the binary mass, mass ratio, and orbital period).

We measure model-independent values of $P$, $\Phi$ and $A_{\rm{IR}}/A$ directly from the optical and WISE light curves. A blackbody spectral fit to the WISE data yields a measurement of the total IR luminosity and also the fraction of the optical/UV luminosity echoed in the IR. Together these determine the size $R_{d}$ and covering fraction $\cos{\theta_{T}}$ of the dust torus, yielding $\theta_{T}$ and $t_{d}/P$. The phase lag $\Phi$ provides a second, independent measurement of $t_{d}/P$ through the dust echo models. Because of the sinusoidal nature of the light curve, the values of $t_{d}/P$ measured from phase lags can only be determined up to addition of an integer multiple of the variability period. The two measurements of $t_{d}/P$ determine whether a consistent dust-echo model is viable and, if they agree, they narrow the possible values of $t_{d}/P$. Hence we may use the inferred values of $t_{d}/P$ and $\theta_{T}$ and determine values of the torus inclination $J$ that yield the observationally determined $A_{\rm{IR}}/A$ (the J dependence of $A_{\rm{IR}}/A$ is discussed in §4.2.1 and §4.2.2). In the Doppler case we choose the binary orbital speed and inclination based on fits to the optical variability data alone (e.g. D’Orazio et al., 2015b). A consistent model yields constraints on the dusty torus opening angle, inclination, and radius.

5.2 Model-Independent Measurements

To extract the phase lag and relative amplitude of IR light curves, we fit separate sinusoids to the V-band, W1, and W2 data.888The UV data shows periodicity consistent with the optical period, but has much fewer data points (D’Orazio et al., 2015b). For all model fitting we maximize a $\chi^{2}$ likelihood using the public Markov Chain Monte Carlo code emcee (Foreman-Mackey et al., 2013). Details can be found in Appendix B.1.

We proceed by assuming that the IR light curves are caused by reverberation of the optical/UV and hence have the same period (J15). Fixing the period to the value observed for PG 1302-102 in the optical, $P=1,884\pm 88$ days, we fit sinusoids of the form $\mathcal{A}\sin{\Omega\left(t-t_{0}\right)}+\mathcal{B}$ to the IR data. The best-fit sinusoids are shown in Figure 13, along with the V-band, W1, and W2 data. The best-fit parameters are recorded in Table 2.

The best-fit sinusoid amplitudes in Table 2 give relative IR to optical amplitudes

[TABLE]

To compare to our analytic solutions we require the quantity $A_{\rm{IR}}/A$. The relevant quantity “$A$” is the variation amplitude of the flux in the frequency range that dominates the dust heating. In the Doppler case we obtain $A$ from knowledge that the dust sees a modulation proportional to the Doppler factor to the fourth power (See Eq. 4), so that the total optical/UV amplitude is related to $A_{V}$ by $A=4/(3-\alpha_{V})A_{V}$. Using $\alpha_{V}=1.1$ (D’Orazio et al., 2015b), we find that $A=2.1A_{V}$. Then averaging the W1 and W2 values we find

[TABLE]

In the isotropic case we have incomplete knowledge of the frequency dependence of source variability. The values of $A$ in the V-band ($\approx 13\%$) and the FUV ($\approx 35\%$) are known to differ significantly, but we know that the slope of the $\lambda F_{\lambda}$ spectrum of PG 1302-102 is increasing in the FUV, levels off in the NUV, and decreases in the optical (D’Orazio et al., 2015b). Hence, the peak optical/UV power is likely emitted somewhere in-between the UV and optical values. To be conservative, we assume that the FUV and optical values bracket the actual effective $A$ by taking the mean and variance of the four values $A_{W1}/A_{V}$, $A_{W2}/A_{V}$, $A_{W1}/A_{FUV}$, and $A_{W2}/A_{FUV}$,

[TABLE]

The best-fit phases relative to the V band are

[TABLE]

The phase lags measured in J15 are

[TABLE]

Our values of the phase lag are consistent with J15’s, though the mean of our W1 phase lag extends out of the range measured in J15. This discrepancy in the phase lags measurements likely emerges from the two different methods used to determine them: J15 cross-correlate the V-band and IR light curves, while we fit sinusoids to each. Our smaller error bars likely derive from imposition of a sinusoid model inherent to our method for estimating $\Phi$. To be conservative, we use smallest and largest values of the phase lag allowed by the $1\sigma$ errors from both measurements to derive a range of values for $t_{d}/P$ from Eq. (15) and (17).

5.3 Constraints from Spectral Information Only

To place further constraints on $t_{d}/P$, and to constrain the covering fraction of the dusty torus $\cos{\theta_{T}}$, we estimate the area and temperature of the emitting dust. These constraints rely only on the observed time-averaged SED, and not on the details of the time-dependent dust echoes. Throughout we assume, as is observed, that PG 1302-102 is an unobscured quasar so that the inferred optical/UV luminosity can be associated with that which heats the dust.

We fit blackbody spectra to the time-averaged fluxes (see Appendix B.2 for details) by varying the dust temperature $T_{d}$, and the total emitting area $4\pi R^{2}_{d}\cos{\theta_{T}}$ of the dust.999We have assumed $Q_{\nu}=1$, but the result does not greatly differ if we fit for the efficiency parameters $\nu_{0}$ and $k$ (see Table 1). We impose the thermal-equilibrium constraint $R_{d}=\sqrt{L/(16\pi\sigma T^{4}_{d})}$, where we estimate the bolometric luminosity as $L=\rm{BC}\times L_{V}$ from the V-band luminosity and from the range of bolometric corrections determined empirically for quasars with similar brightness, $6\leq\rm{BC}\leq 12$ (Richards et al., 2006). Then, given the best-fit value for the total area and the range of possible BC’s, we can solve for the range of possible covering fractions $\cos{\theta_{T}}$ and dust radii. Conceptually, given a best-fit temperature, the range of allowed BC’s sets the range of possible dust sphere radii, then the covering fraction is constrained to match the total observed IR luminosity.

An important source of ambiguity is that we do not know what fraction of detected IR flux in each band is attributable to reverberation of the central source luminosity – as opposed to spatially unresolved emission from unrelated off-nuclear dust heated by star-formation, or any other source of IR emission unrelated to the torus. We here take two approaches to fitting the IR SED:

We attempt to fit all four WISE bands, and find that a blackbody spectrum can be fit to the W1, W2, and W3 bands, but not W4. We are forced to assume that contributions from W4 are primarily not due to reverberation of the central source continuum. 2. 2.

We fit only the W1 and W2 band fluxes, because these are the only two bands for which measurements of variability exist, confirming the same periodicity as in the optical/UV. These two bands should therefore be associated with reverberation of the central source, but without data at more epochs, this remains unclear for the W3 band.

For case (i), using the W1, W2, and W3 fluxes, we find a spectrum with best-fit parameters $T_{d}=577\pm 3$K, $\cos{\theta_{T}}=0.08^{+0.04}_{-0.03}$, and $R_{d}=3.9\pm 1.0$ pc, corresponding to $t_{d}/P=3.1\pm 0.8$.

For case (ii), using only W1 and W2 fluxes, we find a spectrum with best-fit parameters $T_{d}=881\pm 14$K, $\cos{\theta_{T}}=0.05^{+0.02}_{-0.02}$, and $R_{d}=1.9\pm 0.4$ pc, corresponding to $t_{d}/P=1.5\pm 0.3$.

Case (i) yields a solution for dust that is colder, further out, and has a larger covering fraction than in case (ii). This is because there is more flux in the lower-wavelength W3 band, pushing case (1) to a lower temperature but a larger covering fraction (see Figure 18).

A clear way to distinguish between the above two cases, and to remove the ambiguity in the best-fit values for the reverberated IR luminosity, is to continue monitoring the longer-wavelength emission in the W3 (and W4) bands for periodicity.

Table 3 summarizes our measured values for the IR PG 1302-102 quantities. We now use them to constrain dust and source properties of PG 1302-102.

5.4 Interpretation with Dust-Echo Modeling

We now place constraints on dust torus models, based on the predicted properties of their dust echoes.

In Figure 14, we show regions allowed by various constraints in 2D parameter planes: in ($\cos{\theta_{T}}$, $t_{d}/P$) in the top row and in the ($J$, $\cos{\theta_{T}}$) plane in the bottom row. The left panels are for isotropically varying sources and the right panels are for Doppler-boosted sources.

In the top row, we show five distinct constraints:

The regions where the values of $A_{\rm{IR}}/A$ predicted by our models fall within the measured values (Table 3) are shaded blue and green, for two different fixed values of the torus inclination for which the solutions are analytic. The blue regions are for a face-on torus ($J=\pi/2$; Eq. 4.2) and the green regions are for an edge-on dust ring ($J=0$ and $\theta_{T}\approx\pi/2$; Eqs. 24 and 25). There are multiple disjoint allowed regions, because of the oscillatory patterns in $A_{\rm{IR}}/A$ as a function of $t_{d}/P$. 2. 2.

The horizontal red bands show the allowed range of the covering fraction $\cos{\theta_{T}}$ inferred from the IR SED modeling alone. This constraint is independent of $t_{d}/P$. 3. 3.

The hashed, orange regions show the values of the ratio $t_{d}/P$ allowed by fitting the observed phase lags using the dust echo models (Eqs. 15 and 17). The disjointed nature of these orange regions is from the ambiguity in adding an integer number of periods to the phase lag measurements, and also from the half-cycle phase change occurring when $A_{\rm{IR}}/A$ changes sign. Because the sign of $A_{\rm{IR}}/A$ depends on $\cos{\theta_{T}}$ and $J$, we hash the orange regions with the color corresponding to the value of $J$ for which they are calculated: blue ($J=\pi/2$) and green ($J=0$). 4. 4.

The vertical solid orange band spanning the entire y-axis is from the independent measurement of $t_{d}=R_{d}/c$ from the blackbody fits to the average WISE-band fluxes (see Table 3), combined with $P$ measured directly from the sinusoidal fits to the light-curves. 5. 5.

By assuming that the source emits at $\epsilon$ of the Eddington luminosity, we estimate the sublimation radius of dust as a function of black hole mass (using Eq. 2 of Mor & Netzer 2012). Assuming that dust tori cannot exist inside this radius, we obtain the lower limit

[TABLE]

The grey shaded regions to the left are excluded by this constraint (in Figure 14 we use the PG 1302 bolometric luminosity to determine the sublimation region).

The top panels shows that there exist consistent models for periodically reverberated emission from PG 1302-102 which satisfy all five constraints. Graphically these are realised where the red, solid orange, blue-hashed orange, and blue ($J=\pi/2$) or red, solid orange, green-hashed orange and green ($J=0$) regions overlap, and avoid the grey regions. However, these overlaps do not express the full space of solutions, they only take into account specific dust torus models for which we have analytic solutions. To show the full space of allowed solutions, we next numerically calculate $A_{\rm{IR}}/A$ over the range of joint constraints on $t_{d}/P$ and $\cos{\theta_{T}}$ and for the full range of torus inclinations angles $J$.

The regions of allowed values of $t_{d}/P$ and $\cos{\theta_{T}}$, for which we now calculate $A_{\rm{IR}}/A$ as a function of $J$, are highlighted in yellow in the top panels of Figure 14. The bottom panels of Figure 14 calculate $A_{\rm{IR}}/A$ in only these yellow regions, for a grid of torus inclinations spanning from [math] to $\pi/2$ and the allowed range of $\cos{\theta_{T}}$. For each pair of $J$ and $\cos{\theta_{T}}$, we ask if any value of $t_{d}/P$ in the allowed range generates $A_{\rm{IR}}/A$ in the measured range; if so, we assign the region in the bottom panels of Figure 14 the color white, if not that region is colored black. We then draw a black contour excluding the black regions. In this way we discern the limits on the torus inclination in each of our models. Note the we also take into account that each yellow region, in progressing $t_{d}/P$, is allowed only for a specific sign of $A_{\rm{IR}}/A$ (see Eqs. 15 and 17), and this sign alternates between yellow regions (see “+” and “-” labels in the top rows of Figures 14 and 15).

For computational simplicity, to calculate $A_{\rm{IR}}/A$ in the bottom panels of Figure 14, we allow $Q_{\nu}=1$ everywhere and compute the IR variability amplitude when integrating over all frequencies. As Figure 5 shows, this generates a small overestimate for the relative IR amplitude and results in a conservative estimate of the ruled out dust parameters (overestimate of the size of white regions in Figures 14). Practically, we use a $10\times 10$ grid of $J$ and $\cos{\theta_{T}}$, sampling each yellow region uniformly with five values of $t_{d}/P$. With these approximations in mind, we use only the general trends in the bottom panels of Figure 14 to estimate the regions of excluded torus inclination angle.

Figure 15 illustrates the same restrictions on parameter space as Figure 14, but for the constraints found by including only the W1 and W2 band average fluxes in the SED fitting (section II, case (ii) in Table 3).

Figures 14 and 15 reveal that the short-wavelength (\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip 1.0pt\sim}\hss}\raise 1.0pt\hbox{<}}20\mum) IR emission from PG 1302-102 is consistent with reverberation by a periodically variable source. However, only the case which includes W1 and W2 bands is consistent with all constraints assuming an isotropically varying or a Doppler-boosted central source. Case (i) which also includes SED information from the W3 band is consistent only with an isotropically varying source, for a small range of nearly face-on dust tori. The allowed parameter values depend on whether the WISE W3 band is excluded or included in the SED fitting, along with the W1 and W2 bands. Common to both cases is the restriction to a small dust covering fraction (requiring the dusty torus to resemble a thin disc). This constraint arises from fitting the time-averaged WISE SEDs, independently of the time-dependent dust echo modeling.

When we include the W1, W2, and W3 bands in the SED fit (section II, case (i) of Table 3 and Figure 14) we find a dust covering fraction, or equivalently a disc aspect ratio, of 0.05\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip 1.0pt\sim}\hss}\raise 1.0pt\hbox{<}}\cos{\theta_{T}}\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip 1.0pt\sim}\hss}\raise 1.0pt\hbox{<}}0.12, and inner radii of $\sim 2.9\rightarrow 4.9$ pc, but for different disjoint values indicated by the ranges of $t_{d}/P$ listed in §III of Table 3 for isotropic and Doppler-boost models. The bottom panels of Figure 14 show us that, in the isotropic case, the inclination of the dusty disc must be J\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip 1.0pt\sim}\hss}\raise 1.0pt\hbox{>}}0.8 radians ($57^{\circ}$) away from edge on, and may be face on. In the Doppler case, there are no allowed dust torus solutions within our model.

That there are consistent solutions for periodically reverberated emission from an isotropic source, but not a Doppler-boosted source in case (i) may simply be a manifestation of the greater uncertainty on the value of $A_{\rm{IR}}/A$ for the isotropic source (see § 5.2). A better determination of variability amplitudes across the spectrum of PG 1302 will help to resolve this issue.

Including only the bands for which there is evidence for periodic variability, the W1 and W2 WISE bands (section II, case (ii) of Table 3 and Figure 15) we find dusty discs with aspect ratio 0.03\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip 1.0pt\sim}\hss}\raise 1.0pt\hbox{<}}\cos{\theta_{T}}\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip 1.0pt\sim}\hss}\raise 1.0pt\hbox{<}}0.07 and and inner radii of $\sim 1.5\rightarrow 2.3$ pc, but for different disjoint values indicated by the ranges of $t_{d}/P$ listed in §III of Table 3. The bottom panels of Figure 15 show us that, in the isotropic case, the inclination of the dusty disc is largely unconstrained but cannot be face on, J\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip 1.0pt\sim}\hss}\raise 1.0pt\hbox{<}}0.6 ($34^{\circ}$) or 1.1\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip 1.0pt\sim}\hss}\raise 1.0pt\hbox{<}}J\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip 1.0pt\sim}\hss}\raise 1.0pt\hbox{<}}1.5 radians ($63^{\circ}\rightarrow 86^{\circ}$), this is also clear in the top left panel of Figure 15, where the analytic face-on torus solutions, represented by the blue region, do not intersect the yellow highlighted regions. The situation is similar for the Doppler case, the dust inclination is limited to J\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip 1.0pt\sim}\hss}\raise 1.0pt\hbox{<}}0.6 radians ($34^{\circ}$) from edge on.

In both case (i) and case (ii), the Doppler solutions are restricted from having face-on dust inclination angles because a thin, face-on dust disc can not generate large enough reverberated IR amplitudes. This follows because the dust is constrained to be in a thin disc (torus with a large opening angle). When such a disc is oriented nearly face on, the time delay along the line of sight dust structure becomes small and, as discussed in §4.2.2, this results in the relative amplitude of IR variability to also become small for a dust distribution that is symmetric about the origin.

For the isotropically pulsating source, we find that case (i) rules out dust discs that approach edge-on inclinations, while case (ii) rules out those approaching face on. This difference is caused mainly by the different required locations of the dust radius in either case.

In (isotropic) case (i), a more-edge-on ring cannot generate a large enough relative IR amplitude at the required, large values of $t_{d}/P$. A face-on dust torus, however, can generate large relative IR amplitudes for arbitrarily large $t_{d}/P$, as long as the dust torus is thin enough, ($\cos{\theta_{T}}$ is small; see left panel of Figure 8). In (isotropic) case (i), the relative IR amplitude even becomes too large at small $\cos{\theta_{T}}$ as evidenced by a change from allowed to not-allowed face-on dust tori with decreasing $\cos{\theta_{T}}$ (see the top of the bottom left panel of Figure 14).

In (isotropic) case (ii), the smaller radius of the dust torus, and hence, smaller required values of $t_{d}/P$, cause relative IR amplitudes to be too large in the case of a thin, face-on dust torus. At the same time, the relative IR amplitudes generated by an edge-on ring are large enough to meet the measured values of $A_{\rm{IR}}/A$ at these smaller values of $t_{d}/P$.

While the case (ii) spectral fit is more restrictive on the size of the inner dust region $R_{d}$, it is less restrictive in dust inclination angle $J$. This is because it is easier to match the required values of $A_{\rm{IR}}/A$ for the smaller values of $t_{d}/P$ in case (ii).

We have so far ignored the relative motion of the binary with respect to the dusty torus. This motion can introduce additional time-dependence in the light-travel times, with significance on the order of the ratio of binary separation to the size of the IR reprocessing region (the dust). This ratio is

[TABLE]

where we have again estimated the dust radius by the sublimation radius (Eq. 2 of Mor & Netzer 2012) with $\epsilon$ the Eddington fraction. This tells us that the impact of binary orbital motion on dust reverberation models is most important for the lowest mass binaries with the longest periods, and contributes at most at the percent level for the fiducial values taken here for a PG 1302-102-like binary. Because the dependence on mass is weak, and 5 years is a present upper limit on observed binary periods discovered in EM time-domain surveys (e.g. Graham et al., 2015a), it is safe to assume that binary orbital motion is a \mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip 1.0pt\sim}\hss}\raise 1.0pt\hbox{<}}2\% effect. It will nevertheless be important to take this effect into account for longer-period binaries, with orbits still surrounded by a single dust region, binaries emitting at a small fraction of Eddington, and also for modeling reverberation by the closer-in broad-line regions around MBHBs (the subject of future work).

Finally, we reiterate that the constraints for thin dust tori around PG 1302-102 arise from two independent lines of reasoning. The first is based only on the time-averaged SED data and the geometry of the torus, and is independent of the time-dependent dust-echo modeling. The black-body emitting dust must be heated to the observed temperature and emit at the observed overall luminosity. We find that this forces the dust to be far enough from the central source to be cool enough to match the IR SED shape, and to also have a small enough covering fraction to avoid over-producing the IR luminosity (See section II of Table 3). The second requirement is based on the light-curves and the time-dependent dust-echo models. The large measured IR variability amplitudes, and our finding that thinner dusty discs generate larger relative IR amplitudes is also consistent with tori being thin, this can be seen, e.g., from the blue regions in the top panels of Figures 14 and 15: in the allowed range of $t_{d}/P$ values, $\cos{\theta_{T}}$ is restricted to, at thickest, \cos{\theta_{T}}\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip 1.0pt\sim}\hss}\raise 1.0pt\hbox{<}}0.3 in the isotropic case (ii).

6 Discussion and Conclusions

We summarize our key results and discuss their implications for MBHBs.

The phase lag of IR variability relative to optical/UV variability is given by $2\pi t_{d}/P$ radians in the isotropic case and $(1+1/4)2\pi t_{d}/P$ radians in the Doppler case. There is an additional half-cycle phase shift for light curves which are reflection symmetric (sinusoids); whether this phase shift occurs depends on the values of $t_{d}/P$ as well as the dust geometry (see Eq. (4.2)). These differences, for periodic sources, should be considered when relating an IR time lag with the light travel time across the dust reverberation region. 2. 2.

The amplitude of IR variability relative to optical/UV is a function of $t_{d}/P$, the inclination and opening angles of the dusty torus, and also the binary inclination to the line of sight for a Doppler-boosted source (see Figures 8–12). In the isotropic case, the IR amplitude falls to zero for $t_{d}/P\gg 1$, and approaches that of the optical/UV continuum for $t_{d}/P\rightarrow 0$. The Doppler case obtains peak amplitude at t_{d}/P\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip 1.0pt\sim}\hss}\raise 1.0pt\hbox{>}}0.3, depending on the torus properties, and falls to zero for both $t_{d}/P\rightarrow 0$ and $t_{d}/P\gg 1$. The isotropic case also exhibits zero relative IR amplitude when the line-of-sight light crossing time of the dust is equal to an integer multiple the source variability period while the Doppler case exhibits zero relative IR amplitude when the dust line-of-sight light crossing time is approximately an integer plus one quarter cycle of the binary orbital period.

When the dust has the geometry of a face-on or edge-on torus, the IR variability amplitude in the Doppler-boosted case does not exceed $\sim 60\%$ of the optical/UV amplitude (See Figures 10 and 9). As we showed in Figure 12, this is not the case when the torus is misaligned and the binary is nearly face on. We discuss the implications of this in the next bullet point.

While Figure 9 shows that thin, nearly-face-on dust rings can generate large IR variability at any value of $t_{d}/P$ (though in a trade-off for total IR luminosity), we find that for the majority of dust and binary parameters, the largest IR variability amplitudes occurs for smaller values of $t_{d}/P$. Hence we use $t_{d}/P$ as an indicator for IR variability amplitude in Figure 1. In Figure 1, we overlay contours of $t_{d}/P$ to illustrates MBHB parameters conducive to generating significant reverberated IR variability.

We conclude from Figure 1 that, for isotropic sources, long-period optical/UV modulations by emission from low-mass binaries have the potential to generate the largest relative IR modulations. The situation is different in the Doppler case where t_{d}/P\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip 1.0pt\sim}\hss}\raise 1.0pt\hbox{<}}0.1 results in small relative IR variability amplitude. In addition to smaller relative amplitudes, small $t_{d}/P$ also implies weak absolute optical/UV variability amplitudes due to the Doppler boost and hence small absolute IR amplitude as well. Hence, in the Doppler case, either intermediate masses and binary periods, or the thin nearly-face-on dust rings discussed above are favored for detection of reverberated IR variability. 3. 3.

Orphan-IR variability can occur for Doppler-boosted sources which are nearly face on (so we do not see the optical/UV variability), but are surrounded by a dusty torus that is not symmetric across the plane separating the front and back of the Doppler-boosted source. Periodic orphan-IR variability in quasars would be a signature of Doppler-boosted MBHBs at high binary inclination to the line of sight, or to a variety of variable sources completely surrounded by dust. These periodic sources would not yet have been identified in optical searches such as those carried out by Graham et al. (2015a) and Charisi et al. (2016). 4. 4.

Periodic IR variability from PG 1302-102 is consistent with reverberation of a Doppler-boosted source surrounded by a thin disc of dust with inclination within \mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip 1.0pt\sim}\hss}\raise 1.0pt\hbox{<}}34^{\circ} from the line of sight. Consistent solutions for periodic reverberation also exist for an isotropically varying source surrounded by a thin dust disc at any inclination of the disc to the line of sight. Each scenario requires a dusty disc aspect ratio $\sim 0.1$ and inner radius between 1 and 5 pc. The isotropic and Doppler-boost cases each require different inner dust radii in this range, allowing a test to differentiate the source nature if more complete IR SEDs of reverberated emission and phase lag measurements can narrow down the true size.

Because of the small inferred angular size of the reverberated-IR emitting region (\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip 1.0pt\sim}\hss}\raise 1.0pt\hbox{<}} milli-arcsecond), it is not currently possible to image this region in the IR at the distance of PG 1302-102. Though it would be worth determining if imaging of regions further out in the dust structure could lend clues to the dust geometry and orientation at smaller radii.

These scenarios also require that longer wavelength (\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip 1.0pt\sim}\hss}\raise 1.0pt\hbox{>}}22\mum for isotropic source and \mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip 1.0pt\sim}\hss}\raise 1.0pt\hbox{>}}12\mum for a Doppler-boosted source) IR emission from PG 1302-102 is not from reverberation of the central MBHB emission, and so our models predict that IR emission in the WISE W4 and longer wavelength bands will not exhibit periodic variability at the optical period. A clear determination of which IR wavelengths exhibit such periodicity, and are thus due to central source reverberation, will help to constrain the dust+source model for PG 1302-102. Specifically, determination that the W1 through W3 band ($\sim 3\rightarrow 17\mu$m) emission is primarily due to reverberation by the central source while W4 band emission ($22\mu$m) is not, would add evidence against the Doppler-boost nature of the central periodic source of optical/UV emission. The enhanced constraints on dust parameters from narrowing down the reverberated IR SED can be seen by comparing the constraints from the two possible SED fits in Table 3, of which we quote only the weaker, joint constraints above. Further significant improvements would come from measurements of the frequency dependent variability amplitude of PG 1302-102. These would allow a tighter constraint on the relative IR to optical/UV amplitude in the isotropic case.

6.1 Caveats and extensions

• •

We have considered smooth, single-species dust models which are optically thin to their own emission and optically thick to optical/UV emission, and are at a single temperature at each location. Future work should consider dust models which include dust heating by optical/UV over a finite radial extent, IR absorption by dust, clumpy dust (e.g Netzer, 2015, and references therein), and a distribution of dust grain species and sizes, as well as a finite range of temperatures.

• •

The geometry of the dust torus may not be a simple conical cut of a sphere, but may be convex in nature so that optical/UV light could be absorbed by the dust at a range of radii even under the assumption here that the dust is optically thick to optical/UV (e.g. Pozo Nuñez et al., 2014). This could add a colder IR component with significantly different time lag structure and should be considered in future work.

• •

We have considered a dust distribution which is always centered on the emitting source, and reflection symmetric about the origin. Deviation from this symmetry will alter our predictions. For example, the Doppler boost IR amplitudes will not go exactly to zero for $t_{d}/P\rightarrow 0$ if dust is only in front of or only behind the source.

• •

We have considered the case of isotropically variable emission as a proxy for variations due to accretion changes and also as a control with which to compare the Doppler-boosted source models. However, emission from an accretion disc itself is likely spatially anisotropic, emitting less from the disk edges (e.g Namekata & Umemura, 2016, and references therein). Unlike the Doppler-boost case, however, this anisotropy is constant in time. Addition of this effect would introduce another parameter to both Doppler-boost and non-Doppler-boost models: the relative inclination of accretion disc and dust torus.

• •

We assume the binary is on a circular orbit. Some hydrodynamical models of the binary interaction with a gas disc predict the excitation of large binary eccentricities (e.g. Cuadra et al., 2009; Roedig et al., 2011). Binary eccentricity could change the shape and periodicity structure of the optical and hence the IR light curves predicted here (although we note that the nearly sinusoidal light curve in the example PG 1302-102 shows no evidence for eccentricity; D’Orazio et al. 2015b).

• •

If grains can re-form on a timescale shorter than a binary orbital time, the sublimation radius will change periodically with the changing central source flux. The change in dust sublimation due to the changing observed luminosity variations $\delta L$ is $2\delta R_{d}/R_{d}=\delta L/L$. For typical Doppler luminosity variations, this corresponds to a change in the inner dust radius of a few to $\sim 10\%$. In the Doppler case, the largest variations occur for less massive binaries with shorter periods. In the isotropic case, only the amplitude of modulation is relevant. In either case, the dust will emit at a constant (source-frame) temperature, but light travel time lags will acquire extra time dependence.

• •

General relativistic effects such as time delay, lensing, and precession could become important and affect optical/UV continuum and reverberated IR emission differently.

• •

Our predictions for reverberated periodic continuum hold true for any source with a periodic component. Addition of any non-periodic component such as noise, does not changes our findings, but rather adds an additional IR signature on top of the one presented here. We expect that quasi-periodic signatures will create reverberated IR light curves that can be treated as a perturbation from those of the periodic case.

6.2 Conclusion

We have developed models to compute IR light curves that result from the reverberation of a periodic optical/UV source by a surrounding dusty torus. For the first time, we have considered the reverberation of periodic optical/UV sources, with either isotropic luminosity variations, or anisotropic source variability induced by the relativistic Doppler boost. The latter setup resembles a rotating, forward-beamed light-house surrounded by fog. We found a number of differences from the case of reverberation of an isotropically variable source, and applied our findings to the IR emission from the MBHB candidate PG 1302-102. We found that dust tori are constrained to be thin, and that the present data already places non-trivial constraints on the inclination angle of the torus in the Doppler case, requiring the torus to avoid being nearly face on.

The type of models described here can be applied to fit the IR and optical/UV light curves of the growing list of MBHB candidates (Graham et al., 2015a; Charisi et al., 2016; Jun et al., 2015), aiding in their interpretation, and could also aid in identifying new MBHB candidates.

Acknowledgments

We thank Hyunsung Jun for providing the WISE data from his analysis, and Hyunsung Jun, Daniel Stern, Arlin Crotts, Jules Halpern, Andrew MacFadyen, Jeffrey J. Andrews, Adrian Price-Whelan, and James Guillochon for useful discussions. DJD thanks Janna Levin and Pioneer Works for providing a stimulating work space during the completion of this work. The authors additionally thank the anonymous referee for his/her comments. This publication makes use of data products from WISE, a joint project of the University of California, Los Angeles, and the Jet Propulsion Laboratory/California Institute of Technology, funded by NASA. This publication makes use of data products from NEOWISE, which is a project of the Jet Propulsion Laboratory/California Institute of Technology. NEOWISE is funded by the National Aeronautics and Space Administration. Financial support was provided from NASA through Einstein Postdoctoral Fellowship award number PF6-170151 (DJD) and NASA ATP grants NNX11AE05G and NNX15AB19G (ZH). ZH also gratefully acknowledges support from a Simons Fellowship for Theoretical Physics.

Appendix A Analytic Approximation to Doppler relative IR amplitudes

Throughout the text we approximate the Doppler dust echo using a solution where the illuminating source has an effective spectral slope of $\alpha_{\rm{bol}}=4$. This is because this allows a simple, exact analytic formula for the light curve and hence the relative amplitude and phase, which we have shown agrees well with numerical solutions for the more relevant “bolometric” Doppler spectral slope of $\alpha_{\nu}\approx-1$ (see Fig. 4).

Although not used for the calculations in this paper, here we give approximate analytic solutions for the relative amplitudes in the $\alpha_{\rm{bol}}=-1$ case, which may be useful in future studies.

For small binary orbital velocity and when the binary is face on, the approximate solutions presented here approach the $\alpha_{\rm{bol}}=4$ solutions used in the main text.

For a face on torus ($J=\pi/2$), we wish to evaluate the integral

[TABLE]

where $v_{||}$ is a function of $\phi$ and $\theta$ given by Eq. (3.2). Expanding the integrand in a binomial series to second order in $\beta$ (the orbital velocity in units of the speed of light $c$), and choosing a convenient observer position so that $t_{\rm{em}}=t-t_{d}\left(1-\cos{\theta}\right)$, we can solve analytically for the time dependent light curve. Because we already have the optical/UV light curve, we can find the maximum and minimum of both the IR and the UV light curves to find the amplitude of each. Taking the ratio of amplitudes, we then find

[TABLE]

While $A=a$ for the $\alpha_{\rm{bol}}=4$ case and for $\gamma\rightarrow 1$, this is not true in general. For a binary which is also face on ($I=\pi/2$) this solution simplifies considerably to

[TABLE]

which is the $\alpha_{\rm{bol}}=4$ exact solution presented in the main text to first order in $\beta$. Figure 4 of the main text compares the $\alpha_{\rm{bol}}=-1$, $\alpha_{\rm{bol}}=4$, and numerical solutions.

We also find the analog of Eq. (24) for an edge-on dust ring in the Doppler case, but with the broad-band “bolometric” $\alpha_{\rm{bol}}=-1$. Following the same procedure as above, expanding to second order in $\beta$, we find

[TABLE]

where $J_{0}(z)$ and $J_{1}(z)$ are the zeroth- and first-order Bessel functions of the first kind. Figure 16 plots the $I=0$ and $I=\pi/2$ limits of this approximation compared with the numerical evaluation and the $\alpha_{\rm{bol}}=4$ analytic solution, for different values of $\beta$.

We find that

The $\alpha_{\rm{bol}}=-1$, $I=0$ analytic approximation does not match well the numerical solution, while the $\alpha_{\rm{bol}}=-1$, $I=\pi/2$ analytic solution does. 2. 2.

Still, the $\alpha_{\rm{bol}}=-1$ solutions do not differ greatly from the $\alpha_{\rm{bol}}=4$ solution unless $I\rightarrow\pi/2$. 3. 3.

The $\alpha_{\rm{bol}}=-1$, $I\rightarrow\pi/2$ solutions show new behavior. As $I\rightarrow\pi/2$, and for increasing $\beta$, the $\alpha_{\rm{bol}}=-1$ relative amplitudes blow up to infinity starting at $t_{d}/P=0$ and propagating outwards in $t_{d}/P$. This may be because the optical/UV amplitude is approaching zero as the binary approaches face on.

Based on these results, and because we are not concerned with near face-on binaries in the case of PG 1302-102, we use the exact $\alpha_{\rm{bol}}=4$ solutions presented in the main text for vetting PG 1302-102 models.

Appendix B MCMC fits

In each case below we use the emcee code (Foreman-Mackey et al., 2013) to find best-fit model parameters to the data. In the main text and in Tables 2 and 3, we quote our results in terms of a maximum likelihood value $\pm$ the values which exclude the outer $15\%$ of the likelihood.

B.1 Sinusoid fit

For each of the V-band, W1, and W2 time series we first fit a sinusoid by minimizing

[TABLE]

with respect to the three parameters $A$, $t_{0}$, and ${\rm mag_{0}}$ in each band, where $t_{i}$, mag$(t_{i})$, and $\delta$mag$(t_{i})$ are the observation times, magnitudes, and errors. We impose the flat priors $A>0$ and $0\leq t_{0}\leq P$ and fix $\Omega$ to the best-fit value previously found by Graham et al. (2015b) corresponding to a rest-frame period of $P=1473.72$ days (or observed $P_{\rm obs}=1884$ days) from $V$ band data.

Figure 13 shows the best-fit sinusoids with this fixed period. The black data points and error bars are the epoch-averaged magnitudes and standard deviations used in the fit, while the coloured points are the unbinned data. Figure 17 shows the sampling of the posterior distribution.

B.2 Blackbody SED fit

The best-fit blackbody spectral energy distribution is found for three different cases: using the time- averaged fluxes from all four WISE bands, only the fluxes from the W1, W2, and W3 bands, and also using only the W1 and W2 bands. The average WISE fluxes are, $F_{\rm{W1}}=8.83\pm 0.18$ mJy, $F_{\rm{W2}}=13.33\pm 0.26$ mJy, $F_{\rm{W3}}=36.07\pm 0.56$ mJy, and $F_{\rm{W4}}=101.859\pm 2.62$ mJy (1$\sigma$ errors). We minimize

[TABLE]

with respect to two parameters: the dust temperature $T_{d}$ and $X$, a combination of the dust covering fraction $\cos{\theta_{T}}$ and total source luminosity after the bolometric correction,

[TABLE]

where $d$ is the luminosity distance to the source and $L_{V}$ is the V-band luminosity.

The above model flux derives from the assumption that the grains are optically thin to IR. Then the flux from each grain is $\pi B_{\nu}[T_{d}](a_{\rm{eff}}/d)^{2}$ and there are $N=\Sigma_{d}\int^{2\pi}_{0}\int^{\pi-\theta_{T}}_{\theta_{T}}R^{2}_{d}\sin{\theta}\ d\phi\ d\theta=4\pi R^{2}_{d}\cos{\theta_{T}}$ grains. For an optically thick grain distribution $\Sigma\rightarrow(\pi a^{2}_{\rm{eff}})^{-1}$ and we obtain

[TABLE]

which equals $XB_{\nu}[T_{d}]/T^{4}_{d}$ upon imposing the thermal equilibrium constraint

[TABLE]

Once we fit for $T_{d}$ and $X$, we impose an empirically motivated range of the bolometric correction (Richards et al., 2006) which limits the possible range of the covering fraction that can reproduce the total observed IR luminosity. We choose BC$=9\pm 3$ to obtain the range of $\cos{\theta_{T}}$ presented in Table 3.