EVIDENCE FOR PERIODICITY IN 43 YEAR-LONG MONITORING OF NGC 5548

Because ULySS gives us a choice of defining and including components, we adjusted it to simultaneously analyze every component that contributes to the flux in the wavelength region around Hβ. For analyzing the variability of NGC 5548, we defined the model M(x) as follows:

$$\begin{eqnarray}M(x) & = & P(x)\ ([T(x)\otimes G(x)]+C(x)+N(x)\\ & & +B(x)+\displaystyle \sum _{i=1}^{4}{Fe}\,{{\rm{II}}}_{i}(x)+\displaystyle \sum _{j=1}^{n}{S}_{j}(x)),\end{eqnarray} \tag{ 1 }$$

where M(x), represents a bounded linear combination of nonlinear components—stellar template spectrum T(x) convolved with a line of sight velocity broadening function, an AGN continuum model C(x), a sum of narrow N(x), and semi-broad components B(x) of [O iii] emission lines, respectively—as well as a sum of Gaussian/Gauss-Hermit functions S(x) accounting for other AGN emission lines in the analyzed spectral domain and the Fe ii template consisting of four groups of Fe ii lines.

A multiplicative polynomial P(x), which represents a linear combination of Legendre polynomials, was included in a fit to eliminate the overall shape differences between the observed stellar and galactic spectra. The introduction of this polynomial ensures that the results are insensitive to the normalization, Galactic extinction, and the flux calibration of a galaxy and stellar template spectra (Koleva et al. 2008). For simplicity, we assumed a Gaussian velocity broadening function ${\boldsymbol{G}}({\boldsymbol{x}})$, but it is also possible to use the Gauss–Hermite polynomials (Rix & White 1992; van der Marel 1994). The contribution of the components to the total flux can be obtained from their weights, which are determined at each Levenberg–Marquardt (Marquardt 1963) iteration using a bounding value least-square method (Lawson & Hanson 1995).

For the stellar population model we used grid of PEGASE.HR single stellar populations, computed with the Elodie.3.1 library and a Salpeter IMF (Le Borgne et al. 2004). The model M(x) is generated at the same resolution and with the same sampling as the observation, and the fit is performed in the pixel space. The fitting procedure performs the Levenberg–Marquardt minimization (Marquardt 1963). In modeling the integrated spectra of NGC 5548 we added the power law to the stellar population base to represent an AGN featureless continuum (${f}_{\lambda }\sim {\lambda }^{\alpha }$). The spectral index α depends on the continuum slope, and represents the free parameter in the fit. To tie the parameters of the [O iii] lines, we defined two separate components of the model: a narrow component and a semi-broad component. In this way we tied the widths, shifts, and intensities of [O iii] components (the intensity ratio was kept to 3:1). The rest of the emission lines in the domain $\lambda \lambda [4430,6400]$ were fitted with a sum of Gaussians; we fitted He ii 4686 Å with two components, and used four components in the fit of Hβ line: narrow, broad blueshifted, broad redshifted, and very broad component. As proposed by Kovačević et al. (2010), we defined four groups of Fe ii lines to fit the Fe ii multiplets around the Hβ line: Fe ii s, p, and f group, and I Zw1 template. Because it is very difficult to fit the stellar population in the spectra of Sy1 galaxies with broad emission lines, in the first step we fitted the spectra of NGC 5548 in the minimum of activity with all free parameters of the model (defined with Equation (1)), and with multiplicative polynomials of the 15th order (as in Bon et al. 2014). We find that a single stellar population with an age of 7200 Myr and metallicity [Fe/H] = 0.2 fits the best spectra in the minimum. In fitting the rest of the spectra from monitoring campaigns of this galaxy, we fixed age and metallicity at best-fit values from the first step, and used the multiplicative polynomial of the first order to minimize the effects of this polynomial on the fit of emission lines. Free parameters associated with the stellar population were kinematic parameters—mean stellar velocity and dispersion. Examples of best-fit spectra at the minimum and maximum activity are presented in Figure 1.

Zoom In Zoom Out Reset image size

