BASS XXIX: The near-infrared view of the BLR: the effects of obscuration in BLR characterisation

Virial black hole mass ($M_{BH}$) determination directly involves knowing the broad line region (BLR) clouds velocity distribution, their distance from the central supermassive black hole ($R_{BLR}$) and the virial factor ($f$). Understanding whether biases arise in $M_{BH}$ estimation with increasing obscuration is possible only by studying a large (N$>$100) statistical sample of obscuration unbiased (hard) X-ray selected active galactic nuclei (AGN) in the rest-frame near-infrared (0.8-2.5$\mu$m) since it penetrates deeper into the BLR than the optical. We present a detailed analysis of 65 local BAT-selected Seyfert galaxies observed with Magellan/FIRE. Adding these to the near-infrared BAT AGN spectroscopic survey (BASS) database, we study a total of 314 unique near-infrared spectra. While the FWHMs of H$\alpha$ and near-infrared broad lines (He\textsc{i}, Pa$\beta$, Pa$\alpha$) remain unbiased to either BLR extinction or X-ray obscuration, the H$\alpha$ broad line luminosity is suppressed when $N_H\gtrsim10^{21}$ cm$^{-2}$, systematically underestimating $M_{BH}$ by $0.23-0.46$ dex. Near-infrared line luminosities should be preferred to H$\alpha$ until $N_H<10^{22}$ cm$^{-2}$, while at higher obscuration a less biased $R_{BLR}$ proxy should be adopted. We estimate $f$ for Seyfert 1 and 2 using two obscuration-unbiased $M_{BH}$ measurements, i.e. the stellar velocity dispersion and a BH mass prescription based on near-infrared and X-ray, and find that the virial factors do not depend on redshift or obscuration, but for some broad lines show a mild anti-correlation with $M_{BH}$. Our results show the critical impact obscuration can have on BLR characterization and the importance of the near-infrared and X-rays for a less biased view of the BLR.


INTRODUCTION
Supermassive black holes (SMBHs, with black hole masses M BH ∼ 10 5 − 10 9 M ) are ubiquitous in the local Universe, lurking in the spheroid of almost all local galaxies (Kormendy & Richstone 1995). During active accretion phases, the SMBH is no longer dormant but shines as an active galactic nucleus (AGN), due to a surrounding accretion disk of matter, which releases gravitational energy as it infalls toward the central dark attractor. The ultraviolet emission from the inner accretion disk photoionizes nearby clouds located in the broad line region (BLR). Under the hypothesis of a virialized BLR, whose dynamics are dominated by the central SMBH, the M BH can be simply determined from the velocity ∆V BLR of the emitting gas clouds located at distance R BLR in the BLR as M BH = G −1 ∆V 2 BLR R BLR , with G being the gravitational constant. To model the unknown emission-weighted 1 geometry and dynamics of the BLR, the observed width ∆W obs (either the full-width-at-half maximum FWHM or the second moment of the line profile, i.e., the line dispersion σ line ) of a doppler-broadened photoionized element at distance R BLR from the SMBH is used as a tracer of the true velocity in the BLR and a correction factor f known as the virial factor is introduced (Onken et al. 2004): where M vir is the so-called virial product or virial mass. Since it is in practice impossible to spatially resolve the BLR for statistically sized samples, time-resolved observations substitute for spatially-resolved information to estimate the BLR radius, adopting the so-called reverberation-mapping technique (RM; Blandford & McKee 1982). This is gradually changing due to campaigns with the VLTI/GRAVITY instrument, which has opened the path to spatially-resolved observations of the BLR around nearby AGN, allowing high angular resolution (< 0.1 milliarcsec) spectral-spatial interferometric observations of the Paα (in 3C 273, see, Gravity Collaboration et al. 2018) and Brγ (in IRAS 09149-6206, see, Gravity Collaboration et al. 2020) line. However, this powerful technique remains limited to a small sample of AGN. Extensive RM campaigns have found that the radius of the BLR is linked to the AGN luminosity, R BLR ∝ L α AGN (Bentz et al. 2006(Bentz et al. , 2009(Bentz et al. , 2013, where the slope α is consistent with expectations from photoionization (α 0.5). The AGN continuum and the broad-line luminosity (see, e.g., Shen 2013) have been both used as a proxy for the BLR radius, which has allowed an efficient calibration for single-epoch (SE) BH mass estimation.
The virial factor f has been directly inferred only for a limited subsample (<20) of RM AGN with sufficient high quality data available (Pancoast et al. 2014Grier et al. 2017;Williams et al. 2018;Li et al. 2018).
Hence, previous studies have often adopted an ensemble virial factor f that is determined using the M BH − σ relation observed in local samples of quiescent galaxies with dynamically-based BH massses Ho & Kim 2014;Batiste et al. 2017;Yu et al. 2019) M BH,σ = f M vir . (2) However it is still unclear whether the M BH − σ relation is universally followed by all types of galaxies: e.g., barred/unbarred hosts (Graham 2008), early/late type hosts (McConnell & Ma 2013;Sahu et al. 2019), elliptical and classical-/pseudo-bulges (Kormendy & Ho 2013;Saglia et al. 2016;de Nicola et al. 2019). It is also unsettled whether AGN should follow the scaling relations determined by quiescent galaxies Ricci et al. 2017c;Shankar et al. 2019), since the methods used to measure the BH masses in active galaxies are not expected to suffer from the resolution-dependent bias that instead affect dynamical-based BH mass estimates (Bernardi et al. 2007;Shankar et al. 2016). Moreover, the f -factor could change on an object-by-object basis if it depends on some AGN properties, such as the bolometric luminosity, Eddington ratio (λ Edd ), M BH , obscuration or line-of-sight inclination angle θ. The only statistically sound correlation found so far is between f and θ, although this is based on a limited number of RM objects with directly inferred f . This correlation is consistent with expectations of the BLR being a thick disk with clouds moving in a combination of elliptical and inflowing motions (Williams et al. 2018). This f − θ relation is corroborated by statistical studies that used a variety of f -independent BH mass measurements to infer the virial factor, such as the bulge luminosity based BH mass (Decarli et al. 2008), the stellar velocity dispersion σ based BH mass (Shen & Ho 2014) and the accretion disk based BH mass (Mejía-Restrepo et al. 2018). These studies only focused on optical broad-line AGN, for which optical virial-based M BH estimates were available.
While the use of the Hα emission is common practice to derive M BH in statistical samples to study the demography and evolution of the AGN population, it is unclear whether the Hα is completely reliable, particularly in so-called Sy 1.9 (Osterbrock 1981), where the level of extinction due to dust is more relevant than what is usually experienced in optical broad line Seyferts. Indeed Reines & Volonteri (2015), using a sample of 262 broad-line AGN in the nearby Universe (z < 0.055) with Hα broad-lines measured from the Sloan Digital Sky Survey (SDSS) DR8, find that AGN-hosts with SE Hα-based BH masses define a separate M BH − M relation, with a slope similar to that of early-type galaxies with dynamically detected BHs but with a normalization ∼ 1.2 dex lower. Similarly, Koss et al. (2017) show that the SE Hα-based BH masses in Sy 1.9 are undermassive than what is expected from the M BH − σ relation of elliptical/classical-bulges, with the BH mass deviation being more extreme in sources with broad Hα equivalent width EW < 50Å. Additionally, Caglar et al. (2020) find an offset of ∼ 0.6 dex between the SE Hαbased BH masses and those based on the stellar velocity dispersion in a sample of 19 partially obscured local hard X-ray selected Seyferts from the LLAMA sample. This discrepancy is reduced in the LLAMA sample only after accounting for optical extinction, in the Hα measurement, and galaxy rotation, in the σ estimate (Caglar et al. 2020). These results can be either explained with the fact that BH-host scaling relations should be different in active vs inactive galaxies, or that in some cases the Hα-based BH masses could be biased low in presence of extinction (or with a combination of these two effects).
To this end, the rest-frame near-infrared (NIR, 0.8-2.5µm) band allows to deeper probe the physical condition of the BLR gas, being at least a factor of 10 less affected by dust extinction than the rest-frame optical emission (Goodrich et al. 1994;Veilleux et al. 1997;Veilleux 2002). However, NIR ground-based spectroscopic observations are more complex and timeconsuming than optical spectroscopy, due to the lower atmospheric transmission that reduces the observable windows, bright sky background and strong OH sky line emissions. For these reasons, local AGN samples studied so far have remained limited to few objects, usually 50 (Riffel et al. 2006;Glikman et al. 2006;Landt et al. 2008;Mason et al. 2015;Riffel et al. 2015;Onori et al. 2017a).
Determining the BLR properties and f -factors for a less biased AGN sample is of paramount importance to assess the uncertainties and systematics in M BH measurement in active galaxies as a function of obscuration. This is particularly important for Sy 1.9s, where the only optical broad line available is the Hα, and Sy 2s, where sometimes so-called hidden BLR are found in the NIR (∼30% of cases, see, e.g., Veilleux et al. 1997;Riffel et al. 2006;Cai et al. 2010;Smith et al. 2014;Lamperti et al. 2017;Onori et al. 2017a). Such an investigation is only possible by constructing an obscuration unbiased AGN sample, by means of hard X-ray (> 10 keV) AGN selection that is almost unaffected by interviening obscuring material, at least up to N H ∼ 10 23.5 − 10 24 cm −2 Koss et al. 2016), due to the decline of the photoelectric cross-section with increasing photon energy. A sensitive all-sky survey in the ultra-hard X-ray band (14-195 keV), such as the one carried out by the the Burst Alert Telescope (BAT, Barthelmy 2000) on the Neil Gehrels Swift Observatory (Swift/BAT), coupled with NIR spectroscopy and ancillary optical spectroscopic information, is the ideal database to quantify the effects of obscuration on BLR characterisation and thus BH mass estimation. The BAT AGN Spectroscopic Survey (BASS) 2 has for the first time increased the sample size of AGN surveyed with NIR spectroscopy to more than 100 beginning with the DR1 (Lamperti et al. 2017).  (Oh et al. 2018). Colors and symbols according to the legend. Right: histogram of the hard X-ray luminosities.
In this work, we present NIR Magellan spectra of 65 local Seyferts, selected from the 70-month Swift/BAT catalog (Baumgartner et al. 2013), as part of the BASS survey. We complement our near-infared Magellan sample with the NIR BASS database, NIR DR1 (Lamperti et al. 2017), DR2 (den Brok et al., submitted) 3 , to construct the largest sample of local Seyferts with NIR and X-ray spectral information available to date, for a total of 314 unique NIR spectra. We further complement the NIR analysis with optical spectral information on the broad component of the Hα from the optical BASS DR2 (Mejía-Restrepo et al., submitted), to quantify and 2 https://www.bass-survey.com/ 3 We note that the NIR BASS DR2 (den Brok et al., submitted) does not contain the spectral measurements already published in the NIR BASS DR1 (Lamperti et al. 2017) since the data analysis is very consistent between the two different BASS data releases.
compare the BLR characterisation in both optical and NIR. We finally make use of the optical stellar velocity dispersion measurements available in BASS DR2 (Koss et al., in prep.; Caglar et al., in prep.), to infer the individual f factors of our sample. The work is organized as follows: Sect. 2 presents the Magellan data selection and reduction; Sect. 3 describes the Magellan spectral fitting of the most important NIR emission lines; Sect. 4 is devoted to the descriptions of two independent BH mass measurement methods adopted to derive the virial factors f . Results are described in Sect. 5, where we investigate the fraction of hidden BLRs detected in Sy 1.8-1.9-2 (Sect. 5.1), the effects of X-ray obscuration and BLR extinction on the BLR velocity and radius tracers derived from optical (i.e., Hα) and NIR emission lines (He iλ10830Å, Paβ λ12821Å, Paα λ18756Å) in Sy 1-1.9 (Sect. 5.2) and in Sy 1-2 (Sect. 5.4). In Sect. 5.3 we explore the connection between the material repsonsible for BLR extinction and the one absorbing the X-rays. The virial factors of our sample are derived in Sect. 5.5, where we test whether f depends on some parameters, e.g., z, N H and BH mass. Section 6 is devoted to discussions and conclusions, while Sect. 7 briefly summarises our main results. We adopt the concordance cosmological model, Ω M = 0.3, Ω Λ = 0.7 and h = 0.7.

