# The spatially resolved broad line region of IRAS 09149-6206

A. Amorim, M. Bauböck, W. Brandner, Y. Clénet, R. Davies, P. de Zeeuw, J. Dexter, A. Eckart, F. Eisenhauer, N. Förster Schreiber, et al.

► To cite this version:

A. Amorim, M. Bauböck, W. Brandner, Y. Clénet, R. Davies, et al. The spatially resolved broad line region of IRAS 09149-6206. Astronomy and Astrophysics - A&A, 2020, 643, pp.A154. 10.1051/0004-6361/202039067. hal-03970181

HAL Id: hal-03970181

https://hal.science/hal-03970181

Submitted on 3 Feb 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# The spatially resolved broad line region of IRAS 09149-6206

GRAVITY Collaboration: A. Amorim$^{19,21}$, M. Bauböck$^1$, W. Brandner$^{22}$, Y. Clénet$^2$, R. Davies$^1$, P. T. de Zeeuw$^{1,17}$, J. Dexter$^{24,1}$, A. Eckart$^{3,18}$, F. Eisenhauer$^1$, N. M. Förster Schreiber$^1$, F. Gao$^1$, P. J. V. Garcia$^{15,20,21}$, R. Genzel$^{1,4}$, S. Gillessen$^1$, D. Gratadour$^{2,25}$, S. Hönig$^5$, M. Kishimoto$^6$, S. Lacour$^{2,16}$, D. Lutz$^1$, F. Millour$^7$, H. Netzer$^8$, T. Ott$^1$, T. Paumard$^2$, K. Perraut$^{12}$, G. Perrin$^2$, B. M. Peterson$^{9,10,11}$, P. O. Petrucci$^{12}$, O. Pfuhl$^{16}$, M. A. Prieto$^{23}$, D. Rouan$^2$, J. Shangguan$^{1,*}$, T. Shimizu$^1$, M. Schartmann$^1$, J. Stadler$^1$, A. Sternberg$^{8,14}$, O. Straub$^1$, C. Straubmeier$^3$, E. Sturm$^1$, L. J. Tacconi$^1$, K. R. W. Tristram$^{15}$, P. Vermot$^2$, S. von Fellenberg$^1$, I. Waisberg$^{13}$, F. Widmann$^1$, and J. Woillez$^{16}$

*(Affiliations can be found after the references)*

Received 30 July 2020 / Accepted 16 September 2020

## ABSTRACT

We present new near-infrared VLTI/GRAVITY interferometric spectra that spatially resolve the broad Bry emission line in the nucleus of the active galaxy IRAS 09149-6206. We use these data to measure the size of the broad line region (BLR) and estimate the mass of the central black hole. Using an improved phase calibration method that reduces the differential phase uncertainty to 0.05° per baseline across the spectrum, we detect a differential phase signal that reaches a maximum of ~0.5° between the line and continuum. This represents an offset of ~120 µas (0.14 pc) between the BLR and the centroid of the hot dust distribution traced by the 2.3 µm continuum. The offset is well within the dust sublimation region, which matches the measured ~0.6 mas (0.7 pc) diameter of the continuum. A clear velocity gradient, almost perpendicular to the offset, is traced by the reconstructed photocentres of the spectral channels of the Bry line. We infer the radius of the BLR to be ~65 µas (0.075 pc), which is consistent with the radius-luminosity relation of nearby active galactic nuclei derived based on the time lag of the Hβ line from reverberation mapping campaigns. Our dynamical modelling indicates the black hole mass is ~1 × 10$^8$ M$_{\odot}$, which is a little below, but consistent with, the standard M$_{BH}$-σ$_{*}$ relation.

**Key words.** galaxies: active – galaxies: nuclei – galaxies: Seyfert – quasars: supermassive black holes – techniques: interferometric

## 1. Introduction

Massive black holes in the centres of galaxies are a key component of galaxy evolution because of the role that accreting black holes have in the feedback that regulates star formation and galaxy growth (Booth & Schaye 2009; Fabian 2012; Somerville & Davé 2015; Dubois et al. 2016). Knowing their masses is central to our understanding of this co-evolution (Hopkins et al. 2008; Somerville et al. 2008; Kormendy & Ho 2013; Heckman & Best 2014), but measuring their masses robustly is challenging because it requires resolving spatial scales where the black hole dominates the gravitational potential (Ferrarese & Ford 2005). For inactive galaxies, it can be done using spatially resolved stellar (e.g. Thomas et al. 2004; Onken et al. 2014; Saglia et al. 2016) or gas kinematics (Hicks & Malkan 2008; Davis 2014; Onishi et al. 2017; Boizelle et al. 2019). These methods can be difficult to apply to active galaxies displaying broad emission lines because the active galactic nucleus (AGN) itself is so bright, although it has been possible for a few objects (Davies et al. 2006). Instead, the most precise method for measuring black hole masses in AGNs is through megamaser kinematics using very-long-baseline interferometry (Greene et al. 2010; Kuo et al. 2011). However, this requires the megamasers to be orbiting in a Keplerian disk restricted to near-edge-on view, which is very rare (Zhu et al. 2011; van den Bosch et al. 2016). Reverberation mapping (RM) exploits the variability of the black hole accretion to shed light on the size of the broad line region (BLR) and hence leads to a measurement of the black hole mass (Blandford & McKee 1982; Peterson 1993; Peterson

et al. 2004). Despite potential biases (Shankar et al. 2019) and caveats in terms of the virial factor *f* that reflects the unknown geometry and kinematics of the BLR (Onken et al. 2004; Woo et al. 2010; Ho & Kim 2014), this factor has been used successfully for many years.

Spatially resolving the very small size of the BLR, 10$^3$–10$^5$ gravitational radii (Netzer & Laor 1993; Peterson 1993; Netzer 2015), has been a long-standing goal of spectroastrometry (Petrov et al. 2001; Marconi et al. 2003), which has recently become possible with long baseline near-infrared interferometry. GRAVITY, a second-generation Very Large Telescope Interferometer (VLTI) instrument, has greatly improved the sensitivity of earlier efforts and has been able to combine all four of the 8 m Unit Telescope (UT) beams to yield six simultaneous baselines (GRAVITY Collaboration 2017). In GRAVITY Collaboration (2018, hereafter, GC18), we reported the first robust measurements of BLR size and kinematics for 3C 273 by combining differential phase spectra with the Paα emission line profile to model the BLR as a thick rotating disk under the gravitational influence of a black hole of ~3 × 10$^8$ M$_{\odot}$. This value is fully consistent with the result of a subsequent study using 10 year RM data (Zhang et al. 2019). We have now embarked on a programme to make independent estimates of the black hole masses in a sample of AGNs. The aim is not just to verify the masses derived through RM, but to better understand the structure of the BLR, which has implications for inflow (accretion) and outflow processes on small scales.

The classical picture of the BLR is a virialised distribution of clouds, with good evidence that many are rotating systems. This comes from a variety of observations, including variations in the polarisation properties across the broad line

<footnote>* Corresponding author: e-mail: shangguan@mpe.mpg.de</footnote>
profile (Smith et al. 2004, 2005). In addition, a significant minority, ~3%, of BLRs in both radio quiet and radio loud AGNs show broad double-peaked profiles that are well fitted by disk emission (Eracleous & Halpern 1994, 2003; Strateva et al. 2003; Storchi-Bergmann et al. 2017). That relatively few objects show such characteristics may be explained by the apparent relation between line width and shape, measured as $FWHM/\sigma$, reported by Kollatschny & Zetzl (2011). This suggests that rapidly rotating BLRs are flattened systems, while slower rotating BLRs have a more spherical structure due to turbulence. The turbulent velocities of 300–700 km s$^{-1}$ for Hβ and 2000–4000 km s$^{-1}$ for C IV may be indicative of outflowing gas. Disk winds are theoretically appealing and apply in many different situations, including AGNs (Emmering et al. 1992; Murray et al. 1995; Dehghanian et al. 2020); an outflow at a characteristic elevation of 30° above the mid-plane is the basis of the concept proposed by Elvis (2000), a geometry that can explain a variety of observations, including the broad absorption and broad emission lines. Rotating disk winds have been developed on a more physical basis by Everett (2005) and Keating et al. (2012); PG 1700+518 is one example of a rotating outflow that has been observed (Young et al. 2007). It therefore seems likely that the dynamics of the BLR may be a combination of rotation and outflow in a way that depends on the properties of individual objects.

Recently, thanks to high signal-to-noise ratio (S/N) spectra, the non-parametric velocity-resolved RM analyses of tens of bright AGNs have begun to reveal this structure, through qualitative evidence for both inflow and outflow in addition to virial motion (e.g. Denney et al. 2009; Bentz et al. 2010; Peterson 2014; Du et al. 2016a, 2018a; De Rosa et al. 2018; Horne et al. 2020). Furthermore, temporal variation of the BLR geometry and dynamics has been reported for some AGNs (Pei et al. 2017; Pancoast et al. 2018; Xiao et al. 2018), implying that the virial factor of the same source may evolve with time.

Parametric models of the BLR geometry and dynamics have been successfully applied to fewer than two dozen AGNs that have both high S/N spectra and high cadence monitoring (Pancoast et al. 2014a; Grier et al. 2017a; Williams et al. 2018; Li & Wang 2018). These enable one to fit for radial (inflow or outflow) motion of the clouds in addition to rotation, and hence also to derive the virial factors for individual objects. While the number of objects is too small for robust statistics, more than half of the BLRs are dominated by radial motion, with inflow and outflow being equally common. Remarkably, the black hole masses inferred from such dynamical modelling are usually consistent with those measured with the classical RM method, even when the BLR is dominated by outflow (Williams et al. 2018). However, there are degeneracies in the models (Grier et al. 2017a); and the impact of systematic uncertainties, such as the assumption that the ionisation source is point-like and that the BLR structure does not change significantly during the monitoring campaign, are challenging to assess (Pancoast et al. 2014b; Grier et al. 2017a). In particular, assumptions about the geometry of different BLR components may significantly bias the inferred physical interpretation (Mangham et al. 2019). As such, an independent method to constrain the BLR structure is much needed.

This role is fulfilled by optical and near-infrared long baseline interferometry combined with spectroastrometry. In this paper we present an analysis of new GRAVITY observations for IRAS 09149–6206 ($\alpha$ = 09:16:09.39, $\delta$ = -62:19:29.9). Perez et al. (1989) serendipitously discovered IRAS 09149–6206 as an AGN in the IRAS Point Source Catalogue during a search for planetary nebulae after optical spectra

showed characteristic broad Balmer lines. IRAS 09149–6206 is at a redshift of 0.0573 (Perez et al. 1989), which means that the broad Bry line can still be observed in the K-band and thus GRAVITY. Unfortunately, very little archival data exist for IRAS 09149–6206. It is detected in early radio continuum surveys, but not resolved (Murphy et al. 2007; Panessa et al. 2015). Existing optical/NIR images (e.g. Veron-Cetty & Woltjer 1990; Márquez et al. 1999) cannot constrain the morphological type of the host galaxy, and mid-infrared interferometric observations only barely resolve it (Kishimoto et al. 2011; Burtscher et al. 2013; López-Gonzaga et al. 2016). Section 2 describes the observations and data reduction including changes which significantly improve the precision of our differential phase spectra. In Sect. 3, we show that the model-independent photocentre positions reveal the detection of a spatial offset between the hot dust continuum and Bry, and a velocity gradient across the emission line. Building on this, in Sect. 4 we adopt a kinematic cloud model to fit the phase data in order to constrain the BLR size and kinematics. In Sect. 5, we discuss the implications of our measured $M_{BH}$ and place IRAS 09149–6206 on the BLR radius–luminosity relation. The possible interpretations of the offset between the BLR and the continuum is discussed in Sect. 6.

This work adopts the following parameters for a $\Lambda$CDM cosmology: $\Omega_m = 0.308$, $\Omega_\Lambda = 0.692$, and $H_0 = 67.8$ km s$^{-1}$ Mpc$^{-1}$ (Planck Collaboration XIII 2016). Using this cosmology, 1 pc subtends 0.87 mas on sky and $1\mu$s corresponds to 1.37 light day at the redshift of IRAS 09149–6206.

## 2. Observations and data reduction

We observed IRAS 09149–6206 with GRAVITY (GRAVITY Collaboration 2017) using all four UTs, on eight occasions between November 2018 and March 2020, primarily as part of a large programme to spatially resolve the BLR and measure black hole masses in a sample of AGNs<fnref n="1" />. Targets were selected as the brightest type 1 AGNs visible from the VLTI and above the GRAVITY sensitivity limit ($K < 11$). IRAS 09149–6206, in particular, has more than 70% of its total K band flux originating in the nucleus and is also bright and compact in the optical ($R \approx 11.5$). These properties make it an ideal GRAVITY target for observations where the source is phase referenced to itself as well as used for adaptive optics. We adopted the single-field on-axis mode for all of the observations with combined polarisation. The scientific spectra (1.95–2.45 $\mu$m) were taken in medium-resolution mode ($\lambda/\Delta\lambda \approx 500$), with 90 independent spectral elements, extracted and resampled into 210 channels.

All observations followed a similar sequence. After the telescope had pointed to the target and the adaptive optics (MACAO; Arsenault et al. 2003) had closed the loop, the telescope beams were coarsely aligned on the VLTI laboratory camera IRIS (Gitton et al. 2004). The light was then fed to GRAVITY and the internal beam tracking of GRAVITY aligned the fringe-tracking (FT) and science channel (SC) fibres on the target. After the fringe tracker had searched for and found the fringe, the scientific exposures were started. Each set of scientific exposures consisted of ten frames of 30 s integration (NDIT = 10 and DIT = 30 s). Fringe tracking is difficult for faint targets and leads to large phase noise ($\ge 0.5^\circ$). We therefore integrated deeply on-source, with only occasional sky exposures. In order to account for the observatory transfer function (e.g. coherence loss

<footnote n="1">¹ Observations were made using the ESO Telescopes at the La Silla Paranal Observatory, programme IDs 0102.B-0667 and 1103.B-0626.</footnote>
<table><caption>**Table 1.** Exposure time and weather conditions for the observations.</caption><thead><tr><th>Date</th><th>On-source time (min)</th><th>Seeing (′′)</th><th>Coherence time (ms)</th></tr></thead><tbody><tr><td>2018 Nov 20 (*)</td><td>55</td><td>0.65–0.97</td><td>3.3–4.2</td></tr><tr><td>2019 Feb 16 (†)</td><td>75</td><td>0.53–0.95</td><td>6.0–9.1</td></tr><tr><td>2019 Nov 07 (*)</td><td>85</td><td>0.36–0.71</td><td>6.5–10.7</td></tr><tr><td>2019 Nov 09 (*)</td><td>30</td><td>0.42–0.58</td><td>2.4–4.2</td></tr><tr><td>2019 Dec 09 (*)</td><td>45</td><td>0.50–0.71</td><td>1.8–3.9</td></tr><tr><td>2020 Feb 09</td><td>55</td><td>0.47–0.69</td><td>13.2–15.4</td></tr><tr><td>2020 Feb 10 (†)</td><td>75</td><td>0.51–0.99</td><td>7.7–13.2</td></tr><tr><td>2020 Mar 08 (†)</td><td>45</td><td>0.49–0.80</td><td>7.1–9.9</td></tr></tbody></table>

**Notes.** The reported data only include exposures for which we were able to track the fringe for >80% of the exposure time. The seeing and coherence time are based on measurements from the Differential Image Motion Monitor and Multi-aperture Scintillation Sensor on Paranal. (*)K stars are observed as calibrators. (†)B stars are observed as calibrators.

due to vibrations, uncorrected atmosphere, birefringence, etc.), we observed a calibrator star close to the science target. The calibrator data are used to calibrate the complex visibility and closure phase of the FT data, as well as the flux spectrum of the AGN. A calibrator was observed for all runs except 2020 Feb 09. The date, exposure time, and weather conditions for all observations are summarised in Table 1.