All the results presented below are based on the assumption that the luminosity of the [O iii] 5007 line, and all other narrow emission lines, does not change in time (see for example Peterson et al. 2002; Shapovalova et al. 2008). We assumed that the slit width and position angles used in the various campaigns were not affected by the O iii region size in the NGC 5548, which appeared to be very compact (see Schmitt et al. 2003). For that reason we made a simple test to estimate light loss for smallest slit used (see the Appendix). The spectra were scaled to a constant flux of F([O iii]λ5007) = 5.58 · 10⁻¹³ erg s⁻¹ cm⁻², and we constructed host galaxy subtracted continuum light curves of continuum flux measured at 5100 Å and the broad Hβ emission line. The 5100 Å continuum light curve and the Hβ light curve are presented in Figure 2.

Zoom In Zoom Out Reset image size

From the set of 1494 "IAW" spectra, we remove 128 spectra of the Hα region only, 109 spectra with poor spectral resolution, and 48 spectra that did not cover a large enough range to measure the 5100 Å continuum. This left 1209 IAW spectra, which combined with 249 spectra from Crimean Astrophysical observatory, 83 spectra from SAO/GHO observatories, and 7 spectra from Asiago, totaled to 1548 spectra. After careful inspection of the fits, we removed 10 additional spectra, leaving 1538 for further analysis (Table 2).

The zero-point error in the radial velocity was measured from the scatter of the difference between the peak, v_r of [O iii] λ5007 and the narrow component of Hβ, $\delta {v}_{r}={v}_{r}({\rm{H}}\beta )-{v}_{r}({\rm{O}}\,{\rm{III}})\approx 35$ km s⁻¹. We note that there is a known small systematic offset of $\delta {v}_{r}\approx +19$ km s⁻¹, frequently found in [O iii]-strong AGN (e.g., Eracleous & Halpern 2003; Hu et al. 2008; Komossa et al. 2008; Marziani et al. 2016, and references therein). The velocity zero point was set on the peak wavelength of [O iii] λ5007, because [O iii] λ5007 is a very sharp feature of high S/N and is less influenced by the underlying broad Hβ profile.

Typical errors on light and radial velocity curves were estimated from the dispersion of the measured parameters on short time intervals ($\lt 20$ days), in some cases multiple spectra were available (five for the continuum light curve). Because 20 d is a period that is long enough to include possible significant intrinsic variations, we considered the time behavior of the minimum dispersion value around the average of a given parameter (for example, FWHM, HWs) computed over 20 d as a function of time. Typical rms errors for Hβ seem to be around $5\cdot {10}^{-13}$ erg s⁻¹ cm⁻², while for the typical rms scatter for the 5100 Å continuum flux is 10⁻¹⁶ erg s⁻¹ cm⁻² A.

Li et al. (2016) discussed a large part of a similar dataset in a somewhat similar way. However, the main differences between Li et al. (2016) and our analysis was in our unified approach and the quantity of data used. In our study this was achieved using a more robust method and taking into account some important components that were not considered by Li et al. (2016; i.e., the galactic host emission, Fe ii multiplets, He emission, and two components of each [O iii] narrow line in spectral fit). As mentioned above, 1538 spectra were analyzed (with the new spectra covering the last 12 years), instead of only about 850 used by Li et al. (2016), with a large gaps in the last 12 years of their time series.

Light and radial velocity curves were analyzed for possible periodicity using standard methods, such as Lomb–Scargle (Lomb 1976; Scargle 1982), as well as the new method for unevenly sampled data tailored to the specific special conditions of our data series. In addition to light curves, we constructed curves using measurements of different fractional intensities of blue side (blue dots) and red side (red dots) Hβ broad emission line at 25%, 50%, 75%, and 90% as a function of time, which behave as half-width radial velocity curves at these fractional intensities. We also calculated line centroids,¹⁷ and full widths.¹⁸ Half widths and centroid radial velocity curves are presented in Figure 3. Full widths are also calculated, and their curves are presented in Figure 4. We applied different methods to test for possible periodicities using a similar method as for light curves.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Using the Lomb–Scargle (Lomb 1976; Scargle 1982) analysis (LS method), we analyzed light and radial velocity curves with previously removed linear trends. Results can be seen in Table 3.