DATA
Here we present the NIR spectroscopic data (PI: E. Treister, F. Ricci, M. Baloković) 4 obtained at Magellan using the Folded-port InfraRed Echellette (FIRE; Simcoe et al. 2008). The Magellan/FIRE sample was selected from the hard X-ray (14-195 keV) 70-month catalog (Baumgartner et al. 2013) without near-infared coverage in the BASS DR1 (Lamperti et al. 2017), as part of an effort within the BASS collaboration to obtain NIR coverage for an obscuration unbiased census of local accreting SMBHs. The 70-month catalog (Baumgartner et al. 2013) listed 838 AGN, 102/838 were targeted in the BASS DR1 (Lamperti et al. 2017) and 118/838 were observed as part of the DR2 (den Brok et al., submitted) which selected as well additional 50/1016 AGN from the latest 105-month source catalog (Oh et al. 2018). The combined sample totals 314 hard X-ray selected AGN with unique NIR spectra. The FIRE sample was chosen to target the more obscured sources, to estimate the M BH also in obscured Seyfert class AGN (i.e., Sy 1.8-1.9-2). The FIRE sample is composed of 65 targets selected at z 0.2, divided into 52 obscured AGN (i.e., Sy 1.8-1.9-2) and 13 optical broad-line AGN (i.e., Sy 1-1.2-1.5), whose Seyfert classification is defined according to the Osterbrock (1981) standard criteria, using optical spectra collected by BASS (e.g., BASS optical DR2, Mejía-Restrepo et al., submitted). Figure 1  (3) optical Seyfert classification, as defined by Osterbrock (1981); (4) observation date; (5) exposure time; (6) airmass at the midpoint of the observation; (7) J-band Vega mag (from 2MASS extended, e.g., Jarrett et al. 2000); (8) redshift from the [O iii] BASS DR2; (9-10) slit width and aperture of the spectral extraction, in arcsec and kpc; (11) intrinsic 2-10 keV luminosity from Ricci et al. (2017a). (This table is available in its entirety in a machine-readable form in the online journal. A part is shown as guidance for the reader regarding its content.) * Readout mode used was Fowler 4 rather than SUTR.
The FIRE data is complementary to the BASS NIR DR1 and DR2, spanning ∼3 dex of overlap in X-ray luminosity (see, right panel in Fig. 1), and aside from some lower L X and lower-z object of the DR1 and some higher L X and higher redshift source in the DR2, all the three BASS samples should be studying AGN with similar properties. The combined NIR dataset is representative of the parent BAT AGN sample.
The 65 NIR 0.8-2.5 µm spectra were observed using the FIRE instrument in the high-resolution echelle mode in four visiting runs carried out between April 2018 and April 2019. FIRE is a dual-mode IR spectrometer mounted at the Magellan Baade telescope at Las Campanas Observatory (LCO), Chile. Its primary mode employs a combination of a diffraction grating and four prisms to deliver cross-dispersed spectra covering the whole NIR bandpass in a single exposure, with nominal wavelength resolution of R = λ/∆λ ≈ 6000 for a 0. 6-slit width, i.e., ∆v 50 km s −1 . This slit width was adopted for our program, with the exception of two targets, namely BAT 677 and 1085, which were observed with a 0. 45-slit width (i.e., R ≈ 8000, ∆v 37 km s −1 ), since those were expected to have low M BH . Our observations took place the nights UT 2018 April 5, UT 2018 September 30, UT 2019 March 9 and on UT 2019 April 14-15 during gray time. The observations were executed under clear skies with different airmass conditions (see Tab. 1) and visual seeing variying between 1. 5 to 0. 4, with an average of 0. 75.
For each target, the individual spectra were obtained using the nodding technique in a sequence of ABBA acquisitions, with exposure times ranging between 190 to 623 s (multiples of 10.6 s, since sample-up-the-ramp SUTR readout mode was used for all targets but the ones marked with an asterisk in Tab. 1, which were observed using Fowler 4 readout), depending on target magnitude and observing conditions. The acquisition sequence involved a short arc frame (ThAr) just after the target observation in order to correct for telescope flexure and obtain the wavelength solution. For long (>300 s) exposures, OH airglow was used to improve the wavelength calibration. By default, the wavelengths are calibrated in vacuum. Sky and dome (Qz) flats were acquired to correct for detector illumination and pixel gain variations across the slit, respectively. Data was reduced with the Interactive Data Language (IDL) pipeline FireHose v2 package (Gagné et al. 2015), which performs 2D sky subtraction and extracts an optimally weighted 1D spectrum. Nearby A0V stars were observed during the night in order to derive relative flux calibrations. We corrected the atmospheric absorption features (H 2 O, CO 2 , CH 4 and O 2 ) using the software tool molecfit (Smette et al. 2015). Molecfit uses a radiative transfer code to simulate the atmospheric transmission taking into account local weather parameters (temperature, pressure, humidity, etc.), recorded at the LCO/Magellan site.
Flux-calibrated and redshift-corrected NIR spectra of the 65 AGN are shown in Fig. 2 (as a figure set) smoothed using a Savitzky-Golay filter, which preserves the average resolving power. The spectra are ordered by their BAT ID in Fig. 2. Regions of low telluric transmission are plotted in gray. The locations of some of the most intense NIR emission lines are labeled and indicated with dashed purple lines in Fig. 2. The reduced spectra will be available on the BASS survey website 5 .

SPECTRAL MEASUREMENTS
Below we describe the NIR emission line fitting analysis of the FIRE spectra.

Near-infrared Spectral fitting procedure
Our emission-line fitting approach is similar to the one adopted for the BASS NIR DR1 (e.g., Lamperti et al. 2017) and DR2 (den Brok et al., submitted). We use PySpecKit (v0.1.21), an extensive spectroscopic analysis toolkit for astronomy, which uses a Levenberg-Marquardt algorithm (Ginsburg & Mirocha 2011). We fit the Paζ (0.9−0.96 µm), Paγ (1.04−1.15 µm), Paβ (1.15−1.30 µm), and Paα (1.80−2.00 µm) spectral regions separately, to ease the fitting convergence. Each so-defined independent spectral region contains at maximum three lines coming from permitted species in the BLR. There are only few differences with the emission lines considered in Lamperti et al. (2017, see their Tab. A1), which are the following: we fit also the [N i]λ10404Å in the Paγ region and the Br λ18179.1Å, H 2λ18345Å, He iλ18635Å, and [S xi]λ19196Å in the Paα region, while we did not include in the fit [Fe ii]λ9227Å in the Paζ spectral region, since it is a faint iron emission and was not detected in our FIRE observations.
We first deredden the spectra using the Galactic extinction value E(B − V ) (Schlafly & Finkbeiner 2011) as listed in the IRSA Dust Extinction Service 6 , and redshift correct the spectra. We employ a single first-order power-law fit to model and remove the continuum. For each spectral region, we estimate the continuum level using sections of the wavelength range free of emission lines, i.e., excluding 20Å around the narrow lines and 150Å where a broad component was expected, e.g., in permitted species.
As the focus of our investigation is to study the BLR properties, the main goal is to derive the FWHM and flux of broad line species. Emission lines are modeled using a combination of Gaussian profiles, namely one component for all line species associated with the NLR and an additional Gaussian to account for the BLR component in all permitted transitions of each spectral region. The relative central wavelength of the narrow lines are tied together, but are allowed to shift together by a maximum of 500 km s −1 with respect to the systemic redshift listed in Tab. 1. As an initial input value for the width of the narrow lines, we used the best-fit width of the [S iii]λ9531Å line, which is the strongest (and the narrowest, see, e.g., Fig. 3 in Riffel et al. 2013) narrow emission line in the rest-frame NIR wavelength range, with only a minor blending with the Pa on its red side. The width of the narrow components in each spectral region are then tied together.
The threshold between narrow and broad components is set at FWHM = 1200 km s −1 , as consistent with the division defined in BASS NIR DR1 (Lamperti et al. 2017). The central wavelengths of the broad components are allowed to shift up to ≈1000 km s −1 (Shen et al. 2016). Since the widths of the broad components in each spectral region are likely to be similar, we tie them together to avoid unnecessary complexity and degeneracy 7 . As each spectral region is fitted separately, the resulting best-fitting BLR component can be different in each region. The same applies to the NLR estimate, which might vary in different spectral regions.   Gaussians was run twice as in the normal case, in order to discard the components in case those were below the detection threshold. We then correct the measured FWHMs to account for instrumental resolution, even though it is not a substantial correction for broad lines.
To estimate the uncertainties related to emission line measurements in each spectral region, we repeated the fit ten times, adding each time an amount of noise σ randomly drawn from a normal distribution with the deviation equal to the noise level. We computed the median absolute deviation of the ten measurements and we used this value as an estimate of the uncertainty at the one-sigma confidence level (c.l.). We estimate the flux upper limits (at 3σ) on the (undetected) broad line components by assuming a FWHM=4200 km s −1 , which is the average FWHM of the broad line detections in our FIRE dataset. We then visually inspected all the fits and assigned quality flags, following the classification nomenclature of the first BASS paper . Quality flag 1 refers to spectra that have small residuals and very good fits. Flag 2 means that the fits are not perfect, but still acceptable. Flag 3 is assigned to not completely satisfactory fits for high S/N sources due to the presence of either absorption lines, additional components in the fit or structure in the residuals, making the fit decomposition more uncertain. Flag 4 refers to spectra with low S/N and/or strongly affected by telluric residuals, the best-fit NLR and BLR estimates are highly uncertain. Flag 9 refers to spectra where no emission line is detected. The results of the broad emission line fits are presented in Tables A1-A3 for the Paγ, Paβ and Paα broad lines. Tables A4-A5 report the few cases that needed additional components in their spectral fit. Figure 3 shows examples of emission-line fits for a Sy 1, BAT 744 (top panels), with fit quality flags 2 (2,1) in the Paγ (Paβ, Paα) regions, and for a Sy 1.9, BAT 557 (bottom panels), with fit quality flags 2 (1,4) in the Paγ (Paβ, Paα) regions. The remainder of the best-fit models derived for the full FIRE dataset in the spectral regions Paγ, Paβ and Paα are shown as figure set in Fig.  3 (available on the online journal) ordered by increasing BAT number.

