© The Authors 2024

1. Introduction

Many earlier works have shown that supermassive black holes (SMBHs) reside in the centre of most galaxies (e.g. Soltan 1982; Rees 1984; Ferrarese & Ford 2005; Davis 2004; Padovani et al. 2017). The M_BH-σ_* relation (the relationship between the central BH mass and the stellar velocity dispersion of the host galaxy) provides strong evidence for SMBH-host galaxy coevolution (e.g. Magorrian et al. 1998; Gebhardt et al. 2000; Treu et al. 2004; Kormendy & Ho 2013; Caglar et al. 2020). This requires robust measurements of BH masses. There are many ways of doing this, such as spatially resolving stellar kinematics (e.g. Gebhardt et al. 2000; Sharma et al. 2014), measuring the kinematics of (ionised) gas (e.g. Davies et al. 2004a,b; Hicks & Malkan 2008; Davis et al. 2013), utilising megamaser kinematics with very long baseline interferometry (VLBI; e.g. Kuo et al. 2020; Wagner 2013; van den Bosch et al. 2016), and reverberation mapping (RM; e.g. Peterson 1993; Lira et al. 2018; Cackett et al. 2021). Other indirect methods involve scaling relations, relying on single-epoch optical/near-infrared (NIR) spectroscopy to observe several line estimators and to measure their fluxes and widths (e.g. Hα, Hβ, Mg II λ2798, and C IV λ1549; e.g. Kaspi et al. 2000; Greene & Ho 2005; Shen & Liu 2012) as well as coronal lines to estimate the accretion disc temperature (e.g. Prieto et al. 2022).

When the SMBH is ‘active’, it accretes surrounding material and forms an accretion disc. The brightness of the resulting active galactic nucleus (AGN) dilutes the signal from the stars and gas surrounding it (in addition to the gas being subject to other forces than gravity, e.g. radiation pressure, winds, etc.), hindering BH mass measurements via stellar and gas kinematics (Gebhardt et al. 2000). The intense radiation can trigger megamasers and thus allow for this SMBH mass measurement method, but the need for a nearly edge-on line of sight (LOS) allows for only type 2 AGNs (van den Bosch et al. 2016) to be observed.

Reverberation mapping campaigns utilise the fluctuations over time of the continuum flux originating from the accretion disc (Blandford & McKee 1982; Peterson 1993; Peterson et al. 2004). These variations are then echoed or reverberated in the surrounding gas of the broad-line region (BLR). The time lag between the continuum and BLR variations can be measured from their light curves and utilised to estimate the radius of the BLR under the assumption that the measured time delay represents the light-crossing time between the accretion disc and BLR. The mass of the central BH is then calculated via the following formula:

$$\begin{array}{c}\hfill M_{\mathrm{BH}}=\frac{f{R}_{\mathrm{BLR}}{v}^2}{G},\end{array}$$(1)

where G is the gravitational constant, M_BH is the central BH mass, R_BLR is the BLR radius, and v is the full width at half maximum (FWHM) or dispersion (σ) of the line observed with RM. We note that f is the virial factor, which accounts for the unknown geometry and structure of the BLR.

Previous works have revealed that local AGNs exhibit a relationship between the BLR radius measured from AGN emission lines (e.g. Hα and Hβ) and AGN continuum luminosity (usually at 5100 Å), the so-called R-L relation (Kaspi et al. 2000; Bentz et al. 2013). The R-L relation, combined with a measurement of the gas velocity (e.g. the FWHM of the broad emission lines), can then be used to estimate BH masses of distant AGNs using Eq. (1) (e.g. Laor 1998; Vestergaard & Peterson 2006; Shen et al. 2023). However, BH mass estimations via single-epoch spectroscopy are only accurate if the R-L relation accurately represents the whole AGN population. Hence, confirming the relation with a more diverse and broad AGN sample is important. The accretion rate or L/L_edd has been suggested to serve as a ‘third’ parameter in the R-L relation, as highly accreting AGNs seem to have smaller BLR sizes than what is predicted by the canonical R-L relation (Du et al. 2014b, 2015, 2016; Du & Wang 2019). Most of these highly accreting sources are also highly luminous (log λL_λ(5100 Å)/erg s⁻¹ ≳ 44; Du et al. 2018), and hence the deviation from the canonical R-L relation is more prominent at higher luminosities. Du et al. (2018) and Du & Wang (2019) show that using ℛ(Fe II) (the equivalent width ratio of Fe II and Hβ) as a proxy of accretion rate can reduce the scatter of the R-L relation. Not accounting for the accretion rate could then lead to BH mass estimations via single-epoch spectroscopy and the R-L relation to be overestimated by as much as ∼0.3 dex (Alvarez et al. 2020; GRAVITY Collaboration 2023).

Reverberation mapping measurements have concluded that the BLR gas is mostly dominated by Keplerian motion and is therefore orbiting around the central BH (e.g. Gaskell 2000; Denney et al. 2009; Bentz et al. 2010, 2021b; Du et al. 2016; Grier et al. 2017; Williams et al. 2018; Villafaña et al. 2023). The BLR could also have considerable turbulence in the direction perpendicular to the BLR mid-plane, with significant inflow and outflow velocity components (Osterbrock 1978; Ulrich & Horne 1996). Several pieces of evidence of possible inflow and outflow signatures for some BLRs have been found (e.g., Kollatschny & Dietrich 1997; Denney et al. 2009; Bentz et al. 2009, 2010; Du et al. 2016). On the other hand, three BLRs whose kinematics were resolved by GRAVITY were found to be rotation dominated (GRAVITY Collaboration 2018, 2020a, 2021a). Hence, a single description of the BLR has not been established yet.

A crucial limitation of RM is the uncertain virial factor, f. Usually, f is calibrated assuming that AGNs follow a similar M_BH-σ_* relation with those of quiescent galaxies. Although f is different for each AGN and cannot be determined without velocity-resolved RM data, the average f can be calculated for a sample of AGNs. The average virial factor, ⟨f⟩, serves as a scaling factor of RM-derived BH masses to match the said relation. Various works have different values of ⟨f⟩ (e.g. Onken et al. 2004; Woo et al. 2010; Graham et al. 2011; Park et al. 2012; Batiste et al. 2017). Source variability and line width uncertainties could also lead to uncertainties in the measured time lag (see Vestergaard et al. 2011 and references therein). RM campaigns are also usually limited by the luminosity of the AGNs they can target, as more luminous targets should require much longer RM campaigns to measure the expected time lags. Recent RM surveys though are working to expand the high luminosity range (e.g. Woo et al. 2024).

Spatially resolving the BLR was initially difficult to accomplish due to its very small angular size (≲10⁻⁴ arcseconds; Blandford & McKee 1982). However, thanks to GRAVITY, the second-generation NIR beam combiner in the Very Large Telescope Interferometer (VLTI), the sensitivity of NIR interferometry has been exceptionally improved, and the light received by all four 8 m unit telescopes (UTs) has been combined to yield six simultaneous baselines (GRAVITY Collaboration 2017). This has allowed us to not just spatially resolve the BLR, but also obtain measurements of the SMBH mass and investigate the BLR structure at a high-angular resolution. To this end, we carried out an ESO Large Programme to observe the brightest type 1 AGNs, which span four orders of magnitude in luminosity, to spatially resolve their BLR and measure their central BH masses. This programme also aims to further understand the BLR structure and its intricacies, such as inflow and outflow motions that may be present in some systems, and to investigate GRAVITY-derived BLR radius and BH mass measurements and how they compare with scaling relations of local AGNs, such as the R-L and M_BH-σ_* relations. Three BLRs from this sample have already been spatially resolved by GRAVITY (GRAVITY Collaboration 2018, 2020a, 2021a,b). In parallel, several works have focussed on resolving the dust sublimation region (GRAVITY Collaboration 2020b, 2021a), and using the observed dust size to estimate BLR sizes (GRAVITY Collaboration 2020c; Abuter et al. 2024).

In this work, we present an analysis of four new GRAVITY observations of type 1 AGNs, namely Mrk 509, PDS 456, Mrk 1239, and IC 4329A, all in the local Universe (z < 0.2). For the entirety of our paper, we adopt a Lambda cold dark matter (ΛCDM) cosmology with Ω_m = 0.308, Ω_Λ = 0.692, and H₀ = 67.8 km s⁻¹ Mpc⁻¹ (Planck Collaboration XIII 2016). We discuss the targets’ properties in Sect. 2. We describe the observations, data reduction, and the resulting flux and differential phase spectra in Sect. 3. We report the model-independent photocentre positions of each target’s BLR in Sect. 4. We describe our BLR model in Sect. 5, and discuss the importance of differential phase data in our BLR modelling in Sect. 6. We present the results of our BLR model fitting in Sect. 7, while we place the derived BLR sizes and BH masses of our targets in the context of the R-L and M_BH-σ_* relations in Sects. 8 and 9, respectively. We also observed an offset between the BLR and continuum photocentres of our targets and discuss its possible origin in Sect. 10. We also explain how we calculated the virial factors of our targets in Sect. 11. Finally, we present our conclusions and future prospects in Sect. 12.

2. Targets

All four targets were observed with GRAVITY as part of a large programme (and initial pilot projects) that aims to spatially resolve the BLR and measure BH masses for a sample of the brightest (nuclear luminosity of K ≲ 10 mag and V ≲ 15 mag) type 1 (Sy1 or QSOs with broad Paα or Brγ emission lines) AGNs in the local Universe¹. Table 1 summarises the properties of our four targets and the three other targets that are already published.

2.1. Mrk 509

Mrk 509 is a type 1 Seyfert galaxy that shows significant outflow signatures in terms of mildly relativistic (∼0.14–0.2c) FeXXVI Kα and Kβ absorption features (Cappi et al. 2009) detected in its X-ray spectra. These imply a possible varying structure and geometry (Cappi et al. 2009). RM campaigns have estimated the BLR size of Mrk 509 to be R_BLR ∼ 80 ld (Carone et al. 1996; Peterson et al. 1998; Bentz et al. 2009; Shablovinskaya et al. 2023). Its central BH mass was also calculated to be M_BH = 10^7.9 − 10^8.3 M_⊙ (Peterson et al. 1998, 2004).

2.2. PDS 456

PDS 456 is known to be the most luminous radio-quiet AGN in the local Universe (z ≲ 0.3; Torres et al. 1997; Simpson et al. 1999; Bischetti et al. 2019; Yun et al. 2004). Aside from its brightness and proximity, it has been extensively studied due to its ultra-fast outflows (UFOs) of highly ionised Fe (Reeves et al. 2003). These outflows possess high kinetic power (∼20% of its bolometric luminosity; Luminari et al. 2018) and are also radiatively driven, which can be explained by the AGN’s high Eddington ratio (Matzeu et al. 2017). Its BH mass, unfortunately, is not well-studied. No RM measurements are available for PDS 456 because such a luminous target should produce very long time lags (∼several years), requiring decade-long RM campaigns. Therefore, BH mass estimations based on empirical relations from RM are used (Reeves et al. 2009; Nardini et al. 2015). Its central BH mass is estimated to be M_BH ∼ 10⁹ M_⊙. Its BLR size is also implied to be large; GRAVITY Collaboration (2020c) partially resolved the continuum hot dust emission region of PDS 456 and used the measured continuum size from the fringe tracker data and the differential visibility amplitude data from the science channel (SC) to indirectly infer its BLR size. However, their estimate is highly uncertain due to the weak correlation between the SC differential amplitude and the BLR size for very small angular sizes relative to the baseline resolution. GRAVITY Collaboration (2023) estimated its BH mass and BLR radius to be M_BH ∼ 10^8.68 M_⊙ and R_BLR ∼ 150 ld from the relation between the BLR and dust continuum size. It is interesting to note that GRAVITY Collaboration (2023) measured a smaller BH mass than previous works have estimated, and this deviation is suggested as a consequence of its high Eddington ratio (Du & Wang 2019) that has been implied by its strong FeII features and lack of [OIII] emission lines (e.g., Simpson et al. 1999).

2.3. Mrk 1239

Mrk 1239 is a narrow-line Seyfert 1 (NLSy1) galaxy well known for its high polarisation (Goodrich 1989; Smith et al. 2004), relatively redder optical-IR colour compared to other typical NLSy1s, and presence of a radio jet-like structure (Orienti & Prieto 2010; Doi et al. 2015) oriented perpendicular to the polarisation angle. Mrk 1239 has not been a target for RM campaigns, but previous works have estimated its BH mass and BLR size via indirect estimates from various scaling relations. They all point to the same conclusion that Mrk 1239 has a small but uncertain BH mass, M_BH ∼ 10^5.7 − 10^7.0 M_⊙ (e.g. Kaspi et al. 2005; Greene & Ho 2005; Ryan et al. 2007; Du et al. 2014a; Buhariwalla et al. 2020; Pan et al. 2021; Abuter et al. 2024). Its BLR size is also estimated to be R_BLR ∼ 10–20 ld (e.g. Du et al. 2014a; Abuter et al. 2024).

