© The Authors 2025

1. Introduction

The central components of active galactic nuclei (AGNs) include a massive black hole, a compact X-ray emitting region, a relatively cold Keplerian accretion disk, and the broad line region (BLR) (see, e.g., Krolik 1999, for a basic compendium). The Keplerian disk serves as the primary source of continuum emission, while the BLR predominantly produces broad emission lines, along with a minor contribution to the continuum. In terms of spatial extent, the innermost regions of the accretion disk span approximately 10−100 gravitational radii (r_g), whereas the BLR extends outward at distances of order of ∼1000 r_g. Beyond the BLR, a dusty molecular torus is located at even larger radii (∼10⁴ to 10⁵ r_g). Since the earliest studies of active galaxies by Seyfert (1943) and the discovery of quasars by Schmidt (1963), there has been tremendous progress in our understanding of the structure of AGNs (see e.g., Krolik 1999, for a review).

The innermost regions of an AGN are too compact (spanning just a few microarcseconds in the case of the BLR) to be spatially resolved directly. However, their structure can be probed using intrinsic flux variability. Among the various wavelengths, X-ray variability is the fastest, originating from the most compact region near the black hole. The reprocessing of this X-ray emission at longer wavelengths, observed through light echo studies, provides insights into the structure and dynamics of the surrounding material, a technique known as reverberation mapping (RM; Blandford & McKee 1982; Peterson 1993). This technique has been widely applied to AGNs (Sergeev et al. 2005; McHardy et al. 2014; Mudd et al. 2018; Homayouni et al. 2019; Yu et al. 2020; Guo et al. 2022), particularly in studies of continuum variability, which is believed to arise from the irradiated accretion disk.

Although our understanding of AGN structure is still evolving, significant progress has enabled the exploration of AGNs as potential tools for cosmology. Unlike standard candles, AGNs exhibit luminosities that span several orders of magnitude. Nevertheless, the luminosity distance to individual AGNs can be estimated either through geometric methods or by standardizing them in a statistical sense (see, e.g., Czerny et al. 2018, for reviews). Early successful approaches relied on the nonlinear relation between UV and X-ray luminosities (Risaliti & Lusso 2015; Lusso et al. 2025) or the correlation between the BLR radius and the source luminosity (Martínez-Aldama et al. 2019; Cao et al. 2022, 2024, 2025). However, these methods require external scaling and without it, we can only constrain the curvature of the Hubble diagram. As such, they are not directly suitable for determining the Hubble constant.

The Hubble constant can instead be derived using strong gravitational lensing techniques (see Suyu et al. 2024, for a recent review) or through a combination of RM of the BLR with angular size measurements obtained from high-resolution interferometric observations (Li et al. 2022, 2025). While both methods are promising, they demand extensive, dedicated observational campaigns that are currently feasible only for a limited number of AGNs; specifically, nearby, bright sources at low redshift where sufficient spatial resolution can be achieved.

The primary goal of this work is to develop a practical method for measuring the Hubble constant using a single mean spectrum of an AGN combined with photometric RM results at multiple wavelengths. High-quality photometric RM measurements are already available for a number of sources (Fausnaugh et al. 2016; Cackett et al. 2018; Edelson et al. 2024; Prince et al. 2025), and a substantial influx of data is expected from large-scale time-domain surveys, such as the Zwicky Transient Facility (ZTF) and the Vera C. Rubin Observatory’s Legacy Survey of Space and Time (LSST).

The idea of using photometric RM to measure the Hubble constant is not new; first proposed by Collier et al. (1999), it has yet to produce successful results. The method is based on the standard accretion disk model of Shakura & Sunyaev (1973), which predicts a direct relationship between the time delay at a given wavelength and the monochromatic luminosity: intrinsically brighter sources are physically larger and thus produce longer time delays. More specifically, the model anticipates that the time delays (τ) should scale with wavelength (λ) as τ ∝ λ^4/3. While all coefficients in this relation are determined by the model, the viewing angle (i) to the source remains a free parameter and must be estimated independently.

While observational results have roughly confirmed this trend, they have not matched the predicted normalization, as observed delays continue to exceed theoretical expectations (e.g., Cackett et al. 2007; Shappee et al. 2014; Kokubo 2018; Guo et al. 2022; Mandal et al. 2025). This discrepancy, known as the accretion disk size problem, was first recognized in optical microlensing studies of gravitationally lensed quasars (Morgan et al. 2010). Until recently, remained an ongoing challenge in our understanding of AGN accretion disk structure.

Recent intensive RM campaigns (Edelson et al. 2015; De Rosa et al. 2015; Fausnaugh et al. 2016; Pei et al. 2017; Horne et al. 2021; Lu et al. 2022; Prince et al. 2025) have revealed that the observed optical continuum is not solely emitted by the accretion disk. Korista & Goad (2001) first noted that radiation reprocessed by the BLR contributes not only to the formation of emission lines but also to the diffuse continuum emission. This finding was later confirmed by Fausnaugh et al. (2016) and further emphasized by Cackett et al. (2021), Netzer (2022), who demonstrated that the diffuse BLR contribution can significantly dominate the measured continuum time delays. This effect poses challenges for using continuum time delays as a method to determine the Hubble constant, H₀, as originally proposed by Collier et al. (1999). For instance, Cackett et al. (2007) derived an H₀ value of 15 ± 3 km s⁻¹ Mpc⁻¹; however, at the time, not all aspects of the physical mechanisms driving continuum time delays were fully understood.

The contribution of the BLR to the continuum time delay implies that the use of photometric continuum RM of AGNs for measuring the Hubble constant requires not only a reliable accretion disk model, but also a robust model of the BLR structure. Importantly, such a BLR model must be fully determined by the absolute luminosity of the source. However, fully parametric BLR models (e.g., Baldwin et al. 1995; Li et al. 2013; Pancoast et al. 2014; Grier et al. 2017) might not be adequate in this context, as the ability to recover the absolute luminosity of the source could be compromised, unless wavelength-resolved spectroscopic RM data are available. Without such spectral resolution, key constraints on the BLR geometry and its luminosity dependence may be lost, limiting the accuracy of distance measurements based on photometric continuum RM.

We thus propose the use of the Failed Radiatively Accelerated Dusty Outflow (FRADO) model of the BLR (Czerny & Hryniewicz 2011; Naddaf et al. 2021; Naddaf & Czerny 2022). The model assumes BLR formation based on radiation pressure acting on dust in disk regions where the disk atmosphere temperature falls below the sublimation threshold. It was proposed to explain the vertical rise of material responsible for low-ionization lines, such as Hβ. In this study, we thus investigate the reprocessing of central radiation through a combination of two reprocessors: the Keplerian disk, which we modified to include the warm corona in the innermost part, and the BLR, using the BLR response derived from the FRADO model. This gave us access to a three-dimensional (3D) distribution of BLR clouds, which we then combined with a spectral shape model of the BLR generated using the CLOUDY code (version 23.00; Chatzikos et al. 2023). We fit this model to the observed spectrum and time delays of NGC 5548 and evaluated its potential for determining the distance to the source and the Hubble constant. We tested the model for a specific source and treated the results as the pilot study: a base for further model developments and the later use for other objects.

NGC 5548 (z = 0.017175), a Seyfert 1 galaxy, is among the most extensively studied AGNs. Its brightness, relatively close distance, and significant variability made it an early candidate for studying the inner, unresolved structure of an active nucleus (e.g., Clavel et al. 1991; Peterson et al. 1991; Krolik et al. 1991; Rokaki et al. 1993). The recent high-quality photometric monitoring (Fausnaugh et al. 2016), broadband spectral analysis (Mehdipour et al. 2015), measured Hβ time delay from spectroscopic RM (Pei et al. 2017), and wavelength-resolved RM that captured the detailed response of the BLR (Horne et al. 2021) collectively make NGC 5548 an ideal laboratory for testing our methodology.

The structure of this paper is as follows. In Section 2, we describe the data used for fitting. Section 3 outlines our methodology, which involves generating the response function of the accretion disk using the lamp-post model and incorporating the BLR contribution via FRADO and CLOUDY to estimate the combined time delays. Section 4 presents the results obtained from different models, detailing the simultaneous fitting of the time lag and spectral energy distribution (SED) while assessing the strengths and limitations of each approach. Finally, in Section 5, we estimate the luminosity distance based on our best-fitting model. The discussion and conclusions are presented in Sections 6 and 7, respectively.

2. Spectral shape and continuum time delay data of NGC 5548

NGC 5548 (z = 0.017) is one of the most extensively studied AGNs across multiple wavelengths. It was the primary target of the AGN Watch international collaboration (Peterson et al. 2002), which began investigating its emission line variability as early as 1987. Over the years, its optical and UV variability has been the subject of numerous studies (Peterson et al. 1991, 1992, 2002, 2004; Korista et al. 1995). Observations from both ground-based and space-based telescopes have provided critical insights into accretion disk dynamics and emission line variability. Notably, optical variability in both the continuum and emission lines was first reported by de Vaucouleurs & de Vaucouleurs (1972) in the 1970s, while early UV observations commenced in 1988 with the International Ultraviolet Explorer (IUE; Ulrich & Boisson 1983). More recently, intensive multiwavelength monitoring campaigns (Edelson et al. 2015; De Rosa et al. 2015; Fausnaugh et al. 2016; Pei et al. 2017; Horne et al. 2021; Lu et al. 2022) have significantly deepened our understanding of the source.

These studies have revealed that the measured continuum time delays do not follow a simple power-law trend with wavelength. Instead, they show clear evidence of interaction with an extended BLR, where different emission lines originate. Consequently, these delays are often approximated by mean time delays, though some emission line delay profiles exhibit double-peaked features. This complexity opens up new avenues for testing specific models of the BLR structure, as well as for explaining the observed inter-band continuum time delays.

In this study, we utilized the broadband SED and global parameters from Mehdipour et al. (2015), covering the period from 2013 to 2014. The broadband SED, corrected for internal extinction, starlight contamination, the Balmer continuum, and the Fe II pseudo-continuum, was first presented by Mehdipour et al. (2015). Subsequently, Kubota & Done (2018) fit this spectrum to an AGN spectral model, revealing that only the outermost part of the accretion flow corresponds to a standard disk. Their results further indicated that the optical/UV emission primarily originates from a region dominated by the warm corona (see their Figure 4).

Therefore, to systematically model the observed inter-band continuum time delays, we propose a step-by-step approach. We select three representative setups to illustrate how the spectral fit and time-delay measurements depend on the pre-set parameters and details of the geometry.

First, we assumed a standard accretion disk model (model A). This choice is motivated by the delay fitting for NGC 5548 done by Kammoun et al. (2021a). The source was fitted by Kammoun et al. (2021a) without assuming an inner hot flow, only with the coronal height, black hole mass (M_BH), spin, accretion rate, and the X-ray luminosity as free parameters. The color correction was set to 2.4. We modified the approach: we fixed M_BH and the Eddington ratio, but we allow for the inner hot flow and the contribution of the BLR. Our color correction is the same and no warm corona is present. Since we also fit the spectral shape, we included a starlight component in the model.

Next, we considered an alternative scenario (model B) based on the spectral decomposition proposed by Kubota & Done (2018). In this case, the warm corona plays an important role in the spectral fit. We adopted the value of the warm corona temperature from the corresponding paper. To calculate the time delay caused by reprocessing, we assumed that the optical/UV variability arises solely from the reprocessing of X-rays in the outer cold disk. This is motivated by the fact that warm corona with a temperature of approximately ∼10⁶ K is too optically thick to absorb and re-emit incident hard X-rays in the UV/X-ray bands; although it does still contribute to the UV continuum. The model was supplemented with the BLR contribution to the delays and the spectrum, with starlight included.

Finally, we introduced model C, the most general spectral decomposition, based on Figure 5 of Mehdipour et al. (2015). Within its framework, a number of geometrical and physical parameters are optimized to fit the spectral and time delay data. This includes the warm corona temperature and its optical depth, as well as contribution from BLR and starlight.

