A Mote in Andromeda's Disk: a Misidentified Periodic AGN Behind M31

TL;DR

This paper reclassifies a star in M31 as a background AGN, analyzes its spectrum, detects potential binary black hole signatures through periodicity in photometry, and discusses implications for gravitational wave sources.

Contribution

It identifies a misclassified star as a background AGN behind M31 and presents evidence of possible binary supermassive black holes through periodic signals.

Findings

Spectroscopic analysis confirms the AGN nature of J0045+41.

Detected periodicities suggest a binary black hole system.

Potential gravitational wave source within the detection range.

Abstract

We identify an object previously thought to be a star in the disk of M31, J0045+41, as a background $z\approx0.215$ AGN seen through a low-absorption region of M31. We present moderate resolution spectroscopy of J0045+41 obtained using GMOS at Gemini-North. The spectrum contains features attributable to the host galaxy. We model the spectrum to estimate the AGN contribution, from which we estimate the luminosity and virial mass of the central engine. Residuals to our fit reveal a blue-shifted component to the broad H$\alpha$ and H$\beta$ at a relative velocity of $\sim4800$ km s$^{-1}$. We also detect \ion{Na}{1} absorption in the Milky Way restframe. We search for evidence of periodicity using $g$-band photometry from the Palomar Transient Factory and find evidence for multiple periodicities ranging from $\sim80-350$ days. Two of the detected periods are in a 1:4 ratio, which is identical to the predictions of hydrodynamical simulations of binary supermassive black hole systems. If these signals arise due to such a system, J0045+41 is well within the gravitational wave regime. We calculate the time until inspiral due to gravitational radiation, assuming reasonable values of the mass ratio of the two black holes. We discuss the implications of our findings and forthcoming work to identify other such interlopers in the light of upcoming photometric surveys such as the Zwicky Transient Facility (ZTF) or the Large Synoptic Survey Telescope (LSST) projects.

Tables

Table 1. Table 1: Results from § 4 . T 𝑇 T is the period, P S ​ ( T ) subscript 𝑃 𝑆 𝑇 P_{S}(T) is as described in the text, T m ​ i ​ n subscript 𝑇 𝑚 𝑖 𝑛 T_{min} and T m ​ a ​ x subscript 𝑇 𝑚 𝑎 𝑥 T_{max} are the bounds of the period bin in which the FAP is calculated. The last column shows whether the period is detected using a DRW + sinusoidal mean model in celerite .

$T$	$P_S\left(T\right)$	($T_{min},T_{max})$	FAP	Detected with
days		days		celerite?
82.10	0.120592	(78.809,79.280)	$6.98469\times {10}^{-3}$	Yes
117.84	0.139525	(116.428,116.898)	$7.14281\times {10}^{-3}$	Yes
162.04	0.148967	(154.047,154.517)	$7.78917\times {10}^{-3}$	No
202.01	0.212229	(201.071,201.541)	$4.72885\times {10}^{-3}$	Yes
328.03	0.233829	(323.332,323.802)	$3.30188\times {10}^{-3}$	Yes
354.84	0.270498	(351.546,352.016)	$1.01854\times {10}^{-3}$	Yes
409.85	0.248934	(407.974,408.444)	$4.27479\times {10}^{-3}$	No
702.34	0.281859	(699.520,699.990)	$4.59042\times {10}^{-3}$	Yes
867.86	0.300183	(859.400,859.870)	$5.84198\times {10}^{-3}$	No

Table 2. Table 2: Orbital and gravitational properties of proposed orbital periods

$T$	$R$/$\theta$	$v_{orb}$	$t_{GW}$	$h_0$
days	AU/$\mu$arcsec	${10}^3$ km s-1	yr
		($q=1/0.01$)	($q=1/0.01$)	($q=1/0.01$)
82.10	216.02/0.30	14.312/28.341	$3.522\times {10}^2$/$8.982\times {10}^3$	$9.252\times {10}^{-17}$/$3.628\times {10}^{-18}$
162.04	339.90/0.47	11.410/22.594	$2.159\times {10}^3$/$5.505\times {10}^4$	$5.880\times {10}^{-17}$/$2.306\times {10}^{-18}$
328.03	543.93/0.75	9.020/17.860	$1.416\times {10}^4$/$3.610\times {10}^5$	$3.674\times {10}^{-17}$/$1.441\times {10}^{-18}$

Equations

$$E=f-G\cdot c$$

$$R=E^T{\Sigma }^{-1}E$$

$$c=(G^T{\Sigma }^{-1}G{)}^{-1}(G^T{\Sigma }^{-1})\cdot f$$

$$⟨A_V⟩={\tilde{A}}_V{e}^{{\sigma }^2/2}$$

$${\sigma }_A^2={\tilde{A}}_V^2{e}^{{\sigma }^2}(e^{{\sigma }^2}-1)$$

$$PSD(T)=\frac{4{\sigma }^2\tau }{1+4\pi (\tau /T{)}^2}$$

$$S(\Delta t)={\sigma }^2{e}^{-|\Delta t|/\tau }$$

$$\ln \mathcal{L}=-\frac{1}{2}{r}^T{K}^{-1}r-\frac{1}{2}\ln |K|-C$$

