The Host Galaxy of OJ 287 Revealed by Optical and Near-infrared Imaging

Spheroidal galaxies and galactic bulges harbor central supermassive black holes (SMBHs), and there is a fundamental tight relation between the SMBH mass and the bulge luminosity (mass), both in quiescent and active galaxies (Gültekin et al. 2009; Kormendy & Bender 2011; Sanghvi et al. 2014). Methods to measure the central BH mass include stellar and gas kinematics (e.g., Gebhardt & Thomas 2009), megamasers (e.g., Kuo et al. 2011), reverberation mapping (e.g., Kaspi et al. 2000), and emission line widths in quasars (e.g., Decarli et al. 2012; Sanghvi et al. 2014). Given that galaxy mergers are thought to be common, especially at high redshift (e.g., Somerville et al. 2008), it is reasonable to expect that double, even multiple BH systems form in galactic nuclei.

The BL Lacertae object OJ 287 (z = 0.306) presents one of the best cases for a small-separation SMBH binary system (e.g., Valtonen et al. 2008), other examples being SDSS J120136.02+300305.5 (Saxton et al. 2012; Liu et al. 2014) and possibly GSN 069 (Miniutti et al. 2019). Its historical light curve is unique among active galaxies, due to the prominent recurring double-peaked optical outbursts every ∼12 yr. Each outburst consists of two rapid brightenings, separated by ∼1 yr and lasting ∼6 weeks (Sillanpaa et al. 1996; Valtonen & Pihajoki 2013). A notable feature of the outbursts is that they are not strictly periodic, nor is the separation of the two peaks constant from one outburst to the next. This behavior is attributed to strong precession of the secondary's orbit in the model. We also note that the redshift of OJ 287 is secure, based on the detection of several clearly identifiable emission lines (e.g., Nilsson et al. 2010).

The binary BH model of OJ 287 is very successful in reproducing the previous outbursts, and it makes specific predictions about future outbursts. It is extremely unlikely that a single BH model would produce such good light-curve predictions as the binary BH model has done (Valtonen et al. 2011). In the binary BH model, the best-fitting orbital parameters correspond to the mass of the primary BH: M = 1.84 × 10¹⁰ M_⊙, and mass of the secondary BH: m = 1.30 × 10⁸ M_⊙. Notably, the primary SMBH is at the high mass end of the quasar BH mass function (e.g., Vestergaard et al. 2008; Portinari et al. 2012).

If OJ 287 contains an extremely massive primary BH, and if it follows the tight BH–bulge relation of active galaxies (e.g., Dotti et al. 2012), its host galaxy is also expected to be extremely luminous and massive. Puzzlingly, there have only been marginal or nondetections of the host galaxy of OJ 287 based on optical imaging (Benitez et al. 1996; Wurtz et al. 1996; Yanny et al. 1997; Heidt et al. 1999; Urry et al. 2000; Pursimo et al. 2002; Nilsson et al. 2003). The reported host luminosities in these studies are −22.0< M_R < −23.5, much fainter than expected from the BH–bulge relation, M_R ∼ −26). Converting the host galaxy absolute magnitudes into stellar mass by adopting the M/L of a single stellar population originated at z(burst) = 5 and passively evolving down to z = 0 (Decarli et al. 2010), these luminosities translate into host galaxy mass in the range of $\mathrm{log}({M}_{\mathrm{host}}/{M}_{\odot })$ = 11.4–11.9, much smaller than the expected host galaxy mass based on the BH–bulge relation, $\mathrm{log}({M}_{\mathrm{host}}/{M}_{\odot })\sim 13$. Conversely, based on the marginal host galaxy detections in the optical, if OJ 287 follows the BH–bulge relation, the primary BH mass is only $\mathrm{log}({M}_{\mathrm{BH}}/{M}_{\odot })$ = 8–9. Interestingly, this value is within the range of typical BH masses for BL Lac objects (Falomo et al. 2003; Plotkin et al. 2011), but much smaller than required in the binary SMBH model.