To avoid very different sampling, we analyzed rebinned curves. This also eliminated shorter variations that correspond to the light-crossing timescale of the system. We rebinned the radial velocity curves to the 80-day average bins because the variability lags of Hβ to continuum variations were less than 30 days (Kaspi et al. 2000; Peterson et al. 2002; Zu et al. 2011; Koshida et al. 2014). According to Czerny et al. (1999), there are two clearly distinct physical mechanisms of variability in NGC 5548 curves with two different timescales: one less than 30 days and another with an order above 100 days. The short timescale variability is connected to the Comptonization of the soft photons emitted by the innermost part of the accretion disk, while in the long timescales the optical variability is not related to the X-rays (see more in Czerny et al. 1999). Therefore, we assumed that 80-day binning would be long enough to filter short variations and only analyze the longer ones. We found periodicities with very low false alarm probability in the radial velocity curves of half widths measured at 25%, 50%, 75%, and 90% of Hβ line maximum (see Table 3). The obtained periodicities show similar values. We also, searched for periodicity in the radial velocity curves of the full widths at 25%, 50%, and 75% of the maximum intensity of the broad Hβ emission line. The results are presented in Figure 5. There is a significant peak at about 3000 days, which is about half the value detected at half widths radial velocity curves (see Table 3). As argued below, these periodicities are far from being sinusoidal and cover a very small number of repeating periods (a little less then 3). We suspect that the standard LS analysis may not be reliable enough in such cases, and so searched for a more robust method that is more suitable for this type of data. The method is explained below and the numbers listed are based on the standard LS method and should be regarded as "tentative periods." We note that the LS method is known to give spuriously high significance levels to low-frequency periods for "red noise" variability (see Westman et al. 2011; Vaughan et al. 2016; N. Bon et al. 2016, in preparation).

Zoom In Zoom Out Reset image size

It is interesting to note that red and blue side radial velocity curves at 25% and 50% are anticorrelated, while red and blue radial velocity curves at 90% show a positive correlation (Figure 6). As we go toward the top of the line, the anti correlation switches to correlation (see Figure 6). This indicates that the peak of the line is shifting with the red side of the line, while at the base of the line radial velocity curves could be affected by two different kinematic components.

Zoom In Zoom Out Reset image size

To analyze the Hβ line shifts we also obtained a single Gaussian fit of the spectra and constructed a radial velocity curve from the obtained shifts. We fit a sine function, assuming the expected periodicity of 5700 days (see Figure 7), to lead the eye and not to claim a simple sinusoidal periodicity, because there are obvious deviations from the sinusoidal curve in several epochs. One can see some similarity of this radial velocity curve and the one obtained from measurements of red half width at 75% of the line maximum (see Figures 7 and 3), implying that the line shifts are mainly affected by variations on the red side of the line. Also, the similarity of these curves could indicate the same periodicity.

Zoom In Zoom Out Reset image size

As explained, our data are obtained from many different monitoring campaigns with very different sampling patterns. From the light curves one can see that the amount of data is much higher in the second third of the observed interval than the first and third parts; the latter contain less than a few hundred observations each, while the former contained above a thousand spectra. This distribution introduces various biases in standard methods for periodicity analysis of unevenly sampled data, based on sine function decomposition such as Fourier and Lomb–Scargle. The Lomb–Scargle method (LS) is commonly used to detect periodicities in unevenly sampled time series. LS follows the approach of Fourier methods and fits a sinusoid to the data. However, in classic Fourier methods uniform sampling gives rise to the orthogonality properties of the trigonometric functions, which has profound statistical implications. Thus, the Fourier series terms can be shown to be statistically independent under very broad conditions. This is not the case in LS unevenly sampled data. Nevertheless, it is a common practice to interpret the LS periodogram as a kind of Fourier series, relying on the assumption that the uneven sampling times are uniformly distributed. It is possible to neglect this effect if many cycles are included in the analyzed timespan. However, in our case the sampling has an obvious nonuniform nature, where the middle part of the time series is sampled much more densely than the rest. Moreover, we focus our attention on long periods, meaning only two to three cycles in total. Thus, the simplifying assumptions that allow the use of LS are no longer valid, and we must tailor a method to take into account the nature of the specific sampling pattern. Comparing Table 2 and Figure 5 shows that even though Lomb–Scargle analysis gives a very high probability of periodicities, the values of the obtained periods are dispersed. The periodicity result values span between 5500 and 6200 days. Very unevenly sampled data clearly require a method that can recognize repeating patterns, not only sine functions (as in cases of Fourier or Lomb–Scargle analysis). For that purpose we proposed method tailored to this specific case, with the special conditions of very few cycles and specific sampling characteristics.