GRAVITY measures the complex visibilities of six baselines (telescope pairs). The visibility amplitudes measure the angular extent of a structure, whereas the phases provide its position and spatial distribution on the sky (e.g. Buscher & Longair 2015). The GRAVITY synthesis beam (about 3 mas) is usually much larger than the BLR of an AGN at cosmological distance and can only marginally resolve the continuum emission from the hot dust (GRAVITY Collaboration 2018, 2020a). Therefore, the differential phase and differential visibility amplitude are the most important interferometric observables for the BLR. Spectro-astrometry (Bailey 1998) enables us to spatially resolve the BLR, using the continuum as a reference. The differential phase measures the shift in the photocentre on sky at different wavelength channels of the broad line emission with respect to the continuum (see Appendix B). The differential visibility amplitude measures the relative size difference between the broad line emission and the continuum emission (see Sect. 4.4).

## 2.1. Pipeline data reduction

We used the latest version of the GRAVITY pipeline to reduce the data (GRAVITY Collaboration 2017; Lapeyrere et al. 2014). This includes the recent improvements to lower the imprint of the internal source of the calibration unit (Blind et al. 2014), and a better filtering of hot pixels due to cosmic rays. For the FT continuum visibility data, we used the default pipeline settings. However, the low S/N of the FT data and loss of coherence during scientific exposures often caused the pipeline to flag most of the DITs in one scientific exposure. In analysing these data, we found a substantial improvement in residual phase noise (rms scatter) by retaining all data independent of FT S/N or the estimated visibility loss<fnref n="2" /> (GC18; GRAVITY Collaboration 2020a, hereafter GC20a).

<figure><img src="image_2.png" /><figcaption>**Fig. 1.** AGN flux spectra of IRAS 09149–6206 normalised to the continuum. The black vertical bar at 2.2896 µm marks the wavelength of Bry given a redshift $z = 0.0573$ (Perez et al. 1989). We discard the wavelength ranges (≤2.10 µm and ≥2.35 µm) that are significantly affected by the atmospheric absorption.</figcaption></figure>

The pipeline reduces the data using a pixel-to-visibility matrix (Lapeyrere et al. 2014; Lacour et al. 2019). This matrix, encoding the relative throughput, coherence, phase shift, and cross-talk for each pixel, was measured during the day with the GRAVITY calibration unit after every night that GRAVITY observations were performed. Applying the matrix to the detector frames yields the instrument-calibrated complex visibilities. The pipeline then fits the phase of the science channel in each frame with a third-order polynomial to derive the phase reference, which is then subtracted from the phase. This yields the self-referenced complex visibility for every frame. Before averaging the complex visibilities of an individual exposure, we applied the algorithm developed for VLTI/AMBER (Millour et al. 2008) to calculate and remove the average phase using all of the other wavelength channels for each channel. This method produces consistent results and improves our phase errors, typically by about 10%–20%. Initial uncertainties for the combined visibilities are estimated by bootstrapping the different frames in the pipeline. However, this underestimates the uncertainty, so we adopt a different method to estimate the phase uncertainty after the pipeline data reduction.

## 2.2. Normalised profile of the broad Bry line

The emission line profile, normalised to the continuum, is necessary to derive the velocity gradient and model the dynamics of the BLR. However, it is challenging to derive the Bry line profile for IRAS 09149–6206 because the line is red-shifted to ~2.2896 µm where atmospheric water absorption severely affects the red wing. Because late-type stars are usually selected as interferometric calibrators to correct for these telluric features, stellar absorption lines around 2.3 µm complicate the correction despite the use of stellar templates. To address this, we observed B star calibrators for three nights (different stars in different nights; Table 1); and by using the flux spectra from only these nights we were able to accurately recover the Bry line profile. Figure 1 shows that the line profiles from the three nights with early-type calibrators are consistent with each other. We averaged these line profiles, weighted by their statistical uncertainties, to derive the final line profile. The averaged line profile is displayed in both Figs. 3 and 4 to guide the eye, as both the differential phase and differential visibility amplitude scale with the normalised line profile. We do not observe a distinct narrow emission line component because it is weak and its width (Perez et al. 1989) is comparable to our 0.002 µm channel resolution.

<footnote n="2">² We set the thresholds of the pipeline, snr-min-ft and vfactor-min-sc, to zero.</footnote>
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** uv coordinates of the binned data. Each coloured stripe spans the averaged uv coverage of a given baseline over the spectral range. From bin 1 to bin 5, the uv coverage of the 6 baselines rotates clockwise.</figcaption></figure>

The standard deviation of the three spectra around the line region is about 0.002 (normalised units), corresponding to ~5% (i.e. 0.002/0.04) of the line core region. This is consistent with the GRAVITY flux uncertainty we measured for 3C 273 (GC18). It is larger than the statistical uncertainties propagated from the uncertainty of each individual exposure spectrum, which is the result of stacking many exposures together in each night. The systematic uncertainty mainly arises from the variation of the sky absorption and the calibrator data. Therefore, we quadrature summed the systematic uncertainty (the rms) and the statistical uncertainty of the normalised spectral flux in each wavelength channel with a minimum systematic uncertainty of 0.002. The resulting total flux uncertainty is largely uniform across the whole line region. We note that our results are not significantly different if we adopt a slightly different threshold (e.g. 0.0015–0.0025).

## 2.3. The averaged differential phase

The phase of the pipeline reduced visibility shows significant (≥1°) instrumental features across the spectral band, which also varies between exposures. As discussed in Appendix A, we find that the instrumental features primarily consist of two components. A variable component is introduced by the dispersion of the air in the non-vacuum delay line of the VLTI, while a stable component is likely generated inside the cryostat of GRAVITY. Crucially, the variable component has a stable profile as a function of wavelength and only the total amplitude changes with each exposure. We remove these systematic features by fitting and subtracting a simple instrumental phase model – a stable component plus a fixed phase profile with a variable amplitude. Using calibrator star data, we demonstrate that the systematic uncertainty of our flattening method reaches ≤0.05° across most of the wavelength coverage (2.05–2.40 µm) of GRAVITY.

When we remove the residual instrumental features described above, we also mask out the wavelength range (2.25–2.35 µm) where the broad Bry line in IRAS 09149–6206 dominates. The mask does not affect the flattening (see Appendix A). Similarly, the uncertainty of the differential phase is estimated from the rms of the flattened phase at 2.05–2.40 µm, masking the wavelength range of the broad Bry.

In addition, we also fit a second-order polynomial function locally around the Brγ line, avoiding the line core (2.20–2.27 µm and 2.31–2.38 µm), to further flatten the phase around the line. This is necessary because IRAS 09149–6206 shows a slight positive phase gradient across the broad emission lines, including Brγ (2.2896 µm), Brδ (2.0551 µm), and Paα (1.9821 µm)<fnref n="3" />. Without this step, the pipeline would create a phase dip around the blue wing of the Brγ line because it derives the phase reference without considering the scientific signal. Using the calibrator data, we find that the additional flattening procedure does not increase the noise with respect to the systematic uncertainty (~0.05°; Appendix A). Finally, as shown in Fig. 2, we split the data collected across the three years into 5 angular uv bins in order to minimise smearing of the phase signal. We then averaged the differential phase in each spectral bin, weighted by the uncertainty, resulting in total integration times for each bin between 1.3 and 1.8 h. The full set of binned differential phase spectra can be seen in Appendix D while in Fig. 3 we show the averaged differential phase spectra for each baseline. The six baselines have different orientations and uv distances (Fig. 2). The longer baselines are more sensitive to the small scale signal. Therefore, the signals measured by different baselines behave differently.

## 2.4. The differential visibility amplitude

To produce the differential visibility amplitude, we fit and divide out a second-order polynomial from the pipeline reduced SC visibility amplitude while masking the Brγ line region. The uncertainty is estimated from the rms of the resulting normalised continuum channels. We then average the differential amplitudes from each exposure together, weighted by their inverse variance, and the final uncertainty is calculated by propagating the individual exposure uncertainties. Instead of uv binning, we simply average all exposures for a single baseline together, since the signal in the differential visibility amplitude is much weaker than the differential phase. As shown in Fig. 4, the differential visibility amplitude of UT4–UT1 is well above 1 and follows the profile of the emission line. The small dip in the differential visibility amplitude signal around the central wavelength of the line may be due to the narrow Brγ component; but it does not affect our analysis because it is barely significant compared to the noise level.

## 2.5. Continuum visibility from the fringe tracker

The fringe tracker kept a record of fringe measurements at 300 Hz throughout the exposures. There are six spectral channels across K band. The visibility amplitude and closure phases are important to constrain the continuum emission. In Fig. 5, we plot the closure phase distribution for the four triangles. All distributions are consistent with an average closure phase between –1° and 0°.

# 3. Locating the broad line region

The differential visibility amplitude spectra provide direct evidence that the BLR has been unambiguously detected. The signal is especially clear in the UT4–UT1 baseline shown in Fig. 4,

<footnote n="3">³ Our spectrum only covers the red wing of the Paα line. The overlapping line profile and phase signal of Brδ and Paα are significantly affected by atmospheric carbon dioxide absorption at ~2 µm. Therefore, it is very difficult to incorporate the Brδ line in the analysis.</footnote>
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Averaged differential phase spectra (coloured) and the normalised Bry spectrum (grey histogram) for IRAS 09149–6206. Solid and dashed black curves show the differential phase spectra of the best-fit BLR models (see Sect. 4). For the purpose of presentation, we average all of the data in each baseline. BLR modelling, however, was done on the five bins based on uv coverage as shown in Fig. 2. The same Bry spectrum is displayed in all of the panels.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Averaged differential visibility amplitude spectra (coloured) and the normalised Bry spectrum (grey histogram) for IRAS 09149–6206. There is a clear positive differential visibility amplitude signal in the UT4–UT1 spectrum, which follows well the line profile. Solid and dashed black curves show the differential visibility amplitude of the best-fit BLR models (see Sect. 4.4). Several channels with large scatter at >2.3 µm are masked for clarity. The same Bry spectrum is displayed in all of the panels.</figcaption></figure>

where the differential visibility amplitude significantly increases across the BLR line profile. The channels dominated by broad Bry emission display higher visibility amplitude, hence smaller size, than those dominated by the continuum (see Sect. 4.4). Consistent with GC20a, this indicates that the BLR is more compact than the near-infrared continuum, which traces the hot dust distribution around the AGN. IRAS 09149–6206 also shows a strong differential phase signal primarily in the UT4–UT2 and UT4–UT1 baselines. Remarkably, as is apparent in Fig. 3, the signal is also entirely positive and peaks near the line centre, which is different from the “S-shape” seen in 3C273 (GC18) that crosses zero at the line centre.

A differential phase following the shape of the line profile is expected for a constant phase difference between the hot dust continuum and the Bry emission. Both sources are marginally resolved, which is strongly supported by the ~0° closure phase in

Fig. 5. This means that to first order a phase difference measures an offset in photocentre position. The phase signal caused by this offset, we hereafter refer to as the “continuum phase”. By construction, the differential phase data are referenced to the photocentre position of the hot dust continuum. Hence, we have the differential phase, $\Delta\phi_{\lambda} = -2\pi[f_{\lambda}/(1+f_{\lambda})]u \cdot x_{\text{BLR},\lambda}$, where $f_{\lambda}$ is the line flux at wavelength $\lambda$ relative to a continuum level of unity, $u$ is the uv coordinate of the baseline, and $x_{\text{BLR},\lambda}$ is the coordinate of the photocentre w.r.t. the centroid of the continuum (see Appendix B for details). Fitting the photocentre coordinate of each channel to the differential phase data of 30 baselines (6 baselines × 5 angular bins), we can reconstruct the photocentres of the IRAS 09149–6206 Bry line emission.

The result in Fig. 6a shows a systematic offset of the BLR photocentres from the origin by ~120 µas to the east. Moreover, there is clear evidence for a velocity gradient that is
<figure><img src="image_6.png" /><figcaption>Fig. 5. Kernel density estimates for the closure phase distribution from the FT data of IRAS 09149–6206. The closure phase distributions for the four triangles are consistent with 0°.</figcaption></figure>

nearly perpendicular to the offset: The blueshifted channels lie predominantly to the south of the red-shifted channels. While the separation between individual channels is only moderate given the uncertainties, the general gradient from north to south appears robust. Following GC18, we estimate the significance of the offset between the blue and red channels with an F-test, comparing the null hypothesis that the phase signal is produced by a single position of the unresolved BLR with the hypothesis that the blue and red channels are at two distinct positions. Using the same channels as shown in Fig. 6a (6 blue channels <2.2896 µm and 5 red channels >2.2896 µm), we find the red-blue photocentre offset is at >8σ significance. If we use only alternate channels, in order to avoid the impact of possible data correlation, the same method still yields a significance >5σ<fnref n="4" />. According to this simple model, the average photocentre displacement on sky is 33 ± 8 µas, which can be considered as the lower limit of the separation of the photocentres.

This red-blue photocentre displacement, while model-independent, is only a lower limit to the true physical BLR size. To measure the physical size and constrain the BLR kinematics, a model is required. We adopt a flexible model, with a full differential phase (see also Appendix B),

$$ \Delta\phi_{\lambda} = [f_{\lambda}/(1 + f_{\lambda})](\phi_{\text{BLR},\lambda} - 2\pi \mathbf{u} \cdot \mathbf{x}_{\text{o}}), \quad (1) $$

