Modeling time delays from two reprocessors in active galactic nuclei

Context. Continuum time delays from accretion disks in active galactic nuclei (AGNs) have long been proposed as a tool for measuring distances to monitored sources. However, the method faces serious problems as a number of e ﬀ ects must be taken into account, including the contribution from the broad line region (BLR). Aims. In this paper, we model the expected time delays when both the disk reprocessing of the incident X-ray ﬂux and further reprocessing by the BLR are included, with the aim to see whether the two e ﬀ ects can be disentangled. Methods. We used a simple response function for the accretion disk, without relativistic e ﬀ ects, and we used a parametric description to account for the BLR contribution. We included only the scattering of the disk emission by the BLR inter-cloud medium. We also used artiﬁcial light curves with one-day samplings to check whether the e ﬀ ects are likely to be seen in real data. Results. We show that the e ﬀ ect of the BLR scattering on the predicted time delay is very similar to the e ﬀ ect of the rising height of the X-ray source, without any BLR contribution. This brings additional degeneracy for potential applications in the future, when attempting to recover the parameters of the system from the observed time delays in a speciﬁc object. Both e ﬀ ects, however, modify the slope of the delay-versus-wavelength curve when plotted in log space, which opens a way to obtaining bare disk time delay needed for cosmology. In addition, when the disk irradiation is strong, the modiﬁcation of the predicted delay by the BLR scattering and by X-ray source height become considerably di ﬀ erent. In this regard, simulations of the expected bias are also presented.


Introduction
The standard ΛCDM cosmological model is currently under vigorous discussion and testing via a state-of-the-art approach based on current and upcoming astronomical instruments (Freedman 2017).This global model describes the evolution of the entire Universe, so the measurements of the global model parameters should give the same values independently of how and where they are measured.One such parameter is the current (i.e., at redshift zero) expansion rate of the Universe, the Hubble constant, H 0 .However, the local direct measurements based on SN Ia, calibrated predominantly with Cepheid stars, give average values of H 0 of the order of 74 km s −1 Mpc −1 (e.g., H 0 = 74.03± 1.42 km s −1 Mpc −1 , Riess et al. 2019; H 0 = 73.2± 1.3 km s −1 Mpc −1 , Riess et al. 2021).At the same time, measurements of the cosmic microwave background (CMB) properties imply 67.4 ± 0.5 km s −1 Mpc −1 (Planck Collaboration VI 2020) when the standard ΛCDM is used to derive H 0 .In the observational cosmology this discrepancy in the Hubble constant is known as Hubble Tension.The detection of Hubble tension suggests a need of different cosmological model to explain the local universe.However, before going into detail of deriving a different cosmological model, we have to be very sure that the tension really exists -hence, a range of different independent methods is required to probe it.A similar tension between the Planck results and the local measurements shows up for most (and quite numerous) methods (for the most recent comprehen-CNPq Fellow. sive review, see Di Valentino et al. 2021, in particular Fig. 1 therein).However, the disagreement is not at all that clear since the systematic errors in each method are difficult to assess.
Therefore, a relatively simple and direct method, which would not require the involvement of the distance ladder would be extremely useful to fix the problem.One such method is the one based on continuum time delays in accretion disks in active galactic nuclei (AGNs), proposed by Collier et al. (1999; see also Oknyanskij 1999 for a torus-based version of the idea).The method is effectively based on measuring the size of the accretion disk at different wavelengths and comparing it with the classical accretion disk model of Shakura & Sunyaev (1973).This comparison, due to a specific scaling of both the monochromatic flux and the effective temperature with the product of the black hole mass and accretion rate (see e.g., Panda et al. 2018), allows us to get the distance to the source directly, from the measured time delay τ between the two wavelengths and the observed monochromatic flux, f ν , at one of these wavelengths, without any hidden dependence on the black hole mass and accretion rate.The time delay, τ, as predicted by the theory, depends on the wavelength λ as τ ∝ λ 4/3 , and the proportionality coefficient is also strictly predicted by the theory.This coefficient contains the observed flux and the distance, thus offering the possibility to obtain the redshift-independent distance to the source knowing the time delays and observed fluxes.
Observational monitoring of several sources confirmed the expected delay pattern (e.g., Collier et al. 1999;Cackett et al. 2007;Pozo Nuñez et al. 2019;Lobban et al. 2020), specifically, the proportionality: τ ∝ λ 4/3 .However, the observationally determined proportionality constant was frequently deemed as being much too high (by some 40% up to a factor of a few) in comparison with the theory, that is, when the standard cosmology was used (e.g., Collier et al. 1999;Cackett et al. 2007;Lobban et al. 2020;Guo et al. 2022a) and the narrow filters used by Pozo Nuñez et al. (2019) did not help.Some sources have found a good agreement with expectations, particularly if the height of the irradiating source and/or the extended character of the reprocessor are included (Kammoun et al. 2021a); some authors have found a disagreement only at a single wavelength band, close to the Balmer edge (e.g., Edelson et al. 2015;Kammoun et al. 2019;McHardy et al. 2018;Cackett et al. 2018Cackett et al. , 2020;;Hernández Santisteban et al. 2020).Hints of a considerable problem with regard to the contribution from the more distant reprocessing region, namely, the broad line region (BLR), have appeared (Chelouche et al. 2019).In his most recent paper, Netzer (2022) concluded that (for sources analysed by him) most of the reprocessing actually happens in the BLR region, making a continuum time delay longer than the predicted value from the disk itself.A similar suggestion was made earlier by Lawther et al. (2018) for the case of NGC 5548.If this is true, it is very difficult to disentangle the BLR time delay with the continuum disk time delay -hence, these objects cannot be used for the cosmological purpose.The source of contamination comes, apart from strong emission lines, also from broadband spectral features such as the Fe ii pseudo-continuum and the Balmer continuum (Wills et al. 1985).The problems of reconciling the disk size with the standard model also appeared in microlensing studies (Rauch & Blandford 1991;Mosquera et al. 2013), but the presence of the additional reprocessing medium most likely solves this problem as well.
In the present paper, we model the combined reprocessing by two media: an accretion disk and the extended BLR, with the aim to find a way to disentangle efficiently these two effects.Finally, we aim to use these results to reconstruct the disk time delay.