2.4. IC 4329A

IC 4329A is a type 1 Seyfert galaxy that is well-observed in X-ray (e.g., Madejski et al. 1995; Delvaille et al. 1978; Piro et al. 1990; Nandra & Pounds 1994). Several RM campaigns were implemented to measure the BLR size and BH mass of IC 4329A. For instance, Winge et al. (1996) presented the first RM data for IC 4329A, which was reanalysed by Peterson et al. (2004). However, the poor quality of the light curves raised caution about their estimated BH masses and BLR sizes. Wandel et al. (1999) calculated a BLR size and BH mass of $R_{\mathrm{BLR}}=1.4_{-2.9}^{+3.4}$ ld and M_BH ∼ 10^7.4 M_⊙ after re-analysing the spectrum taken by Winge et al. (1996). Bentz et al. (2023) presented the latest RM campaign for IC 4329A, and together with BLR model fitting, they acquired a BLR size of $R_{\mathrm{BLR}}=14.2_{-3.7}^{+7.2}$ ld and a BH mass of $M_{\mathrm{BH}}=4_{-2}^{+10}\times {10}^{7.0}{M}_{\odot }$. Abuter et al. (2024) also estimated the BLR size of IC 4329A via dust size measurements, yielding R_BLR ∼ 17.4 ld.

3. Observations and data reduction

3.1. Observations

The GRAVITY observations were done over several nights between August 2017 and July 2021. We used the single-field on-axis mode with combined polarisation for all observations. Each observation sequence follows that described in GRAVITY Collaboration (2020a). First, the telescopes pointed to the target and closed the adaptive optics (MACAO, Arsenault et al. 2003) loop. The light was then propagated to GRAVITY, where the fringe tracking (FT) and science channel (SC) fibres align on the target via internal beam tracking of GRAVITY. The exposures were acquired after the fringe tracker found the fringes and began tracking. The integration time for each exposure frame is DIT = 30 s, and the number of frames for each set of exposures is NDIT = 12. All observations were taken in MEDIUM resolution. The number of obtained exposures varied among objects. A calibrator star (usually an A- or B-type star) close to the target was also observed and used to calibrate the flux spectrum of the AGN. This ensures that atmospheric and vibrational effects, coherence loss, and birefringence are all accounted for. We adopted the same pipeline data reduction as GRAVITY Collaboration (2020a) to calculate the complex visibilities from the raw data, enabling us to extract quantities correlated to the physical properties of our targets. Table 2 shows the date, exposure time and weather conditions (seeing and coherence time) during observations.

We use all available GRAVITY data for three of the four targets. However, we only use the observations of Mrk 509 from 2021. This is because earlier (2017 and 2018) observations of Mrk 509 show differential phase errors that are ∼50% larger than those from 2021, which benefit from the factor two better throughput of the science channel spectrometer after the grating upgrade at the end of 2019 (Yazici et al. 2021).

3.2. Differential phase spectra

To produce the differential phase curve on each baseline for all the targets, we first selected exposures with fringe tracking ratio (percentage of utilised time that fringe tracking was working) greater than or equal to 80%, removing 16, 17, 8, and 33 exposures for Mrk 509, PDS 456, Mrk 1239, and IC 4329A, respectively. We followed the method from GRAVITY Collaboration (2020a) in removing instrumental features from the differential phase curves, estimating their uncertainties in each channel, and polynomial flattening. To check whether our uncertainties match the observed noise in the differential phase spectra and are a good representation of the actual dispersion of the differential phase spectra, we compare the average uncertainty in the whole spectra and the dispersion of the differential phase in regions beyond the wavelength range where the broad emission line is expected to be found. Our comparison suggests that only Mrk 1239 shows a significant discrepancy: the uncertainties in its differential phase spectra (which are ∼0.06–0.08° on average for all baselines) are lower compared to the dispersion of the differential phase values at around 2.18–2.20 μm and 2.22–2.24 μm (which are ∼0.08–0.14° on average). Hence, we adopted the dispersion of the differential phase values at these wavelength ranges as the representative uncertainty/error bar of all channels in each baseline. However, we emphasise that we also tried adopting the uncertainties in Mrk 1239’s differential phase spectra as its representative uncertainty, and we conclude that this does not greatly affect our results.

The definition of differential phase in the partially resolved limit that is appropriate for our AGN observations (e.g. GRAVITY Collaboration 2018, 2020a) is as follows:

$$\begin{array}{c}\hfill \mathrm{\Delta }{\varphi }_{\lambda }=-2\pi \frac{f_{\lambda }}{1+f_{\lambda }}\mathit{u}·{\mathit{x}}_{\mathrm{BLR},\lambda },\end{array}$$(2)

where Δϕ_λ is the differential phase measured in a certain wavelength channel λ, f_λ is the normalised flux in that channel, u is the uv coordinate of the baseline, and x_BLR, λ denotes the BLR photocentres. The differential phases are referenced to the continuum photocentre, which is placed at the origin. If the centre of the BLR coincides with the continuum photocentre, the differential phase signal reflects the kinematics of the BLR. Meanwhile, if the BLR to continuum offset is much larger than the BLR size, the x_BLR, λ becomes a constant approximately over different channels. The differential phase signal will be proportional to $\frac{f_{\lambda }}{1+f_{\lambda }}$ and resemble the profile of the broad line (see more discussion in Sect. 10). Hereafter, we call the phase signal generated by the global offset between the BLR and the continuum emission the “continuum phase”. Such a global offset is observed in IRAS 09149-6206 and NGC 3783 primarily due to the asymmetry of the hot dust emission (GRAVITY Collaboration 2020a,b). More details about the derivation of Eq. (2) are shown in Appendix B of GRAVITY Collaboration (2020a).

Figure 1 shows the differential phase spectra of our 4 targets averaged for each baseline. Most of the baselines of the targets show signals in their differential phase curve that are either negative (Mrk 509 and PDS 456) or positive (Mrk 1239 and IC 4329A) and mostly coincide with the peak of their line profiles. As for IRAS 09149-6206 (GRAVITY Collaboration 2020a), this indicates that their differential phase spectra show signatures of a continuum phase. The signals are also relatively small with absolute values ≲0.6°. We can, therefore, conclude that their BLRs are either intrinsically small or viewed at a very low inclination angle.

Fig. 1. Averaged differential phase spectra (coloured, based on their respective baselines labelled in the top row) and the normalised Brγ (for Mrk 509, Mrk 1239, and IC 4329A) and Paα (for PDS 456) spectra (black histogram) of Mrk 509 (first row), PDS 456 (second row), Mrk 1239 (third row), and IC 4329A (fourth row). The curves show negative signals in Mrk 509 and PDS 456 with peaks ranging from ∼ − 0.4° to −0.5°, and positive signals in Mrk 1239 and IC 4329A with peaks ranging from ∼0.5°–0.6°.

3.3. Normalised profiles of the broad Brγ and Paα emission lines

Modelling the dynamics and deriving the velocity gradient of a target’s BLR requires a well-measured nuclear emission line profile (Brγ for Mrk 509, Mrk 1239, and IC 4329A, and Paα for PDS 456). In general, the stacking and calibration of the line profiles are as follows: after fitting for and removing the Brγ absorption line of the calibrator star, each spectrum was normalised by fitting a third-degree polynomial to the continuum before stacking. We then followed the method of extracting the final line profile from GRAVITY Collaboration (2020a) wherein we considered the statistical uncertainty (i.e. the rms from each individual spectrum’s uncertainties) and the systematic uncertainty (i.e. mainly caused by variations in the calibrator data and sky absorption) of the spectra in producing the final flux error of the spectra. Lastly, the narrow components of the emission lines were also accounted for and removed following previous work (GRAVITY Collaboration 2020a). Removing the narrow components from the flux spectra of our targets is important because the flux spectra, together with the differential phase spectra, are fitted together with our BLR model (to be discussed in Sect. 5). Equation (2) assumes that all of the flux in the line profile is only partially resolved. Narrow line emission that enters the GRAVITY fibres occurs on large scales and is over-resolved and thus does not contribute to the differential phase. We note that the narrow component removal is only applied to the flux spectra and does not affect the interferometric data. We use the GRAVITY spectra of Mrk 509, PDS 456, and IC 4329A, while we use the spectrum from the Apache Point Observatory’s (APO) TRIPLESPEC instrument of Mrk 1239. The APO spectrum of Mrk 1239 has a higher spectral resolution than that of GRAVITY, and we decided to use it to better remove the narrow component and characterise the line profile. Our comparison between the two Mrk 1239 spectra confirms that they are both consistent in terms of line shape and flux levels. Figure 2 shows the utilised K-band spectra of all targets in this work.

Fig. 2. Average AGN flux spectra (black steps with grey error bars) of (a) Mrk 509 and (b) PDS 456 taken by GRAVITY, (c) Mrk 1239 taken by APO/TRIPLESPEC, and (d) IC 4329A taken by GRAVITY, all normalised to the continuum. The name of the emission lines for each target is shown in each panel.

3.3.1. GRAVITY spectra: Mrk 509, PDS 456, and IC 4329A

The flux spectra (1.95–2.45 μm) of Mrk 509, PDS 456, and IC 4329A were taken with GRAVITY in medium resolution mode (R = λ/Δλ ≈ 500) with 90 independent spectral elements, which were extracted and resampled (and interpolated, if needed) into 210 channels. We extracted the final spectra of Mrk 509, PDS 456, and IC 4329A by weight-averaging 3, 10, and 3 GRAVITY spectra, respectively. Due to the low number of individual spectra to average for Mrk 509, we cannot simply estimate the flux errors in each channel based on multiple observations. Hence, we decided to multiply the flux errors by a factor of 3 to balance the weight of the flux spectrum and differential phase spectrum during BLR modelling. We also did the same for IC 4329A, but we decided to remove the factor of the flux error as it does not greatly affect the results of our analyses. To remove the Mrk 509’s narrow Brγ emission component, we used the optical spectrum of Mrk 509 taken by the Kitt Peak National Observatory (KPNO) 2.1 m telescope (Shang et al. 2005). We utilised the optical blue spectrum (0.33–0.55 μm), adopting [O III]λ5008 Å as our narrow line template. For PDS 456, we assume all of its light originates from the BLR. Simpson et al. (1999) showed very little [O III] emission in PDS 456’s optical spectrum, which strongly suggests negligible narrow emission in its broad Paα line. Hence, no narrow-line component removal was performed. For IC 4329A, no narrow line template was used to fit its narrow component. Although there is available VLT/SINFONI K-band (1.93–2.47 μm) spectrum for IC 4329A, it was taken last Jun. 2019, which is too distant in time for variability to be neglected. Therefore, we simply fit the Brγ line profile with 6 Gaussian components and subtracted the component describing the narrow component at the central wavelength of Brγ. We note that for the flux spectrum of IC 4329A, any number of Gaussian components below 6 is insufficient to provide a good fit, while above 6 does not provide a better fit.

3.3.2. APO and XSHOOTER spectra: Mrk 1239

We acquired 3 Mrk 1239 spectra (0.94–2.47 μm) with the TRIPLESPEC instrument at the Apache Point Observatory (APO) 3.5 m telescope on two nights of 27th and 29th Dec., 2020. Since the GRAVITY interferometric observations were done ∼1–3 months after the APO spectroscopic observations, variability of spectral shape and intensity should be negligible as confirmed by Pan et al. (2021) after checking multi-epoch observations of Mrk 1239 in various wavelength ranges. The resolution of our Mrk 1239 flux spectrum is R = 3181, almost 6 times that of GRAVITY’s, making it a better choice for our analysis. We used a 1.1″ × 45″ slit for our APO observations. The spectra were reduced using the modified version of the Spextool package (Cushing et al. 2004). An A0V star was also used for telluric correction and flux calibration (Vacca et al. 2003).

To get a narrow line template suitable for removing the narrow component of Mrk 1239’s Brγ emission line, we observed Mrk 1239 with X-shooter in Dec. 2021 (PI: Shangguan, programme ID 108.23LY.001). The VIS (0.56–1.02 μm, R = 18 400, using slit dimensions of 0.7″ × 11″) spectrum of Mrk 1239 was reduced using the X-shooter pipeline version 3.5.3 running under the EsoReflex environment version 2.11.5. From this spectrum, we utilise [SII]λ6716 Å as a template for the narrow component.

We implement two methods to study the BLR of our targets: photocentre measurements (a model-independent method that reveals the BLR’s possible velocity gradient and offset with respect to the hot dust) and BLR modelling (reveals BLR properties and central BH mass). More details about each method are discussed in Sects. 4 and 5.

4. Measuring the BLR photocentres