M BH MEASUREMENTS
Below we briefly describe the two independent M BH estimates adopted in order to understand if the NIR view of the BLR gives consistent BH mass estimates compared to the more often adopted Hα-based BH mass estimate, in case of Sy 1 up to Sy 1.9 AGN, and the σ -based BH mass estimates, which include Sy 2 AGN lacking broad Hα line and strong AGN continua but with NIR broad lines. Combining Eqs. 1-2, we evaluate the individual virial factors f of our sample as: where M BH,σ is the BH mass estimated from the optical stellar velocity dispersion (Sect. 4.1), and M vir,line is the virial-based BH mass estimated for each near-infared broad line considered in our study (see Sect. 4.2).

BH mass estimate from the M BH − σ scaling relation
The relation between M BH and the bulge stellar velocity dispersion is probably the most fundamental and most used BH-host scaling relation due to its instrinsic small scatter (∼ 0.3 dex Kormendy & Ho 2013;van den Bosch 2016;Saglia et al. 2016;Shankar et al. 2016;de Nicola et al. 2019;Marsden et al. 2020), and lack of strong redshift evolution, at least until z ∼ 1 (e.g., Shen et al. 2015;Sexton et al. 2019 Saglia et al. 2016;de Nicola et al. 2019). We adopt the M BH − σ relation proposed by Kormendy & Ho (2013), that is calibrated on dynamical-based M BH measurements of local ellipticals and classical-bulges, and for which an average virial factor f has already been estimated (Ho & Kim 2014;Yu et al. 2019). More recent estimates of the average virial factor determined by Batiste et al. (2017) and Yu et al. (2019) are consistent, within their (large) uncertainties, with the values found by Grier et al. (2013) and Ho & Kim (2014).
The M BH −σ relation of elliptical and classical bulges from Kormendy & Ho (2013) is amongst the highest relations in normalization in the M BH vs σ plane up to date (see, e.g., Fig. 2 in Ricci et al. 2017c). This means that at given velocity dispersion it predicts the largest SMBH masses, at least until predicted M BH ≈ 10 9 M . In other words, the f -factors derived using Eq. 3 are the largest that can be possibly predicted with the currently calibrated M BH −σ relations, even considering the most recent updates (Saglia et al. 2016;de Nicola et al. 2019). We will discuss later how this assumption affects the virial factors f and our analysis.

BH mass estimate from the virial method
The virial M BH estimate implicitly assumes that the motion of clouds in the BLR is virialized. Thanks to the R − L relation established by RM campaigns (see, e.g., Bentz et al. 2013), Eq. 1 can be rewritten using easily-accessible quantities, like broad emission line or continuum luminosity. We adopted the mixed virial BH mass estimator put forward by Ricci et al. (2017b), for a similar mixed M BH virial estimator see also Bongiorno et al. (2014); La Franca et al. (2015). The virial BH mass estimator proposed by Ricci et al. (2017b) is based on either optical (e.g., Hβ or Hα) or NIR (Paα, Paβ or He iλ10830Å 9 ) FWHM and on hard X-ray luminosities, either 2-10 or 14-195 keV. Thus, we can compute the virial BH mass M vir,line using the 14-195 keV luminosity L X and any single broad line reliably detected, i.e., He i, Paβ and Paα. The use of the observed BAT luminosity does not affect our analysis since the majority of our sample has N H < 10 23.5 cm −2 , and the observed 14-195 keV luminosity is almost unaffected by X-ray absorbing columns up to at least N H ∼ 10 23.5 − 10 24 cm −2 Koss et al. 2016).
We can also compute the mass with the FWHM of all the reliable detections and with the weighted 10 average, in presence of more than one reliable detection. We call this case N IR 11 in all subsequent analysis. The mixed virial mass is thus M vir,line ∝ F W HM (line) 2 × L 0.5 X , with line being He i, Paβ, Paα, and N IR samples. The statistical uncertainties on the virial BH mass estimate are then the combination of the errors on the broad FWHM and on the X-ray luminosity L 14−195 keV . For simplicity, we assumed a 5% uncertainty on the hard X-ray luminosity, while for the FWHM we used the uncertainties determined from the spectral fit. The aforementioned statistical uncertainty does not take into account the intrinsic spread of virial BH mass estimates, which is of the order of ∼ 0.5 dex (McLure & Jarvis 2002;Ricci et al. 2017b).
As for the majority of the virial BH mass estimators, the relations in Ricci et al. (2017b) were calibrated against a sample of local RM AGN with Hβ measurements, and therefore an average virial factor f was adopted. This f changes depending on whether the RM calibrating virial masses are based on the FWHM or the second moment of the line profile σ line (see, e.g., Onken et al. 2004;Collin et al. 2006) such that The relations calibrated by Ricci et al. (2017b) adopt as a calibration sample the RM virial masses based on the Hβ line dispersion measured from the rms spectra, σ Hβ,rms (Ho & Kim 2014), and in particular the value f σ = 4.31 ± 1.05, derived in Grier et al. (2013) by requiring that RM AGN reproduce the M BH − σ relation found in quiescent galaxies by Woo et al. (2013). Alternatively, Ricci et al. (2017b) also calibrated BH mass estimators considering the bulge morphology, proposing BH mass estimators also for classical-bulges, f σ = 6.3 ± 1.5, and pseudo-bulges, f σ = 3.2 ± 0.7. These f σ values were determined in Ho & Kim (2014) to put RM AGN virial BH masses on the M BH − σ relation of classical-bulges, given by Kormendy & Ho (2013), and the one followed by pseudo-bulges, determined by Ho & Kim (2014). For more details on the BH mass estimators, see Table 4 in Ricci et al. 2017b.
Since the chosen reference BH-σ relation is the one from Kormendy & Ho (2013) (see Sect. 4.1), we adopt the BH mass relation based on the f σ = 6.3, that is 10 The weights are the square of the inverse of the measured broadline FWHM uncertainties, 1/σ 2 i 11 The term N IR in italic should not be confused with the general term NIR that indicate the 0.8-2.4 µm wavelength range. relation b3 in Tab. 4 of Ricci et al. (2017b), in the classical bulge case with a = (7.96 ± 0.02).
In order to more easily compare with literature works about the f -factor (see, e.g., Collin et al. 2006;Decarli et al. 2008;Shen & Ho 2014;Mejía-Restrepo et al. 2018), we convert the individual f -factors from f σ to f F W HM using Eq. 4 and assuming the FWHM/σ line ratio for a single Gaussian case, i.e., √ 8 ln 2 ≈ 2.355. This is equivalent to rescale the f σ adopted as 6.3/(2.355) 2 . We denote this rescaled average virial factor as f 0 . The symbol f 0 is hence the equivalent to writing f F W HM , and it is the average value adopted to convert the virial BH mass M vir,line to the BH mass M (line), i.e., M (line) = f 0 × M vir,line . With this notation, the f derived in Sect. 5.5 using Eq. 3 will be, strictly speaking, a f F W HM . From now on, we will drop the suffix FWHM and simply write f when referring to f F W HM .
When optical broad Hα mass estimates are available from the BASS DR2 (Mejía-Restrepo et al., submitted), we compare them to the NIR-based virial BH mass estimated here (see, Sect. 5.2). These optical BH masses are computed with the prescription proposed by Greene & Ho (2005, see, e.g., their Eq. 6), based on the FWHM and luminosity of the broad Hα, with virial factor of unity (Mejía-Restrepo et al., submitted).