Method
We performed a set of numerical simulations that allows us to see whether the adopted geometry for the disk and the BLR region can be recovered in measurements of the time delays.We created artificial light curves to mimic the incident radiation, assumed a set of parameters describing the disk and BLR reprocessing, and, finally, calculated the time delays using interpolated crosscorrelation function (ICCF) method to see whether the geometry can be determined and, in particular, the conditions under which the time delay related to the accretion disk alone can be recovered in such simulations.

Incident light curve simulation
According to the general picture used in the description of the disk reprocessing, we assume that variable X-ray emission is responsible for the variability of the disk emission (Rokaki et al. 1993).
We modeled the lightcurve using the algorithm of Timmer & Koenig (1995;hereafter, the TK method), which is based on the adopted shape of the power spectrum.We modeled the X-ray power spectrum as a broken power law, with two breaks and three slopes.The higher frequency break has been relatively well studied for a number of AGNs.For a given set of source parameters, the relation between the black hole mass, the Eddington rate, and the position of the break is given by McHardy et al. (2006).Older and newer measurements are roughly consistent with this law (e.g., Czerny et al. 1999Czerny et al. , 2001;;Markowitz 2010).For the slopes, we assumed −2, −1, and 0, respectively.For the lower frequency break, we assumed that its timescale is by a factor of 100 longer than the short timescale break.This is somewhat arbitrary, since the long timescale trends are not well measured.Alternatively, we could use a bending power law as, for example, in Georgakakis et al. (2021).The normalization of the power spectrum is then adjusted to the required level of the source variability.This model is certainly better than a damped random walk (DRW), corresponding to a single break and slopes −2 and 0, proposed by Kelly et al. (2009) for modeling the optical variability of AGNs.As shown by Yu et al. (2022), more advanced models are needed for precise description of quasar variability in Stripe 82.For our purposes, TK method is satisfactory, as it broadens the frequency range in comparison to DRW.

Accretion disk reprocessing
Our description of the disk reprocessing is relatively simple.For the geometry, we assumed a simple lamppost model that represents the X-ray corona, namely, one that is frequently adopted for the disk reprocessing in compact X-ray binaries and AGNs (e.g., Martocchia & Matt 1996;Miniutti & Fabian 2004;Niedźwiecki et al. 2016), and the disk height is neglected.We do not include general relativity (GR) effects and we assume perfect thermalization of the incident radiation.Once the X-ray photon hits the disk all the radiation get absorbed by the disk which increases the disk temperature locally.Thus, in our model, we do not consider any energy-dependent reflection as in Kammoun et al. (2019) but all the incident emission is absorbed by the disk and gets reprocessed.
The response of the disk to the variable irradiation was recently studied in much more detail by Kammoun et al. (2019Kammoun et al. ( , 2021b)).Their approach included numerous effects, such as the assumption of the Kerr metric in the description of the disks, full GR treatment of the photon propagation, and energy-dependent reflection of the X-rays by the disk surface.We do not aim to achieve such a detailed approach here.Instead, we adopt the rather simple approach of Shakura & Sunyaev (1973) for the disk and geometrical optics, where a perfect thermalization of the incident X-ray flux is considered.In creating our software, we aim to check qualitatively the effect of the second reprocessor on the measured time delay and we do it not only with the use of the transfer function approach, but also with the use of artificial light curves simulating the data cadence and quality.In such cases, we simply carry out a direct calculation of a sequence of disk reactions to variable irradiation by X-ray flux, where the incident flux can vary with arbitrarily large amplitudes.
Our numerical scheme is thus different from the one used by Kammoun et al. (2019Kammoun et al. ( , 2021b)).We first set the radial grid covering the accretion disk between R isco and R out by a variable sample size.The radial bin size is quasi-logarithmic, to allow for a proper resolution at the inner parts.Specifically, we define this grid by the following formula, dR = 0.085 * ( R R isco ) 0.85 , and for each "R," we also increase the grid step in angular (φ) direction by dφ = 1.5700N div .For a given R and φ, we calculate the Cartesian (x, y) coordinate in the disk plane, and the surface element, ds = R * dR * dφ * R g * R g .The disk height is neglected, that is, we assume z = 0.For a given x, y coordinate, we calculate the total delay τ total (x, y).This delay is the sum of time τ d (r) taken by photon to reach the given disk location from corona, located at the height, H, along the symmetry axis to accretion disk, and the A147, page 2 of 14 time to reach observer after disk reprocessing, τ do (x, y).For convenience, we define this last time delay with respect to the plane crossing the equatorial plane at R out , φ = 0, and perpendicular to the direction towards the observer.This delay depends on the inclination angle i, the corona height, the position on the accretion disk, and the black hole mass.All the delays are calculated with respect to the photon generated from the corona.
Once we have the total delay τ(r) for a differential disk elements, we calculate the new temperature from the total flux, F, the sum of flux generated from the non-irradiated disk, and irradiation.We thus used the following expression: The local temperature at each moment is calculated from the local flux assuming a blackbody emission: where σ B is the Boltzmann constant.We do not include the color correction, which would affect the normalization of the time delay but does not affect the basic trends with the parameters that are the topic of our study.However, the inclusion of the color correction indeed makes all the delays longer, which is critical for actual data fitting.From Planck's formula, we can generate the entire spectrum for the given differential area only to store it in a photon table, which is a 2D matrix P(t , λ), using a determined time delay appropriate for the disk position and wavelength.Every flux element of the photon table has a unique delay and a wavelength λ.We chose the values of λ from 1000-10 000 Å for the simulation, using a logarithmic scale grid for selection.The time delay step is either linear or logarithmic, depending on the need.
We assume that the incident light curve is provided with equal bin size, ∆t.We treat such a light curve as a histogram, which is a step function with ∆t.Thus, the time bin size of the photon table matches the resolution of the incident light curve and the time span must be long enough to cover the duration of the incident lightcurve, t irrad .For each incident light curve bin, we include this uncertainty of ∆t of the photons arrival when locating the flux contribution into the photon table.We first performed the loop with respect to incident light curve bins for a given location at the disk, and then repeated the process for all the elements of the disk surface, thus creating the lightcurve expected from the irradiated disk, equivalent to: for a very dense grid.We can calculate the response function of the disk, ψ d , which is a very useful concept if the response of the medium is linear (see e.g., Eq. (9) of Peterson 1993).In this case, we replace the light curve L X (t) with a very short impulse of duration of one second, and we normalize the result based on the incident bolometric luminosity: We do not subtract the flux from the non-irradiated disk, since we do not aim to linearize the equation and, in general, ψ(t, λ) depends on the parameters, including L x , which is a function of time.When generating the response function, we used two types of time bins: linear and logarithmic -as the linear time bin failed to capture the smallest delays created by the disk elements close to the black hole.
To create the disk light curves from the incident X-ray light curves, we generally do not use the concept of response function as defined by Eq. ( 4), but we calculate the result directly from Eq. ( 3).In this approach, the amplitude of the irradiating flux can be arbitrarily large in comparison with the locally dissipated radiation flux.