There is a reported near-infrared (NIR) detection of the host galaxy by Wright et al. (1998), who found for the host galaxy to be extremely bright, M_K = −28.9 ± 0.6. Transforming this to our cosmology we obtain M_K = −28.5 ± 0.6. With this luminosity, OJ 287 would fall exactly on the host galaxy luminosity—SMBH mass relation. However, the implied color of the host galaxy, R–K = 6.1 , would be much redder than in a typical elliptical galaxy (R–K = 2.7). We note that the data were taken with one of the first generation NIR imagers (IRCAM 62 × 58 px) and that they could not resolve the host galaxy in the R, J, or H bands with observations under similar seeing conditions and with similar exposure times as in the K band. Furthermore, the relatively poor spatial resolution (03–06 pixels and 09–13 seeing) and poor definition of the PSF severely limit their accuracy.

In this paper we present new observations, which enable us to detect the host galaxy of OJ 287 both in the optical and in near-IR. This enables us to constrain the host galaxy luminosity further and to shed light on some of the conflicting results obtained so far. In all calculations we use the cosmology H₀ = 70 km s⁻¹ Mpc⁻¹, Ω_Λ = 0.67, and Ω_M = 0.33.

OJ 287 went into a deep minimum in 2017 November, fading to R ∼ 15.8 (S. Zola 2017, private communication; Komossa et al. 2020), which is about 1 mag fainter than on average between the flares. We exploited the favorable AGN/host brightness ratio by obtaining a deep SDSS i-band image at the Gran Telescopio Canarias (GTC) on the night of 2017 December 2 when OJ 287 was already coming up from the minimum, but still about R ∼ 15.2. We obtained 77 exposures of 1.5 s with the OSIRIS instrument with a pixel scale of 0.254 arcsec/pix, totaling 115.5 s of exposure. The images were bias subtracted and flat-fielded with twilight flats and registered using stars in the field of view. The images were then coadded and deviations from a uniform background were removed by fitting a low-order polynomial to the background pixels and subtracted. The coadded image has an FWHM of 078. Calibration of this image was obtained through stars 15 and 16 in González-Pérez et al. (2001) using SDSS DR14 published magnitudes. The zero-points derived from these two stars differed by 0.01 mag.

K_s-band observations were obtained at the 2.5 m Nordic Optical Telescope (NOT), La Palma, during the night of 2012 December 18. At this epoch OJ 287 was also at a relatively low state at R ∼ 15.4. We used the 1024 × 1024 pixel NOTCam detector in imaging mode with a pixel scale of 0.235 arcsec px⁻¹, giving a field of view of ∼4 × 4 arcmin². The seeing during the observations was ∼08. The images were acquired by dithering the target across the array in a random grid within a box of ∼20'', and taking a 14.4 s exposure at each position. Individual exposures were then coadded to produce the final frame. A total of 237 exposures of 14.4 s were acquired, which provided a total exposure time of 3412 s.

Data reduction was performed using IRAF¹³ and followed the procedure described in Kotilainen et al. (2007) and Hyvönen et al. (2007). In each image, bad pixels were corrected for using a mask made from the ratio of two sky flats with different illumination levels. Sky subtraction was performed for each science image using a median averaged frame of all the other temporally close frames in a grid of eight exposures. Flat fielding was performed using normalized median averaged twilight sky frames with different illumination levels. Finally, images were aligned to subpixel accuracy using field stars as reference points and combined after removing spurious pixel values to obtain the final reduced coadded image.

For photometric calibration, we used stars 16 and 17 in González-Pérez et al. (2001), imaged simultaneously with the science observations. The zero-points derived from these two stars deviate 0.03 mag from each other, which is an ∼1σ difference considering the error bars in González-Pérez et al. (2001). A small systematic error may remain due to the deviation of the OJ 287 spectrum from a stellar spectrum and from the fact that González-Pérez et al. (2001) list K-band magnitudes and we used the K_s-band filter.

We analyzed the images by fitting two-dimensional surface brightness models to the observed light distribution of OJ 287. The models consisted of two components, an unresolved AGN nucleus and a host galaxy. The former was described by three parameters, x–y position and the flux in mJy, and the latter by four parameters, x–y position, flux, and effective radius. Given that the host galaxy is very weak and difficult to detect, we abandoned attempts to characterize it in more detail. Therefore, the host was assumed to be a early-type galaxy with Sérsic index equal to 4 and zero ellipticity. These assumptions are quite reasonable since BL Lac nuclei are almost exclusively found in giant ellipticals with Sérsic indices concentrated around n = 4 (e.g., Falomo et al. 2000; Nilsson et al. 2003). On the other hand, given the special nature of OJ 287, it may not be hosted by an "ordinary" galaxy. We will discuss this point later.