RESULTS
To help the reader, we here present a brief summary of the topics in the results section. In Sect. 5.1 we introduce the sample of Sy 1.8-1.9-2 types for which we reliably detect NIR broad lines. We compare in Sect. 5.2 the NIR and optical, i.e., Hα, broad line widths to understand if the ability of measuring the BLR velocity and radius changes with increasing obscuration/extinction. We then examine how the medium responsible for obscuring the BLR is related to the material absorbing the X-rays in Sect. 5.3. By investigating the effect of X-ray obscuration on Hα and NIR broad line luminosities we demonstrate that the use of Hα (Sect. 5.2.2) and NIR (Sect. 5.4) broad line luminosities in SE BH mass determinations induce a bias in the BH mass measurements when the column density is above a certain threshold. Finally in Sect. 5.5 we compare two independent obscuration-unbiased BH mass measurements, derive the virial factors and examine whether they depend on some additional parameters. . Hα vs NIR BLR velocities for BASS targets, from left to right: He i, Paβ, Paα and N IR broad lines samples, respectively. The N IR sample is the collection of reliably detected NIR broad-line FWHM, being the weighted average FWHM in case there is more than one reliable broad-line detection. The colored circles mark the Sy 1-1.5 (blue), Sy 1.8 (cyan), and Sy 1.9 (orange) types. Only targets with reliable Hα and NIR detection in each line, as defined in Sect. 5.1 and 5.2.1 are shown. The black filled squares are the more robust optical measurements, having < 10% uncertainties in the broad FWHM (Hα). Pearson correlations coefficients r computed for the more robust subsample (black squares), without separating between Sy subclasses, are reported in each panel.

Broad line detection in reddened Seyferts
presented in this work is composed of 33 Sy 2, 18 Sy 1.9 and 1 Sy 1.8, or a total of 52 reddened Sy types. As stated in the introduction, in ∼30% of cases Sy 1.8 to 2 with narrow emission lines in the rest-frame optical spectra have been found to exhibit broad hydrogen and helium recombination lines, i.e., Paα, Paβ and He i. We focus mainly on these three emission lines since all higher order Paschen transitions are strongly blended with emission lines from other elements, and He i, even though blended with Paγ, is among the most intense rest-frame NIR lines (see, e.g., Fig. 9 in Riffel et al. 2006), and therefore it is possible to detect and deblend faint broad components relatively easily. The rest-frameNIR also contains the hydrogen Brackett series, the strongest emission being the Brγ λ21661Å. The Brγ is isolated, located in a region less affected by low atmospheric transmission, and redder than Paα, thus effects of dust extinction are expected to be even lower. However, its intensity is rather low compared to the Paschen and helium lines (Brγ/Paβ ≈ 0.16 in case B with T=10 4 K and density n = 10 6 cm −3 , see e.g., Osterbrock & Ferland 2006), thus we do not investigate its properties in this work. Based on the fact that the He i, Paβ and Paα lines have been found to have similar FWHMs in Sy 1 and in intermediate types (Landt et al. 2008;Ricci et al. 2017b;Lamperti et al. 2017;Onori et al. 2017a, see also Sect. 5.2 and in particular Fig. 4), we can also examine the collection of nearinfrared broad-line measurements, that we call N IR.
When more than one observation was present, we chose the one with the highest S/N evaluated on the continuum, in each spectral region. Combining the DR1, DR2 and FIRE observations the total parent sample is thus composed of 235 obscured BASS AGN, divided into 168 Sy2, 60 Sy 1.9 and 7 Sy 1.8.
Considering only the most reliable measurements, i.e., only the cases where the spectral fit quality is good (quality flag equal to 1 or 2) and the relative uncertainty on the measured broad FWHM is < 10%, the final reliably detected broad lines are 63/235 (27 +4 −3 %) combining the number of detections of at least one of the above transitions. The broad-line detection rate is expected to be consistent in the three BASS samples since the targeted AGN should be similar in their properties (see, Fig. 1) 12 .

Optical and near-infrared view of the BLR in Sy
1-1.9 It is important to understand whether the measured Hα broad line width and luminosity are accurately trac- 12 The FIRE/BASS sample shows a higher detection rate of hidden BLRs compared to the other two BASS works, even though consistent within 2σ given the poissonian uncertainties. This might be caused by the higher resolving power of this work (R∼6000) with respect to the BASS DR1 (most observations have R=800, ∆v = 375 km s −1 , see, Tab. 4 in Lamperti et al. 2017) and with the better observing conditions, resulting in slightly higher S/N, of this FIRE dataset in comparison to the DR2 (average seeing 0. 7 in the BASS/FIRE sample vs 1 in the DR2, den Brok et al., submitted). Columns: (1)-(2) sample size and type, (3) average FWHM of the broad Hα, (4) average broad near-infrared line, (5)-(6) Pearson correlation coefficient and probability. ing the underlying BLR motion and radius, particularly for intermediate Sy 1.8 and 1.9 where the extinction is higher than in normal optical broad-line AGN. The NIR offers a view of the BLR complementary to the optical, as it can help penetrate into the innermost and fastest moving BLR clouds. Together with X-ray ancillary data provided by BASS , NIR spectroscopy can help to understand biases or systematics in BLR measurements, e.g., the BLR velocity and radius, obtained from the Hα (BASS DR2 Mejía-Restrepo et al., submitted). Indeed, a large community effort is being spent to find intermediate mass black holes (IMBHs, M BH < 10 5 M , Greene & Ho 2004, 2007Reines et al. 2013;Moran et al. 2014;Baldassare et al. 2017;Mezcua 2017;Mezcua et al. 2018;Chilingarian et al. 2018;Martínez-Palomera et al. 2020), as this population has a huge impact on several aspects concerning BH seed formation and BH growth at high redshift (see, e.g., Volonteri et al. 2008;Treister et al. 2013;Pacucci et al. 2018;Inayoshi et al. 2020). One of the most commonly used lines in the local Universe to this end is Hα (Reines & Volonteri 2015;Baldassare et al. 2016), since it is about three times more intense than Hβ, and it is in practice the only line available in the rest-frame optical in case of Sy 1.9. Figure 4 reports the measured broad Hα from BASS DR2 (Mejía-Restrepo et al., submitted) compared to the measured broad near-infrared lines from BASS NIR DR1, DR2 and FIRE targets, separated into Sy 1-1.5 (blue circles), Sy 1.8 (cyan circles), and Sy 1.9 (yellow circles) where both optical and NIR BLR components are detected. We recall that we consider the collection of reliably detected broad-line measurements as the N IR line sample (in italic, see footnote in Sect. 4.2), being the weighted average FWHM in case there is more than one near-infrared broad-line detection. The sample of reliable near-infrared, as defined in Sect. 5.1, and Hα broad lines 13 contains 60, 47, 50 and 86 AGN with He i, Paβ, Paα and N IR broad lines, respectively. The sample with robust optical measurements, i.e., with < 10% uncertainty in the measured FWHM Hα, contains 42, 31, 33 and 60 AGN with He i, Paβ, Paα and N IR broad lines, respectively (see black filled squares in Fig.  4). The measurements show some dispersion around 13 From the BASS Hα DR2 database we excluded BAT 557 whose observation does not have nightly flux calibrations, thus resulting in higher uncertainties on the flux measurement (e.g., Koss et al. 2017).

BLR velocity estimates from Hα and near-infrared lines
the one-to-one locus, with the Paschen-line FWHM being on average higher than the Hα one, which might be driven by the less extincted view in the NIR with respect to the optical band, though the averages velocities measured from Hα and all the NIR lines considered are consistent (< 1σ) within the uncertainties of the mean (see Tab. 2). We found that for the whole sample (see 'all' in Tab. 2), the Pearson coefficients are high (>0.66) and statistically significant, as also confirmed by a Student's T-test. This means that the Hα and NIR measurements of the BLR velocities are statistically describing the same velocity field. The intrinsic spreads with respect to the Hα are 0.14, 0.10, 0.15 and 0.13 dex for He i, Paβ, Paα and N IR, respectively. Note that these correlations between optical and NIR broad lines do not depend on our fitting approach, since we are comparing broad lines that come from different spectral regions that are fit independently.
In order to verify whether some particular classes of Sy are strongly influencing the result, we split the sample into Sy 1-1.5 and Sy 1.8-1.9 and compare the more secure optical and NIR detections (i.e., the black symbols in Fig. 4 and the 'more robust' cases in Tab. 2). The results are confirmed: the optical and NIR BLR measurements give consistent estimate of the BLR velocity fields. Additionally, the FWHM measured in Sy 1-1.5 vs Sy 1.8-1.9 come statistically from the same parent population in both NIR and optical. Similar statistical conclusions are derived even adopting the Kolmogorov-Smirnov test and a non-parametric test like the Kendall τ and bootstrapping to consider measured errors on the FWHM (e.g., pymccorrelation python package, Privon et al. 2020). Therefore, we can conclude that once a broad Hα line is securely detected, the Hα width is in agreement within the intrinsic scatter, with the broadline width as measured in the NIR. This scatter might be partially explained by AGN variability since the spectroscopic measurements of the optical and NIR are often obtained from non-simultaneous observations.
We examine in Fig. 5 the difference in the measured BLR velocities as a function of the line-of-sight column density N H measured in the X-rays . In order to verify the dependencies between the y and x variables, we adopt a standard forward regression using the linmix err routine of Kelly (2007), that employs a fully Bayesian approach, and we fit linear log-log relations log(y/y 0 ) = α + β log(x/x 0 ) .
The linmix err routine can include censored data (i.e., upper limits) on the y-variable. In our analysis, all the times upper limits are on the x-variable, being the N H . The N H upper limits are objects with N H < 10 20 cm −2 ,  In all the subsamples inspected, the Bayesian best-fit linear regressions (solid black lines in Fig. 5 with the 1σ c.l. regions in green) are rather flat, indicating no significant effect of N H on the difference in the line widths due to N H . We note that the best-fit slope β found for He i has an opposite sign than what is found for the Paβ and Paα samples, though all slopes are statistically consistent with zero (p-value >0.1). The reason for this might be the ionization structure in the BLR, since the ionization potential of He i is much higher (≈6 times) than the one of Paα, Paβ and Hα. We can conclude that there is a lack of correlation between the NIR-to-optical BLR velocity measurement ratio and the column density N H , at least until N H 10 23 cm −2 . At higher N H , the sample is too small to derive meaningful conclusions. This result nicely complement what can be seen using the larger statistical sample of optical BASS DR2 (Mejía-Restrepo et al., submitted) having both Hα broad line and N H measurements ). In particular, the FWHM(Hα) -N H plane shows a decrease in the optical Hα width in Sy 1.9 at N H > 10 23 cm −2 , while below this column density the FWHM -N H distribution is rather flat (see, e.g., their Fig. 9).
Observationally, the line-of-sight X-ray column density N H has shown a correlation with occultation events in the X-rays, due to gas clumps located in a dust-free region or at the inner edge of the dusty torus (see, e.g., Risaliti et al. 2007Risaliti et al. , 2011Maiolino et al. 2010;Markowitz et al. 2014;Ricci et al. 2016;Zaino et al. 2020). More prevalent eclipsing events are linked to higher covering factors, placing at least some of the X-ray obscuring material in the dust-free BLR (see also Schnorr-Müller et al. 2016). In fact, Gandhi et al. (2015) have argued that narrow iron Kα line-emitting material could reside in between the BLR and the putative dusty torus. In addition, some works have suggested that the accretion disk could partially contribute to the observed column density in Compton-thick AGN, as they have large inclination angles (e.g., Masini et al. 2016;Ramos Almeida & Ricci 2017). Those observations imply the presence of gas inside the sublimation radius. However, the X-ray column density N H is a line-of-sight local measurement, that is very unlikely to strongly obscure the full BLR. Therefore, we calculated the optical-to-NIR broad line flux ratios as a way to measure the extinction experienced by the broad Hα with respect to the broad NIR emission lines. Figure 6 shows only a mild trend, statistically consistent with being flat given the dispersion (p-value >0.05), when comparing the BLR velocity ratio with the broad Hα to near-infrared broad line fluxes.
From Figs. 5 and 6 we can then conclude that neither the X-ray obscuration nor the BLR extinction significantly affect the measurement of the width of the broad Hα line and of the NIR broad lines, at least among the total local hard X-ray selected AGN sample up to columns N H 10 23 cm −2 , once the broad Hα line detection is secure.

