Chandra Observations of Spikey: A Possible Self-lensing Supermassive Black Hole Binary System

This work examines the recent X-ray observations of the active galactic nucleus KIC 11606854 (nicknamed “Spikey”) by the Chandra space telescope. Based on previous observations of a symmetric flare in the system’s light curve by the Kepler space telescope, Spikey has been proposed to be a self-lensing supermassive black hole binary system in which the more massive black hole gravitationally lenses the accretion flow of its smaller companion. The recent Chandra observations (2020 March–May) correspond to the time when the next pulse was expected to occur and were separated in enough time to observe the apparent relativistic Doppler boosting effect from the high orbital velocities of the black holes. We model the expected self-lensing plus Doppler boosting light curve using our wavelength-dependent extended source self-lensing model combined with our relativistic orbital motion code. This orbital motion code is capable of modeling the expected apsidal precession for Spikey, which can be used to predict future pulses. We show that the expected signal was undetectable in the Chandra data as the intrinsic X-ray variability associated with the system was large relative to the changes expected by self-lensing and Doppler boosting. Expected flux increases in more favorable wavelengths were also calculated using our wavelength-dependent self-lensing model, revealing a relationship between the observing wavelength and measured orbital inclination.


Introduction
analyzed the light curves of 21 type 1 active galactic nuclei (AGNs) located within the field of view of the Kepler space telescope. The AGN KIC 11606854 possessed two rather peculiar features in its light curve: the quasar exhibited an almost sinusoidal modulation in its relative flux as well as a very sharp pulse that lasted roughly 15 days. Smith et al. (2018) attempted to fit a t −5/3 model for a possible tidal disruption event but were unable to find a good fit to the data. The phenomenon remained unexplained until Hu et al. (2020) proposed that this flare was the result of a gravitational lensing event between two supermassive black holes (SMBHs) in orbit around one another. The pulse would correspond to the more massive black hole lensing the accretion flow of its smaller companion during eclipsing events, resulting in a large flux increase in the system's light curve, whereas the modulation would be due to the relativistic Doppler boosting effect brought about by high orbital velocities. Given that the smaller (secondary) black hole has a higher orbital velocity than its larger counterpart (roughly 4% the speed of light), Doppler boosting is driven by the secondary's motion, where the apparent magnitude of the system decreases as the secondary moves away from the observer and increases as it moves toward the observer.
A lensing event that occurs between two components in a binary system is called "self-lensing" and, historically, all selflensing events that have been observed come from white dwarf stars lensing their sunlike companions (Kruse & Agol 2014;Masuda et al. 2019). From the maximum flux increase observed in the light curve, as well as the timescale of the event, one can determine the mass of the compact object. This method of mass determination has been discussed in the literature (Maeder 1973;Witt & Mao 1994;Marsh 2001;Rahvar et al. 2011), although the focus was primarily on stellar mass lensing events. While a few studies (D'Orazio & Di Stefano 2018; Ingram et al. 2021) have modeled how possible SMBH self-lensing systems would manifest themselves, until now, there have been no observed self-lensing events or known candidate black holes of any mass. This would make Spikey (a nickname given to the system by Hu et al. 2020) a unique source of self-lensing and potential discoveries regarding the accretion process and evolution of such a system would be extensive.
A self-lensing event in a system such as Spikey should occur once per orbit. Due to the Doppler boosting effect, one can determine the orbital period of the system in the observer's frame of reference (Spikey is located at a redshift of z = 0.926). Hu et al. (2020) estimated a period of 805 days, corresponding to a rest-frame orbital period of 418 days. Since the Kepler pulse data from 2011, (so three pulses could have occurred by now), there has been no publicly released data with a high enough signal-to-noise ratio to possibly detect the pulse. This lasted until 2020, when the Chandra space telescope carried out several observations of the system before, during, and after the anticipated pulse. We analyze the data set from the Chandra observations, estimate the timing of the pulse, and we discuss the modeled X-ray signal from these observations due to the Doppler boosting and self-lensing effects.

Observations and Data Analysis
The source was observed 5 times by Chandra/ACIS-S during 2020 March-May, each time for nearly 20 ks ( Table 1). The raw data were reduced using CIAO-4.13, following standard processing trends including chandra_reprocess. The reprocessing takes care of bad pixel and GTI corrections, while barycenter corrections were performed using the axbary command. Source extraction was performed using a 2 5 circular region and the background was subtracted using a 2 5 inner and 5″ outer radius annular region. Our main motivation was to look for flux enhancement as previously observed in the Kepler data for this source (Smith et al. 2018). In this pursuit, we extracted all the light curves using dmextract, in four energy bands (broad band B = 0.5-7 keV, hard band H = 2-7 keV, medium M = 1.2-7 keV, and soft S = 0.5-1.2 keV). Source flux and mean count rates for all the observations were calculated using the srcflux command in broad band.
Source and background spectra with region-specific response files were created using the specextract command for each observation. Finally, the spectral data and response files were fitted in Sherpa.

Light Curves: Flux and Hardness Ratio Variability
The source has been detected in all five observations, and the resulting long-term light curve (binned at 4 ks) is plotted in Figure 1. The source flux varies within and between all the observations, so we created additional light curves binned at 1 ks to look for intraobservation variability (Figure 2). The maximum count rate appears during the first observation (ObsID 22557) and decreases with each subsequent data set. Furthermore, in the last three observations, the interobservation variability is less than that of the first two observations. The hardness ratio evolution is shown in Figure 3. We observe the hardness ratio decrease in the second observation and an increase in the third observation before decreasing again for the remainder of the observations. This sudden increase in the hardness ratio during the third observation is interesting given that this is where the projected self-lensing flare was expected to occur. This ratio shows an anticorrelation with X-ray luminosity, which decreases continuously throughout the observations (Figures 1 and 3).
Neil Gehrels Swift observations analyzed by Hu et al. (2020) showed decreasing flux measurements that were consistent with the overall trend caused by Doppler boosting that was expected during that time frame. However, in the observations made by Chandra, we find that the flux is not only decreasing when it was expected to be increasing due to Doppler boosting, but it is decreasing relative to the Swift observations made in    58,913.32, 58,951.89, 58,964.59, 58,975.63, and 58,997.30, respectively. the prior year. The dimmest Swift observation measured a flux of 4.5 ± 0.9 × 10 −13 erg cm −2 s −1 in the 0.5-0.7 keV energy range, whereas the brightest Chandra observation found a flux of 3.303 ± 0.254 × 10 −13 erg cm −2 s −1 . This points to a longterm secular trend of decreasing flux in the system.

Spectral Changes
ACIS-S pulse height spectra with sufficient counts were used to detect spectral changes over the five observations. The spectra were binned in such a way that they had at least 30 counts per bin and the PI channels above 5.5 keV were ignored. Then, an absorbed power-law model was used to fit all of the observations. The results of the spectral fitting are listed in Table 2, and the energy spectra are shown in Figure 4.
It was expected that, due to Doppler boosting increasing throughout the five observations (Figure 9), the slope of the energy spectrum would flatten with each observation, resulting in a lower photon index. In contrast, we see that the fifth observation has the steepest slope and, thus, the highest photon index of any observation. This is further evidence that spectral variability overwhelms any effect expected to result from Doppler boosting.

X-Ray Variability in PSD Profiles
The power spectrum density (PSD) profiles represent the variability distribution as a function of frequency. The PSDs were generated from light curves binned at 6 s using the standard process (Priestley 1981;Vaughan et al. 2003;Ingram & Done 2011). The periodograms were created using the powspec tool in heasoft (v6.29). We fitted the PSDs in frequencies 0.003 Hz with a power-law model using the X-ray spectral fitting packagexspec (v12.12) package (following the procedure described in Appendix 1 of Ingram & Done 2011). All of the periodograms were normalized throughout using (rms/mean) 2 (Vaughan et al. 2003).
The power-law indices have been used for comparison between different observations in each band. We have used a high-frequency range where the red noise decreases steeply compared to the low-frequency counterparts (Uttley et al. 2002). The results from this fitting are reported in Table 3, and the PSD profiles are provided in Figure A1 of Appendix A.

Extended Source Self-lensing Model
In order to predict the profile of a self-lensing pulse in an eccentric system such as Spikey, one would need to be able to keep track of the relative orbital velocities between the two  objects as well as their binary separation at any given time. In order to do this, we developed a code in the program Mathematica that uses a system of differential equations derived from the Paczynski-Wiita potential (Paczyńsky & Wiita 1980) that models the orbits of these two black holes. We chose this potential due to its accuracy in modeling the relativistic orbital precession expected in Spikey. Four differential equations are used to keep track of the black holes' positions and acceleration in Cartesian coordinates as a function of time. The four equations are as follows: and R s is the total binary mass Schwarzschild radius. The double primed variables are the accelerations in the x and y directions of the more massive and less massive black holes (labeled one and two, respectively) and the unprimed variables represent their positions. This system of differential equations is solved using Mathematica's built-in numerical differential equation solver called "NDSolveValue," which returned four interpolating functions that allow one to track the black holes' motion as a function of time. Using the binary parameters estimated in Hu et al. (2020), we have plotted out the orbits of the two black holes over the course of 42 yr in the observer's frame (see Figure 5). We can see the orbital precession over this time frame, which results in the predicted self-lensing pulse occurring ≈2.41 days earlier per orbit as observed on Earth. Hu et al. (2020) modeled the self-lensing pulse from Spikey using a point-source model where the smaller black hole's accretion disk (referred to as the source) is treated as a single point. This represents a reasonable approximation depending on the size of the source accretion disk. However, it has been noted that the apparent size of the source accretion disk is highly wavelength dependent; in particular, the accretion disk will appear smaller at shorter wavelengths and larger at longer wavelengths. This is because material flowing closer to the black hole is much hotter due to the increased orbital speeds and, thus, increased friction (Lynden-Bell 1969;Shakura & Sunyaev 1973). This implies that the effective size of the source accretion disk will vary depending on the wavelength the system has been observed. The choice of a point source effectively eliminates the chromatic dependency of the selflensing pulse due to the change in apparent source size. Therefore, we have decided to abandon the point-source model in favor of an extended source model where we treat the source accretion disk as an object of finite, wavelength-dependent size.
In an extended source self-lensing model, one of the most important parameters for determining the expected amplification is the value P å , the ratio between the radius of the source (the background accretion disk) and the radius of the lens (known as the Einstein radius). We adopt the P å equation from D' Orazio & Di Stefano (2018) to estimate the ratio between the radius of the source accretion disk and the radius of the lens as a function of the emitted wavelength. The equation is as Note. The model used here was a single power law. The PSDs for each ObsID and each band are provided in Appendix A, Figure A1.
where η is the accretion efficiency factor, ò is the disk accretion rate in units of the Eddington rate, n a is the binary separation in units of the Schwarzschild radius of the total mass, and M is the total binary mass. The final term was added to approximate the surface area of the accretion disk that is visible when the binary system and accretion disk are inclined by f and J degrees, respectively. Figure 6 helps break down the geometry of this self-lensing system. It is worth noting that the specification of λ emitted instead of λ is necessary for Spikey given that, due to its location at a redshift of z = 0.926, the maximum wavelength emitted will not be the same as the wavelength observed. The emitted wavelength can be calculated using In order to estimate the amplification during the lensing event, one must also compute the angular separation in units of the angular Einstein radius as a function of time with the equation given that where u 0 is the closest angular approach, a is the orbital separation at eclipse, and f is the inclination of the orbit with respect to the lens-observer line of sight (Rahvar et al. 2011).
The variable t E represents the Einstein time, the time it takes the lens to move the length of its Einstein radius (R E ). The Einstein radius is found from the equation Next, we substitute the values of u(t) and P å into the equation for amplification as a function of time: This analytical solution, originally solved by Witt & Mao (1994), plots the relative flux of a lensing event for a circular source of constant surface brightness. Now, an accretion disk does not have constant surface brightness, and it will only appear circular when viewed directly face on. This is where the addition of that final term in the P å equation comes into play.
The final term will adjust the corresponding P å value depending on the inclination of the accretion disk J. Assuming the orbital inclination is f = 0°, a disk inclination of J = 90°r esults in this term vanishing, leaving P å unchanged. If the value of J is smaller, it causes the value of P å to shrink as less of the disk is visible. At J = 0°, one finds that P å = 0, resulting in a point source. We adopted the same Doppler boosting model from that of Hu et al. (2020) given by Figure 6. The geometry of a self-lensing event in an SMBH binary system seen from two vantage points. On the left, a face-on view of a lensing event where the foreground black hole, with the blue accretion disk, gravitationally lenses the background companion. The surface area of the background black hole's accretion disk is increased (a secondary image forms as well), resulting in an increase in its apparent magnitude. The right panel gives a side view of this system where the orbital inclination and the inclination of the background accretion disk are more apparent. This figure was created using SpaceEngine (Romanyuk 2016). The software allows for real-time simulations of gravitational lensing using the potential derived by Paczyńsky & Wiita (1980).
where α is the spectral index (measured through observation), v is the orbital velocity of the secondary black hole, f is the orbital phase, and I is the orbital inclination.

Modeling the Spikey System
From the profile of the pulse observed by Kepler (and an orbital period of T = 418 days), Hu et al. (2020) determined the total mass of the pair to be 3 × 10 7 M e , the inclination angle i = 82°, and the eccentricity e = 0.524. We chose to model the expected X-ray pulse using these values and with an observed wavelength of λ = 6 nm. Hu et al. (2020) also reported an average spectral index of α = 0.6 for Spikey as observed by the Neil Gehrels Swift Observatory from two observations. This matched what was seen during the Chandra observations (the range we measured from the energy spectra was α = 0.6 − 0.8), so we chose to use α = 0.6 in our model. In Figure 7, we compare the expected Chandra light curve with that of our model Kepler light curve. There are two immediate features that are noticeable. First, the expected flux increase, as caused by Doppler boosting, is higher for the Chandra model than what was simulated for the Kepler model (a 6% flux increase compared to a 2% increase). Second, the self-lensing pulse of the Chandra model is smaller than that of its Kepler counterpart (a 5% increase compared with a 7% increase). It was expected that the smaller the apparent size of the source accretion disk, the larger the amplification. Given that the disk would appear smaller in the shorter X-ray wavelength, one may have expected the Chandra observation to produce a larger flare during the lensing event. However, in the case of an inclined system such as Spikey, it is possible for longer-wavelength observations to result in larger amplifications. This can be seen in Figure 8, where the P å value corresponding to the maximum amplification of a pulse increases with increasing orbital inclination.

Discussion
We have plotted the Chandra data as compared to our model's light curve in Figure 9. The X-ray flare was expected to be rather similar to that observed by the Kepler telescope, but its absence is not entirely surprising. For the third Chandra observation, we see a relative flux dispersion of 0.23, a full 0.17 larger than the expected pulse itself. Comparing the Chandra data with those of the Kepler observation, one can see the large difference in dispersion between the two data sets. The relative flux dispersion in the Kepler data during its pulse is only 0.01. It is known that X-ray variability in quasars is high relative to other wavelengths, and the presumption that there are here two SMBHs in orbit around one another would likely increase this variability even more (Graham et al. 2015). While the lack of a flare is disappointing, it does not conflict with the self-lensing explanation. Radio observations by Kun et al. (2020) did support the original prediction for a self-lensing system (although no pulse would be expected at such radio wavelengths). Spikey will have to be observed again in a wavelength in which variability is lower (preferably in the infrared).
We showed in Section 4 that shorter wavelengths do not necessarily correspond to higher amplification. Assuming the estimated orbital inclination of f = 8°, an observation in the infrared (≈3300 nm) should correspond to a ≈20% relative flux increase, the maximum such increase for the system given  Maximum amplification for four different P å values as a function of the lens-centered orbital inclination f. As expected, when the inclination is close to zero, a point source will produce the maximum flux increase. However, it can be seen that as inclination increases, larger P å values will begin to produce higher magnifications. Figure 9. The Chandra data as compared with the expected light curve from our self-lensing plus Doppler boosting model (SL+DB). The purple dashed curve indicates the location of the pulse if orbital precession were not taken into consideration. The black dots represent the Kepler data, which has been shifted in time to align with the projected Chandra lensing event. The intrinsic X-ray variability in the system is too high and would have masked any self-lensing pulse that may have occurred. the estimated binary parameters. Even with some variability, a 20% increase could lead to a significant detection. In fact, multiple wavelength observations of self-lensing pulses can be used to better constrain the orbital inclination of the system. Figure 10 shows how the maximum flux increases measured at each observing wavelength give very distinct curves depending on the inclination of the system. It also shows, unexpectedly, how longer wavelengths can correspond to higher flux increases. In edge-on systems (f = 0°), a point source would result in the highest amplification possible for the system. In short, the smaller the source accretion disk, the higher the amplification. However, in inclined systems such as Spikey, observations in longer wavelengths (thus larger source disks) can lead to higher amplifications than their smaller accretion disk counterparts. In an inclined system, more of a larger disk would be passing through the midpoint of the lens than in the case of a point source, resulting in higher amplification. While this result defied initial expectations, it has been shown in previous studies that this is, indeed, possible for systems given appropriate conditions (Heyrovský & Loeb 1997). Witt & Mao (1994) showed that, in general, larger source sizes will result in larger amplification in the cases of P å  2u and lower amplifications when P å  2u. Nevertheless, the size of the accretion disk will be limited by tidal truncation, resulting in a wavelength of peak emission (Papaloizou & Pringle 1977). This effectively means there is a limiting wavelength where one can increase the apparent size of the accretion disk. One can approximate this wavelength of peak emission using the relation described in D'Orazio & Di Stefano In the case of Spikey, the emitted wavelength at the outermost edge of the source accretion disk is 815 nm corresponding to an observed wavelength of 1569 nm. As we can see in Figure 10, this wavelength results in a modest amplification, but an increase nonetheless (10% instead of 8%). It may also be worth investigating pulses at longer wavelengths ∼2000-3000 nm as a test of the model's results. An SMBH binary such as Spikey is rare, but with the aid of self-lensing, it may possible to test the interaction of accreted material between the two black holes and how this accretion process changes with time.

Conclusion
A self-lensing pulse from Spikey (the putative self-lensing binary SMBH candidate) in the soft X-ray band should result in a 5% flux increase lasting a few days due to the pulse itself, accompanied by a slower 6% modulation at the orbital period due to Doppler boosting, similar to what was seen by Kepler at visible wavelengths. We analyzed a series of 5 × 20 ks Chandra observations evenly spaced during a 3 month window scheduled to include a self-lensing pulse predicted for 2020 March-May. The third observation is closest in time to the predicted peak of the lensing pulse, which was projected using our Paczysnki-Wiita orbital motion code. We also examined published Swift results (Hu et al. 2020) from a previous epoch that were reported to show a decline consistent with the Doppler boosting timescale. We characterized the timing and spectral properties of the X-ray data, finding that flux variability on both inter-(weeks) and intraobservation (<20 ks) timescales exceeds the magnitude of the anticipated lensing and Doppler boosting signals. We do find a statistically significant continuous long-term decline in flux. The decline in flux seen by Swift has continued so that an alternative explanation is required, one that does not involve Doppler boosting. Subtracting a quadratic fit to the declining Chandra fluxes reveals a dip in the residual of the third observation, where the pulse was expected. The size of this dip is 4.35% below the quadratic fit to the nonlensed observations (Appendix B), which is inconsistent with predictions of the lensing pulse's height ( Figure B1). Simultaneously, this residual is consistent with the stochastic variability seen throughout the data set.
Doppler boosting is known to produce shifts in color (or apparent temperature) of a source, and a shift in the observed variability timescale. Accordingly, we examined the energy spectra ( Figure 4) for changes in photon index and the powerdensity spectra ( Figure A1) for changes in shot-noise slope. Small changes are found in both parameters (Tables 2 and 3), but no coherent behavior is seen. We conclude that flux variability in soft X-ray observations of Spikey is dominated by variability inherent to the accretion process. Future efforts to place meaningful constraints on Doppler boosting and selflensing in Spikey and other binary SMBHs at X-ray wavelengths will require higher effective area (to suppress Poisson noise to the level required) and higher observing cadence in order to pick the lensing pulse from the shot-noise and secular variations (inherent variability on longer timescales) in AGN light curves.
It was shown that observations at longer wavelengths can result in higher possible flux increases and that multiwavelength observations can be used to test theories related to the tidal truncation of accretion disks in such SMBH binary systems. Observations in longer wavelengths may reveal the expected pulse or, if there is a lack of a pulse, open the door to other explanations of the flare observed by Kepler. Figure 10. Maximum amplification for three different inclination angles as a function of the observed wavelength. By varying the value of P å through multiple wavelength observations of Spikey, one could determine the inclination of the orbit based on the shape of the corresponding curve. The largest flux differential occurs in the infrared. However, tidal truncation would likely cut off the source accretion disk at an observing wavelength of λ = 1569 nm (denoted with the red dashed line). Detrending the Chandra Observations The five observations made by Chandra showed a general downward trend in the count rate with each observation. We decided to remove this trend from the data set by applying a fit to the mean count rate of each observation and subtracting the model from the data. This was done in an attempt to search for evidence of a pulse in the third observation once the trend was subtracted. We created two fits, a linear and a quadratic fit, both of which were created using Mathematica's built-in Fit function. We used the mean from four of the five observations, ignoring the third observation as this is where the self-lensing pulse was projected to occur. Our fits, as well as the resulting residuals from the detrending process, are shown in Figure B1.
In the case of the residuals produced by the linear fit, there does not appear to be any trend consistent with Doppler boosting or self-lensing. Interestingly, for the residuals produced from the quadratic fit, the third observation does indeed show the greatest deviation from the model curve, though its deviation is negative. This is the opposite of what would have been expected from a self-lensing pulse, but it is still consistent with the stochastic variability seen throughout the data set.
To compare the observed deviation from the quadratic model in the third observation with that of self-lensing, we took the model as equivalent to the unlensed baseline flux level in Figure 7 and then we calculated the relative deviation as (CR [3]-Model[3])/Model[3]), which yields −4.35%. The same procedure was repeated using the flux inferred from the spectral fits, resulting in a relative deviation of −9.63% ± 5.20%. Figure B1. The top panel showcases the Chandra observations of Spikey with linear (purple) and quadratic (red) fits through the data. These models were subtracted from the mean count rate of each observation to remove the downward trend observed over the five observations. The fit does not include the third observation (black) as this is where the self-lensing pulse was projected to occur. The bottom panel shows the residuals from this detrending process.