We start our analysis of the BLR structure by measuring the photocentre positions of the line emission. This model-independent method provides a direct representation of the differential phases in the various baselines. By measuring the BLR photocentres, we can confirm whether there is an overall velocity gradient in the BLR and investigate the offset between the BLR and hot dust photocentres, which previous works have also measured (GRAVITY Collaboration 2020a, 2021a).

We first select the spectral channels where the Brγ/Paα emission dominates. This criterion gives us 9, 8, 7, and 9 spectral channels to fit for Mrk 509, PDS 456, Mrk 1239, and IC 4329A, respectively. The photocentre displacement for each channel can then be calculated using Eq. (2). We also estimate the astrometric accuracy of each target (listed in Table 2) by averaging the errors of each photocentre position.

The left column of Fig. 3 shows the BLR photocentres of the 4 targets as a function of wavelength. The photocentres show a velocity gradient in the BLR, indicated by the separation between the blueshifted and redshifted channels. To investigate the significance of this gradient, we fit all blueshifted and redshifted channels, assuming a single photocentre for each side. We call this our “2-pole” model. We also fit all spectral channels into one photocentre, which we call our “null” model. The results of the 2-pole/null model fitting are shown in Fig. 3, middle column. The resulting 2-pole model fittings reveal the separation between the red and blue poles (velocity gradient) of Mrk 509, PDS 456, Mrk 1239, and IC 4329A to be 48, 20, 30, and 37 μas, respectively. These correspond to 0.034, 0.064, 0.013, and 0.012 pc, respectively. An F-test was utilised to confirm the significance of the velocity gradient, with the null model as the null hypothesis and the 2-pole model as the alternative hypothesis. Mrk 509, PDS 456, and IC 4329A reveal p-values that are ≪0.05, indicating 6.7σ, 6.9σ, and 6.6σ significance, respectively. However, the red-blue photocentre offset of Mrk 1239 has a p-value of 0.88, meaning we cannot reject the null hypothesis of a single photocentre describing the BLR. In this case, the confidence in the offset between the red and blue photocentres and the resulting position angle begs a more detailed model fitting.

Fig. 3. Best-fit BLR photocentres of our four targets. The columns from left to right show the photocentres from the data, the “2-pole” and “null” photocentre models (described in Sect. 4), and the photocentres from the best-fit BLR model of each target. The red cross refers to the continuum photocentre. The direction of the radio jet is shown as a green solid line. The black arrow shows the position of the best-fit “null” model fitted from the best-fit BLR model. The ellipses around each centroid refer to 68% (1σ) confidence intervals. The colours refer to the corresponding velocity of each spectral channel. The (significant) separation of the red- and blueshifted photocentres confirm that we have resolved the BLRs of our targets. For Mrk 1239 and IC 4329A, the perpendicular alignment between the radio jet and red-blue photo centres from the “2-pole” model fitting of the data and the straight alignments of their model photo centres indicate that their BLR is rotation-dominated. On the other hand, for Mrk 509 and PDS 456, the radio jet and red-blue photocentres are not closely perpendicular, and the model photocentres are curved, which suggests that their BLRs are radial motion-dominated.

The BLR photocentres of our targets also reveal a systematic offset to the continuum photocentre (red cross in Fig. 3). This offset, which we call the “BLR offset”, is measured as the distance between the continuum photocentre and the BLR null photocentre, and is about 28, 17, 137, and 34 μas (0.020, 0.054, 0.057, and 0.011 pc) for Mrk 509, PDS 456, Mrk 1239, and IC 4329A, respectively. The offsets are mostly close to perpendicular to the BLRs’ velocity gradient, which is similar to what was observed with IRAS 09149-6206 (GRAVITY Collaboration 2020a). To avoid the impact of possible data correlation between neighbouring channels during photocentre measurements, we re-binned the spectra by a factor of 2. The same method still yields a significance > 4.5σ (except Mrk 1239), regardless of the inclusion of the bluest/reddest channels that are furthest from the photocentres of the other channels.

On the right column of Fig. 3, we also show the model photocentres, which are produced from the mock differential phase and flux spectra of the best-fit BLR model (see Sect. 7). These model photocentres show the average photocentres of the BLR clouds of the best-fit BLR model in each wavelength channel and provide another way to see the cloud distribution, which can be linked to the differential phase data. These are discussed in more detail in Sect. 7, where they are used to explain the observed differential phases of our targets.

Since we were able to constrain the BLR velocity gradient of most of our targets, we can compare the alignment of the red and blue photocentres with the jet alignment, as revealed by radio observations. The orientation of the radio jet is shown as the green line in the middle 2-pole/null model fitting results. The radio jet orientations are taken from Ulvestad & Wilson (1984) for Mrk 509, Yang et al. (2021) for PDS 456, Doi et al. (2015) and Orienti & Prieto (2010) for Mrk 1239, and Colbert et al. (1996) for IC 4329A. It is clear that the red-blue photocentre orientations of Mrk 1239’s and IC 4329A’s BLRs are perpendicular to their radio jets, indicating that their BLRs are rotating, similar to what is found for 3C 273 (GRAVITY Collaboration 2018). However, Mrk 509 and PDS 456 reveal that the red-blue photocentre orientations of their BLRs are closely aligned to their radio jets.

To conclude our photocentre fitting analysis, we find significant (> 5σ) velocity gradients between the blueshifted and redshifted channels of all our targets except Mrk 1239. These velocity gradients are perpendicular to the BLR offsets. Our comparison between the red and blue photocentres and the orientation of the radio jets shows that the BLRs of Mrk 1239 and IC 4329A are rotating. However, Mrk 509 and PDS 456 show otherwise. Photocentre fitting cannot explain this phenomenon, let alone fully constrain the physical BLR size, as the photocentre offsets average over the emission of the channel bandpass and will therefore underestimate the BLR size. Therefore, a flexible model fit to the full differential phase spectra must be utilised to investigate our targets’ BLR structure and kinematics.

5. Description of the BLR model

We follow the BLR model described in Pancoast et al. (2014) (henceforth called the Pancoast model). This model defines the BLR as a collection of non-interacting clouds encircling the central SMBH. This is a convenient way to model the kinematics and structure of the BLR, but we do not mean to imply that the BLR specifically comprises discrete clouds. The physical structure of the BLR is described by several parameters that dictate the clouds’ position and motion. Table 2 of GRAVITY Collaboration (2020a) summarises the parameters included in the model and their possible ranges of values. We use the same priors as GRAVITY Collaboration (2020a), with a few exceptions, which are discussed later. Hence, we refer the readers to GRAVITY Collaboration (2020a) for the full list of parameters in the model.

In the model, R_BLR is the average (emissivity-weighted) BLR radius, and β controls the radial cloud distribution, which can be Gaussian (0 < β < 1), exponential (β = 1), or heavy-tailed and steep inner profile (1 < β < 2). The fractional inner radius is defined as R_min/R_BLR, where R_min is the minimum or inner radius of the BLR. The BLR can be inclined relative to the plane of the sky by an inclination angle (i) defined such that face-on is i = 0° and edge-on is i = 90°. The half-opening angle of the disc, θ₀, describes the thickness of the BLR: θ₀ = 90° produces a spherical BLR, while θ₀ = 0° produces a thin disc. Finally, the BLR is rotated within the plane of the sky by the position angle, PA. The parameter γ controls the angular concentration of BLR clouds relative to θ₀. When set to 1, it pertains to the uniform case with the clouds equally distributed as a function of angular height. Increasing γ > 1 redistributes the BLR clouds more and more to the outer faces of the disc, with γ = 5 as the maximum possible value.

Asymmetries can also be introduced in the BLR model through two effects. We can introduce anisotropic emission from each individual BLR cloud via the parameter κ. We assign a weight “w” to each cloud that defines how much BLR emission is directed into the LOS:

$$\begin{array}{c}\hfill w=0.5+\kappa cos\varphi \end{array}$$(3)

where κ ranges from −0.5 to +0.5. When κ > 0, it results in preferential emission from the BLR’s near side (side closer to the observer), while κ < 0 results in preferential emission from the BLR’s far side. The parameter ϕ describes the angle between the LOS of the observer and the LOS of the BLR cloud to the central ionising source. Anisotropic emission from the near side could be caused by BLR clouds situated closer to the observer obstructing the emission from the gas farther away. On the other hand, preferential emission from the far side of the BLR could be caused by self-shielding within individual BLR clouds, resulting in emission only at the back towards the central continuum source.

Finally, the transparency of the midplane can be modelled with the parameter ξ, which ranges from 0 to 1. When ξ = 1, the clouds are evenly distributed on both sides of the equatorial plane, while ξ = 0 means that the emission from the clouds behind the equatorial plane is obscured. The physical cause of possible BLR mid-plane opacity is not well-understood (GRAVITY Collaboration 2020a), and it should be interpreted cautiously.

Two parameters control the effect of radial motion in the BLR clouds: f_ellip and f_flow. f_ellip is the fraction of BLR clouds that have bound elliptical orbits. The rest (1-f_ellip) have much more elongated orbits dominated by radial motion. f_flow is a flag that controls whether these elongated orbits are inflowing (0 < f_flow < 0.5) or outflowing (0.5 < f_flow < 1). The parameter f_flow acts as a “binary flag” and is thus not a continuous variable. We also fit the central wavelength of the line, the line peak flux, and the BLR centre. However, we do not focus our discussion on these parameters.

We consider two variations of the Pancoast model in our work. One is termed the circular model. This model sets κ = 0 and ξ = 1 to ensure that all clouds emit isotropically and are uniformly distributed above and below the BLR mid-plane. No angular asymmetry is considered in this model; that is, γ = 1. All BLR clouds are then subject to circular Keplerian rotation (i.e., f_ellip = 1). The second model is the elliptical/radial model. In this model, angular asymmetry (γ) and all parameters corresponding to the asymmetrical properties of the BLR (κ, ξ) are fitted. Inflowing/outflowing clouds (f_ellip and f_flow) are fitted as well (together with the Keplerian clouds). We fit the circular and elliptical/radial models to our 4 targets, but we only show the best-fit results for each target. We discuss the reason for our choice of BLR model during the fitting of each target in Sect. 7.

For the exceptions compared to GRAVITY Collaboration (2020a), we fix the angular and radial standard deviations of circular and radial orbit distributions and the standard deviation of turbulent velocities to zero because they do not affect the model fitting. We limit the parameter space of the inclination angle to be Uniform(cos i(0, π/4)). We tested the fitting for all targets with a much larger range, Uniform(cos i(0, π/2)), and found that all targets except Mrk 509 and IC 4329A prefer a smaller inclination. Hence, we decided to use Uniform(cos i(0, π/2)) for the parameter space of the inclination angle of Mrk 509 and IC 4329A, while we adopt the smaller parameter space, Uniform(cos i(0, π/4)), to the rest of the targets to avoid multiple peaks in the posterior distribution. More discussion about the inclination angle of Mrk 509 and IC 4329A is presented later (in Sects. 7.1 and 7.4, respectively).

We fit the models to the total flux line profile and differential phase curves from all baselines simultaneously. The chosen wavelength channels for fitting each target cover the wavelength range where Brγ/Paα emission flux is relevant. All models utilise 2 × 10⁵ clouds, which are randomly given parameter values based on their model distributions.

The interferometric data and prior information about the source and model are used to infer the best-fit model parameters via Bayes’ theorem. We use the Python package, dynesty (Speagle 2020), to implement posterior distribution sampling (random walk) of 2000 live points for each BLR fitting. dynesty allows the nested sampling algorithm to be utilised, which is powerful in estimating Bayesian evidence and dealing with complex models with presumably multimodel posteriors. The Bayesian evidence or marginal likelihood, Z, measures how well the probability distributions are constrained. The larger the Z, the stronger the constraint to the probability distributions. The Bayes factor, or the ratio of the Bayesian evidences, is then calculated to identify which model is more apt to fit the data.

We prefer to fit two models (inflowing and outflowing elliptical/radial model) separately and compare their Bayes factor, or the ratio of the outflowing model evidence to the inflowing model evidence, due to the nature of f_flow parameter as a binary flag. We do this by fitting an elliptical/radial model with f_flow restricted to a value between 0–0.5 (for inflowing clouds, representing an inflow elliptical/radial model) or 0.5–1.0 (for outflowing clouds, representing an outflowing elliptical/radial model). If the Bayes factor > 1, the source prefers outflowing radial motions; otherwise, it prefers inflowing radial motions. For all of our fittings, we use the nesting sampling algorithm (NestedSampler) as the sampler and random walk (rwalk) as the sampling method.

In this work, we report the best-fit values and their 68% (1σ) credible intervals of each parameter. We also show the model-inferred photocentre fitting results in Fig. 3, right column, where we simulate the differential phase and flux spectra from the best-fit BLR models of our targets. We observe that our model-inferred photocentre fitting results are consistent with what we get from fitting the photocentre from our data.

6. Importance of differential phase data in BLR fitting