$$\begin{array}[]{ccc}t_{GW}&=&\frac{5}{256}\frac{c^{5}}{G^{3}}\frac{R^{4}}{(M_{1}+M_{2})(M_{1}M_{2})}\\ &=&\frac{5}{256}\frac{c^{5}}{G^{3}}\frac{R^{4}}{M_{tot}^{3}}\frac{(1+q)^{2}}{q}\end{array}$$

$$h_{0}=\frac{4G}{c^{2}}\frac{qM_{tot}}{(1+q)^{2}D_{L}(z)}\bigg{(}\frac{2\pi GM_{tot}}{c^{3}T}\bigg{)}^{2/3}$$

Full text

galaxies: active — galaxies: individual: LGGS J004527.30+413254.3 — quasars: supermassive black holes

††journal: ApJ††software: aplpy v1.1.1 (Robitaille & Bressert, 2012), Astropy v1.2.1 (Astropy Collaboration et al., 2013), Ciao v4.7 (Fruscione et al., 2006), celerite v0.3.0 (Foreman-Mackey et al., 2017), corner v2.0.1 (Foreman-Mackey, 2016), emcee v2.2.1 (Foreman-Mackey et al., 2013), gatspy v0.3 (VanderPlas & Ivezić, 2015; VanderPlas, 2015), IRAF v2.16 (Tody, 1986), Matplotlib v2.0.0 (Hunter, 2007), NumPy v1.11.3 (Van Der Walt et al., 2011), PySynphot v0.9.8.4 (Lim et al., 2015), and Sherpa v4.7 (Freeman et al., 2001)

1 Introduction

Active Galactic Nuclei (AGN) are among the most luminous persistent sources of radiation in the Universe, capable of outshining their host galaxies when in a quasar state. They are hosts to supermassive black holes (SMBHs) and are found throughout the history of the universe from redshift $z\sim 7$ onward (Mortlock et al., 2011). With the advent of surveys like the Sloan Digital Sky Survey (SDSS, York et al. 2000), the number of cataloged AGN has increased by many orders of magnitude.

As incredibly powerful sources of ionizing radiation, AGN drive and regulate the evolution of the stars, gas, and dust of their host galaxies. The major merger of two gas-rich galaxies can trigger intense dust production and star formation, while the increased accretion onto the central black hole of one or both galaxies can increase its luminosity, triggering outflows and regulating star formation (Sanders et al., 1988), leaving behind a massive, gas-poor elliptical remnant. Such mergers appear to be not only frequent, but the primary means by which both SMBHs and galaxies are grown (Kauffmann & Haehnelt, 2000). If both galaxies in a merger contain SMBHs, simulations indicate that the black holes themselves can merge over $\sim$Gyr timescales (Volonteri et al., 2003; Tremmel et al., 2017). At early times, the SMBHs in a merger will appear as dual or offset AGN (depending on the accretion rate of both black holes, Comerford et al. 2015). As their orbits decay, the black holes can form a supermassive black hole binary (SMBHB), which could be observed as an apparently single AGN that displays periodic variability. We present here spectroscopic and time-domain analyses of an AGN behind M31 that has been previously misidentified as a red supergiant, a globular cluster, and an eclipsing binary. We find evidence for the periodic variability of the AGN and discuss the implications of its misidentification in light of forthcoming large photometric surveys.

1.1 J0045+41

As a part of a search for red supergiant X-ray binaries — a still-theoretical class of exotic stellar binary system — we used the single-epoch photometry of the Local Group Galaxy Survey (LGGS, Massey & Olsen 2003; Massey et al. 2006, 2007), which covers M31, M33, the Magellanic Clouds and 7 dwarf galaxies in the Local Group, to assemble a statistical sample of Local Group red supergiants (RSGs). We used the method of Massey (1998) to reduce contamination from the far more prevalent foreground M-dwarfs by taking advantage of the separation of the two populations in $B-V$ vs. $V-R$ space. After creating our sample (and ensuring our results agreed with Massey et al. 2009 in M31, where we find 437 candidate RSGs), we searched the Chandra Source Catalog (CSC, Evans et al. 2010) for X-ray sources within 10” of the LGGS RSGs. This search yielded one close match.

LGGS J004527.30+413254.3 ($\alpha=00^{\rm h}45^{\rm m}27^{\rm s}.30$, $\delta=+41^{o}32^{\prime}54^{\prime\prime}.31$, Figure 1), which we will refer to as J0045+41 hereafter, is a bright ($V\approx 19.9$) object of previously-unknown nature in the disk of M31. Vilardell et al. (2006) classified J0045+41 as an eclipsing binary with a period of $\sim 76$ days. While the observed variability is of order 1 magnitude in $B$ and $V$, their data are poorly sampled in phase. On the other hand, Kim et al. (2007) included J0045+41 in a catalog of candidate globular clusters, and it has also been included in catalogs of M31 globular clusters as recently as 2014 (Wang et al., 2014). The LGGS photometry was consistent with the color and brightness of a typical 12-15 $M_{\odot}$ RSG in M31, with an inferred effective temperature of $\sim$3500 K and bolometric magnitude of -6.67 (following Massey et al. 2009). However, the best SED fit to photometry from the Panchromatic Hubble Andromeda Treasury (PHAT, Dalcanton et al. 2012) using the Bayesian Extinction And Stellar Tool (BEAST, Gordon et al. 2016) is a 300 M*⊙, 105 K “star”, extincted by $A_{V}\sim$4 magnitudes, which we exclude as being unphysical. This discrepancy is likely due to the broader wavelength coverage of the PHAT dataset, as well as the fact that the BEAST performs a complete SED fit, whereas our RSG selection criteria are purely based on color and magnitude cuts to select for bright, red objects roughly consistent with the photometric properties of RSGs. Furthermore, the object appears extended in the PHAT images (though its radial profile appears similar to that of other nearby stars; see Figure 1), implying that J0045+41 may be a background AGN or quasar. Given the angular size of M31 at optical wavelengths ($\sim 10$ deg2) and the typical surface density of bright quasars on the sky ($\sim 18$ deg-2*, Richards et al. 2002), we expect $\sim 180$ sources in the entirety of M31 to actually be background AGN.