Second reprocessor
The BLR is a secondary medium responsible for the reprocessing of the irradiating flux (Lawther et al. 2018;Korista & Goad 2019;Chelouche et al. 2019;Netzer 2022).This medium is the source of the emission lines and the emission line time delays have been measured by many authors since many years using the reverberation mapping technique (e.g., Blandford & McKee 1982;Peterson 1988Peterson , 1993;;Kaspi et al. 2000;Peterson et al. 2004;Bentz et al. 2013;Du et al. 2014;Grier et al. 2017;Martínez-Aldama et al. 2019;Du & Wang 2019;Panda et al. 2019;Zajaček et al. 2021).This process has mostly been modeled and measured as an independent process, neglecting the disk reprocessing.However, as pointed recently by many authors (Korista & Goad 2001, 2019;Netzer 2022;Guo et al. 2022b), BLR is also the source of diffuse continuum, including Balmer continuum and scattering, which also vary -this, together with some level of line contamination affects the measured continuum time delays.
In this paper, we focus on the simplest aspect of the BLR, which is electron scattering of the photons by the inter-cloud BLR medium.Such a process does not imprint any characteristic features as a function of wavelength but can effectively modify the predicted net time delay between the two continuum bands.The schematic illustration of the geometry is shown in Fig. 1.We include this effect through the Thomson scattering approximation.This means that the scattering does not change the photon frequency.
In order to model the scattering effect, we used a simple analytical parametrization of the response of the second reprocessor, ψ BLR .There are no simple direct observational determinations of such a response, but we can look for suitable parameters by looking at the measured responses for Hβ lines (e.g., Grier et al. 2013;Xiao et al. 2018;Du et al. 2018) and assume that the inter-cloud medium follows a similar distribution.Thus, for the shape of this function, we assumed one symmetric function in the form of a Gaussian and two asymmetric functions: half-Gaussian and an exponential decay.The peak location of the function in the time axis is given by T peak .The non-zero contribution is included starting at a specific time delay t 0 up to t max .The width of the Gaussian, σ, or the decay timescale of the exponential (t decay ) are the parameters of the model.Exemplary shapes are presented in Fig. 2. We assume the same shape of the ψ BLR for all the disk photons, independently from their wavelength and location on the accretion disk.
In the case of the presence of the second reprocessor, we introduce a parameter f BLR which weights the relative contribution from the two regions.The parameter f BLR accounts for the fraction of photons scattered by the inter-cloud medium and introduces an additional time delay.The factor (1 − f BLR ) represents the fraction of disk photons that reaches the observer without any scattering or additional time delay.In case of f BLR different from zero the total time delay is a combination of two processes: initial delay due to light travel time from the X-ray source to the disk and the extra time delay due to scattering in extended BLR inter-cloud medium.
The assumption of extended BLR medium modeled ψ BLR implies that the final result is a convolution of these two effects: The value of f BLR can vary from 0 (no BLR reprocessing) to 1 (no disk reprocessing).The usually expected value is rather in the range from 0.1 to 0.3, if measured by an estimated solid angle of the BLR.However, fully ionized medium can have a different spatial distribution than the BLR denser clouds.
If we aim to calculate the effect of reprocessing from the long X-ray lightcurve, we apply the response function directly to our photon table P(t , λ), using the same parameter f BLR : a fraction of (1− f BLR ) remains unchanged and the fraction of f BLR is smeared by ψ BLR (t).