6.1. Resulting posterior distribution from BLR fitting with and without phase data

As discussed in the previous section, the differential phase and flux spectra are our key inputs for BLR fitting. The differential phase measures the broad emission line’s photocentre shift at different wavelength channels with respect to the continuum, providing spatial information. On the other hand, the flux spectrum gives information on the distribution of gas velocities in the BLR (see Raimundo et al. 2019, 2020, and references therein), as well as the inclination angle of the BLR, which is degenerate with the observed BH mass (Rakshit et al. 2015; GRAVITY Collaboration 2018). The work of Pancoast et al. (2014), from which our model is derived, performed fitting of broad emission lines (usually Hβ). In this section, we investigate the importance of adding the differential phase spectrum in fitting the broad emission line with our BLR model.

First, we fit our four targets with and without their differential phase spectra, assuming the same priors as the fitting with the differential phase spectra, and then compare their posterior distributions. We focus on the posterior distributions of the BLR radius, log R_BLR [ld], and BH mass, log M_BH [M_⊙] since these are the parameters that cannot be constrained with only a line profile. Figure 4 shows the comparison between the posterior distributions of our BLR fitting with (orange line) and without (blue line) differential phase spectra for our 4 targets. It is clear that without the phase data, our BLR fitting fails to constrain these parameters. Our simple comparison shows that the phase data is crucial for unveiling more accurate estimates of our targets’ BLR size and BH mass and more precise pictures of their BLR geometry and kinematics. This still holds true even for cases when there is no differential phase signal detected.

Fig. 4. Posterior distributions of BLR radius (left column) and BH mass (right column) resulting from BLR fitting without differential phase data (blue histogram) and with differential phase data (orange histogram). Each row corresponds to each target; from top to bottom: Mrk 509, PDS 456, Mrk 1239, and IC 4329A. The histograms are normalised such that they integrate into 1. It is clear that the spatial information from the differential phase significantly constrains both the BLR radius and BH mass.

6.2. Expected differential phase signals of Mrk 1239 and IC 4329A

In the previous subsection, we highlight the importance of differential phase data in our BLR model fitting even if there is no phase signal detected. We investigated further by estimating the expected phase signals for Mrk 1239 and IC 4329A if only their flux profile were available for fitting. We focus on these two objects because, among our 4 targets, they have very weak signals that are still below their noise level (∼0.05° and ∼0.10° for Mrk 1239 and IC 4329A, respectively). To constrain the possible phase signals of Mrk 1239 and IC 4329A without the phase data, we ran our BLR fitting for these two objects without their differential phase spectra while fixing their BH masses into their possible maximum and minimum values from previous works. For Mrk 1239, the maximum BH mass is log M_BH [M_⊙] = 7 from Pan et al. (2021), while the minimum BH mass is log M_BH [M_⊙] = 5.9 from Kaspi et al. (2005) and Greene & Ho (2005). For IC 4329A, the maximum BH mass is log M_BH [M_⊙] = 7.6 from Bentz et al. (2023), while the minimum BH mass is log M_BH [M_⊙] = 6.8 from Kaspi et al. (2000). The expected averaged phase signals are shown in Fig. 5 as a 1σ shaded region. The chosen baselines to average are UT4-UT3, UT4-UT2, and UT4-UT1 for Mrk 1239, and UT4-UT2, UT4-UT1, and UT3-UT2 for IC 4329A. To create Fig. 5, the PA of each object were also fixed to their best-fit values (see best-fit results in Sect. 7). Since there is a degeneracy between the BH mass, BLR size, and inclination angle (Rakshit et al. 2015; GRAVITY Collaboration 2018), the expected phase signal will then reflect the possible BLR sizes and inclination angles for these targets; that is, the larger the expected phase signal is, the bigger the inferred BLR size and inclination angle are.

Fig. 5. Expected differential phase signals for Mrk 1239 (left panel) and IC 4329A (right panel) when their line profiles were fitted with our circular model assuming their maximum (red shaded region) and minimum (blue shaded region) BH masses were fixed during fitting. The assumed BH masses for Mrk 1239 and IC 4329A are taken from previous literature as discussed in the text. The shaded regions refer to 1σ uncertainty. The chosen baselines to average are the same as the ones used to create Fig. 6: UT4-UT3, UT4-UT2, and UT4-UT1 for Mrk 1239; UT4-UT2, UT4-UT1, and UT3-UT2 for IC 4329A. The black points refer to the maximum and minimum observed phase signals at the line region of Mrk 1239 and IC 4329A after averaging the differential phase spectra from the same chosen baselines. The horizontal lines show the wavelength range of the maximum and minimum observed phase signals.

Figure 5 shows the range of differential phase signals that we are expecting based on what was previously known. The expected phase signals for the two objects span a large range when the maximum BH mass is assumed, indicating the wide range of possible BLR radius and inclination angles. The resulting absolute phase signal goes as high as ∼0.09° for Mrk 1239 and ∼0.10° for IC 4329A. The opposite is true when the minimum BH mass is assumed; the resulting absolute phase signals are very small (∼0.005° for Mrk 1239, and ∼0.02° for IC 4329A) and cover a very narrow area in the differential phase parameter space. The areas covering the expected phase signals for Mrk 1239 and IC 4329A when the maximum BH masses are assumed are much larger than the typical signals shown in the differential phase spectra of these objects (see Fig. 1). This emphasises the importance of phase data in our fitting; without the phase data, the constraining power of our fitting is lessened, and the resulting posterior distributions of their BLR radius and inclination angle span a much wider range. The rather weak signals of Mrk 1239 and IC 4329A have a discriminating power, allowing us to limit the possible values of their BLR radii and inclination angles. Another way to look at this is that the small phase signals of Mrk 1239 and IC 4329A rule out higher values of inclination angle because these would lead to stronger phase signals.

The expected phase signals of Mrk 1239 and IC 4329A could also be used to infer the possible BH masses of these targets by comparing them with their observed phase signals. In Fig. 5, we plot the average differential phase on the line region channels located at the left and right of the line centre. These represent the average observed phase signals of Mrk 1239 and IC 4329A. For Mrk 1239, its average observed phase signal runs from −0.04° to +0.02°. Based on the expected phase signals in Fig. 5, this indicates that Mrk 1239 may not have as low of a BH mass as Kaspi et al. (2005) measured, and instead, it may have a large BH mass close to what Pan et al. (2021) measured. The average observed phase signal of IC 4329A (from −0.06° to +0.02°) also suggests that IC 4329A may not possess such a low BH mass similar to what Kaspi et al. (2000) measured. This indicates that IC 4329A may have a large BH mass similar to the value Bentz et al. (2023) measured. The BLR fitting results discussed next section will confirm if these points are true or not.

7. BLR modelling results

Figure 6 summarises our BLR modelling results. We compare the differential phase spectra from our data with the ones predicted by the model, using the weighted average of the differential phase spectra from the longest baselines and/or the baselines that show the strongest or most prominent BLR phase signal after continuum phase removal (note that all the baselines were used for the actual BLR fitting; we only do this step for clearly visualising the averaged phase of each target). We show the 16-to-84 percentile range (1σ confidence interval) of the model phase and flux spectra as a red-shaded region by randomly sampling the posterior distribution of our fitting results 100 times.

Fig. 6. Summary of BLR fitting results for (a) Mrk 509, (b) PDS 456, (c) Mrk 1239, and (d) IC 4329A. The observed averaged differential phase spectra with the continuum phase signal removed are shown in blue, while the median differential phase spectra from the best-fit BLR model are shown in red (see text for the details of averaging the differential phase spectra for each target). Above each panel, the observed line profile (black) is shown together with the median best-fit model line profile (red). The central wavelengths of the line profiles are shown by the grey vertical dashed lines. The red-shaded region shows the 16-to-84 percentiles of the averaged differential phase spectra and flux spectra from the best-fit model, which was created by getting randomly selected samples of each parameter in 100 different instances and then recreating the averaged differential phase spectra. Mrk 509 and PDS 456 clearly show an asymmetrical signal typical for the Pancoast model, while Mrk 1239 and IC 4329A show the S-shape signal typical for the Keplerian model.

The differential phase signal will vary over different baselines due to their different orientation with respect to the observed target. Therefore, for visualisation purposes, we have the freedom to choose which baselines we should average to compare the data and model fitting results. For Mrk 1239, the averaged differential phase spectrum is derived from averaging the differential phase spectra from three baselines: UT4-UT2, UT4-UT1, and UT4-UT3. For IC 4329A, the chosen baselines are similar to that of Mrk 1239, except that UT3-UT1 is chosen instead of UT4-UT3. For Mrk 509 and PDS 456, we averaged the spectra from all baselines. The observed differential phase includes both the continuum phase and the BLR differential phase. We are only interested in the latter. Hence, we subtract the former based on our model fitting. Mrk 509 and PDS 456 show a strong asymmetric signal in their averaged BLR differential phase. Mrk 509 also shows a slight asymmetry in its flux spectrum. On the other hand, Mrk 1239 and IC 4329A reveal weak BLR signals, which are below their respective noise. Considering this, we decided to fit the elliptical/radial model for Mrk 509 and PDS 456 data, while we fit the circular model for Mrk 1239 and IC 4329A data. The low signal-to-noise ratio (S/N) data of Mrk 1239 and IC 4329A cannot be well-constrained by the elliptical/radial model, and even if we fit them with the elliptical/radial model, the resulting BH mass and BLR size are not much different given our uncertainties, similar to what is found for IRAS 09149-6206 in GRAVITY Collaboration (2020a).

Tables 3 and 4 show the best-fitting parameters of the elliptical/radial and circular models for our 4 targets. We discuss each target’s BLR modelling results (i.e., averaged differential phase spectra, best-fit parameter values, and their interpretation) in the following subsections.

7.1. Mrk 509

Mrk 509 is best fit with the elliptical/radial model, which better reproduces its asymmetric flux spectrum. The Bayesian evidence from the elliptical/radial Pancoast model, which is ∼100× greater than that from the circular Pancoast model, strongly indicates the need for including both asymmetric emission and radial motion. Figure A.1 shows the corner plot of the resulting posterior distributions of all fitted parameters. The inferred inclination angle is ∼69°. Although it is quite high for a face-on target, we note that with the presence of radially moving clouds, the definition of the inclination angle is not as straightforward as that of a BLR possessing only rotation-dominated clouds.

The inferred BLR radius for Mrk 509 is log R_BLR [ld] = 2.29${}_{-0.26}^{+0.01}$. This is slightly larger than what is measured previously by other works (Peterson et al. 1998). However, the best-fit BH mass of Mrk 509 is log M_BH [$M_{\odot }]=8.{00}_{-0.23}^{+0.06}$, consistent with previous estimates (e.g., Peterson et al. 1998; Grier et al. 2013). The inferred inner BLR radius is log R_BLR, min [ld] = 0.93${}_{-0.39}^{+0.40}$.

Our BLR model for Mrk 509 prefers only ∼30% of the clouds to be in circular orbits, meaning the majority (∼70%) of the clouds have significant radial (highly elliptical) motion. The resulting Bayes factor after fitting an inflowing and outflow radial model is 32, indicating a strong preference for the outflow radial model over inflow. This supports the result that outflows dominate the BLR of Mrk 509.

Mrk 509 shows an asymmetric line profile with a blueshifted shoulder, and a negatively peaked differential phase signal shifted redwards of the Brγ central wavelength (see Fig. 6a). Such profiles can only be produced by a BLR possessing asymmetric properties. Our best fit angular distribution parameter is γ ∼ 4.1, indicating that most of the BLR clouds are on the outer faces of the BLR disc, and fewer clouds are located in the midplane of the BLR. This is physically consistent with what one would expect for an outflowing BLR. Our best fit BLR model prefers a BLR with low mid-plane transparency (ξ ∼ 0.21) and a moderate preference for the emission to originate from the far side of the clouds (κ ∼ −0.18). Figure 7a shows the edge-on view of Mrk 509, emphasising the lack of clouds seen below the midplane. We note that this does not necessarily mean that there are actually fewer clouds below the midplane since we are only modelling the observed broad-line emission. The high inclination angle allows the observer to look directly into the edge of the disc, causing the observed differential phase on the blueshifted side to be small. On the other hand, the preferential emission on the far side (where the redshifted BLR clouds are mostly found, as shown in Fig. 7) creates a larger phase on the redshifted side of the line centre. In addition, the midplane obscuration affects more redshifted clouds than blueshifted clouds due to the inclination angle, and therefore produces an increase on the blue wing of the flux spectrum. Therefore, the best-fit BLR model of Mrk 509 (e.g., high inclination angle, large thickness, presence of BLR asymmetry) can recreate the asymmetric flux and differential phase spectra of Mrk 509. The combination of the effects of midplane obscuration and a significant fraction of outflowing BLR clouds causes the photocentres to be curved around the centre (right column, first row in Fig. 3), producing a strong dip in the averaged differential phase signal. This also explains why the red-blue photocentres of Mrk 509’s BLR are not perpendicular to its radio jet, as the red-blue photocentres are dominated by the clouds dominated by radial motion and not by tangential/Keplerian motion.