The disk parameters used for models A, B, and C are summarized in Table 1, with more detailed descriptions provided in Section 4.4. To test these models, we utilized the continuum time delay data from Fausnaugh et al. (2016), obtained during the STORM campaign from December 2013 to August 2014. This campaign integrated optical light curves data from 16 ground-based observatories, across the B, V, R, I filters, as well as the SDSS −u, g, r, i, z filters, coupled with ultraviolet data from the Hubble Space Telescope (HST) and Swift instruments. The interband time delays were measured relative to the HST light curve at 1367 Å, providing a crucial dataset for evaluating the proposed disk reprocessing scenarios.

3. Method

Understanding the origin of continuum time delays in AGNs is crucial for probing the structure of the central engine and its surrounding medium. However, fitting these delays solely as a result of accretion disk reprocessing has often led to the so-called disk-size problem, as discussed above. Motivated by previous studies (Korista & Goad 2001; Netzer 2022; Jaiswal et al. 2023; Beard et al. 2025), we extended this framework by considering the reprocessing of central flux by both the accretion disk and the BLR medium. Given our twofold motivation, we adopted a theoretical modeling approach that accounts for all relevant aspects. This allows us not only to test the BLR model we employed, but also to explore the determination of the Hubble constant, as the relatively small number of model parameters enables us to treat the luminosity distance as an unknown parameter.

Building on this foundation, our model of an active nucleus comprises multiple key components: a standard accretion disk in the outer parts of the flow, a warm Comptonizing corona, and a hot corona in the inner parts of the flow, the BLR region, and the contribution of the host galaxy to the total spectrum. In the following sections, we examine each of these components in detail.

3.1. Hot and warm corona

The hot corona, responsible for hard X-ray emission, was modeled by Mehdipour et al. (2015). In our approach, we adopted the proposed power-law shape and normalization for the observed spectrum when fitting the warm corona. We did not model the hot corona; instead, the power-law component was treated as part of the broadband spectrum that irradiates the BLR clouds, as discussed in the next section. It is also included in the fit to the global spectrum, in the X-ray region. However, we neglect the contribution of hard X-rays to the optical/UV part of the spectrum.

Unlike Mehdipour et al. (2015), we modeled the warm corona with a different approach. In our framework, the inner and outer radii of the warm corona are treated as parameters. The outer radius of the warm corona coincides with the inner radius of the standard (outer) disk. The Comptonization process in the warm corona is computed using the analytical formulae from Sunyaev & Titarchuk (1980) and is parameterized by the optical depth and electron temperature. Since the warm corona covers the inner disk, we do not assume a single temperature for the soft photons. Instead, we determine the temperature of the underlying cold disk at each radius and adjust it by the Comptonization amplification factor. These calculations were performed iteratively, ensuring that the soft photon temperature spans a continuous range, while preserving the total energy from the accretion flow. The coronal parameters (optical depth and the temperature) were assumed to be radius-independent, as we did not have firm predictions for these parameters from the warm corona physics. This methodology is closely aligned with the approach of Kubota & Done (2018), while the computational code used for this purpose was originally presented in Czerny et al. (2003). The warm corona can, thus, be fully modeled, with its spectral shape and normalization determined by global fitted parameters, including the optical depth, temperature, the disk radius at which the warm corona develops, and the accretion rate. These parameters are incorporated in the spectral fitting, but not included separately in BLR reprocessing where the incident radiation is made at the basis of the broadband SED, instead of a separate modeling of the disk and the warm corona. This choice allows us to avoid the computational complexity of calculating the reprocessed emission using the CLOUDY code (Chatzikos et al. 2023) at every step of the simulation. This is stated more explicitly later on in this work.

Additionally, we did not include irradiation of the warm corona by hard X-rays, as the high temperature and large optical depth of the warm corona would cause the incident radiation to be predominantly scattered rather than absorbed and reprocessed. Consequently, no significant thermalization was expected, although some high-ionization lines, such as iron lines in the X-ray spectrum, could still be present (Petrucci et al. 2020; Ballantyne et al. 2024).

3.2. Cold disk reprocessing

The irradiating flux is thermalized only in the outer, cold regions of the standard accretion disk. To further explore this process, we simulate the lamp-post model to generate the disk’s response function, assuming a pulse of short duration. In this simulation, we excluded corrections for general relativity (GR) and energy-dependent reflection effects, as introduced by Kammoun et al. (2023, 2024), Papoutsis et al. (2024), to maintain a relatively simple code suitable for data fitting. As a result, the incident radiation from the corona was assumed to be fully absorbed by the cold disk, leading to a localized increase in temperature. We briefly comment on the implications and limitations of this simplification in Section 6.

To accurately model the disk reprocessing, we employed a non-uniform radial grid extending from the innermost stable circular orbit (r_isco) to the outer disk radius (r_out), ensuring sufficient resolution at smaller radii. The radial grid spacing follows the relation $\mathrm{d}r=0.085\times {(\frac{r}{r_{\mathrm{isco}}})}^{0.85}$, where r is the radial distance from the central source in units of gravitational radius, r_g. For each radial position, the azimuthal angle, ϕ, is sampled using the relation $d\varphi =\frac{1.57}{N_{\mathrm{div}}}$, where N_div = 3800 ensures sufficient resolution across the entire radial range. The surface element for each (r, ϕ) coordinate is defined as ds = r ⋅ dr ⋅ dϕ in units of r_g². To facilitate further calculations, we converted (r, ϕ) into Cartesian coordinates, assuming a geometrically thin disk with negligible height.

For each (x, y) coordinate, we computed the total time delay, τ_total(x, y), as the sum of two components: the delay from the corona to the disk, τ_d(r), and the delay from the disk to the plane intersecting the equatorial plane at (r = r_out, ϕ = 0), denoted as τ_do(x, y). The delay is influenced by key parameters, including the lamp-luminosity (L_x), M_BH, accretion rate (Ṁ), coronal height (h), and the inclination angle of the system (i). By default, we assumed a cold Keplerian disk extending down to the ISCO at 6 r_g, with an outer radius of 10⁴ r_g. However, we also considered cases with different inner and outer radii, which can be set through the parameter r_in.

The non-irradiated disk is characterized by its flux, F_{non − irradiated}, and temperature, T_disk, as defined in Equations (1) and (2), respectively. When the contribution from the corona is included, the resulting total flux and temperature are given by Equations (3) and (4). For a specified delay, we use the combined flux from both the disk emission and irradiation, F_{disk + irradiation}, as described in Equation (3). We then converted this total flux into an effective temperature, T_eff, according to Equation (4). This temperature is subsequently used to compute the blackbody emission.

$$\begin{array}{cc}& F_{\mathrm{non}-\mathrm{irradiated}}(r,t+{\tau }_{\mathrm{d}}(r))=\left(\frac{3GM\dot{M}}{8\pi r^3}(1-\sqrt{\frac{r_{\mathrm{in}}}{r}})\right),\hfill \end{array}$$(1)

$$\begin{array}{cc}& T_{\mathrm{disk}}(r,t+{\tau }_{\mathrm{d}}(r))={\left[\left(\frac{3GM\dot{M}}{8\pi r^3{\sigma }_B}(1-\sqrt{\frac{r_{\mathrm{in}}}{r}})\right)\right]}^{\frac{1}{4}},\hfill \end{array}$$(2)

$$\begin{array}{cc}& F_{\mathrm{disk}+\mathrm{irradiation}}(r,t+{\tau }_{\mathrm{d}}(r))=\left(\frac{3GM\dot{M}}{8\pi r^3}(1-\sqrt{\frac{r_{\mathrm{in}}}{r}})\right)+\left(\frac{L_x(t)h}{4\pi r^3}\right),\hfill \end{array}$$(3)

$$\begin{array}{cc}& T_{\mathrm{eff}}(r,t+{\tau }_{\mathrm{d}}(r))={[\left(\frac{3GM\dot{M}}{8\pi r^3{\sigma }_B}(1-\sqrt{\frac{r_{\mathrm{in}}}{r}})\right)+\left(\frac{L_x(t)h}{4\pi r^3{\sigma }_B}\right)]}^{\frac{1}{4}}·\hfill \end{array}$$(4)

Additionally, we accounted for a color correction factor, f_c, to the disk temperature, which raises the temperature of the disk atmosphere to T_c, where f_c = 1 corresponds to no color correction applied, as described below (Shimura & Takahara 1995):

$$\begin{array}{c}\hfill T_{\mathrm{c}}=f_{\mathrm{c}}{T}_{\mathrm{eff}}.\end{array}$$(5)

To compute the full spectrum for a specific differential area, we used Planck’s formula and store the results in a photon table, represented as a 2D matrix denoted as P(t′,λ). This involves applying the calculated time delay corresponding to the disk position and wavelength. Each element in the photon table represents a unique delay and wavelength λ, with the λ values ranging from 1000 to 10 000 Å and selected using a logarithmic scale grid for simulation. We also introduce a color correction as a free parameter of the model.

Next, we generate the response function by sending a very short light pulse of Δ_t = 0.05 days duration and normalizing the result by the incident bolometric luminosity. The pulse duration is chosen to balance accuracy and temporal resolution: a very short pulse can introduce significant errors in the temperature enhancement, while a very long pulse tends to smear the signal and reduce time resolution, especially at the shortest wavelengths. The response function for a bare disk is given by the following equation,

$$\begin{array}{c}\hfill \psi (t,\lambda )=\frac{1}{{\mathrm{\Delta }}_t{L}_x}\frac{1}{f_{\mathrm{c}}^4}{\displaystyle {\int }_{S_{\mathrm{disk}}}}{B}_{\lambda }(T_{\mathrm{c}}(r(t\prime -{\tau }_{\mathrm{do}}(x,y))))\mathrm{d}s.\end{array}$$(6)

The time delay in the disk has been described previously by Jaiswal et al. (2023). Our adopted method is indeed simpler than the approaches presented by Kammoun et al. (2021a,b), which include corrections for GR and disk albedo; however, it is sufficiently accurate for the current purpose of this pilot study.

3.3. BLR structure and reprocessing

To introduce the contribution from the BLR, we first describe its properties. In this study, we did not parameterize the structure of the BLR using arbitrary numbers or functions; instead, we calculated it based on the FRADO model, which was qualitatively proposed by Czerny & Hryniewicz (2011) and further developed in several subsequent studies (Czerny et al. 2017; Naddaf et al. 2021, 2023; Naddaf & Czerny 2022, 2024).

The basic idea behind the model is to associate the low-ionization parts of the BLR that are responsible for emission lines, such as Hβ, Mg II, and Fe II, with dust-driven massive winds. In stellar environments, such winds are typically denser, possess higher optical depths, and exhibit significantly slower outflow velocities. Similarly, in AGNs, dust-driven winds are expected to be launched from the outer disk, where the effective temperature falls below the dust sublimation threshold. This temperature constraint naturally defines the inner radius of the BLR. The nature of the wind depends on M_BH and Eddington ratio: for sufficiently massive black holes and high accretion rates, the wind escapes, contributing to the BLR outflow structure; otherwise, it forms a failed wind, in which the material ultimately falls back onto the disk. Additionally, irradiation from the central parts of the disk plays an increasingly important role in shaping the dynamics and ionization state of the wind as the cloud elevates above the disk plane.

Specifically, we utilized the code from Naddaf et al. (2021), which provides a detailed numerical description of the wavelength-dependent cross-section for dusty particles. According to this model, radiatively dust-driven pressure lifts the clouds from the disk surface, while preserving their angular momentum derived from the Keplerian motion of the disk surface. As a cloud is lifted, it becomes increasingly illuminated by the inner parts of the disk. However, if the cloud becomes too hot, the dust evaporates, allowing the cloud to continue its motion along a ballistic orbit. The global parameters of the model, such as M_BH, Eddington rate, and metallicity, govern the behavior of the BLR. For lower black hole masses, lower Eddington ratios, and lower metallicities, the clouds form a failed wind, while in the opposite case, a fraction of the clouds may form an escaping wind. Hence, we emphasize that the inner and outer radii of the BLR, along with the statistical distribution of the clouds, are governed by global parameters, such as M_BH and the Eddington ratio.