Time delay measurements
We considered two mathematical approaches to determine the expected time delay.In the first method, we use a single flare event (not a delta function but of final duration of 1 s), we determined the response function by the disk for the corresponding parameters and we eventually added the response function of the second reprocessor.In this case (see Sect. 2.4.1), the whole light curve is not created, so this method is most accurate but does not adequately represent how the time delay is measured in the actual observational data.In the second approach, we created a realistically sampled light curve (incident and in each of the bands), and the time delay is measured using methods comparing the two light curves (see Sect. 2.4.2).Independently from the method, all delays are always measured with respect to the X-ray flare event.

Single-event delay
In the case of a single event, we constructed the standard response function ψ for the accretion disk with the lamp post geometry (see Sect. 2.2), combined it with properly normalized response function from the second reprocessor, and calculated the expected time delay using the formula (Koratkar & Gaskell 1991): Computations in this case require much denser time grids in comparison with the computations for long X-ray light curves (see Sect. 2.2), since the onset of the reprocessed flare has to be well resolved in this case and that happens very close to the black hole.

Realistic light curves
In this case, we simulated the entire incident radiation curve, with realistic sampling, and determined the observed continuum curves at selected wavelengths.Three methods are most frequently adopted: interpolated cross-correlation function, ICCF (Peterson 1993;Sun et al. 2018), JAVELIN (Zu et al. 2011(Zu et al. , 2013(Zu et al. , 2016)), and χ 2 (Czerny et al. 2013;Bao et al. 2022).In the present paper, we concentrate on the first one (ICCF), which brings rather stable results.

Mkn 110 as a motivation for the adopted parameters in simulations
Complex variability was recently discussed in detail for the source Mkn 110 (Vincentelli et al. 2021).So, in order to put our simulations onto realistic footing, we predominantly focus on parameters well representing this source and the light curve duration and spacing characteristic for this source.Mkn 110 (Markaryan 1969)  that the cold disk down to 1.24r g is needed (Mantovani et al. 2016), although we note that such a fit was achieved for unlikely inclination of 80 deg.A more recent study of combined data from XMM-Newton and NuSTAR gave more conservative values of ∼20r g (Porquet et al. 2021).The level of polarization in the optical band in this source is low (∼0.5%,Afanasiev et al. 2019), but the low polarization does not necessarily imply the low optical depth of the scattering medium (e.g., Śniegowska et al. 2022) and the presence of very highly ionized medium is revealed through the presence of Fe xxvi emission line (Mantovani et al. 2016).The viewing angle is estimated at 18.0 ± 3.1 deg.(Afanasiev et al. 2019).The black hole mass in this source was estimated to be 2×10 7 M (Bentz & Katz 2015) and also adopted by Vincentelli et al. (2021).Older measurements have claimed a higher value, 1.2 × 10 8 M , and even larger values has been determined from the recent polarization method, log M = 8.32 ± 0.21 M (Afanasiev et al. 2019).The source has been monitored with Neil Gehrels Swift Observatory (Gehrels et al. 2004), with a cadence roughly once per day, for about 200 days, and we adopted this setup in our simulations.A detailed disk reverberation is studied by Vincentelli et al. (2021), using good cadence data from Neil Gehrels Swift Observatory and LCO network and they found the variability on two different time scales in Mkn 110.The variability time scale below ten days is mostly consistent with accretion disk reverberation with a maximum two-day lag between the shortest wavelength (W2 band) and longest z-band.On the longer time scale they found that the g-band lags the hard X-ray (from BAT) by ten days and a similar lag was also observed between the z-and g-band which is not consistent with the disk reverberation.The author proposed that the longer time scale and higher time delay can be due to the contamination from the diffuse emission of the BLR.Vincentelli et al. (2022) also discussed the effect of X-ray luminosity on the lag spectra.During the low X-ray luminosity state, they do not observe any u-band excess and negative time lag excess as frequently seen in many AGNs.However, at high a X-ray luminosity state, the u-band excess is visible, which is dedicated to the diffuse BLR emission.However, these authors also argue that the excess lag in X-ray can also be explained by moving the corona to farther distances.In the same sense, our study also sees the possibility to disentangle the BLR and corona height contribution in order to explain the disk reverberation.

Results
We presents the result of our project aimed at testing the time delay in the presence of the two reprocessors using simulated light curves.This allows us to test under which circumstances, if at all, the time delays from the disk alone can be recovered.

Response function of the disk and single event delay
We first calculated one example of the response function, that is, the result of the reprocessing of the delta signal from the X-ray lamp by the disk, without any presence from the second reprocessor.Such computations require much denser grids in space and time as specified in Sect.2.2 to adequately see the onset of the radiation.For this exercise, we used the accretion disk model with the following parameters: M = 10 8 M , L/L Edd = 1, lamp height of H = 5R g , the lamp luminosity of L X = 10 40 erg s −1 , and the viewing angle of i = 30 deg.All the delays are calculated with respect to the corona.Although we can generate a response function for any wavelength, we usually store and show     the response function only for nine values of λ, from 1258.92 Å (response1) to 7961.59 Å (response9), adopting a constant logarithmic step.We show the results in Fig. 3.The overall shape is similar to the responses derived by Kammoun et al. (2021b), although we do not have GR corrections.We comment more quantitatively on this issue in Appendix A.
We calculated the centroid times of all response functions using Eq. ( 6).The results (see Fig. 4) are consistent with the simple analytic formula of Collier et al. (1999).We also compared the normalization of the best fit τ ∝ λ 4/3 trend with the expectations from Collier et al. (1999).Their formula contains an unspecified factor X which accounts for the peak contribution to the total emission from a given radius through a scaling of X = hc/(kT λ), which they estimated to be of the order of 3-4.Our value, derived from numerical computations, gives X = 2.47, much lower than the factor suggested by Collier et al. (1999), but higher than the semi-analytical relation (X = 1.65) proposed by Siemiginowska & Czerny (1989).

Response function from two reprocessors and single event delay
We go on to calculate the time delay from Eq. ( 6) as a function of the wavelength, for a combined disk and BLR effect and for several values of the parameter f BLR .
A147, page 5 of 14   The exemplary response functions for the two reprocessor setup is shown in Fig. 5.We plot the shape for the shortest wavelength only, but for three shapes of the BLR response illustrated in Fig. 2.
We see that in the two-reprocessor setup, the combined response has two peaks and the shape of the new response depend on the adopted description of the BLR.In Fig. 6, we plot the response functions for representative wavelengths, selecting half-Gaussian (Fig. 2, middle panel) that represents the BLR.We see that the deep valley between the two peaks becomes more shallow when we go towards longer wavelengths and, finally, the two-peak structure disappears.We note that the dependence of the time delay on wavelengths at the longest time delay is not due to wavelength-dependent effect in BLR itself (as we assume Thomson scattering in the inter-cloud medium), but it is due to the fact that those are photons generated at large disk radii, with the delay generated between the X-ray source and their origin and the net effect is a convolution, as given by Eq. ( 5).
Next, for the same parameters, we calculated the time delays from two reprocessors.The results are shown in Fig. 7.For f BLR = 0 the usual τ ∝ λ 4/3 is recovered.Now we add an offset in the Eq. ( 5), which corresponds to the fraction of BLR.The time delay for different contribution of BLR is estimated and in the linear scale the delays at higher BLR contribution are just shifted with respect to f BLR = 0 (lower panel of Fig. 7).However, the result appears to be very interesting when we plot the time delays in log-log space (upper panel of Fig. 7).With increasing f BLR , we obtain more shallower relations than the standard one.So, not only does the time delay become longer overall due to the extra scattering in the BLR region, but also the slope of the relation changes.In the extreme case, when the disk contribution becomes small and most of the photons are actually scattered, the delay only weakly depends on the wavelength.However, some dependence on the wavelength would remain even for 100% BLR contribution since the photons (before going on to the observer) are first reprocessed in the disk with a wavelengthdependent delay.It is important to note that the change of the slope is actually seen only if we use log-log plot.It is also interesting to show the dependence of the time delay at selected wavelengths, but as a function of BLR contribution.Netzer ( 2022) argued (partially following Lawther et al. 2018) that the delay from two reprocessors is a linear combination of the two time delays (disk and BLR) weighted with the flux contribution.We plot the expectations from our model in Fig. 8.The plot supports the claim of the linear dependence of the time delay on the BLR contribution for all wavelengths if this contribution comes from scattering.The dependence is indeed linear, with the time delay plot shifted up with the increase in the wavelength.
The contribution from the scattering by the inter-cloud medium in BLR can thus easily account for too large disk sizes claimed from the data.In addition, the change of the delay shape shown in Fig. 7 might, in principle, reveal this effect in the data.However, Kammoun et al. (2021a) were able to fit well the observational data for a number of sources at the expense of postulating large height of the illuminating source.Thus, our more general model -which includes the arbitrary lamppost height and arbitrary contamination by the disk photon scattering -might be degenerate with respect to these two parameters and, in the data fitting in the future, we will not be able to discriminate among them.In order to address this problem in advance, in our simulations, we calculated the effect of these two parameters for a range of lamppost source luminosities.
First, we assumed a very faint lamppost of 10 40 erg s −1 , and we repeated the disk delay computations for several height values and compared the results with the expected time delay for small height but with BLR contribution.We see that a change of the lamppost height leads to a flattening of the delay curve, as shown by Kammoun et al. (2021b).In particular, Fig. 18 of their paper shows that the effect is very similar to the introduction of the BLR scattering.
The curvature introduced to the delay curves in lowluminosity case is actually very similar in the case of an increased height or some BLR scattering (see Fig. 9, upper panel).It is a serious source of the degeneracy in the future data fitting, although the relative importance of the two effects depends strongly on the parameters.For example, from our standard model (σ = 10 days, L X = 10 40 erg s −1 , the contribution of BLR just ∼7% gives the same effect as moving the lamp height from 5 to 100R g , and smaller values of σ further reduce this factor. Next, we repeated the simulations for the incident X-ray flux to 30% of the disk bolometric luminosity.Such high luminosity (L X ∼ 3 × 10 46 erg −1 ) brings different result: 9% BLR contribution delay is only matched at the smallest wavelength, not at longer wavelengths.Thus, the overall curvature becomes different (see Fig. 9, lower panel), which opens up the possibility to differentiate the effect of the height from the effect of the BLR scattering in the data.This shows that for X-ray-bright sources, we can differentiate between the lamppost height and the BLR contamination if the data is of sufficient quality, but it might be much more difficult for X-ray weak sources.

Time delays from light curves
In the case of observational data, we do not have a direct insight into the response function; however, some techniques allow us to recover it from the data.The observed light curves depend not only on the system parameters, but also on the sampling, while the measured time delay depends not only on the light curves but on the method to determine the time delay.Therefore, we repeated our experiment from Sects.3.1 and 3.2, using the long artificial light curves with different adopted power and data sampling.As a standard, we adopted the frequency break timescale of 75 days, sampling of one day, and the ICCF method.Computations are based on ten statistical realizations of the process, allowing us to show the errors representing the dispersion, namely, a likely error in a single measurement.

Role of the frequency break in the power spectrum in disk time delay measurement
First, again for test purposes, we calculated the examples of X-ray light curves representing different intrinsic timescales assumed in parametrization of the power spectrum.We adopted the time step of one day, duration of 200 days, and the irradiating X-ray lightcurve was calculated as described in Sect.2.1.We assumed the level of X-ray variability of set by normalized dispersion of 0.3 in the whole light curve of duration of 10 8 s.
The examples of the light curve for three values of the high frequency break are shown in Fig. 10.We see that a small value for the timescale corresponding to the frequency break gives much sharper values of the curve peaks and much higher variability amplitude in a period of 200 days (i.e., much shorter than the whole curve duration).All three curves were obtained from the same value of the parameter initializing random generator -for a better comparison.The X-ray curve seems more smooth when the timescale corresponding to high frequency break is longer, but otherwise, the geometry of the system is not affected.In order to check whether this indeed could affect the measured time delay, we calculated the delay for the three values of the frequency breaks, using ICCF.The results are shown in Fig. 11.
Comparing the time delays obtained from the light curves to the delays calculated directly from the response function (Fig. 4, upper panel), we see that for a black hole mass, 10 8 M , the time delays are relatively well recovered only in the case of a 75-day characteristic time variability (middle panel of Fig. 11).When variability is faster, the numerically calculated time delays are systematically much too short in comparison with expectations.If the variability timescales are longer, the numerical delays are marginally consistent with expectations within an error, but they again locate themselves systematically below the expected values.At the shortest wavelengths, even the optimum variability timescale underestimates the delay, but this is directly caused by Fig. 11.Time delay as a function of the wavelength calculated with ICCF method from the artificial light curves for the frequency break corresponding to timescales of 7.5 day (upper panel), 75 days (second panel), and 750 days (third panel).Model parameters: black hole mass 10 8 M , Eddington ratio = 1.0, corona height = 5r g , inclination angle = 30 degree, and luminosity is on the order of 10 45 erg s −1 .Black points and continuous line represents the delays expected from the disk response function.the adopted one-day sampling which is not enough to resolve the innermost part of the disk.Increasing the incident X-ray flux to ∼3 × 10 46 erg s −1 does not improve the results.We discuss the issue later.
This trend to obtain numerically the time delays which are shorter than expected is rather interesting and potentially important for actual data analysis.The optimum characteristic variability of 75 days is within the duration of the total lightcurve of ∼150 days, so characteristic peaks are a few and apparently well sampled.If the actual timescale is much longer than the total observing time, we may have no strong features to rely on for A147, page 8 of 14 10 4 2 × 10 3 3 × 10 3 4 × 10 3 6 × 10 3 [Å] 10 1 10 0 Delay[days] binsize=1.0day binsize=0.5day binsize=0.25day response function delay Fig. 12.Time delay as a function of the wavelength calculated with ICCF method from the artificial light curves for different bin size.Black curve is expected delay from the response function.Model parameters: black hole mass 10 8 M , Eddington ratio = 1.0, corona height = 5r g , inclination angle = 30 degrees, and luminosity is on the order of 10 45 erg s −1 .time delay measurement.Indeed, if we use the light curves of the duration of 1000 days (again, with a one-day sampling) and the remaining parameters unchanged, the agreement between the numerical results and predictions is much better.Additionally, the asymmetry in the time reprocessing by the disk can also contribute if there are only very few strong peaks in both curves.On the other hand, if the characteristic timescale is much shorter than the total observing run and the sampling rate is not very dense the curve is too noisy.It might thus be recommended to check the characteristic timescale in the data (e.g., using the structure function) and compare it to the derived time delay in order to additionally discuss the potential bias in the time lag determination.

Bin size effect in the disk time delay measurements
As we show in Fig. 11, the prediction of the delays from the simulated light curves systematically underestimate the delay by ∼30% at 4000 Å, and the effect is stronger at the shortest wavelengths.
In order to see whether this is the result of inadequate sampling, we repeated the analysis for just one frequency break corresponding to 75 days, but for an increased data sampling and keeping the total length of the curves unchanged -effectively increasing the number of observational points.The effect is shown in Fig. 12.Indeed, with the denser sampling the numerical light curve time delay was systematically approaching the expected response of the disk.Already, the sampling of 0.5 day was enough to measure well the time delay at 2000 Å and longer wavelengths, for the adopted black hole mass of 10 8 M .The shorter wavelengths ∼1000 Å would require still denser sampling, as even 0.25 of a day would still underestimate the delay almost by a factor of 2.

Time sampling and the black hole mass
The sampling rate of the lightcurve must be adjusted to the source parameters.Our previous discussion focused on 10 8 M black hole mass.However, if we increase the black hole mass by a factor of 10, the delays are still not well recovered at the 10 4 2 × 10 3 3 × 10 3 4 × 10 3 6 × 10 3 [Å] 10 0 10 1 Delay[days] Fig. 13.Time delay from a bare disk as a function of the wavelength calculated with ICCF method from the artificial light curves.Model parameters: black hole mass 10 9 M , Eddington ratio = 1.0, corona height = 5r g , inclination angle = 30 degrees, luminosity is of the order of 10 45 erg s −1 , and the frequency break corresponding to 75 days.
shortest wavelengths as shown in Fig. 13.However, at longer wavelengths, the delay is comparable to the value expectations based on the response function and the slope is well recovered (see Fig. 4).Thus, for a larger mass, a one-day sampling is fully adequate.

The influence of stochastic approach on time delay with BLR contribution
In the case of the analytical (response function) approach, the presence of the additional scattering in the BLR resulted in a simple shift in the net time delay, as expected previously (e.g., Netzer 2022;Lawther et al. 2018).However, the stochastic light curve approach may not preserve such a simple trend in all parameter range.We selected the timescale break of 75 days, since it was working relatively well for the disk delay, and we adopted a one-day sampling which was adequate at longer wavelengths.The resulting light curves thus become distorted (smoothed and shifted) by the BLR (as illustrated in Fig. 14).
Smoothing is clearly stronger when the width of the BLR response function is larger.We now calculate the time delay using those stochastic curves.We noticed that the numerically calculated time delay does not increase with f BLR up to the critical moment when f BLR crosses (unrealistic) value of 50%.This is in contrast to expectations based on the response function approach.The peak in the combined response function is still due to the disk for smaller values of f BLR , which is apparently confusing in terms of the numerical method.We see the same effect for the other values of the width of the Gaussian.It may indicate that in real data analysis, we would actually be recovering the disk delay, independently of the BLR contamination.
However, the symmetric Gaussian shape for the BLR response is unlikely, so we repeated the same analysis for half-Gaussian shape.We found that for half-Gaussian shape the BLR contamination shows more similarity between the time delay predicted by response function and the stochastic prediction, although the determined lags are always below the ones expected from the combined response function.
We finally checked, in a systematic way, how the adopted width of the BLR response function affects the delay.This time  we performed ten simulations for each parameter set and the errors mark the dispersion.We illustrate the complex trend in the time delay with the change of the response model for BLR, σ, and f BLR in Fig. 15.We see that for a Gaussian shape the departure from the linear trend of the rise of the expected delay with the importance of the BLR contamination is strong.But most of the other shapes also predicted similar trend -the initial rise was slower than expected and only after crossing rather unrealistic level of BLR contribution (above 50%), the time delay flipped to values close to the BLR time delay.For example, when the departure between the measured time delay and the linear time delay is determined at 50% of the BLR contribution we see a delay longer by 33.88% (Gaussian), 28.85% (half-Gaussian), 15.01%(half-Gaussian2), 35.03% (half-Gaussian3), 25.62% (half-Gaussian4), and 0.66% (half-Gaussian5).We refer to the caption of Fig. 15 for the model parameters.Thus, no departure is seen for very wide asymmetric BLR response profile while narrow asymmetric or symmetric response show a considerable departure from a linear trend.
We see from the performed simulations that the measured time delay depends on the light curve properties as well as BLR response, and the dispersion in a single measurement is considerable.This means that when an actual reverberation mapping campaign is performed, corresponding to a single realization of our process, some modeling adjusted to the observational setup and source properties is useful for estimating the possibility of the systematic bias in the measured time delays.

Discussion
We studied the wavelength-dependent time delay of optical photons originating from the X-ray photons generated in the lamppost geometry above the AGN accretion disk and reprocessed by the surroundings.We included the photon thermalization and re-emission in the accretion disk, but we also allowed for an additional scattering of the generated optical photons by the intercloud medium of the BLR.Such a scattering does not change the photon energy but introduces an additional time delay with respect to the arrival of the primary X-ray emission as well as with respect to the unscattered optical photons.We constructed the response functions for the combined effect of the accretion disk and stud-ied the time delays analytically, but we also constructed simulated X-ray light curves and their reprocessing.
The results based on the response function computations give a very smooth dependence on the model parameters.The most interesting result of this study is the modification of the time delay by the rising contribution of BLR scattering.This effect is difficult to distinguish from the effect of rising the height of the lamppost, without postulating any contribution from the BLR.In noticing the difference in the curvature of the time delay pattern is practically impossible even with high-quality data, if the mean incident X-ray flux is small in comparison with the disk bolometric luminosity.On one hand, this degeneracy between the lamppost height and the BLR contribution can account for surprisingly large heights obtained from data fitting.Kammoun et al. (2021a) successfully modeled the time delay in seven nearby AGNs, but the derived height of the lamppost ranged from 11.2R g (for Mkn 509, maximally rotating black hole) to ∼75R g for NGC 7469, independently from the spin.This is not consistent with many of the fits of the X-ray spectra that require low lamppost heights to model the relativistically broadened Kα line (e.g., Parker et al. 2014;Jiang et al. 2019;Walton et al. 2021).We have an independent insight into the geometry of the X-ray reprocessing from the measurement of the Kα line delays, and they rise with the black hole mass from ∼100 to 1000 s for mass increasing from 10 6 M to 10 8 M (Kara et al. 2016); for 10 8 M , this implies a geometrical delay of ∼2R g .However, in the case of NGC 7469, the Kα line is broad (broadening velocity about 2700 km s −1 ) but not relativistically distorted (Mehdipour et al. 2015), so it can come from the outer disk and/or BLR, so a large height is not in contradiction with the X-ray spectrum.
High values of the irradiating flux allow us to differentiate the delay curve shape caused by the increase of the lamppost height and by the BLR scattering.The question of whether such high values -namely, up to 30% of the disk bolometric luminosity -are possible is directly related to the question of the origin of the irradiating flux.Hard X-ray emission, as argued by Kubota & Done (2018), contributes less than 2% to the bolometric luminosity of bright AGNs, not containing inner ADAF.On the other hand, soft X-ray excess can contain much higher fraction of the total flux.The lamppost model is more likely to represent better the hard X-ray emission while the geometry of the soft X-ray emission is still under debate, but it is most likely a warm corona (e.g., Czerny et al. 2003;Różańska et al. 2015;Petrucci et al. 2020).However, studies of other geometries besides the lamppost is beyond the scope of the present paper.
In general, a potential data fitting of the time delay faces a number of degeneracies.As demonstrated by Kammoun et al. (2021a), independent information about the black hole mass and accretion rate would reduce it considerably; usually, estimates of the black hole mass, based on line widths, are available.Knowing the monochromatic flux, we can also estimate the accretion rate in a way that only weakly depends on the black hole spin.The viewing angle remains, however, an issue, since the monochromatic flux roughly depends on cos(i).However, a dusty or molecular torus limits the available viewing angles to between 0 and ∼70 (e.g., Prince et al. 2022, and the references therein).The new degeneracy between the lamppost height, H X , and the BLR contribution, f BLR , creates an additional issue.When the height is small and the high-quality X-ray data are available, we can independently estimate its height, but no such estimate is possible if the height is large and the relativistic distortion of the line is not strong.Notes.The corresponding plot is shown in Fig. 15; τ is the delay measured in days between X-ray light curve and 7961.59Å light curve, and ∆τ is error in delay measured in days.
Perhaps, in the future, a more careful modeling of the disk reprocessing plus BLR scattering may help to ease the problem.In our simple code, indeed the effect of the height and the effect of BLR show similar trend for low X-ray luminosity; whereas, in Kammoun et al. (2021b), the disk height results in a convex-shape plot of the time delay versus wavelength, while in our simple model the pattern is concave both for height and BLR contribution.We think that the shape should actually be concave, and the convex shape results from too small outer radius adopted in the computations (see Appendix A).Repeating the calculations of the disk plus BLR scattering using full GR, X-ray reflection, and color correction to the temperature may reveal a systematic difference in the system reaction to these two parameters.In this case, the data fitting should not be done just using a power law part of the delay curve, as in Kammoun et al. (2021a), but the fits should include the full wavelength-dependent model with the curvature.Also, studies of the same source at different flux levels are very helpful in disentangling the lamp height and BLR effect, as argued by Vincentelli et al. (2022).
In the present study, we did not include the re-emission by the BLR clouds.Such emission has clear spectral signatures, A147, page 11 of 14 including the prominent Balmer edge (Lawther et al. 2018;Korista & Goad 2019;Chelouche et al. 2019;Netzer 2022).This effect is also important but, in principle, it is easier to include it later (in the data fitting), since the prominent Balmer edge should fit the corresponding drop in the time delay.In numerical computations of the reprocessed BLR component with the use of Cloudy (Ferland et al. 2017) or equivalent code, the effect of scattering is included but only for (usually) constant density clouds, not accounting for the inter-cloud medium.Thus, the scattering effect can be more difficult to disentangle in the real data.Both broad band data, possibly dense in the wavelength (e.g., coming from specially designed narrow-band filter photometry, as used by Pozo Nuñez et al. 2019), but also very dense cadence is essential, as we can see from our experiments with the artificial light curves.Also, the broad wavelength coverage is very important since it allows to determine the shape of the relation more accurately and to improve the disentangling of the contributions from the disk and BLR.Finally, there are two other possible effects that could modify the delay obtained for the disk continuum: the disappearance of the inner cold disk and the disk winds, as argued by Zdziarski et al. (2022) -an insight into this issue could be expected from a fitting of the broadband spectra of the studied objects.

Conclusions
The results of our modeling of the X-ray reprocessing by the accretion disk, with the additional scattering of disk photons in the BLR region are as follows: -For low-irradiating X-ray flux, the lamppost height and BLR contribution through scattering are degenerate; -For high-irradiating flux, there is a difference in the curvature in delay versus wavelength plot that allows us to distinguish between the two effects -if the wavelength coverage is broad enough; -The time delay rises linearly with the BLR contribution in the description, which uses the response function; -When stochastic incident light curves are used, the time delay is aptly recovered only if the time-step of the curve is considerably denser than the characteristic variability timescale (set by the high-frequency break in the power spectrum) and when the total duration of the light curve is much longer than this timescale; -In numerical stochastic incident light curves, this linear dependence is perturbed and the time delay rise is initially slower than linear, then rising rapidly with the BLR contribution; -Our modeling shows that the results of the time delay based on a single observational campaign should be supplemented with simulations in order to identify the potential bias in measuring the time delays.

Fig. 1 .
Fig. 1.Geometry of the reprocessing by the extended disk and extended BLR.We include only the scattering from the inter-cloud medium.

Fig. 4 .
Fig. 4. Time delays calculated from the disk response functions from Fig. 3. Parameters as in Fig. 3.

Fig. 8 .
Fig. 8. Time delay as a function of the BLR contribution for exemplary wavelengths.For BLR profile, we used half-Gaussian as shown in the middle panel Fig. 2. Model parameters are as in Fig. 7.

Fig. 9 .
Fig. 9. Time delay as a function of the wavelength for models with different lamppost height and for no BLR and for BLR contribution.For the BLR profile, we used a half-Gaussian, as shown in the middle panel of Fig. 2. Model parameters (top): black hole mass 10 8 M , L/L Edd = 1, inclination angle 30 degree, and luminosity of the corona is 10 40 erg s −1 .Model parameters (bottom): black hole mass 10 8 M , L/L Edd = 1, inclination angle 30 degree, and luminosity of the corona is 3.78 × 10 46 erg s −1 .

Fig. 10 .
Fig.10.Different X-ray light curve generated for different breaking time.Red, blue, and green are for breaking times of 7.5 days, 75 days, and 750 days, respectively.

Fig. 14 .
Fig.14.Light curve for λ = 1258.92Å generated for different BLR contribution.We used half-Gaussian for BLR profile σ = 5 days and we used an X-ray light curve for breaking time 75 days.Parameters: black hole mass 10 8 M , Eddington ratio = 1.0, height: h = 5r g , and viewing angle: i = 30 degrees.