Fig. 7. Edge-on views of the best-fit outflow model for (a) Mrk 509 and (b) PDS 456. For clarity, the PA for both configurations is adjusted to 180°, and the BLR centre is positioned to the origin. The LOS at the best-fit inclination angle is shown as a black arrow. The midplane, represented as a green line, is defined as the line passing through the origin and is perpendicular to the LOS at i = 0°. The regions above (closer to the observer) and below (farther to the observer) the midplane are also labelled. The colours of each cloud refer to their LOS velocity, while their sizes reflect the weight given to each cloud towards the total emission; the larger the size of the cloud, the more this cloud contributes to the broad-line emission. The edge-on view of the BLRs of Mrk 509 shows that its asymmetry is due to a highly obscuring midplane; for PDS 456, a preference for broad-line emission to come from the far side of the BLR.

7.2. PDS 456

Similar to Mrk 509, PDS 456 is best fit with the elliptical/radial model. This model better fits PDS 456’s asymmetric differential phase spectrum, which is difficult to achieve with the circular model. The Bayesian evidence for the elliptical/radial model is ∼14× greater than that from the circular Pancoast model fitting. Figure A.2 shows the resulting corner plot of its BLR model fitting. It has a small inclination angle (i ∼ 13°), and the opening angle (θ₀ ∼ 42) suggests that its BLR is moderately thick. This indicates that the system, despite the thickness, is unobscured. The BLR cloud distribution also has a steep profile (β ∼ 1.83). Our BLR modelling with PDS 456 infers a BLR radius of log R_BLR [ld] = 2.49${}_{-0.38}^{+0.08}$; the largest among our sample. The best-fit BH mass of PDS 456 is log M_BH [$M_{\odot }]=8.{23}_{-0.49}^{+0.01}$. Our BH mass is a factor of 10 smaller than that measured by Reeves et al. (2009) and Nardini et al. (2015) via scaling relations. This is similar to what Abuter et al. (2024) concluded via dust size measurements as they also measured a smaller BH mass than previous works, albeit 0.5 dex higher than our result. This also suggests that the discrepancy in the BH mass between our work and that estimated by Abuter et al. (2024) is just due to the assumed virial factor to calculate the latter. The inferred BLR inner radius is log R_BLR, min [ld] $1.{16}_{-0.38}^{+0.42}$.

Our modelling results with PDS 456 emphasise that almost half of the BLR system is dominated by radial motions, as the best-fit value of f_ellip ∼ 0.52. Letting f_flow as a free parameter causes the BLR fitting results to slightly prefer inflowing radial motions (f_flow < 0.5). However, the Bayes evidence of fitting with f_flow fixed to inflow and outflow is equal to 1.818, which means that PDS 456 does not show a preference between inflow and outflow radial motions. Therefore, we conclude that the data cannot distinguish between the inflow and outflow model, and we choose to fix f_flow > 0.5 for PDS 456 as other observations reveal outflowing signatures at both smaller and larger scales than the BLR (Reeves et al. 2003; Luminari et al. 2018; Bischetti et al. 2019).

PDS 456 also exhibits a negative averaged differential phase spectrum almost centred at the Paα central wavelength, similar to Mrk 509. Its phase signal is more asymmetric compared to that of Mrk 509, as it shows a more extended redward tail. Again, we look at the best-fit asymmetry parameters and the model photocentres to explain this phase signal. The best-fit value of ξ is ∼0.75, indicating that the midplane is mostly transparent. However, the best-fit value of κ is ∼ − 0.44, suggesting that the BLR clouds emit preferentially from the far side. The effects of these parameters are evident in Fig. 7b, where the edge-on view of its BLR is shown. Here, the number of observable BLR clouds above and below the midplane are very similar, as expected for the best-fit ξ. However, the clouds below the midplane are displayed with larger sizes compared to those above the midplane, which reflects the best-fit κ. The preferential emission of the far side of the BLR, where the redshifted clouds reside, creates a larger phase on the redshifted side of the line centre. In addition, similar to Mrk 509, the combination of the effects of outflowing BLR clouds and their preferential emission cause the locus of model photocentres of PDS 456 to look curved (right column, second row, in Fig. 3). The null model BLR photocentre of PDS 456 is located above the BLR photocentres, which creates the negative differential phase signal. Similar to that of Mrk 509, the dominance of radial motion in the BLR clouds causes the red-blue photocentres to be shifted less perpendicular to its radio jet. In conclusion, the cause of the asymmetric signal in PDS 456’s average differential phase spectrum is due to (1) the BLR clouds’ outflowing radial motions and (2) the preference of the Paα emission to originate from the far side of the BLR.

7.3. Mrk 1239

The red- and blueshifted photocentres of Mrk 1239 show an insignificant separation (see Fig. 3). The photocentre fitting results of Mrk 1239 suggest that the differential phase is dominated by the continuum phase, and the BLR differential phase is very moderate (low S/N). However, the perpendicular orientation of the BLR velocity gradient and continuum photocentre and the radio jet is consistent with the picture of a rotating BLR. Despite the relatively weak signal of Mrk 1239, we argue that this aids in adding further constraint to our BLR fitting compared to when we only fit its flux spectrum (as proven in Sect. 6). Given the low signal, we choose to fit the data with the simpler, circular model. The flux and differential phase spectra of Mrk 1239 are well-fit with the circular model as shown in Fig. 6c. While there is no clear differential phase signal for Mrk 1239, the combined fitting of the emission line profile and differential phase allows for meaningful constraints on the BLR size and SMBH mass as shown in the corner plot (see corner plot in Fig. A.3). While Fig. 6c does not show a clear differential phase signal, the very lack of a signal combined with the impressively low noise level still allows the model to converge and constrain the key BLR and SMBH properties of Mrk 1239.

The measured BLR PA of Mrk 1239 is ∼197°, similar to what is measured for the red-blue photocentres of Mrk 1239, even if their separation is insignificant. The best-fitting parameters indicate that the data favours a face-on disc (i ∼ 11°) with a thickness of θ₀ ∼ 42°. This geometry is consistent with the results of Pan et al. (2021), who suggest the observed broad emission of Mrk 1239 is actually reprocessed through scattering by polar dust. Effectively, we are observing from the LOS of the polar dust, which has a face-on viewing angle.

The cloud distribution follows an exponential profile (β ∼ 1.21), while the mean radius is measured to be log R_BLR [ld] = 1.77${}_{-0.74}^{+0.03}$. The inferred BH mass is log M_BH [$M_{\odot }]=7.{47}_{-0.92}^{+0.15}$. Our inferred BLR radius is consistent within uncertainties with those published previously (Du et al. 2014a; Abuter et al. 2024), and so is our BH mass (Kaspi et al. 2005; Du et al. 2014a; Abuter et al. 2024). Lastly, the inferred BLR inner radius is log R_BLR, min [ld] = $0.{64}_{-0.28}^{+0.27}$. However, we caution that the BLR size of Mrk 1239 is marginally constrained, as its lower error is relatively larger than its upper error. This means that the BLR size and, therefore, the BH mass of Mrk 1239 has a higher chance of being lower than our reported best-fit values.

As discussed in Sect. 6.2, the average observed differential phase signal of Mrk 1239 suggests that its BH mass may be similar to the Pan et al. (2021) measurement. Indeed, from the modelling, the best-fit BH mass is log M_BH [M_⊙] = 7.47, close to the log M_BH [M_⊙] = 7 measurement of Pan et al. (2021).

7.4. IC 4329A

Even though the noise level of individual differential phase spectra in each baseline is ∼0.1°, the S-shape signal is not as obvious as that of Mrk 1239 (after removing the continuum phase). Nevertheless, the radio jet of IC 4329A is shown to be perpendicular to the red-blue photocentre orientation of its BLR, indicating the rotation of its BLR and we therefore fit the data with the simpler circular model. As for Mrk 1239, the combination of the emission line profile and differential phase data strongly constrains the BLR structure and SMBH mass. Similar to Mrk 1239, IC 4329A does not show a clear phase signal (Fig. 6d). As discussed in Sect. 6, this lack of clear phase signal acts as a strong constraining factor in IC 4329A’s BLR and SMBH properties. Without the differential phase spectrum, the BLR properties of IC 4329A will not be well constrained.

The cloud distribution of its BLR (β ∼ 1.81) suggests it has a steep inner radial profile. The inclination angle is i ≈ 54°; similar to Bentz et al. (2023), the high inclination angle of IC 4329A suggests that the AGN system and its host galaxy disc are misaligned significantly. The opening angle is θ₀ ≈ 54°, indicating that the BLR is thick. We resolved the BLR of IC 4329A, with a BLR radius of log R_BLR [ld] = 1.13${}_{-0.23}^{+0.10}$ and a BH mass of log M_BH [$M_{\odot }]=7.{15}_{-0.26}^{+0.38}$. The inferred BLR inner radius is log R_BLR, min [ld] = $0.{64}_{-0.28}^{+0.27}$.

Bentz et al. (2023) recently published results of a new RM campaign for IC 4329A which spanned about 5 months during 2022. They used CARAMEL (Bentz et al. 2021b, 2022) to fit their RM data with the BLR model from Pancoast et al. (2014), which is similar to what we use. Our best-fit BLR size, BH mass, i, and θ₀ are all consistent with RM measurements from Bentz et al. (2023). Indeed, the average observed phase signal of IC 4329A spans the same range as the expected phase signal of IC 4329A with its BH mass fixed to what Bentz et al. (2023) measured.

8. The GRAVITY-AGN radius-luminosity relation

One of the main objectives of our GRAVITY-AGN large programme is to place AGN properties derived by GRAVITY in the context of AGN scaling relations previously derived by other methods, namely the R-L and M_BH-σ_* relations. In this section, we show the first non-RM-derived R-L relation. To start, we have the BLR radius and BH mass measurements of our 4 targets, and adding the measurements of the 3 previously published targets (IRAS 09149-6206 from GRAVITY Collaboration 2020a, 3C 273 from GRAVITY Collaboration 2018, and NGC 3783 from GRAVITY Collaboration 2021a), gives us 7 AGNs to plot in the R-L space, now spanning a wide range of luminosity (log λL_λ(5100 Å) ∼ 43.0–46.5). We fit these data points with a power-law relation similar to that of Bentz et al. (2013):

$$\begin{array}{c}\hfill \log (R_{\mathrm{BLR}}/\mathrm{ld})=K+\alpha \log (\lambda L_{\lambda }/{10}^{44}\mathrm{erg}{\mathrm{s}}^{-1}).\end{array}$$(4)

We employed the LINMIX algorithm (Kelly 2007), a package that utilises a hierarchical Bayesian approach to linear regression with measurement errors and a parameter, σ², which encapsulates the intrinsic random scatter around the regression. We calculate the 16th, 50th, and 84th percentile of the resulting posterior distributions of each parameter to get their best-fit values and their 1σ uncertainties. We use the optical AGN luminosities of our targets listed in Table 1 and the best-fit BLR sizes of our targets to produce our GRAVITY-AGN R-L relation. However, we emphasise that Mrk 1239 is known to be a physically complicated object based on various observed properties (e.g., Balmer decrement, polarisation) leading to the high obscuration of its nucleus (Goodrich 1989; Doi et al. 2015). Its observed broad lines and optical continuum are also suggested to originate from scattering by polar dust around the BLR (Pan et al. 2021). Therefore, we cannot directly use the observed λL_λ(5100 Å) of Mrk 1239. We opt to use the extinction-free λL_λ(5100 Å) of Mrk 1239 from Pan et al. (2021).

Figure 8 shows the result of our R-L relation fit compared with that of Bentz et al. (2013). The best-fit values for our GRAVITY-observed AGNs are $K=1.{69}_{-0.23}^{+0.23}$, $\alpha =0.{37}_{-0.17}^{+0.18}$, and ${\sigma }^2=0.{23}_{-0.13}^{+0.48}$. The large uncertainties in our best-fit parameters are inevitable due to our small sample size. By observing more AGNs with GRAVITY, we can increase the sample size and, therefore, better constrain the R-L relation independent of RM (see Sect. 12 for more discussion of these future prospects). We emphasise that our approach allows us to extend the R-L relation to much higher luminosities without the need for decade-long time baselines to observe their expected year-long time lags. Our results suggest that the GRAVITY R-L relation has a flatter slope (∼0.4) compared to that of the canonical R-L relation from Bentz et al. (2013), which has a slope of ∼0.5. Woo et al. (2024) recently showed a similar conclusion as this work by focussing on reverberation mapping of high-luminosity AGNs (10⁴⁴ L_⊙ < λL_λ(5100 Å) < 10⁴⁵ L_⊙, with a few sources possessing λL_λ(5100 Å) > 10⁴⁵ L_⊙).