J0045+41 is separated by $\sim 1.18^{\prime\prime}$ (4.45 pc at the distance of M31) from an X-ray source in the Evans et al. (2010) catalog. This source, CXO J004527.3+413255 ($\alpha=00^{\rm h}45^{\rm m}27^{\rm s}.30$, $\delta=+41^{o}32^{\prime}55^{\prime\prime}.46$), is bright ($F_{X}=1.98\times 10^{-13}$ erg s*-1* cm*-2*) and has hardness ratios from Evans et al. (2010) that are consistent with an unabsorbed X-ray binary or AGN. To confirm this, we fit a spectrum from the publicly available ChandraPHAT dataset (Obs. ID 17010, Williams 2014) with an absorbed power law model (xstbabs * powlaw1d) in Sherpa (Freeman et al., 2001). We use the atomic cross-sections from Verner et al. (1996), and abundances from Wilms et al. (2000). The spectrum is binned to ensure each bin has a minimum of five counts, and we fit the background-subtracted spectrum from 0.3 to 8 keV. The best-fit ($\chi^{2}_{red}=0.34$) model has a neutral H column density $N_{H}=1.7\times 10^{21}$ cm*-2* and a power-law slope $\Gamma=1.5$. The spectrum and fit are shown in Figure 2. The value of $N_{H}$ derived from the fit corresponds to an extinction of $A_{V}<1$, which would be surprising if CXO J004527.3+413255 was a background AGN or quasar seen through the disk of M31, as we would expect a significantly higher column density. In addition, using the ChandraPHAT data, Williams et al. (in prep) derive improved source locations and positional errors, resulting in a much better alignment between CXO J004527.3+413255 and J0045+41 (see Figure 1)

To conclusively determine the nature of J0045+41, we decided to obtain optical spectrophotometry. We discuss our observations and data reduction in §2. We present the spectrum, use it to classify J0045+41 as an AGN, identify key features, and analyze it in §3, and search for evidence of periodicity using archival data in §4. We conclude with a discussion of our results and their implications in §5.

2 Observations and Data Reduction

We obtained a longslit spectrum of J0045+41 using the Gemini Multi-Object Spectrograph (GMOS) on Gemini-North (Hook et al., 2004). Four 875 second exposures were taken 2016 July 5 using the B600 grating centered on 5000 Å, and four 600 second exposures were taken 2016 July 9 using the R400 grating centered on 7000 Å, with a blocking filter to remove 2nd-order diffraction. Two of each set of exposures were offset by +50 Å to fill in the gaps between the three CCDs in GMOS. We followed the standard GMOS-N reduction pipeline using the gemini package in IRAF (Gemini Observatory & AURA, 2016). Flux calibration was performed using HZ 44 (Oke, 1990) as a standard star for both sets of observations. The final reduced spectrum is continuous from $\sim$4000 to $\sim$9100 Å at a resolution of $R\sim 1688\>({\rm blue})/1918\>({\rm red})$.

3 Spectrum and Analysis

The optical spectrum is shown in Figure 3. It shows the broad emission lines characteristic of an AGN. We use Ca II H & K, the Fe I/H$\gamma$/[O III] G band, [O III] $\lambda 5007$, Mg I $\lambda\lambda$5192,5197, Na I D, and He I $\lambda 7067$ to determine that J0045+41 is at $z\approx 0.215$. We also detect Na I D doublet absorption in the rest frame of the Local Group; however our data are not of sufficient resolution to distinguish Milky Way from M31 absorption. Both H$\alpha$ and H$\beta$ are broad, with full widths at half maximum of $\sim 10^{4}$ km s*-1*. The centers of broad H$\alpha$ and H$\beta$ are slightly blueshifted ($z\approx 0.21$) relative to the rest of the spectrum, which may be indicative of an outflow or motion of the central engine relative to the host galaxy.

Mistaking a blue AGN for a red star might seem unsurprising given that it is seen through the disk of M31. However, the low amount of extinction implied from the fit to the X-ray spectrum seems inconsistent with an object seen through an entire galactic disk. In Figure 4, we show our spectrum of J0045+41 compared with the composite Sloan Digital Sky Survey (SDSS, York et al. 2000) quasar template spectrum from Vanden Berk et al. (2001) as well as a template Seyfert 2 spectrum from PySynphot (a Python implementation of Synphot destributed by Space Telescope Science Institute, Lim et al. 2015), both redshifted to $z=0.215$ and reddened by 1 (top) and 2 (bottom) magnitudes of extinction in $V$ using a standard Cardelli et al. (1989) $R_{V}=3.1$ extinction law. While hardly a robust fit, this comparison serves to illustrate that either a larger value of extinction is required to reproduce the overall spectral shape of J0045+41 with a pure QSO template or that many of the spectral features — e.g., the apparent break in spectral slope at $\sim 5500$ Å and the presence of strong absorption lines in the spectrum – are intrinsic to the host galaxy.