where $\mathbf{x}_{\text{o}}$ is the coordinate of the origin of the BLR with respect to the centroid of the continuum. This velocity-independent photocentre offset could for example result from asymmetric structure in the continuum, or a physical offset between the BLR and hot dust continuum. The kinematic model described in the next section provides the velocity-dependent phase $\phi_{\text{BLR},\lambda}$ of the BLR itself.

## 4. Rotation versus outflow in the broad line region

The GRAVITY measurements of the line and phase profile for 3C 273 were modelled with a simple model consisting of a symmetric distribution of clouds in circular orbits, which fitted those

data very well GC20a. In that particular case, the model is supported by the symmetric profile of the broad Paα emission line, as well as the fact that the orientation of the radio jet is almost perfectly perpendicular to the gradient of the photocentres. With more limited knowledge available for IRAS 09149–6206, it is not clear a priori whether the velocity gradient we observe reflects rotational or radial motion of the BLR. As such, we adopt a flexible BLR model with a generalised prescription of the BLR dynamics (Pancoast et al. 2014b). The specific code we use here was developed by Stock (2018) and already adopted in the analysis of 3C 273 (GC18).

Following a general description of the model and its various parameters in Sect. 4.1, we compare the fits for two different implementations, allowing in both cases for an offset between the continuum and the BLR as discussed above. For the first case, in Sect. 4.2, we reduce the model to circular orbits as was applied in the case of 3C 273. For the second case, in Sect. 4.3, we apply the full model, allowing for radial motions. Section 4.4 then checks the consistency of both models with the observed differential visibility amplitudes. Finally, a comparison of the fits in Sect. 4.5 shows that the goodness of fit does not indicate a preference based on the data alone, and it is instead the astrophysical implication that leads to a preference of one model.

### 4.1. The generalised BLR model

The generalised BLR model was developed by Pancoast et al. (2014b, P14), with the original purpose to model the spectra and light curves from RM campaigns. The P14 model describes the BLR as a large number of non-interacting clouds and includes a large number of parameters, summarised in Table 2, that define the position and motion of each of those clouds. Below, we briefly describe how these affect the geometry and dynamics of the model.

The first set of parameters defines the locations of the clouds. Their distances from the black hole are given as

$$ r = R_S + F R_{\text{BLR}} + g (1 - F) \beta^2 R_{\text{BLR}}, \quad (2) $$

where $R_S = 2G M_{BH}/c^2$ is the Schwarzschild radius, $R_{BLR}$ is the mean radius, $F = R_{\min}/R_{BLR}$ is the fractional inner radius, $\beta$ is the shape parameter, and $g = p(x|1/\beta^2, 1)$ is drawn randomly from a Gamma distribution

$$ p(x|a, b) = \frac{x^{a-1} e^{-x/b}}{\Gamma(a) b^a}, \quad (3) $$

where $\Gamma(a)$ is the gamma function. Using a shape parameter in this way provides enough flexibility to reproduce several qualitatively different radial distributions, namely a Gaussian ($0 < \beta < 1$), exponential ($\beta = 1$), or heavy-tailed ($1 < \beta < 2$) profile. The angular distribution of the clouds is given by

$$ \theta = \arccos(\cos \theta_o + (1 - \cos \theta_o) \times U^\gamma), \quad (4) $$

where $\theta_o$ is the angular thickness of the distribution (defined as the angle between the mid-plane and the upper edge of the distribution) and $U$ is a random number drawn uniformly between 0 and 1. Setting $\gamma > 1$ concentrates more clouds closer to the maximum angular height $\theta_o$. The structure is viewed at an inclination angle $i$ (where $i = 0^\circ$ is defined to be face-on) and rotated in the plane of the sky so that the line of nodes is at position angle PA (measured east of north). A weight is assigned to each cloud to represent the relative strength of its emission, and is defined as

$$ w = 0.5 + \kappa \cos \phi, \quad (5) $$

<footnote n="4">⁴ The result remains robust independent of whether one includes the bluest channel that is furthest offset from the photocentres of the other channels.</footnote>
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Best-fit centroids to the differential phases of (a) observed data of IRAS 09149–6206, (b) best-fit Keplerian model described in Sect. 4.2, and (c) best-fit outflow model described in Sect. 4.3. The colour code represents the wavelength of the channels around the centre of the Bry line. The coloured ellipses around each centroid in panel a represents the 68% (1σ) credible intervals of the uncertainty. The red plus sign at the origin represents the photocentre of the continuum. The black arrows in panels b and c indicate the origin of the BLR according to the inferred offset of the models.</figcaption></figure>

<table><caption>**Table 2.** Parameters of the BLR model with definitions, priors, and units where appropriate (all angles are in radians).</caption><thead><tr><th></th><th></th><th></th></tr></thead><tbody><tr><td>R<sub>BLR</sub></td><td>Mean radius of the BLR</td><td>LogUniform(10<sup>-4</sup>, 10 pc)</td></tr><tr><td>F</td><td>Minimum radius of the BLR in units of R<sub>BLR</sub></td><td>Uniform(0, 1)</td></tr><tr><td>β</td><td>Unit standard deviation of BLR radial profile</td><td>Uniform(0, 2)</td></tr><tr><td>θ<sub>o</sub></td><td>Angular thickness measured from the mid-plane</td><td>Uniform(0, π/2)</td></tr><tr><td>i</td><td>Inclination angle</td><td>Uniform(cos i(0, π/3))</td></tr><tr><td>PA</td><td>Position angle of the line of nodes on sky (east of north)</td><td>Uniform(0, 2π)</td></tr><tr><td>κ</td><td>Anisotropy of the cloud emission</td><td>Uniform(-0.5, 0.5)</td></tr><tr><td>γ</td><td>Clustering of the clouds at the edge of the disk</td><td>Uniform(1, 5)</td></tr><tr><td>ξ</td><td>Mid-plane transparency</td><td>Uniform(0, 1)</td></tr><tr><td>M<sub>BH</sub></td><td>Black hole mass</td><td>LogUniform(10<sup>5</sup>, 10<sup>10</sup> M<sub>☉</sub>)</td></tr><tr><td>f<sub>ellip</sub></td><td>Fraction of clouds in bound elliptical orbits</td><td>Uniform(0, 1)</td></tr><tr><td>f<sub>flow</sub></td><td>Flag for specifying inflowing or outflowing orbits</td><td>Uniform(0, 1)</td></tr><tr><td>θ<sub>e</sub></td><td>Angular location for radial orbit distribution</td><td>Uniform(0, π/2)</td></tr><tr><td>σ<sub>ρ,circ</sub></td><td>Radial standard deviation for circular orbit distribution</td><td>LogUniform(0.001, 0.1)</td></tr><tr><td>σ<sub>Θ,circ</sub></td><td>Angular standard deviation for circular orbit distribution</td><td>LogUniform(0.001, 1)</td></tr><tr><td>σ<sub>ρ,radial</sub></td><td>Radial standard deviation for radial orbit distribution</td><td>LogUniform(0.001, 0.1)</td></tr><tr><td>σ<sub>Θ,radial</sub></td><td>Angular standard deviation for radial orbit distribution</td><td>LogUniform(0.001, 1)</td></tr><tr><td>σ<sub>turb</sub></td><td>Normalised standard deviation of turbulent velocities</td><td>LogUniform(0.001, 0.1)</td></tr><tr><td>λ<sub>emit</sub></td><td>Central wavelength of the emission line</td><td>Norm(2.2896, 0.002 μm)</td></tr><tr><td>f<sub>peak</sub></td><td>Peak flux of the normalised line profile</td><td>Uniform(0.05, 0.065)</td></tr><tr><td>(x<sub>o</sub>, y<sub>o</sub>)</td><td>Offset of the origin of the BLR</td><td>Uniform(-1, 1 mas)</td></tr></tbody></table>

**Notes.** The priors for most parameters are specified here in one of two ways. Uniform(min,max) denotes uniform sampling over the range specified. LogUniform(min,max) indicates that the logarithm of the parameter is sampled uniformly over the logarithm of the range. The prior for the inclination angle (i) is set between 0 and π/2 so that cos i is uniformly sampled between 0 and 1. The prior for the central wavelength of the emission line (λ<sub>emit</sub>) follows a Gaussian distribution centred at 2.2896 µm and with a standard deviation 0.002 µm.

where κ is a parameter in the range (-0.5, 0.5) reflecting any anisotropy of the emission, and φ is the angle between the line of sight from the cloud to the observer and to the central ionising source. The mid-plane transparency is modelled using the parameter ξ which controls the fraction of clouds located behind the equatorial plane. If ξ = 1 then the clouds are evenly distributed on either side of the equatorial plane, while ξ = 0 means that all the clouds are in front of it.

The remainder of the parameters define the kinematics of the clouds, under the assumption that their motions are governed entirely by the gravitational potential of the black hole. A fraction, f<sub>ellip</sub>, of the clouds are put on bound elliptical orbits. The rest are placed on much more elongated orbits, which are dominated by radial motion. A single parameter, f<sub>flow</sub>, is used

as a binary switch to control whether the radial motion is inflow (0 < f<sub>flow</sub> < 0.5) or outflow (0.5 < f<sub>flow</sub> < 1).

For the bound elliptical orbits, radial and tangential velocities, v<sub>r</sub> and v<sub>φ</sub>, are drawn randomly from a distribution centred on the point {0, v<sub>circ</sub>} in the v<sub>r</sub>-v<sub>φ</sub> plane (see Fig. 2 in Pancoast et al. 2014b), where v<sub>circ</sub> = √G M<sub>BH</sub>/r is the circular velocity. The distribution itself follows an ellipse in the v<sub>r</sub>-v<sub>φ</sub> plane that is defined as

$$ \frac{v_r^2}{2v_{\text{circ}}^2} + \frac{v_{\phi}^2}{v_{\text{circ}}^2} = 1 \qquad (6) $$

and is defined to be a Gaussian with standard deviation σ<sub>Θ,circ</sub> along the ellipse and σ<sub>ρ,circ</sub> perpendicular to it.
Because the ellipse naturally connects points around {0, $v_{\text{circ}}$} corresponding to circular orbits, with points around {$\pm\sqrt{2}v_{\text{circ}}$, 0} corresponding to highly elongated orbits, the orbits dominated by radial motion can be defined in a similar way. In this case, there is an additional parameter $\theta_e = \arctan(|v_{\phi}/v_i|)$ that defines the angular location around the ellipse where the distribution is centred. If $\theta_e$ is close to $\pi/2$, the orbit distribution is centred around {0, $v_{\text{circ}}$} exactly as for the bound elliptical orbits and there is very little inflow or outflow. As $\theta_e$ approaches 0, the centre of the distribution shifts to {$\pm\sqrt{2}v_{\text{circ}}$, 0} where orbits are dominated by radial motion at the escape velocity. The distribution is defined around the point on the ellipse corresponding to $\theta_e$ as a Gaussian with standard deviation $\sigma_{\theta_e, \text{radial}}$ along the ellipse and $\sigma_{\rho, \text{radial}}$ perpendicular to it. The units of $\sigma_{\rho, \text{circ}}$ and $\sigma_{\rho, \text{radial}}$ are given in terms the circular velocity of the clouds, while $\sigma_{\theta_e, \text{circ}}$ and $\sigma_{\theta_e, \text{radial}}$ are given as a fraction of $\pi$. The last velocity component is a random velocity $v_{\text{turb}}$ describing the macroturbulence. This is randomly sampled from a Gaussian distribution with standard deviation $\sigma_{\text{turb}}v_{\text{circ}}$, and added to the line-of-sight velocity of each cloud.

The full relativistic Doppler effect and gravitational redshift are also taken into account in generating the spectrum and differential phases. For each cloud, the intrinsic line width is assumed to be negligible. The wavelength of the line is shifted from $\lambda_{\text{emit}}$ to $\lambda_{\text{obs}}$ (both in the observed frame) by

$$ \lambda_{\text{obs}} = \lambda_{\text{emit}} \frac{1 - \frac{v_{\text{los}}}{c}}{\sqrt{1 - \frac{v^2}{c^2}}} \frac{1}{\sqrt{1 - \frac{R_i}{r}}}, \quad (7) $$

where $v_{\text{los}}$ is the total line-of-sight velocity and $v$ is the total velocity. Finally, we bin the clouds in the spectral channels according to their $\lambda_{\text{obs}}$ and sum their weights to derive the normalised line profile. The profile is then scaled by $f_{\text{peak}}$, so that it has the same normalisation as the continuum. The projected coordinates perpendicular to the line of sight are averaged in each bin according to the cloud weights to derive the photocentre of each spectral channel. The differential phase is then calculated with Eq. (1) with,

$$ \phi_{\text{BLR},i} = -2\pi \frac{f_i}{1 + f_i} \mathbf{u} \cdot \left( \frac{\sum_i w_i \mathbf{x}_i}{\sum_i w_i} \right), \quad (8) $$

where $w_i$ and $\mathbf{x}_i$ are the weight and coordinate of the $i$th cloud of the BLR with $\lambda_{\text{obs}}$ within the wavelength channel $\lambda$.

Our prior assumptions on all these parameters, as given in Table 2, generally follow the choices of Pancoast et al. (2014b). The main exception is the inclination angle, which we require to be below 50°, because the BLR of a type 1 AGN is expected to be relatively face-on<fnref n="5" /> For the prior on $\lambda_{\text{emit}}$, we adopt a Gaussian distribution centred at 2.2896 µm based on the redshift measured by Perez et al. (1989). The standard deviation of the Gaussian distribution is 0.002 µm, which is the width of the spectral channel and the equivalent $\sigma$ of the line spread function for the GRAVITY science spectrograph. Three additional parameters that we include are the peak flux $f_{\text{peak}}$ of the emission line

normalised to the continuum, and the offset $\{x_o, y_o\}$ for the origin of the BLR. The prior range adopted for $f_{\text{peak}}$ is 0.05–0.065. The offset of the BLR, $\{x_o, y_o\}$, is allowed to vary by up to –1 to 1 mas in both directions.

## 4.2. Model with circular Keplerian rotation

The P14 model described above can easily be reduced to circular Keplerian rotation. Setting $\kappa = 0$ and $\xi = 1$ ensures that the clouds are equally weighted (i.e. emitting isotropically) and are distributed uniformly above and below the equatorial plane. Instead of using Eq. (4) to determine the initial angular distribution of the clouds, the circular Keplerian model distributes them uniformly between 0 and $\theta_o$ (GC18). To produce bound circular orbits, we set $f_{\text{ellip}} = 1$, $\sigma_{\rho, \text{circ}} = 0$, and $\sigma_{\theta_e, \text{circ}} = 0$. Finally, setting $\sigma_{\text{turb}} = 0$ ensures that there is no additional turbulence. We refer to this model as the “Keplerian model”, hereafter, for simplicity.