Fig. 8. Logarithm of the BLR radius as a function of log λLλ(5100 Å) (R-L relation). The blue line shows the best-fit line derived from all 7 GRAVITY-observed AGNs, and the blue shaded region is its 1σ confidence interval. The AGN luminosity is taken from Table 1. The dashed line represents the best-fit radius luminosity relation from Bentz et al. (2013). The 3 published AGNs (IRAS 09149-6206 from GRAVITY Collaboration 2020a, 3C 273 from GRAVITY Collaboration 2018, and NGC 3783 from GRAVITY Collaboration 2021a) are shown as a red star, red diamond, and red square, respectively. Our 4 targets are shown in differently coloured filled squares with 1σ errors. Our best-fit R-L (power-law) relation has a slope of $\alpha =0.{37}_{-0.17}^{+0.18}$, which is lower but, given the current large uncertainty, still consistent within uncertainties with the slope of R-L relation fit from Bentz et al. (2013; $\alpha =0.{533}_{-0.033}^{+0.035}$). If we follow the prescription from photoionisation models, the BLR sizes of 3C 273 and PDS 456 will be adjusted (shown as open symbols), and the fitted GRAVITY R-L relation is shown as the green line.

We also consider the relative time lags of the emission lines with respect to the Hβ time lag based on the predictions of photoionisation models. These models (e.g., Netzer 2020; Korista & Goad 2004) predict that the BLR size is dependent on certain parameters such as the optical depth of the observed emission line and the variability and the photoionisation flux of the central engine. Kuhn et al. (2024) recently confirmed the consistency of observed relative time lags of several emission lines (e.g., Hα, He I) with respect to Hβ time lag with what the radiation pressure confined (RPC) BLR model by Netzer (2020) described. According to this model, the time lags of Brγ and Paα are about 1.0–1.2 and 1.5–1.7 times larger than that of Hβ. It is clear, therefore, that BLR sizes derived from Paα are greatly affected, and thus, an adjustment must be taken into consideration. We show the adjusted BLR sizes of 3C 273 and PDS 456 (the objects with Paα line profiles) as open symbols in Fig. 8. If we consider these adjusted BLR sizes in our GRAVITY R-L relation fit, we get $K=1.{66}_{-0.24}^{+0.24}$, α = 0.30${}_{-0.19}^{+0.19}$, and σ² = 0.27${}_{-0.16}^{+0.54}$. We, therefore, conclude that the resulting slope of the R-L relation will be shallower if we consider prescriptions from photoionisation models and adjust BLR sizes measured from Paα lines. As for the effect of this “BLR stratification” phenomenon on the rest of the BLR model parameters, future work (Kuhn et al. 2024) on NGC 3783 has shown that the geometry beyond the radial distribution of the line emission (e.g., inclination angle, opening angle, black hole mass) are consistent among the different emission lines (Hα, Hβ, Paα, Paβ, HeI).

8.1. Effect of Eddington ratio in our results

An important aspect to discuss with the GRAVITY R-L relation fit is the Eddington ratio of our targets since previous works have suggested that it is the “third” parameter in the R-L relation, and it can explain the smaller observed BLR sizes compared to what the canonical R-L relation expects (Du et al. 2018; Du & Wang 2019) especially for PDS 456 whose best-fit BH mass and BLR size are much smaller than previous estimates. In this work, we calculate the Eddington ratio as λ_Edd = L_bol/L_Edd, where L_Edd = 1.3 × 10³⁸ (M_BH/M_⊙) erg s⁻¹. We also calculate the dimensionless accretion rate, $\dot{\mathcal{M}}$, based on the standard disc model (Shakura & Sunyaev 1973) shown in Eq. (4) of Du et al. (2018). We calculate the former to easily facilitate comparison with other literature, while we calculate the latter to compare with the super-Eddington accreting massive black hole (SEAMBH) sample of Du et al. (2018). We use the best-fit inclination angle of each object to calculate $\dot{\mathcal{M}}$. However, due to the presence of outflowing radial motion in Mrk 509 and PDS 456, their best-fit inclination angles may not be easily as straightforward as the rest of the GRAVITY-observed AGNs, and therefore, the values of their $\dot{\mathcal{M}}$ should be taken with caution.

For Mrk 509, IC 4329A, NGC 3783, and IRAS 09149-6206, we use the 14–195 keV luminosity from the BAT AGN Spectroscopic Survey (Koss et al. 2017), which is directly taken from the 70-month Swift-BAT survey (Baumgartner et al. 2013), to calculate L_bol through the calibration in (Winter et al. 2012). For 3C 273 and PDS 456, we decided to use the bolometric luminosity calculated from the 5100 Å luminosity and the bolometric correction of Trakhtenbrot et al. (2017), which was also used by GRAVITY Collaboration (2020c). This is because 3C 273 is most likely dominated by the jet in the X-ray regime (Dermer et al. 1997; Vasudevan & Fabian 2007), and PDS 456 exhibits strong X-ray variability, causing the observed L_{2 − 10 keV} to be 0.2% of its bolometric luminosity (Reeves et al. 2009). We also use the 5100 Å luminosity and bolometric correction of Trakhtenbrot et al. (2017) to get the bolometric luminosity of Mrk 1239 from the L_λ(5100 Å) measurement of Pan et al. (2021). This bolometric luminosity is comparable to (only slightly larger than) the bolometric luminosity estimated from WISE W3 (12 μm). Both of these L_bol values are much larger than the bolometric luminosity estimated from the 2–10 keV X-ray luminosity (Jiang et al. 2021), which is likely because of the high obscuration that is not accounted for.

Table 5 lists the calculated λ_Edd and $\dot{\mathcal{M}}$ of our 7 targets based on GRAVITY-derived properties, together with the L_bol used to calculate λ_Edd in this work. The λ_Edd values calculated from GRAVITY-derived values are consistent with those taken from the literature, except PDS 456. As for $\dot{\mathcal{M}}$, we see that the range of the $\dot{\mathcal{M}}$ of our targets are similar to that of Du et al. (2018) (about 10⁻³–10³). PDS 456 shows a relatively large Eddington ratio compared to the previous estimate assuming a bolometric luminosity of ∼10⁴⁷ erg s⁻¹ (Reeves et al. 2000) and a BH mass of log M_BH/M_⊙ = 9.24 (Nardini et al. 2015). However, Yang et al. (2021) also find evidence that the accretion rate of PDS 456 may exceed the Eddington rate based on the observed X-ray wind velocities on PDS 456.

Our Eddington ratio calculations suggest that the BLR sizes of our targets are consistent with the picture that the accretion rate plays a role in the observed deviation from the canonical R-L relation, as shown by Du et al. (2018) and Du & Wang (2019). These previous works have shown that ℛ(Fe II), which is thought to be an indicator of accretion rate, is correlated with the deviation of AGNs from the canonical R-L relation. However, recent results from Woo et al. (2024) undermine these conclusions, as they found the aforementioned correlation to be weaker after adding more luminous sub-Eddington AGNs. Nevertheless, our sample size is limited, and therefore, more AGN observations with GRAVITY will help provide independent evidence to further test this hypothesis.

8.2. Effect of used luminosity in the R-L relation

We also investigate what happens to the R-L relation if we use a different indicator of AGN luminosity other than λL_λ(5100 Å), the standard luminosity for reporting the R-L relation. It is possible λL_λ(5100 Å) is not a good measure of a target’s ionising luminosity at high luminosities if, for example, the SED shape changes significantly. In this case, the bolometric luminosity (measured from X-ray observations) might be considered. If we use the log L_bol listed in Table 5 to plot the R-L relation, we get a best-fit slope of $\alpha =0.{49}_{-0.30}^{+0.30}$ (in our discussion, we disregard the effect of relative time lags with respect to Hβ for simplicity). This suggests that the apparent decrease in the slope of the R-L relation might instead be caused by a non-linear UV-to-optical relationship. Similarly, Netzer (2019) presents a luminosity-dependent bolometric correction. Using this to convert our log λL_λ(5100 Å) values to L_bol, we find the slope of the R-L_AGN relation to be 0.48${}_{-0.21}^{+0.22}$. Therefore, while we are limited to our small sample size, our large luminosity range has revealed that a changing SED shape could instead explain the apparent departure of the R-L relation from R ∼ L^0.5. This will be further assessed with a much larger sample of GRAVITY-observed AGNs.

Woo et al. (2024) also recently investigated the possibility of log λL_λ(5100 Å) not being a good representative of the AGN luminosity of their targets. Their tests show the systematic change of the SED slope between the UV and optical wavelengths may be partly responsible for the deviation from the 0.5 slope of the R-L relation, similar to what our test above suggests.

9. The black hole mass-stellar velocity dispersion relation

Since one of the ultimate goals of our work is to acquire precise BH masses, we are also interested in where our targets lie in the M_BH-σ_* relation. Figure 9 compares our data points with M_BH-σ_* relations from other studies. We focus our comparison with Greene et al. (2020), wherein they fit the M_BH-σ_* relations for early- and late-type galaxies. This sample comprises various subsamples of galaxies from different works (e.g., Kormendy & Ho 2013; Greene et al. 2016; Saglia et al. 2016; Krajnović et al. 2018; Thater et al. 2019; den Brok et al. 2015; Nguyen et al. 2018, 2019). Upper limits of galaxies with smaller stellar velocity dispersions were estimated via stellar dynamics (see Table 4 of Greene et al. 2020 and references therein). The stellar velocity dispersions of our targets are drawn from the literature. For Mrk 509, Mrk 1239, and IC 4329A, we adopt the measurements from Oliva et al. (1999) derived from Si 1.59 μm, CO (6,3) 1.62 μm, and CO(2,0) 2.29 μm. We use the average σ_* if there are two or more available measurements for each object, and the typical errors of their measured σ_* are about ±20 km s⁻¹. For PDS 456, we use the calculated dynamical mass from Bischetti et al. (2019) to derive its tentative stellar velocity dispersion. They acquired CO(3-2) emission maps of PDS 456, from which they measured the dynamical mass of PDS 456 within 1.3 kpc. The inferred dynamical mass is 1.0 × 10¹⁰ M_⊙. We can estimate the stellar velocity dispersion via virial theorem: ${\sigma }_{\ast }=\sqrt{\frac{\mathit{GM}}{R}}$, where G is the gravitational constant, M is the dynamical mass measured at R = 1.3 kpc. With this method, we estimate σ_* ∼ 182 km s⁻¹ for PDS 456. However, we caution readers that this estimate is uncertain, as an unknown geometric correction factor C is usually added (especially if the object is not a disc; Binney & Tremaine 2008). In a variety of galactic mass distributions, C = 6.7 could be an appropriate value to assume (e.g., Förster-Schreiber et al. 2009). In addition, the most accurate way of measuring the stellar velocity dispersion of a target is through measuring equivalent widths of stellar absorption lines (e.g., Oliva et al. 1999). However, this is yet to be done for PDS 456. We also caution that the σ_* of IRAS 09149-6206 used in this work is based on the [OIII] line and therefore is very uncertain (GRAVITY Collaboration 2020a).

Fig. 9. Plot showing BH mass vs stellar velocity dispersion (MBH-σ* relation). The symbols of our 7 AGNs are similar to that of Fig. 8. The figure is redrawn from Fig. 3 of Greene et al. (2020), but with the addition of the 4 targets introduced in this work. Early- and late-type galaxies are shown as red circles and blue hollow squares, respectively, while BH mass upper limits of some low-mass late-type galaxies are shown as inverted blue triangles. All these galaxies’ measurements are taken from Greene et al. (2020) and referenced therein. The best-fit MBH-σ* relation for early- and late-type galaxies are shown as red and blue-shaded regions, while the best-fit relation for all galaxies is shown as a grey-shaded region. The best-fit lines are surrounded with their 1σ confidence interval, and all best-fit values are also taken from Greene et al. (2020).

All of our objects lie close to the standard relation, showing similar scatter to other samples. However, Mrk 1239 and IC 4329A are shown to be below the local M_BH-σ_* relation. As shown in Table 5, these two objects also have relatively high Eddington ratios. These results are consistent with previous works showing highly accreting BHs to be below the local BH-galaxy scaling relations (Ding et al. 2022; Zhuang & Ho 2023). However, we find 3C 273 and PDS 456, the sources with the highest Eddington ratios from our sample, to be placed consistently with the M_BH-σ_* relation of early-type galaxies and above the M_BH-σ_* relation of late-type galaxies. We note that the σ_* of PDS 456 is very uncertain, but it is not the case for 3C 273. Nevertheless, the two sources are the most luminous AGNs at low redshift, which by selection requires a high BH mass and high accretion rate.

10. Possible origin of the spatial BLR offset and its relation with AGN luminosity