We assumed a universal value for the dust sublimation temperature of 1500 K, representative for all grain species and sizes. It was suggested by Barvainis (1987) based on observational data for PG quasars. It is frequently used for zero-order approximations, although it is well known that actually the sublimation temperature depends on the chemical composition as well as size of the dust grains (Draine & Lee 1984; Baskin & Laor 2018). However, calculating the range of evaporation temperature of each grain is time consuming and not appropriate for a pilot study. In addition, the chemical composition of grains in AGN is not firmly established, and the presence of the graphite with usual grain size distribution is highly unlikely (e.g., Czerny et al. 2004; Gaskell et al. 2004). We further discuss this issue in Section 6.3.

With knowledge of the global parameters of the source, we can determine the 3D locations of statistically representative clouds within the BLR. The code calculates the 3D trajectories of clouds launched from the disk. As described in Naddaf et al. (2021), clouds are then located at these trajectories proportionally to time they spent at each part of the trajectory. This information is then combined with the emissivity law for specific emission lines, particularly low-ionization lines such as Hβ, Mg II, or Fe II. In our model, we assume that the emissivity weight of each cloud is influenced by its distance from the central disk: clouds closer to the center receive more radiation, resulting in a higher emissivity weight, while clouds farther away receive less radiation, leading to a lower weight, as defined in Equation (7),

$$\begin{array}{c}\hfill w\propto d^{-2},\end{array}$$(7)

where d is the distance of a cloud from the center. This scaling of the emissivity, which inversely depends on the square of the distance, neglects the local efficiency of converting incident flux into BLR spectrum. However, since our model does not yet predict the local density, this simplification is reasonable at this stage of development. Thus, we fixed the local density of the cloud at 10¹¹ cm⁻³ and the column density at 3 × 10²³ cm⁻².

By supplementing the cloud distribution with emissivity, we can determine the emission line profile for any viewing angle relative to the observer, as demonstrated in Naddaf & Czerny (2022). Next, we can focus on calculating the response function ψ_BLR(t) for the BLR based on this distribution. This is accomplished by calculating the time delays for each representative cloud as observed, using a method analogous to that applied for disk emission. At this stage, we assumed that the emissivity of the clouds is universal (albeit wavelength-dependent), determined through photoionization calculations for a representative cloud.

3.4. Photoionization calculations using CLOUDY

To determine the BLR spectral shape, we performed photoionization calculations using the code CLOUDY, version C23.00 (Chatzikos et al. 2023). In these computations, we simplify the 3D cloud distribution by replacing the entire ensemble of clouds with a single representative cloud, positioned at the radius inferred from the time delay. For the analysis at this established distance, we adopted an incident radiation luminosity of log L (erg/s) = 44, and a BLR distance from the central source of log r (cm) = 16 (Dalla Bontà et al. 2020), a constant hydrogen gas density of log n_H (cm⁻³) = 11, and a column density of log N_H (cm⁻²) = 23.5, following the prescriptions of Korista & Goad (2001), Panda et al. (2022). The adopted distance in our representative case corresponds to the measured time delay of the Hβ line in NGC 5548. The assumed bolometric luminosity is also representative of this source, based on the redshift-inferred distance within the framework of standard cosmology. In this pilot study, we do not explore a wide range of gas densities or column densities in detail, as photoionization models, particularly single-zone models, still do not perfectly represent the full complexity of BLR properties (e.g., Netzer 2020; Pandey et al. 2025; Floris et al. 2025). We assumed a metallicity five times higher than solar, as this was favored by comparisons between the FRADO model and quasar spectra (Naddaf et al. 2021). Independent studies of quasar metallicity also support super-solar values (e.g., Śniegowska et al. 2021), and in the specific case of NGC 5548, high metallicity was advocated by Netzer (2020) based on the line emissivity analysis. For simplicity, we neglected turbulence in the medium.

The shape of the incident SED is taken from the file “NGC5548.sed” available in the CLOUDY database, which corresponds to the SED derived by Mehdipour et al. (2015). Their study estimated the starlight contribution to NGC 5548 using HST observations, providing a well-constrained incident SED that we adopt for our calculations.

This is clearly an oversimplification, in two aspects. Each cloud receives the radiation of the same shape, only scaled down with the distance as specified in Equation (7). This radiation includes both the hard X-ray radiation as well as the warm corona and the disk emission. Apart from the spectral shape, we also did not differentiate here between the arrival time from the hot corona, and warm extended corona. Accounting for such differences would require us to model the BLR response separately to the hard X-ray emission from the hot corona and to the emission from the warm corona, followed by merging several resulting transfer functions. As this is a pilot study, we have not undertaken that level of complexity at this stage, although it is feasible in principle, given that our spectral decomposition separates the contributions from all three components. For the current analysis, we adopted a fixed SED shape, which is a simplification that is unlikely to significantly affect our results. The hard X-ray power law from the hot corona is consistent with SED. The warm corona was fitted in way that allowed us to reproduce the spectrum. A slight error on the time delay, coming from the fact that some photons should come directly from the hot corona (instead of the disk) should not be significant when taking into account the much larger size of the BLR, as compared the disk.

The resulting BLR emissivity profile, ϵ(λ), was then incorporated into the final delay model. When we later accounted for the assumption of an unknown luminosity distance, the irradiating flux was adjusted accordingly to ensure consistency with the observed hard X-ray flux.

Since CLOUDY does not compute higher-order Balmer lines, we incorporated a scaled component into the spectrum, following the method of Kovačević et al. (2014). This modification improves the spectral fit, as demonstrated in Pandey et al. (2025), although it only affects the wavelength range between 3646 Å and 4000 Å in the rest frame. The shape of this component is taken directly from Pandey et al. (2025); however, the emissivity still appears too low to fully capture the gradual decline observed beyond the Balmer edge in the data. More advanced modeling may be required to improve the fit in that region.

When we performed our computations of the model C for a fixed luminosity distance, the value of the incident flux in CLOUDY computations was fixed to correspond to the observed luminosity and the time delay to the BLR, as detailed above. However, later on (as described in Section 5) when the luminosity distance is treated as unknown, we scaled the incident bolometric luminosity to the fitted accretion rate, which is a function of the luminosity distance.

3.5. Starlight contribution

The starlight contribution in NGC 5548 was estimated by Mehdipour et al. (2015) and we adopted the same level in some of our models. However, in other models, we allowed for flexibility in adjusting the starlight level. To represent the starlight profile, we used the template of an Sa galaxy from Kinney et al. (1996), which is available in the Kinney-Calzetti Spectral Atlas¹.

3.6. Combined time delay and NGC 5548 data fitting

After analyzing each component of the central regions assumed in our model for NGC 5548, we derived the final response function by incorporating both disk reprocessing and BLR contributions, as expressed in the following equation,

$$\begin{array}{cc}\hfill \psi (\lambda ,t)=& (1-f_{\mathrm{BLR}}){\psi }_{\mathrm{d}}(\lambda ,t)\hfill \\ \hfill & +f_{\mathrm{BLR}}ϵ(\lambda ){\displaystyle {\int }_{t_0}^{t_{\max }}}{\psi }_{\mathrm{d}}(\lambda ,t-t\prime ){\psi }_{\mathrm{BLR}}(t\prime )\mathrm{d}t\prime ,\hfill \end{array}$$(8)

where the parameter f_BLR defines the BLR fraction contributing to the continuum time delay.

As mentioned in Section 3.4, we did not separately treat the photons reaching the BLR directly from the hot corona and/or from the warm corona and from the disk. Photons from the hot and warm corona were included in BLR computations as the incident continuum (i.e., they were included in the ϵ(λ)), but their arrival times were the same as for the photon disks under consideration. This was done for the efficiency of computations when fitting the data. Replacing the incident continuum adopted for photoionization modeling with three separate spectral components would additionally require repeating photoionization computations at every computational step.

We generate response functions for 100 wavelengths, covering the range from 1000 Å to 10 000 Å. Using these response functions, we could then compute the time delay at each wavelength using the following equation,

$$\begin{array}{c}\hfill \tau (\lambda )={\displaystyle \frac{\int t\psi (t,\lambda )\mathrm{d}t}{\int \psi (t,\lambda )\mathrm{d}t}}·\end{array}$$(9)

To fit the data of NGC 5548, we utilized all the parameters listed in Table 1, integrating the contributions from both the disk and the BLR described in Equation (8).

When searching for the optimal solution, we impose constraints from both the observed SED and the measured time delays. Specifically, we incorporated all 17 available time delay measurements and selected 17 representative points, evenly distributed across the entire observed SED range. To systematically evaluate and compare different solutions, we integrated these constraints into a unified criterion, expressed as

$$\begin{array}{c}\hfill {\chi }^2={\displaystyle \sum \frac{{(O_i-E_i)}^2}{{(\delta O_i)}^2}},\end{array}$$(10)

where O_i = observed data, δO_i = is the error in observed data, and E_i = value from the model. We assumed an error in the spectral data of 7%, while the errors of the delay measurements were taken from Fausnaugh et al. (2016).

4. Results

In this section, we concentrate on testing the FRADO model against the spectral and time-delay data for NGC 5548 and summarize the results obtained from our model fitting. We fixed M_BH and the bolometric luminosity in our analysis. We assumed that the distance to the source is known and corresponds to the source redshift. The issue of the Hubble constant determination is addressed separately in Section 5.

4.1. BLR properties from FRADO

The global model parameters adopted for NGC 5548, combined with the assumption of a dust sublimation temperature of 1500 K, allow us to uniquely determine the structure of the BLR. Notably, the FRADO model is specifically applicable to low-ionization lines such as Hβ, Mg II, and Fe II. In this framework, the cloud positions within the BLR are determined by the FRADO model, with parameters relevant to NGC 5548: M_BH = 5 × 10⁷ M_⊙ (Netzer 2022), an Eddington ratio of 0.02 (Crenshaw et al. 2009), and an assumed metallicity that is five times solar.

We present the distribution of clouds forming the BLR in NGC 5548 in Figure 1. The upper panel presents the overall 3D structure, revealing a complex cone-like configuration where no clouds escape to infinity due to the low Eddington rate. The lower panel provides a cross-section of this distribution, with the black line marking the thickness of the underlying Keplerian disk. This disk structure accounts for the effects of radiation pressure and opacity, as described by Różańska et al. (1999), both of which play a crucial role in shaping its complex geometry. A notable feature is the abrupt change in disk height, corresponding to the transition from the inner radiation pressure-dominated region to the outer gas-dominated region, which significantly impacts the cooling efficiency of the disk and its overall structure.

Fig. 1. Upper panel: 3D plot of cloud positions from FRADO model for NGC 5548 with MBH = 5 × 107 M⊙, L/LEdd = 0.02, metallicity Z = 5 in solar units. The axes are in units of rg. Bottom panel: Cross-section of cloud positions for y > −300 rg and y < 300 rg. The black line indicates the thickness of the Keplerian disk. Clouds form a geometrically thin complex layer above the disk.

In this model, clouds are launched from the disk surface but remain relatively close to it, with the ratio of vertical height (z) to radial distance (r) staying below 5%. This limited height is a direct consequence of the low Eddington ratio, giving the BLR a geometry that closely resembles a puffed-up disk surface. The z/r ratio in the model peaks near the inner edge of the BLR, which is defined by the dust sublimation temperature reached at approximately ∼2260 r_g for the adopted M_BH, Eddington rate, and metallicity. At this radius, the Keplerian velocity is ∼6300 km s⁻¹ for this M_BH, implying that the full width at half maximum (FWHM) of the emission lines could reach up to twice this value. This result based on FRADO model is consistent with the Hβ FWHM of 10 161 ± 587 km s⁻¹ reported by Pei et al. (2017). However, contributions from larger radii act to moderate this broadening.

The position of the clouds and their velocity field allowed us to construct a 2D velocity-delay map (see Figure 2, upper panel) which we consider as representative for Hβ line, assuming a viewing angle of 40°, which is likely representative of NGC 5548 (Pancoast et al. 2014; Cappi et al. 2016; Wildy et al. 2021; Horne et al. 2021). Our velocity–delay map closely resembles the predictions of a simple flat Keplerian disk model, as the vertical velocities are minimal, no outflows are included, and no selective shielding is assumed. It shows some similarity in the overall shape to the map derived observationally by Horne et al. (2021) (see their Figures 5 and 7).