We started from the working assumption that the data between JD = 47508 and JD = 52175 are well sampled and of better quality than the rest (see e.g., Figure 2). Those two problems prevent proper use of conventional period search techniques, because they usually assume many cycles and homogeneous data quality and sampling. Our new method treats the best-sampled part as fixed, and tests only the noise and sparsely sampled parts on either end to see how well they fit the hypothesis of a period. We first linearly interpolate the best-sampled part. Then we produce a periodogram by trying many periods. For each trial period we fold every point we test into the best-sampled part, and calculate its deviation from the interpolation (Similar to the phase dispersion minimization method; see Stellingwerf 1978; Plavchan et al. 2008). We sum those squared deviations and obtain a significance score that is similar to the χ² score for the periods. If there is a periodicity it should come up as a lower than the usual value of this score.

We look for a minimum (it is basically a ${\chi }^{2}$, but without normalization by the errors), and found it at a period of 5676 days (see Figure 8). But we still need to test how statistically significant the minimum is. To do this we consider the central segment as a reference we do not change, and test whether the points in the first and third segments fit the central part after phase folding. For this we carry out the following permutation test. In each iteration we randomly reshuffle the values of the points in the first and third segments (the measurements, not the sampling times). Thus we "ruin" their time dependency, but keep their "window function." We then recalculate the periodogram, and look for the minimum value. This determines the fraction of random reshufflings that will produce a value lower than the value obtained with the actual light curve. Only 2 out of 10,000 values were lower, which means the results have a significance p-value of around 0.0002.

Zoom In Zoom Out Reset image size

The Figure 8 shows the actual periodogram (in black) against the background of 1000 periodograms out of the 10,000. The structure of the problem imposes some structure on the periodograms, but the result stays significant because the 5676 periods the χ² trough look quite displaced from their randomly produced counterparts.

We repeated the analysis on the radial velocity curves shown in Figure 3, and show the results in the bottom panel of Figure 8. The periodicity looks very significant on the radial velocity curve measured at red side 75% of the line maximum; the lowest point (i.e., best fit) is at a period of 5725 days, which is essentially identical to the 5676 day period in the flux. As judged from the new method, none of the other radial velocity curves shown in Figure 3 shows a statistically significant periodicity.

It is important to emphasize that real radial velocity periodicities may be present in the Hβ profile, but the noisy half-width radial velocity curves do not allow us to detect them unambiguously.

If we assume that the periodic component is sine-like with a periodicity of 5676 days, the variance (see e.g., Nandra et al. 1997; Nikolajuk et al. 2004) of the data carried by the periodic component (with a secular trend) is about 19%. We consider both the variance in the data and the secular component under the assumption of long timescale (15.7 years) sinusoidal variations.