Hα line intensity suppression with NH in Sy 1.9
Rather than affecting the width of the Hα broad line, increasing obscuration might simply diminish the entire line intensity, and therefore the flux, of the broad Hα. Such line suppression would preferentially affect the rest-frame optical lines more than the NIR broad lines. Figure 7 investigates whether there is an increasing trend between the optical/NIR BLR extinction and X-ray column obscuration. This correlation is only marginally significant for the Paα sample (pvalue 0.05), while it is not significant for the He i and Paβ (p-value >0.10 for both). We find similar results employing the bootstrapping and point perturbation method with the pymccorrelation package. This observed trend, though only weakly seen for the Paα, is consistent with results from the BASS optical DR2 (Mejía-Restrepo et al., submitted, see their Fig. 8), where the broad Hα to 14-150 keV X-ray luminosity ratio shows a sharp decrease of 1 dex at N H 10 22 cm −2 for Sy 1.9 AGN, while this ratio remains rather constant up to N H ≈ 10 23 cm −2 for Sy 1-1.5 types.
A trend of increasing Hα line suppression as a function of N H implies that M BH measurements based on the Hα broad line luminosity is biased in presence of obscuration. Figure 8 shows the ratio between the virial (i.e., f 0 = 1) Hα-based BH mass, [M vir,Hα = M vir (F W HM (Hα), L Hα )], and the virial NIR+L X -based BH mass for each NIR line, [M vir,line = M vir (F W HM (line), L X )], as a function of the X-ray obscuration (left) and broad lines flux ratios, used as proxy for the BLR extinction (right). The difference between the BH mass estimates increases with increasing obscuration and extinction, more prominently for Sy 1.9 AGN, indicating that the Hα-based BH mass is biased low in presence of obscuration. The trend expected if the NIR-based BH mass was underestimating the true virial M BH at a given Hα-based BH mass would be positive with increasing N H and increasing BLR extinction, while the opposite trend is observed.
By dividing the sample in two bins of obscuring column and running a Student's T-statistic test on the values in the two bins, we find a statistically significant decrement of ∆ = 0.23±0.14 (0.44±0.17, 0.46±0.13) dex with p-value of 4 × 10 −2 (1.7 × 10 −4 , 1.3 × 10 −3 ) in the BH mass difference between Hα and He i (Paβ, Paα) in the two groups at low and high N H . There is a significant bias when the Hα broad line luminosity is used in SE BH mass estimates in presence of obscuration N H > 10 21 cm −2 , as can bee seen in the two shaded Figure 8. Difference between the virial BH mass estimates from the optical and NIR as a function of obscuration (left) or broad line fluxes (right), proxy for the BLR extinction. As previously done in Fig. 6 and 7, the x-axes in the right panels show the direction of increasing BLR extinction (towards the right). Each row reports a NIR line, from top to bottom: He i, Paβ and Paα. Symbols are plotted with the same color scheme as in Fig. 4. The gray solid lines mark the identity relation. The shaded rectangles in the left panels show the average BH mass difference (together with the error on the mean) computed in two bins of NH , low NH in green and high NH in magenta. Red stars show the difference between the BH mass measurements for Sy 1.9 when the Mvir,Hα is estimated using the same prescription as in the NIR, that is to use the hard X-ray luminosity instead of the Hα broad line luminosity. The vertical panels show the histograms of the BH mass difference in the case of completely Hα-based (black dashed histogram) or when the hard X-ray luminosity is adopted instead of the broad Hα line luminosity (solid red histogram). In the latter case, the distribution of the mass difference is more consistent with scatter around zero. bins in the left panels of Fig. 8 that report the average BH mass difference in the N H bin together with the uncertainty on the mean. Indeed the presence of a bias for the Hα-based BH masses can be seen in the dashed black histogram of the BH mass differences, that are skewed to high (negative) values.
When the same prescription is adopted to measure the optical BH masses, i.e., using the hard X-ray luminosity instead of the broad Hα luminosity as a proxy of the BLR radius (see the red stars for Sy 1.9), the difference between the BH masses derived in the optical and NIR is consistent with scatter around zero offset (see red solid histograms in Fig. 8), with no obvious trend versus obscuration or BLR extinction. Therefore, we can conclude that a mixed BH mass estimate based on the hard X-ray luminosity, as put forward by Ricci et al. (2017b, see also Bongiorno et al. 2014La Franca et al. 2015), can overcome the biases due to extinction and obscuration of the optical Hα broad line in Sy 1.9 and in galaxy-dominated AGN, where the UV/optical continuum emission can be diluted by the host starlight. Whereas, as long as a broad line is reliably detected and the hard X-ray luminosity is higher than what is expected from the emission coming from X-ray binaries and hot extended gas in the host galaxy, the mixed BH mass estimator can be used to get an unbiased BH mass measurement. We note that adopting L X instead of the UV/optical luminosity may introduce a weak dependence on M BH and/or L Edd (see Liu et al. 2021;Lusso et al. 2010, albeit there is significant scatter in the correlations). However, the extremely tight (∼0.2 dex) non-linear L X -L U V relation has been extensively used to compute Hubble diagrams for quasars, after accounting for flux-limit related biases and testing for additional systematics (Lusso et al. 2020, and references therein). Further investigation on possible L Edd dependencies will be addressed in a forthcoming paper (F. Ricci et al. in prep).

