A novel approach to study the variability of NGC 5548

Understanding the properties of the continuum radiation and broad emission lines of active galactic nuclei provide significant information not only to model the radiation mechanism and constrain the geometry and kinematics of the broad-line region (BLR) but also to probe the central engine of the sources. Here we investigate the multifractal behaviour of the H$\beta$ emission line and the 5100 {\AA} continuum flux light curves of NGC 5548. The aim is to search for multi-scaling signatures in the light curves and check if there is a possible nonlinear relationship between them. To this end, we use a multifractality analysis technique called Multifractal Detrended Moving Average (MFDMA) analysis. We detect multifractal (nonlinear) signatures in the full monitoring and densely sampled period of the H$\beta$ line and 5100 {\AA} continuum light curves of NGC 5548, possibly indicating the presence of complex and nonlinear interaction in the 5100 {\AA} continuum and H$\beta$ emission line regions. Moreover, the degree of multifractality of H$\beta$ line is found to be about twice that of the 5100 {\AA} continuum. The nonlinearity of both emissions could be generated when the broad-line region reprocesses the radiation from the central compact source. Finally, we found that anti-persistent long-range temporal correlation is the main source of the multifractality detected in both light curves.


Introduction
In the study of active galactic nuclei (hereinafter AGNs), understanding the nature of the central engine, the geometry and kinematics of the broad-line region (hereinafter BLR), and the connection between the continuum and the BLR remain the most important questions. Several studies to uncover the working mechanism of the central engine (e.g., Netzer & Shields 2007;Carol et al. 2008;Neronov et al. 2008), to investigate the variability and constrain the geometry and kinematics of the BLR (e.g., Chiang & Murray 1996;Sulentic et al. 2000;Denney et al. 2009;Gaskell 2009;Liu et al. 2011;Chelouche & Daniel 2012;Grier et al. 2013Grier et al. , 2017Braibant et al. 2017;Hutsemékers et al. 2017;Yong et al. 2017), and to learn the correlation between the continuum and the BLR (e.g., Yee 1980;Kaspi et al. 2005;Punsly et al. 2018) have been published. It is now believed that AGNs are powered by mass accretion onto a super-massive black hole (e.g., Dopita et al. 1998;Marconi et al. 2004;Ricci et al. 2011 and the references therein). Most, if not all, AGNs consist of mainly four important regions: the central engine (super-massive black hole), the continuum region, the BLR, and the narrow-line region (NLR). The study of broad emission lines of AGNs is very important not only to constrain the structure and kinematics of the BLR but also to understand the nature of the central continuum and the engine powering these sources. The broad emission lines we observe in the UV/optical regime are generated by photoionization in the outer regions of an accretion disk (hereinafter AD) that surrounds the central black hole (Kollatschny & Zetzl 2013). Since the BLR is close to the center of the AGN, it responds to continuum variability. It has been recognized that knowledge about the variations in the broad emission lines of quasars provides us with valuable information that can be used to constrain the size of the BLR and the continuum region (e.g., Lewis & Ibata 2004). Furthermore, the correlation between the continuum and emission line properties of AGNs gives the sizes of line-emitting regions and provides us with potentially valuable information to probe the unobservable ionizing continuum region and central engine in AGNs (Cherepashchuk & Lyutyi 1973;Gaskell & Sparke 1986). The technique based on the lag between the intrinsic variability of the continuum and of the broad emission lines, reverberation mapping, has been applied in studying the geometry and kinematics of the BLR (e.g., Peterson 1993Peterson , 2006Gondhalekar et al. 1994;Bon et al. 2018;Lira et al. 2018). Though it has been applied largely following its success, the reverberation mapping technique of determining the size of BLR continues to have several limitations (e.g., Guerras et al. 2013 and the references therein).
Multifractality analysis has been applied to different AGN light curves (e.g., see Bewketu Belete et al. 2018Belete et al. , 2019aBelete et al. , 2019b, and here we study the multifractal, and thus the nonlinear, behavior of the Hβ emission line and 5100 Å continuum flux of the Seyfert type I galaxy NGC 5548. Type I Seyfert galaxies show broad emission lines because their BLR is not hidden by a dusty torus (e.g., see Dopita et al. 1998;Pozo Nuñez et al. 2015, and references therein). Significant variations in the intensities of the broad emission lines are identified as one of the features of Seyfert I galaxies (Robinson & Perez 1990 and the references therein), and Lyutyi (1973) also showed that the optical continuum of NGC 5548 is highly variable. NGC 5548 was among the first Seyfert galaxies to be studied (Bon et al. 2016;Li et al. 2016) and is one of the nearby active galaxies that have deserved more attention (Kollatschny & Zetzl 2013;Mehdipour et al. 2015). In recent years, spectral variability and other properties of this object have been intensively studied by the Space Telescope and Optical Reverberation Mapping campaign (e.g., De Rosa et al. 2015;Edelson et al. 2015;Fausnaugh et al. 2016;Goad et al. 2016;Mathur et al. 2017;Pei et al. 2017;Starkey et al. 2017).
The first reverberation mapping of NGC 5548 (Peterson & Gaskell 1986) has shown that Hβ responded to continuum changes with a delay of only a few weeks. According to Sergeev et al. (2007), inspection of individual broad Hβ profiles over a 30 yr period reveals that the broad emission line profiles can undergo dramatic changes (from a typical singlepeaked profile centered near the systemic redshift of the galaxy to profiles that show prominent blue/red peaks or two peaks). Gilbert & Peterson (2003) have found a nonlinear relationship between the Hβ emission line and the optical continuum flux of NGC 5548. Denney et al. (2009) have seen evidence for outflowing, in-falling, and virialized BLR gas motions based on the optical continuum and broad Hβ emission line variations observed in the nuclear regions of selected Seyfert galaxies including NGC 5548. Using the 43 yr  monitoring data of NGC 5548 (Bon et al. 2016;Li et al. 2016) revealed ∼5700 day periodicity in the continuum light curve, the Hβ light curve, and the radial velocity curve of the red-wing of the Hβ line. Orbiting dusty and dust-free clouds, a binary black hole system, tidal disruption events, and the effect of an orbiting star periodically passing through an AD are explained to be the possible physical mechanisms for the periodicity. Additional studies of NGC 5548 were performed by Bon et al. (2018), Mehdipour et al. (2015), Bentz et al. (2009Bentz et al. ( , 2007, Peterson et al. (1991Peterson et al. ( , 1992Peterson et al. ( , 1999, and Netzer et al. (1990).
Our aim is to characterize the multifractal (nonlinear) behavior of the Hβ emission line and the 5100 Å continuum flux of NGC 5548 and search for any possible nonlinear relationship, or similarity/difference in the degree of nonlinearity between them-how the emission line responds to continuum variations. This provides us with potentially valuable information to model the structure and dynamics of both the continuum and its associated BLRs, which in turn can be used as input to better understand the physics of the central engine and accretion. This work is organized as follows. Section 2 includes the data and analysis method we used. The results are given in Section 3. Section 4 presents a discussion of the results. Finally, the main conclusions are summarized in Section 5.

Light Curves
We use the log-term  light curves of the Hβ emission line and 5100 Å continuum flux of the Seyfert type I galaxy NGC 5548 3 (see Figure 1; Bon et al. 2016;Li et al. 2016).

Framework for Analyzing Light Curves
In the last few decades, stochastic processes have been used to describe the variability of quasars, e.g., computing structure functions (SFs) or simulating light curves. A popular model to account for quasar variability, mainly in the optical band, is the so-called damped random walk (DRW; e.g., Kozłowski et al. 2010;Zu et al. 2011Zu et al. , 2013MacLeod et al. 2012;Ivezić & MacLeod 2014;Kasliwal et al. 2015;Guo et al. 2017). The SF for the DRW model is given by (MacLeod et al. 2010, and reference therein): where τ is the characteristic time, ¥ SF = s 2 , and σ is the variability amplitude. This SF has two asymptotic regimes: The characteristic timescale τ and ¥ SF are the two key parameters, and the behavior of SF∝| | D b t is related to a power spectral distribution PSD∝f −δ , where δ=2β+1. Thus, in Equation (2), one deals with the standard asymptotic limit of a random walk at short values of Δt, i.e., a Wiener process with δ=2, while the damped nature manifests at long values of Δt, when the SF and PSD become constant (δ=0). In addition, the power-law index δ can be also related to the generalized Hurst exponent from a multifractal analysis (see the discussion below Equation (9)). This allows us to compare multifractality-based power-law indexes and those predicted by the DRW model. Naturally, complex systems generate time series that exhibit fluctuations over a range of timescales and/or broad distributions of the values (Kantelhardt 2009). These natural fluctuations are often found to follow a scaling relation over several orders of magnitude, from which one can characterize the time series and the generating complex systems by a set of scaling exponents (fractals). This characterization of the time series and the systems can be used to make comparisons with other systems or/and models. Moreover, determining the nature of scaling between fluctuations in a variable of interest and whether or not some types of power-law exponents are present for various statistical moments at different scales is very significant to better characterize the behavior of the time series and the system generating the time series. Since the dynamics of such complex systems cannot be described by a single power exponent (monofractal), it is necessary to apply a different formalism that can explore the set of power exponents (fractals/multifractals) to fully characterize the systems. Particularly, a multifractal analysis approach is applicable when different scaling exponents are required to describe different parts of the series. There are different techniques that have been developed to characterize the properties of stationary and nonstationary fractal and multifractal time series. Some of them are Hurst's rescaled-range analysis (R/S), fluctuation analysis (FA), autocorrelation function analysis, spectral analysis, wavelet transform module maxima (WTMM) method, detrended fluctuation analysis (DFA), multifractal detrended fluctuation analysis, detrending moving average (DMA), and multifractal detrended moving average (MFDMA) algorithm (Kantelhardt 2009;Gu & Zhou 2010, and the references therein).
The standard FA is a fractal analysis technique for stationary time series and is based on random walk theory (Kantelhardt 2009, and the references therein). For example, for a time series x(t), t=1, 2, 3, ..., N, the global profile is given by Equation (3) and we study how the fluctuations of the profile, in a given segment size s, increase with s. One can determine the square-fluctuations of the profile scale with s, by dividing each record of N elements into N s =int(N/s) nonoverlapping segments of size s starting from the beginning and another N s nonoverlapping segments of size s starting from the end of the considered series. Furthermore, the mean fluctuation function can also be estimated by averaging the square-fluctuations over all subsequences (see Kantelhardt 2009 for the detail). Since FA is unable to analyze nonstationary fractal time series, a different technique called DFA was introduced by Peng et al. (1994). The establishment of DFA was basically to detect longrange (auto) correlations in nonstationary time series. The formalism DFA is also based on random walk theory and basically represents a linear detrending version of FA (Kantelhardt 2009). But the technique DFA presents abrupt jumps in the detrended profile, which is the limitation of the method. To fix this drawback, several modifications and extensions have been introduced. DMA is one of these modifications. Furthermore, DMA was extended to MFDMA to analyze multifractal time series and multifractal surfaces (Gu & Zhou 2010).
Here we apply the backward (θ=0) one-dimensional MFDMA algorithm briefly discussed in Gu & Zhou (2010). In multifractality analysis the most crucial parameters to describe the structural behavior of a time series x(t) are (i) a qth-order fluctuation function F q (n), from which one can calculate the local Hurst exponent h(q), (ii) the multifractal scaling exponent function τ(q), and (iii) the multifractal spectrum function f (α). We obtain these parameters following the procedures explained in Gu & Zhou (2010): 1. Consider a time series x(t), t=1, 2, 3, ..., N. The time series is then reconstructed as a sequence of cumulative sums given by: where N is the length of the data. 2. We calculate the moving average function˜( ) y t of Equation (1) in a moving window using the relatioñ where n is the segment (window) size. The lowest segment size (n min ) is chosen to be 10 and the maximum (n max ) is around 10% of T, where T is the length of the time series. In addition, for our case, θ is adopted as zero, referring to the backward moving average. See Bewketu Belete et al. (2018) and Gu & Zhou (2010) for a detailed explanation. 3. We remove the trend from the reconstructed time series y (t) using the function˜( ) y t and the residual sequence ò(t) as follows: 1 . The residual time series ò(t) is subdivided into N n disjoint segments with the same size n given by In this sense, the residual sequence ò(t) for each segment is denoted by ò v , where ò v (i)=ò(l+1) for 1in and l=(v−1)n.
4. Then we calculate the root-mean-square function using the relation: n v 1 2 1 2 5. We determine the overall fluctuation function F q (n) as: where q is the statistical moment or power argument, which takes a set of zero-mean numbers. For our case it is chosen to be in the interval [−5 5].
When q=0, according to L'Hôspital's rule, we have: We determine the power-law relation between the overall fluctuation function F q (n) and the segment size n by using the relation: where h(q) is the slope (the scaling exponent), to be calculated from the log (F q (n)) versus log(n) plot using the least square fitting technique. This relationship is applicable if the time series has long-range power-law correlations, or if F q (n) increases for larger values of n as a power law. For a monofractal series, a series with a single exponent, the Hurst exponent (Hurst 1951) is sufficient to characterize the behavior of the series at different scales, and thus H=h(q), where H is the generalized Hurst exponent. For a multifractal time series, the time series scales differently at different timescales producing a range of exponents, and thus H is not sufficient to describe the behavior of the series. Therefore, for a multifractal time series ¹ H h(q), and the scaling exponent h(q) can be thought of as "local Hurst exponents" instead (Hampson & Mallen 2011). The local Hurst exponent defines a continuum between a noise-like time series and a random-walk-like time series. For a stationary time series, such as fractional Gaussian noise (fGn), h(q=2) will be between 0 and 1, and h(2)=H (Feder 1988;Sadegh Movahed et al. 2006). For a nonstationary time series (i.e., a time series with timedependent variance) with a random walk or fractional Brownian motion (fBm), the local Hurst exponent h(q) > 1 (Sadegh Movahed et al. 2006;Ihlen 2012), and related to H as H=h(q=2)−1 (Sadegh Movahed et al. 2006;Hu et al. 2009). The H values 0<H<0.5 and 0.5<H<1 indicate antipersistence and persistence long-range correlations, respectively, whereas, H=0.5 indicates an uncorrelated time series. The temporal or spatial fluctuation of complex systems can be characterized by a power-law decaying power spectral density PSD∝f −δ , where δ is the power spectrum scaling exponent. According to random fractal theory, δ can be related to the generalized Hurst exponent H as δ=2H+1 (Sadegh Movahed et al. 2006;Hu et al. 2009). For a fGn process, δ<1, whereas, δ>1 for fBm (Hampson & Mallen 2011, and reference therein). 7. Using the calculated h(q), we determine the Renyi scaling exponent τ(q) as follows: where D f is the fractal dimension of the geometric support of the multifractal measure. For our case, D f =1, because we are applying the MFDMA for one-dimensional time series analysis. 8. We determine the singularity strength function (the Hölder exponent) α(q), and finally estimate the multifractality spectrum function f (α) as follows: The Δα=α max −α min (the range of singularity strength α) is the measure of the range of multifractal singularity strength and denotes the range of exponents present in a time series. If there is strong nonlinearity in the scaling exponent curve, i.e., a different slope for the negative and positive q values, then we can have wider Δα, which reflects strong multifractality behavior in the data considered (Kantelhardt et al. 2002;Telesca et al. 2004). A physical system is said to be intermittent if it has a wide multifractal spectrum or concentrates into smallscale features with a large magnitude of fluctuations enclosed by extended areas of weaker fluctuations (Moffatt 1994;Frisch 1995;Monin & Yaglom 2007). In addition to the width, the symmetry in the shape of α, defined as A=(α max −α 0 )/(α 0 −α min ), where α 0 is the value of α when f (α) assumes its maximum value, provides valuable information about the singularities present in the time series. The asymmetry presents three shapes: asymmetry to the rightskewed (A>1), left-skewed (0<A<1), or symmetric (A = 1). A brief explanation of this can be found in Bewketu Belete et al. (2018) and the references therein.

Results
Here we analyze the multifractal signatures present in the light curves of the Hβ line and 5100 Å continuum of NGC 5548 using the MFDMA procedures. First, we determine the scaling relationship between the overall fluctuation function F q (n) and the scale n for each light curve, Figure 2, to understand how the fluctuation function scales at each scale n. To make sure the nature of the relationship existed between the fluctuation function F q (n) and the scale n, we calculate the slope, usually known as the local Hurst exponents h(q), from the log-log plot of F q (n) versus n, for each moment q using the least square fitting technique, Figure 3. In addition, we determine the classical scaling exponent function τ(q) using Equation (10) and present the results obtained for each light curve in Figure 3. Furthermore, we estimate the multifractal spectrum function f (α) by using Equations (11) and (12) for each light curve, Figure 3. The main purpose of estimating the multifractal spectrum function is mainly to determine the range of exponents present in the light curves and to quantify the degree of multifractality detected, which in turn can be used to identify the light curve with a strong multifractal signature.
Since the data between day 47,508 and day 52,175 are more densely sampled and of better quality than the rest (Bon et al. 2016;Li et al. 2016), we repeat the same analysis for this time segment separately. The results obtained, the fluctuation functions F q (n), the local Hurst exponents h(q), the Renyi scaling exponents τ(q), and the multifractality spectrum functions f (α) are given in Figures 4 and 5, respectively.

Discussion
The results obtained are discussed as follows. The local Hurst exponent h(q) describes the scaling behavior of the scales with large and small fluctuations for positive and negative values of q, respectively (Tanna & Pathak 2014). For a multifractal time series, the fluctuation function behaves differently for negative and positive values of q at the smallest segment sizes (Ihlen 2012). In addition, for a multifractal time series, the local Hurst exponent h(q) shows nonlinear dependency on the moment q and decreases monotonically when q increases (Tanna & Pathak 2014). In contrast, for monofractal time series, h(q) has a constant value, i.e., the local Hurst exponents do not change with moment q (Kantelhardt et al. 2002). As shown in Figures 2 and 3, for both light curves, F q (n) is different for positive and negative values of q at the smallest scales, and the local exponent h(q) has a nonlinear relationship with the moment q and decreases monotonically when q increases, revealing the presence of a multifractal signature both in the Hβ line and 5100 Å continuum light curves of NGC 5548. The difference in the nature of the relationship between the fluctuation function F q (n) and the scale n, or the difference in the behavior of the local Hurst exponents between the light curves indicates the presence of a different degree of multifractality. Judging only from the value of the local exponent at the second moment h(q=2), this is 1.0255 for Hβ and 1.0308 for 5100 Å, so we are dealing with h (q=2)>1 nonstationary time series. Using the relation H=h(q=2) − 1 (Sadegh Movahed et al. 2006), we find that H=0.0255 (Hβ) and H=0.0308 (5100 Å), indicating the presence of antipersistent long-range correlations that tends not to continue in the same direction but to turn back on itself giving a less smooth signal (Hampson & Mallen 2011). In addition, the power spectrum scaling exponent δ=2H+1= 1.0510 (Hβ) and δ=1.0616 (5100 Å), which are consistent with random-walk-like time series or fBm processes (see details in Section 2.2). The measured power-law index of the PSD is a kind of average of the two asymptotic values for the DRW model, and thus, the observed time series may be in rough agreement with DRWs for some choices of τ and δ. However, we note that this last statement is exclusively based on the second moment h(q=2), and one should properly check the DRW predictions for all moments h(q). A multifractal analysis of DRW-based simulated light curves is out of the scope of this paper, although we verified that SFs for both time series over a wide range of time lags have complex behaviors, which cannot be reproduced for DRW models.
To verify the discussion based on the fluctuation function and the local Hurst exponents, we determine the scaling exponent function τ(q). Similarly, for monofractal time series τ(q) is independent of q whereas it shows a nonlinear behavior for a multifractal time series (Kantelhardt et al. 2002;Ihlen 2012;Tanna & Pathak 2014). For the light curves considered here, the nonlinear functionality between the moment q and the scaling exponent function τ(q) confirms the presence of multifractal signature in the light curves. The difference in the behavior of τ(q) between the light curves again shows that the degree of multifractality, and thus nonlinearity between the light curves is different. Once the presence of multiscaling (multifractal) behavior is confirmed, it is important to determine how strong the detected multifractal signature is. To this end, we calculate the width (singularity strength) of the multifractal spectra by using Δα=α max −α min . The width value for the Hβ line and 5100 Å continuum are 0.4507 and 0.2002, respectively. Since width of the multifractal spectrum is a measure of degree of multifractality (Ashkenazy et al. 2003), the wider the width, the stronger (richer) the multifractality (Kantelhardt et al. 2002). Hence, comparing the width values, the degree of multifractality of the Hβ line is stronger than that of 5100 Å.
It has been indicated that the width and shape of the multifractal spectrum are related to temporal variation of the local Hurst exponents. The symmetric spectrum originated from the leveling of the qth-order local Hurst exponent both for negative and positive q values (Bewketu Belete et al. 2018 and the references therein). The leveling of qth-order local Hurst exponents reflects that the q-order fluctuation is insensitive to the magnitude of local fluctuations. A multifractal spectrum will be found with right truncation if the multifractal structure is sensitive to the small-scale fluctuation with large magnitude, whereas the multifractal spectrum will be found with left truncation if the structure is sensitive to the local fluctuation with small magnitudes. For our case, the multifractal spectra of both the light curves are found to be left-side truncated, indicating that the multifractal structure is sensitive to the local fluctuation with small magnitudes and the existence of smallscale intermittency in the continuum and BLR. In addition, we take into account the effects of the observational noise (flux uncertainties) to the width Δα and shape of the spectra. We generated simulated light curves of the continuum and Hβ emissions at epochs equal to those of observation, modifying the observed fluxes by adding random quantities. These additive random numbers were realizations of normal distributions around zero, with standard deviations equal to the measured uncertainties. Our analysis method was then applied to the simulated curves to produce distributions of spectral shapes and widths. Despite the two spectral shapes (Hβ line and 5100 Å continuum) do not change when flux uncertainties are taken into account, we found appreciable differences in the Similarly, for the time segment between day 47,508 and day 52,175, the behavior of all functions, the fluctuation functions, the local Hurst exponents, the Renyi scaling exponents, and the multifractal spectra functions, reveal the presence of multifractal, and thus nonlinear signatures in this part of the light curves. The calculated width values for this time segment of the Hβ emission line and the 5100 Å continuum flux light curves are 0.3390 and 0.1733, respectively. Here also the degree of multifractality of the Hβ emission line is stronger than that of 5100 Å continuum. Similarly, to take the flux uncertainty in this particular period into account, we generated simulated light curves of the continuum and Hβ emissions at epochs equal to observation in this time interval. Though the two spectral shapes remain the same when flux uncertainties are taken into account, we found differences in the spectral width. These differences allowed us to infer 1σ errors: 0.34±0.09 (Hβ) and 0.17±0.05 (5100 Å). This may imply that observational noise could affect the width of the multifractal spectrum, but it may not affect the shape of the spectrum. Bewketu Belete et al. (2018) have shown that the widening in the width of a multifractal spectrum could be due to noise contamination. Now, we discuss possible origins of multifractality detected in the light curves. There are two commonly known origins of multifractality: (i) a broad probability density function for the values of the time series, or/and (ii) different long-range correlations for small and large fluctuations (Kantelhardt et al. 2002). The simplest way to identify the origin of multifractality in the time series is by analyzing the corresponding shuffled fluxes, i.e., in order to destroy the temporal correlations in the original data, observed fluxes are randomly shuffled at the observing epochs. If the observed multifractality is due to different long-range correlations for small and large fluctuations, then there will be no difference in the fluctuation function F q (n) of the shuffled data for positive and negative q at the smallest segment size, and consequently, the local Hurst exponent h(q) will be independent of the moment q, having a constant value h shuff (q)=0.5 at all q. In addition, the scaling exponent τ(q) of the shuffled data will be a linear function of q, and as a result Δα shuff =0. If the multifractality is due to the probability density function, h shuff (q)=h ori (q), and as a result Δα shuff =Δα ori , this is from the fact that the shuffling process does not affect the probability density function. If both sources of multifractality are present, Δα shuff <Δα ori . The fluctuation functions F q (n) of the corresponding shuffled time series are given in Figure 6. In addition, we present the local Hurst exponent h(q), the scaling exponent function τ(q), and the multifractal spectrum function f (α) of the shuffled series in comparison to the corresponding original light curve of the Hβ emission line (Figure 7) and the 5100 Å continuum (Figure 8). For the light curves considered here, h shuff (q) is almost constant, but it differs from 0.5. Additionally, Δα shuff =Δα ori . Thus, although both mechanisms appear to play a role in the observed multifractality, we find that the long-range correlations contribute more to the observed multifractality, because the shuffling procedure greatly reduced the multifractal signature present in the original time series.

Summary and Conclusions
We perform multifractality analysis for the full record and densely sampled period of the Hβ emission line and 5100 Å continuum light curves of the Seyfert type I galaxy NGC 5548 using the backward (θ=0) MFDMA procedures. We also repeat the same analysis for the corresponding shuffled time series to identify the origin of multifractality. First, we calculate the overall fluctuation function F q (n) for each light curve, from which we estimate the local Hurst exponents h(q) using the least square fitting technique. Second, knowing the local Hurst exponents, we determine the Renyi scaling exponent function  τ(q). At last, we estimate the Hölder exponent α(q) and multifractal spectra function f (α), from which we calculate the width of the spectra of the light curves.
We detect multifractality, and thus nonlinear signatures in the full record and densely sampled period of the Hβ line and 5100 Å continuum light curves of NGC 5548. The Hβ light curves (full record and densely sampled period) exhibit stronger multifractality than those of the 5100 Å continuum, i.e., Δα (Hβ)∼2Δα (5100 Å). The presence of different antipersistent long-range temporal correlations for small and large fluctuations is identified to be the main source of the observed multifractal (nonlinear) signatures in both light curves.
The detection of multifractal signature in the light curves may indicate the presence of complex and nonlinear interactions in the continuum and Hβ emission regions, which cannot be described by linear models. For example, a period where the BLR weakly responded to changes in the continuum was observed (Goad et al. 2016). The observed multifractal signature could be due to turbulence in the continuum and emission line regions, because turbulent dynamics produces multifractal and intermittent structures (e.g., see Yordanova et al. 2004;Leonardis et al. 2013). In general, it is the variation in the flux that results in multifractal structure in a time series, and therefore the observed multifractality signature could be due to different physical mechanisms. Regarding the optical continuum, its variability may be generated in several ways: AD intrinsic instabilities or flares in a central irradiating source (CIS) that are reprocessed in the AD (e.g., Collier & Peterson 2001), reprocessing of the CIS-AD (nuclear) variations in the BLR (e.g., Korista & Goad 2001), and others. Korista & Goad (2001) reported that the diffuse continuum emission from the BLR is important for NGC 5548, because the diffuse-to-nuclear continuum ratio is of about 25% at 5100 Å (see Figure 2 in that paper). There is also recent observational evidence for this type of emission in NGC 5548 (e.g., Edelson et al. 2015;Fausnaugh et al. 2016). Thus, signatures of nonlinearity might be produced during reprocessing in the BLR clouds. The nonlinearity of the continuum is less significant than that of the Hβ emission line, which may be plausibly interpreted as a dilution effect, i.e., all the Hβ emission is generated in the BLR after reprocessing of the CIS-AD light, while only a fraction of the continuum comes from BLR clouds.
The presence of long-range temporal correlations shows the existence of periodic and/or long timescale physical phenomena. For example, the observed long memory could be due to the mechanism responsible for the ∼5700 day periodicity in the Hβ emission line and continuum flux light curves of NGC 5548, such as orbiting dusty and dust-free clouds, a binary black hole system, tidal disruption events, or the effect of an orbiting star periodically passing through the AD (Bon et al. 2016;Li et al. 2016). Recently, it has been remarked that AGN variability is a complex phenomenon that requires advanced