The faintness of the host galaxy and the inevitable errors in the PSF model make it very difficult to reliably establish any offset ≤FWHM between the host and the AGN nucleus. We thus assume that the AGN nucleus is centered on the host galaxy. The total number of parameters is therefore 5: AGN/host x, AGN/host y, AGN flux, host galaxy flux, and host galaxy effective radius.

Prior to fitting, we performed three critical steps: masking, sky background subtraction, and PSF construction. First, we masked out any pixels with visible contribution from targets other than OJ 287. Then we determined the sky background by measuring the average pixel brightness in ∼10 rectangular clean areas around OJ 287 and then computing their average. We then checked the obtained sky value by plotting line and column cuts near OJ 287, looking also for any slope in the background that might have remained from previous steps.

Obtaining a good PSF was complicated by the fact that OJ 287 was ∼2 times brighter than any star in the field of view in both bands. This means that using a stellar image directly as the PSF model does not work very well; the noise in the PSF image gets amplified when the PSF is scaled to the brightness of OJ 287 and both the fit and the residuals become very noisy. Analytical models are noise-free and therefore do not have this problem. However, stellar images give more accurate PSF models than analytical models, so they are preferred, if possible.

Since the PSF is often a function of the position on the detector, we constructed PSF models of several stars in the field of view. More specifically we used stars 15, 16, and 17 in the i band and 16 and 17 in the K band. We first fitted an analytical model consisting of two Moffat profiles to the stars. Subtracting this model revealed some asymmetric features plus the diffraction pattern from the secondary mirror support. These residuals were filtered using a wavelet-based multiscale filtering technique (Murtagh & Starck 1998), which recovers any significant signal in the residual image with high fidelity. This filtered image was then summed to the analytical model. In this way we obtain a low-noise PSF image, which accurately preserves the deviations from perfect symmetry.

The model fits were made with an affine-invariant Monte Carlo Markov Chain (MCMC) ensemble sampler (Goodman & Weare 2010). The posteriori probability distribution was first sampled using this sampler with multiple walkers. The number of walkers was equal to three times the number of free parameters. We used Gaussian priors that were not very restrictive, but nevertheless rejected unreasonable parameter values, like, for example, parameters that would place the model several pixels away from the observed position. The host galaxy priors we especially carefully constructed so that the priors would not dominate the fit over data. Typically 5000–10,000 iterations were made during each fit. After the sampler had finished, we checked that the success rate for most of the walkers was close to 50%. Then the marginalized histograms for each parameter were extracted and the best-fit parameter values and their errors were derived by fitting a Gaussian to these histograms.

The following noise model was used to compute the likelihood: a pixel with a brightness of N ADUs had an associated uncertainty

$$\begin{eqnarray}&&\sigma =\sqrt{{\sigma }_{N}^{2}+{\sigma }_{\mathrm{sky}}^{2}+{\sigma }_{\mathrm{PSF}}},\end{eqnarray} \tag{ 1 }$$

where

$$\begin{eqnarray}&&{\sigma }_{N}=\displaystyle \frac{\sqrt{{\left({gN}\right)}^{2}+{R}^{2}}}{{gs}}\end{eqnarray} \tag{ 2 }$$

is the photon and readout noise term, where g is the effective gain, R is the effective readout noise, and s is a constant (∼1.2 in this case) taking into account the smoothing caused by interpolation when the images were registered. The term

$$\begin{eqnarray}&&{\sigma }_{\mathrm{sky}}=\displaystyle \frac{{\sigma }_{\mathrm{bg}}}{c},\end{eqnarray} \tag{ 3 }$$

where σ_bg is the background rms and c is a constant parameterizing our estimate of the background accuracy. Here we use c = 10. The last term σ_PSF represents the PSF uncertainty and it is obtained from

$$\begin{eqnarray}&&{\sigma }_{\mathrm{PSF}}=\sqrt{{\sigma }_{\mathrm{sky}}^{2}+{\left(d\cdot N\right)}^{2}},\end{eqnarray} \tag{ 4 }$$