To decompose the spectrum into host and AGN spectra, we follow Vanden Berk et al. (2006). We use the first five galaxy eigenspectra and the first ten QSO eigenspectra derived from a Principal Component Analyses (PCA) of SDSS galaxy and quasar samples (Yip et al., 2004a, b) as a set of basis spectra, which we redden using the Cardelli et al. (1989) extinction law, redshift to $z=0.215$, and fit to our spectrum of J0045+41 as follows. If the measured fluxes are represented by a column vector, $f$, then the residuals between the data and the basis spectra fit is simply

[TABLE]

where $G$ is a matrix whose columns are the redshifted and reddened basis spectra interpolated to the values of the observed wavelengths in our spectrum and $c$ is a column vector containing the coefficients for each basis spectrum. Taking the errors on each point into account, the scaled residual at each point can be represented by the scalar

[TABLE]

where $\Sigma$ is the covariance matrix and $E^{T}$ denotes the matrix transpose. It can be shown that the coefficients that minimize $R$ are given by

[TABLE]

In order to estimate a suitable value of $A_{V}$ to use when reddening the basis spectra, we redden the spectra with integer values of $0\leq A_{V}\leq 10$ mag. Some of these fits are shown in Figure 5. While the basis spectra sufficiently fit the spectrum for $0\leq A_{V}\leq 2$ mag, at higher values, the basis spectra are unable to reproduce the spectral shape, especially in the blue. Going forward, we adopt $A_{V}=1$ mag. Dalcanton et al. (2015) mapped the dust extinction in M31 at a resolution of 25 pc using the PHAT dataset. They model the probability distribution of $A_{V}$ in each pixel with a log-normal distribution, parametrized by the median extinction, $\tilde{A}_{V}$ and the dimensionless width, $\sigma$, such that the mean extinction $\langle A_{V}\rangle$ is

[TABLE]

and the variance in the extinction $\sigma_{A}^{2}$ is

[TABLE]

Dalcanton et al. (2015) also include the fraction of stars in each pixel that are reddened, $f_{red}$. In the pixel containing J0045+41, $f_{red}=0.206$, $\tilde{A}_{V}=0.72$, and $\sigma=0.28$. The latter two values correspond to $\langle A_{V}\rangle=0.75$, $\sigma_{A}=0.21$, consistent with our estimate of $A_{V}$. Spectral modeling at higher resolution would further constrain the extinction along the particular line of sight towards J0045+41.

The galaxy and AGN components of this fit are shown in the top panel of Figure 6. The bottom panel shows the dereddened rest-frame luminosity spectrum of each component. The luminosity of the underlying AGN component is $L_{\lambda}=3.46\times 10^{39}$ erg s${-1}$ Å*-1* at 5100 Å. The derived host galaxy spectrum appears similar to an early type galaxy. This is unsurprising as the hosts of low-luminosity AGN (like J0045+41) tend to be early type (Kauffmann et al., 2003). If the periodicity (discussed in Section 4) arises from a SMBH binary formed through the major merger of two late type AGN hosts, it would also be unsurprising that the resulting host is an early type galaxy.

With the underlying contribution to the spectrum from the central engine now known, it is possible to estimate the mass of the SMBH (Shen et al., 2008). We use the full width at half maximum of H$\beta$ ($1.11\times 10^{4}$ km s*-1*), the continuum rest frame luminosity from the quasar at 5100 Å, and the H$\beta$ virial mass estimator coefficients from McLure & Dunlop (2004) to calculate $\log(M/M_{\odot})=8.30$. We use the bolometric correction from Runnoe et al. (2012) to calculate the bolometric luminosity, from which we determine the Eddington ratio $\Gamma\equiv L_{bol}/L_{Edd}=0.007$. This small value for $\Gamma$ may indicate that the accretion flow is radiatively inefficient (Casse & Keppens, 2004).

4 Potential Periodicity

4.1 Searching for Periodicity Using the Supersmoother Algorithm

Though the light curve in Vilardell et al. (2006) is sparsely sampled, the suggestion of a $\sim$76 day period in J0045+41 prompted further investigation. While continuum emission from AGN is well-known to be stochastically variable due to a variety of phenomena associated with the central engine and surrounding environment, periodicities in the variability have long been predicted as a signature of SMBHBs (e.g., Bogdanović et al. 2008). A short-period SMBH system would be well within the gravitational wave regime. We investigated the reported periodicity using data from the Palomar Transient Factory (PTF, Law et al. 2009). PTF observed J0045+41 in both $g$ and $r$, though the $g$-band data cover a broader range in time, and thus we focus our analyses solely on those data. These data are shown in Figure 7.

AGN continuum variability is well fit by a damped random walk (DRW) process (Kelly et al., 2009), described by a characteristic timescale ($\tau$) and long-term rms variability ($\sigma$ or $SF_{\infty}=\sqrt{2}\sigma$). The power spectral distribution (PSD) of a DRW process (Charisi et al., 2016) is

[TABLE]

and the covariance function is

[TABLE]

where $\Delta t$ is the time between two observations.