Dust extinction toward the BLR and gas absorption in the X-rays
The central and right panels of Figs. 6 and 7 report the expected value of the Hα/Paβ and Hα/Paα flux ratio for case B recombination at T = 10 4 K and n = 10 6 cm −3 (Osterbrock & Ferland 2006). We see a wide range of these ratios, with the Paschen lines ruling out case B recombination in most of the X-ray obscured objects, (see also Soifer et al. 2004;Glikman et al. 2006;Riffel et al. 2006;La Franca et al. 2015).
If anyhow we assume that the Paschen lines respect the case B recombination, we can use those broad-line flux ratios to quantify the amount of extinction toward the BLR A V (BLR). Then we could check whether the decline by dust extinction observed in Fig. 7 is consistent with the gas absorption, by assuming the Galactic dust-to-gas ratio.
We compute the dust extinction toward the BLR A V (BLR) following Domínguez et al. (2013, see their Eq. 3): we assumed the reddening curve k(λ) from Calzetti et al. (2000), and that the relationship between the nebular emission line color excess and the Paschento-Hα decrement is given by log (P a/Hα) obs (P a/Hα) case B where E(Hα − P a) is the E(Hα − P aβ) and E(Hα − P aα) excess, the k(λ Hα ) and k(λ P a ) are the Calzetti et al. (2000) reddening curve evaluated at the Hα, Paβ and Paα wavelengths and (Pa/Hα) are the broad-line flux ratio Paβ/Hα and Paα/Hα. We also assume that A V (BLR) = 3.1 × E(B − V ). We then calculate the A V (N H ) from the Galactic dust-to-gas ratio, N H /A V = 2.69 × 10 21 cm −2 (Nowak et al. 2012). The comparison between the two independent A V estimates is presented in Fig. 9. Figure 9 shows a separation with Seyfert subclasses: most of the Sy 1s lie above the 1:1, thus the A V (N H ) is somewhat underestimating the BLR extinction, whereas Sy 1.9s have A V (BLR) ∼ 1 − 5 and their A V (BLR) are mostly below the 1:1, thus the A V derived from dust extinction is less than what is expected from gas absorption. This result is in agreement with what  Fig. 4, with the addition of Sy 2 as red circles. The dashed areas show the average luminosity ratio (together with the error on the mean) computed in two bins of NH , low NH in green and high NH in magenta (see text for more details). The filled squares are the averages in quantiles of log NH for Sy 1-1.5 (black) and Sy 1.8-1.9-2 (dark red), with bars reporting the error on the mean of the luminosity ratio and the bin size in log NH . Figure 11. The FWHM-NH plane for BASS AGN having reliable broad detected NIR lines as reported in the labels. Symbols plotted as in Fig. 10.
has been found in Shimizu et al. (2018, see, their Fig. 3), where the A V (BLR), estimated from the broad Hαto-hard X luminosity ratio, in Sy 1.9 is lower than what is expected from the Galactic dust-to-gas ratio, particularly for A V (BLR) > 3. We note that we adopt a single Galactic dust-to-gas ratio of 2.69 × 10 21 cm −2 , while in the literature it spans a range from 1.79 × 10 21 cm −2 (Predehl & Schmitt 1995) to the value we assume. Furthermore, there is evidence that in AGN the dust-to-gas ratio is not Galactic (Maiolino et al. 2001). If we adopt the average dust-to-gas ratio from Maiolino et al. (2001), e.g., 1.1 × 10 22 cm −2 , then the derived A V (N H ) would be a factor 1.1/0.269 4.1 smaller than what is shown in Fig. 9. Figure 9 suggests that the extinction towards the BLR and line-of-sight X-ray N H are distinct, i.e. they are due to separate obscuration medium. This difference might arise because the BLR extinction is a "global" diffuse measurement while the N H is a more "local" line-ofsight measurement, or if there is gas within the BLR. Alternatively, it might suggest that the dust-to-gas ratio is somehow different from the Galactic value for BASS AGN. Given the large spread in Galactic dust-to-gas ratio and the uncertainty in what is the A V /N H in AGN environments, we cannot distinguish between these two possibilities.
We might then ask whether the same effect of line luminosity suppression is similarly experienced by Paschen and helium lines in the NIR as observed for the broad Hα line. Figure 10 shows the ratio of the broad Hα (left) and NIR emission lines to the BAT 14-195 keV hard X-ray luminosity as a function of the X-ray column density for all BASS targets with reliable NIR broad-line detections and with available N H . There is a decreasing trend of the NIR to X-ray luminosity ratio with increasing X-ray obscuration, as similarly observed for the Hα. Binning the data in quantiles of log N H and splitting the sample into Sy 1-1.5 (black squares) and Sy 1.8-1.9-2 (dark red squares), Sy 1-1.5 AGN exhibit a roughly constant average L(line)/L X across the N H range probed by our sample, while the L(line)/L X of Sy 1.8-.9-2 types seems to decrease as N H grows. This effect seems particularly more prominent at shorter wavelengths. Running a Student's T-statistic test we find a statistically significant decrement of ∆ = 0.54 ± 0.15 (0.54 ± 0.12, 0.46 ± 0.14, 0.27 ± 0.12) dex with p-value of 6.8 × 10 −6 (4.9 × 10 −5 , 2.1 × 10 −3 , 3.3 × 10 −2 ) in the Hα (He i, Paβ, Paα) to hard X-ray luminosity ratio between the two groups divided at log(N H /cm −2 ) = 21. 25 (21.75, 21.45, 21.85), as can be seen in the two shaded rectangles in Fig. 10 that report the average luminosity ratio in the two N H bins together with the uncertainty on the mean. The decrement observed in the Hα to X-ray luminosity can explain the bias observed in Fig. 8. We note that the decrement in the Hα to X-ray luminosity happens at a column density level that is slightly lower than what observed in the other NIR lines, which might be explained considering that dust attenuation diminishes when moving to longer wavelengths. Near-infrared line luminosity, when available, should be preferred to the Hα broad line luminosity as a proxy for the BLR radius when estimating M BH using single-epoch techniques.
Finally, the last question we ask is the following: is NIR really penetrating into the BLR or is the NIR BLR velocity estimate as good as the Hα (e.g., see Figs. 4, 5 and 6) just because both are biased in the same way with increasing column density? Figure 11 can help address this, showing the FWHM-N H plane for all the BASS Sy AGN having reliable NIR broad line detection and available N H . The Hα BLR velocity estimate is also reported for the subsample of the optical BASS DR2 with N H and both NIR and Hα reliable broad-line detections. When binning in quantiles in log N H , the average FWHM measured in the NIR in Sy 1.8-1.9-2 remains quite constant across the N H range, with only a slight decrease in the highest N H bin (log(N H /cm −2 ) 23.5); the error bars are fully consistent with a constant trend versus obscuration. This behaviour is in agreement with the results presented in Onori et al. (2017b), who used a much smaller sample of local hard X-ray BAT selected obscured AGN (N=17). Also among Sy 1-1.5 types the average FWHM is rather constant, given the uncertainties, until log(N H /cm −2 ) 22, while in the last N H bin there is an upturn. We note that a similar behaviour is observed as well in the optical BASS DR2 in the FWHM(Hβ) -N H plane (see, Fig. 9 of Mejía-Restrepo et al., submitted). This upturn behaviour in Sy 1-1.5 might be related just to small sample statistics, since type 1 AGN with high N H are quite rare. However taken at face value, we can speculate that it might be related to a transition in Sy type when the FWHM is smaller than 4000 km s −1 and N H 10 22 cm −2 . At high X-ray columns and high inclinations, the optical photons might in some cases still find their way through the obscuring medium by experiencing multiple scatterings. If this hypothesis were plausible, the fraction of optical polarized light in this particular sample of Sy 1 with high N H should be higher than what is usually found in normal Sy.
From Figs. 10-11 we conclude that 1) the Hα and near-infrared FWHMs are not affected by obscuration, Figure 12. Histograms of the normalized virial factor f /f0, i.e. BH mass deviation MBH,σ /M (line), for each line subsample, as labeled in each panel. Total sample in black, Sy 1-1.5 in blue, Sy 1.8 in cyan, Sy 1.9 in orange and Sy 2 in red. The shaded area defines the IQR, that is the region enclosing the 25th to the 75th percentile of the total distribution. The 50th percentile is overplotted with a solid green vertical line.
at least up to log(N H /cm −2 ) ≈ 23.5, and 2) the nearinfrared line luminosities are not as strongly extincted as the Hα and can be used in SE BH mass measurements until log(N H /cm −2 ) ≈ 21.5 − 22, depending on the specific line chosen.

The virial factor f
We finally verify if the BH mass estimates obtained from two independent methods, i.e., the NIR+L X based BH mass M (line) and the σ -based BH mass M BH,σ , are consistent and, if not, what are the quantities possibly driving the difference.
The resulting sample is composed of 60, 37, 39 and 88 local BASS Sy having optical σ (Koss et al., in prep.; Caglar et al., in prep.) and He i, Paβ, Paα and N IR robustly detected broad lines available, respectively. We were able to gather optical σ measurements inside BASS DR2 for ∼62-68% of the sample with near-infrared broad line detections, as there are only 30, 23, 24 and 41 sources without available optical σ for the He i, Paβ, Paα and N IR sample, respectively. We show in Fig. 12 the histograms of the normalized virial factor f /f 0 which is basically the ratio M BH,σ /M (line). Each panel reports the distribution of f /f 0 separately for each line subsample, color coded according to the Sy classification as shown in the legend. The vertical green line marks the 50th percentile, and the shaded area encloses the 25th to 75th percentiles of the whole sample (black histogram), also known as interquartile range IQR. The histograms are all roughly centered around 0, meaning that the normalized virial factor is of the order of unity, with a distribution ranging from -1 to 1. Some sources show high f /f 0 values, e.g., 1 < log f /f 0 < 2 and −2 < log f /f 0 < −1 , being considered as outliers. These sources have the ratio of BH masses deviating from each other more than one order of magnitude. We discuss the possible reasons why these sources are outliers in Appendix B. In all the subsequent analysis and plots, these most deviating objects are reported as black crosses, and are omitted when performing linear regression fits. Table B6 lists the broad line measurements of the samples and ancillary quantities.
In Fig. 13   As can be seen from the color gradient in Fig. 13, some of the scatter between M BH,σ and M (line) estimates might be driven by the observed FWHM; smaller FWHMs (lighter pink) are located above the 1:1 locus while broader lines (darker red) are below it. This difference in the BH mass estimates, i.e., the virial factor f , is expected in case the FWHM gets broadened with inclination (see, e.g., Collin et al. 2006;Decarli et al. 2008;Shen & Ho 2014;Mejía-Restrepo et al. 2018).
We now investigate whether the scatter between the two different mass estimates depends on additional variables by examining the ratio between the two mass estimates, M BH,σ /M (line), as a function of either redshift (left panels in Fig. 14), X-ray column N H (center panels in Fig. 14) or BH mass estimate based on the virial assumption M (line) (right panels in Fig. 14). In all panels, the solid gray line shows the 1:1 zeropoint and the dashed gray lines are indicating the ±1σ scat- In order to quantify the presence of possible correlations between f /f 0 and other physical quantities, we perform a forward regression Bayesian fit as outlined in Eq. 6, where x is either redshift, column density or the NIR+L X virial-based BH mass estimate M (line) and y/y 0 is f /f 0 . Table 3 reports the best-fit parameters obtained using the IDL routine linmix err, and solid black lines in Fig. 14 show the log-log linear best fit regression with 1σ c.l. In each panel the Pearson correlation coefficient is also reported.
At fixed aperture, with increasing redshift a bigger part of the host (bulge) is sampled in 1D spectra, possibly producing an increase on the measured stellar velocity dispersion and therefore an enhancement of M BH,σ , particularly relevant in late-type systems (Falcón-Barroso et al. 2017). The BAT AGN are typically found in massive spirals with strong bulges (with high concentration index) in between spirals and ellipti-  Columns are: (1) sample of each line; (2-3) the zero point and slopes of the best-fit relations; (4-5) the Pearson correlation coefficients with the related probability; (6) the probability of the slople being different from zero; (7) same test with the null hypothesis of β = −1, carried out only for the f /f0 vs M (line) correlation.
cals (Koss et al. 2011). As such, the aperture corrections go in different directions and it is not straightforward to determine what is the appropriate aperture correction to apply (this issue is further explored by Caglar et al., in prep.). Moreover, line-of-sight velocity dispersion could be broadened due to the disc rotation (Bennert et al. 2011;Har et al. 2012;Kang et al. 2013;Caglar et al. 2020), and this effect should be higher outside the spheroid. In the left panels of Fig. 14 we can see that there is no strong gradient of M BH,σ along the x-axis, while there is some gradient along the y-axis, meaning that at each redshift there is a range of measured M BH,σ . This implies that M BH,σ has a negligible dependence on redshift. As a matter of fact, the log-log relation (solid black line) has always a flat slope in all cases consistent with zero (see Tab. 3) and there is no clear correlation between f /f 0 and redshift. Therefore we can conclude that aperture effects are most probably not important for this sample. As a matter of fact, the optical slit width is typically ≈ 1. 5 in BASS DR2, and at z = 0.1 the sampled region would be of 1.85 kpc/arcsec ×1. 5 2.8 kpc. The average effective radius in SDSS late-type galaxies (Sa, Sb, Sc) is ∼2.7 kpc (Oohama et al. 2009), therefore the spectral extraction would still be roughly within the bulge or spheroid, even at the highest redshift probed in this work.
If obscuration was biasing the NIR+L X virial-based BH masses, the sample should exhibit a gradient with M (line) as a function of the X-ray column, and also we should see a positive correlation between the f /f 0 and the column density. Our sample shows at most the opposite behaviour, a mild anti-correlation between f /f 0 and N H , and no clear gradient in M (line) at increasing N H (center panels in Fig. 14). The derived best-fit relations cannot be statistically distinguished from a relation with a zero slope, thus again there is no statistical evidence for a dependence with X-ray obscuration. This is indeed not surprising since we show in Sect. 5.4 that the FWHM measured from NIR lines is not affected by obscuration until high X-ray column densities log(N H /cm −2 ) 23.5.
Finally, we explore whether there is an anti-correlation between f /f 0 and M (line), which could be inherited by the adopted definition of f , such that f ∝ M (line) −1 (we recall that M (line) = f 0 × M vir,line , see Sect. 4.2). The right panels in Fig. 14 show that indeed there is an anti-correlation, but the best-fit slopes are always flatter than −1. The correlations between f and M (line) are statistically different from the −1 slope expected simply by definition (see Tab. 3, p-value < 7 × 10 −8 ), while the statistical dependence is significantly different from a zero slope only in the He i and N IR cases. Similar conclusions are reached when boostrapping and point perturbation analysis are adopted. We note that this f −M (line) mild anticorrelation might as well be driven by a more fundamental relation with the FWHM, since M (line) ∝ F W HM (line) 2 , and already from Fig. 13 it is evident that there is a dependency of the scatter between the two M BH estimates and FWHM. A thorough investigation of the f − F W HM dependence will be carried out in a separate publication (F. Ricci et al. in prep.).