The flux and differential phase spectra of IRAS 09149–6206 are fit reasonably well by the Keplerian model, as shown in Fig. 3 for the averaged bins (see also Fig. D.1 for the individual bins). The most prominent phase signals in the UT4–UT2, UT4–UT1, and UT3–UT1 baselines are primarily due to the continuum phase produced by the offset between the BLR and the centre of the continuum emission. For this model, the BLR is ~120 µas east of the continuum centre. The cloud distribution of the model is shown in Fig. 7a, and the corresponding photocentres in Fig. 6b are qualitatively consistent with those reconstructed from the data. The best fitting parameters summarised in Table 3 indicate that a rather face-on disk with $i \approx 21^\circ$ is favoured by the data, which is consistent with the inclination found from an upcoming dust reverberation study (Hönig et al., in prep.). A low inclination is consistent with the Seyfert 1 classification and, as one would expect, low inclinations are generally inferred when fitting RM data of other objects. The disk is also very thick, with $\theta_o \approx 71^\circ$. Grier et al. (2017a) have suggested that, when fitting RM data, the thickness $\theta_o$ is always similar to the inclination angle $i$, due to degeneracy between the two quantities. Interferometry breaks this degeneracy. Our derived values for these parameters are significantly different, and the posterior distributions in Fig. D.2 show no particular coupling. The mean radius of $R_{\text{BLR}} = 65$ µas corresponds to 0.075 pc. The posterior distributions in Fig. D.2 indicate that there is some degeneracy between BLR radius, black hole mass, and the inclination angle, which is consistent with previous studies (Rakshit et al. 2015; GC18). Although the inferred line centre $\lambda_{\text{emit}} = 2.2892$ µm corresponds to a velocity offset of $\Delta v_{\text{BLR}} = -56$ km s$^{-1}$ with respect to the systemic velocity from the [O III] rotation curve (Perez et al. 1989), the uncertainty of $\Delta v_{\text{BLR}}$ is large enough that the modelled line centre is fully consistent with it.

In order to display the phase signal specific to the BLR, we subtract the best-fit continuum phase from the data in the three longest baselines (UT4–UT2, UT4–UT1, and UT3–UT1) as shown in Fig. 7c, and then average them. The residual BLR signal, shown in Fig. 7b, exhibits the expected S-shape profile for a rotating structure. Based on the analysis in Appendix A, and taking into account that several baselines were combined, the uncertainty of this phase is expected to be below 0.03° per spectral channel. As such, even though the ~0.1° BLR signal is several times weaker than the continuum phase signal, it remains a significant detection. We note that the S-shape profile is due to specifically fitting with the Keplerian model. The significance of the BLR phase signal is constrained model-independently with the reconstructed photocentres (Sect. 3). Finally, we also need to

<footnote n="5">⁵ Allowing $i$ to vary over the full range of 0°–90° leads to multiple modes in the posterior distribution for the outflow model described in Sect. 4.3 and a preference for a highly inclined geometry inconsistent with the expectation for a Seyfert 1. We therefore apply a more restrictive prior to the inclination angle. Our results are unaffected if we increase the upper limit of $i$ to, for example, 60°. We choose $i < 50^\circ$ because the posterior probability density starts to rise towards higher values of $i$.</footnote>
<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Panel a: cloud distribution of the best-fit Keplerian model. Each circle represents one cloud, colour-coded by the line-of-sight velocity. The zero velocity is defined at the best-fit $\lambda_{\text{emit}}$. The red plus sign at the origin represents the photocentre of the continuum. The orientation of the model is consistent with the observed photocentre gradient. Panel b: observed averaged differential phase from UT4–UT2, UT4–UT1, and UT3–UT1 after removing the continuum phase signal (blue) compared to the averaged differential phase from the best-fit BLR model (red). These baselines were chosen since they contain the strongest S-shape signal. Above, the observed line profile (black) is compared with the model line profile (red). The residual of the spectrum subtracting the model line profile is displayed in the lower panel. The left y axis corresponds to the averaged differential phase, while the right y axis corresponds to the line profile. Panel c: differential phase data and the best-fit models (solid lines) of the three baselines that show the strongest signal of the BLR component (dashed lines). The phase in panel b is calculated by averaging the phases of these three baselines after subtracting the best-fit continuum phases (dotted lines).</figcaption></figure>

<table><caption>**Table 3.** Inferred maximum a posteriori value and central 95% credible interval for the modelling of the spectrum and differential phase of IRAS 09149–6206.</caption><thead><tr><th>Parameters</th><th>Circular Keplerian model</th><th>Outflow model</th></tr></thead><tbody><tr><td>R<sub>BLR</sub> (µas)</td><td>65<sup>+30</sup><sub>-39</sub></td><td>50<sup>+38</sup><sub>-11</sub></td></tr><tr><td>R<sub>min</sub> (µas)</td><td>8.3<sup>+13.1</sup><sub>-7.4</sub></td><td>6.8<sup>+7.5</sup><sub>-5.0</sub></td></tr><tr><td>β</td><td>1.17<sup>+0.30</sup><sub>-0.33</sub></td><td>1.27<sup>+0.17</sup><sub>-0.32</sub></td></tr><tr><td>θ<sub>o</sub> (°)</td><td>71<sup>+16</sup><sub>-27</sub></td><td>61<sup>+20</sup><sub>-20</sub></td></tr><tr><td>i (°)</td><td>21<sup>+20</sup><sub>-8</sub></td><td>35<sup>+13</sup><sub>-10</sub></td></tr><tr><td>PA (° E of N)</td><td>130<sup>+29</sup><sub>-34</sub></td><td>219<sup>+27</sup><sub>-37</sub></td></tr><tr><td>κ</td><td>...</td><td>-0.32<sup>+0.44</sup><sub>-0.17</sub></td></tr><tr><td>γ</td><td>...</td><td>1.27<sup>+2.40</sup><sub>-0.20</sub></td></tr><tr><td>ξ</td><td>...</td><td>0.05<sup>+0.30</sup><sub>-0.04</sub></td></tr><tr><td>Offset (µas)</td><td>(121.6<sup>+6.5</sup><sub>-9.7</sub>, 4.5<sup>+8.4</sup><sub>-8.6</sub>)</td><td>(109.2<sup>+11.5</sup><sub>-21.8</sub>, -13.9<sup>+12.3</sup><sub>-17.4</sub>)</td></tr><tr><td>log (M<sub>BH</sub>/M<sub>☉</sub>)</td><td>8.06<sup>+0.41</sup><sub>-0.57</sub></td><td>7.70<sup>+0.41</sup><sub>-0.18</sub></td></tr><tr><td>f<sub>ellip</sub></td><td>1</td><td>0.19<sup>+0.35</sup><sub>-0.17</sub></td></tr><tr><td>f<sub>flow</sub></td><td>...</td><td>0.71<sup>+0.26</sup><sub>-0.20</sub></td></tr><tr><td>θ<sub>e</sub></td><td>...</td><td>5.0<sup>+22.9</sup><sub>-4.3</sub></td></tr><tr><td>log σ<sub>turb</sub></td><td>0</td><td>-1.87<sup>+0.80</sup><sub>-1.04</sub></td></tr><tr><td>Δv<sub>BLR</sub> (km s<sup>-1</sup>)</td><td>-56<sup>+78</sup><sub>-67</sub></td><td>380<sup>+208</sup><sub>-356</sub></td></tr><tr><td>χ<sup>2</sup><sub>t</sub></td><td>1.39</td><td>1.38</td></tr><tr><td>ln K</td><td>0</td><td>7.1 ± 0.2</td></tr><tr><td>Δ AIC</td><td>0</td><td>-12.6</td></tr><tr><td>Δ BIC</td><td>0</td><td>44.0</td></tr></tbody></table>

**Notes.** $\Delta v_{BLR}$ is the difference between the velocity derived from the best-fit $\lambda_{\text{emit}}$ and the systemic velocity based on the [O III] line (Perez et al. 1989). $\chi_t^2$ is the reduced $\chi^2$ of the models with the best-fit parameters. The Bayes factor, AIC, and BIC are relative to the Keplerian model.

consider the fit to the spectral line profile. The Keplerian model is only able to generate a symmetric line profile, which also means that $\lambda_{\text{emit}}$ is close to the wavelength of the peak of the line profile. However, the observed line profile of IRAS 09149–6206 is slightly asymmetric. As such, although the line profile is reasonably well matched by the model, a difference between the model and data, especially in the blue wing (e.g. 2.24–2.29 µm), is apparent in Fig. 7b. This issue is addressed in the next Section.

### 4.3. Model including radial motion

To be able to fit the asymmetries in the emission line profile, we apply the full P14 model. The best fitting parameters of the P14 model (Table 3) include $f_{\text{ellip}} \approx 0.2$ indicating that the majority of the clouds are on orbits with a dominant radial component, $\theta_e \approx 5^\circ$ indicating that the orbits are sufficiently elongated that the radial motion is very close to the escape velocity, and $f_{\text{flow}} > 0.5$ indicating that this radial motion is outwards. Together, these indicate that, although it is not required a priori, the configuration of the model preferred by the data is very much dominated by outflow. As such, hereafter we call this the “outflow model”.

As before, the phase signals shown in Fig. 3 (see also Fig. D.3 for the individual bins) are dominated by the continuum phase. The BLR offset of ~110 µas, which can be seen in Fig. 8a, is statistically consistent with that of the Keplerian model. The modest difference is due to the different BLR phase signals (Fig. 8b). The orientation and gradient of the photocentres in Fig. 6c, are also consistent with the data. The outflow model indeed better fits the line profile compared to the Keplerian model. No systematic residual is seen in Fig. 8b, especially in the blue wing. Because the two models fit the differential phase data equally well, as shown in Fig. 3, the total goodness of fit for the two models are nearly indistinguishable (Sect. 4.5).

The mean radius of $R_{BLR} = 50$ µas is slightly smaller than, but statistically consistent with, that of the Keplerian model. The cloud radial distribution given by $\beta = 1.27$ prefers to be exponential or heavy-tailed, with a small inner radius of $R_{min} = 6.8$ µas. The model has $\text{PA} \approx 219^\circ$, which is 90° different from Keplerian model. The reason is simply that the BLR kinematics, and hence the orientation of the velocity gradient, are now dominated by radial motion rather than rotation. This model also prefers anisotropic emission from the clouds, with $\kappa = -0.32$ indicating that line emission is stronger from the inner illuminated side of the clouds and hence the far side of the distribution, similar to many of the results inferred from RM data. The parameter $\sigma_{\text{turb}} = 0.013$ implies that additional macroscopic turbulence is not significant. As for the Keplerian model, $\theta_0 \sim 61^\circ$ indicates that the cloud trajectories are distributed over a wide range of angles from the mid-plane. We
<figure><img src="image_9.png" /><figcaption>Fig. 8. Panel a: cloud distribution of the best-fit outflow model. The symbols and lines are the same as for Fig. 7, except that in panel a the sizes of the circles scale with the weight of the cloud.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 9. Edge-on (panel a) and face-on (panel b) views of the best-fit outflow model, except that the PA is adjusted to 180° and the BLR is moved back to the origin of the plot for clarity. In both panels, the colour coding represents the line-of-sight velocity for each cloud. Following Pancoast et al. (2014a,b), Δx is along the line of sight, Δy is the direction of the right ascension, and Δz is the direction of declination.</figcaption></figure>

also note, similarly to the Keplerian model, that the thickness and inclination have somewhat different values, and there is no evidence for degeneracy between them in the posterior distributions shown in Fig. 4. Indeed, except for the greater thickness of the BLR in IRAS 09149–6206, the configuration is rather similar to that inferred from RM data for Arp 151 (Pancoast et al. 2014a), Zw 229–015 (Williams et al. 2018), or Mrk 142 (Li et al. 2018). Finally, the best fit line centre is $\lambda_{emit} = 2.2923 \mu m$ which corresponds to an offset $\Delta v_{BLR} = 380 \text{ km s}^{-1}$ from the systemic velocity. A discussion of whether this is physically plausible is deferred to Sect. 4.5.

One of the most important aspects of this model is that the differential phase of the BLR component is very different to the S-shape seen in the Keplerian model. It is clear from Fig. 8b that the continuum subtracted phase data fitted by the BLR model is dominated by a positive signal on the red-shifted side of the line profile. Figure 8c illustrates the decomposition of the BLR phase and continuum phase components. The asymmetric BLR phase signal is produced by two main effects: (i) The BLR kinematics are dominated by outflow (as discussed above), and (ii) $\xi \approx 0$ means that the mid-plane is opaque. The impact of this second effect is discussed below in the context of the distribution and motions of the clouds.

The edge-on and face-on views of the cloud distribution presented in Fig. 9 can shed more light on the role of the mid-plane obscuration in generating the positive phase signal. The edge-on view clearly shows that the cloud distribution extends far above the mid-plane because $\theta_0 \approx 60^\circ$; while $\xi \approx 0$ means that the mid-plane obscuration is so strong that there are few clouds below it. As discussed above, our model is significantly dominated by

outflowing clouds. The blue-shifted clouds are on the near side towards the observer, while the red-shifted clouds are on the far side. In addition, the anisotropy parameter $\kappa \approx -0.3$ means that the weighting applied to clouds on the far side is much larger than for the near side. This effect is indicated in the figure by the size of the circles representing the clouds: On the side nearer the observer, the circles are much smaller than those on the far side. The mid-plane obscuration and inclination angle together mean that, as is apparent in Fig. 8a, the blue-shifted clouds are distributed fairly symmetrically around the centre of the BLR. This means that their photocentre is close to the BLR centre and the corresponding differential phases are close to $0^\circ$. In contrast, the red-shifted clouds are primarily located to the north-east of the BLR centre, so the corresponding red-shifted channels show significant differential phase signal. As a result, as shown in Fig. 6, the origin of the BLR coincides with the blue-shifted channels instead of with the channel associated with the line peak.

## 4.4. Model prediction of the differential visibility amplitude

The measured differential visibility amplitude is useful to provide an independent check of the BLR model fits. The differential visibility amplitude ($V_{diff}$) is derived by normalising the total visibility amplitude with the visibility amplitude of the continuum $V_c$. In each spectral channel $\lambda$, this is

$$V_{\text{diff}}(\lambda) = \frac{1 + f_{\lambda} V_{\text{BLR}}(\lambda)/V_c(\lambda)}{1 + f_{\lambda}}, \qquad (9)$$