Previously in Sect. 4, we find a spatial offset between the BLR and the continuum photocentre for all of our targets. We call this the “BLR offset”, and it ranges from ∼30–140 μas for our four targets. Upon investigation, we find a strong positive (Pearson correlation coefficient = 0.81) correlation (p = 0.026) between the BLR offsets (henceforth called R_off) and optical AGN luminosity of all GRAVITY-observed AGNs. We show this in Fig. 10a, where we also compare the data with R-L relation taken from Bentz et al. (2013), the fitted R-L relation in this work, and the dust continuum R-L relation derived from our GRAVITY data (Abuter et al. 2024). Our comparison with these relations suggests two things about the BLR offsets: (1) that they are all within the dust sublimation radius and (2) that they are of the same scale as BLR sizes. Due to the latter, we can rule out the BLR structure as a possible origin of the observed BLR offsets.

Fig. 10. Spatial offset between the BLR and continuum photocentres and its possible physical representation and their connection with optical AGN luminosity. (a) The BLR offset and (b) the flux ratio of two sides of the hot dust (Wr/Wl) as a function of optical AGN luminosity, assuming that the hot dust emits asymmetrically as explained in Sect. 10. The 4 targets introduced in this work are shown in circles, while the previously published AGNs (3C 273, IRAS 09149-6206, and NGC 3783) are shown in triangles. The error bars are 1σ errors calculated through the propagation of uncertainties from the centroid position of the null photocentre fitting. For comparison, the relation between hot dust size and AGN luminosity (Abuter et al. 2024) is shown by the red dashed line, while the relation between the BLR radius and AGN luminosity (the so-called R-L relation) based on Bentz et al. (2013) is shown by the blue dashed line. We also show the fitted R-L relation in this work depicted by the green dashed line. The Pearson correlation coefficients and p-values for the left and right plots are 0.81 and 0.026, and −0.15 and 0.73, respectively.

We include the BLR offsets of the 3 previously published GRAVITY AGNs in Fig. 10a. The photocentre of IRAS 09149-6206 is reported in GRAVITY Collaboration (2020a). We newly measured the photocentre offsets of 3C 273 and NGC 3783. For 3C 273, the reduction method in GRAVITY Collaboration (2018) removed the continuum phase signal, so we reduce the data again with our new method based on that of GRAVITY Collaboration (2020a). For NGC 3783, we already reduce the data with our new method, so we simply perform photocentre and BLR fitting again². We measure a new BLR offset of 22.8 ± 2.1 μas and 71.6 ± 6.2 μas for 3C 273 and NGC 3783, respectively. This does not affect the derived BLR differential phase of the two targets.

Given the luminosity dependence of the offsets, we propose that the offset between the BLR and continuum photocentres results from asymmetric K-band emission from the hot dust. This asymmetry can then be simply modelled as the hot dust having a side with a brighter flux and a side with a fainter flux. Figure 11 shows the schematic diagram for an easy visualisation of the model. This asymmetry can be produced by several factors, for example, the presence of a parcel of dusty cloud/s with significant flux located within the hot dust (similar to that of GRAVITY Collaboration 2021a) or the dusty clouds are optically thick at ∼2 μm, so one preferentially sees the side illuminated by the AGN. Some of the emitting hot dust structures of AGNs could also have irregularities or clumpy regions, or the edge of a foreground dust lane could coincide with the LOS to the nucleus (GRAVITY Collaboration 2020a). If the asymmetric K-band emission from the hot dust is caused by the coincidence of the LOS to a BLR irregularity, a connection between the BLR offsets and inclination angles might be implied. However, the Pearson correlation coefficient and p-value of the log R_BLR, off [pc] and i [°] are −0.42 and 0.35, respectively, indicating the absence of a significant correlation between the two. Even if we exclude Mrk 509 and PDS 456, the targets whose best-fit inclination angles may not represent their “true” inclination angle due to the presence of outflowing radial motion in their BLRs, the resulting Pearson correlation coefficient and p-value are −0.39 and 0.51, respectively. Therefore, the asymmetric K-band emission from the hot dust could not be due to the coincidence of the LOS to any BLR irregularities.

Fig. 11. 2D schematic diagram showing our explanation for the observed offset between the BLR and hot dust. The BLR (orange disc) is rotating, has a size of RBLR, and its photocentre is marked as a red cross and is assumed to be the same as the position of the central BH. For simplicity, the hot dust (with a radius of Rdust) is assumed to have two sides: the left side being the fainter one and the right side being the brighter one. This causes the continuum photocentre (blue dot) to shift to the brighter side of the hot dust. The observed offset between the BLR and hot dust photocentres is labelled Roff.

The next objective is to determine how bright (in terms of flux) the brighter side of the hot dust is relative to the fainter side. Following our model, we investigate this by deriving a formula based on the concept of the centre of mass that will give us the flux ratio between the two sides. We arrive at the following relation (see Appendix B for more discussion):

$$\begin{array}{c}\hfill \frac{W_{\mathrm{r}}}{W_{\mathrm{l}}}=\frac{2+\pi R_{\mathrm{ratio}}}{2-\pi R_{\mathrm{ratio}}}\end{array}$$(5)

where $R_{\mathrm{ratio}}=\frac{R_{\mathrm{off}}}{R_{\mathrm{dust}}}$, R_dust is the dust sublimation radius derived from the dust radius – luminosity relation from Abuter et al. (2024), W is the flux, and the subscripts l and r refer to the left and right sides of the hot dust, respectively.

We show W_r/W_l as a function of optical AGN luminosity in Fig. 10b. The average value of W_r/W_l is 1.29 ± 0.01, meaning that the brighter side of the hot dust has a ∼30% higher flux than its fainter side. The Pearson correlation p-value of W_r/W_l vs. log λL_λ(5100 Å) is ∼0.73, suggesting insignificant correlation between the two quantities.

The absence of a significant relationship between the BLR offsets and inclination angle also purports that the “tilted torus” model (Lawrence & Elvis 2010) could not explain the presence of BLR offsets. Other possible reasons include variable obscuration on the BLR and smaller scales. Dehghanian et al. (2019) proposed that a “holiday” period (a period when the covering factor of the LOS obscurer varies) can explain the decorrelation of the emission lines and continuum variations observed on NGC 5548. However, it would be too much of a coincidence to conclude that all the targets were observed during their “holiday” periods so that their BLR photocentres are offset with respect to their continuum photocentres.

11. Virial factors

Finally, we calculate the virial factors of our 4 targets. Dalla Bontá et al. (2020) argued that the dispersion or σ (i.e., square root of the second moment of the line profile) is better than the FWHM in calculating the virial factor due to the latter introducing a bias in the BH mass scale. Yu et al. (2019) and Wang et al. (2019) also purported a similar conclusion as they found the σ is insensitive to the inclination angle, and hence, the BLR geometry doesn’t affect the resulting BH mass. Hence, we follow this prescription and measure the σ of the model line profile. Afterwards, the values of the relevant model parameters are drawn randomly from the sampled posterior parameter space created during BLR model fitting. We use the formula $f_{\sigma }=G{M}_{\text{BH}}/(R_{\text{BLR}}{v}_{\sigma }^2)$ to calculate the virial factor.

Figure 12 shows the virial factors of our sources as a function of optical AGN luminosity. The virial factors of Mrk 509, PDS 456, Mrk 1239, and IC 4329A are 1.10${}_{-0.18}^{+0.16}$, 1.17${}_{-0.38}^{+0.27}$, 6.39${}_{-3.37}^{+3.32}$, and 2.98${}_{-1.53}^{+1.23}$, respectively. The virial factor of 3C 273 was taken from GRAVITY Collaboration (2018), while the virial factors of NGC 3783 and IRAS 09149-6206 are calculated from the posterior distributions of their parameters. The average virial factor of all 7 AGNs is ⟨f_σ⟩ = 3.04 ± 0.64. This is lower than that calculated by previous works that assumed their AGN sample follows the same M_BH-σ_* as those of quiescent galaxies (Onken et al. 2004; Woo et al. 2010; Park et al. 2012; Batiste et al. 2017), which is ⟨f_σ⟩∼5. It is not clear whether we should expect a match here, since the host galaxy type plays a role in the determination of the M_BH-σ_* relation. And, our average virial factor is consistent with that of Grier et al. (2017) and Williams et al. (2018), who also calculated the individual virial factors of their AGN sample and got ⟨f_σ⟩∼3. We also note that we find lower f_σ values for Mrk 509 and PDS 456, which intriguingly are the two targets with significant radial motion in their BLR.

Fig. 12. Virial factor (fσ) as a function of optical AGN luminosity. The legends are similar to that of Fig. 10. The Pearson correlation p-value is ∼1.0, indicating an insignificant correlation between the two quantities.

The resulting Pearson correlation coefficient between log f_σ and log λL_λ(5100 Å) is ∼ − 0.02, and the p-value is ∼0.96, indicating insignificant correlation. Previous works (Villafaña et al. 2023; Williams et al. 2018) reported an positive correlation, although non-significant (p-values ≫ 0.05), between the two quantities. They also reported a possible positive correlation between log f_σ and log M_BH/M_⊙. However, calculating the Pearson correlation coefficient and p-value of these two quantities in our work suggests otherwise (−0.27 and 0.56, respectively). Hence, we conclude that there is no correlation between the virial factor and optical luminosity or BH mass for our targets. However, the uncertainties are large, and bigger samples will be needed to confirm this result.

12. Conclusions and future prospects

We investigate the broad Brγ line-emitting regions of Mrk 509, Mrk 1239, and IC 4329A, and the broad Paα line-emitting region of PDS 456. To study the kinematics and properties of their BLRs, we performed photocentre and BLR model fitting. Our results support many of the assumptions and scaling relations used and derived by RM, and they provide an independent method that can be extended to high-z and large BH mass:

1. Most of the AGNs observed by GRAVITY can be well described by a thick, rotating disc of clouds, consistent with RM assumptions. However, two of our targets (Mrk 509 and PDS 456) show evidence of outflowing BLR clouds, while the other two targets (Mrk 1239 and IC 4329A) have relatively weaker differential phase signals that limit the constraints of their BLR sizes and BH masses. Nevertheless, we prove from our analyses that the differential phase is crucial in providing stronger constraints on the BLR kinematics.

2. By adding the other three AGNs that were previously observed with GRAVITY (3C 273, IRAS 09149-6206, and NGC 3783), we derived a new R-L relation based on GRAVITY-derived BLR sizes only. We derived a slope and intercept of $\alpha =0.{37}_{-0.17}^{+0.18}$ and $K=1.{69}_{-0.23}^{+0.23}$, respectively. Our results are consistent with works showing shallower R-L relation slopes, although the effect of the Eddington ratio in this relation cannot be fully realised with our relatively small sample.

3. Most of the GRAVITY AGNs are consistent within the scatter of the standard M_BH-σ_* relation. Mrk 1239 and IC 4329A, two of our targets with relatively high Eddington ratios, are placed below the M_BH-σ_* relation consistent with previous works, while 3C 273 and PDS 456, the targets possessing the highest Eddington ratios among our sources, are still consistent within the M_BH-σ_* relation, suggesting that their high luminosity may play a role in their difference compared to other high Eddington-ratio sources such as Mrk 1239 and IC 4329A.

4. We find a significant correlation between the offset between the BLR and continuum photocentre and optical AGN luminosity, and the offsets of all GRAVITY AGNs are found to be of a similar scale as their BLR sizes. This correlation is surmised to be due to asymmetric K-band emission of the hot dust, and a simple model was created to estimate that this emission is ∼30% brighter on one side than the other for our sample.

5. Lastly, we calculated the virial factors of our four targets and found an average virial factor (based on the dispersion of the model line profiles) of ⟨f_σ⟩ = 3.04 ± 0.64 for all GRAVITY-observed AGNs. This is consistent with previous works that calculated the individual virial factors of their targets. Further increasing the observed AGNs with GRAVITY will greatly help in providing better hints in answering some of the remaining questions unanswered by this work. With the advent of GRAVITY+, the limits of GRAVITY will be pushed even further, providing wide-field off-axis fringe-tracking, newer adaptive optics systems on all the UTs, and laser guide stars. With GRAVITY+, observing hundreds of AGNs and even those at high redshift will be achievable (GRAVITY+ Collaboration 2022). An example of this feat is the first dynamical mass measurement of a z ∼ 2 quasar with GRAVITY-Wide (Abuter et al. 2024). It is expected, therefore, that we may also be able to shed light on the R-L relation at high redshift and provide more accurate BH mass measurements for galaxies at the cosmic noon.

Acknowledgments

This research has used the NASA/IPAC Extragalactic Database (NED), operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration. This research has used the SIMBAD database, operated at CDS, Strasbourg, France.

References

Appendix A: Corner plots of BLR model fits

The posterior probability density distributions of the fitted parameters for all of our targets are shown in Figs. A.1 to A.4 as corner plots.