Fig. 2. Upper panel: Velocity–delay map of NGC 5548 constructed from FRADO model with inclination angle i = 40°. Color coding represents the number density of BLR clouds in each velocity-delay bin. The white dashed line shows the virial envelope corresponding to Keplerian disk-like rotation for v2 × τ = constant with MBH = 5 × 107 M⊙. Bottom panel: BLR response function generated using the cloud positions from the model shown in Figure 1, assuming an inclination angle of i = 40°.

However, the main inconsistency lies in the narrower velocity extent observed in our velocity-delay map, which shows emission up to ∼ ± 4000 km s⁻¹, while Horne et al. (2021) reported observed Hβ emission extending to ∼ ± 8000 − 10 000 km s⁻¹, although the majority of the response is concentrated within ∼ ± 5000 km s⁻¹. Note that, FRADO inherently assumes that the BLR originates beyond the dust sublimation radius, which is appropriate for low-ionization lines like Hβ and thus naturally predicts narrower velocity distributions, consistent with emission from outer, slower-moving regions of the BLR. While, the broader velocity distributions observed by Horne et al. (2021) are shaped by contributions from the inner, dust-free high-ionization zones (dominated by CIV, He II lines) of the BLR, where gas orbits closer to the black hole at higher Keplerian speeds.

Nevertheless, the exact comparison is difficult since the observational Hβ map in Horne et al. (2021) seems to be strongly contaminated by He II emission. Notably, while the red wing of the observed Hβ shows a Keplerian structure, the blue wing appears significantly weaker, suggesting substantial He II contamination. Our model, in its present form, cannot explain this asymmetry. However, such asymmetries in the Hβ line are well-documented and known to evolve over time in a quasi-periodic fashion (Shapovalova et al. 2004). Bon et al. (2016) explained such phenomena as caused by the presence of the secondary black hole perturbing the BLR via the precession of the outer part of the disk or, alternatively, the presence of spiral waves in the disk. All such mechanisms will give asymmetry between blue and red-shifted wing of the flat BLR coexisting with the disk. In the future we could include the warped disk into FRADO model, but this is beyond the scope of the current paper.

Next, we derived the BLR response function by projecting the map onto time axis. The resulting shape, shown in Figure 2, lower panel, exhibits two distinct peaks: one at shorter time delays, around ∼2.6 days, corresponding to the clouds on the same side as the observer, and the other peak arising from the clouds on the opposite side of the black hole. The shape of our modeled response function does not perfectly replicate the observed one reported by Horne et al. (2021). In their analysis using the MEMECHO software for Hβ (see Figure A.1 for comparison), the projected response exhibits two peaks: a prominent one at approximately 2 days and a second, much smaller peak at around 25 days. While our first peak aligns well with theirs at ∼2 days, our second peak appears closer in time, at ∼11 days, and is stronger than the second peak, reported by Horne et al. (2021). Despite these differences, it is important to emphasize that our model relies solely on global parameters, without introducing arbitrary values such as inner or outer BLR radii. This absence of free parameters enhances the significance of the observed similarity between the modeled and empirical response functions. Additionally, we assess the consistency between the BLR response function recovered from our FRADO model and the Hβ emission-line response function derived from spectroscopic RM data by Horne et al. (2021) in Appendix A.

By integrating the BLR response profile derived from the FRADO model, we calculate a mean delay of 8.101 days for a viewing angle of 30° and 8.105 days for a viewing angle of 40°. This delay is measured relative to the X-rays. To compare it with delays at other wavelengths, such as the commonly used 5100 Å, it is essential to account for the net time delay between the X-rays and the specific wavelength of interest.

The measured time delay of the Hβ line relative to the 5100 Å continuum during the 2014 campaign (January to June) was reported by Pei et al. (2017) as $4.{17}_{-0.3}^{+0.36}$ days, which is shorter than predicted by the radius–luminosity relation. This Hβ delay ranks among the shortest ever recorded for NGC 5548. In contrast, the mean response function from the 1998−2001 campaign peaked around 20 days (Cackett & Horne 2006). Horne et al. (2021) later explored the nature of this shorter delay in detail, finding that the response function exhibited a secondary peak at approximately ∼25 days but began with a high value near to zero time delay. Their observed map also supports the primary response originating from the near side of the BLR, facing the observer.

To compare the predicted and observed time delays, we need to add the delay of the 5100 Å continuum relative to X-rays to the observed data. Since the exact time delay at 5100 Å was not measured by Fausnaugh et al. (2016), we averaged their measurements around this wavelength and included the hard X-ray time delay, also provided in their study. This results in an estimated observed time delay for Hβ relative to X-rays of 6.83 ± 0.53 days. Our predicted delay is slightly longer. This discrepancy may stem from the long tail in the delay, which the model accounts for but may not be fully captured in the observations. Alternatively, the delay could be shorter if there is intervening material between the black hole and the BLR, potentially shielding part of the most distant flow. Another possibility is that the dust temperature in the BLR model, as discussed by Panda et al. (2020) and Naddaf et al. (2021), differs from the assumed 1500 K.

Additionally, the model predicts an upper limit for the covering factor of the BLR. If the cloud distribution is not transparent, the region intercepts all radiation emitted within inclinations greater than the aspect ratio z/r of the cloud distribution. As shown in Figure 1, the model suggests that less than 4% of the radiation can be intercepted by the BLR, as the cloud height is generally small compared to the radius. This is notably lower than the 10% to 30% covering factor typically expected from BLRs (e.g., Baldwin et al. 1995; Korista & Goad 2019; Panda 2021; Panda et al. 2022). In contrast, the yearly variations observed in the line profile of NGC 5548 suggest that the disk can be distorted and/or precessing (e.g., Shapovalova et al. 2009, see their Figure 2). Such a departure from the disk model strictly perpendicular to the symmetry axis may affect the irradiation of the launched clouds and this way to modify the BLR emissivity.

The distribution of BLR clouds enables us to roughly estimate the shape of the Hβ line predicted by the model. In this estimation, we neglect the vertical cloud velocities and only consider the projected rotational velocities towards the observer. The resulting line shape, shown in Figure 3, reveals a two-peak structure that is even more pronounced than that observed in the data. The FWHM predicted by our model approaches 8000 km s⁻¹. In comparison, the measured FWHM of the Hβ line from the rms spectrum was 10 161 ± 587 km s⁻¹ (Pei et al. 2017). Notably, the FWHM values have varied over the years, indicating the dynamic nature of the BLR. For instance, Peterson et al. (2004) reported an Hβ FWHM of 7202 ± 392 km s⁻¹, which is in good agreement with our model estimate. These variations suggest that our model captures a plausible state of the BLR, consistent with the historical range of observed line widths.

Fig. 3. Line profile generated by the FRADO model using the parameters from Figure 1.

4.2. The BLR emissivity

To model the emission characteristics of photoionized clouds, we perform CLOUDY simulations that account for both line and continuum contributions. During our computations, we consider both line and continuum emissions, incorporating radiation from both inward and outward directions to ensure equal visibility of the illuminated and un-illuminated sides of the clouds. Among the most prominent features that emerge are the Balmer and Paschen edges, as illustrated in Figure 4. The depths of these edges are governed by the specific model parameters adopted (see the most recent study by Pandey et al. 2023).

Fig. 4. Emissivity profile (ratio of the reprocessed to the incident continuum) of a BLR cloud for the adopted parameters: log nH[cm−3] = 11, log L [erg s−1] = 44, and the BLR distance of 1016 cm.

The reprocessing calculations are performed for a single representative cloud, using the full shape of the SED as the incident radiation, as described in Section 3.4. This approach is admittedly a significant simplification, since in reality each cloud is exposed to a different incident spectrum depending on its location and orientation. This spectral variation is properly accounted for in the computation of the radiation force. However, generating comprehensive tables of local emissivity for the full range of incident spectra lies beyond the current scope of this work. Adopting such an extension is essential in future studies to improve the physical realism of the model.

4.3. Exemplary shape of the response function

After deriving the BLR response function using FRADO model and the emissivity profile from CLOUDY for NGC 5548, we compute the disk response function and the combined disk-BLR response function using Equations (6) and (8), respectively. These results are presented for nine different wavelengths in the upper and bottom panels of Figure 5. The disk response function (upper) displays a single-peaked structure, whereas the combined response function (bottom) exhibits a bimodal shape, as expected, with the first peak corresponding to direct emission from the disk, and the second arises from delayed reprocessing by the BLR. Moreover, the overall time delay associated with the combined response function is significantly longer than that of the disk alone, consistent with the findings of Netzer (2022).

Fig. 5. Upper panel: Response functions of the disk at different wavelengths. Bottom panel: Combined response functions of the disk and BLR across the same wavelengths. These representative response functions are obtained for model C. Parameters used: MBH = 5.0 × 107 M⊙, Eddington ratio = 0.015, LX = 9.68 × 1043 erg s−1, height: h = 48.29 rg, Rin = 94.87 rg, Rout = 10 000 rg, and viewing angle: i = 40°.

4.4. Time delays in representative models

This paper serves as a pilot study, focusing on three distinct global setups rather than exploring the full parameter space. In this section, we assume that the source distance can be determined from its redshift, enabling a unique conversion of bolometric luminosities to fluxes. We define these setups as model A, model B, and model C, which were shortly introduced in Section 2. Each setup is characterized by different global parameters. In all cases, M_BH was fixed; noting that in models A and C, the M_BH is derived from observational data (Netzer 2022), whereas in model B, it is adopted from the corresponding model reference (Kubota & Done 2018). All models are supplemented with the BLR components. A detailed discussion of each model follows in the subsequent sections.

4.4.1. Model A

To account for disk atmosphere effects, model A employs a standard accretion disk framework with a fixed color correction. It uses well-established literature values for M_BH and Eddington ratio, and does not include the warm corona component. With only a few free parameters, the model incorporates the BLR contribution alongside the accretion disk emission to fit the observed continuum time delays and spectral energy distribution. We minimize the number of free parameters by relying on the basic parameter values from the literature, including the M_BH, Eddington rate, color correction, and viewing angle (see Table 1 for the adopted values). Specifically, the adopted color correction of 2.4 comes from the studies of Kammoun et al. (2021a, 2023). The model does not include the warm corona, and the cold disk extends down to the ISCO (6 r_g). The BLR model is then directly determined based on the assumed M_BH, Eddington rate, and a fixed metallicity of 5, but the covering factor is left as a free parameter of the model. The other key parameter is the coronal height. The resulting predictions are presented in Figure 6.

Fig. 6. Model A. Upper panel: Observed delay from Fausnaugh et al. (2016) is represented in black. The delay calculated from the disk response function alone is shown in red, while the delay calculated from the combined response function of the disk and BLR is shown in brown. Bottom panel: Observed SED from Mehdipour et al. (2015) is shown in blue. The disk SED is represented in red, while green illustrates the disk SED combined with the emission lines. The starlight contribution is shown in black color. The final SED, incorporating contributions from the disk, star, and emission lines, is shown in purple. The selected points used for estimating the χ2 value are represented by brown squares. Parameters: MBH = 5.0 × 107 M⊙, Eddington ratio = 0.02, LX = 1.26 × 1044 erg s−1, height: h = 20 rg, Rin = 35 rg, Rout = 10 000 rg, color correction = 2.4, and viewing angle: i = 40°.

The model reproduces the time delay quite well, the Balmer edge at ∼3600 Å is well reproduced, the region of the Paschen jump is not so well fitted. The covering factor and the coronal height optimizing the fit of the time delay are given in Table 1. Below, we outline the key advantages and limitations of this model.

Model advantages:

• three free parameters: starlight level, inner disk radius, and BLR fraction;

• the Hβ delay is approximately recovered;

• the double peak shape is approximately recovered;

• the time delay is well-modeled.

Model limitations:

• The optical/UV SED is inconsistent with the observed data.