where $V_{BLR}$ is the visibility amplitude of the line emission of the BLR. At the wavelength of the line emission, one will find
$V_{\text{diff}} > 1$ if $V_{\text{BLR}}/V_c > 1$, this is if the BLR emission is more compact than the continuum emission. This can be seen in Fig. 4, in particular for the UT4–UT1 baseline. Similar results have been reported for 3C 273 (GC18) and PDS 456 (GC20a).

The visibility amplitude of the continuum emission of IRAS 09149–6206 has already been studied in GC20a. Using the visibility amplitude from the fringe tracker channel and the differential visibility from the science channel, we derived consistent full width at half maximum (FWHM) sizes $0.54 \pm 0.05$ mas and $0.64 \pm 0.06$ mas, respectively, for a circular Gaussian profile. This indicates that the continuum emission is only marginally resolved. For consistency, we adopt a Gaussian profile with FWHM = 0.6 mas to calculate

$$V_c(\lambda) = \exp\left(-\frac{(\pi FWHM)^2(u^2 + v^2)}{4 \ln 2}\right). \qquad (10)$$

Following Waisberg et al. (2017), we calculate $V_{\text{BLR}}$ from the Keplerian and outflow models from the second moment of the source emission,

$$V_{\text{BLR}}(\lambda) \approx 1 - \frac{2\pi^2}{\mu_{00,\lambda}} (u^2 \tilde{\mu}_{20,\lambda} + v^2 \tilde{\mu}_{02,\lambda} + 2uv \tilde{\mu}_{11,\lambda}), \qquad (11)$$

where $\mu_{00,\lambda} = \sum_{\lambda} w_i$ is the total intensity of the BLR line emission (zero-order moment), and summing up the weight ($w_i$ for the ith cloud) of all of the clouds that belong to each spectral channel $\lambda$. In addition $\tilde{\mu}_{pq,\lambda} = \sum_{\lambda} w_i(l_i - l_{c,\lambda})^p(m_i - m_{c,\lambda})^q$ is the relative moment of the line emission with respect to the photo-centre of the spectral channel

$$(l_{c,\lambda}, m_{c,\lambda}) = \left( \frac{\sum_{\lambda} w_i l_i}{\mu_{00}}, \frac{\sum_{\lambda} w_i m_i}{\mu_{00}} \right). \qquad (12)$$

The marginally resolved approximation here is the same one used in calculating the BLR differential phase. It is valid when the BLR is not more than marginally resolved ($2\pi u \cdot x \ll 1$ or $R_{\text{BLR}} \ll 1$ mas). We calculated the differential visibility amplitude for both models, and plot them over the data in Fig. 4 taking into account the instrumental line spread function. The predicted differential visibility amplitudes for the two models are almost indistinguishable because the derived BLR size $R_{\text{BLR}}$ is almost the same at ~60 µas in both cases. The clear signal in the UT4–UT1 baseline is consistent with the compact BLR size inferred from the differential phase data. We note that although this is sensitive to the adopted size of the continuum emission, it produces qualitatively similar results for continuum FWHM in the range 0.54–0.64 mas. Additional details are discussed in Appendix E.

### 4.5. Comparing the models

In the preceding sections, we have discussed the Keplerian and outflow models individually from a phenomenological perspective. The offset and the orientation of the velocity gradient of the models are qualitatively consistent with those derived directly from the observed data of IRAS 09149–6206. We have also shown that the size of the BLR in both cases is consistent with the differential visibility amplitudes, which provide a direct comparison to the measured size of the continuum. Here we try to compare them, both in terms of a statistical perspective and with reference to the literature.

We calculate the reduced $\chi^2(\chi^2_r)$, Bayes factor (K), Akaike information criterion (AIC) and Bayesian information criterion (BIC), in order to compare the goodness of fits for the two models. These quantities are described in more detail in Appendix C,

and their values are given in Table 3 (in terms of the difference of the outflow model with respective to the Keplerian model) together with the relevant set of best fitting parameters. The $\chi^2_r$ is almost exactly the same for both models. In contrast, the other quantities do indicate a preference. A high Bayes factor and negative $\Delta$AIC favours the outflow model. However, a negative $\Delta$BIC favours the Keplerian model. This difference reflects the different approaches that the criteria adopt, in penalising the free parameters versus the prior information. Using mock data generated from the inferred parameters with the Keplerian and outflow models (see Appendix E), we confirm that the Bayes factor seems to provide a more reliable model selection for this case.

Unfortunately, the literature provides no useful information about the spatially resolved radio emission nor the gas kinematics that might be able to shed light on the interpretation of the BLR kinematics. In addition, the early interferometric measurements (Kishimoto et al. 2011; Burtscher et al. 2013; López-Gonzaga et al. 2016) barely resolve the mid-infrared continuum emission.

Nevertheless, we can compare the models from the point of view of the physical interpretation. The key parameter that generates the asymmetric line profile for the outflow model is the mid-plane transparency, with $\xi \sim 0$ indicating that it is opaque. In contrast, the mid-plane of the Keplerian model is, by our definition, fully transparent. It is not immediately clear what might physically cause an opaque mid-plane. Absorption by gas is not possible because any given cloud could only absorb a narrow range of velocities; and so it would have to be due to extinction by dust. But to obscure Brγ emission requires an optical depth $\tau \sim 1$ at 2.2 µm. For a standard gas-to-dust ratio, this would require a column of $N_H > 10^{22} \text{cm}^{-2}$ of dusty gas that is nominally located inside the dust sublimation radius. This might be possible for the BLR concept proposed by Baskin & Laor (2018) in which, because the emission from the accretion disk is anisotropic, large ($\ge 0.3$ µm) graphite grains can survive close to the disk plane at radial scales associated with the BLR. Indeed, the model requires this, since it purports that the BLR may be a failed dusty wind – failed because while dust opacity allows a wind to be launched, the dust sublimates once clouds move upwards. In the context here, the key point is that it does imply a possible physical source for the mid-plane opacity without affecting the optical/UV emitting inner accretion disk. However, this model would also imply that much of the hot dust should exist on the same spatial scales as the BLR, while our interferometric data (especially Fig. 4), indicate that the hot dust is much more extended. As such, we would argue that it is difficult for this model to explain such a high mid-plane opacity in a way that is consistent with the data.

An alternative possibility arises because, for a BLR dominated by clouds on radial outflowing trajectories, it is not entirely clear whether the parameters, such as the inclination angle and the disk thickness, should still retain their original physical interpretation. The best-fit outflow model tends rather to mimic a polar outflow, although intrinsically limited by its disk-like construct. However, a truly polar outflow would be inconsistent with the broad concept described by Elvis (2000) in which, because the outflow originates in a disk, it has a major rotational component (as observed in PG 1700+518 by Young et al. 2007) and is directed at an angle significantly offset from the polar direction. Consistency with physically motivated rotating wind models (Everett 2005; Keating et al. 2012; Mangham et al. 2017) is difficult to achieve for the outflow model here, which is strongly polar. On the other hand, the Keplerian model is plausibly
consistent with a disk wind because of the wide angle $\theta_o$ above and below the mid-plane over which its cloud trajectories are distributed.

Another impact of the opacity of the mid-plane in the outflow model is to reduce the number of red-shifted clouds, which means that the resulting line emission is dominated by the blue-shifted clouds. The model, therefore, prefers a large $\lambda_{emit} \approx 2.2923 \mu m$. This moves the line profile from the model back so that it matches the wavelength of the observed line profile. The best fitting $\lambda_{emit}$ corresponds to a ~380 km s⁻¹ shift with respect to the systemic velocity, which is slightly below the 2-$\sigma$ lower boundary of the probability distribution (see Table 3 and Fig. D.4). This deviation is moderately significant given our spectral resolution. The main issue is that the systemic velocity has been derived from the symmetry of the rotation curve of the [O III] line that includes significant regions at radii $\ge$3 kpc, where its shape is flat (Perez et al. 1989) and is a reliable method. The outflow model therefore requires that the black hole is offset by ~380 km s⁻¹ from the expected velocity. While black hole recoil in merging systems has been discussed extensively in the literature and a significant minority are expected to have recoil speeds exceeding 500 km s⁻¹ (Schnittman & Buonanno 2007), there have been few convincing cases for such candidates (Komossa et al. 2008; Eracleous et al. 2012; Komossa 2012). Moreover, this particular case would require a remarkable coincidence that the black hole recoil velocity exactly matches the outflow velocity of the BLR clouds. This is another argument against the outflow model.

We can try to avert this problem by fixing $\lambda_{emit} = 2.2896 \mu m$ so it exactly matches the systemic velocity. Doing this results in a similar configuration to that in Fig. 9, except that the weighting of the blue clouds in the near side is pushed to its lower limit at $\kappa \approx -0.5$. And although the profile of the line wing is still matched reasonably well, the line core can no longer be fit as well as in Fig. 8b. Thus, although the flexibility of the P14 model means it is still able to fit the interferometric data well, the original rationale for trying the outflow model, with its increased number of parameters, is lost. A similar situation occurred when we fixed $\xi = 1$ so that we avoid the mid-plane transparency problem. The line profile then becomes completely symmetric just as in the Keplerian model and again the rationale for using the outflow model is lost.

Taking all these arguments together, we strongly favour the results of the Keplerian model as the most likely for IRAS 09149–6206, and caution against over-interpreting the best-fit parameters of the outflow model. Nevertheless, we emphasise that the main results – the photocentre offset and gradient, the BLR radius, and the black hole mass – inferred from the Keplerian model and the outflow model are statistically consistent. Despite these similarities in the key parameters, the two models are expected to show different characteristic features in velocity-resolved RM data (Peterson 2014), and it would certainly be interesting to compare our results with dynamical modelling of high quality RM data for IRAS 09149–6206.

**5. Black hole mass and the radius-luminosity relation**

One of the parameters in the BLR model is the black hole mass. For the Keplerian model, Table 3 shows that the best fitting value is $1.1 \times 10^8 M_\odot$. The posterior probability distributions in Fig. D.2 show that it is correlated with the mean radius $R_{BLR}$ and the inclination angle $i$, and this is likely what drives the large formal uncertainty of +0.3/-0.4 dex. Interestingly, Fig. D.4 shows

<figure><img src="image_11.png" /></figure>

**Fig. 10.** Comparison of the BLR radii for IRAS 09149–6206 and 3C 273 measured with GRAVITY to those based on the RM time lag. The black circles are the RM measurements collected from the literature by Du & Wang (2019), which mainly include data from Bentz et al. (2013) and the SEAMBH campaign (see text for details). The blue squares are based on Hβ time lags from the SDSS-RM project Grier et al. (2017b). The dashed line is the best-fit radius–luminosity relation from Bentz et al. (2013).

that the correlation is weaker for the outflow model, and that the black hole mass derived is rather similar at $0.5 \times 10^8 M_\odot$. This is likely because the outflow velocities are linked to the local Keplerian speed in the P14 model. As a result, as noted previously, one should be cautious when interpreting the black hole mass derived from the outflow model.

The only other report of the black hole mass is by Koss et al. (2017) for single-epoch estimates (which can have large uncertainties) based on the broad Balmer lines. Using the method of Trakhtenbrot & Netzer (2012) with the Hβ line, these authors estimate a mass of $3.8 \times 10^8 M_\odot$. Similarly, using the method of Greene & Ho (2005) for the Hα line they find a mass of $2.4 \times 10^8 M_\odot$. Although slightly less than these, our new value is fully consistent within the uncertainties of the methods used. There is no reported measurement of the stellar velocity dispersion in the literature, making it difficult to place this object on the $M_{BH}-\sigma_*$ relation. However, as has been commonly done, we can use the width of the [O III] line. From the numbers reported by Perez et al. (1989), we can estimate the dispersion as 250 km s⁻¹. This puts IRAS 09149–6206 only a factor of two to three below the relation as defined by Gültekin et al. (2009), which is within the scatter. Although the object is a factor of five below the full relation of McConnell & Ma (2013), this offset is reduced when one considers the relation they find for late-type galaxies. Good agreement is also found if we adopt the correlations for the late-type galaxy sample from Greene et al. (2020). This object therefore does not appear to be unusual in terms of its black hole mass.

We calculated the virial factor, $f_{FWHM} \equiv GM_{BH}/(R_{BLR}v_{FWHM}^2)$, by randomly drawing the BLR model parameters from the posterior parameter space sampled from our fitting procedure. The FWHM of the model line profile was used to calculate the velocity, $v_{FWHM}$. From the Keplerian model, $f_{FWHM} = 0.59^{+0.67}_{-0.26}$, while $f_{FWHM} = 0.31^{+0.25}_{-0.09}$ from the outflow model, again statistically consistent with each other. The typical $f_{FWHM}$, based on calibration against the $M_{BH}-\sigma_*$ relation, is
~1.3 (e.g. Onken et al. 2004; Woo et al. 2010; Ho & Kim 2014). The difference in the virial factor explains most of the difference between our BH mass and those from the single-epoch estimate. Remarkably, Ho & Kim (2014) found $f_{FHWM} \approx 0.5$ for AGNs with pseudobulges. Unfortunately, high resolution imaging is not available to reveal the bulge properties of IRAS 09149–6206. Although the average $f_{FHWM}$ from BLRs with dynamical modelling is about 1 (Williams et al. 2018), the inferred value of $f_{FHWM}$ for individual AGNs shows a wide distribution. In particular, the AGNs with similar BLR structure from RM dynamical modelling, Arp 151 (Pancoast et al. 2014a) and Mrk 142 (Li et al. 2018), also show a comparably low $f_{FHWM}$.

The BLR size of IRAS 09149–6206 is robustly measured from our data. In Fig. 10, we compare the BLR radius of IRAS 09149–6206 and 3C 273 measured with GRAVITY to the radius–luminosity (R–L) relation of the RM results. Du & Wang (2019) provide the most recent collection of 75 Hβ time lags from various reverberation campaigns during the past two decades<fnref n="6" />. The collection primarily includes 41 AGNs from Bentz et al. (2013) and 25 AGNs with high accretion rate from the SEAMBH campaign (Du et al. 2014, 2015, 2016b, 2018b). We also compare with 44 Hβ time lags from the SDSS-RM campaign (Grier et al. 2017b), which is based on a homogeneously selected quasar sample (Shen et al. 2015). The BLR radius and 5100 Å continuum luminosity of IRAS 09149–6206 are 89 light days and $8.32 \times 10^{44}$ erg s⁻¹ (Koss et al. 2017), respectively. Those of 3C 273 are 146 light days (GC18) and $8.43 \times 10^{45}$ erg s⁻¹ (Zhang et al. 2019), respectively. Both IRAS 09149–6206 and 3C 273 show good consistency with the RM results. Moreover, IRAS 09149–6206 is very close to the best-fit relation from Bentz et al. (2013).