We can readily exclude two mechanisms that could give rise to a significant periodicity in a single BH system, as discussed in Bon et al. (2012): geodetic precession and disk-self warping induced by radiation pressure (Pringle 1996). The first of these occurs on timescales that are in general much longer than the observed periodicity (Begelman et al. 1980). As for the latter, Graham et al. (2015a) showed that warped disks around single BHs are not favored in an AGN context. In particular, the expected periodicity in such a scenario is much longer than the period found for NGC 5548 (e.g., a BH mass of 10⁸ M_⊙ gives a period between ${10}^{2.2-6.9}$ years; Pringle 1996; Graham et al. 2015b). Mechanisms that could produce orbital periodicity are reperturbation in a very large central disk that extends far beyond the inner accretion disk whose dimensions are known from recent reverberation mapping to be of order 2 ld (De Rosa et al. 2015; Edelson et al. 2015; Fausnaugh et al. 2016). A hot spot rotating at a distance of about 15 ld, which is not observed during a fraction of the orbit, can explain the periodic continuum variations. Emissivity perturbation by spiral arms in a large central disk can result in a double-peaked profile and emission line shape variability (see Chakrabarti & Wiita 1994; Jovanović et al. 2010; Lewis et al. 2010). Spiral arms can be triggered by the close passage of a massive object, such as a massive star cluster or super massive BH (see Lewis et al. 2010, and references within), by gravitational instabilities (see e.g., Flohic & Eracleous 2008, and references therein), or by an object passing through the extended disk (Chakrabarti & Wiita 1993). Fragmented spiral arms can account for emission line shape variation on relatively short time and radial velocity changes occurring on the dynamical timescale (Jovanović et al. 2010; Lewis et al. 2010). Such substructures in a nonuniform central disk could cause excess emission moving across the line profile (Jovanović et al. 2010; Lewis et al. 2010; Goosmann et al. 2014). We do not consider such models to be plausible explanations for NGC 5548 because of the huge central disk that is not observed in this source (the hot spot that, to explain the periodic variation should emit a sizable fraction of the luminosity of the small inner disk), and because the overall emission line spectrum of NGC 5548 is very similar to those observed in thousands of AGN of similar luminosity; these scenarios cannot explain all the population properties.

A TDE when a star is disrupted by a BH creates an appearance somewhat like an AGN for a limited duration. The tidal radius (i.e., the distance from the BH at which a star is tidally disrupted) can be written as ${r}_{t}\approx 1.5\times {10}^{13}{({M}_{\mathrm{BH},7})}^{\tfrac{1}{3}}{({M}_{\star ,\odot })}^{\beta }\,\mathrm{cm}$ (Hills 1975; Komossa 2015), where M_BH is the mass of the BH in units of 10⁷ M_⊙, and the mass of the star is in solar mass. The exponent β is ≈2/3 or 1/6 depending on whether the star is a main sequence star M_⋆, has mass ≲1 M_⊙, or belongs to the upper main sequence with M > 2 M_⊙ (Torres et al. 2010). For the NGC 5548 BH mass and for a main sequence star, the tidal radius is extremely small: a few gravitational radii. A main sequence star may therefore orbit at 2–15 ld without suffering a TDE. In the absence of a pre-existing accretion disk, tidal disruption causes a luminous flare with a short rise and a longer-lasting decline (e.g., Rees 1990), as observed in several cases (e.g., Komossa & Bade 1999; Gezari et al. 2012, and references therein). A partially stripped star in an orbit around the SMBH can cause repeated accretion events each time its orbit passes near pericenter, and may thus produce a semi-periodic signal in the light curve (e.g., Hayasaki et al. 2015). However, in the case of NGC 5548, this is an unlikely explanation if we were to interpret all its properties in the context of a TDE. For instance, its optical narrow emission lines imply the presence of a classical narrow-line region and therefore a much longer-lived AGN. Similar arguments hold for the possibility of causing a semi-periodic light curve from a TDE in a binary SMBH system (Liu et al. 2009). Another possibility is that we have a TDE in addition to the permanent accretion disk of a long-lived AGN, and the TDE contributes extra accretion during each pericenter passage. This would boost the accretion onto the BH in a periodic fashion, but material must be added to the accretion flow very close to the event horizon because the viscous time of accreting through the disk is very long. Perhaps such material is added to the central part, which periodically increases the X-ray emission from the disk corona; this in turn illuminates the central disk and boosts the optical light emitted from its surface. X-ray illumination is very important in NGC 5548 (Kaastra et al. 2014; Mehdipour et al. 2015, 2016), so such possibilities cannot be excluded. A detailed discussion of such a scenario is beyond the scope of this paper. We note, however, that TDEs are very rare events (Rees 1990; Magorrian & Tremaine 1999, one event per inactive galaxy every 10⁴–10⁵ years) even though rates can be higher in AGN (Karas & Šubr 2007) and in the presence of SMBBHs (e.g., Ivanov et al. 2005; Chen et al. 2011). In any case, the chances of seeing such an event in only one nearby galaxy analyzed are very small, and we therefore consider it very unlikely. Similarly, Landt et al. (2015) concluded that a TDE is a very unlikely explanation in NGC 5548.