In model A, the discrepancy between the observed and modeled SED is substantial in the blue-shaded part and cannot be resolved by adjusting any model parameters without making the time delay fit worse. The required lamp height leads to a very blue spectrum in the UV due to the color correction applied. The time delay is primarily driven by the disk response, with the BLR contribution appearing only near the Balmer edge. This is contrary to the findings of Netzer (2022). This discrepancy arises because the disk contribution is underestimated due to the fixed color correction being set too high. While this increases the time delay by pushing the fixed-color emission outward, it simultaneously shifts the disk spectrum toward the far-UV. Consequently, we do not include the SED fit from model A in the χ² minimization, as further adjustments would not improve the fit.

The starlight contribution in our decomposition is 6.74 × 10⁻¹⁵ erg s⁻¹ cm⁻² Å⁻¹, comparable to reported value of 6.2 × 10⁻¹⁵ erg s⁻¹ cm⁻² Å⁻¹ (Mehdipour et al. 2015). Furthermore, this model fails to reproduce both soft and hard X-rays, as it includes only a cold disk. As a result, there is an inconsistency between the shape of the incident continuum assumed in the CLOUDY modeling and the continuum produced by the model itself.

4.4.2. Model B

In model B, we adopt the framework of Kubota & Done (2018), where the innermost part of the accretion disk is replaced by a hot corona. In our setup, the standard accretion disk extends from 151 r_g to 282 r_g, and we do not apply any color correction since the disk in this region remains cold. The global system parameters and the inclination angle are consistently taken from Kubota & Done (2018). The free parameters in the model are the BLR contribution, f_BLR, and the level of starlight. The coronal height is fixed at the value of the inner radius of the warm corona, namely, 43 r_g.

The results of this model are presented in Figure 7. As in model A, we employ the same BLR response function (see Section 4.1). However, this approach may introduce minor inconsistencies due to slight differences in M_BH and accretion rate. While the time delay fit remains broadly satisfactory, the SED is again not well reproduced. The reduced inner disk radius results in a steep decline in the disk spectrum at shorter wavelengths.

Fig. 7. Model B. Upper panel: Observed delay from Fausnaugh et al. (2016) is represented in black. The delay calculated from the disk response function alone is shown in red, while the delay calculated from the combined response function of the disk and BLR is shown in brown. Bottom panel: Observed SED from Mehdipour et al. (2015) is shown in blue. The disk SED is represented in red, while green illustrates the disk SED combined with the emission lines. The starlight contribution is represented in black, while the warm corona is shown in orange. The final SED, incorporating contributions from the disk, star, warm corona, and emission lines, is shown in purple. The selected points used for estimating the χ2 value are represented by brown squares. Parameters: MBH = 5.5 × 107 M⊙, Eddington ratio = 0.02, LX = 1.247 × 1044 erg s−1, height: h = 43 rg, Rin = 151 rg, Rout = 282 rg, color correction = 1.0, and viewing angle: i = 45°.

In this configuration, the total time delay is predominantly influenced by the BLR, with only a minimal contribution from the accretion disk. Consequently, the contribution from BLR to the time delay is significantly higher than in model A. Additionally, the amount of starlight in this case is 5.39 × 10⁻¹⁵ erg s⁻¹ cm⁻² Å⁻¹ (less than the corresponding value in model A).

Model advantages:

• only two free parameters: starlight level, and BLR fraction;

• the Hβ delay is approximately recovered;

• the double peak shape is approximately recovered;

• the time delay is reasonably well modeled.

Model limitations:

• Although the high-energy portion of the SED is reasonably well recovered, this model fails to accurately fit the 3000−4000 Å wavelength range, which includes the Balmer jump, and shows some inconsistencies at longer wavelengths in the optical regime. This is related to the limited flexibility due to a small number of free parameters;

• The modeled time-delay spectrum is predominantly influenced by the BLR, with negligible contribution from the accretion disk, in contrast to observed AGN time-delay spectra that exhibit significant contributions from both the accretion disk and the BLR (Mandal et al. 2025).

Fig. 8. Model C. Upper panel: Observed delay from Fausnaugh et al. (2016) is represented in black. The delay calculated from the disk response function alone is shown in red, while the delay calculated from the combined response function of the disk and BLR is shown in brown. Middle panel: Observed SED from Mehdipour et al. (2015) is shown in blue. The disk SED is represented in red, while green illustrates the disk SED combined with the emission lines. The starlight contribution is represented in black, while the warm corona is shown in orange. The final SED, incorporating contributions from the disk, star, warm corona, and emission lines, is shown in purple. The selected points used for estimating the χ2 value are represented by brown squares. Parameters used: MBH = 5.0 × 107 M⊙, Eddington ratio = 0.015, LX = 9.68 × 1043 erg s−1, height: h = 48.29 rg, Rin = 94.87 rg, Rout = 10 000 rg, and viewing angle: i = 40°. Bottom panel: X-ray-UV-optical SED fitting. The warm corona is represented by the magenta line, the hard X-ray power law is shown in dashed gold, the green line represents the sum of the two, and the black line displays the optical data.

4.4.3. Model C

This model represents a generalization of model B. Rather than fixing certain parameters to the values derived by Kubota & Done (2018), we treat them as free parameters and allow the data to constrain them through fitting. First, we allow for the variation of the radii dividing the zones of hot corona, warm corona, and outer disk. We assume that the outer radius of the cold disk is large. Then, we varied the warm corona parameters (optical depth, temperature) as well as the Eddington rate, to achieve the best fit to the broadband spectrum (including X-rays), and to time delays. Model parameters are listed in Table 1.

The time delays recovered from this model, along with the observed inter-band delays from Fausnaugh et al. (2016), are presented in Table 2. The corresponding fitted SEDs, are displayed in the upper, middle, and bottom panels of Figure 8, respectively. The Balmer edge is well reproduced in the lag-spectrum. Paschen edge is expected to be less prominent, as shown by Pandey et al. (2023) at the basis of CLOUDY simulations.

Additionally, the broadband SED is now well represented, with the starlight continuum of 6.76 × 10⁻¹⁵ erg s⁻¹ cm⁻² Å⁻¹, which is nearly identical to the value of 6.2 × 10⁻¹⁵ erg s⁻¹ cm⁻² Å⁻¹ reported by Mehdipour et al. (2015). The middle panel of the figure provides a detailed view of the optical/UV portion of the spectrum, consistent with previous presentations. Furthermore, the bottom panel displays the full X-ray-UV-optical spectral range on a logarithmic scale, allowing for clearer visualization of the broadband fit. In this representation, both the hard X-ray component and the contribution of the warm corona to the X-ray band are clearly visible. This contrasts with model A, which lacked a warm corona and therefore failed to reproduce the observed soft X-ray data. In comparison, model B, by construction from Kubota & Done (2018), is consistent with the X-ray band and thus its X-ray spectral range is not plotted in this figure. However, it is important to note that Kubota & Done (2018) did not fit the spectroscopic optical data, and as a result, the discrepancy between their model and the actual optical observations was not apparent in their analysis.

The disk contribution to the time delay in this model falls between the contributions in models A and B. Overall, the time delay is reproduced much more accurately compared to the other two models, providing a better fit to the observed data. The contribution of the corona to the spectrum is dominating the part of the spectrum below ∼2000 Å.

Model advantages:

• limited number of parameters (still larger than in model B);

• the Hβ delay is approximately recovered;

• the double peak shape is approximately recovered;

• the time delay is well modeled;

• we can successfully reproduce the Balmer jump in the observed SED and our modeled SED is consistent with the observed spectrum.

In conclusion, this model offers a significant improvement in both time delay and SED fitting, along with a more realistic representation of the BLR contribution to the observed continuum time delays and a better overall alignment with the observed spectral features.

5. Distance estimation from the generalized model C

Having considered the final fit from model C as successful, we can now assess whether our model improves the determination of the Hubble constant for this source, as done in Cackett et al. (2007). To begin, we first explore the analytical estimates and subsequently evaluate the method that ought to be employed for actual data fitting.

5.1. Analytical estimate of H₀

Collier et al. (1999) proposed a method for determining H₀ based on continuum time delays, under the assumption that the inter-band time delays are a result of a classical accretion disk structure, as given by the following equation,

$$\begin{array}{c}\hfill H_0=89.6{\left(\frac{\lambda }{{10}^4}\right)}^{3/2}\left(\frac{z}{0.001}\right){\left(\frac{\tau }{\mathrm{day}}\right)}^{-1}{\left(\frac{f_{\nu }/cosi}{\mathrm{Jy}}\right)}^{1/2}{\left(\frac{\chi }{4}\right)}^{4/3}\left[\frac{\mathrm{km}}{\mathrm{s}\mathrm{Mpc}}\right],\end{array}$$(11)

where τ is the time delay measured at the wavelength λ, z is the source redshift, f_ν is the measured flux at the frequency ν corresponding to the wavelength λ, and χ is the parameter which is used to account for systematic discrepancies due to conversion of annulus temperature to the corresponding wavelength λ at a given radius. According to Wien’s law, χ = 4.97 (Netzer 2022); however, under the assumption of a flux-weighted radius, χ is typically taken to be 2.49 (Edelson et al. 2024). In our analysis, considering the AGN accretion disk SED as a multicolored blackbody spectrum, we adopted χ = 4.

We can test the model using this formula, but it is important to use the disk delay rather than the total combined disk plus BLR time delay. The disk delay from model C is τ₅₁₀₀ = 1.008 days at λ = 5100 Å relative to the X-ray source. For a coronal height of h = 0, this delay becomes τ₅₁₀₀ = 0.898 days (i.e., 1.008 − 0.110 days), and this value is therefore used in Equation (11). We convert the flux from Lu et al. (2022), F_λ = 5.22 × 10⁻¹⁵ erg s⁻¹ cm⁻² Å⁻¹, to f_ν = 4.590 × 10⁻²⁶ erg s⁻¹ cm⁻² Hz⁻¹. Assuming an inclination angle of i = 40° for the source, we derive H₀ of 47.54 km s⁻¹ Mpc⁻¹. This represents a considerable improvement over the value of H₀ = 15 km s⁻¹ Mpc⁻¹ obtained by Cackett et al. (2007) for the same source. The difference is largely due to the fact that we are using only the disk delay. If we had used the total time delay (disk + BLR) of 2.00 days, we would have obtained a much lower value of H₀ = 21.34 km s⁻¹ Mpc⁻¹, closer to the result from Cackett et al. (2007).

However, adopting χ = 2.49 yields a derived H₀ of 25.27 km s⁻¹ Mpc⁻¹ for the disk delay, or 11.34 km s⁻¹ Mpc⁻¹ when including both the disk+BLR delay. Since there is no universally accepted value of χ, Equation (11) inherently carries additional uncertainty associated with this parameter.

5.2. Direct method to measure the distance to the source and H₀ from combined spectral and time delay fitting

Model C in the version described in Section 4.4.3 is based on the known distance which in turn constraints the accretion rate and the bolometric luminosity of the source. Therefore, if we indeed want to measure the distance to the source we cannot rely on quantities that were determined knowing its distance. We thus must allow the Eddington rate to cover a broad range, reflecting the possibly broad range of the luminosity distance of the source. However, we keep M_BH fixed, as it is derived from spectral fitting – using the FWHM/line dispersion (σ) of the Hβ line and the time delay measured from RM, making it independent of the luminosity distance. We also fix the viewing angle.

In this section, we present a series of computations using a grid of pre-selected luminosity distances. For each luminosity distance, we fit several model parameters by minimizing the χ² value. The parameters we vary when fitting the spectrum include the Eddington ratio (ṁ), the cold disk inner radius (r_in), the height of the corona (h), the optical depth of the warm corona (τ_{op − depth}), the temperature of the hot corona (T_e), the covering factor of the BLR (f_BLR), and the normalization of the starlight (A_star).

The incident radiation in the CLOUDY computations is determined by the choice of the luminosity distance, as the incident radiation must match the observed hard X-ray flux. Therefore, for each value of the luminosity distance, we use a different spectral shape for the BLR component. Figure 9 illustrates exemplary shapes corresponding to different normalizations of the incident flux. It is important to note that the distance to the BLR and the internal parameters of the clouds remain unchanged throughout these computations.