DISCUSSION AND CONCLUSIONS
In this work we present 0.8-2.5 µm NIR spectroscopic observations of 65 local BASS selected AGN obtained at Magellan/FIRE, splitted into 13 Sy 1-1.5 and 52 Sy 1.8-1.9-2 types. We fit four NIR spectral regions (0.9-0.96 µm, 1.04-1.15 µm, 1.15-1.30 µm, and 1.80-2.00 µm) to study the most characteristic NIR BLR properties, i.e., the BLR velocities and radii, estimated from the most prominent NIR hydrogen (Paα and Paβ) and helium (He i10830Å) transitions. We combine our NIR FIRE sample with the whole BASS NIR database (DR1; Lamperti et al. 2017, and DR2;den Brok et al., submitted), finding NIR broad emission lines in 64/235 Seyferts 1.8-1.9-2. The results of this analysis confirm the possibility of using the NIR band to deeper probe the BLR conditions also in optical narrow-line AGN. The line showing the highest success rate of broad-line detection in reddened Seyferts is He i (43/235), followed by Paα (24/235) and Paβ (20/235), suggesting a possible correlation between the NIR line intensity and the higher ability of isolating faint BLR components reliably in Sy 1.8-1.9-2 types.
We then complement the NIR BLR view with the one obtained from the optical, namely the Hα, from BASS DR2 (Mejía-Restrepo et al., submitted). In this way we constructed the largest NIR sample of hard X-ray selected local Sy having robust ancillary data available from BASS. In terms of L X and z, our sample is representative of the BASS AGN sample, even though the fraction of Seyfert types is not (in the BASS DR1 there is an equal portion of Sy 1 and Sy 2; see, e.g., Koss et al. 2017), since our program was aimed at detecting the BLR in the most elusive Seyferts 1.8-1.9-2 (235/314 objects, ≈75%).
With this obscuration unbiased sample in hand, we investigate the presence of possible systematics in the BLR characterisation taking advantage of the X-ray, Hα and NIR spectral information in Sy 1 up to Sy 1.9 (and Sy 2, without the Hα but with NIR broad line detection). We verify that the FWHM measured from Hα and NIR lines are consistent in Sy 1 up to Sy 1.9. The same re-sults are found even when splitting the sample according to the Sy subclassification. The Hα and NIR FWHM broad line measurements statistically describe the same velocity field in the BLR. Thus, once a broad Hα line is reliably detected, its FWHM is in agreement with the NIR FWHM within some scatter (≈ 0.10 − 0.15 dex, depending on the specific line). We then demonstrate that the Hα and NIR BLR velocity estimates are not significantly affected by neither obscuration nor BLR extinction, as measured by the flux decrement of the broad Hα flux to the NIR broad line flux, at least until log(N H /cm −2 ) ≈ 23, as also consistently found in the optical BASS DR2, where the FWHM(Hα)-N H distribution remains constant up to log(N H /cm −2 ) ≈ 23 (Mejía-Restrepo et al., submitted).
Rather than the Hα FWHM, it is the entire Hα broad line intensity that gets suppressed with increasing obscuration, implying that the optical Hα photons are obscured by a uniform screen placed outside the BLR. The ratio of broad Hα-to-NIR broad lines flux shows a mild decreasing trend with N H for the Paα sample, while the statistical evidence is weaker for He i and Paβ. We quantify the amount of dust extinction towards the BLR and the level of gas absorption measured in the X-rays, finding a clear separation in Seyferts subclasses. This has been shown already by several other works (Burtscher et al. 2016;Schnorr-Müller et al. 2016;Shimizu et al. 2018). The comparison of two different A V estimates suggests that either the material obscuring the BLR is distinct from the one producing the absorption in the X-rays, or that the dust-to-gas ratio in local hard X-ray selected AGN environments is not Galactic (see, e.g., Maiolino et al. 2001). There are some limitations in the approach adopted to estimate the A V (BLR), the most notable being the assumption of case B recombination, which might not hold for Paschen lines (Soifer et al. 2004;Glikman et al. 2006;Riffel et al. 2006;La Franca et al. 2015). Additionally, we stress that the broad-line ratios Paschen-to-Hα, not only depend on dust extinction but also on collisional effects (i.e., by the ionization parameter U and particle density n) taking place at the high densities of the BLR, 10 10−11 cm −3 (see, e.g., Schnorr-Müller et al. 2016). Therefore the observed dispersion in broad line flux measurements is likely due to collisional effects in addition to extinction.
We then explore if the line intensity suppression with increasing X-ray column is similarly experienced by NIR line as found in the Hα. We find a decrement in the broad lines to hard X-ray luminosity ratio of 0.54 ± 0.15 dex for Hα, decreasing at longer wavelengths to 0.27 ± 0.12 dex in case of Paα, occurring at a level of obscuring column density log(N H /cm −2 ) =21.25 for Hα, moving to higher N H levels going to longer wavelengths, up to log(N H /cm −2 ) = 21.85 for Paα. The Hα broad line intensity suppression with increasing N H induce a bias in SE Hα-based BH masses of ≈0.2-0.45 dex, depending on the specific line chosen to compare the Hα with. The NIR line luminosity should be preferred to the Hα line luminosities, when available, to estimate M BH using single-epoch relations (Kim et al. 2010;La Franca et al. 2015). Notwithstanding, in presence of substantial obscuration log(N H /cm −2 ) 22 the NIR line luminosity can be underestimated, particularly in Sy 1.9-2. Thus to overcome these shortcomings, it is preferable to use a more unbiased luminosity as proxy to estimate the BLR radius, and a mixed NIR+L X approach as put forward by Ricci et al. (2017b, see also, Bongiorno et al. 2014La Franca et al. 2015) is a more unbiased estimate of the virial BH mass. The NIR FWHMs seem as well to be much less susceptible to obscuration, at least up to columns log(N H /cm −2 ) ≈ 23.5.
We finally evaluate the consistency between σ -based (from the optical BASS DR2, Koss et al., in prep.; Caglar et al., in prep.) and mixed NIR+L X virial based BH mass estimates, finding that the two quantities agrees within ∼0.4-0.57 dex scatter. The scatter expected in this case is at least the combination of two factors: i) there is an inherent scatter in the M BH − σ relation of the order of 0.3 dex, and ii) the intrinsic scatter of virial-based M BH estimates is 0.4-0.5 dex. Thus, combining the two log-normal contributions, the expected scatter is of the order of 0.5-0.58 dex, which is consistent with what we find. We note however that there are some works suggesting that the M BH − σ intrinsic scatter might be higher if the host morphology is not properly taken into account (e.g., Falcón-Barroso et al. 2017), thus we caution the reader that the M BH based on σ for single BASS targets might be characterised by a systematic uncertainty >0.3 dex (e.g., Gültekin et al. 2009), as there is an inherent difficulty in applying morphology-based corrections to the BAT AGN sample that is composed primarily by massive spirals and lenticulars with strong bulges, with a high fraction of disturbed mergers (Koss et al. 2010(Koss et al. , 2011, adding further complexity to aperture corrections for stellar kinematics. The ratio of two independent BH mass estimates is here used to derive the virial factor f , for the first time for a less biased sample, since we also consider the Sy 1.8-1.9-2 showing BLR components in the NIR. The BHσ scaling relation adopted to derive M BH,σ is the one proposed by Kormendy & Ho (2013) calibrated for elliptical and classical bulges. Still, by normalizing the virial factor f with the ensemble virial factor f 0 we es-sentially minimize the effects due to the choice of the BH-σ scaling relation. As long as the BH-σ relation for pseudo-bulges and/or late-type galaxies differ from the one of elliptical/classical-bulges by only a normalization term, as it occurs for the M BH − σ relation of pseudobulges proposed by Ho & Kim (2014), the different normalization is absorbed in f 0 . If instead the M BH − σ relation of pseudo-bulges or late-type galaxies has a different slope than the one observed for classical-bulges, then this effect is not absorbed in the f 0 normalization.
Examining the distributions of normalized virial factors f /f 0 divided according to the Seyfert subclassification, our data do not show evidence of different M BH between type 1 and type 2 AGN, as instead claimed by Onori et al. (2017b) and Ricci et al. (2017c). These works used a small sample of BAT-selected local Seyferts 1.8-1.9-2 with NIR broad lines and compare their M BH to the well-known and best studied BATselected type 1 AGN with RM-based BH masses. By matching type 1 and obscured Seyferts in hard X-ray luminosity and stellar velocity dispersion they find that Sy 1.8-1.9-2 are undermassive than RM AGN. This BH mass difference is mainly driven by small sample statistics, since i) their small obscured Seyfert sample do not span the same FWHM range of type 1, being limited to FWHM 3600 km s −1 , and ii) the type 1 control sample is biased to higher FWHM, having FWHM 1600 km s −1 , while in this work, thanks to the additional spectra obtained with FIRE and Xshooter, we much better populate the FWHM-L X plane for both AGN populations.
Dependencies of the normalized virial factor f /f 0 with redshift, obscuration and virial-based BH mass estimate have been explored, only finding a mild dependency with M (line). This mild anti-correlation might be actually due to a more fundamental relation with the FWHM. Indeed, using a sample of about 600 local SDSS optical broad line AGN, Shen & Ho (2014) show that the virial factor, estimated as M BH,σ /M (Hβ), is anticorrelated with the FWHM(Hβ), and claim that this effect is a by-product of the line broadening due to inclination of clouds moving in a disk-like BLR (see, e.g., Collin et al. 2006). Similarly, an earlier attempt by Decarli et al. (2008) demonstrated that there is an anticorrelation between the virial factor f , computed as the ratio of the BH mass estimated from the BH-bulge luminosity M BH − L bul relation and the virial Hβ-based BH mass, and the broad Hβ FWHM, and this inclination effect might bias single-epoch virial-based BH masses. At higher redshift, using accretion-disk models to get a BH mass estimate independent of the virial assumption, Mejía-Restrepo et al. (2018) have further demonstrated that this f -FWHM anticorrelation is found not only for the Hβ line but also in Hα, Mg ii and C iv, and can be explained with inclination effects on the measurements of the velocities of clouds moving in a thin disk BLR. We aim to further investigate this possibility with a more detailed analysis on this issue in a separate paper, exploring possible BLR parameters space consistent with our data (F. Ricci et al. in prep.).
Our findings corroborate that the use of a mixed virialbased NIR BH mass estimator gives a less biased view of the BLR properties, and that the use of a single ensemble average virial factor f 0 for the whole AGN population, might bias our understanding about the AGN demographics and evolution, since we showed that f is at least related (though weakly) to M (line) in two broad lines samples (He i and N IR line data sets). Hopefully in the near future, the astronomical community will take advantage of the James Webb Space Telescope (JWST), that will essentially allow to explore the rest-frame NIR properties of AGN and galaxies also at higher redshift. Indeed, the Paβ line at z ∼ 2 − 3 can be observed in the 1 − 5 µm wavelength range with NIRSpec on board JWST. Additionally, being a satellite, the limitations due to ground-based observations will be overcome and more efficient NIR observation will be carried out, opening the path to the construction of large statistical samples of AGN with rest-frame NIR coverage at higher redshifts.

