A Bayesian Analysis of SDSS J0914+0853, a Low-Mass Dual AGN Candidate

TL;DR

This paper introduces BAYMAX, a Bayesian tool for analyzing Chandra X-ray data to determine if an AGN is single or dual, and applies it to a low-mass dual AGN candidate, finding it most likely a single AGN.

Contribution

The paper presents BAYMAX, a novel Bayesian framework for classifying AGNs as single or dual sources using X-ray data, especially effective at small separations.

Findings

BAYMAX analysis indicates SDSS J0914+0853 is most likely a single AGN.

The probability of it being a dual AGN with flux ratio >0.3 and separation >0.3'' is very low.

BAYMAX will aid in systematic classification of dual AGN candidates.

Abstract

We present the first results from BAYMAX (Bayesian AnalYsis of Multiple AGN in X-rays), a tool that uses a Bayesian framework to quantitatively evaluate whether a given Chandra observation is more likely a single or dual point source. Although the most robust method of determining the presence of dual AGNs is to use X-ray observations, only sources that are widely separated relative to the instrument PSF are easy to identify. It becomes increasingly difficult to distinguish dual AGNs from single AGNs when the separation is on the order of Chandra's angular resolution (<1''). Using likelihood models for single and dual point sources, BAYMAX quantitatively evaluates the likelihood of an AGN for a given source. Specifically, we present results from BAYMAX analyzing the lowest-mass dual AGN candidate to date, SDSS J0914+0853, where archival Chandra data shows a possible secondary AGN ~0.3'' from the primary. Analyzing a new 50 ks Chandra observation, results from BAYMAX shows that SDSS J0914+0853 is most likely a single AGN with a Bayes factor of 13.5 in favor of a single point source model. Further, posterior distributions from the dual point source model are consistent with emission from a single AGN. We find the probability of SDSS J0914+0853 being a dual AGN system with a flux ratio f>0.3 and separation r>0.3'' to be very low. Overall, BAYMAX will be an important tool for correctly classifying candidate dual AGNs in the literature, and studying the dual AGN population where past spatial resolution limits have prevented systematic analyses.

Figures

Tables

Table 1. Table 1: Symbols

Symbol	Definition
(1)	(2)
($x_i,y_i$)	Sky coordinate of photon $i$
$E_i$	Energy of photon $i$, in keV
$n$	Total flux (counts) of given source
$\mu$	Central position of given source in sky coordinates (2D; $\mu =[{\mu }_x,{\mu }_y]$)
$k$	Number of Chandra observations being modeled
$\Delta x_K$	Translational astrometric shift in $x$ ($K=[1,\dots ,k-1]$)
$\Delta y_K$	Translational astrometric shift in $y$
$\Delta r_K$	Radial astrometric shift ($\Delta r_K=\sqrt{{\left(\Delta x_K\right)}^2+{\left(\Delta y_K\right)}^2}$)
$\Delta {\varphi }_K$	Rotational astrometric shift
$f$	Flux ratio between secondary and primary source (/)
$M_j$	Given model being analyzed by BAYMAX
${\theta }_j$	Parameter vector for $M_j$, i.e. [$\mu$, $f$, $\Delta x_K$, $\Delta y_K$, $\Delta {\varphi }_K$] .

Table 2. Table 2: Posterior Results for θ 1 subscript 𝜃 1 \theta_{1} and θ 2 subscript 𝜃 2 \theta_{2}

Parameter	Best-fit Value
(1)	(2)
Single Point Source Model
${\mu }_x$	4074.6 $\pm$ 0.1
${\mu }_y$	4063.6 $\pm$ 0.1
$\alpha$	138.70
$\delta$	+8.89
$\Delta x_1$	$-$14.8″$\pm$ 0$\stackrel{\prime \prime }{.}$1
$\Delta y_1$	$-$25.8″$\pm$ 0$\stackrel{\prime \prime }{.}$1
$\Delta r_1$	29.8″$\pm$ 0$\stackrel{\prime \prime }{.}$1
Dual Point Source Model
${\mu }_{x,P}$	4074.6 $\pm$ 0.1
${\mu }_{y,P}$	4063.6 $\pm$ 0.1
${\mu }_{x,S}$	4074.6 $\pm$ 1.5
${\mu }_{y,S}$	4063.4 $\pm$ 1.5
${\alpha }_P$	138.70
${\delta }_P$	+8.89
${\alpha }_S$	138.70
${\delta }_S$	+8.89
$r$	0$\stackrel{\prime \prime }{.}$15 $\pm$ 0$\stackrel{\prime \prime }{.}$15
$\log f$	$-$1.6 $\pm$ 0.4
$\Delta x_1$	$-$14.8″$\pm$ 0$\stackrel{\prime \prime }{.}$1
$\Delta y_1$	$-$25.8″$\pm$ 0$\stackrel{\prime \prime }{.}$1
$\Delta r_1$	29.8″$\pm$ 0$\stackrel{\prime \prime }{.}$1

Equations

$$\underset{\text{Posterior odds}}{\underbrace{\frac{P(M_2\mid D)}{P(M_1\mid D)}}}=\underset{\text{Bayes factor}}{\underbrace{\frac{P(D\mid M_2)}{P(D\mid M_1)}}}\times \underset{\text{Prior odds}}{\underbrace{\frac{P(M_2)}{P(M_1)}}},$$

$$BF=\frac{\int P(D\mid {\theta }_2,M_2)P({\theta }_2\mid M_2)d{\theta }_2}{\int P(D\mid {\theta }_1,M_1)P({\theta }_1\mid M_1)d{\theta }_1},$$

$$\begin{array}{rl}\mathcal{L}& =P(x,y\mid \mu ,E)=\prod _{i=1}^{n}P(x_i,y_i\mid \mu ,E_i)\\ & =\prod _{i=1}^n\frac{M_{1,i}({\theta }_1{)}^{D_i}}{D_i!}\exp (-M_{1,i}({\theta }_1)),\end{array}$$

$$\text{mbox}ln\mathcal{L}=\sum _{i=1}^n{D}_i\text{mbox}ln{M}_{j,i}({\theta }_j)+\text{mbox}constant,$$

$$\mu =\mathcal{U}(a,b),$$

$$\mu =\mathcal{N}({\mu }_m,{\sigma }^2),$$

Full text

galaxies: active — galaxies – X-rays — galaxies: interactions

††software: CIAO (v4.8; Fruscione et al. 2006),

XSPEC (v12.9.0; Arnaud 1996),