Fig. 9. Examples of the BLR emission profile for three values of the logarithm of the bolometric luminosity log L [erg s−1] = 44.0 (black line), 43.8 (blue line), and 44.2 (red line) corresponding to the three values of the luminosity distance, DL = 70, 56, and 84 Mpc, respectively.

The dependence of the BLR shape on the incident flux is not monotonic. However, the changes are relatively small, as the range of the considered flux amplitude variation is limited. It is important to emphasize that we do not modify the spectral shape of the incident radiation, as it was derived from observations rather than the model. During spectral fitting, we construct the final model spectrum based on absolute luminosity and compare it to the observed flux, taking into account the adopted distance as

$$\begin{array}{c}\hfill F_{\lambda }=\frac{L_{\lambda }}{4\pi D_{\mathrm{L}}^2}·\end{array}$$(12)

The value of χ² is calculated at several wavelengths, excluding the intense emission lines. Soft X-ray data points, taken from Mehdipour et al. (2015), are also included in the calculation. We assume a 10% error in the spectral fitting process. For the time delay, the χ² value is based on the errors complementing the time delay measurements (Fausnaugh et al. 2016).

The results of the combined delay and spectra fitting is shown in the top panel of Figure 10 in luminosity distance space. Now, in the χ² − D_L plane, each point in the sequence corresponds to the same Eddington ratio and other parameters in both the time delay and spectral fitting. The minimum occurs at D_L = 74 Mpc in our grid, though the minimum is relatively shallow. The best fits for both the spectrum and the time delay at this grid are also presented in Figure 8.

Fig. 10. From top to bottom: (i) Best-fit results for the simultaneous fitting of spectrum and the time delay as a function of the luminosity distance, marginalized over other model parameters; (ii) the same plot but as a function of the resulting Hubble constant; (iii) best-fit results for the spectrum alone as a function of the luminosity distance, marginalized over other parameters; and (iv) the corresponding plot for the time-delay fit alone.

The accretion rate predictably increases along the sequence, from 0.008 to 0.038, while the covering factor remain roughly the same, at ∼30%. The results are weakly sensitive to the lamp height, which remains close to ∼25 r_g, as well as to the inner radius of the cold disk, which is located at ∼90 r_g in all fits.

The value of the luminosity distance can be directly translated into the Hubble constant using the expression H₀ = zc/D_L, given that the redshift to the source is small. Thus, the same result, but as a function of the implied H₀, is shown in the second panel from the top in Figure 10. The distance of D_L = 74 Mpc corresponds to a Hubble constant of 66.8 km s⁻¹ Mpc⁻¹. Interpolating between the grid points we find the minimum location at 73.96 Mpc. Since the minimum is shallow, the error is large, leading to our constraint of D_L as $D_{\mathrm{L}}=74.0_{-10.0}^{+2.4}$ Mpc. This translates to Hubble constant constraints: $H_0=66.9_{-2.1}^{+10.5}$ km s⁻¹ Mpc⁻¹.

The error estimate presented here is necessarily approximate, as it is based on a qualitative assessment of the effective number of independent degrees of freedom (d.o.f.) in the model. Our dataset includes 17 measurements of time delays and we sample the spectrum in the same number of wavelengths. The model involves seven free parameters, in addition to four parameters that are kept fixed, namely: M_BH, viewing angle, local cloud density, and column density. Each model realization was computed for a fixed luminosity distance. This setup yielded 23 d.o.f. We simply determined the 1 sigma error of the reduced χ² as χ²/d.o.f. ≤ χ_min²/d.o.f. + 1/23. We did not perform tests to check whether all these measurements were statistically independent.

A more robust approach to error estimation would be to use a bootstrap resampling method, which could account for eventual correlations. Unfortunately, our current implementation is not optimized for such computationally intensive analysis, which would require at least 100 realizations to obtain reliable error estimates. Code optimization is underway, as we view this methodology as promising for future studies.

We note that in this analysis, we did not vary all the parameters involved in model C when calculating the spectral and time delay fits, as described above. Therefore, the derived constraint should be interpreted as an indication that the methodology is working, with further progress expected to reduce the error.

The value is well within the error range of the expected values (Abdalla et al. 2022; Freedman et al. 2025). The measurement error may decrease with improved modeling of the continuum in this source. Additionally, combining measurements from more sources will help further reduce the error. This value is higher than the 47.54 km s⁻¹ Mpc⁻¹ obtained from the analytical formula, which suggests that the proper methodology – avoiding forced power-law solutions and the extra χ parameter – is essential. To illustrate the necessity of combining delay and spectral fitting, we also present the results obtained when only one type of data (either delays or spectra) is used independently.

In experiments where only the time delay is fitted, we vary four parameters (ṁ, r_in, h, f_BLR), while the other three (T_e, τ_{op − depth}, A_star) are fixed at the values determined by the spectral best-fit, as these parameters have a negligible impact on the time delay. This strategy is motivated by the fact that spectral fitting is computationally efficient, whereas time delay fitting is considerably more time-consuming. Moreover, when fitting the SED, all parameters must be optimized simultaneously.

The χ² values in the luminosity distance space for the SED-only fit are shown in the second panel from the bottom of Figure 10. The key parameter that changes with D_L is the accretion rate. The formal minimum of χ² occurs at D_L = 58 Mpc, with the Eddington ratio varying from ∼0.008 at D_L = 50 Mpc to ∼0.034 at D_L = 102 Mpc. Along with this change in the accretion rate, the coronal height increases from 16.5 r_g to 53.0 r_g. The covering factor of the BLR, f_BLR, remains roughly the same, ∼0.30. The stellar contribution also decreases by a factor of 4 as D_L increases. Despite these variations, the optical depth of the hot corona remains roughly unchanged at 20, while the electron temperature increases slightly with D_L, from 6.6 × 10⁶ K to 8.0 × 10⁶ K. The reduced χ² is of order of 1, suggesting that the adopted spectral errors of 7% are realistic.

The χ² values for the case where we fit only the delay data set in the lag-spectrum is shown in the bottom panel of Figure 10. The dependence of χ² on D_L is relatively weak, exhibiting an almost monotonic decrease in χ² values with increasing luminosity distance. This sequence exhibits very different properties compared to the spectral sequence. The Eddington ratio remains constant at 0.03, independent of D_L, since there is no direct coupling between the delay, the distance, and the source luminosity. The lamp height is monotonically increasing with D_L. Additionally, the weak variation in the irradiation flux does not lead to a substantial change in the model structure.

Therefore, the combined fitting of both the SED and the time lag-spectrum is essential. Spectral fitting enforces a change in the best-fit Eddington ratio with D_L, which, in turn, affects the delay fitting. A small Eddington ratio implies a smaller source size and shorter time delays, as expected from standard disk.

6. Discussion

The aim of this paper is two-fold: first, to test the FRADO model for the formation of the BLR and second, to attempt the estimation of the Hubble constant from such a self-consistent model by determining the redshift-independent distance to the source NGC 5548. While this is solely a pilot study, our results demonstrate that the method holds significant potential.

The model itself involves very few parameters. Here, M_BH has been fixed and we have not yet tested its dependence on the results. However, it is worth noting that M_BH can be considered independent from cosmology since it is based on line kinematic width and the measured time delay. The key parameters, namely, the Eddington ratio, BLR contribution, and cold disk inner radius, can be fitted using the observed mean spectrum and the continuum time delays measured in the continuum–RM campaign reported in Fausnaugh et al. (2016). The best fit favors a redshift-independent distance of ∼74 Mpc to the source. This distance, when converted to the Hubble constant, implies a value of $H_0=66.9_{-2.1}^{+10.5}$ km s⁻¹ Mpc⁻¹. Additionally, the model successfully reproduces the Hβ time delay. One hidden parameter in the FRADO model is the dust sublimation temperature, T_sub, which we adopt as 1500 K. Variations in this temperature would directly influence the favored distance, roughly following the relationship r ∝ (T_sub/1500 K)^−4/3.

The multiwavelength time delays in NGC 5548 have also been modeled in other studies. The original work by Fausnaugh et al. (2016), which provided the delay measurements, employed only power-law parametric modeling. More recently, Kammoun et al. (2023) fit the NGC 5548 time delay using a new model that incorporates relativistic effects, such as the X-ray reflection component and albedo, while also exploring the influence of black hole spin and coronal height. Their SED fitting yielded a higher Eddington ratio (0.05) than our solution, along with a higher M_BH, resulting in a significantly greater bolometric luminosity for the source compared to our model. However, we cannot analyze this discrepancy in detail, as their best-fit SED was not explicitly presented. While Kammoun et al. (2023) strongly advocated for the inclusion of color correction, they did not account for the BLR contribution in their model. Notably, in both of these studies, the distance to NGC 5548 was determined using redshift measurements.

The recent study by Netzer et al. (2024) has several key elements in common with our work. However, their BLR model is not based on a dynamical framework (as ours is), requiring additional arbitrary parameters to define the dynamics and geometry. At the same time, their study addresses important factors, such as obscuration, which are not considered in our analysis. Despite successfully fitting both the spectra and time delays for their target, Mrk 817, they did not attempt to interpret their results in terms of a redshift-independent distance or derive an estimate for the Hubble constant.

Constructing a complete model of the BLR remains a significant challenge, with progress being made along several fronts. For instance, Long & Dexter (2025) analyzed the STORM data for NGC 5548 and, by combining a few representative geometries with table models that span a range of cloud ionization properties, concluded that none of the simplest models adequately reproduce the observed data. This highlights the limitations of overly simplified assumptions. Additionally, parametric approaches, such as those implemented in the BELMAG code (Rosborough et al. 2024, 2025), exemplify the complexity and richness required to capture the full diversity of BLR structures.

Our approach does not yet incorporate a detailed treatment of the ionization levels of individual clouds, as our objective is more modest and thus we do not aim to reproduce all features of the available data for the NGC 5548 source. Instead, we begin by establishing the dynamics based on the FRADO model, which alone enables us to roughly reproduce key observables such as the Hβ time delay, the transfer function, and the mean spectrum (apart from intensities of all lines), as well as the overall time-delay spectrum.

The geometry adopted in our model naturally leads to double-peaked emission lines. This is consistent with the observed double-peaked structure of the Hβ line, particularly in the rms profile (e.g., Lu et al. 2016; Shapovalova et al. 2009; Long & Dexter 2025). Moreover, other emission lines also exhibit double-peaked profiles (Horne et al. 2021) and the velocity-resolved time delays further support this structure (Long & Dexter 2025). Our model recovered transfer function similarly displays a double-peaked profile, consistent with an axially symmetric geometry dominated by Keplerian motion. In contrast, Horne et al. (2021) obtained an almost single-peaked transfer function for Hβ, which may indicate the presence of cloud emission anisotropy, an effect that was previously considered by Ferland et al. (1979), as well as possible selective shielding of parts of the BLR by outflowing material located within the region.

There remains considerable room for advancement in the coming years, particularly in refining the physical details of BLR modeling. Furthermore, the potential application of continuum time-delay measurements to cosmology provides an additional strong motivation for continued progress in this area. Taken together, these comparisons highlight the unique strengths of our approach. Unlike previous studies that either relied on power-law fitting, incorporated relativistic effects without addressing the BLR contribution, or did not pursue an independent distance determination, our method integrates spectral fitting and time delay modeling within a self-consistent framework. By leveraging the FRADO model, we obtained a redshift-independent distance to NGC 5548. While further refinements (e.g., improved continuum modeling and broader sample studies) are needed to reduce the uncertainties, our results demonstrate the potential of this methodology in offering an alternative approach to cosmic distance measurements.

6.1. Covering factor

Although we refer to our disk-plus-BLR model as self-consistent, not all of its components are fully self-consistent yet. In principle, for a given M_BH and Eddington ratio, the BLR covering factor should be determined by the model. Some tests of the predicted covering factor have been performed for relatively high Eddington ratios (Naddaf et al. 2023; Naddaf & Czerny 2024). However, at lower Eddington ratios, the clouds are located too close to the disk, leading to inconsistencies.

The aspect ratio of the cloud distribution, as shown in Figure 1, implies a covering factor of approximately 5%, whereas the value favored by time delay and spectral fitting is 36%. While the BLR location is well represented, the vertical motion of the clouds appears insufficient when considering only radiation pressure from dust. This suggests that line driving likely plays a role, but it is not yet included in the FRADO model.