There are several scenarios involving binary BH systems in AGNs. These can be divided into two broad groups: One is where only one of the BHs has an accretion disk and a BLR associated with it. The other is the case where both BHs are accreting through their own disks. For roughly equal BH mass (the only case considered here), the average separation of the two is of order 20 ld, and hence there is only one dusty toroidal structure around the two. In general, such systems are thought to be the end result of a galaxy merger, where the two BHs from the two galaxies are at the final stage of merging (see Milosavljević & Merritt 2001; Merritt & Milosavljević 2005). Earlier studies of NGC 5548 suggest some evidence for a merger 0.6–1.0 Gyr ago (more details in Tyson et al. 1998; Steenbrugge et al. 2005; Slavcheva-Mihova & Mihov 2011). Some simulations performed to characterize the SMBBH systems show the formation of a circumbinary disk, inside of which the two BHs are accreting matter and forming mini-accretion disks (e.g., Bogdanović et al. 2008, 2009; Hayasaki et al. 2008; MacFadyen & Milosavljević 2008; Cuadra et al. 2009; Smailagić & Bon 2015). Further out, the circumbinary disk cools and may form a torus. Blending BLRs has also been investigated (e.g., Shen & Loeb 2010; Popović 2012), or in case of a high-mass-ratio system, only one shifting BLR may be seen. A second BH can give rise to a host of phenomena that can yield periodic signals (see Katz 1997; Sillanpaa et al. 1988; Bogdanović et al. 2008; Bon et al. 2012; Kun et al. 2014; Graham et al. 2015a, 2015b). A detailed scenario of this type has been investigated in several papers by Bogdanović and collaborators (Bogdanović et al. 2008, 2009) who simulated a high-mass-ratio system with nearly identical time interval and periodicity as found here. These involve disk disruption, the formation and destruction of spiral arms in the gas between the BHs, and more. Some features of this model are appealing, especially those corresponding to periodic changes in the velocity curve of part of the gas, including cases where only one side of the profile is affected. Unfortunately, there have been no attempts to use the results of the dynamical simulation to calculated the resulting emission line spectrum, line profile, and time variations in the systems. While we are not in a position to look into this in detail, we note that the gas configuration in this model, and the gas properties, may be very different from what is known from many years of study of NGC 5548. Moreover, as noted earlier, the broadline spectrum of this source is very similar to the spectra of thousands of other type-I AGN. We thus consider it less plausible that NGC 5548 contains such a binary BH system.

A binary BH system where only one of the BHs carries its own disk and BLR (although the outskirts of the BLR must be disturbed by the "naked" BH) is perhaps easier to explain. Obscuration in this scenario is very inefficient. However, we note that gravitational lensing of the luminous disk around the primary BH by the second BH can enhance the continuum emission by a factor of order 1.11 over a period of a few hundred days; the effects over BLR emission would be even smaller, by a factor of less then 1%. For example, if both BHs have a mass of 5 × 10⁷ M_⊙ and their separation is 20 ld, the size of the Einstein ring is ≈0.48 ld. This should be compared with the size of the disk at 5100 Å (≈2 ld). Such an enhancement is achromatic, which gives the immediate prediction that other continuum wavelengths would show an identical change of amplitude during the passage. Such a scenario cannot explain the periodic Hβ variations or its periodic light curve. Moreover, 15 ld is well inside the BLR, so dynamical changes in the line-emitting gas must be considered too.

Using a radial velocity test for supermassive BBHs for broad, double-peaked emission lines (Liu et al. 2016) assuming equal mass components, the line peaks should be at about 2300 km s⁻¹. At some epochs, very small moving peaks in emission lines are identified corresponding to such velocities (see Shapovalova et al. 2004). However, it is very clear that the red and blue wings of Hβ respond to the same continuum variability at roughly the same time, so the gas emitting them is approximately at the same distance from the primary accretion disk. Recently, Li et al. (2016) proposed a BBH scenario in NGC 5548, which was a result of their two Gaussian decomposition model, fitted into 150 day averaged spectra. We tried to test these claims using a series of two Gaussian decomposition models with different types of constraints (e.g., a constant intensity ratio with a significant width difference of each Gaussian, forcing them to fit different parts of the line: narrower fits for the core and wider for the wings). We found that such a configuration could result in Gaussian components that switch sides and cross from blue to red and vice versa. Unfortunately, in every case tested (assuming constant initial parameters of component shifts), there were additional crossings in radial velocity curves that ruined the expected periodicity. We also note that in our modeling with a single Gaussian fit to the spectra (Figure 7), we obtained a radial velocity curve that appeared somewhat similar to the one of the 75% red half width (see Figures 7 and 3). This could imply the possibility of the same periodicity and could indicate the presence of a high-mass-ratio system. According to this analysis, we again find a BBH hypothesis to be unlikely, except for a in high-mass-ratio system, which we can neither disprove nor support.