nestle (https://github.com/kbarbary/nestle),

PyMC3 (Salvatier et al., 2016),

SAOTrace (http://cxc.harvard.edu/cal/Hrma/Raytrace/SAOTrace.html),

MARX (v5.3.3; Davis et al. 2012)

1 Introduction

Given that almost all massive galaxies are thought to harbor nuclear supermassive black holes (SMBH; Kormendy & Richstone 1995) and that classical heirarchical galaxy evolution predicts that later stages of galaxy evolution are governed by mergers (e.g., White & Rees 1978), galaxy mergers provide a favorable environment for the assembly of active galactic nuclei (AGNs) pairs (Volonteri et al., 2003). The role galaxy mergers play in triggering AGN and/or AGN pairs remains unclear (e.g. Hopkins & Quataert 2010; Kocevski et al. 2012; Schawinski et al. 2012; Hayward et al. 2014; Villforth et al. 2014, 2017; Capelo et al. 2015), however both observations and simulations agree that AGN activity should increase with decreasing galaxy separation (e.g. Koss et al. 2012; Blecha et al. 2013; Ellison et al. 2013; Goulding et al. 2018; Capelo et al. 2017; Barrows et al. 2017).

“Dual AGNs” are usually defined as a pair of AGN in a single galaxy or merging system (with typical separations of $\sim$1 kpc), while a “binary AGN” is a pair of AGNs that are gravitationally bound with typical separations $\lessapprox$100 pc (see Begelman et al. 1980 for a summary of the main merging phases of SMBHs). Understanding the specific environments where dual AGNs occur provides important clues about black hole growth during the merging process. Additionally, as progenitors to SMBH-binary mergers, the rate of dual AGNs is intimately tied to gravitational wave events detectable by pulsar timing arrays and space-based interferometry (see Mingarelli 2019 and references within). Thus, dual AGNs offer a critical way to observe the link between galaxy mergers, SMBH accretion, and SMBH mergers.

The frequency of galaxy mergers in our observable universe implies that dual AGN should be relatively common. In particular, assuming a dynamical friction timescale of $\sim$1 Gyr we expect the galaxy merger fraction with separations $\leq$ 1 kpc by $z=0.1$ to be between $\sim$6–10% (Hopkins et al., 2010). However, these estimated merger fractions don’t take into account the AGN duty cycle, and observations of nearby AGN have shown that at separations $\leq$ 1 kpc the fraction of dual AGN may be much higher (Barrows et al., 2017). Yet, very few dual AGN with separations $<$ 1 kpc have been confirmed. Such systems become difficult to resolve with Chandra beyond $z\geq 0.05$, where separations on the order of 1 kpc approach Chandra’s angular resolution (where the half-power diameter is $\sim$0$\farcs$8 at $\sim$1 keV). For example, the closest dual AGN candidate identified using two resolved point sources with Chandra is NGC 3393 (Fabbiano et al., 2011) with a projected separation of $\sim$150 pc ($\sim$0$\farcs$6; however see Koss et al. 2015 for a critical analysis of the X-ray emission). Thus, many indirect detection techniques have been developed to search for evidence of dual and binary AGNs, primarily relying on optical spectroscopy and photometry.

Perhaps the most popular method of finding dual AGN candidates is via double-peaked narrow line emission regions (e.g., Zhou et al. 2004; Gerke et al. 2007; Comerford et al. 2009; Liu et al. 2010; Fu et al. 2012; Comerford et al. 2012, 2013). Double-peaked narrow lines can be a result of a dual AGN system during the period of the merger when their narrow line regions (NLR) are well separated in velocity. The system can display two sets of narrow line emission regions, such as [O III], where the separation and width of each peak will depend on parameters such as the distance between the two AGN. However the optical regime alone is insufficient in confirming a dual AGN candidate because of ambiguity in interpretation of the observed double-peaked narrow line regions. For example, bipolar outflows and rotating disks can also can produce the double-peaked emission feature (see, e.g., Greene & Ho 2005; Rosario et al. 2010; Müller-Sánchez et al. 2011; Smith et al. 2012; Nevin et al. 2016). Indeed, follow-up observations using high-resolution imaging and spatially resolved spectroscopy have found that many double-peaked dual AGN candidates are most likely single AGN (Fu et al., 2012; Shen et al., 2011; Comerford et al., 2015). Dual or binary AGN candidates can be confirmed using high resolution radio imaging (see Rodriguez et al. 2006; Fu et al. 2015; Müller-Sánchez et al. 2015; Kharb et al. 2017); however an absence of radio emission does not necessarily mean an absence of AGNs (as only $\sim 10\%$ of AGN are radio-loud), while a detection of two radio nuclei can have multiple physical explanations (such as star-forming nuclei). Nuclei can only be classified as AGN at radio frequencies if they are compact and have flat or inverted spectral indices (see, e.g., Burke-Spolaor 2011; Hovatta et al. 2014).

1.1 X-ray observations of dual AGN candidates

The most robust method of confirming the presence of dual AGNs is to use X-ray observations. Due to the relatively few possible origins of emission above $10^{40}$ erg s*-1* besides accretion onto a SMBH (Lehmer et al., 2010), X-rays are one of the most direct methods of finding black holes, especially with Chandra’s superb angular resolution. Unlike the optical regime, X-rays are less sensitive to absorption from the dusty environments of merger remnants. Currently, many analyses searching for dual AGN candidates using Chandra observations implement the Energy-Dependent Subpixel Event Respositioning (EDSER) algorithm (Li et al., 2004). EDSER improves the angular resolution of Chandra’s Advanced CCD Imaging Spectrometer (ACIS) by reducing photon impact position uncertainties to subpixel accuracy, and in combination with Chandra’s dithering can resolve sub-pixel structure down to the limit of the Chandra High Resolution Mirror Assembly. However, thus far it has only been used to make images and qualitatively analyze them for dual point sources. In the absence of corroborating evidence from other data, the reliance on visual interpretation of dual AGNs with separations comparable to Chandra’s resolution leads to both false negatives and false positives. This issue is worse in the low-count regime ($<$ 200 counts), where even dual AGN with larger separations ($>0\farcs{5}$) but low flux ratios are not clearly distinct.

We have developed a PYTHON tool BAYMAX (Bayesian AnalYsis of AGNs in X-rays) that allows for a quantitative and rigorous analysis of whether a source in a given Chandra observation is more likely composed of one or two point sources. This is done by taking calibrated events from a Chandra observation and comparing them to the expected distribution of counts for single or double source models. The main component of BAYMAX is the calculation of Bayes factor, which represents the ratio of the plausibility of observed data, given two different models. Values $>1$ or $<1$ signify which model is more likely (see Section 2 for explicit details). Further, BAYMAX returns the maximum likelihood values for the parameters of each model. In this paper we introduce our tool BAYMAX and present its analysis on the Chandra observations of the lowest-mass dual AGN candidate SDSS J0914+0853. Here we specifically highlight BAYMAX’s capabilities with respect to the Chandra observations of SDSS J0914+0853. We are using a subset of BAYMAX’s full capabilities, i.e., analyzing an on-axis source, assuming identical spectra for both the primary and secondary AGN, and the background contribution is deemed negligible. As well, false positives are only analyzed for regions in parameter space (such as count number, separation, and flux ratio) that are specific to SDSS J0914+0853. Our following paper (Foord et al. 2019 (in prep)) will expands upon the explicit details of BAYMAX, including its capabilities of correctly identifying dual AGN as a function of observed flux, angular separations, off-axis angle, and flux ratios. In this paper, we restrict our discussion to BAYMAX’s abilities on our observations of SDSS J0914+0853.

1.2 SDSS J0914+0853

SDSS J0914+0853 was originally identified by Greene & Ho (2007) as one of $\sim$200 low-mass SMBH based on “virial” black hole mass estimates, where the velocity dispersion and radius of the broad line region (BLR) were estimated from H$\alpha$ emission line characteristics. The system is at $z=0.14$ ($D_{L}=661$ Mpc and $D_{A}=509$ Mpc for a $\Lambda$CDM universe, where $H_{0}=69.6$, $\Omega_{M}=0.286$, and $\Omega_{\Lambda}=0.714$) and is a low-mass ($M_{BH}=10^{6.3}M_{\odot}$), low-luminosity AGN. SDSS J0914+0853 was observed by Chandra as part of a Cycle 13 program targeting low-mass AGNs (Proposal ID:13858, PI:Gültekin). These data were taken to investigate the fundamental plane in the low-mass regime and thus are on-axis (Gültekin et al., 2014). Analyzing the 15 ks Chandra data with EDSER, the archival Chandra exposure shows a possible secondary source 0$\farcs$3 away from the primary. Possible contamination from an ultraluminous X-ray source (ULX) is very low; following the methodology in Foord et al. (2017a) we calculate the number of expected ULXs with $L_{\geq}10^{41}$erg s*-1* to be $<10^{-3}$ within a radius of 0$\farcs$3 from the center of the galaxy. If the emission is found to most likely originate from two point sources, it will be the lowest-mass dual AGN discovered, and analysis of this system paves the way for a better understanding of the role of mergers and AGN activity in low-mass systems. In particular, dual AGNs in low-mass galaxies with low luminosities are the perfect testbed for discerning between competing models for the connection between galaxy mergers and AGN activity. It has been argued that mergers can trigger high-luminosity AGN but not low-luminosity AGN, which are triggered by stochastic processes (Hopkins & Hernquist, 2009; Treister et al., 2012). A competing hypothesis is that there is no correlation between AGN luminosity and mergers (e.g.,Villforth et al. 2014). Since dual AGNs most likely arise from mergers, the presence of a low-luminosity dual AGN in SDSS J0914+0853 would show that low luminosity AGNs can arise from mergers. However, effects due to pileup and artifacts from the Point Spread Function (PSF) cannot be ruled out at a high statistical confidence for the low-count ($\sim$250 counts between $2$–$7$ keV) image. At 10%, the pile-up fraction is relatively small, but combined with asymmetries in the Chandra PSF (Juda & Karovska 2010), it could produce a spurious dual AGN signature. Thus, a statistical analysis is necessary before a discovery can be confirmed.

We aim to unambigiously determine the true nature of SDSS J0914+0853. As stated above, the existing Chandra data cannot do this because (i) the pile-up introduces systematic uncertainties in the EDSER processing, (ii) the existing exposure is relatively shallow (15 ks), and (iii) potential PSF artifacts can produce spurious dual AGN signatures. To help determine the true nature of SDSS J0914+0853, we received a new observation (Proposal ID:19464, PI:Gültekin) that addresses all three of the above points. In particular, the observation (i) uses the shortest possible frame time with a subarray, thereby eliminating pileup, (ii) goes 3$\times$ deeper with a 50 ks exposure, and (iii) uses a substantially different roll angle so that any PSF artifacts will not appear in the same location on the sky. With a total of $\sim$723 counts between $2$–$7$ keV (combining both datasets), BAYMAX is able to statistically analyze the likelihood that SDSS J0914+0853 is a dual AGN for separations $>0\farcs 3$ and flux ratios $>0.1$.

The remainder of the paper is organized into 5 sections. Section 2 introduces Bayesian inference, focusing on the specific components of Bayes factor and how BAYMAX calculates the likelihood and prior densities. In Section 3 we analyze the Chandra observations of SDSS J0914+0853, including both a photometric and spectral analysis. In Section 4 we present our results when running BAYMAX on the Chandra observations of SDSS J0914+0853, and in Section 5 we discuss the sensitivity and limitations of BAYMAX across parameter space and how they affect our results. Lastly, we summarize our findings in Section 6. Please see Table 1 for a list of symbols used throughout this paper.

2 Methods

2.1 Bayesian Inference

BAYMAX is capable of statistically and quantitatively determining whether a given observation is better described by a model composed of one or two point sources based on a Bayesian framework. A Bayesian approach combines all available information (using prior distributions and likelihood models) to infer the unknown model parameters (posterior distributions). Bayes Theorem implies:

[TABLE]

where the posterior odds represents the ratio of the dual point source model ($M_{2}$) vs. the single point source model ($M_{1}$) given the data $D$; the Bayes factor ($BF$) quantifies the evidence of the data for $M_{2}$ vs. $M_{1}$, and the prior odds represents the prior probability ratio of $M_{2}$ vs. $M_{1}$. Specifically, the Bayes factor is the ratio of the marginal likelihoods:

[TABLE]

representing the ratio of the plausibility of observed data $D$, given two different models, and parameterized by the parameter vectors $\theta_{2}$ and $\theta_{1}$. Values $>$1 or $<$1 signify whether $M_{2}$ or $M_{1}$ is more likely (see Jeffreys 1935 and Kass & Raftery 1995 for the historic interpretations of the strength of a $BF$ value; we analyze our data to define a “strong” $BF$ value in Section 5.) In this paper, we assume that $M_{2}$ and $M_{1}$ are equally probable, so that $P(M_{2})=P(M_{1})=0.5$ and the Bayes factor directly represents the posterior odds. Thus calculating the Bayes factor can be broken into two components — the likelihood density, $P(D\mid\theta_{j},M_{j})$, and the prior density, $P(\theta_{j}\mid M_{j})$.

2.2 Data Structure and Modeling the PSF

In this section we will focus on the likelihood density implemented in BAYMAX. Each reprocessed Chandra level-2 event file tabulates the directional coordinates ($x_{i}$, $y_{i}$) and energy $E_{i}$ for each detected photon, where $i$ indexes each detected photon (See Table 1 for a summary of notation). The detector itself records the pulse height amplitude (PHA) of each event, which is roughly proportional to the energy of the incoming photon. In the reprocessed files, the energy $E_{i}$ is calculated from the event’s PHA value, using the appropriate gain table. Thus, BAYMAX takes calibrated events ($x_{i}$, $y_{i}$, $E_{i}$) from reprocessed Chandra observations and compares them to new simulations based on single and dual point source models.

We characterize the properties of the Chandra PSF by simulating the PSF of the optics from the High Resolution Mirror Assembly (HRMA) via ray tracing simulations. The two primary methods to simulate the HRMA PSF are SAOTrace111http://cxc.harvard.edu/cal/Hrma/Raytrace/SAOTrace.html and the Model of AXAF Response to X-rays (MARX, Davis et al. 2012). While the MARX model uses a slightly simplified (and faster) description of the HRMA, differences between SAOTrace and MARX simulations are minimal for on-axis simulations. For our PSF analysis below, we find consistent parametric fits between an SaoTrace generated PSF and one generated by MARX – in particular the root-mean-square error between the two fits is on the order of $\sim 0.1\%$

Thus, our Likelihood models for single and dual point sources are created by parametrically modeling the Chandra PSF using high count simulations created by MARX-5.3.3. To translate the PSF model to an event file, the HRMA ray tracing simulations are projected on to the detector-plane via MARX. Ray tracing simulations generated by both MARX and SAOTrace will have roughly the correct total intensity, but small deviations in the overall shape. Specifically regarding MARX – the PSF wings are broader than observations while the PSF core is narrower than observed (Primini et al., 2011). These discrepancies can be reduced by blurring the PSF when projecting it to the detector-plane via the AspectBlur parameter. This parameter is used to account for the uncertainty in the determination of the aspect solution (such as effects from pixel quantization and pixel randomization), as well as the uncertainty in the instrument and dither models within MARX. The best value should be considered carefully for each unique observation222http://cxc.harvard.edu/ciao/why/aspectblur.html. For MARX generated simulations on ACIS-S, we expect the AspectBlur parameter to have values between $0\farcs 25-0\farcs 28$. For our PSF analysis we set AspectBlur to $0\farcs 28$. We note that value used for AspectBlur does not represent the accuracy at which we can centroid.

For a given observation, a user-defined source model is input to MARX to generate X-ray photons incident from a single point source centered on the observed central position of the AGN ($\mu_{obs}$, defined as the coordinates where the hard X-ray emission from the AGN is estimated to peak). Because we do not model the spectral parameters of the system (see Section 1), we are only interested in modeling the spatial distribution of a photon due to its energy $E_{i}$ and our PSF does not depend on the spectral shape of our model. Each simulation uses the observation-specific detector position (RA_Nom, Dec_Nom, Roll_Nom) and start time (TSTART). We set the number of generated rays (NumRays) to $1\times 10^{6}$ and the read-out strip is excluded by setting the parameter ACIS_Frame_Transfer_Time to 0.

We model the PSF as a summation of 2D Gaussians, where the amplitude and standard deviation of each Gaussian is energy-dependent. In general, the PSF may be any function which is unique to a given observation and can be quickly evaluated. For both the 15 ks and 50 ks observation, we fit a variety of possible functions to the PSF. Using the Bayesian Information Criterion as a diagnostic for model comparison, we find that a summation of three circular concentric 2D Gaussians yields the best-fit for both Chandra observations (specifically, we find the PSF wings are best-modeled by the broadness of a Gaussian component versus the addition of a Lorentzian component). We model the PSF for each observation individually, however we find that the best-fit parameters for each PSF model are consistent with one-another within the $1\sigma$ error bars (which is not surprising, given that both sources were observed on-axis for ACIS-S albeit in different Chandra cycles).

Each photon is assumed to originate from a single or dual point source system. For example, for a single point source, the probability that a photon observed at location $x_{i}$,$y_{i}$ on the sky with energy $E_{i}$ is described by the PSF centered at $\mu$ is $P(x_{i},y_{i}\mid\mu,E_{i})$, i.e., the energy dependent PSF. For $n$ total events, the total probability is the product of the probability for each detected photon, i.e., the likelihood density is:

[TABLE]

where we use the Poisson likelihood, appropriate given that Chandra registers each event individually. Here, $M_{1,i}(\theta_{1})$ is the probability for event $i$ given our PSF model, and $D_{i}$ is the data value for event $i$. For a dual point source the total probability is $P(x,y\mid\mu_{P},\mu_{S},E,n_{S}/n_{P})$, where $\mu_{P}$ and $\mu_{S}$ represent the location of the primary and secondary AGN. The ratio of the fluxes (or, total counts) between the secondary and primary is represented by $n_{S}/n_{P}=f$, where $0\leq f\leq 1$. We note that our analysis on SDSS J0914+0853 does not include fitting for the spectral models. Using the archival data we find consistent hardness ratios between the candidate primary and secondary AGN, where we use circular and non-overlapping apertures centered on their apparent locations. Thus, we assume that the spectra are the same spectral shape as that for the entire system, but with different normalizations. Future analyses with BAYMAX will include fitting for different spectral shapes.

Because Bayes factor represents the ratio of likelihood densities, and we use the same data across both models, our calculations become simplified. We are left with:

[TABLE]

where $M_{j,i}(\theta_{j})$ is calculated for either a single ($j=1$) or dual ($j=2$) point source model via our parametrically-fit PSF, for each detected event $i$.

2.3 Prior Distributions

BAYMAX requires user input regarding (i) the number of datasets and (ii) the prior distributions for each parameter. Regarding point (i), SDSS J0914+0853 has $k=2$ observations and thus the parameter vector $\theta_{1}=[\mu,\Delta x_{1},\Delta y_{1}]$ while $\theta_{2}=[\mu_{{P}}$, $\mu_{{S}},{\log{f}},\Delta x_{1},\Delta y_{1}]$. Here, $\mu=(\mu_{x},\mu_{y})$ is the central sky $x$,$y$ positions of the AGN; $\Delta x_{1}$ and $\Delta y_{1}$ account for the translational components of the relative astrometric registration for the $k-1$ observation; and $\log{f}$ is the log of the flux ratio where $f=n_{S}/n_{P}$. The relative astrometric registration adds an uncertainty that must be taken into account in order to avoid spurious dual AGN signals that can be generated from slight mismatches between two or more observations. We take this into account by including the astrometric registration of multiple observations as a set of parameters to be marginalized over. For SDSS J0914+0853, we find that the rotational component of the relative astrometric registration is expected to be very small ($\Delta\phi_{1}<1^{\circ}$) and including the parameter does not affect our results. Thus, we only include $\Delta x_{1}$ and $\Delta y_{1}$, which are analyzed for the shallower observation (i.e., relative to the 50 ks exposure).

Regarding point (ii), BAYMAX can incorporate any user-defined function to describe the prior distributions for each parameter. For SDSS J0914+0853, the prior distributions of $\mu$ for both $M_{1}$ and $M_{2}$ are described by a continuous uniform distribution333The continuous uniform distribution $\mathcal{U}(a,b)$ is a probability distribution where all values between the minimum $a$ and maximum $b$ are equally probable.:

[TABLE]

where we constrain all $\mu$ values to be between $a=\mu_{obs}-2$ and $b=\mu_{obs}+2$. Thus, the 2D parameter space for possible $\mu_{x}$ and $\mu_{y}$ is a 4$\times$4 sky-pixel box ($\approx$ 1.98″$\times$1.98″) centered on the observed central X-ray coordinates of SDSS J0914+0853. Further, the prior distributions of $\Delta x_{1}$ and $\Delta y_{1}$ are also described by a uniform distribution with $a=\delta\mu_{obs}-3$ and $b=\delta\mu_{obs}+3$, where $\delta\mu$ represents the difference between the observed central X-ray coordinates of the two observations (in practice, $\delta\mu$ is expected to be small; however because the most recent observation of SDSS J0914+0853 was taken in a subarray mode, the difference between the aimpoints of the two observations is $\approx$15 sky pixels). For $M_{2}$, the prior distribution for $\log{f}$ is also described by a uniform distribution (and thus $f$ is described by a log uniform distribution), where $a=-2$ and $b=0$. The range for the prior distribution of $\log{f}$ covers possible values expected for “major mergers” (with mass ratios $>1/3$), while accounting for a large range of possible Eddington fractions between the two black holes. In general, informative priors can be incorporated if prior information is available. For example, we might set the prior distributions of $\mu_{P}$ and $\mu_{S}$ to Gaussian distributions centered on coordinates that are better constrained by other observations (such as spectroastrometric [O III] observations, or complementary IR photometry).

2.4 Calculation of Bayes Factor

Bayesian inference can be divided into two categories: model selection and parameter estimation. In this section we review how we address each component in BAYMAX.

Computing the marginal likelihood is challenging, as it involves a multi-dimensional integration over all of parameter space. Only over the last $\sim$20 years have the advances in computational power allowed Bayesian inference to become a more common technique for model selection. In addition, general numerical methods based on Markov chain Monte Carlo (MCMC; e.g., see Metropolis et al. 1953) have been developed, allowing one to conduct Bayesian inferences in an efficient manner, with few constraints on dimensionality or analytical integrability.

To calculate the marginal likelihood, BAYMAX implements a sampling technique called nested sampling (Skilling, 2004). In nested sampling, the marginal likelihood is rebranded as the “Bayesian evidence”, denoted by $Z$. Here, $Z=\int P(D\mid\theta_{j},M_{j})P(\theta_{j}\mid M_{j})d\theta_{j}$. Nested sampling transforms the multi-dimensional integral to a one dimensional integral by introducing the prior mass $X$, defined as $X=\int_{\mathcal{L}(\theta)>\lambda}P(\theta_{j}\mid M_{j})d\theta_{j}$; here the integral extends over the regions of parameter space contained within the iso-likelihood contour $\mathcal{L}(\theta)=\lambda$ and at any given time has a value $0<X<1.$ For a step-by-step explanation of nested sampling we refer the reader to Skilling (2004), Shaw et al. (2007), Feroz & Hobson (2008), and Feroz et al. (2009). A direct effect of the nested sampling methodology is sparsely sampling in low likelihood regions and densely sampling where the likelihood is high. To calculate the evidence, BAYMAX uses the PYTHON package nestle444https://github.com/kbarbary/nestle. The package provides a pure-PYTHON implementation of nested sampling, where prior mass space can be sampled via different techniques. In particular, we use multi-ellipsoidal sampling (by setting method=‘multi’; see Mukherjee et al. 2006; Shaw et al. 2007; Feroz & Hobson 2008).

The two parameters that affect the accuracy of $Z$ are the number of “active points” and the stopping criterion dlog$Z$. The number of active points represent how many points in prior mass space one is sampling at a given time (roughly analogous to the number of walkers in an MCMC run). The stopping criterion determines when the nested sampling loop terminates — when the current largest sampled likelihood does not increase by more than the stopping criterion value, the sampling will end555Specifically, at a given iteration $i$ where the current evidence is $Z_{i}$ and the estimated remaining evidence in the likelihood landscape is $Z_{est}$, if $\mbox{log}(Z_{i}+Z_{est})-\mbox{log}(Z_{i})<\mbox{dlog}Z$, the sampling will terminate. For our analysis of SDSS J0914+0853, we use 500 active points (generally, a lower limit on the number of active points is 2$N_{dim}$, and due to using the multi-ellipsoidal method we add additional points to characterize each mode well) and set the stopping criterion dlog$Z=0.1$. By continually increasing the number of active points and decreasing the stopping criterion, the estimated evidence and its accuracy should converge. We find no significant difference in results when increasing our active points above 500 and using dlog$Z<0.1$, and conclude that these values have properly sampled the likelihood space.

Lastly, for parameter estimation BAYMAX uses PyMC3 (Salvatier et al., 2016), which uses gradient-based MCMC methods for sampling. Specifically, we use PyMC3’s built in Hamiltonian Monte Carlo (HMC) sampling method. HMC uses the gradient information from the likelihood to much more quickly converge than normal Metropolis-Hastings sampling. In general, HMC is more powerful for high dimensionality and complex posterior distributions (see, e.g., Betancourt et al. 2014).

In Figure 1 we show the results when we analyze two simulations with BAYMAX. The simulations have been reprocessed using EDSER, and binned by $2/3$ of the native pixel size. We simulate a single and a dual AGN system via MARX using the same telescope configuration as our new 50 ks observation. Both simulations have $n=700$ photons between $2$–$7$ keV, and the same $2$–$7$ keV spectrum as SDSS J0914+0853. We do not include a background contribution in our simulations (see Section 3). For the dual AGN simulation, each AGN has the same spectra but with normalizations such that the flux ratio $f=n_{S}/n_{P}=0.8$ while the separation between the two AGN is 0$\farcs$4. As evident, it is difficult to visually distinguish whether a given simulation is actually composed of one or two sources. Using the methodology presented above, BAYMAX favors the correct model for both simulations: for the single AGN simulation BAYMAX estimates a $BF$ of $17\pm 1.6$ in favor of the single point source model, while for dual AGN simulation BAYMAX estimates a $BF$ of $25\pm 1.5$ in favor of the dual point source model (error bars have been determined by running BAYMAX multiple times on each simulation, see Section 4). Although the hard X-ray emission appears quite similar between the two simulations, the joint posterior distributions are significantly different from one another. Specifically, for the dual AGN simulation the joint posterior distribution is more tightly concentrated around the true values. BAYMAX is able to recover the true separation and flux ratio within the 68% credible interval. However, for the single AGN, the separation and flux ratio are consistent with 0 at the 99.7% confidence level. We note that this particular joint-distribution shape (“L” shape) is consistent with a single AGN, where at very large flux ratios the dual AGN candidate is likely to have $r=0$, and at very large separations the dual AGN candidate is likely to have $\log{f}=-2$. More specifically, the dual point source model places each AGN at the same location and/or with arbitrarily low $f$, effectively consistent with one point source.

3 Data Analysis

3.1 X-ray Data

SDSS J0914+0853 was originally targeted to study low-mass AGNs and their relation to the plane of black hole accretion. The quasar was placed on the back illuminated S3 chip of the Advanced CCD Imaging Spectrometer (ACIS) detector, with an exposure time of 15 ks (Obs ID: 13858). We received a new 50 ks exposure at a roll angle significantly different from the previous observation, and using the smallest subarray (1/8) on a single chip to get the shortest standard frame time (Obs ID: 19464). This was done to (1) place the PSF artifact in a different location, (2) avoid pileup, and (3) receive $\sim$ 2–3 times more counts. We re-reduced and re-analyzed the archival data to ensure a uniform analysis between the two datasets.

We follow a similar data reduction as described in previous X-ray studies analyzing AGN (e.g., Foord et al. 2017a, b), using the Chandra Interactive Analysis of Observations (CIAO) v4.8 (Fruscione et al., 2006). Both datasets are analyzed with the energy-dependent sub-pixel event repositioning algorithm (EDSER; Li et al. 2004), which can be included in the standard CIAO reprocessing command chandra_repro with the parameter pix_adj=EDSER. For each observation, we first evaluate the aspect solutions of the reprocessed level-2 event files to ensure the Kalman lock was stable at all times. Further, we inspect the event detector coordinates as a function of time and find that they followed the instrument’s dither pattern, indicating no aspect-correction based degradation of the PSF.

We then correct for astrometry, cross-matching the Chandra detected point-like sources with the Sloan Digital Sky Survey Data Release 9 (SDSS DR9) catalog. The Chandra sources used for cross-matching are detected by running wavdetect on the reprocessed level-2 event file. We require each matched pair to be less than 2″ from one another and have a minimum of 3 matches. The 15 ks observation meets the criterion for an astrometrical correction; we find 8 matches between the Chandra observation and the SDSS DR9 catalog, resulting in a shift less than 0$\farcs$5. The 50 ks observation was taken in a subarray, and thus does not meet these criterion (however, because BAYMAX takes into account astrometric shifts between observations this step will not affect our final results, see Section 4). Background flaring is deemed negligible as neither dataset contain intervals where the background rate is 3$\sigma$ above the mean level.

We then rerun wavdetect on filtered 0.5 to 7 keV data to generate a list of X-ray point sources. We use wavelets of scales 1, 1.5, and 2.0 pixels using a 1.5 keV exposure map, and set the detection threshold significance to $10^{-6}$ (corresponding to one false detection over the entire S3 chip). We identify the quasar as an X-ray point source 0$\farcs$4 (Obs ID: 13858) and 0$\farcs$7 (Obs ID:19464) from the SDSS-listed optical center (2″corresponds to 95% of the encircled energy radius at 1.5 keV for ACIS). Counts are extracted from a 2″ radius circular region centered on the X-ray source center, where we use a source-free annulus with an inner radius of 20″ and outer radius of 30″ for the background extraction. We compare the estimated background contribution from the datasets to the Chandra blank-sky files. Here, the blank-sky files are properly scaled in exposure time and have matching WCS coordinates, dimensions, and energies. We find consistent results, where within $2$–$7$ keV we expect $\lesssim$ 1 and 1.5 background counts within a 4$\times$4 sky-pixel box ($\approx$ 1.98″$\times$1.98″) centered on the quasar for the archival and new dataset, respectively.

Our reduced data are shown in Figure 2. Here, both exposures have been reprocessed using the EDSER algorithm and are binned by a tenth of the native pixel size. The archival data appear to have sub-pixel structure, indicating a possible secondary AGN $\sim$ 0$\farcs$3 away from the primary, however our new observation is inconsistent with this picture. Although the X-ray emission may be slightly extended in the East-West direction (North is up, while East is left), we find minimal photometric evidence supporting the presence of a secondary AGN.

3.2 Spectral Fitting

The quasar’s net count rate and flux value are determined using XSPEC, version 12.9.0 (Arnaud, 1996). All errors evaluated in this section are done at the 95% confidence level, unless otherwise stated, and error bars quoted are calculated with Monte Carlo Markov chains via the XSPEC tool chain. We implement the Cash statistic (cstat; Cash 1979) in order to best assess the quality of our model fits.

Both spectra have an excess of flux at soft X-ray energies ($<1$ keV) with respect to the power-law continuum, while the 15 ks data appear to catch the source in a higher flux state in the soft X-ray band (see Fig. 3). Both of these behaviors are seen in AGN with a “soft excess” component (e.g., Lohfink et al. 2012, 2013; see Miniutti et al. 2009; Ludlam et al. 2015 for examples of soft X-ray excess in low-mass AGN candidates), an excess in emission above the extrapolated $2$–$10$ keV flux that is detected in over 50% of Seyfert 1s (Halpern, 1984; Turner & Pounds, 1989; Piconcelli et al., 2005; Bianchi et al., 2009; Scott et al., 2012). The physical origin of the soft excess remains uncertain; the shape is suggestive of a low-temperature high optical depth Comptonization of the inner accretion disk, however the temperature of this component appears to be constant over a wide range of black hole masses (and thus inferred accretion disk temperatures; see Gierliński & Done 2004; Crummy et al. 2006). The two most popular explanations for the soft excess are blurred ionized reflection from the inner parts of the accretion disk (e.g., Fabian et al. 2002, 2005; Gierliński & Done 2004; Crummy et al. 2006) and Comptonization components (such as partial covering of the source by cold absorbing material, see Boller et al. 2002; Tanaka et al. 2004).

Indeed, we find a statistically better fit when including an absorbed redshifted blackbody component to account for this soft excess (phabs$\times$zphabs$\times$(zpow + zbbody)). We fix the Galactic hydrogen column density (the photoelectric absorption component phabs) to $4.21\times 10^{20}$ cm*-2* (Kalberla et al., 2005), and the redshift to $z=0.14$. In Figure 3, we show the X-ray spectrum of both observations, along with the best-fit XSPEC models. For the archival dataset, we find best-fit values for intrinsic $N_{H}=3.38^{+0.10}_{-0.10}\times 10^{20}$ cm*-2*, power law component $\Gamma=2.01^{+0.11}_{-0.12}$; and blackbody component $kT=0.10^{+0.03}_{-0.05}$ keV. For our new dataset, we find best-fit values for intrinsic $N_{H}=4.07^{+0.10}_{-0.10}\times 10^{20}$ cm*-2*, power law component $\Gamma=2.51^{+0.10}_{-0.12}$; and blackbody component $kT=0.11^{+0.05}_{-0.04}$ keV.

However, because our analysis with BAYMAX is restricted to the $2$-$7$ keV photons from SDSS J0914+0853, our results are not affected by the soft emission component in the spectrum. In particular, although we detect variability between the two observations in the low-energy band, the $2$–$10$ kev fluxes are consistent with one another (at the 99.7% C.L.) when we fit each spectra independently between $2$–$7$ keV with an absorbed redshifted power law. For the 15 ks observations we calculate a total observed $2$–$10$ keV flux of $3.20^{+0.90}_{-0.80}\times 10^{-13}$ erg s*-1* cm*-2*, while for the 50 ks observation we calculate a total observed $2$–$10$ keV flux of $2.23^{+1.0}_{-0.49}\times 10^{-13}$ erg s*-1* cm*-2* s*-1*. This corresponds to a rest-frame $2$–$10$ keV luminosity of $1.83^{+0.31}_{-0.40}\times 10^{43}$ erg s*-1* and $1.25^{+0.35}_{-0.21}\times 10^{43}$ erg s*-1* at $z=0.14$ (assuming isotropic emission).

4 Results

Analyzing the 15 ks Chandra data with the EDSER option enabled, SDSS J0914+0853 appears to be an interesting dual AGN candidate. When binned, the data show a possible secondary source 0$\farcs$3 away from the primary (see Fig. 2). Although a possibly interesting result, classifying the source based on a qualitative analysis runs the risk of a false positive. A statistical analysis is necessary before a discovery can be confirmed. With an abundance of photons, and a robust model of the Chandra PSF, in the following section we aim to unambigiously determine the true nature of SDSS J0914+0853. We first analyze each observation individually using BAYMAX, and then combine the two (yielding a total of $n=723$ counts between $2$–$7$ keV).

We restrict our analysis to photons with (i) energies between $2$–$7$ keV and (ii) contained within a 4$\times$4 sky-pixel box (1.98″$\times$ 1.98″) centered on the nominal X-ray coordinates of the AGN. This corresponds to $\sim$95% of the encircled energy radius for the $2$–$7$ keV photons. Because we expect $\lesssim$ 1 and 1.5 background counts within this region for the archival and new dataset, each photon is assumed to originate from either one ($M_{1}$) or two ($M_{2}$) point sources, with no background contamination. The asymmetric PSF feature is within this extraction region, and sits approximately 0$\farcs$7 from the center of the AGN (see Fig. 2). Within $2$-$7$ keV, there are $6$ and $12$ photons that spatially coincide with the feature for the 15 and 50 ks observations, respectively. We mask the feature in both exposures before running BAYMAX.

We run BAYMAX with the initial conditions for the parameter vectors $\theta_{1}$ and $\theta_{2}$ as stated in Section 2. When running BAYMAX on our 15 ks and 50 ks observation individually, $k=1$ and thus we exclude the $\Delta x_{1}$ and $\Delta y_{1}$ from $\theta_{1}$ and $\theta_{2}$. Further, we run BAYMAX with the initializations for nestle as described in Section 2, with 500 active points and dlog$Z=0.1$

Our 15 ks observation has a total of $n=251$ counts between $2$–$7$ keV, while our 50 ks observation has a total of $n=472$ counts between $2$–$7$ keV. Using BAYMAX, we calculate a Bayes factors (defined as the ratio of the evidence for the dual point source model to the single point source model) of $\frac{Z_{2}}{Z_{1}}=0.154$ and $\frac{Z_{2}}{Z_{1}}=0.102$ for the 15 ks and 50 ks observations, respectively. This represents a Bayes factor of $\approx 6.5$ and $\approx 9.8$ in favor of a single point source. The relative magnitudes of the $BF$ values are not surprising – because the 15 ks observation has fewer counts than the 50 ks observation we expect there to be less evidence in favor of a given model. Indeed, using the definitions presented in Kass & Raftery (1995), both of these $BF$ values are considered “positive” against the dual point source model. Further, the posterior distributions for $\theta_{2}$ are consistent across both datasets: the best-fit locations for $\mu_{P}$ and $\mu_{S}$ are consistent with one-another at the 95% confidence interval, and the joint posterior distributions have shapes consistent with a single point source (consistent with the “L” shape seen in Fig. 1).

Given that the individual analyses on each dataset favor the same model, and that we can treat the two spectra as the same between $2$–$7$ keV, we increase our statistical power and combine both datasets. This yields a total of $n=723$ counts between $2$–$7$ keV. Although we analyze the two observations jointly, we emphasize that the likelihoods for each observation are calculated independently of one another, and are a function of their respective PSF models. Here, $k=2$ and $\Delta x_{1}$ and $\Delta y_{1}$ are included in parameter vectors for each model. We use BAYMAX to calculate a Bayes factor $\frac{Z_{2}}{Z_{1}}=7.40\times 10^{-2}$. This represents a Bayes factor of $\approx$13.5 in favor of a single point source.

To test the impact of the MCMC nature of nested sampling, we run BAYMAX multiple times on the combined datasets. We find consistent results, with a spread in ln$BF$-space well-described by a Gaussian distribution centered at ln$BF=2.6$ with standard deviation of 0.2. This Bayes factor strongly supports that the single point source model best describes the X-ray emission from SDSS J0914+0853. In Table 2 we list the best-fit values (defined as the median value of the posterior distributions) for parameter vector $\theta_{1}$.

We examine the posterior distributions for $\theta_{2}$ to better understand our results. In Figure 4 we show the combined $2$–$7$ keV dataset ($\approx$ 65 ks, where the photons associated with the 15 ks exposure have been spatially shifted by the most-likely $\Delta x_{1}$ and $\Delta y_{1}$) with the best-fit $x$ and $y$ sky-coordinates for the primary and secondary AGN ($\mu_{P}$ and $\mu_{S}$), as well as the joint posterior distribution for the separation, $r$, and log flux ratio, $\log{f}$, parameters. Here, $r=\sqrt{(\mu_{x,P}-\mu_{x,S})^{2}+(\mu_{y,P}-\mu_{y,S})^{2}}$. Spatially, the best-fit locations for $\mu_{P}$ and $\mu_{S}$ are consistent with one-another at the 95% confidence interval. Further, the joint posterior distribution has a shape consistent with a point source — the median values of the marginal posterior distributions are $r=0.15\pm 0.5$ and $\log{f}=-1.6\pm 0.4$, at very large flux ratios ($\log{f}\rightarrow{0}$) the separation is consistent with 0, and at very large separations ($r\rightarrow{2\arcsec}$) the flux ratio is consistent with 0. The best-fit values for all the parameters in parameter vector $\theta_{2}$ are listed in Table 2.

We investigate the influence of our prior distributions on our results. In particular, the Bayesian evidence automatically implements Occam’s razor — the simpler model will be more easily favored than the more complicated one, unless the latter is significantly better at explaining the data. For our analysis, this means that the dual point source model needs enough data to overcome the inherent bias that BAYMAX has towards favoring the single point source model. Whenever the prior distribution is relatively broad compared with the likelihood function, the prior has fairly little influence on the posterior. Thus, we re-run BAYMAX with Gaussian prior distributions for $\mu$, $\mu_{P}$, and $\mu_{S}$:

[TABLE]

where $\mu_{m}$ and $\sigma^{2}$ represent the mean and variance of the distribution. For $\mu_{P}$ and $\mu_{S}$, we set $\mu_{m}$ to the nominal X-ray positions of the potential primary and secondary AGN and set $\sigma$ to the observed separation between the two ($\sim 0\farcs 3$), given the 15 ks archival observation. For $\mu$ we set $\mu_{m}$ to the nominal X-ray position of the AGN, and similarly set $\sigma$ to 0$\farcs$3. BAYMAX calculates a Bayes factor of 10.8 $\pm$ 1.2, consistent within the errors of our previous measurement. Further, the posterior distributions returned by BAYMAX are consistent with those listed in Table 2. We conclude that using sharper priors (comparable to the sharpness of our likelihoods), has no effect on our results.

5 Discussion

Our results support the hypothesis that the low-mass dual AGN candidate SDSS J0914+0853 is instead a single AGN. Individually, we find $BF$ values of 6.5 and 9.8 in favor of a single point source for the 15 ks and 50 ks observations. When we combine the two datasets for a joint analysis, we find a $BF\sim 13.5$ in favor of a single point source, and the posterior distributions are consistent with this model. Further, the prior distributions do not appear to have a great influence on our posteriors, reflecting that the data should be sufficient to favor the correct model, even when accounting for the Bayesian bias. In the following section we discuss the significance of our results by analyzing BAYMAX’s capabilities across a range of parameter space for both the single and dual point source models. Assuming that SDSS J0914+0853 is indeed a dual AGN system, we investigate how the $BF$ determined by BAYMAX depends on parameters $r$ and $f$. In particular, we aim to understand where in parameter space BAYMAX loses sensitivity for simulations with a comparable number of counts as our observations.

5.1 BAYMAX’s Sensitivity Across Parameter Space

The first step is to investigate how well BAYMAX can classify a sample of simulated single AGN, i.e., our frequency of false-positives. This measurement will allow us to better define a range of Bayes factors that we can consider “strongly” support the dual point source model. We simulate 100 single AGN via MARX, assuming the same telescope configuration and spectrum as our new dataset. Further, each simulation has 700 photons between $2$–$7$ keV. We analyze each simulation with BAYMAX and find that only $2$ are misclassified as a dual AGN with $BF>3$ (with the largest $BF=3.5$). Thus, we define a $BF>3$ in favor of a dual AGN as “strong evidence”, while anything below this cut is classified as inconclusive.

We then run BAYMAX on a suite of simulated dual AGN systems, generated via MARX. The simulations were created with the same assumptions as listed above. Each simulation has 700 photons between $2$–$7$ keV, and each simulated AGN has the same $2$–$7$ keV spectrum as SDSS J0914+0853, but with normalizations proportional to their flux ratio. We simulated systems with separations that range between $0\farcs{3}$–$0\farcs{5}$ and flux ratios that range between $0.1$-$1.0$. For each $r$–$f$ point in parameter space, we evaluate 100 simulations with randomized position angles between the primary and secondary. Our results are shown in Figure 5, where we plot the logarithm of the mean $BF$ for each point in parameter space. Consistent with expectations, BAYMAX favors the dual point source model more strongly as the separation and flux ratio of a given dual AGN simulation increases, where we can expect $BF$ on the order of $\approx 10^{7}$ for systems with $r\geq 0\farcs{5}$. We enforce a cut of $BF>3$, where only $BF$ above this value are classified as strongly in favor of the dual point source model. We find that we are sensitive to most flux ratios where $r\geq 0\farcs{35}$, and for the smallest separations ($r\leq 0\farcs{35}$) BAYMAX is capable of identifying the correct model when $f\geq 0.8$

5.2 A Quasi-Frequentist Approach

Our analysis is intended to be a fully Bayesian inference, however some readers may find a frequenstist interpretation more intuitive. In the following section, we describe a potential interpretation of our results using a quasi-frequentist perspective.

On average, for separations below 0$\farcs$35, BAYMAX will not necessarily favor the correct model for a dual AGN system. For SDSS J0914+0853, we estimate a possible separation of 0$\farcs$3, given the shallower Chandra observation. However, the strength of the Bayes factor in favor of a single AGN has its own significance. From a frequentist perspective, we may ask what is the probability of measuring a $BF\geq 13.5$ in favor of a single AGN if the system is dual AGN. In this specific scenario, our “null hypothesis” is that SDSS J0914+0853 is a dual AGN system and our $p$-value represents the probability of measuring a $BF\geq 13.5$ in favor of a single AGN. Using our suite of dual AGN simulations, we analyze the probability of measuring a $BF\geq 13.5$, as a function of $r$ and $f$. Across all of parameter space, we find $p\leq 0.05$ and thus reject the null hypothesis at a 95% confidence level. If we set our $p$-value threshold to $p<0.03$, we find that only for the smallest flux ratios ($f<0.2$) can the null hypothesis not be rejected for $r<0.3$ (see Fig. 1). Thus, the probability of SDSS J0914+0853 being a dual AGN system with a (1) flux ratio $f>0.3$, (2) separation $r>0\farcs{3}$, and (3) measured $BF=13.5$ in favor of a single AGN, is very low.

We find that for observations with 700 counts BAYMAX is sensitive to a large region of $r$–$f$ parameter space, such that if SDSS J0914+0853 were a dual AGN, we expect different results. Our results and discussion highlight the importance of a robust, quantitative analysis of dual AGN candidates that are classified by their X-ray emission. Most candidate dual AGN are discovered via indirect detection methods, such as narrow-line optical spectroscopy or optical/IR photometry. However, directly detecting the X-ray emission unambigiously associated with a AGN is necessary for confirmation. For candidate AGN with separations on the order of Chandra’s resolution ($<1\arcsec$), receiving observations with sufficient counts, paired with a robust model of the Chandra PSF will allow for the most accurate analysis. In particular, we may expect that most dual AGN candidates should have separations $<1\arcsec$, as at a distance of 200 Mpc ($z\approx 0.05$) the physical-to-angular scale becomes 1.0 kpc/arcsec. Given the small number of currently confirmed dual and binary AGN, tools such as BAYMAX will be important for a precise measurement of the dual AGN rate, and as a result, an improved physical understanding of the evolution of SMBHs and their activity.

6 Conclusions

In this work, we present the first analysis by BAYMAX, a tool that uses a Bayesian framework to statistically and quantitatively determine whether a given observation is best described by one or two point sources. BAYMAX takes calibrated events from a Chandra observation and compares them to simulations based on single and dual point source models. BAYMAX determines the most likely model by the calculation of the Bayes factor, which represents the ratio of the plausibility of the observed data $D$, given the model $M_{j}$ and parameterized by the priors. We present the results of BAYMAX analyzing the lowest-mass dual AGN candidate SDSS J0914+0853, which was originally targeted as a dual AGN based on shallow archival Chandra imaging. The 15 ks exposure appears to have a secondary AGN $\sim$ 0$\farcs$3 from a primary AGN We received a new 50 ks Chandra exposure, with (i) a shorter frame time to avoid pileup and (ii) a different roll angle, with the aim of unambiguously determining the true accretion nature of the AGN. The main results and implications of this work can be summarized as follows:

Analyzing our new 50 ks observation, we find (by visual analysis) that the $2$–$7$ keV emission more closely resembles that of a single point source. Both spectra have an excess of flux at soft X-ray energies ($<1$ keV) with respect to the power-law continuum, while the 15 ks observation appear to catch the source in a higher flux state in the soft X-ray band. Both of these behaviors are seen in AGN with a “soft excess” component, and we fit our spectra with an absorbed redshifted powerlaw and blackbody (phabs$\times$zphabs$\times$(zpow + zbbody)). For the archival dataset, we find best-fit values for intrinsic $N_{H}=3.38^{+0.10}_{-0.10}\times 10^{20}$ cm*-2*, power law component $\Gamma=2.01^{+0.11}_{-0.12}$; and blackbody component $kT=0.10^{+0.03}_{-0.05}$ keV. For our new dataset, we find best-fit values for intrinsic $N_{H}=4.07^{+0.10}_{-0.10}\times 10^{20}$ cm*-2*, power law component $\Gamma=2.51^{+0.10}_{-0.12}$; and blackbody component $kT=0.11^{+0.05}_{-0.04}$ keV. 2. 2.

We find that the $2$–$7$ keV emission is consistent between the two observations, and fit the spectra in this energy-range with an absorbed redshifted powerlaw. For the 15 ks observations we calculate a total observed $2$–$10$ keV flux of $3.20^{+0.90}_{-0.80}\times 10^{-13}$ erg s*-1* cm*-2*, while for the 50 ks observation we calculate a total observed $2$–$10$ keV flux of $2.23^{+1.0}_{-0.49}\times 10^{-13}$ erg s*-1* cm*-2* s*-1*. This corresponds to a rest-frame $2$–$10$ keV luminosity of $1.83^{+0.31}_{-0.40}\times 10^{43}$ erg s*-1* and $1.25^{+0.35}_{-0.21}\times 10^{43}$ erg s*-1* at $z=0.14$ (assuming isotropic emission). 3. 3.

We use BAYMAX to analyze the 15 ks and 50 ks observations both individually, as well as combined, restricting our analysis to photons with energies between $2$–$7$ keV. Using BAYMAX we calculate a Bayes factor in favor of the single point source model of $\approx$ 6.5 and 9.8 for the 15 ks and 50 ks observations, respectively. When combining the two observations, we calculate a Bayes factor of 13.5 in favor of a single point souce. To test the impact of the MCMC nature of nested sampling, we run BAYMAX multiple times on the combined datasets. We find consistent results, with a spread in ln$BF$-space well-described by a Gaussian distribution centered at ln$BF=2.6$ with standard deviation of 0.2. 4. 4.

Our posterior distributions for both the single and dual point source model further support that SDSS J0914+0853 is a single AGN. Spatially, the best-fit locations from $\mu_{P}$ and $\mu_{S}$ are consistent with one-another within the 68% error level. Further, the joint posterior distribution has a shape expected from a single point source — the median values of the marginal posterior distributions are $r=0.15\pm 0.5$ and $\log{f}=-1.6\pm 0.4$, at very large flux ratios ($\log{f}\rightarrow{0}$) the separation is consistent with 0, and at very large separations ($r\rightarrow{2\arcsec}$) the flux ratio is consistent with 0. 5. 5.

We investigate the influence of our prior distributions, by running BAYMAX with Gaussian prior distributions for $\mu$, $\mu_{P}$, $\mu_{S}$. BAYMAX calculates a Bayes factor in favor of a single point source of $10.8\pm 1.2$, consistent within the errors of our initial measurement. Further, the posterior distributions returned by BAYMAX are consistent with those listed in Table 2. 6. 6.

We investigate how the Bayes factor determined by BAYMAX depends on the separation and flux ratio of a given dual AGN system. We find that for Chandra observations with at least 700 counts between $2$–$7$ keV, BAYMAX is capable of strongly favoring the correct model for most flux ratios when $r\geq 0\farcs{35}$. For the smallest separations ($r\leq 0\farcs{3}$), BAYMAX is capable of identifying the correct model when the flux ratio $f\geq 0.8$. 7. 7.

From a quasi-frequentist perspective, we estimate the probability of measuring a $BF\geq 13.5$ in favor a single AGN, using a null hypothesis that SDSS J0914+0853 is actually a dual AGN. Across all of parameter space ($0\farcs{3}<r<0\farcs{5}$ and $0.1<f<1.0$), we find $p\leq 0.05$ and can reject the null hypothesis at a 95% confidence level. Thus, the probability of SDSS J0914+0853 being a dual AGN system with a (1) flux ratio $f>0.3$, (2) separation $r>0\farcs{3}$, and (3) measured $BF=13.5$ in favor of a single AGN, is very low.

We have shown through various analyses that there is an absence of evidence supporting SDSS J0914+0853 as a dual AGN system. Specifically, BAYMAX estimates a Bayes factor strongly in favor of a single AGN and posterior distributions for a possible separation and flux ratio between a primary and secondary AGN are consistent with 0. Moving forward, statistical analyses with BAYMAX on Chandra observations of dual AGN candidates will be important for a robust measurement of the dual AGN rate across our visible universe. Lastly, our Bayesian framework will eventually be capable for more general analyses, such as evaluating binary active star systems.

A.F. and K.G. acknowledge support provided by the National Aeronautics and Space Administration through Chandra Award Number TM8-19007X issued by the Chandra X-ray Observatory Center, which is operated by the Smithsonian Astrophysical Observatory for and on behalf of the National Aeronautics Space Administration under contract NAS8-03060. We also acknowledge support provided by the National Aeronautics and Space Administration through Chandra Award Number GO7-18087X issued by the Chandra X-ray Observatory Center, which is operated by the Smithsonian Astrophysical Observatory for and on behalf of the National Aeronautics Space Administration under contract NAS8-03060. A.F. thanks Abderahmen Zoghbi for helpful discussion regarding PyMC3. This research has made use of NASA’s Astrophysics Data System.