where d is a constant, whose value was determined during the fit process (see below). This PSF error model, while simplified, captures one crucial aspect, namely that far from the PSF center the sky uncertainty dominates the PSF error, whereas close to the PSF core, other factors dominate.

The fit proceeded as follows. First the value of d was determined in both bands by fitting a pure point-source model to a star in the field. We used stars 15 and 17 in the SDSS i band and star 17 in the K_s band. This fit had four free parameters, x, y, flux, and d. This fit yielded d = 0.05 and 0.06 for the SDSS i and K_s bands, respectively. Then we fitted OJ 287 with the five parameter model keeping the d value constant.

The results of model fitting are summarized in Table 1 and Figure 1. There is a small but significant excess over the PSF in both bands, visible in both the 1D profiles and in the 2D residuals. This excess is most visible in the lower portion of the middle panel, which shows the residuals after subtracting the scaled PSF. Since the host galaxy is finite in size, we expect these residuals to tend toward zero at distances of a few effective radii, about 3''–5'' in the case of OJ 287. In the K band the excess starts to turn down at r ∼ 4'', but in the r band there is no clear turning.

Zoom In Zoom Out Reset image size

However, none of the stars show any hint of an excess. Furthermore, far from the center, the error bars are dominated by the uncertainty of the sky level, which affects all points similarly. For instance, if the sky level is increased by 1σ, the residual plots would follow the lower ends of the error bars. Thus within 1σ sky uncertainty a turning down is possible in both bands. Thus we conclude that there is an excess of light around OJ 287 with respect to the stars in the vicinity, but this excess is rather weak and difficult to characterize.

Modeling the excess as a de Vaucouleurs profile, we obtain the values in Table 1. Two things can be read from this table. First, changing from one PSF star to another changes the host galaxy properties more than indicated by the error bars derived from the posteriori. This is especially true for the host galaxy magnitude and it is a clear indication that the dominant source of error is the PSF error.

Second, the host galaxy magnitude is less constrained in the K band. This is not due to the K-band observation having a less favorable AGN/host flux ratio. Both images were taken with similar AGN brightness and since the host galaxy is brighter in the K band than in the i band, the core/host flux ratio is more favorable in the K band. The reason for worse performance in the K band is most likely due to the fact that the PSF stars are relatively fainter in the K band than in the i band. Thus when the PSF is scaled to the brightness of OJ87, the inevitable errors in the PSF model are amplified more in the K band than in the i band.

For the discussion below we convert the magnitudes to fluxes and take the mean of the results in Table 1 in both bands. The error bars are set to the smallest number that keeps all results within the error bars. Thus in what follows we use f_host(i) = (140 ± 40)μJy and f_host(K) = (930 ± 410)μJy.

In Figure 2 we show the rest-frame broadband flux of the host galaxy as a function of wavelength. We also overplot an early-type galaxy spectrum at zero redshift. To transfer the data to the same redshift, we need to apply the K-correction K(z) and evolution correction e(z), both of which depend on the filter and galaxy type. We assume here that the host galaxy is an early-type galaxy and obtain K(z) from Chilingarian et al. (2010). The K-correction has almost no dependence on galaxy type in the K band and also in the SDSS r band, which we use as a proxy for the R band, the dependence is mild, except for the starburst galaxies. There appears to be very little star formation going on in OJ 287. Nilsson et al. (2010) measured a narrow ${\text{}}{H}_{\alpha }$ line flux of 4 × 10⁻¹⁶ erg s⁻¹ cm⁻², which transforms to a H_α line luminosity L(H_α) = 1.6 × 10⁴⁰ erg s⁻¹. Using the formula for the star-forming rate SFR = 4.6 × 10⁻⁴² L(H_α) M_⊙ yr⁻¹ (Twite et al. 2012), we obtain SFR = 0.07 M_⊙ yr⁻¹. This has to be considered as a generous upper limit since the narrow H_α line could be produced mostly by the AGN activity. We thus do not consider the starburst models in Chilingarian et al. (2010) and obtain K(z) = 0.45, 0.25, and −0.5 for the R, SDSS i, and K bands, respectively. For the evolution correction, we use e(z) = 0.2 for both bands. This value was determined by running the Pegase 3 code (Fioc & Rocca-Volmerange 2019) for a scenario with a single starburst 10 Gyr ago and passive evolution thereafter. The vertical scaling of the template in Figure 2 was determined by integrating the template over the respective bandpasses and finding the scaling factor that minimizes the chi-squared between our i- and K-band observations and the template fluxes.