Towards the high luminosity end of the R–L relation, the RM results also show large scatter and tend to drop below the best-fit relation from Bentz et al. (2013), which we refer to as the “standard” R–L relation. A study of high accretion rate AGNs finds a deviation from this relation that is primarily driven by accretion rate (Du et al. 2015; Du & Wang 2019). These authors proposed that the size of the BLR is reduced at high accretion rate due to the anisotropic emission of the “slim” accretion disk (Abramowicz et al. 1988). Interestingly, the BLR radii measured from the SDSS-RM campaign also lie mostly below the R–L relation, although the accretion rates of these AGNs are not high by the standard of the SEAMBH AGNs (Grier et al. 2017b; Fonseca Alvarez et al. 2020). With the bolometric luminosity $\sim 2 \times 10^{45}$ erg s⁻¹ (GC20a), the Eddington ratio of IRAS 09149–6206 is $\sim 0.2$. Following Eq. (2) of Du et al. (2015) and adopting $i \approx 21°$ from our best-fit Keplerian model, the dimensionless accretion rate of IRAS 09149–6206 is $\dot{M} \approx 4.1$. This is about the limit beyond which AGNs start to deviate from the R–L relation. With $i \approx 12°$ (GC18), the dimensionless accretion rate of 3C 273 is $\sim 23.8$. The accretion rate estimates are quantitatively consistent with IRAS 09149–6206 being close to the R–L relation, while 3C 273 slightly deviates from it.

In closing, we note that different emission lines may trace regions at different radii of the BLR (e.g. Gaskell & Sparke 1986; Clavel et al. 1991; Peterson et al. 1991; Kaspi et al. 2000; Bentz et al. 2010; Grier et al. 2017b). For example, the Hα lag is expected to be larger than the Hβ lag due to radial stratification and optical depth effects (e.g. Rees et al. 1989; Korista & Goad 2000, 2004; Bottorff et al. 2002; Netzer 2020). Zhang et al. (2019) found the Hβ time lag of 3C 273 is fully consistent with the BLR radius measured by GRAVITY from the Paα

<figure><img src="image_12.png" /><figcaption>Fig. 11. Cartoon illustrating a possible cause for the observed offset between the near-infrared continuum photocentre and the BLR. Dust continuum is assumed to follow a ring that is centred on the BLR and with radius based on the size measurement of GC20a. Brightness variation along this ring causes the offset between BLR and continuum photocentre. Variation in brightness of the hot dust distribution along a ring is supported by the resolved observations of NGC 1068 reported by GRAVITY Collaboration (2020b).</figcaption></figure>

emission, after de-trending the contamination from the jet emission (Li et al. 2020). Considering that Hβ and Paα both come from the $n=4$ level of hydrogen, Wang et al. (2020) argued that the two lines are likely to originate from similar regions of the BLR. They also estimated the possible size difference (~13%) of the Hβ and Paα emission regions based on the difference of their FWHM. The size difference of Hβ and Bry for IRAS 09149–6206 is less clear than that of 3C 273, as the two lines come from different upper levels. However, we find the FWHM of Bry is ~3350 km s⁻¹, which is very close to the Hβ FWHM ~3500 km s⁻¹ (Perez et al. 1989). Following Wang et al. (2020), we estimate a size difference between the Hβ and Bry emitting regions of $\le 10\%$.

## 6. Origin of the spatial offset between BLR and continuum photocentre

Our models (Table 3) place the continuum photocentre outside the bulk of the BLR cloud distribution, at a physical scale ~0.14 pc. While there could be many possible explanations for such an offset, we argue that it is consistent with the simplest explanation: Both BLR and dust are centred on the black hole but there is a modest level of asymmetry in the K-band emission which is arising from hot dust on scales larger than the BLR, near the sublimation radius. An asymmetry might arise from differential brightness between the near and far sides, from clumpiness or irregularities in the emitting dust structure, or possibly also if the edge of a foreground dust lane crosses the line of sight to the nucleus. The first of these is illustrated in the cartoon of Fig. 11, where, inspired by the asymmetric ring-like dust emission of NGC 1068 (GRAVITY Collaboration 2020b) we assume a dust ring with radius based on the K-band size measurement for IRAS 09149–6206 by GC20a. A linear brightness gradient is imposed to reproduce the shift between geometric centre of both BLR and dust ring, and the continuum photocentre. Such

<footnote n="6">⁶ See Table 1 of Du & Wang (2019) for more details.</footnote>
asymmetries will lead to nonzero continuum closure phases, but these remain small for compact sources. For the specific configuration of Fig. 11, closure phases <1° are expected, consistent with the measurements shown in Fig. 5. More generally, in the marginally resolved limit, the maximum closure phase on VLTI triangles is approximately, 1° (FWHM/0.6 mas)³, where the source FWHM is scaled to that measured from FT data for IRAS 09149–6206.

The Fig. 11 scenario is clearly not unique, but one of the simple and plausible ways to create an offset between BLR and continuum photocentre, while staying consistent with all GRAVITY observations. Other plausible ways include a tilted view at non-planar (e.g. bowl-shaped) dust emission. And finally, more exotic explanations for an offset are not ruled out, such as the recoil option discussed in Sect. 4.5.

## 7. Conclusion

With 7.8h on-source integration of GRAVITY, we successfully spatially resolve the broad Bry emission line region of IRAS 09149–6206. This is the second source, following 3C 273 (GC18), for which near-infrared interferometric observations directly constrain the size of the BLR and enable an estimate of the mass of the central black hole. With an improved phase calibration method, the differential phase can be uniformly calibrated to systematic uncertainty ~0.05° for each baseline. This enables us to robustly resolve the BLR of the nearby AGN with the broad Bry line. The main results are summarised as follows.

- We obtain a ~0.5° differential phase signal on two baselines, which is measured from the Bry emission line with a peak flux that is ~6% of the continuum. The differential visibility amplitude of the BLR is ~0.8% above the continuum, indicating that the BLR is much more compact than the continuum emission. The closure phase of the continuum emission is ~0°, consistent with the continuum being only marginally resolved.

- The model-independent reconstruction of photocentres reveals that the BLR is offset to the east of the photocentre of the continuum by ~120 µas. While the offset dominates the differential phase signal, the photocentres display a significant blue-red velocity gradient in a north-south direction, indicating that we are resolving the kinematics of broad Bry emission.

- We model the interferometric data with (1) a simplified BLR model including only clouds on circular orbits and (2) a generalised dynamical model that allows for radial motions, and which is widely used in analysing AGN RM data. Both models provide an adequate fit to the data. We argue against the outflow model because there are several difficulties associated with its physical interpretation and implication, and caution is needed when interpreting the parameters from the fit. Based on the favoured Keplerian model, and with 95% credibility intervals, we report a radius for the BLR of $65_{-39}^{+30}$ µas or $89_{-53}^{+41}$ light days, and a black hole mass $1.1_{-0.84}^{+1.80} \times 10^8 M_\odot$.

- The BLR radii measured by GRAVITY (Bry size for IRAS 09149–6206 and Pa$\alpha$ size for 3C 273) are quantitatively consistent with the radius-luminosity relation based on H$\beta$ RM of AGNs.

*Acknowledgements.* We thank the referees for their careful reading of the manuscript and their suggestions that have helped to improve it. This research has made use of the NASA/IPAC Extragalactic Database (NED), which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration. J.D. was supported in part by NSF grant AST 1909711. A.A. and P.G. were supported by

Fundação para a Ciência e a Tecnologia, with grants reference UIDB/00099/2020 and SFRH/BSAB/142940/2018. SH acknowledges support from the European Research Council via Starting Grant ERC-StG-677117 DUST-IN-THE-WIND. J.S. thanks the helpful discussions with Jian-Min Wang, Pu Du, Yan-Rong Li, and Yu-Yang Songsheng.

## References

Abramowicz, M. A., Czerny, B., Lasota, J. P., & Szuszkiewicz, E. 1988, ApJ, 332, 646

Akaike, H. 1973, in *Proceedings of the 2nd International Symposium on Information Theory*, eds. B. N. Petrov, & F. Csaki (Budapest: Akademiai Kiado), 267

Arsenault, R., Alonso, J., Bonnet, H., et al. 2003, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, eds. P. L. Wizinowich, & D. Bonaccini, Proc. SPIE, 4839, 174

Bailey, J. A. 1998, in Optical Astronomical Instrumentation, ed. S. D'Odorico, SPIE Conf. Ser., 3355, 932

Baskin, A., & Laor, A. 2018, MNRAS, 474, 1970

Bentz, M. C., Walsh, J. L., Barth, A. J., et al. 2010, ApJ, 716, 993

Bentz, M. C., Denney, K. D., Grier, C. J., et al. 2013, ApJ, 767, 149

Blandford, R. D., & McKee, C. F. 1982, ApJ, 255, 419

Blind, N., Eisenhauer, F., Haug, M., et al. 2014, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Proc. SPIE, 9146, 91461U

Boizelle, B. D., Barth, A. J., Walsh, J. L., et al. 2019, ApJ, 881, 10

Booth, C. M., & Schaye, J. 2009, MNRAS, 398, 53

Bottorff, M. C., Baldwin, J. A., Ferland, G. J., Ferguson, J. W., & Korista, K. T. 2002, ApJ, 581, 932

Burtsher, L., Meisenheimer, K., Tristram, K. R. W., et al. 2013, A&A, 558, A149

Buscher, D. F., & Longair, F. B. M. 2015, Practical Optical Interferometry (Cambridge: Cambridge University Press)

Ciddor, P. E. 1996, Appl. Opt., 35, 1566

Clavel, J., Reichert, G. A., Alloin, D., et al. 1991, ApJ, 366, 64

Colavita, M. M., Swain, M. R., Akeson, R. L., Koresko, C. D., & Hill, R. J. 2004, PASP, 116, 876

Davis, T. A. 2014, MNRAS, 443, 911

Davies, R. I., Thomas, J., Genzel, R., et al. 2006, ApJ, 646, 754

De Rosa, G., Fausnaugh, M. M., Grier, C. J., et al. 2018, ApJ, 866, 133

Dehghanian, M., Ferland, G. J., Kriss, G. A., et al. 2020, ApJ, 898, 141

Denney, K. D., Peterson, B. M., Pogge, R. W., et al. 2009, ApJ, 704, L80

Du, P., & Wang, J.-M. 2019, ApJ, 886, 42

Du, P., Hu, C., Lu, K.-X., et al. 2014, ApJ, 782, 45

Du, P., Hu, C., Lu, K.-X., et al. 2015, ApJ, 806, 22

Du, P., Lu, K.-X., Hu, C., et al. 2016a, ApJ, 820, 27

Du, P., Lu, K.-X., Zhang, Z.-X., et al. 2016b, ApJ, 825, 126

Du, P., Brotherton, M. S., Wang, K., et al. 2018a, ApJ, 869, 142

Du, P., Zhang, Z.-X., Wang, K., et al. 2018b, ApJ, 856, 6

Dubois, Y., Peirani, S., Pichon, C., et al. 2016, MNRAS, 463, 3948

Elvis, M. 2000, ApJ, 545, 63

Emmering, R. T., Blandford, R. D., & Shlosman, I. 1992, ApJ, 385, 460

Eracleous, M., & Halpern, J. P. 1994, ApJS, 90, 1

Eracleous, M., & Halpern, J. P. 2003, ApJ, 599, 886

Eracleous, M., Boroson, T. A., Halpern, J. P., & Liu, J. 2012, ApJS, 201, 23

Everett, J. E. 2005, ApJ, 631, 689

Fabian, A. C. 2012, ARA&A, 50, 455

Ferrarese, L., & Ford, H. 2005, Space Sci. Rev., 116, 523

Fonseca Alvarez, G., Trump, J. R., Homayouni, Y., et al. 2020, ApJ, 899, 73

Gaskell, C. M., & Sparke, L. S. 1986, ApJ, 305, 175

Gitton, P. B., Leveque, S. A., Avila, G., & Phan Duc, T. 2004, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, ed. W. A. Traub, Proc. SPIE, 5491, 944

GRAVITY Collaboration (Abuter, R., et al.) 2017, A&A, 602, A94

GRAVITY Collaboration (Sturm, E., et al.) 2018, Nature, 563, 657

GRAVITY Collaboration (Dexter, J., et al.) 2020a, A&A, 635, A92

GRAVITY Collaboration (Pfuhl, O., et al.) 2020b, A&A, 634, A1

Greene, J. E., & Ho, L. C. 2005, ApJ, 630, 122

Greene, J. E., Peng, C. Y., Kim, M., et al. 2010, ApJ, 721, 26

Greene, J. E., Strader, J., & Ho, L. C. 2020, ARA&A, 58, 257

Grier, C. J., Pancoast, A., Barth, A. J., et al. 2017a, ApJ, 849, 146

Grier, C. J., Trump, J. R., Shen, Y., et al. 2017b, ApJ, 851, 21

Gültekin, K., Richstone, D. O., Gebhardt, K., et al. 2009, ApJ, 698, 198

Heckman, T. M., & Best, P. N. 2014, ARA&A, 52, 589

Hicks, E. K. S., & Malkan, M. A. 2008, ApJS, 174, 31

Ho, L. C., & Kim, M. 2014, ApJ, 789, 17

Hopkins, P. F., Hernquist, L., Cox, T. J., & Kereš, D. 2008, ApJS, 175, 356
Horne, K., De Rosa, G., Peterson, B. M., et al. 2020, ApJ, submitted [arXiv:2003.01448]

Hurvich, C. M., & Tsai, C.-L. 1989, Biometrika, 76, 297

Kaspi, S., Smith, P. S., Netzer, H., et al. 2000, ApJ, 533, 631

Keating, S. K., Everett, J. E., Gallagher, S. C., & Deo, R. P. 2012, ApJ, 749, 32

Kishimoto, M., Hönig, S. F., Antonucci, R., et al. 2011, A&A, 536, A78

Kollatschny, W., & Zetzl, M. 2011, Nature, 470, 366

Komossa, S. 2012, Adv. Astron., 2012, 364973

Komossa, S., Zhou, H., & Lu, H. 2008, ApJ, 678, L81

Korista, K. T., & Goad, M. R. 2000, ApJ, 536, 284

Korista, K. T., & Goad, M. R. 2004, ApJ, 606, 749

Kormendy, J., & Ho, L. C. 2013, ARA&A, 51, 511

Koss, M., Trakhtenbrot, B., Ricci, C., et al. 2017, ApJ, 850, 74

Kuo, C. Y., Braatz, J. A., Condon, J. J., et al. 2011, ApJ, 727, 20