Previous searches for periodicities in AGN lightcurves commonly use Lomb-Scargle periodograms (Liu et al. 2016; Charisi et al. 2016; Zheng et al. 2016). Lomb-Scargle periodograms detect periodicities in irregularly-sampled lightcurves by fitting sinusoids to the data (Lomb, 1976; Scargle, 1982). It is important to note that sinusoidal variability is expected if the periodicity arises due to the relativistic Doppler boost of the emission of the secondary component of a steadily-accreting binary (see D’Orazio et al. 2015). However, the predicted periodicity from SMBHBs is not necessarily sinusoidal if caused by periodic episodes of accretion (e.g., Farris et al. 2015). Furthermore, Vaughan et al. (2016) show that the behavior generated by red noise processes can be well fit by a sinusoid over a few ‘cycles’. Therefore the statistical significance of previously-reported detections using Lomb-Scargle periodogram analysis may be overestimated.

To provide a robust assessment of periodicities in the lightcurve of J0045+41, we utilize the Supersmoother algorithm (Reimann, 1994), which uses a non-parametric periodic model to test the strength of signals at various periods. Using the implementation in the gatspy Python package (VanderPlas & Ivezić, 2015), we calculate the periodogram of the $g$-band data on a linearly spaced grid of 2000 periods between 60 and 1000 days — we are unlikely to see periods shorter than 60 days (see Charisi et al. 2016), and our data do not cover more than two cycles of a signal with more than a 1000 day period. The periodogram is shown in Figure 8. As expected by a DRW signal, the power appears to rise to a constant level at long periods. However, there do appear to be real peaks superimposed onto the expected DRW behavior.

4.2 Estimating the Significance of Measured Peaks

To check that the measured power of the true signal ($P_{S}(T)$) is not attributable to a DRW process, we generate simulated DRW lightcurves, following the prescription of MacLeod et al. (2010), and compare the distribution of the periodograms of the simulated lightcurves to $P_{S}(T)$. While it is possible to calculate the DRW parameters, $\sigma$ and $\tau$, from the estimated mass of J0045+41, we choose to instead estimate those parameters by fitting the lightcurve directly, thus incorporating the distribution of possible values. We implement (7) as a kernel function in celerite (Foreman-Mackey et al., 2017), a Python package for Gaussian process computations, which calculates the likelihood, $\mathcal{L}$, of a DRW with given $\sigma$ and $\tau$:

[TABLE]

where $r$ is a vector of the observed data minus the mean, $K$ is the covariance matrix incorporating the photometric errors and the DRW covariance function, and $C$ is a constant proportional to the number of measurements (for a discussion of Gaussian processes and the derivation of this likelihood function, see Rasmussen & Williams 2006). We then use emcee (Foreman-Mackey et al., 2013), an affine-invariant MCMC Python package, to fit for $\sigma$, $\tau$, and the mean magnitude $\langle g\rangle$ by sampling the posterior distribution. We use 32 walkers, and, after discarding 500 burn-in steps, record 3000 samples per walker for a total of 96,000 samples. A corner plot of these samples is shown in Figure 9.

Drawing the value of $\sigma$, $\tau$ and $\langle g\rangle$ from the posterior distribution of samples, we generate 96,000 DRW lightcurves. The lightcurves are sampled at the same times as the PTF observations and have identical photometric errors. The final points in the simulated lightcurve are then drawn from a Gaussian distribution with the magnitude of the raw point as the mean, and standard deviation equal to the photometric error. We then calculate periodograms for each simulated DRW lightcurve on the same grid of periods as $P_{S}$. The mean ($P_{DRW}$) and standard deviation ($P_{\sigma}$) of the simulated periodograms are plotted as $P_{DRW}\pm P_{\sigma}$ along with $P_{S}$ and the theoretical DRW PSD with $\sigma=0.2$, $\tau=200$ days (scaled to match the values returned by Supersmoother) for comparison in the left panel of Figure 10. Much of the structure in the true periodogram is matched by the simulated periodograms, but not in the theoretical PSD. This is likely due to the irregular sampling of the PTF lightcurve, which is reflected in the simulated lightcurves. However, some of the peaks in the true periodogram do not appear in the DRW noise.

To identify periods with power in excess of the DRW noise, we search for peaks in $\sigma=(P_{S}-P_{DRW})/P_{\sigma}$. $\sigma(T)$ is plotted in the right panel of Figure 10, with the ten peaks with largest $\sigma$ indicated by blue triangles. As Supersmoother only returns values between 0 and 1 when it calculates the periodogram — and thus the values are not normally distributed — $\sigma$ as a statistic is meaningless by itself. We instead want to estimate the false-alarm probability (FAP) of each peak. Traditional estimates of significance (see Horne & Baliunas 1986, for example) assume that the null hypothesis is pure white noise. Because the background noise is dependent on the period, we split the grid of periods into $N_{trial}=100$ bins with 20 periods each. In each bin, we find the period $T$ associated with the largest value of $P_{S}$. We then calculate the number of simulated periodograms that have at least one point with power greater than $P_{S}(T)$ ($N_{DRW}(>P_{S}(T))$) within the period bin. The FAP is thus $N_{DRW}(>P_{S}(T))$ divided by the number of simulated DRW periodograms ($N_{DRW}=96,000$) times $N_{trial}$, which accounts for the fact that there are $N_{trial}\times N_{DRW}$ ‘chances’ to randomly generate a peak with more power than the true peak (the look-elsewhere effect).

4.3 Distinguishing Periodicity from Systematics