Zoom In Zoom Out Reset image size

The observed fluxes are consistent with an early-type galaxy spectrum, corrected for passive evolution from z = 0.306 (Figure 2). The scaled template corresponds to a galaxy with M_R = −22.5, which is within the range of BL Lac host galaxy luminosities $\langle {M}_{R}\rangle$ = −22.8 ± 0.5 found by Sbarufatti et al. (2005). The K-band luminosity of the scaled template is M_K = −25.2 with (R–K) = 2.7. For the effective radii we find r_e = 5 ± 2 kpc in the SDSS i band and r_e = 4 ± 3 kpc in the K_s band. Both are consistent with sizes found for BL Lac host galaxies, although the errors are so high that we can only rule out a very large galaxy. Therefore, the host galaxy of OJ 287 does not look special in any particular way in comparison to BL Lac host galaxies in general.

We next study how our results relate to the SMBH mass $1.84\times {10}^{10}$ M_⊙ indicated by the binary black hole model. In Figure 3 we plot the host galaxy luminosity—SMBH mass data from Graham & Scott (2013) together with some recently found very massive SMBHs and OJ 287, assuming for the latter the mass indicated by the binary model and the luminosity of our scaled template. The line in Figure 3 shows the fit to the "core-Sérsic" galaxies of the Graham & Scott (2013) sample, i.e., essentially to galaxies brighter than M_K = −23.5. The SMBH in OJ 287 is overmassive by 1.13 dex, roughly 2.6 times the vertical rms scatter, 0.44 mag, in this plot, compared to SMBHs found in luminous nearby bulges studied by Graham & Scott (2013).

Zoom In Zoom Out Reset image size

The central SMBH in OJ 287 is thus at the high mass end of SMBHs with a similar host, but not significantly deviating from the distribution. The BH mass is basically a result of applying general relativity to the timing of the flares and it is very accurately determined (e.g., Valtonen et al. 2008). Furthermore, such overmassive SMBHs, challenging the coevolution between SMBHs and their host galaxy, have recently been found in an increasing number of galaxies, e.g., NGC 4486B (McConnell et al. 2011), NGC1277 (van den Bosch et al. 2012), NGC1281 (Yıldırım et al. 2016) , NGC1332 (Rusli et al. 2011), NGC 1271 (Walsh et al. 2015), and SDSS J151741.75-004217.6 (Lasker et al. 2013). However, to our knowledge no such detections have ever been made in supposedly spheroidal, massive host galaxies of bright AGN.

According to Ferré-Mateu et al. (2015), such galaxies form their SMBHs early (z ∼ 2) and then experience a lack of merger activity and subsequent galaxy growth, preventing them from reaching the present-day scaling relation. They are compact (R < 2 kpc) and consist of a purely old stellar population (t > 10 Gyr), and may be the local analogs (relics) of z ∼ 2 massive galaxies (e.g., Trujillo et al. 2007; van Dokkum et al. 2008). Interestingly, these galaxies, including OJ 287, are at the expected location of the M(BH)-Mbulge scaling relation at high redshift (e.g., Decarli et al. 2010; Portinari et al. 2012). Another way to produce an undermassive host galaxy is via tidal stripping (see e.g., Graham & Scott 2015), but this applies only to the case of stellar mass loss due to tidal interaction with a more massive galaxy, which is not the case for OJ 287.

Based on observations made with the Nordic Optical Telescope, operated by the Nordic Optical Telescope Scientific Association at the Observatorio del Roque de los Muchachos, La Palma, Spain, of the Instituto de Astrofisica de Canarias. J.K. acknowledges financial support from the Academy of Finland, grant 311438. A.J.C.T. acknowledges support from the Spanish Ministry Projects AYA2012-39727-C03-01 and 2015-71718R. I.H.P. acknowledges support from NRF 2018R1A2A1A05022685. S.Z. was supported by NCN grant No. 2018/29/B/ST9/01793. The work is partly based on the observations made with the Gran Telescopio Canarias (GTC), installed in the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofisica de Canarias, in the island of La Palma.