Incorporating line driving, such as by combining FRADO with the QWIND model of Risaliti & Elvis (2010), is highly complex. The improved version of the QWIND software (Quera-Bofarull et al. 2023) is publicly available, providing a strong foundation for future work. However, it still lacks dust physics, which is crucial for determining the radial structure of the BLR at low luminosities. Future developments integrating these elements will be essential for achieving a fully self-consistent model.

6.2. Relativistic effects and X-ray albedo

Currently our model does not account for the effects of GR, which are critically important for phenomena within the inner 10 r_g. In these inner regions, spacetime curvature significantly influences both the propagation of photons from the lamppost (a point-like X-ray source) to the accretion disk and their subsequent travel from the disk to the observer. Additionally, the energy dissipation profile in the cold disk depends on the black hole spin (see e.g., Novikov et al. 1973).

However, the lamppost geometry itself is a considerable simplification, and in the case of NGC 5548, the cold accretion disk does not extend all the way to the ISCO. As a result, the reprocessing primarily occurs at larger radii, where GR effects are less dominant.

Another important consideration is the disk albedo, which is a complex function of the density in the disk atmosphere and the ionization parameter, both of which are likely to vary with radius. This issue is best addressed in systems that exhibit broad, relativistically smeared X-ray Fe Kα lines, where model assumptions can be more directly tested (e.g., Ballantyne et al. 2001; Dovčiak et al. 2004; García et al. 2013). Even models of reprocessing by the warm corona were studied recently by Petrucci et al. (2020).

For the case of continuum time delays, GR effects have been specifically studied for NGC 5548 by Kammoun et al. (2024). However, these models do not include contributions from the BLR and fail to reproduce features such as the Balmer and Paschen edges in the time-delay spectrum.

In the future, it may become possible to integrate these various components into a single model. However, full GR-based computations are computationally expensive and require the inclusion of additional free functions (radial profile of the ionization parameter, disk geometrical height, hard X-ray source geometry) to model the reprocessing accurately. Overall, as demonstrated by Kammoun et al. (2021b), the disk response is indeed sensitive to both relativistic effects and albedo. Nonetheless, those models do not incorporate an inner warm corona and therefore cannot be used for direct comparison in our case.

6.3. Dust sublimation temperature

In this work, we adopt a universal dust sublimation temperature of 1500 K, which plays a crucial role in the FRADO model as it sets the onset of BLR cloud launching. This choice directly influences the analytical expression for the Hubble constant derived in Equation (11). Given the temperature profile in the standard accretion disk that underlies this formula, our inferred values for the Hubble constant scale approximately as:

$$\begin{array}{c}\hfill H_0\propto {\left(\frac{T_{\mathrm{sub}}}{1500\mathrm{K}}\right)}^{4/3}.\end{array}$$(13)

However, the dust sublimation temperature is not firmly established. While Barvainis (1987) originally recommended 1500 K as a representative value, more recent papers recommend 1500 K for standard mixture of silicate/graphite dust (e.g., Mor & Netzer 2012), but for pure graphite, higher values such as 1800 K are proposed (e.g., Mor & Netzer 2012). A broader, density-dependent range was presented by Baskin & Laor (2018), whose Figure 1 shows that at densities around 10¹¹ cm⁻³, sublimation temperatures for silicate and graphite grains could reach 1600 K and 2000 K, respectively.

Nonetheless, graphite is unlikely to be present in the nuclear parts of AGN, as the 2175 Å feature is generally not seen in absorption (Laor & Draine 1993; Gaskell et al. 2004), which rather suggests the presence of amorphous carbon grains (e.g., Czerny et al. 2004). However, there are no reliable estimates for the sublimation temperature of amorphous carbon in AGN environments. Stellar measurements from AGB stars give a temperature of 1170 ± 60 K, very close to the sublimation temperature of silicate in these surroundings.

An alternative constraint was proposed by Vazquez et al. (2015), who argued for dust temperatures as high as 1800 K based on dust time-delay measurements rather than direct spectral analysis. However, their interpretation may be affected by an incomplete subtraction of the accretion disk contribution to the observed dust-delays (Mandal et al. 2024).

Given these uncertainties, we adopt 1500 K as a conservative and widely used value, but we acknowledge a plausible error margin of ±100 K. Based on the scaling relation above, this translates into a ∼9% uncertainty in the inferred Hubble constant, representing an important systematic contribution to our results.

The presence of dust, particularly of dust with atypical properties, may additionally imply problems with establishing the intrinsic SED. For example, Gaskell et al. (2023) argue that, based on their extinction curve, the UV flux of NGC 5548 may be underestimated by a factor of seven. Tests of such issues can be performed in the future.

6.4. Other fixed parameters as sources of systematic errors

It would be computationally time-consuming to move the other fixed parameters to the parameters fitted to the data, and to study the degeneracy between all these parameters. However, we can easily discuss their expected importance looking at the analytical formula for the Hubble constant of Collier et al. (1999) (see Equation (11)). Those extra fixed parameters are: M_BH, inclination angle, gas density, and column density. The black hole mass does not enter into Equation (11). The monochromatic flux, f_ν, is not sensitive on M_BH as long as the optical spectrum does not show strong departure from a power law (Shakura & Sunyaev 1973). This happens only in the far UV. The presence of the warm corona complicates the reasoning. Also the shape of the structure function shows some dependence on M_BH (see e.g., Naddaf et al. 2025, for most recent study), and not just on local monochromatic flux which is actually fitted. The onset of the BLR does not depend on M_BH, as it is set by the dust sublimation temperature, but the vertical velocities increase with M_BH. Overall, it should not lead to significant effect. The viewing angle is important, as it is present in Equation (11). For the viewing angle of 30°, the Hubble constant would be lower by 6%, and for 50° it would be higher by 8%. Therefore, this effect is important for any specific source. On the other hand, if the method is later applied to a large sample of objects at different redshifts, the viewing angle effects will average out, as there is no strong evidence for a redshift-dependent trend in the viewing angle or the dusty torus opening angle (e.g. Prince et al. 2022; Rałowski et al. 2024). However, assessing the influence of the local density and column density of the BLR clouds is considerably more challenging, as both parameters impact the ionization state and the effective optical depth of the region. These quantities can be effectively constrained by analyzing the shape of the delay enhancement near the Balmer and Paschen edges. Moreover, any reduction in effective optical depth can be compensated for by adjusting the covering factor, which is treated as a free parameter in the model. We certainly plan to refine and improve the model in future work.

6.5. FRADO versus LOC models as a method to reproduce the BLR transfer function

FRADO is a dynamical model which aims to reproduce the physical origin of the clouds, their location and their velocities. The model is not yet fully matured, but there are clear prospects for the model further development. For the dynamics, the addition of the line driving, mentioned above, should be included as its importance is increasing with the height of the clouds. Next, the emissivity of each clouds should be calculated individually, according to its location. Finally, there are two even more complex problems. One is a good realistic description of the shielding of clouds by themselves, important at low heights from the disk. The second one is determination of the density distribution of material based on thermal instability of the medium, as described already by Krolik et al. (1981). The wind is launched as a continuous medium, with clouds forming out of this wind in short thermal timescale.

Locally optimally emitting cloud (LOC) models, on the other hand, by their concept (power law distribution of clouds in radius and in density, with the limiting densities and limiting radii, and power law slopes as parameters) are better suited for accurately reproducing the line ratios observed in AGN spectra. The additional parameters are still required if the distribution is flat, concentrated around the accretion disk, instead of a spherical distribution as in the original model of Baldwin et al. (1995). This model is mature, but its parameters must be controlled by the data. However, to derive the BLR response, additional assumptions are required. Lawther et al. (2018) combined a spherically symmetric approach with the continuity equation to constrain the radial distribution of clouds, and were thereby able to predict the expected time delays for all emission lines. Additionally, using LOC model, Korista & Goad (2019) pioneered the study of the disk continuum contamination by BLR continuum. However, it remains to be seen whether LOC models can be applied in the same way as FRADO to simultaneously fit both the spectra and time delays, which would then enable a determination of the Hubble constant. This goal has not been achieved so far.

6.6. Other distance measurements to the source

The heliocentric redshift of NGC 5548 (z = 0.017175 ± 2.30 × 10⁻⁵; de Vaucouleurs et al. 1991) corresponds to a distance of 79.05 ± 5.54 Mpc. It is different when the motion of the source is measured relative to the Cosmic Microwave Background (CMB)². The proper motion of the source is thus not negligible, as indicated by the difference between its heliocentric and CMB-based velocities reported by NED (5149 ± 7 km s⁻¹ vs. 5359 ± 16 km s⁻¹).

In addition to redshift-based estimates, several independent distance measurements exist for this source. Two estimates rely on the Tully-Fisher method: an older measurement of 103 Mpc (Bottinelli et al. 1984) and a more recent one that closely aligns with the redshift-based value (82.2 ± 16.4 Mpc; Robinson et al. 2021). Meanwhile, the AGN continuum time lag method yielded a significantly higher distance of 341 ± 62 Mpc (Cackett et al. 2007). In contrast, dust reverberation mapping provided a distance of 92.5 ± 1.5 Mpc (Yoshii et al. 2014), which, despite its small uncertainty, is not consistent with the redshift-based value.

Applying a velocity correction from the heliocentric frame to the CMB frame shifts the redshift from 0.017175 to 0.017876. This adjustment would, in turn, modify our estimated value of the Hubble constant to $69.6_{-2.2}^{+11.0}$ km s⁻¹ Mpc⁻¹. However, given the relatively large uncertainties involved, achieving a highly precise correction for the proper motion is not crucial in this context.

6.7. The Hubble tension

The recent measurements of the Hubble constant generally highlight a persistent tension between early- and late-Universe determinations (see Abdalla et al. 2022, for a review). A key issue in this debate is the determination of measurement uncertainties. For example, studies by Riess et al. (2016) and Planck Collaboration XIII (2016) suggest a significant discrepancy, while the baryon acoustic oscillation (BAO) measurements by Wang et al. (2017) are consistent with both within their reported errors. However, the overall tension remains unresolved, as indicated by the most recent results from Riess et al. (2024), Uddin et al. (2024), and Planck Collaboration VI (2020).

In response to this discrepancy, several independent methods have been proposed to confirm or address the tension (e.g., see de Jaeger et al. 2022; Wong et al. 2024; Trefoloni et al. 2024; Liu et al. 2024; Bargiacchi et al. 2025, for some of the latest examples). Notably, some local measurements yield lower values of H₀ (e.g., Grillo et al. 2024). The latest findings from the clustering of galaxies, quasars, and Lyman-α forest tracers, based on the first year of observations from the Dark Energy Spectroscopic Instrument (DESI Data Release 1; DESI Collaboration 2025), support a lower Hubble constant of 68.40 ± 0.27 km s⁻¹ Mpc⁻¹ for a flat ΛCDM model, with remarkably small uncertainty. Their results, incorporating constraints from Big Bang nucleosynthesis (BBN) and acoustic angular scale measurements from the CMB, suggest H₀ = 68.52 ± 0.62 km s⁻¹ Mpc⁻¹, while a joint analysis of CMB anisotropies and CMB lensing from Planck+ACT yields H₀ = 67.97 ± 0.38 km s⁻¹ Mpc⁻¹. Most recent measurements from the Dark Energy Survey (DES) experiment yields an H₀ value of $67.{81}_{-0.86}^{+0.96}$ km s⁻¹ Mpc⁻¹ (DES Collaboration 2025).

Although our method lacks the precision needed to decisively address this tension, it represents a significant improvement over the results obtained by Cackett et al. (2007) for this object. This progress is both substantial and encouraging.

Interestingly, distance-ladder-based measurements tend to yield higher values of H₀ compared to single-step methods, with Perivolaropoulos (2024) reporting an average of H₀ = 69.0 ± 0.48 km s⁻¹ Mpc⁻¹. Our target source, at a redshift of z = 0.017175, lies near the boundary of z ∼ 0.01, where some studies suggest a possible transition in cosmological behavior (e.g., Gavas et al. 2025; Verde et al. 2024). However, given the current uncertainties in our results, we cannot yet distinguish between competing interpretations. Future efforts should focus on reducing measurement errors and extending the analysis to AGN at both lower and higher redshifts to further explore this issue.