SUMMARY
With the largest and least-biased sample of local hard X-ray selected AGN with restframe NIR spectroscopy, we characterise three key BLR properties directly related to virial-based BH mass measurements, i.e., the BLR velocity, estimated from broad line widths FWHM, the average BLR radius, from broad line and X-ray luminosities, and BLR geometry/dynamics, enclosed in the virial factor f . Our findings are the following: 3. the material obscuring the BLR is distinct from the one responsible of the X-ray absorption, or the dust-to-gas ratio is not Galactic in hard X-ray selected local Seyfert environments; 4. NIR broad lines are suppressed with increasing obscuration, but the decrement is smaller than the Hα one (0.54 ± 0.12, 0.46 ± 0.14 and 0.27 ± 0.12 for L(He i)/L X , L(P aβ)/L X and L(P aα)/L X , respectively) and occurs at slightly higher N H levels, log(N H /cm −2 ) = 21.75, 21.45 and 21.85 for He i, Paβ and Paα, respectively. Even NIR broad line luminosities should not be used in SE M BH estimates when N H 10 22 cm −2 , and a less biased BLR radius proxy should be used, as the hard Xray luminosity L X ; 5. using two obscuration-unbiased BH mass estimates, one based on the σ and the other based on the mixed near-infared+L X virial mass, we show that the two BH mass measurements agree with each other with an intrinsic scatter of 0.4-0.57 dex; 6. we quantify the virial factors f as the ratio of these two independent BH mass measurements and verify that Sy 1 and Sy 1.8-1.9-2 types have the same distribution of virial factors and that our virial factors are not biased with z or N H but show a mild anti-correlation with M (line). This last finding might be driven by a more fundamental anticorrelation with the observed FWHM expected due to inclination effects.

ACKNOWLEDGMENTS
We thank the anonymous referee for the useful comments that improved our manuscript. FR thanks the LCO Instrument & Operations Support Specialists that supported our runs at LCO using Magellan/FIRE, and in particular G. Prieto for sharing the record of local weather data used to run the molecfit correction on our FIRE observations. We acknowledge support from FONDECYT Postdoctorado 3180506 (FR) and 3210157 (ARL); from PRIN MIUR 2017 project "Black Hole winds and the Baryon Life Cycle of Galaxies: the stone-guest at the galaxy evolution supper", contract #2017PH3WAT (     Columns: (1) BAT ID, (2) quality fit flag, (3) noise measured in the continuum in erg/cm 2 /s/Å, (4)-(5)-(6) Paβ flux (in units of 10 −14 erg/cm 2 /s), FWHM and velocity shift, respectively. Values in columns (4) preceded by '<' indicate upper limits. The upper limits on the broad line fluxes have been computed using a FWHM=4200 km s −1 (5). (This table is available in its entirety in a machine-readable form in the online journal. A part is shown as guidance for the reader regarding its content.) Columns: (1) BAT ID, (2) quality fit flag, (3) noise measured in the continuum in erg/cm 2 /s/Å, (4)-(5)-(6) Paα flux (in units of 10 −14 erg/cm 2 /s), FWHM and velocity shift, respectively. Values in columns (4) preceded by '<' indicate upper limits. The upper limits on the broad line fluxes have been computed using a FWHM=4200 km s −1 . (This table is available in its entirety in a machine-readable form in the online journal. A part is shown as guidance for the reader regarding its content.) Columns: (1) BAT ID, (2) quality fit flag (same as in Tab  Columns: (1) BAT ID, (2) quality fit flag (same as in Tab B. F/F 0 OUTLIERS There are some targets showing an outlier nature in the f /f 0 distributions that have been omitted in our analysis. These sources are five in total, namely: IDs 670, 677, 1131, 1465 and 1470. The outliers are mostly (4/5) located above the 1:1 line in Fig. 13, meaning that either the M BH,σ is overestimated, or that the NIR BH mass is underestimated. The former case would indicate that the measured optical stellar velocity dispersion is unreliable for those sources, due to AGN contamination of the stellar absorption features, aperture effects or to host galaxy rotation contamination, while the latter case would happen if there is substantial obscuration along the line of sight, such that the highest velocity clouds of the BLR are not observable, being embedded in some obscuring medium. These four outliers are IDs 670, 1131, 1465 and 1470.
For BAT IDs 670, 1131, 1465 and 1470 the NIR+L X based BH mass is somehow underestimated. Those are part of the NIR DR2 (den Brok et al., submitted): those NIR spectra show only broad He i line, as Paα is located in a heavily affected telluric region and Paβ is only detected, marginally, as a narrow line. All the other higher-order Paschen lines are missing as well. For ID 670, the broad He i transition might not be (fully) tracing the BLR, since a blueshifted outflow in the [O iii] has been detected, with v max = 1072 +80 −990 km s −1 (Rojas et al. 2020).
For IDs 1131, 1465 and 1470 (Sy 1.8, 2 and 2, respectively), such information is not available since they were not part of the sample studied in Rojas et al. (2020). However, the optical BASS DR2 spectrum of BAT 1131 appears to have some ionized outflows in the [O iii], as the residuals from a single Gaussian fit are rather prominent, blueshifted and asymmetric. Therefore the broad component detected in the He i might not be describing the BLR. Additionally, this source, a Sy 1.8, is partic-ularly AGN-dominated in the optical spectrum, making the stellar velocity dispersion measurement a bit tricky.
Finally, we note that for ID 1470 the Galactic extinction is fairly important, E(B − V ) = 0.48 (corresponding to an A J = 0.354 mag according to the IRSA dust database 14 ), while for 1465 is quite relevant, E(B − V ) = 3.03 (A J = 2.150 mag according to the IRSA database), probably affecting the NIR spectrum measurement.
BAT 677 is part of the FIRE sample, and has a good Paβ fit (FWHM=2401±8 km s −1 , fit quality 2). The M BH,σ is computed from a rather small velocity dispersion measurement σ = 41 ± 5 km s −1 performed on the CaT region. We note that the BH-σ relation is not calibrated for σ < 65 km s −1 , thus extrapolating to such low velocity dispersion might introduce a higher error budget on the M BH,σ estimation. We finally note that if instead of σ = 41 we adopt σ = 68 km s −1 , as measured from the Ca ii and Mg i absorptions, the resulting M BH,σ is about an order of magnitude higher (M BH,σ = 10 6.44±0.20 M ) and it is more consistent with the virial NIR+L X based estimate.
Therefore, for the former four sources the NIR+L X based M (line) is not reliable, either due to the possible presence of outflows in ionized material disguised as BLR components (670, also possibly 1465 and 1470) or high Galactic extinction in our line of sight (1465,1470). For BAT 1131 both the σ and the NIR measurements are not trustable, for strong AGN contamination in the optical spectrum and for possible outflows in ionized material. Finally for BAT 677 the M BH,σ is more uncertain being located in the low-velocity dispersion extrapolation of the BH-σ scaling relation. Columns are: (1) BAT ID; (2) optical Seyfert classification; (3-6) near-infrared broad line measurements, from either BASS DR1 (Lamperti et al. 2017), DR2 (den Brok et al., submitted) or from this work, i.e., FIRE spectra; (7) logarithm of the mixed NIR+LX -based BH mass, calculated using the FWHM(N IR); (8) logarithm of the virial factor computed as the ratio of the σ -based BH mass with velocity dispersions from the BASS DR2 (either from Koss et al., in prep., or Caglar et al., in prep.) and M (N IR); (9) column density derived from X-ray spectral fitting . Upper limits on the X-ray columns are denoted with <. (This table is available in its entirety in a machine-readable form in the online journal. A part is shown as guidance for the reader regarding its content.)