The above process results in a number of periods that correspond to local minima in FAP vs. $T$, shown as blue triangles in the right side of Figure 10. Between the sampling of the lightcurve and the algorithm used to generate the DRW lightcurves, it is possible that some of these detections are only arising due to artificial suppression of the DRW noise. To determine this, we use the same algorithm to simulate white noise lightcurves ($\tau\rightarrow 0$, with $\sigma$ and $\langle g\rangle$ drawn from the DRW samples in Figure 9), and calculate the average ($P_{WN}$) and standard deviation ($P_{\sigma,WN}$) of the periodograms. $P_{WN}\pm P_{\sigma,WN}$ is shown in purple on the left panel of Figure 10. It is clear that $P_{WN}$ and $P_{\sigma,WN}$ are roughly constant over the range of tested periods, and thus that none of the detected periodicities arise due to suppression of the DRW noise.

It is also possible that the period detected at $T=354.8$ days is due to approximately yearly systematic variations in observing conditions — e.g., airmass, observability, weather, etc. — at Palomar Observatory, and that the period at $T=708.5\approx 2\times 354.8$ is an alias of the same effects. This appears to be reflected in Figure 12, where the phase-sampling of both the $g$ and $r$ band data is nearly identical at these periods. Because J0045+41 is nearly at the detection limit of PTF, it is certainly possible that those systematics can masquerade as real effects; our discussion of these results comes with the major caveat that the $\sim$yearly periodicity may not be real. However, even discounting the $354.8$-day period, there is a secondary peak at $328$ days that is unlikely to be a result of these yearly systematics.

Finally, if these periods are real, they should be detectable by other means. We add a sinusoidal mean model to our implementation of the DRW kernel within celerite, and simultaneously sample the posterior distribution of the model parameters — mean, amplitude, period, and phase — and the DRW parameters as described above using emcee, using double the number of walkers, and restricting the period of the sinusoid to lie between 60 and 1000 days. As discussed above, a sinusoidal model is not necessarily an accurate one; however, the periods revealed by this analysis should be similar to the periods found above. A histogram of the posterior distribution of the period is shown in Figure 11, with the periods with local minima in FAP indicated by blue triangles. It is clear that at least some of the peaks found — namely at $T=82.1,117.8,202.0,328.0,354.8,\>{\rm and},708.3$ days — are retrieved. The phase-folded, mean-subtracted data and the best-fit Supersmoother model at the six periods detected with celerite, along with the phase-folded $r$-band are shown in Figure 12. Table 1 contains the period $T$, the value of $P_{S}(T)$, the bounds of the period bin containing $T$ ($T_{min}$ and $T_{max}$), the estimated FAP, and whether a strong peak in the celerite posterior appears at a similar period.

The period of $\sim 82.1$ days (FAP$\sim 0.007$) is similar to Vilardell et al. (2006) who find a period of $\sim 76$ days. We plot the PTF data, the historical data from Vilardell et al. (2006) (offset by a constant for clarity), and the best-fit Supersmoother model folded on the period found by Vilardell et al. (2006) in Figure 13. None of the structure in the Vilardell et al. (2006) data is seen in the PTF data or the Supersmoother fit; however, with so few observations, it is possible that the true period detected by Vilardell et al. (2006) is closer to that detected in the PTF data. Unfortunately, the historical data are only available phase-folded, and we are unable to include them in our analysis of other periods.

5 Discussion and Conclusion

One possible interpretation of a periodic signal in an AGN is that it is due to the orbital motions of a SMBHB, formed through a major galaxy merger. Though small, the number of $z<1$ candidate SMBHBs discovered is consistent with this model (Volonteri et al., 2009). The detected periodicities of J0045+41 are thus quite interesting. Most intriguingly, the $\sim 82.1$ day period is almost exactly in a 1:4 ratio with the $\sim 328$-day period. It is possible that either of these peaks is an alias of the other, as the observed periodogram is the convolution of the true periodogram with the Fourier transform of the sampling function (Roberts et al., 1987; Charisi et al., 2015). However, multiple periodicities beyond the orbital period are predicted to occur in SMBHBs at similar period ratios as a result of interactions in the circumbinary disk (MacFadyen & Milosavljević, 2008; Shi et al., 2012; Farris et al., 2014).

In particular, MacFadyen & Milosavljević (2008) found that the periodogram of the accretion rate in their simulation displayed significant peaks at frequencies approximately generated by the formula $\omega=\frac{2}{9}K\Omega_{bin}$, where $\Omega_{bin}$ is the binary orbital angular frequency, and $K=1,2,6,7,8,9,10$. We search for the orbital period, $T_{bin}$, that generates a set of periods closest to the first five observed periods (discounting the 354.8 and 708.5 day periods). We find that $T_{bin}=169.29$ ($\Omega_{bin}=3.7\times 10^{-2}$ day*-1*) creates periods that match quite well with the two shortest periods, though it underpredicts the 202 day period by $\sim 75$ days, and overpredicts the 328 day period by $\sim 50$ days. Finally, Farris et al. (2014) find that, for varying mass ratios and simulation setups, periodic variations in the accretion rate onto one or both black holes can arise at frequencies with the same 1:4 correspondence as the 82 and 328 day periods. These occur at $1/4\Omega_{bin}$ and $\Omega_{bin}$. This points to the 82.1 day period being the orbital period of the binary. Farris et al. (2014) also find frequencies arising at $\Omega_{bin}$ and $2\Omega_{bin}$. Interestingly, we do detect a period with FAP$\sim 0.008$ at $162\approx 2\times 82\approx\frac{1}{2}\times 328$ days. While we do not detect a strong peak in the celerite posterior around this period, this hints that the orbital period may also be 162 or 328 days.