6.8. Prospects for continuum reverberation from LSST

The standard six-color photometry provided by the Legacy Survey of Space and Time (LSST), soon to be available from the Vera Rubin Observatory, may not achieve the same level of precision in time delay measurements as combined SWIFT/ground-based monitoring. However, there remain promising opportunities for detecting continuum time delays, particularly for AGNs with larger black hole masses with high luminosities (Pozo Nuñez et al. 2023, 2024; Panda et al. 2024).

For the most massive quasars, however, the contribution of the BLR to the observed time delay may be less significant. This is suggested by results for PG 2308+098, a quasar with M_BH of 10^{9.6 ± 0.1} M_⊙, as reported by Kokubo (2018). To better understand this effect, we can also investigate a potential scaling relation between BLR continuum contamination and emission line intensities, which may provide further insight into the role of the BLR in time delay measurements.

7. Conclusions

Our pilot study is aimed at testing the FRADO model for the BLR structure and estimate the Hubble constant based on the mean spectrum data and continuum time delay measurements for NGC 5548. The study has been highly successful, yielding several key outcomes, which are detailed below.

• The time delay of the Hβ line is approximately recovered, even though our model does not allow for direct adjustments of the inner and outer radii of the BLR. Instead, these parameters emerge naturally from the computations of radiation pressure acting on dust and the dust sublimation temperature.

• By using the redshift-based distance of NGC 5548 and its known luminosity, we were able to successfully model the overall spectrum, incorporating contributions from both the accretion disk and the BLR. Additionally, we accurately reproduce the time-delay pattern, emphasizing the crucial role played by the Balmer and Paschen edges.

• Relaxing the assumption of a known redshift-based luminosity, we modeled the spectrum and time-delay data across a range of Eddington ratios. For each value, we optimized the source distance, ultimately identifying a best-fit Eddington ratio and corresponding distance. This, in turn, provided an estimate of the Hubble constant.

• Although the uncertainty in the Hubble constant remains relatively large ($66.9_{-2.1}^{+10.5}$ km s⁻¹ Mpc⁻¹), our method demonstrates the potential for refinement. Further improvements, both for this source and for similar determinations using other AGNs with dense monitoring, could lead to significantly more precise measurements.

Our results validate the viability of using AGN continuum time delays and spectral modeling as an independent method for measuring cosmic distances and estimating the Hubble constant. While the current uncertainty is still substantial, the methodology holds promise for future applications. Expanding the sample to include additional AGNs with well-sampled continuum–RM data will be crucial for reducing uncertainties and improving the accuracy of this approach.

Acknowledgments

We thank the anonymous referee for valuable comments and suggestions. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (grant agreement No. [951549]). VKJ acknowledges the OPUS-LAP/GA CR-LA bilateral project (2021/43/I/ST9/01352/OPUS 22 and GF23-04053L). MHN also acknowledges the financial support by the University of Liege under Special Funds for Research, IPD-STEMA Program. SP is supported by the international Gemini Observatory, a program of NSF NOIRLab, which is managed by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the U.S. National Science Foundation, on behalf of the Gemini partnership of Argentina, Brazil, Canada, Chile, the Republic of Korea, and the United States of America. BC acknowledges the support from COST Action CA21136 – Addressing observational tensions in cosmology with systematics and fundamental physics (CosmoVerse), supported by COST (European Cooperation in Science and Technology). FPN gratefully acknowledges the generous and invaluable support of the Klaus Tschira Foundation. We thank Dr. Missagh Mehdipour for sharing the broadband SED data points with us.

References

Appendix A: H₀ based on observationally determined response function

We employed an alternative method to validate the results obtained from the FRADO model. In rare cases, NGC 5548 being one notable example, the transfer function for Hβ was directly determined from the data using densely sampled spectroscopic RM. For this purpose, we use the results presented by Horne et al. (2021). These observations were conducted in 2014, roughly during the same period with the photometric time-delay measurements. The study provides response functions for various emission lines, but we focus on those corresponding to the Hβ line, as they are most representative of the regions in the BLR where the Balmer and Paschen continua originate and play a dominant role in shaping the observed time delays. Specifically, we consider two versions of the Hβ response functions presented by Horne et al. (2021): one that encompasses the full Hβ line and another centered narrowly around 4863 Å. A comparison between our model-derived BLR response function and the two response functions from Horne et al. (2021) is shown in the upper panel of Figure A.1.

Both observed response functions lack a double-peak structure, which, however, is present in the case of the CIV line (see Figure 4 of Horne et al. 2021). More importantly, these response functions exhibit a significantly longer tail, with this extended tail being notably stronger in the full-Hβ response profile than in the narrower 4863 Å case. This long tail is also clearly visible in the 4830 Å bin shown in Figure 9 of Horne et al. (2021), suggesting that the extended delay may arise from contamination by Fe II emission. On average, the optical and UV Fe II delays tend to be longer than those of typical low-ionization lines such as Mg II and Hβ (Zajaček et al. 2024).

Having an observationally derived transfer function significantly reduces our dependence on the FRADO model summarized in Section 3.3. However, we still require CLOUDY computations to determine the wavelength-dependent efficiency of re-emission by the BLR, as well as a disk/corona model.

Therefore, we performed a fitting of both the mean spectrum and the wavelength-dependent time delay as before, for a set of values of the luminosity distance. For each luminosity distance D_L, the model parameters were refitted accordingly.

In the case of the entire Hβ transfer function, we were unable to tightly constrain the luminosity distance (blue data points in the bottom panel of A.1); the monotonic trend continues up to 98 Mpc, with marginal rise at 102 Mpc. We did not perform computations beyond this distance. However, for the transfer function centered around 4863 Å, a broad minimum is observed from 82 Mpc till 94 Mpc (red data points), which is shifted compared to the 74 Mpc value obtained using the theoretical FRADO model for the transfer function (black data points). Interestingly, the overall fit is the best for the FRADO model. The resulting constraint on the Hubble constant H₀ from the 4863 Å bin is $60.3_{-9.1}^{+2.0}$ km s⁻¹ Mpc⁻¹, generally lower than the FRADO-derived range but only marginally outside the 1 sigma error value due to relatively large uncertainty, further supporting the applicability of the FRADO model in such analyses.

Fig. A.1. Upper panel: Response function from the FRADO model (black), total Hβ response (blue) presented in Figure 4 of Horne et al. (2021), and the response at 4863 Å (red) shown in Figure 9 of the same paper. Lower panel: χ2 values as a function of luminosity distance for the three response functions are shown in their respective colors.

All Tables

All Figures

	Fig. 1. Upper panel: 3D plot of cloud positions from FRADO model for NGC 5548 with MBH = 5 × 107 M⊙, L/LEdd = 0.02, metallicity Z = 5 in solar units. The axes are in units of rg. Bottom panel: Cross-section of cloud positions for y > −300 rg and y < 300 rg. The black line indicates the thickness of the Keplerian disk. Clouds form a geometrically thin complex layer above the disk.
In the text

Fig. 2. Upper panel: Velocity–delay map of NGC 5548 constructed from FRADO model with inclination angle i = 40°. Color coding represents the number density of BLR clouds in each velocity-delay bin. The white dashed line shows the virial envelope corresponding to Keplerian disk-like rotation for v2 × τ = constant with MBH = 5 × 107 M⊙. Bottom panel: BLR response function generated using the cloud positions from the model shown in Figure 1, assuming an inclination angle of i = 40°.

In the text

	Fig. 3. Line profile generated by the FRADO model using the parameters from Figure 1.
In the text

	Fig. 4. Emissivity profile (ratio of the reprocessed to the incident continuum) of a BLR cloud for the adopted parameters: log nH[cm−3] = 11, log L [erg s−1] = 44, and the BLR distance of 1016 cm.
In the text

Fig. 5. Upper panel: Response functions of the disk at different wavelengths. Bottom panel: Combined response functions of the disk and BLR across the same wavelengths. These representative response functions are obtained for model C. Parameters used: MBH = 5.0 × 107 M⊙, Eddington ratio = 0.015, LX = 9.68 × 1043 erg s−1, height: h = 48.29 rg, Rin = 94.87 rg, Rout = 10 000 rg, and viewing angle: i = 40°.

In the text

Fig. 6. Model A. Upper panel: Observed delay from Fausnaugh et al. (2016) is represented in black. The delay calculated from the disk response function alone is shown in red, while the delay calculated from the combined response function of the disk and BLR is shown in brown. Bottom panel: Observed SED from Mehdipour et al. (2015) is shown in blue. The disk SED is represented in red, while green illustrates the disk SED combined with the emission lines. The starlight contribution is shown in black color. The final SED, incorporating contributions from the disk, star, and emission lines, is shown in purple. The selected points used for estimating the χ2 value are represented by brown squares. Parameters: MBH = 5.0 × 107 M⊙, Eddington ratio = 0.02, LX = 1.26 × 1044 erg s−1, height: h = 20 rg, Rin = 35 rg, Rout = 10 000 rg, color correction = 2.4, and viewing angle: i = 40°.

In the text

Fig. 7. Model B. Upper panel: Observed delay from Fausnaugh et al. (2016) is represented in black. The delay calculated from the disk response function alone is shown in red, while the delay calculated from the combined response function of the disk and BLR is shown in brown. Bottom panel: Observed SED from Mehdipour et al. (2015) is shown in blue. The disk SED is represented in red, while green illustrates the disk SED combined with the emission lines. The starlight contribution is represented in black, while the warm corona is shown in orange. The final SED, incorporating contributions from the disk, star, warm corona, and emission lines, is shown in purple. The selected points used for estimating the χ2 value are represented by brown squares. Parameters: MBH = 5.5 × 107 M⊙, Eddington ratio = 0.02, LX = 1.247 × 1044 erg s−1, height: h = 43 rg, Rin = 151 rg, Rout = 282 rg, color correction = 1.0, and viewing angle: i = 45°.

In the text

Fig. 8. Model C. Upper panel: Observed delay from Fausnaugh et al. (2016) is represented in black. The delay calculated from the disk response function alone is shown in red, while the delay calculated from the combined response function of the disk and BLR is shown in brown. Middle panel: Observed SED from Mehdipour et al. (2015) is shown in blue. The disk SED is represented in red, while green illustrates the disk SED combined with the emission lines. The starlight contribution is represented in black, while the warm corona is shown in orange. The final SED, incorporating contributions from the disk, star, warm corona, and emission lines, is shown in purple. The selected points used for estimating the χ2 value are represented by brown squares. Parameters used: MBH = 5.0 × 107 M⊙, Eddington ratio = 0.015, LX = 9.68 × 1043 erg s−1, height: h = 48.29 rg, Rin = 94.87 rg, Rout = 10 000 rg, and viewing angle: i = 40°. Bottom panel: X-ray-UV-optical SED fitting. The warm corona is represented by the magenta line, the hard X-ray power law is shown in dashed gold, the green line represents the sum of the two, and the black line displays the optical data.

In the text

	Fig. 9. Examples of the BLR emission profile for three values of the logarithm of the bolometric luminosity log L [erg s−1] = 44.0 (black line), 43.8 (blue line), and 44.2 (red line) corresponding to the three values of the luminosity distance, DL = 70, 56, and 84 Mpc, respectively.
In the text

Fig. 10. From top to bottom: (i) Best-fit results for the simultaneous fitting of spectrum and the time delay as a function of the luminosity distance, marginalized over other model parameters; (ii) the same plot but as a function of the resulting Hubble constant; (iii) best-fit results for the spectrum alone as a function of the luminosity distance, marginalized over other parameters; and (iv) the corresponding plot for the time-delay fit alone.

In the text

	Fig. A.1. Upper panel: Response function from the FRADO model (black), total Hβ response (blue) presented in Figure 4 of Horne et al. (2021), and the response at 4863 Å (red) shown in Figure 9 of the same paper. Lower panel: χ2 values as a function of luminosity distance for the three response functions are shown in their respective colors.
In the text