This category includes several possibilities for obscuration of the central disk and part of the BLR by a moving object at a distance corresponding to a 15.7 year period. We consider two different possibilities: an object inside the BLR moving around the primary BH, and an object that is part of an outflowing wind moving around the polar axis of the disk. Both scenarios correspond to a situation where the length of obscuration is a small fraction of the period, perhaps a year or less. This situation is consistent with the new scheme presented in Section 4, which confirms the periodicity but does not show whether it is sinusoidal or corresponds to only one short event of dimming or enhancing the radiation of the central source. A large dust-free cloud moving around the primary BH, inside the BLR, at a distance corresponding to a period of 15.7 years, can cause the periodic obscuration of the central continuum source. This explanation is appealing because the required distance, about 15 ld, is exactly in the middle of the range of the multi-year RM size of the Hb line, and coincidental agreement between the two distances is unlikely. However, there are various difficulties to this scenario. The obscuring material must be thick enough and large enough to occult a large fraction of the central accretion disk, which is 2 ld in radius (see e.g., Fausnaugh et al. 2016). This dimension is larger than the typical size of BLR clouds that are of order 0.1 ld across, assuming that the density and column density are ≈10¹⁰ cm⁻³ and ≈10¹⁴ cm⁻², respectively (e.g., Netzer 2013). A dust-free cloud must be Compton thick to block the 5100 Å continuum because the ionized column of such gas is only of order ${10}^{22-23}\,{\mathrm{cm}}^{-2}$. Thus both the radial and lateral dimensions of such a cloud are orders of magnitude larger than those considered typical of the BLR (a collection of clumps adding up to the required dimensions is just as difficult to explain). Obscuration by a dusty gas cloud is easier to explain because a very small column density, corresponding to A_V < 1, is all that is required to absorb much of the radiation at 5100 Å. However, 15 ld is well within the dust sublimation radius for this source, and such grains will not survive in this environment. One can consider a very large dusty cloud in a spiraling elliptical orbit (e.g., Netzer & Marziani 2010) where, in this case, part of the dust is not sublimated because it is shielded from the central source radiation during a big part of the orbit. We did not explore this possibility in detail, but consider it problematic because efficient shielding of the dust from the radiation at wavelengths longer than the Lyman or Balmer continuum edges is hard to explain. Explaining the observed periodic variations in L(Hβ), and the period in its velocity curve, is even more challenging. Obscuration by dusty gas spiraling along the polar axis as part of a large scale disk wind is an alternative explanation. Disk winds have been considered for years as a general scenario to explain both the BLR and the dusty torus structure around the BH (see e.g., Elvis 2000; Elitzur & Shlosman 2006; Elitzur 2008; Czerny & Hryniewicz 2011; Netzer 2013, 2015; Arav et al. 2015; Elitzur & Netzer 2016). This geometry, which is sketched in Figure 9, allows a period of 15.7 years at a distance much larger than 15 ld from the BH, because material is rotating around the polar axis of the system with a small opening angle of about 30° (as suggested earlier for this object; see e.g., Rokaki et al. 1993; Kaastra et al. 2014) or less. For example, a dusty cloud at a distance of 50 ld, which is well inside the dust sublimation radius, can obscure the central source every 15.7 years for a few hundred days. The required column density is small, corresponding to A_V < 1 mag, and the lateral dimension is of order the disk size, ∼2 ld (while we talk about "a cloud," this may well be a collection of clumps). Such a cloud can also obscure part of the BLR, although this fraction is much smaller because of the much larger dimensions of the BLR. Obscuration by dusty material must result in wavelength dependent reddening of the central continuum. In principle, the obscuration may last only a few hundred days every 15.7 years, and the available spectroscopy is not good enough to exclude such short-term reddening events. Another difficulty with spiraling out material is that the outward motion of the gas will result in a continuous change of radius and periodicity; this change may be small over the 43 years of observations considered here. These details are beyond the scope of the present work. The geometry of the torus considered here is very different from the simple tori considered in earlier studies (Netzer 2015, and references therein). This opens a range of possibilities that may be related to the periodicity discussed here, in particular the obscuration by line of sight dust. Large structures in the nonuniform torus wall, scattering by the torus dust, and so on should be investigated in detail.