If we assume that any of these three periods is the orbital period of a SMBHB in a circular Keplerian orbit, and that the virial mass derived in §3 is the total mass of the two black holes $M_{tot}$, then the semimajor axis of the orbit ranges from 216 to 544 AU (or 0.3 to 1 microarcseconds at the angular diameter distance of J0045+41, which is unresolvable using current radio interferometric arrays). Such a separation would be well within the regime where loss due to gravitational radiation is significant. We can approximate the time for two circularly orbiting black holes to inspiral due to gravitational radiation using equations (5.9) and (5.10) from Peters (1964):

[TABLE]

where $R$ is the semimajor axis of the orbit, $M_{1},\>M_{2}$ are the masses of the individual black holes and $q\equiv M_{2}/M_{1}$. $t_{GW}$ ranges between $\sim 350$ yr (for the shortest period, with $q=1$) to 360 kyr (for the longest period, with $q=0.01$).

The gravitational waves produced by SMBHBs are expected to be detectable at the nHz frequencies probed by pulsar timing arrays (PTAs, Foster & Backer 1990). The amplitude of the dimensionless gravitational strain ($h_{0}$) of a SMBHB with mass ratio $q$ at redshift $z$, assuming a circular orbit with period $T$ can be expressed as

[TABLE]

where $D_{L}(z)$ is the luminosity distance (Thorne, 1987). The expected strain of a SMBHB with the derived mass and orbital period of J0045+41 would range from $\sim 10^{-16}$ (for the shortest detected period, with $q=1$) to $\sim 10^{-18}$ (for the longest period, with $q=0.01$). These results, in addition to the expected orbital velocity of the secondary black hole (see below) are summarized in Table 2. While the latter strain would be orders of magnitude below the stochastic background of gravitational radiation from all SMBHBs at that period ($h\approx 10^{15}$ at $T=1$ yr, Shannon et al. 2013), the background falls off at higher frequencies as fewer sources are expected to be inspiraling at shorter and shorter periods, and the signal from a $\sim 80$ day SMBHB would be detectable above the background (Moore et al., 2015). Indeed, the signal would be just shy of the anticipated sensitivity — $\sim 6\times 10^{-16}$ (Lazio, 2013) — of the Square Kilometer Array (SKA, Dewdney et al. 2009). While this is an exciting finding, it is important to note that there are a number of other possible interpretations of a periodic signal, e.g: a long-lived or periodically-generated hot spot in the accretion disk, geodetic precession, and self-warping of the disk (see Bon et al. 2017 for a concise review).

Even if it is not a SMBHB, J0045+41 is an interesting object. For one, it appears to be probing a relatively extinction-free region of the ISM in M31. The detection of the Na I D doublet is promising, and follow-up optical and infrared observations at higher spectral resolution may disentangle absorption from M31 and from the Milky Way, and reveal more about the dynamics of the ISM along the line of sight towards J0045+41. The spectrum is well fit by a mixture of the galaxy and quasar eigenspectra from Yip et al. (2004a, b) redshifted to $z=0.215$ and reddened by an $A_{V}=1.0\pm 1.0$ mag Cardelli et al. (1989) extinction law. However, H$\alpha$ and H$\beta$ both have a blueshifted broad component. Indeed, the residuals to the fit shown in Figure 6 appear to be Gaussian. Fitting these residuals with a Gaussian profile shows that this component is at $z=0.196$, a $\sim 4800$ km s*-1* difference from the host redshift. This shift may be due to an outflow from the central engine, a hot spot in the blueshifted side of the accretion disk, or the blending of the broad lines of each SMBH component; as the less massive SMBH moves towards us, we would see its broad lines blueshifted, which would explain the excess of blue flux in the broad lines (Shen & Loeb, 2010). Indeed, a similar binary model has been used to explain SDSS J092712.65+294344.0, which also appears to have blueshifted broad lines relative to the narrow lines in the spectrum (Dotti et al., 2009; Bogdanović et al., 2008). At the short periods found in §4, orbital velocities are expected to be $\sim 10^{4}$ km s*-1* (depending on the assumed mass ratio), so this blueshift would be consistent with the orbital velocities for all of the periods in Table 2, for any value of the mass ratio. Follow-up spectroscopy on a cadence of a few months would be able to search for or exclude periodic changes of the H$\alpha$ and H$\beta$ profiles relative to the narrow lines over time, which would help point to an explanation.

To search for any objects similar to J0045+41 in color space, we used PySynphot (a Python implementation of Synphot distributed by Space Telescope Science Institute, Lim et al. 2015) to generate synthetic photometry from our spectrum in $g$, $r$, $i$, and $z$ — there was not enough signal in $u$ to synthesize a magnitude. We then downloaded photometry of all low-redshift ($z<1$) SDSS quasars from Data Release 13 (SDSS Collaboration et al., 2016) within 0.1 magnitudes of J0045+41 in $g-r$ vs. $r-i$ vs. $i-z$ color space. These quasars are shown in color space in Figure 14. Each point is colored by the assumed value of the extinction in $g$. Of these 446 objects, only 197 of them have redshifts that are positive — implying the remaining objects are not plausibly quasars. Indeed, the spectra of many of the ‘quasars’ in this sample are quite clearly cool stars. Some of these objects are simply misidentified; however, many are flagged with a Z_WARNING: NOT_QSO by the SDSS pipeline. While this is helpful for reducing contamination of the quasar sample, it illustrates than many objects of interest fall through the cracks of classification algorithms (see Dorn-Wallenstein & Levesque 2017 for further discussion).