Lacour, S., Dembet, R., Abuter, R., et al. 2019, A&A, 624, A99

Lapeyrere, V., Kervella, P., Lacour, S., et al. 2014, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Proc. SPIE, 9146, 91462D

Li, Y.-R., & Wang, J.-M. 2018, MNRAS, 476, L55

Li, Y.-R., Songsheng, Y.-Y., Qiu, J., et al. 2018, ApJ, 869, 137

Li, Y.-R., Zhang, Z.-X., Jin, C., et al. 2020, ApJ, 897, 18

López-Gonzaga, N., Burtscher, L., Tristram, K. R. W., Meisenheimer, K., & Schartmann, M. 2016, A&A, 591, A47

Mangham, S. W., Knigge, C., Matthews, J. H., et al. 2017, MNRAS, 471, 4788

Mangham, S. W., Knigge, C., Williams, P., et al. 2019, MNRAS, 488, 2780

Marconi, A., Maiolino, R., & Petrov, R. G. 2003, Ap&SS, 286, 245

Márquez, I., Durret, F., & Petitjean, P. 1999, A&AS, 135, 83

McConnell, N. J., & Ma, C.-P. 2013, ApJ, 764, 184

Millour, F., Petrov, R. G., Vannier, M., & Kraus, S. 2008, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Proc. SPIE, 7013, 70131G

Murphy, T., Mauch, T., Green, A., et al. 2007, MNRAS, 382, 382

Murray, N., Chiang, J., Grossman, S. A., & Voit, G. M. 1995, ApJ, 451, 498

Netzer, H. 2015, ARA&A, 53, 365

Netzer, H. 2020, MNRAS, 494, 1611

Netzer, H., & Laor, A. 1993, ApJ, 404, L51

Onishi, K., Iguchi, S., Davis, T. A., et al. 2017, MNRAS, 468, 4663

Onken, C. A., Ferrarese, L., Merritt, D., et al. 2004, ApJ, 615, 645

Onken, C. A., Valluri, M., Brown, J. S., et al. 2014, ApJ, 791, 37

Pancoast, A., Brewer, B. J., Treu, T., et al. 2014a, MNRAS, 445, 3073

Pancoast, A., Brewer, B. J., & Treu, T. 2014b, MNRAS, 445, 3055

Pancoast, A., Barth, A. J., Horne, K., et al. 2018, ApJ, 856, 108

Panessa, F., Tarchi, A., Castangia, P., et al. 2015, MNRAS, 447, 1289

Pedregosa, F., Varoquaux, G., Gramfort, A., et al. 2011, J. Mach. Learn. Res., 12, 2825

Pei, L., Fausnaugh, M. M., Barth, A. J., et al. 2017, ApJ, 837, 131

Perez, E., Manchado, A., Pottasch, S. R., & Garcia-Lario, P. 1989, A&A, 215, 262

Peterson, B. M. 1993, PASP, 105, 247

Peterson, B. M. 2014, Space Sci. Rev., 183, 253

Peterson, B. M., Balonek, T. J., Barker, E. S., et al. 1991, ApJ, 368, 119

Peterson, B. M., Ferrarese, L., Gilbert, K. M., et al. 2004, ApJ, 613, 682

Petrov, R. G., Malbet, F., Richichi, A., et al. 2001, C.R. Phys., 2, 67

Planck Collaboration XIII. 2016, A&A, 594, A13

Raimundo, S. I., Pancoast, A., Vestergaard, M., Goad, M. R., & Barth, A. J. 2019, MNRAS, 489, 1899

Raimundo, S. I., Vestergaard, M., Goad, M. R., et al. 2020, MNRAS, 493, 1227

Rakshit, S., Petrov, R. G., Meilland, A., & Hönig, S. F. 2015, MNRAS, 447, 2420

Rees, M. J., Netzer, H., & Ferland, G. J. 1989, ApJ, 347, 640

Saglia, R. P., Opitsch, M., Erwin, P., et al. 2016, ApJ, 818, 47

Schnittman, J. D., & Buonanno, A. 2007, ApJ, 662, L63

Schwarz, G. 1978, Annal. Stat., 6, 461

Shankar, F., Bernardi, M., Richardson, K., et al. 2019, MNRAS, 485, 1278

Shen, Y., Brandt, W. N., Dawson, K. S., et al. 2015, ApJS, 216, 4

Smith, J. E., Robinson, A., Alexander, D. M., et al. 2004, MNRAS, 350, 140

Smith, J. E., Robinson, A., Young, S., Axon, D. J., & Corbett, E. A. 2005, MNRAS, 359, 846

Somerville, R. S., & Davé, R. 2015, ARA&A, 53, 51

Somerville, R. S., Hopkins, P. F., Cox, T. J., Robertson, B. E., & Hernquist, L. 2008, MNRAS, 391, 481

Speagle, J. S. 2020, MNRAS, 493, 3132

Stock, M. R. 2018, Master's Thesis München Technische Universität München

Storchi-Bergmann, T., Schimoia, J. S., Peterson, B. M., et al. 2017, ApJ, 835, 236

Strateva, I. V., Strauss, M. A., Hao, L., et al. 2003, AJ, 126, 1720

Thomas, J., Saglia, R. P., Bender, R., et al. 2004, MNRAS, 353, 391

Trakhtenbrot, B., & Netzer, H. 2012, MNRAS, 427, 3081

van den Bosch, R. C. E., Greene, J. E., Braatz, J. A., Constantin, A., & Kuo, C.-Y. 2016, ApJ, 819, 11

Veron-Cetty, M. P., & Woltjer, L. 1990, A&A, 236, 69

Waisberg, I., Dexter, J., Pfühl, O., et al. 2017, ApJ, 844, 72

Wang, J.-M., Songsheng, Y.-Y., Li, Y.-R., Du, P., & Zhang, Z.-X. 2020, Nat. Astron., 4, 517

Williams, P. R., Pancoast, A., Treu, T., et al. 2018, ApJ, 866, 75

Woo, J.-H., Treu, T., Barth, A. J., et al. 2010, ApJ, 716, 269

Xiao, M., Du, P. Lu, K.-K., et al. 2018, ApJ, 865, L8

Young, S., Axon, D. J., Robinson, A., Hough, J. H., & Smith, J. E. 2007, Nature, 450, 74

Zhang, Z.-X., Du, P., Smith, P. S., et al. 2019, ApJ, 876, 49

Zhu, G., Zaw, I., Blanton, M. R., & Greenhill, L. J. 2011, ApJ, 742, 73

1 Max Planck Institute for Extraterrestrial Physics (MPE), Giessen-bachstr. 1, 85748 Garching, Germany

2 LESIA, Observatoire de Paris, Université PSL, CNRS, Sorbonne Université, Univ. Paris Diderot, Sorbonne Paris Cité, 5 Place Jules Janssen, 92195 Meudon, France

3 I. Institute of Physics, University of Cologne, Zülpicher Straße 77, 50937 Cologne, Germany

4 Departments of Physics and Astronomy, Le Conte Hall, University of California, Berkeley, CA 94720, USA

5 Department of Physics and Astronomy, University of Southampton, Southampton, UK

6 Department of Physics, Kyoto Sangyo University, Kita-ku, Japan

7 Université Côte d'Azur, Observatoire de la Côte d'Azur, CNRS, Laboratoire Lagrange, Nice, France

8 School of Physics and Astronomy, Tel Aviv University, Tel Aviv 69978, Israel

9 Department of Astronomy, The Ohio State University, Columbus, OH, USA

10 Center for Cosmology and AstroParticle Physics, The Ohio State University, Columbus, OH, USA

11 Space Telescope Science Institute, Baltimore, MD, USA

12 Univ. Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France

13 Department of Particle Physics and Astrophysics, Weizmann Institute of Science, Rehovot 76100, Israel

14 Center for Computational Astrophysics, Flatiron Institute, 162 5th Ave., New York, NY 10010, USA

15 European Southern Observatory, Casilla 19001, Santiago 19, Chile

16 European Southern Observatory, Karl-Schwarzschild-Str. 2, 85748 Garching, Germany

17 Sterrewacht Leiden, Leiden University, Postbus 9513, 2300 RA Leiden, The Netherlands

18 Max Planck Institute for Radio Astronomy, Auf dem Hügel 69, 53121 Bonn, Germany

19 Universidade de Lisboa – Faculdade de Ciências, Campo Grande, 1749-016 Lisboa, Portugal

20 Faculdade de Engenharia, Universidade do Porto, Rua Dr. Roberto Frias, 4200-465 Porto, Portugal

21 CENTRA – Centro de Astrofísica e Gravitação, IST, Universidade de Lisboa, 1049-001 Lisboa, Portugal

22 Max Planck Institute for Astronomy, Königstuhl 17, 69117 Heidelberg, Germany

23 Instituto de Astrofísica de Canarias (IAC), 38200 La Laguna, Tenerife, Spain

24 Department of Astrophysical & Planetary Sciences, JILA, University of Colorado, Duane Physics Bldg., 2000 Colorado Ave, Boulder, CO 80309, USA

25 Research School of Astronomy and Astrophysics, Australian National University, Canberra, ACT 2611, Australia
# Appendix A: Calibration of the pipeline reduced differential phase

<figure><img src="image_13.png" /><figcaption>**Fig. A.1.** Variation of the differential phase of GRAVITY data. The solid and dashed curves display the differential phases of two observations of calibrator stars, whose differential phase should be consistent with 0°. The vertical bar in the lower left corner of each panel indicates the typical uncertainty, enlarged by a factor of three for clarity.</figcaption></figure>

As shown in Fig. A.1, the pipeline reduced differential phases of the calibrator stars show, at different times, considerable non-zero variation as a function of wavelength. To investigate this, we collected from the ESO archive<fnref n="7" /> on-axis observations of calibrator stars with the same spectral resolution and polarisation configuration. Before the exchange of the GRAVITY science detector during the intervention in October 2019, we found 122 good files<fnref n="8" />. We found 23 good files after the intervention, up to March 2020. The calibrator data reveal a clear systematic phase variation. A principle component analysis<fnref n="9" /> (PCA) shows that there is only one variable component dominating the whole variation of the differential phase in all of the baselines. We refer to this variable component as C1 (Fig. A.2), which stays the same both before and after the intervention of GRAVITY.

Using PCA, we also find the variability is introduced by the dispersion of the air in the non-vacuum delay line. The amplitude of the C1 component is well correlated with the light path difference in the delay line (Fig. A.3). In addition, we also find stable instrumental phase features, that are likely to come from the Fabry-Pérot effect of the elements in the cryostat of GRAVITY. As shown in Fig. A.4, these features changed after the GRAVITY intervention, which is expected since the cryostat was opened and many elements were adjusted.

We therefore flatten the pipeline reduced AGN data by fitting and subtracting an instrumental phase model. This includes the stable phase feature as well as the C1 component, the amplitude of which is determined from the fit. The data before and after the intervention are fit separately using the appropriate stable phase feature. We measure the systematic uncertainty of the differential phase after removing the instrumental features using the calibrator data. As Fig. A.5 shows, the calibrator differential phases, observed after the intervention, are reduced and stacked in the same way as the AGN data. The stacked differential phases are fully consistent with zero, and have an rms close to, or better than, 0.05°. The phase at ≤2.05 µm and ≥2.40 µm suffers strongly from the absorption bands of carbon dioxide and water in the atmosphere, therefore the phase variation increases at those wavelengths. For calibrator data taken before the intervention, we find the rms of the stacked differential phases is quantitatively similar to that in Fig. A.5. We therefore conclude that the systematic uncertainty of the differential phase is 0.05°. We also find that the accuracy of the flattening is insensitive to the choice of the wavelength mask (e.g. across an emission line), as long as the mask is only a reasonable fraction of the entire spectrum. This is because we only have one free parameter in the fit, the amplitude of the C1 component, which has an unvarying profile.

The pipeline uses a third-order polynomial to fit the phase of the complex visibility ($R + iI$) of each DIT to derive the self-reference phase ($\phi_{self}$). The corrected visibility ($R' + iI'$) of each DIT is,

<footnote n="7">⁷ http://archive.eso.org/wdb/wdb/eso/gravity/form</footnote>

<footnote n="8">⁸ We require the fringe tracking ratio of the data better than 99.9%.</footnote>

<footnote n="9">⁹ We adopted `sklearn.decomposition.PCA` function (Pedregosa et al. 2011).</footnote>
<figure><img src="image_14.png" /><figcaption>Fig. A.2. First principle components (C1) from the PCA of all of the baselines (grey curves), before and after the GRAVITY intervention, separately, are very consistent with each other. We fit them together with an eighth-order polynomial (red curve). The strong variations beyond 2.4 µm are due to the resonances of water vapour (Colavita et al. 2004).</figcaption></figure>

$$ 
\begin{aligned}
R' &= \cos(\phi_{\text{self}})R - \sin(\phi_{\text{self}})I, \\
I' &= \sin(\phi_{\text{self}})R + \cos(\phi_{\text{self}})I.
\end{aligned} 
\qquad (A.1) 
$$

It would be very convenient if we could avoid the a posteriori flattening by adopting a better self-reference phase model that takes into account the dispersion of the air. Unfortunately, we find that the models of the refractive index of the air are not accurate enough, in terms of their functional form, for our accuracy requirement (Ciddor 1996). We also tested whether a higher order polynomial fit can remove the effect of the dispersion of the air. However, even fitting with a seventh order polynomial still leaves a residual larger than 0.1°, and using such a high order polynomial is risky because it may also fit the scientific phase signal. Therefore, we prefer to apply the separate phase flattening method described above after the default pipeline data reduction.

<figure><img src="image_15.png" /><figcaption>Fig. A.3. Amplitude of the C1 component as a function of the length difference between the non-vacuum delay lines in each baseline.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>**Fig. A.4.** Stable phase features for each baseline, before (solid) and after (dashed) the intervention of GRAVITY.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Fig. A.5.** Differential phase of the calibrator stars, flattened and stacked in the same way as the AGN. The 23 calibrator files after the GRAVITY intervention in late 2019 are used, which is quantitatively close to the results using the 122 files before the intervention. The rms of the phase at 2.05 to 2.40 µm is reported on the upper right corner of each panel.</figcaption></figure>
## Appendix B: The differential phase of AGN BLRs

In this section, we derive the generalised differential phase of the AGN including the continuum phase, under the assumption that the BLR and continuum are marginally resolved. The complex visibility is defined as,

$$ \tilde{V} = V(\lambda) e^{i\phi_{\lambda}} = \frac{\iint I_{\lambda}(\boldsymbol{\sigma}) e^{-2\pi i \boldsymbol{u} \cdot \boldsymbol{\sigma}} \, d\boldsymbol{l} \, d\boldsymbol{m}}{\iint I_{\lambda}(\boldsymbol{\sigma}) \, d\boldsymbol{l} \, d\boldsymbol{m}}, \quad (B.1) $$