Zoom In Zoom Out Reset image size

The main difficulty in the previous scenarios is the lack of a clear connection between the periodic continuum variations, the periodic Hβ flux variations, and the periodic Hβ velocity variations. A more logical explanation is a mechanism that causes a periodic enhancement of the ionizing continuum, which in turn causes a periodic Hβ intensity variation (reverberation) that causes time-dependent changes in the mean emissivity radius of the Hβ line. In a virialized cloud system, an increase in the Hβ emission region is correlated to smaller gas velocities, which is in agreement with the observed line width variations. About 2/3 of the bolometric luminosity in NGC 5548 is due to far-UV and soft X-ray emission (Gaskell 2008). The X-ray source is probably illuminating the central disk causing much of the optical-UV continuum variations (e.g., Uttley et al. 2003; Edelson et al. 2015; Fausnaugh et al. 2016). Accretion events close to the primary BH, which enhance the emitted X-ray flux in a periodic fashion, can trigger the entire chain of events from optical continuum variations to Hβ flux variations to Hβ velocity variations. As suggested above, a pericenter passage of a partially tidally disrupted star is one possibility for such an event. An orbiting G2-like object (see Witzel et al. 2014) or some other stellar-size object with an orbit in a collision with the disk are other possibilities that could lead to such events. The energetics of these kinds of pericentric encounters depend on the consistency of the medium (the accretion disk or the torus) that is being crossed with their orbit. One can also imagine one or more gas clouds from the BLR in very eccentric orbit as the origin of such events. However, the mass of a single BLR cloud is tiny, which makes it very unstable against tidal forces, and its mass supply to the BH is many orders of magnitude below what is required to produce a significant X-ray flare.

Periodic variations could be caused by an orbiting object perturbing the accretion disk while passing through it (e.g., Syer et al. 1991; Chakrabarti & Wiita 1993; Armitage et al. 1996; Šubr & Karas 1999; Kieffer & Bogdanović 2016; Nguyen & Bogdanovic 2016; Pihajoki 2016). Such a scenario was proposed for NGC 5548 in Shapovalova et al. (2004) and was connected to the appearing and shifting bumpy features in the red wing of Hβ. If we assume a stellar-mass object passing through the disk at a radius smaller than 15 ld, it could cause a perturbation in the disk and produce shocks (e.g., Chakrabarti & Wiita 1993).

Such a collision could heat the disk (see e.g., Kieffer & Bogdanović 2016, raising the temperature above 10⁷ K), producing periodic optical and X-ray emission. It is not our intention to attempt such calculations for NGC 5548, only to mention that the disk dimension and the orbit eccentricity are likely to be the limiting factors in such cases. The hot spots caused by impacts may not need to be present for the entire orbital time, but would be made periodically with each collision. Knowing that in our Galaxy there are a number of central stars that are on highly eccentric, randomly inclined orbits (i.e., the S0 stars; Eckart & Genzel 1997; Ghez et al. 1998; Gillessen et al. 2009), and some show similar periodicities (e.g., for S0-2 the periodicity is 15.2 years and for S0-14 it is 38 years, Zucker et al. 2006; for S0-102 it is 11.5 years, see Meyer et al. 2012), it is conceivable that one could find objects on inclined orbits that could cross the accretion disk of NGC 5548. In the case of a star passing through the part of the disk responsible for the optical emission at about 2 ld radius, the eccentricity of the orbit would be about 0.7.