Of the true quasars in the sample, none are extincted by more than 1.5 magnitudes in $g$. It is likely that these quasars (and the AGN component of J0045+41) are intrinsically red as described by Richards et al. (2003). These quasars may have been reddened by dust intrinsic to the host galaxy, or have excess red flux due to synchrotron emission with an optical turnover. Higher resolution spectroscopic follow-up would allow for more detailed fitting of J0045+41 to determine if a red quasar template yields a better fit.

The confusion of stars and quasars represents a unique problem for purely photometric surveys, such as the upcoming LSST project (Ivezic et al., 2008). Stars and higher redshift ($z>2.2$) quasars are well separated in color space. However, at lower redshifts, the two color loci appear closer and closer. The difference between the two populations is most apparent in $u$-band flux and $u-g$ color; indeed the $u$ filter was designed in part to leverage the difference between power-law spectra and spectra with strong Balmer decrements (Fukugita et al., 1996; Stoughton et al., 2002; Richards et al., 2002). Thus, in any single-visit catalog, the colors of the lower-redshift, low-luminosity, and intrinsically red AGN are the hardest population to distinguish from stars. LSST will visit most of its survey area $\sim 50$-180 times in each filter over 10 years. Peters et al. (2015) demonstrated that it is possible to use variability in addition to colors to distinguish stars from AGN with high ($>90\%$) completeness. However, the accuracy of classifications in the lowest redshift bins studied drops to $\sim 80\%$. While the number of quasars at low redshift is small, this highlights the importance of developing accurate classification algorithms for objects similar to J0045+41. Forthcoming work will focus on distinguishing between stars and quasars in the low-redshift, low-luminosity, red regime.

J0045+41 is an exciting and unique object. It represents an extreme end of color space in which photometric classification methods fail. Both the simple selection methods (described in §1) and more sophisticated machine learning algorithms are unable to correctly classify objects in this regime. Finding these intrinsically red AGN is important, as they are still poorly understood. The evidence of multiple periodic signals in the photometric lightcurve of J0045+41 is compelling, and warrants more dedicated spectroscopic observations at higher spectral resolution and deeper photometric observations sampled at a higher rate. Such observations would be crucial to confirm the presence of a SMBHB in J0045+41. They would also allow for the confirmation of the periods that we detected. The photometric data will soon be attainable in the form of the Zwicky Transient Facility (ZTF, Bellm 2014), a next-generation transient survey that will see first light this year.

The authors thank Jessica Werk, Julianne Dalcanton, Ben Williams, Jake VanderPlas and Scott Anderson for their valuable advice and feedback on this work. We wish to thank the anonymous referee for their extremely helpful comments. Based on observations (Program ID GN-2016A-FT-30) obtained at the Gemini Observatory (processed using the Gemini IRAF package), which is operated by the Association of Universities for Research in Astronomy, Inc., under a cooperative agreement with the NSF on behalf of the Gemini partnership: the National Science Foundation (United States), the National Research Council (Canada), CONICYT (Chile), Ministerio de Ciencia, Tecnología e Innovación Productiva (Argentina), and Ministério da Ciência, Tecnologia e Inovação (Brazil). The authors thank the Gemini-North support staff. The scientific results reported in this article are based in part on data obtained from the Chandra Data Archive and the SDSS archive. Funding for the Sloan Digital Sky Survey IV has been provided by the Alfred P. Sloan Foundation, the U.S. Department of Energy Office of Science, and the Participating Institutions. SDSS-IV acknowledges support and resources from the Center for High-Performance Computing at the University of Utah. The SDSS web site is www.sdss.org. SDSS-IV is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS Collaboration including the Brazilian Participation Group, the Carnegie Institution for Science, Carnegie Mellon University, the Chilean Participation Group, the French Participation Group, Harvard-Smithsonian Center for Astrophysics, Instituto de Astrofísica de Canarias, The Johns Hopkins University, Kavli Institute for the Physics and Mathematics of the Universe (IPMU) / University of Tokyo, Lawrence Berkeley National Laboratory, Leibniz Institut für Astrophysik Potsdam (AIP), Max-Planck-Institut für Astronomie (MPIA Heidelberg), Max-Planck-Institut für Astrophysik (MPA Garching), Max-Planck-Institut für Extraterrestrische Physik (MPE), National Astronomical Observatories of China, New Mexico State University, New York University, University of Notre Dame, Observatário Nacional / MCTI, The Ohio State University, Pennsylvania State University, Shanghai Astronomical Observatory, United Kingdom Participation Group, Universidad Nacional Autónoma de México, University of Arizona, University of Colorado Boulder, University of Oxford, University of Portsmouth, University of Utah, University of Virginia, University of Washington, University of Wisconsin, Vanderbilt University, and Yale University. This work was facilitated though the use of advanced computational, storage, and networking infrastructure provided by the Hyak supercomputer system at the University of Washington. This work made use of the following facilities and software:

Gemini-North (GMOS)