Fig. A.1. Posterior probability distribution of the fitted parameters from the elliptical and radial model fitting Mrk 509 data. The blue lines and cross points refer to the maximum a posteriori value. The vertical dashed lines represent the 68% (1σ) credible intervals.

Fig. A.2. Similar to Fig. A.1 but for PDS 456.

Fig. A.3. Posterior probability distribution of the fitted parameters from the circular model fitting Mrk1239 data. The legends from Fig. A.1 still apply.

Fig. A.4. Similar to Fig. A.3 but for IC 4329A.

Appendix B: Derivation of the flux ratio of two sides of the hot dust

First, let P_cont be the distance of the hot dust photocentre to the origin:

$$\begin{array}{c}\hfill P_{\mathrm{cont}}=\frac{\mathrm{\Sigma }{w}_i{r}_i}{\mathrm{\Sigma }{w}_i}\end{array}$$(B.1)

where w_i is the weight/flux of a dust clump at r_i. We define P_BLR as the photocentre position of the BLR. In our simple model, the central BH is located at the origin of our coordinate system, and P_BLR is also situated at the origin. Considering the left and right side of the BLR, we have

$$\begin{array}{c}\hfill P_{\mathrm{cont}}^{\prime }=\frac{\mathrm{\Sigma }{w}_l{r}_l+\mathrm{\Sigma }{w}_r{r}_r}{\mathrm{\Sigma }{w}_l+\mathrm{\Sigma }{w}_r}\end{array}$$(B.2)

where the subscripts l and r refer to the left and right sides of the hot dust, respectively. The summations must be converted into definite integrals with the range of angles encompassing the left and right sides as their respective limits:

$$\begin{array}{c}\hfill \begin{array}{c}{P}_{\mathrm{cont}}=\frac{1}{{\int }_{\pi /2}^{3\pi /2}{w}_l\mathrm{d}\theta +{\int }_{-\pi /2}^{\pi /2}{w}_r\mathrm{d}\theta }(W_{\mathrm{l}}{R}_{\mathrm{dust}}{\int }_{\pi /2}^{3\pi /2}cos\theta \mathrm{d}\theta \hfill \\ +W_{\mathrm{r}}{R}_{\mathrm{dust}}{\int }_{-\pi /2}^{\pi /2}cos\theta \mathrm{d}\theta )\hfill \end{array}\end{array}$$(B.3)

where W_l and W_r are the total flux of the left and right sides of the hot dust, respectively. R_dust is the average size of the hot dust (we disregard the effect of a possible large discrepancy between the inner and outer radius of the hot dust ring for simplicity). The equation simplifies to:

$$\begin{array}{c}\hfill P_{\mathrm{cont}}=\frac{2R_{\mathrm{dust}}(W_r-W_l)}{\pi (W_l+W_r)}.\end{array}$$(B.4)

If we re-define R_off = P_cont (since P_BLR coincides with the origin and the continuum photocentre simply represents the actual offset), we get $\frac{R_{\mathrm{off}}}{R_{\mathrm{dust}}}$ = $\frac{2}{\pi }\frac{W_{\mathrm{r}}-W_{\mathrm{l}}}{W_{\mathrm{r}}+W_{\mathrm{l}}}$. To calculate W_r/W_l, we have:

$$\begin{array}{c}\hfill \frac{W_r}{W_l}=\frac{2+\pi R_{\mathrm{ratio}}}{2-\pi R_{\mathrm{ratio}}}\end{array}$$(B.5)

where $R_{\mathrm{ratio}}=\frac{R_{\mathrm{off}}}{R_{\mathrm{dust}}}$.

All Tables

All Figures

Fig. 1. Averaged differential phase spectra (coloured, based on their respective baselines labelled in the top row) and the normalised Brγ (for Mrk 509, Mrk 1239, and IC 4329A) and Paα (for PDS 456) spectra (black histogram) of Mrk 509 (first row), PDS 456 (second row), Mrk 1239 (third row), and IC 4329A (fourth row). The curves show negative signals in Mrk 509 and PDS 456 with peaks ranging from ∼ − 0.4° to −0.5°, and positive signals in Mrk 1239 and IC 4329A with peaks ranging from ∼0.5°–0.6°.

In the text

	Fig. 2. Average AGN flux spectra (black steps with grey error bars) of (a) Mrk 509 and (b) PDS 456 taken by GRAVITY, (c) Mrk 1239 taken by APO/TRIPLESPEC, and (d) IC 4329A taken by GRAVITY, all normalised to the continuum. The name of the emission lines for each target is shown in each panel.
In the text

Fig. 3. Best-fit BLR photocentres of our four targets. The columns from left to right show the photocentres from the data, the “2-pole” and “null” photocentre models (described in Sect. 4), and the photocentres from the best-fit BLR model of each target. The red cross refers to the continuum photocentre. The direction of the radio jet is shown as a green solid line. The black arrow shows the position of the best-fit “null” model fitted from the best-fit BLR model. The ellipses around each centroid refer to 68% (1σ) confidence intervals. The colours refer to the corresponding velocity of each spectral channel. The (significant) separation of the red- and blueshifted photocentres confirm that we have resolved the BLRs of our targets. For Mrk 1239 and IC 4329A, the perpendicular alignment between the radio jet and red-blue photo centres from the “2-pole” model fitting of the data and the straight alignments of their model photo centres indicate that their BLR is rotation-dominated. On the other hand, for Mrk 509 and PDS 456, the radio jet and red-blue photocentres are not closely perpendicular, and the model photocentres are curved, which suggests that their BLRs are radial motion-dominated.

In the text

Fig. 4. Posterior distributions of BLR radius (left column) and BH mass (right column) resulting from BLR fitting without differential phase data (blue histogram) and with differential phase data (orange histogram). Each row corresponds to each target; from top to bottom: Mrk 509, PDS 456, Mrk 1239, and IC 4329A. The histograms are normalised such that they integrate into 1. It is clear that the spatial information from the differential phase significantly constrains both the BLR radius and BH mass.

In the text

Fig. 5. Expected differential phase signals for Mrk 1239 (left panel) and IC 4329A (right panel) when their line profiles were fitted with our circular model assuming their maximum (red shaded region) and minimum (blue shaded region) BH masses were fixed during fitting. The assumed BH masses for Mrk 1239 and IC 4329A are taken from previous literature as discussed in the text. The shaded regions refer to 1σ uncertainty. The chosen baselines to average are the same as the ones used to create Fig. 6: UT4-UT3, UT4-UT2, and UT4-UT1 for Mrk 1239; UT4-UT2, UT4-UT1, and UT3-UT2 for IC 4329A. The black points refer to the maximum and minimum observed phase signals at the line region of Mrk 1239 and IC 4329A after averaging the differential phase spectra from the same chosen baselines. The horizontal lines show the wavelength range of the maximum and minimum observed phase signals.

In the text

Fig. 6. Summary of BLR fitting results for (a) Mrk 509, (b) PDS 456, (c) Mrk 1239, and (d) IC 4329A. The observed averaged differential phase spectra with the continuum phase signal removed are shown in blue, while the median differential phase spectra from the best-fit BLR model are shown in red (see text for the details of averaging the differential phase spectra for each target). Above each panel, the observed line profile (black) is shown together with the median best-fit model line profile (red). The central wavelengths of the line profiles are shown by the grey vertical dashed lines. The red-shaded region shows the 16-to-84 percentiles of the averaged differential phase spectra and flux spectra from the best-fit model, which was created by getting randomly selected samples of each parameter in 100 different instances and then recreating the averaged differential phase spectra. Mrk 509 and PDS 456 clearly show an asymmetrical signal typical for the Pancoast model, while Mrk 1239 and IC 4329A show the S-shape signal typical for the Keplerian model.

In the text

Fig. 7. Edge-on views of the best-fit outflow model for (a) Mrk 509 and (b) PDS 456. For clarity, the PA for both configurations is adjusted to 180°, and the BLR centre is positioned to the origin. The LOS at the best-fit inclination angle is shown as a black arrow. The midplane, represented as a green line, is defined as the line passing through the origin and is perpendicular to the LOS at i = 0°. The regions above (closer to the observer) and below (farther to the observer) the midplane are also labelled. The colours of each cloud refer to their LOS velocity, while their sizes reflect the weight given to each cloud towards the total emission; the larger the size of the cloud, the more this cloud contributes to the broad-line emission. The edge-on view of the BLRs of Mrk 509 shows that its asymmetry is due to a highly obscuring midplane; for PDS 456, a preference for broad-line emission to come from the far side of the BLR.

In the text

Fig. 8. Logarithm of the BLR radius as a function of log λLλ(5100 Å) (R-L relation). The blue line shows the best-fit line derived from all 7 GRAVITY-observed AGNs, and the blue shaded region is its 1σ confidence interval. The AGN luminosity is taken from Table 1. The dashed line represents the best-fit radius luminosity relation from Bentz et al. (2013). The 3 published AGNs (IRAS 09149-6206 from GRAVITY Collaboration 2020a, 3C 273 from GRAVITY Collaboration 2018, and NGC 3783 from GRAVITY Collaboration 2021a) are shown as a red star, red diamond, and red square, respectively. Our 4 targets are shown in differently coloured filled squares with 1σ errors. Our best-fit R-L (power-law) relation has a slope of $\alpha =0.{37}_{-0.17}^{+0.18}$, which is lower but, given the current large uncertainty, still consistent within uncertainties with the slope of R-L relation fit from Bentz et al. (2013; $\alpha =0.{533}_{-0.033}^{+0.035}$). If we follow the prescription from photoionisation models, the BLR sizes of 3C 273 and PDS 456 will be adjusted (shown as open symbols), and the fitted GRAVITY R-L relation is shown as the green line.

In the text

Fig. 9. Plot showing BH mass vs stellar velocity dispersion (MBH-σ* relation). The symbols of our 7 AGNs are similar to that of Fig. 8. The figure is redrawn from Fig. 3 of Greene et al. (2020), but with the addition of the 4 targets introduced in this work. Early- and late-type galaxies are shown as red circles and blue hollow squares, respectively, while BH mass upper limits of some low-mass late-type galaxies are shown as inverted blue triangles. All these galaxies’ measurements are taken from Greene et al. (2020) and referenced therein. The best-fit MBH-σ* relation for early- and late-type galaxies are shown as red and blue-shaded regions, while the best-fit relation for all galaxies is shown as a grey-shaded region. The best-fit lines are surrounded with their 1σ confidence interval, and all best-fit values are also taken from Greene et al. (2020).

In the text

Fig. 10. Spatial offset between the BLR and continuum photocentres and its possible physical representation and their connection with optical AGN luminosity. (a) The BLR offset and (b) the flux ratio of two sides of the hot dust (Wr/Wl) as a function of optical AGN luminosity, assuming that the hot dust emits asymmetrically as explained in Sect. 10. The 4 targets introduced in this work are shown in circles, while the previously published AGNs (3C 273, IRAS 09149-6206, and NGC 3783) are shown in triangles. The error bars are 1σ errors calculated through the propagation of uncertainties from the centroid position of the null photocentre fitting. For comparison, the relation between hot dust size and AGN luminosity (Abuter et al. 2024) is shown by the red dashed line, while the relation between the BLR radius and AGN luminosity (the so-called R-L relation) based on Bentz et al. (2013) is shown by the blue dashed line. We also show the fitted R-L relation in this work depicted by the green dashed line. The Pearson correlation coefficients and p-values for the left and right plots are 0.81 and 0.026, and −0.15 and 0.73, respectively.

In the text

Fig. 11. 2D schematic diagram showing our explanation for the observed offset between the BLR and hot dust. The BLR (orange disc) is rotating, has a size of RBLR, and its photocentre is marked as a red cross and is assumed to be the same as the position of the central BH. For simplicity, the hot dust (with a radius of Rdust) is assumed to have two sides: the left side being the fainter one and the right side being the brighter one. This causes the continuum photocentre (blue dot) to shift to the brighter side of the hot dust. The observed offset between the BLR and hot dust photocentres is labelled Roff.

In the text

	Fig. 12. Virial factor (fσ) as a function of optical AGN luminosity. The legends are similar to that of Fig. 10. The Pearson correlation p-value is ∼1.0, indicating an insignificant correlation between the two quantities.
In the text

	Fig. A.1. Posterior probability distribution of the fitted parameters from the elliptical and radial model fitting Mrk 509 data. The blue lines and cross points refer to the maximum a posteriori value. The vertical dashed lines represent the 68% (1σ) credible intervals.
In the text

	Fig. A.2. Similar to Fig. A.1 but for PDS 456.
In the text

	Fig. A.3. Posterior probability distribution of the fitted parameters from the circular model fitting Mrk1239 data. The legends from Fig. A.1 still apply.
In the text

	Fig. A.4. Similar to Fig. A.3 but for IC 4329A.
In the text