where $V(\lambda)$ and $\phi_{\lambda}$ denote the visibility amplitude and phase as functions of wavelength ($\lambda$); $I_{\lambda}(\boldsymbol{\sigma})$ is the source intensity distribution; $\boldsymbol{\sigma} = (l, m)$ coordinate on the sky; $\boldsymbol{u} = \mathbf{B}/\lambda = (u, v)$ is the baseline vector. When the source is only marginally resolved, namely $2\pi \boldsymbol{u} \cdot \boldsymbol{\sigma} \ll 1$, we have $e^{-2\pi i \boldsymbol{u} \cdot \boldsymbol{\sigma}} \approx 1 - 2\pi i \boldsymbol{u} \cdot \boldsymbol{\sigma}$ and $e^{i\phi_{\lambda}} \approx 1 + i\phi_{\lambda}$. Therefore, from Eq. (B.1) and setting $V(\lambda) \approx 1$,

$$ \phi_{\lambda} \approx -2\pi \boldsymbol{u} \cdot \iint I_{\lambda}(\boldsymbol{\sigma}) \boldsymbol{\sigma} \, d\boldsymbol{l} \, d\boldsymbol{m} = -2\pi \boldsymbol{u} \cdot \boldsymbol{x}_{\lambda}, \quad (B.2) $$

where $\boldsymbol{x}_{\lambda}$ is defined as the photocentre of the source at the wavelength $\lambda$. The phase is proportional to the photocentre of the source projected onto the baseline. For three baselines that form a closed triangle, meaning $\boldsymbol{u}_1 + \boldsymbol{u}_2 + \boldsymbol{u}_3 = 0$, their closure phase is naturally, $\phi_{1,\lambda} + \phi_{2,\lambda} + \phi_{3,\lambda} = -2\pi(\boldsymbol{u}_1 + \boldsymbol{u}_2 + \boldsymbol{u}_3) \cdot \boldsymbol{x}_{\lambda} = 0$.

The observed AGN emission consists of two components, the continuum emission from the hot dust and the line emission from the BLR. Therefore, the observed complex visibility is,

$$ \tilde{V} = \frac{f_{c,\lambda} \tilde{V}_c + f_{BLR,\lambda} \tilde{V}_{BLR}}{f_{c,\lambda} + f_{BLR,\lambda}}, \quad (B.3) $$

where $f_{c,\lambda}$ and $\tilde{V}_c$ are the spectral flux and complex visibility of the continuum emission as functions of the wavelength, while $f_{BLR,\lambda}$ and $\tilde{V}_{BLR}$ are the flux and complex visibility of the BLR. Under the marginally resolved assumption, the differential phase is

$$ \Delta\phi_{\lambda} = \phi_{\lambda} - \phi_{c} = \frac{f_{\lambda}}{1 + f_{\lambda}} (\phi_{BLR,\lambda} - \phi_{c}), \quad (B.4) $$

where $\phi_c$ is the phase of the continuum, which is not expected to vary with different spectral channels; $\phi_{BLR,\lambda}$ is the phase of the BLR as a function of wavelength; and $f_{\lambda} \equiv f_{BLR,\lambda}/f_{c,\lambda}$ is the emission line flux normalised by the continuum. Taking Eq. (B.2) for the continuum and BLR separately into Eq. (B.4), we obtain

$$ \Delta\phi_{\lambda} = -2\pi \frac{f_{\lambda}}{1 + f_{\lambda}} \boldsymbol{u} \cdot (\boldsymbol{x}_{BLR,\lambda} - \boldsymbol{x}_{c}), \quad (B.5) $$

where $\boldsymbol{x}_{BLR,\lambda}$ and $\boldsymbol{x}_c$ are the photocentres of the BLR and the continuum.

When we use the BLR model to fit the data, we calculate the differential phase based on Eq. (B.4). The BLR clouds contribute to $\phi_{BLR}$, while the offset of the BLR contributes to the continuum phase,

$$ \Delta\phi_{\lambda} = \frac{f_{\lambda}}{1 + f_{\lambda}} [\phi_{BLR,\lambda} - 2\pi(ux_{o} + vy_{o})], \quad (B.6) $$

where $\{x_o, y_o\}$ is the vector of the offset of the BLR with respect to the photocentre of the continuum.

## Appendix C: Bayesian inference

We infer the optimal model parameters ($\mathbf{\Theta}$) based on the interferometric data ($D$) and our prior knowledge of the source and the model ($I$), according to Bayes theorem,

$$ p(\mathbf{\Theta}|D, I) = \frac{p(\mathbf{\Theta}|I)p(D|\mathbf{\Theta}, I)}{p(D|I)}, \quad (C.1) $$

where $p(\mathbf{\Theta}|D, I)$ is the posterior probability density function of the model parameters. The prior, $p(\mathbf{\Theta}|I)$, is provided by our prior knowledge about the probability distribution of the model parameters. The evidence, $p(D|I)$, is useful to compare across different models using the Bayes factor,

$$ K = \frac{p(D|I_{\text{model 2}})}{p(D|I_{\text{model 1}})}, \quad (C.2) $$

assuming the prior knowledge of model 1 and model 2 are equivalent. The likelihood function, $\mathcal{L} = p(D|\mathbf{\Theta}, I)$ is defined assuming a Gaussian probability distribution,

$$ \ln \mathcal{L} = - \frac{1}{2} \sum_{i}^{n} \left( \frac{(f_i - \tilde{f}_i(\mathbf{\Theta}))^2}{\sigma_{f,i}^2} + \ln(2\pi\sigma_{f,i}^2) \right) \quad (C.3) $$

$$ - \frac{1}{2} \sum_{i}^{n} \left( \frac{(\phi_i - \tilde{\phi}_i(\mathbf{\Theta}))^2}{\sigma_{\phi,i}^2} + \ln(2\pi\sigma_{\phi,i}^2) \right), \quad (C.4) $$

where $f_i$ and $\sigma_{f,i}$ are the normalised line flux and its uncertainty, respectively, at the $i$th spectral channel; $\phi_i$ and $\sigma_{\phi,i}$ are the differential phase and its uncertainty, respectively$^{10}$ and $\tilde{f}_i(\mathbf{\Theta})$ and $\tilde{\phi}_i(\mathbf{\Theta})$ are the normalised line flux and differential phase from the BLR model.

We sample the posterior distribution and calculate the Bayesian evidence for each of the two models using the nested sampling method implemented in the Python package, `dynesty` (Speagle 2020). The nested sampling method has been extensively used for RM studies of the BLR (e.g. Pancoast et al. 2014a; Grier et al. 2017a; Li & Wang 2018; Li et al. 2018, 2020; Williams et al. 2018; Raimundo et al. 2019, 2020; Wang et al. 2020), due to its efficiency in sampling the posterior even when there are complex multi-modal structures. The best-fit model parameters and the uncertainties (Table 3) are estimated, respectively, from the maximum a posteriori value and the 95% (2$\sigma$) credible interval of the posterior probability density distributions. The Bayesian evidence is used to compare the probability of the two models with the Bayes factor. We also calculate the $\chi^2$, Bayesian information criterion (BIC, Schwarz 1978), and Akaike information criterion (AIC, Akaike 1973; Hurvich & Tsai 1989),

$$ \chi_r^2 = \sum_i^N \left( \frac{(f_i - \tilde{f}_i(\mathbf{\Theta}))^2}{\sigma_{f,i}^2} + \frac{(\phi_i - \tilde{\phi}_i(\mathbf{\Theta}))^2}{\sigma_{\phi,i}^2} \right) / (N - k), \quad (C.5) $$

$$ \text{BIC} = k \ln N - 2 \ln \mathcal{L}_{\text{max}}, \quad (C.6) $$

$$ \text{AIC} = 2k - 2 \ln \mathcal{L}_{\text{max}} + \frac{2k(k+1)}{N-k-1}, \quad (C.7) $$

where $N$ is the number of data points, $k$ is the number of free parameters, and $\mathcal{L}_{\text{max}}$ is the maximum likelihood based on the samples from the fit.

<footnote>$^{10}$ For simplicity, and with reference to Eq. (B.4), we write $\phi$ rather than $\Delta\phi$ here.</footnote>
**Appendix D: Plots of the BLR model fits and results**

The fits to the flux spectra and the differential phases for each baseline and angular *uv* bin are displayed in Figs. D.1 and D.3 for the Keplerian model and the outflow model respectively. The

models were fit to the 5 bins of averaged differential phase based on *uv* coverage (Fig. 2) in order to minimise smearing and ensuring the signal is still visible. The posterior probability density distributions of the physical parameters that are important for our discussion are plotted in Figs. D.2 and D.4.

<figure><img src="image_18.png" /><figcaption>**Fig. D.1.** Keplerian model best-fit spectrum and phases. The *first row* displays the same normalised spectrum. The black solid line is the best-fit model. The rows are the differential phases of different baselines. Each column displays the phases of each *uv* bin. In each panel, the solid line is the total phase signal, which sums up the differential phase of the BLR (dashed line) and the continuum phase (dotted line) due to the offset between the BLR and the photocentre of the continuum.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>**Fig. D.2.** Posterior probability distribution of selected parameters from the Keplerian model. The blue lines and the cross points represent the maximum a posteriori.</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>Fig. D.3. Outflow model best-fit spectrum and phases. The arrangements and symbols are the same as Fig. D.1.</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>**Fig. D.4.** Posterior probability distribution of selected parameters from the outflow model. The blue lines and the cross points represent the maximum a posteriori.</figcaption></figure>

## Appendix E: Additional tests of the BLR modelling

In order to reduce the number of parameters, we considered a simplified version of the P14 model in which we fix $\sigma_{\rho,\text{circ}} = \sigma_{\Theta,\text{circ}} = 0$, so that bound orbits are always circular, and $\sigma_{\rho,\text{radial}} = \sigma_{\Theta,\text{radial}} = 0$ so that orbits with radial motion are always at the same location on the phase space ellipse defined by Eq. (6). In addition we set $\sigma_{\text{turb}} = 0$ so that there is no additional macroturbulence. The outcome is that the inferred parameters of this simplified model are entirely consistent with the full P14 model. This and further tests suggest that these five technical parameters ($\sigma_{\rho,\text{circ}}, \sigma_{\Theta,\text{circ}}, \sigma_{\rho,\text{radial}}, \sigma_{\Theta,\text{radial}},$ and $\sigma_{\text{turb}}$) are less important in terms of defining the model. This is consistent with the fact that these parameters are largely ignored by Pancoast et al. (2014a) and other works that use the model to fit RM data. In the main text here, we therefore also do not discuss the specific values of these technical parameters.

We interpret the dominant differential phase signal as the continuum phase, which is produced by the offset between the BLR and the photocentre of the continuum emission. However, in Sect. 4.3 on the outflow model, we show that the BLR alone is also able to produce an asymmetric differential phase with an all positive signal shape. It is therefore interesting to test whether the outflow model is able to fit the data without the need for a continuum phase, so that the BLR would lie at the centre of the continuum emission. We are indeed able to fit the differential phase data reasonably well. However, the inferred radius of the BLR is very large, with $R_{\text{BLR}} \approx 200$ µas. This is easily understood because the best-fit dynamical model in Sect. 4.1 has a mean radius $R_{\text{BLR}} \approx 60$ µas, which corresponds to a phase signal $\sim 0.15^\circ$ shown in Fig. 8b. When the continuum phase is fixed to be zero, the BLR size has to increase by a factor of $\sim 3$ in order to fit the much larger $\sim 0.5^\circ$ phase signal. However, in GC20a we found that the size of the continuum emission
for IRAS 09149-6206, measured as a Gaussian FWHM, is only 0.54-0.64 mas. Thus its BLR radius should be <100 µas, since it must be significantly smaller than the continuum. Similarly, following the method presented in Sect. 4.4, we find the differential visibility amplitude of the large BLR model to be about a factor of two lower than the other two models shown in Fig. 4. Therefore, the differential visibility amplitude data strongly disfavour the large BLR model.

We also attempted to include the differential visibility amplitude data in the fit. However, as the differential visibility amplitude is very sensitive to the relative size between BLR and the continuum, the inferred BLR size strongly depends on the assumed FWHM of the continuum emission. Using FWHM = 0.6 mas, our inferred BLR sizes are only 21$^{+22}_{-9}$ and 30$^{+24}_{-18}$ µas for the Keplerian and outflow models, respectively. GC20a estimated the FWHM of the continuum emission to be 0.54 ± 0.05 mas by directly fitting the visibility amplitude of the fringe tracker data (the “FT size” in their Table 2). From the differential visibility amplitude of the science channel, they also obtained a size of 0.64 ± 0.06 µas (the “SC size”). Under the marginally resolved limit, the latter method only yields the correct continuum size when the BLR is exactly a point source; for better resolved sources, the derived quantity is the quadrature difference between the continuum emission size and the BLR size (Waisberg et al. 2017). Therefore, we should expect the SC size to be slightly smaller than the FT size. This is true for 3C 273 and PDS 456 but not for IRAS 09149-6206. Although the problem can be easily explained by the uncertainty, this indicates that the continuum size of IRAS 09149-6206 is still quite uncertain.

Therefore, we exclude the differential visibility amplitude data from our primary BLR fittings and use it only as a consistency check.

For our last test, we check whether the input parameters can be recovered by fitting mock data generated by the Keplerian model (mock 1) and the outflow model (mock 2). We also test whether the Bayes factor, BIC, and AIC provide a reliable judgment on the model fits. We use the best-fit parameters of the Keplerian model and the outflow model to generate mock data, adopting Gaussian noise based on the measured uncertainty. In the next step, we fit both sets of mock data with both the Keplerian model and the full P14 model. The input parameters are usually recovered if the same model is used to generate and to fit the mock data. Focusing on the mock 2 data, the best-fit Keplerian model results are similar to those in Figs. D.1 and 7. The averaged phase, after subtracting the best-fit continuum phase, still shows clear S-shape profiles. To compare the two models fitting the mock 2 data, the Bayes factor and AIC prefer the P14 model, while the BIC incorrectly selects the Keplerian model. For the fits to mock 1 data, the P14 model can naturally fit the data and provide the parameters consistent with the input, which is as expected. The Bayes factor determines the Keplerian model and the outflow model are equivalent. The AIC incorrectly prefers the outflow model, while the BIC prefers the Keplerian model. In a nutshell, these tests demonstrate that: (1) our models infer the parameters self-consistently when the mock data are generated from them and (2) The Bayes factor provides a reasonable judgement on the preferred model.