BASS XXVIII: Near-infrared Data Release 2, High-Ionization and Broad Lines in Active Galactic Nuclei

We present the BAT AGN Spectroscopic Survey (BASS) Near-infrared Data Release 2 (DR2), a study of 168 nearby ($\bar z$ = 0.04, $z$<0.6) active galactic nuclei (AGN) from the all-sky Swift Burst Array Telescope X-ray survey observed with Very Large Telescope (VLT)/X-shooter in the near-infrared (NIR; 0.8 - 2.4 $\mu$m). We find that 49/109 (45%) Seyfert 2 and 35/58 (60%) Seyfert 1 galaxies observed with VLT/X-shooter show at least one NIR high-ionization coronal line (CL, ionization potential $\chi$>100 eV). Comparing the emission of the [Si vi] $\lambda$1.9640 CL with the X-ray emission for the DR2 AGN, we find a significantly tighter correlation, with a lower scatter (0.37 dex) than for the optical [O iii] $\lambda$5007 line (0.71 dex). We do not find any correlation between CL emission and the X-ray photon index $\Gamma$. We find a clear trend of line blueshifts with increasing ionization potential in several CLs, such as [Si vi] $\lambda$1.9640, [Si x] $\lambda$1.4300, [S viii] $\lambda$0.9915, and [S ix] $\lambda$1.2520, indicating the radial structure of the CL region. Finally, we find a strong underestimation bias in black hole mass measurements of Sy 1.9 using broad H$\alpha$ due to the presence of significant dust obscuration. In contrast, the broad Pa$\alpha$ and Pa$\beta$ emission lines are in agreement with the $M$-$\sigma$ relation. Based on the combined DR1 and DR2 X-shooter sample, the NIR BASS sample now comprises 266 AGN with rest-frame NIR spectroscopic observations, the largest set assembled to date.

Active galactic nuclei (AGN) are accreting, supermassive black holes (SMBHs) located in the center of certain galaxies. They can be among the most luminous, nontransient objects in the known universe (Bañados et al. 2018). While AGN spectra have been extensively analyzed in many wavelength regimes from radio to gamma rays, the rest-frame near-infrared (NIR) wavelength regime (0.8 − 2.4 m) has, to date, only been sparsely studied. Early works include studies of large samples (27 sources, Glikman et al. 2006;47 sources, Riffel et al. 2006; 23 sources, Landt et al. 2008). Over the past few years, studies have increased the number of sources investigated (50 sources, Mason et al. 2015;41 sources, Onori et al. 2017; 102 sources, Lamperti et al. 2017;40 sources, Müller-Sánchez et al. 2018). Spectroscopic NIR observations are advantageous because the NIR wavelengths are less susceptible to interstellar dust extinction by up to a factor 10 as compared to the optical regime (Goodrich et al. 1994;Veilleux et al. 1997;Veilleux 2002), allowing more obscured AGN to be studied (e.g. Lamperti et al. 2017). The NIR band also contains a wealth of emission lines that can help to characterize the ionization structure of the material that may eventually feed the accreting SMBH.
Hydrogen Pa ( = 1.8751 m) and Pa ( = 1.2818 m) are prominent emission lines that are regularly found in the NIR regime. Previous studies have used these lines to derive black hole mass estimates ( BH ) based on their width and strength (Kim et al. 2010;Landt et al. 2013;Kim et al. 2015;La Franca et al. 2015;Ricci et al. 2017e). In certain sources, broad NIR line components have been detected in galaxies that lack broad H or H (e.g., Goodrich et al. 1994;Veilleux et al. 1997;Smith et al. 2014;Lamperti et al. 2017). This is explained by dust obscuration within the host galaxy. Consequently, the Paschen lines provide an additional way to derive black hole masses for obscured AGN (e.g. Ricci et al. 2017e). Furthermore, NIR [Fe ] emission lines can be used to study physical characteristics, as they give important clues on the detailed structure of the emitting gas and they constitute important cooling lines (Riffel et al. 2013;Marinello et al. 2016). In addition, several high-ionization coronal lines (CLs; ionization potential > 100 eV) can be found in the NIR spectral region, such as [Si ] Mazzalay et al. 2010) or the mid-IR region (e.g. [Ne ] 14.3 m; see Sturm et al. 2002). Because of their high-ionization potential (IP), CLs are hard to produce in starburst regions (Marco & Prieto 2005). While type II supernovae can also cause CL emission (Komossa et al. 2009), the lines are generally weak and short-lived (Izotov & Thuan 2009). Since CLs mostly survive only very close to a hard ionization source, they are generally unique tracers of AGN. A proposed mechanism for producing these lines is a strong, central source of intense ionizing continuum in the energetic ultraviolet (EUV) and soft X-ray bands that photoionizes the species (Shields & Oke 1975, Rodríguez-Ardila et al. 2011. Another proposed mechanism is shocks of high-velocity gas clouds that interact with the narrow-line region (NLR) gas (Osterbrock & Parker 1964, Oke & Sargent 1968. These shocks heat the gas to high temperatures ≥ 10 6 K (Oliva 1997). With greater sensitivity of observations, however, emission mechanisms such as shocks could produce detectable CL emission in the absence of AGN, though only rarely in some of the highest star formation mergers in nearby luminous infrared galaxies (Rich et al. 2011). Finally, both mechanisms may occur simultaneously to explain the observed line ratios Geballe et al. 2009;Rodríguez-Ardila et al. 2011). If photoionization is the main generator of these emitting species, a hard radiation field is needed in order to consistently match up the levels of ionizing photons required to produce CL emission (Oliva 1997). This is consistent with the EUV and soft X-rays seen in many AGN, suggesting that CL emission scales with the AGN X-ray emission.
The main interest in NIR CLs mainly derives from the fact that they may be used to detect AGN in dusty environments because of the lowered effect of extinction in the NIR. In the UV to IR regime, the dominant source of obscuration is dust, while high columns of gas are the most important cause of extinction in the X-ray (see the review by Hickox & Alexander 2018). Theoretical arguments indicate that the accretion rate onto SMBHs peaks during the period when the AGN is obscured by dust and gas (e.g. Hopkins et al. 2009). Furthermore, hard X-ray observations show that a large fraction of SMBHs are located in gas-rich (e.g., Koss et al. 2013Koss et al. , 2021, dusty nuclei of galaxies (e.g. , and a large fraction are obscured by high columns of gas (e.g. Brandt & Alexander 2015;Kocevski et al. 2015;Koss et al. 2016;Ricci et al. 2017a). This is further highlighted by the fact that recent NuSTAR observations have found an increasing number of nearby, low-luminosity, Compton-thick AGNs (e.g. Annuar et al. 2015;Ricci et al. 2016;Annuar et al. 2017). Finding and correctly identifying obscured AGN has implications for observational cosmology. As a majority of the AGN population is obscured, a complete census of all sources, obscured and unobscured, is needed to correctly constrain the evolution of SMBH growth over cosmic times. With the advent of the James Webb Space Telescope (JWST), it will be possible to perform infrared spectroscopic observations with an unprecedented sensitivity (Gardner et al. 2006). NIR CLs thus provide several advantages for the identification of AGN activity.
In this work, we investigate CL emission from AGN selected above 10 keV from the Burst Array Telescope (BAT) on the Neil Gehrels Swift Observatory. We examine the properties of CLs in the largest sample of AGN with NIR spectra to date with the goal of learning about the physical mechanisms behind their production. An additional goal is to determine the rate of appearance of such lines in the NIR to determine their variability as a robust tracer of AGN activity. For the distance calculations in this work, we use the concordance cosmological model with Ω = 0.3, Ω Λ = 0.7 and 0 = 70 km s −1 Mpc −1 .

Sample
The BAT AGN Spectroscopic Survey (BASS) project is a collaborative effort to characterize a complete survey of local hard X-ray selected AGN Ricci et al. 2017b), based on the Swift-BAT all-sky survey. This 105 month Swift-BAT all-sky survey has identified 1632 objects, of which 1105 (68%) are AGN (Oh et al. 2018). Due to the hard X-ray (14−195 keV) AGN selection, the sample is nearly unbiased with respect to obscuration up to Compton-thick AGN Koss et al. 2016) and very faint AGN due to X-ray flux limits (e.g Ichikawa et al. 2017). For the second data release(DR2), Very Large Telescope (VLT) X-shooter observations in queue mode were obtained for 269 AGN over several semesters (098.A-0635, 099.A-0403, 0101.A-0765, 0102.A-0433, and 0103.A-0521; these were carried as filler programs), focusing on Type 1.9 or Type 2 AGN or newly identified AGN. A key goal of the high spectral resolution was to measure black hole masses from velocity dispersions in Type 1.9 or Type 2 AGN (Mejía-Restrepo et al. 2022;Koss et al. 2022a), but the NIR arm also provides access to less obscured features. The median seeing was 1. 0 based on the Differential Image Motion Monitor in the band with a standard deviation of 0. 7. A summary and information on the individual observations can be found in Table A.1 in appendix A. From the X-shooter DR2, we only selected nearby AGN ( < 0.5) and excluded beamed AGN (Paliya et al. 2019) to avoid sources with differential beaming of the X-ray emission. Additionally, we have included 10 archival observations of BAT AGN in our sample (e.g. from 086.B-0135) fulfilling the conditions mentioned above (low redshift and no beamed AGN). Our final sample totals 168 unbeamed AGN (Table A.1), of which 110/168 (66 %) are Seyfert 2, 28/168 (17 %) are Seyfert 1.9 and 30/168 (18 %) are Seyfert 1 -1.8 type AGN with broad H . The final sample is biased toward Seyfert 1.9 and Seyfert 2 AGN compared to BAT-detected AGN, which show equal fractions of Type 1 and Type 2 AGN ). Depending on the spectral setup of the instrument, this includes sources with either > 0.3 (if the full NIR range of 9940 -24 790 was covered) or > 0.1 (for which we only have limited NIR coverage of 9940 -21 010). Figure 1 (left) shows the hard X-ray versus redshift plane of the sample of AGN for which NIR spectra were obtained. Figure 1 (right) shows the redshift distribution of our sample.
For completeness, we include all sources from BASS NIR data release 1 (DR1) (Lamperti et al. 2017) in our analysis. This sample consists of 102 NIR spectra of nearby AGN from several observation programs. Most of the sources (55/102) were observed from the 2.2 m NASA Infrared Telescope Facility telescope, with resolution of R = 800 -1000. Seven of the 102 sources were taken with the Florida Multi-object Imaging Near-IR Grism Observational Spectrometer (FLAMINGOS) at the Kitt Peak 4m telescope. Additional sources were taken from archival data from Gemini. The DR1 sample shows a bias toward Seyfert 1 galaxies (∼68% are Seyfert 1−1.9), due to the setup of the archival surveys. We refer the reader to Lamperti et al. (2017)  The sample of additional DR2 data of reduced spectra will be made public on the BASS survey website. We note that the BASS follow-up with X-shooter is ongoing; in this study we use X-shooter observations taken through 2019 October 13. The additional X-shooter observations taken since 2019 October 13 will be presented in later BASS releases and are part of ongoing European Southern Observatory (ESO) programs. The additional data will include other NIR spectroscopy efforts within BASS that are currently ongoing, including follow-up of 65 BASS AGN with Magellan/FIRE (Ricci et al. 2022) and with Palomar/Triplespec (M. Balokovic, in preparation).

Observational Setup
The observations were all carried out with X-shooter, a multiwavelength (0.3 -2.5 m) echelle spectrograph with medium spectral resolution = 4000 − 18 000 (D'Odorico et al. 2006;Vernet et al. 2011). It has three spectroscopic arms, each equipped with optimized optics, dispersive elements, and detectors. Two dichroics are used to split the incoming light into the three arms for efficient observation of all three arms simultaneously. The NIR arm has a wavelength coverage of 1 -2.5 m and includes the traditional atmospheric bands J, H, and K.
The bulk of the observations were carried out between 2017 and 2018. A summary of all observations is listed in Table A.1. Two spectral setups were chosen for the NIR arm: for 84/168 (50%) we obtained full coverage between 0.994 and 2.479 m, while the other 84/168 (50%) had more limited coverage between 0.994 and 2.101 m. The slit width was set to 0. 9, giving a spectral resolution of ∼ 5400 (note that in one archival observation, a slit width of 0. 4 was used; see Table A.1). For most observations the typical total integration time was set to 480 s (for 60/168) or 960 s (for 71/168). To remove the thermal background, sky emission lines, and detector artifacts, the science targets were observed in two positions on the slit in an ABBA nodding sequence with a 5 nod throw. In one archival source a STARE observation was used.
We also obtained an independent estimate of the spectral resolution with the penalized pixel fitting method (ppxf; Cappellari & Emsellem 2004;Cappellari 2017) by fitting stellar absorption lines to individual stars that were observed with the 0.9 slit with a default pipeline extraction of 4 along the slit. To measure the X-shooter spectral resolution, we followed the approach of Gonneau et al. (2020), which was also used for measuring resolutions in the other BASS DR2 optical spectra (Koss et al. 2022a). We use the PHOENIX theoretical spectral library (Husser et al. 2013) as templates, which have much higher resolution ( ∼ 500, 000) than the observations. We fit the 1.45-1.78 m and 2.285-2.38 m regions, respectively, to target stellar absorption features in the CO bandheads in the and bands. In five different stars, we measured =20 ± 1 km s −1 . This corresponds to = 6150, or an FWHM of 0.00026 m at 1.6 m, slightly better than the nominal instrumental resolution listed in the manual.

Data Reduction
The spectra were first reduced using the standard pipeline in the ESO reflex software (Freudling et al. 2013). Pipeline v2. 9.3 was used in all of the sources presented in this paper. We used the default parameters for the creation of the calibration frames. We used the xsh_scired_slit_nod recipe to transform the science and flux-standard frames into flat-fielded, rectified, and wavelength-calibrated 2D-order spectra. The standard 4 extraction region along the slit was used for each spectrum. One of the X-shooter spectrophotometric standard stars was selected for the flux calibration. We corrected atmospheric absorption features that contaminated the spectra (H 2 O, CO 2 , CH 4 and O 2 ) using the software tool molecfit (v1. 5.9;Kausch et al. 2015;Smette et al. 2015). Molecfit uses a radiative transfer code to simulate the atmosphere adopting the observed atmospheric parameters including the ambient temperature, pressure, mirror temperature, and outside humidity.
For the software to work properly the observed spectrum needs to have distinctive, but not saturated telluric features for correction and should avoid intrinsic emission or absorption features from the AGN. With molecfit, no observation time needs to be allocated to telluric standard stars, and because molecfit simulates the atmosphere, small atmospheric changes over a night are better accounted for (Ulmer-Moll et al. 2019).

NIR Emission Line Measurements
For the emission line fitting, the software tool PySpecKit (v0. 1.20;Ginsburg & Mirocha 2011) was used following the procedure of Lamperti et al. (2017). The software is an extensible, spectroscopic toolkit. The fitting procedure relies on the Levenberg-Marquardt algorithm. For the modeled emission lines, a single Gaussian profile is used, or the combination of two Gaussian profiles if the second is detected above 2 above the standard deviation in amplitude above the noise. Before we fit the spectra, we first correct for Galactic extinction, using the built-in deredden function, which takes the − value into consideration (values from Schlegel et al. 1998). The following physical quantities are fitted: the width line and height/amplitude line , as well as the wavelength position obs line of the Gaussian profile. For the final line FWHM measurement, we subtract the instrumental dispersion (i.e., ∼56 km s −1 ) in quadrature.
In order to facilitate the fitting procedure, the NIR spectrum is split into smaller wavelength regions to best fit the varying continuum. The separately fitted regions are (see Table 1 ] spectral region is fitted separately and not included in the Pa region is that the spectra are cut at 1 m, meaning that depending on the redshift, part of the region 9400 -10 000 might be in the NIR arm and part of it in the VIS arm. By separating the [S ] region, issues from the separation of spectra and flux calibrations are minimized. The emission lines we fit in the NIR regime are described in appendix B. The first step in fitting the emission lines is determining the continuum of the spectrum. To allow for more flexibility, especially in telluric-corrected regions with possible residuals, a fourth-order polynomial is fitted to the spectrum. Fitting the AGN continuum using a fourth-order polynomial has been done in several previous studies (e.g Krajnović et al. 2007;Raimundo et al. 2013;Zeimann et al. 2015;Husemann et al. 2020). The continuum level is estimated individually for each of the specific spectral regions (described in Table 1). Emission lines and heavily affected telluric regions are masked.
In certain cases, the continuum shape is irregular. Either the intrinsic continuum or the telluric correction residuals cause an irregular continuum shape and a spline fit is used to estimate the continuum level. In 123 spectral regions for 88/168 (52%) AGN mainly due to strong telluric residuals a spline fit is applied to correct for the continuum. An example where lines are heavily affected by tellurics and the spline fit is applied is presented in Figure D.2 (see Pa region in bottom panel). In Appendix G, we provide more details on the spline fit.
[  The emission lines are fit using Gaussian profiles. We distinguish between narrow lines (FWHM < 1200 km s −1 ) and broad lines (FWHM > 1200 km s −1 ). A broad component is only allowed for the hydrogen recombination lines (Pa , Pa , Pa , and Br ), the strong He lines, and the [S ] 9531 emission line, which is the strongest narrow line in the NIR wavelength range. For [S ] 9531 we also use a third, blueshifted component, which is empirically motivated. Such a blueshifted component is for example also seen in the bright [O ] emission lines (e.g. Rojas et al. 2020). The other NIR lines we fit do not show evidence of significant blue shifted narrow-line components.
As a first step, the Pa region is used to set constraints on the width and offset of the other emission lines. The relative velocity centers of the narrow lines are tied together and the width of the strongest narrow line is used to constrain the width of the other narrow lines in velocity space (with an allowed difference of 200 km s −1 for narrow and 500 km s −1 for broad lines). If no line is found in the Pa range, the Pa range is used instead to constrain parameters. The broad lines are similarly tied together, if detected in the Pa or Pa region. The broad component's centroid wavelength can be shifted by a larger amount. This is empirically motivated by a study of a large sample of AGN looking at shifts of H with respect to the systemic redshift (Shen et al. 2016) finding shifts up to 1000 km s −1 with a mean velocity shift of 109 km s −1 . For high-ionization lines, the allowed offset (of the line's position and width) is set to 400 km s −1 , motivated by observations that these lines tend to be blueshifted.
For a detection, the amplitude line has to be above a certain threshold · , where is the noise level of the surrounding continuum and is the targeted threshold limit. For the determination of the noise level, a window of 0.015 m toward the blue and red of where the line is expected to be, while masking the line itself, is used to calculate the rms value. The threshold is set to = 3 with a width set to the FWHM of other more prominent emission lines. Thus the sample is equivalent width rather than flux limited. For nondetected [Si ] 1.9640 emission lines, we determine upper limits using UL = 3 .
All fits are inspected visually to see whether the lines were fit well. In 14 cases manual intervention is needed for a highionization line because a residual is fitted instead of an actual emission line. In 16 cases, the emission line needs to be fitted manually because of complications with the surrounding noise of the telluric correction.
Errors in the fitted parameters are estimated by performing 20 Monte Carlo simulations drawn from a normal distribution with a standard deviation equal to the noise level in the spectrum. The full table with the measurements is described in Appendix C.

Black Hole Mass Estimation
For narrow-line sources the black hole mass measurements used in this paper are obtained from velocity dispersion measurements using the Ca H + K, Mg , or Ca triplet (around 0.845 -0.87 m) absorption lines using the -relation (Kormendy & Ho 2013). The method is described in detail in Koss et al. (2022b), which is part of this special ApJ series. For broad-line sources, black hole masses are obtained from Balmer lines (mostly H ; see Mejía-Restrepo et al. (2022) for a description of methods). In a future paper, the CO bandheads in the NIR + bands will be used to estimate the mass. For 138/168 (82%) of the sample, black hole mass estimations are available from either Paschen lines or optical velocity dispersions or both. In certain cases, broad NIR emission line components are detected while optical broad Balmer lines are not. For these cases, we use the width and strength of so-called NIR hidden broad lines to estimate the black hole mass and compare the result with the values from other methods described above. Specifically, we use the Pa and Pa -based prescriptions from Kim et al. (2010). We scale down these mass prescriptions by −0.13 dex, to bring them into agreement with the virial factor of = 1 used throughout the BASS/DR2 analyses. Although there is a range of relevant virial factors discussed in the literature, generally in the range ≈ 0.7−1.1 (e.g., Greene & Ho 2005;La Franca et al. 2015;Woo et al. 2015;Yong et al. 2016;Mejía-Restrepo et al. 2018, and references therein), we stress that the differences between them are much smaller than the scatter that dominates the resulting black hole mass estimates in our present analysis (see below). The resulting BH prescriptions are therefore

Ancillary Measurements
In addition to the NIR line measurements, we use X-ray data as well as the [O ] 5007 emission line, which is located in the optical rest-frame regime. The [O ] 5007 observations are from the same X-shooter spectrum; hence instrumental offsets and differences from the NIR in emitting regions are minimized. We note that we did not account for aperture effects in the different slit sizes of the [O ] 5007emission (1.6 in the UVB arm) and CL emission in the NIR (0.9 ). However, as shown by Berney et al. (2015), such aperture effects are negligible even in the more extended [O ] 5007emission, since the emitting region is very concentrated. Likewise for the more compact NIR CL emission, adaptive optics integral-field units (IFU) studies have found the emission to extend to as much as 150 pc (Müller-Sánchez et al. 2011), which would correspond to missing extended emission only in < 0.007 AGN with a 0.9 NIR slit, which represents only 6/168 AGN of our sample suggesting aperture effects are very minimal for the very nearest of our AGN. The [O ] 5007 emission line measurements are presented in a companion paper (Oh et al. 2022). They are detected in the optical data of the X-shooter observations used in our study and have been corrected for Galactic extinction in the same manner. The intrinsic X-ray luminosity and column density H are determined using X-ray observations from Swift-BAT in combination with soft X-ray telescopes such as XMM-Newton, Suzaku, Chandra, and Swift-XRT (see Ricci et al. 2017b for a description of the models). The Swift-BAT telescope provides the observed 14 -195 keV flux.
Throughout the BASS/DR2 analyses, a virial factor of = 1 is used for virial BH estimates that rely on the FWHM of broad emission lines. If one uses the respective line velocity dispersion ( ) instead, this choice would correspond to = 5.5, assuming a Gaussian line profile. The Paschen line prescriptions in Kim et al. (2010) are calibrated against H -based BH estimates from Greene & Ho (2005), which in turn assume = 0.75. Kim et al. (2010)    Additionally, spectral fitting of AGN-specific models to the combined X-ray spectra provides intrinsic luminosities and column density estimates for 116/168 (69%) AGN. X-ray spectral fitting of all 105 month sources will be included in a future release (C. Ricci et al. in preparation). As shown in Koss et al. (2016a) and Ricci et al. (2017b), for BAT observations, the observed flux significantly underestimates the intrinsic flux for H > 10 24 cm −2 which only affects a small number of sources (only 7.6% of the full BASS sample are Compton-thick AGN; Ricci et al. 2015). Because we do not yet have intrinsic flux measurements for the complete sample, we will use the intrinsic 14 -195 keV flux for X-ray luminosity measurements for sources of the 70 month sample (116/168), and use the observed 14 -195 keV flux for the remaining (52/168) sources. In practice the observed BAT 14 -195 keV flux is significantly different (i.e., >20%) for Compton-thick AGN which are rare in the Swift sample (i.e., 7.6%; see Ricci et al. 2017c), which would only correspond to ∼3 sources in our sample of 52 observed 14 -195 keV fluxes. For the derived intrinsic X-ray luminosity the error is < 0.1 dex (Lanz et al. 2019), unless the AGN are Compton-thick, for which the typical errors are 0.4 dex . The typical uncertainty for the observed X-ray luminosity is ∼0.25 dex (Ricci et al. 2017b).
In this study, when talking about the "hard X-ray" flux we are refering to the 14 -195 keV X-ray flux. 4. RESULTS

CLs
If CLs are an efficient tracer of AGN activity they should be detectable in all bright nearby AGN detected in Swift-BAT. Figure 2 (right) shows the distribution of the number of detections per spectrum. In 49/109 (45%) Seyfert 2 spectra, at least one CL is detected. For 3/109 (4%) Seyfert 2 spectra, five or more CLs are detected in a single spectrum. In 19/30 (63%) Seyfert 1 -1.8 and in 16/28 (57%) Seyfert 1.9 spectra at least one CL is detected. The uncertainties of the detection rates are estimated using binomial proportion confidence intervals. The probability confidence interval is set to 1 .
The sample size is smaller than 168 because we exclude here spectra that do not cover the [Si ] 1.9640 emission line due to a combination of the galaxy's redshift and wavelength coverage.  Figure 3 shows the average number of simultaneously detected CLs binned by redshift. The gray area indicates the range in number of detections for each of the different sources in a given bin. A trend can be seen such that we have fewer simultaneous CL detections with increasing redshift, for Seyfert 1-1.8, Seyfert 1.9, and Seyfert 2 galaxies. This decrease is due to a number of factors: less spectral coverage for higher redshifts, the shift of CLs into heavy telluric absorption regions, and generally weaker line fluxes due to increased distance (Rodríguez-Ardila et Lamperti et al. 2017).

Comparison of CLs and X-Ray Luminosity
Naively speaking, CL emission is expected to be driven by soft X-ray and far-UV high-energy photons (>100 eV), which ionize the CL species (Done et al. 2012). So as a first step, we check the correlation between high-ionization and X-ray emission. We use B . The result can be seen in Figure 4. We fit the data using an ordinary least squares (OLS) bisector. In addition, the Python module Linmix is used for regression analysis. The package uses hierarchical Bayesian regression (HBR), which can take the upper flux limits into account. Table 2 lists the regression fit parameters. We note that a positive correlation will be induced due to the correlated axes in a luminosity-luminosity plot. However, in the subsequent analysis, we will mainly focus on and describe the quality of the regression using the scatter around the regression because it is the same in a luminosity-luminosity or a flux-flux plot.
For the relation of [Si ] 1.9640 and X−ray (14-195 keV), the scatter is = 0.37 dex (Figure 4, left). This scatter takes only detections into account. Consequently, the actual intrinsic scatter is most likely higher. We can get a conservative estimate of the lower limit of the scatter assuming the nondetections are 3 detections. This shows that the scatter is actually > 0.42 dex. Using the Python Linmix package we estimate the intrinsic scatter taking nondetections into account. Because the module runs a Markov Chain Monte Carlo (MCMC), we further estimate the uncertainties in the intrinsic scatter. For comparison of [Si ] 1.9640 and X-ray emission, we get intr. = 0.30 ± 0.10 dex. The Pearson correlation coefficient is pear = 0.9 ( pear = 3 × 10 −21 ) showing a strong correlation. As expected when considering flux values rather than luminosities the correlation is more moderate with pear = 0.74 ( pear = 6 × 10 −11 ). In Figure 4 (right), the correlation of [Si ] 1.4300 emission with X-ray emission is shown, again using an OLS bisector and an HBR to fit the data. For the relation of [Si ] 1.4300 and X−ray (14-195 keV), the scatter of the detections is = 0.39 dex. Again assuming the nondetections to be detections, the lower limit of the scatter is estimated to be > 0.4 dex. The intrinsic scatter estimate from the MCMC method is intr. = 0.30 ± 0.20 dex. The Pearson correlation with the hard X-ray luminosity is pear = 0.85 ( pear = 3 × 10 −17 ). For the flux, the correlation is more moderate with pear = 0.61 ( pear = 2 × 10 −6 ). In order to quantitatively investigate whether the correlations of [Si ] 1.9640 and [Si ] 1.4300 with X-ray luminosity differ significantly, the Fisher -test is used, based on the two Pearson correlation coefficients of the luminosity correlation. The two-tailed -value is 0.3, indicating the two distributions are not significantly different.
As a further step concerning the comparison of the correlation of [Si ] 1.9640 and [Si ] 1.4300 with the X-ray emission, we only include sources that show emission from both CLs simultaneously in their spectrum. Figure 5 shows    The scatter of the [Si ] 1.9640-X−ray relation is lower than that of the [O ] 5007-X−ray relation. However, this scatter is a lower limit as it only takes detections into account and the actual intrinsic scatter might be higher. In light of what we find it would seem that CLs are a better proxy for AGN power. In Figure 6 we show the ratio of [O ]     In Figure 8, the emission of [Si ] 1.9640 and the X-ray is shown for both DR1 and DR2. Seyfert 1 galaxies show a statistically higher luminosity. The median luminosity and the 16 th and 84 th percentile ranges of the Seyfert 1 sample are log [Si VI] /erg s −1 Sy1 = 39.8 +0.7 −0.7 and those for the Seyfert 2 galaxies are log [Si VI] /erg s −1 Sy2 = 39.2 +0.6 −0.5 . Furthermore, we find that the scatter is smaller for the [Si ] 1.9640 emission with the hard X-ray for Seyfert 2 galaxies (0.37 dex) than for Seyfert 1-1.9 galaxies (0.45 dex). The scatter across the full sample (Seyfert 1-2 galaxies) is 0.47 dex. Table 3 provides a summary of the scatter and regression parameters for the combined DR1 and DR2 set.
Applying the -test to investigate whether indeed Seyfert 1 and Seyfert 2 luminosity values of the [Si ] 1.9640 emission line differ, we get a -value of = 3 × 10 −5 . We therefore can reject the null hypothesis that the distributions are equal. Figure 9 shows a histogram of the luminosity distributions of Seyfert 1-1.9 and Seyfert 2 galaxies. We see that the [Si ] 1.9640 emission from Seyfert 1-1.9 galaxies is shifted toward brighter values.

Scaling Relations with Black Hole Mass
Theoretical calculations predict a tight dependence between CL emission and the mass of the central black hole, BH (Cann et al. 2018

CL FWHM and Offset
For the FWHM analysis, we take into account how the velocity of the CL-emitting gas clouds depends on the black hole mass (see Netzer 2013 where is the gravitational constant, is the distance from the black hole, Δ line is the velocity measure from the line profile, and ( ) is the virial factor, which takes into account the unknown geometry and orbital structure of the CL region (CLR). Therefore, We Because we look at the ratio, we use for the velocity measure the FWHM determined from our line fitting with PySpeckit. Furthermore, the virial factor cancels out, as we assume similar inclinations between orbits of the CLR. There are seven sources for which this is the case. Figure 11 (left) shows the median of line / [Si ] for the seven sources and errors based on the standard deviation. We see that CLs with higher IP tend to be closer in.
The CL velocity offset is another interesting parameter to analyze, as it can give further information about the kinematics of the CLR. Line offsets are calculated relative to the NLR's velocity offset, as determined from looking at the Pa or He emission line. Figure 11 (right) shows the average mean velocity offset for seven spectra that show all CLs simultaneously. There is a trend toward increasing blueshift with decreasing IP. Figure 12 illustrates that significant shifts are robustly seen, even by eye, in our data, and cannot be explained by poor statistics or a low signal-to-noise ratio (S/N). We focus as an example on the [Si ] 1.4300 and [Si ] 1.9640 CLs because they have the highest detection rates of all the high-ionization lines. For example, Figure 12 shows the spectrum of BAT 1092. It can be clearly seen that the [Si ] 1.9640 line has a systematic blueshift, while such a shift is less clear for [Si ] 1.4300 , which shows a similar offset to the NLR region (as indicated by the dashed line). In Figure 13, we see that such a blueshift is systematically observed for [Si ] 1.9640throughout our sample. In the figure, we color-coded the targets by their respective hydrogen column density. However, we do not find any clear trend with the column density and the magnitude of the velocity offset.

Hidden Broad Lines
Our sample consists of 110 sources that are classified as Seyfert 2 galaxies based on the lack of a broad Balmer line component (FWHM > 1200 km s −1 ) in the optical spectrum. For 59 cases, we have Pa or Pa measurements together with line-of-sight X-ray column density measurements. Figure 14 shows the distribution of the hydrogen column density as a function of the NIR Paschen emission line FWHM. If both the Pa and Pa emission lines have broad components, the average is taken. The vertical dashed line indicates the separation of AGN into Seyfert 1 and Seyfert 2 based on the emission line width and the horizontal dashed line that based on the hydrogen column density (AGN with log( /cm −2 ) >21.9 are considered to be Seyfert 2 galaxies . This is consistent with the fact that the bulk of Seyfert 2 galaxies have a narrow FWHM (< 1200 km s −1 ) and a high column density log( H /cm −2 ) = 21.1 − 25.4 (mean column density: log( H /cm −2 ) = 23.4 ± 0.9). In the following, we only take into account those sources that have velocity dispersion BH   we are interested to see whether the Paschen lines can be used for mass estimations. In Table 4 the Seyfert 1.9 and 2 sources that match these criteria are summarized. In 59 Seyfert 2 sources a Paschen line is found and in 6/59 (10%) a broad Paschen component is detected, but there is no detection of broad H components. These are so-called hidden broad lines. Theses sources show a column density in the range of log( H /cm −2 ) = 22.3 − 23.8 (mean: log( H /cm −2 ) = 23.1. As an example, the spectrum of LEDA 157443 (BAT ID 597) is shown in Figure F Additionally for the 15 Seyfert 1.9 sources considered, 7/15 (47%) show a broad Paschen component despite a column density above log( H /cm −2 ) = 21.9 . These cases have column densities in the range of log( H /cm −2 ) = 22.0 − 23.0 (mean: log( H /cm −2 ) = 22.5). Figure 15 (left) compares the black hole mass estimates using the FWHMs of Paschen and Balmer lines and the stellar velocity dispersion. Figure 15 (right) compares the mass estimates when using the H and the Pa emission line, demonstrating an offset from the 1:1 locus. In Figure 16, the ratio between the FWHMs of H and Pa is compared with the column density. The ratio FWHM(Paschen)/FWHM(H ) might increase with column density. However, this is based on the few Sy 1.9 objects having both H and Paschen broad lines, and more data are needed to definitely understand whether this ratio changes with H ((see, e.g., Ricci et al. 2022)).

Detection of Coronal Lines
A necessary condition for a line to be an efficient tracer of AGN activity is that it should be detected in a large number of targets. We see a trend that with increasing IP, the fraction of detected lines decreases (see Figure 2 left Lamperti et al. (2017) found a higher detection rate of CLs for Seyfert 1 -1.9 galaxies than for Seyfert 2 galaxies (they find a rate of 53% for Seyfert 1-1.9 and 20% for Seyfert 2 galaxies with at least one CL). Despite the large bias toward Seyfert 2 galaxies in the DR2 sample, we detect a larger absolute number of CLs compared to Lamperti et al. (2017) because we have a B .
larger sample, and our VLT observations are of higher quality in terms of spectral resolution and sensitivity, which is essential for deblending the lines. While [Si ] 1.4300 has a fairly high detection rate (30-40%), [Si ] 1.9320 is not detected in our sample. A possible explanation for the nondetection could be a loss of spatial resolution, because of the increased distance to the sources (Rodríguez-Ardila et al. 2011) or a generally poor resolution of the instrument. CLs are thought to be produced in the nuclear region, between the broad-line region (BLR) and the NLR (Rodríguez-Ardila et al. 2006) and so they lose contrast as more nearby continuum stellar light from the host galaxy is included in the spectral aperture (Mazzalay & Rodríguez-Ardila 2007,Mazzalay et al. 2013). Additionally, the nondetection of [Si ] 1.9320 can be explained by the difference in IP and critical density.
[Si ] 1.9320 has a higher IP than [Si ] 1.9640 and is thus likely produced closer to the black hole (as also inferred from our discussion of CLR constraints; see section 5.2). Rodríguez-Ardila et al. (2011) suggest a density gradient toward the center of the AGN and the critical density of [Si ] 1.9320 is lower than the local density of where the emission is produced. As a consequence, the emission of [Si ] 1.9320 might be suppressed due to collisional deexcitation.
Concerning the strength of the [Si ] 1.9640 emission line, looking at the distributions of observed luminosities (see Figure 9), we see that most [Si ] 1.9640 lines are located in the luminosity range of X−ray = 10 39 − 10 41 erg s −1 for Seyfert 1 galaxies and X−ray = 10 38.5 − 10 40 erg s −1 for Seyfert 2 galaxies. Based on the Lamperti et al. (2017) NIR data in DR1, which consists mainly of Seyfert 1 galaxies, and the DR2 sample here, we find that the average flux of Seyfert 1 CLs is higher than that for Seyfert 2 galaxies, indicating that torus obscuration might play a role. We do not find that Seyfert 1 sources show a higher CL detection rate than Seyfert 2 sources. Of the NIR high-ionization lines, the To further expand on torus obscuration, we investigate in Figure 17 whether there is a correlation between the dusty torus covering factor and the high-ionization lines vs. the X-ray (14-195 keV) flux. Ichikawa et al. (2019) measured the dusty (IR) covering factor in several sources in our sample. They computed the geometrical covering factor by assuming the two-phase torus model (described in Stalevski et al. 2016). Ricci et al. (2017d) found a trend of the covering factor with the Eddington ratio. Upon including the upper limit in our regression analysis with the Linmix module, we find evidence of a positive correlation between the covering factor and the high-ionization line flux. Including the upper limits in the regression analysis is important because they place strong constraints on emission at lower covering factors. This trend coincides with the assumption that CL emission originates from a layer close to the torus. In X-ray heated wind, the CL emission becomes more efficient (Pier & Voit 1995), driving the correlation.
For tracing potentially obscured AGN, [Si ] 1.9640 is the most promising CL in terms of detection rate and line detection. Lamperti et al. (2017) Cann et al. (2018) propose might be more prominent in the case of IMBHs.

Comparison of CL and X-Ray Luminosities
Looking at the scatter between the line to that of the IR [Si ] line. Their study found a linear correlation between X-ray emission and CL emission. They found a tighter correlation for Seyfert 1 galaxies and claimed that the scatter is mainly introduced by Seyfert 2 galaxies. Lamperti et al. (2017), however, found that there is no tighter correlation between [Si ] 1.9640 emission and CL hard X-ray flux when looking at Seyfert 1 galaxies compared to Seyfert 2 galaxies. We find, using the full DR1+DR2 sample, that the scatter for Seyfert 2 is tighter ( Sy1 = 0.45 dex and Sy2 = 0.37 dex). Therefore the scatter in the CL emission compared to X-ray emission is not primarily caused by obscuration. This is also evident from Figure 6, as we see that the scatter is larger for the [O ]/[Si ] line ratio, which traces the IP, than for the [Fe ] line ratio, which traces obscuration (Riffel et al. 2006). We note that aperture effects, while playing a role, do not fully explain the scatter in the figure, as even after excluding the most redshifted sources (sources with > 0.1 or > 0.3 depending on the spectral coverage of the X-shooter setup used), the scatter is around 1 dex. Furthermore, also metallicity cannot be the primary cause of the large scatter of the -axis, as the majority of BAT AGN host galaxies have a stellar mass > 10 10 ) and have a constant metallicity gradient.
Another factor is the physical state of the gas in the emitting media. For example, the electron gas density ( ) influences the strength of the CL emission. Rodríguez-Ardila et al. (2011) estimated the CL-emitting region to have a density straddling typical values for the NLR and BLR (10 8 − 10 9 cm −3 ). Using detailed IFU and spectrograph studies of Seyfert 2 galaxies, Rodríguez-Ardila et al. (2017b,a) also found that high values ( > 10 5 cm −3 ) are very likely required. However, Landt et al. (2015) contradicted this, finding the CL gas density is low with ≈ 10 3 cm −3 . To fully understand the influence of the CL gas density on high-ionization emission, more detailed studies of the gas conditions are necessary.
A further potential explanation for the scatter is AGN variability. The NIR and X-ray observations are not made contemporaneously, leading to increased scatter. Looking at Figure 7, where the luminosity of [Si ] 1.9640 is compared with the [Si ] 1.4300 luminosity (thus avoiding scatter caused by differences in observation time), the scatter is not found to be significantly smaller than when looking at the comparison of [Si ] 1.9640 or [Si ] 1.4300 emission with the X-ray emission. So AGN variability is unlikely to account for the scatter. For the comparison of the two CLs, the scatter is = 0.40, while for the comparison of [Si ] 1.9640 with X−ray the scatter is = 0.36 and for [Si ] 1.4300 with X−ray the scatter is = 0.4 Another aspect is radius-dependent variability based on varying distances of the emission regions to the ionizing source. If CLs indeed originate from a region between the BLR and the NLR, then we expect them to be more correlated with the X-ray emission than, for example, with[O ] 5007, which is simply due to the light traveling time; regions that are located further away than the typical X-ray variability timescales will show less variability.
Since the detection frequency and the flux of CLs in Seyfert 2 galaxies are (on average) lower than those in Seyfert 1 galaxies (see Figure 9), obscuration could play a role in CL detection. This strengthens the argument that the CLs are produced in the region between the BLR and NLR, and are thus affected by obscuration by the torus. Thus the CLRs seem to be an extended region in accordance with previous studies (e.g. Rodríguez-Ardila et al. 2006, Landt et al. 2015. In fact, the [O ] 5007 emission shows indeed a larger scatter and lower correlation with the X-ray luminosity, when compared with the [Si ] 1.9640 emission, though the difference is not significant. To further investigate the correlation of [Si ] 1.9640 with X-ray emission, we consider the Eddington ratio ( bol / Edd ) dependence (see Figure 18 top left). Previous studies by Oh et al. (2017Oh et al. ( , 2019 found a correlation between the Eddington ratio and narrow-line emission. The correlation is most likely caused by X-ray heating processes or removal of material by an energetic outflow. B . Looking at Figure 18, there is no clear dependence visible for the correlation with the Eddington ratio ( bol / Edd ), at least for the high-ionization CL considered here.
In Figure 18 (right), we compare the [Si ] 1.9640 emission with the power-law photon index Γ. X-ray spectra can be described to first order by a power law, parameterized by the photon index Γ (for a more detailed description see Ricci et al. 2017b). In a previous study, Rodríguez-Ardila et al. (2011) claimed to have found a linear correlation between Γ and CL emission, while we find none. They postulated that CL emission is predominantly present in sources with a soft excess (when Γ ≥ 2.5). We do not see any correlation when looking at the data. Furthermore, we observe hardly any sources with Γ ≥ 2.5. This could result from bias in our sample, which are predominantly Seyfert 2. Another reason why we do not find any correlation could be that we have a larger sample than they had (we have = 36, Rodríguez-Ardila et al. (2011) had = 13), we cover a larger energy-range (0.5−150 keV against 0.1−2.4 keV) and the photon index is estimated using a more sophisticated model that includes higher-energy photons. It is also possible the soft excess at low energies may be important in the correlation. The column density is also shown in Figure 18 (right). We also do not observe any column density dependence, again indicating that obscuration seems not to cause the large scatter, nor any obvious bias.
Another influencing factor in CL emission could be the central mass of the black hole. Comparing the mass of the black hole ( BH ) with CL emission, only a weak correlation ( pear = 0.48 for [Si ] 1.9640) is found (see Figure E.1). This is also seen in Figure 18, where higher-mass sources are located more frequently in the range of higher [Si ] 1.9640 luminosity values. The mass range covered by our measurements reaches log BH / = 6.5 − 9. A reason for the weak correlation is also the fact that only a small luminosity range is covered. Based on our data, the luminosity of CLs is not a good indicator of black hole mass. However, CLs might be a good tracer of IMBHs according to Cann et al. (2018), who found a correlation between BH and CL ratios. According to their theoretical models, for masses log BH / < 6, the ratio between the fluxes of [Si ] 1.9640 and [Si ] 1.4300 changes by over seven orders of magnitude. Unfortunately however, our data points do not cover this mass range and go only to log BH / > 6.5. The drop of over seven orders of magnitude is explained by the interplay of black hole mass, ionization parameters, and physical properties of the gas, in which, at low masses, the effective number of ionizing photons is a strong function of black hole mass for a fixed Eddington ratio based on standard disk theory. This drop, however, also has implications for the search for IMBHs using, for example, the [Si ] 1.9640 emission line. If the calculations are correct, the ratio peaks precisely in the range log BH / = 6 − 8. This would mean that the high detection fraction of [Si ] 1.9640 may potentially be a selection effect. However, the drop for sources with log BH / > 8 is not seen in our data. The fact that we do not see the predicted drop of the CL emission ratio at log BH / > 8 suggests that the sources have a strong UV emission Figure 19. Outflow movement of gas clumps. Clumps closer to the bicone axis are more highly ionized. In a Seyfert 2 configuration, the projected velocity (i.e. the velocity observed by us) is indeed more blue shifted for species further away from the axis, because the velocity component toward the observer is larger.
even at high masses (that are capable of ionizing the species). It must be noted that the theoretical calculations are based on fixed parameters, such as bol / Edd = 0.1. The BASS sample covers a broad range of bol / Edd at every BH (see also Koss et al. 2017), providing a broader range in parameter space than the Cann et al. (2018) models. We also note that in the BAT sample, high-mass sources have lower Eddington ratios (Ricci et al. 2017d). Such lower Eddington ratios may change the UV ionizing spectrum (Lusso et al. 2010) altering the relation used in Cann et al. (2018). We also investigate the connection between the S/N and the scatter of the CL emission with X-ray luminosity. We separate the data into high and low S/N and compare the scatter. We find that the scatter stays constant irrespective of the S/N cut applied and we see that the points follow the same distribution according to the 2D Kolmogorov-Smirnov test described in Fasano & Franceschini (1987).
So to answer the question about CL correlation with X-ray emission, we indeed see the trend that with increasing X-ray emission, the emission of the CLs [Si ] 1.9640 and [Si ] 1.4300 increases. The scatter is smaller, but comparable with that of [O ] 5007. Will this suffice for CLs to be considered efficient tracers of AGN activity? This is not evident from our analysis. What we have shown is that high-ionization lines are detectable, even with the challenges of telluric absorption. Furthermore, they show a relation to other properties of the accreting system, such as the mass BH or the X-ray luminosity. However, the scatter is still quite large ( ∼ 0.4 dex). An additional advantage is that in obscured AGN (Seyfert 2 galaxies) the CLs are also detected. To fully answer the question about the efficiency of production, we need a wider range of luminosities and Edington ratios, including, e.g., galaxies that are likely to host IMBHs (i.e. not the typical BASS sources).

Constraining the Geometry of the CLR
The link between the FWHM and the IP of CLs has already been intensively studied in previous works. Giannuzzo et al. (1995) postulated that the CLR might occupy different regions in different galaxies based on the wide range of CL FWHMs. Later studies (Reunanen et al. 2002, Rodríguez-Ardila et al. 2011 have found a correlation between the FWHM and IP up to some IP for certain cases. Up to energies of 200 -300 eV, we also observe an increase of the FWHM with increasing IP. If we take further high-ionization lines into account, we see a drop in the FWHM again. However, because the highest ionizing species (>400 eV) are relatively weak, it is difficult to make conclusions. Rodríguez-Ardila et al. (2011) attribute the finding that the increase of the FWHM with IP is only seen up to 300 eV to the combined effect of the electron density gradient toward the center and a spatial extension of the emitting material. The critical density can give insights into possibly why the highest-ionizing species are not observable. One of the common lines detected, [Si ] 1.9640, has a critical density of = 2.5 × 10 9 cm −3 . This sets an upper limit, because certain higher-ionization lines, such as [S ] 19196, have an order of magnitude lower critical density ( = 3.2 × 10 8 cm −3 ). If they are produced closer to the center and there is an increase in density toward the center, [S ] 19196 emission might be suppressed due to collisional deexcitation. The lower critical densities for lines with IPs higher B .
than 350eV are a possible explanation for why the detection frequency drops sharply at high IP (see Figure 2). The FWHM measurements tell us approximately how far from the central ionizing source a certain line is produced within the CLR for a given black hole. Different ionization species will be dominant in different regions, with a mild dependence on BH . Murayama & Taniguchi (1998) proposed that high-density clumps that are radially moving outward produce the CLs. The high density clumps are separated into various segments of different dominating ionizing species. This model cannot explain the difference in FWHM with increasing IP, as the FWHM should be similar, because the high-density clump moves as a whole. The model by Fischer et al. (2017) can already better explain the finding, because it allows for different velocity dispersions in the infalling dust spirals.
Another finding is that the offset of the line peak, i.e. the bulk motion of the emitting material with respect to the observer, is systematically blueshifted for the [Si ] 1.9640 emission line (see Figure 13), while for higher-IP species, the offset is actually redshifted with respect to [Si ] 1.9640 (i.e., they always have a smaller bulk offset). This can be understood in terms of the general geometry of the CLR: Incoming or accreting gas is mainly ionized as it enters the bicone axis (according to Fischer et al. 2017). Due to the hard radiation field or outflows, the material can be accelerated outward. If we observe an AGN in a Type 2 orientation, and the more highly ionized gas is closer to the bicone axis, the higher-ionization lines could be observed at lower apparent radial velocities (the flow being more directed along the plane of the sky). This could explain why all [Si ] 1.9640 emission lines are blueshifted with respect to the NLR, while [Si ] 1.4300 emission is both red-and blueshifted (with respect to the NLR). Outflows have also been found in previous IFU studies looking at the NLR and CLR (Müller-Sánchez et al. 2011). The simplified concept is illustrated in Figure 19, similar to the illustration in Murayama & Taniguchi (1998). The central AGN is shielded by a dusty torus (Marinucci et al. 2016;Ramos Almeida & Ricci 2017) from the observer in the case of galaxies with large covering factors, which is generally the case in Seyfert 2 galaxies. The clouds move outward along different ionization cones, which leads to different observed (relative) velocity shifts/offsets of the CL species. The highest-ionization lines move along a narrower cone closer to the bicone axis, while CLs with a lower ionization move in a wider cone. In a Type 2 configuration, clumps moving in wider cones have a larger velocity component along the line of sight. This could explain why they are more blueshifted (e.g., [Si ] 1.9640 as opposed to [Si ] 1.4300; see Figure 12). Especially in the case of a single cone (or because the second cone is more heavily obscured), if the axis points slightly away from the observer, it might explain why some of the highest-ionization lines are even slightly redshifted as compared to the narrow-line emission. A further indication that this phenomenon is an orientation effect can be seen in Figure 13. A trend can be seen where the sources with lower column density show that both lines are blueshifted, meaning a more face-on look into the center of the bicone. Also, there is a component of gas close to the AGN, essentially at the apex of the bicone. This usually has higher ionization (see Fischer et al. 2017) than the rest of the NLR and does not seem to fit the flow pattern of the more extended (in situ) gas. This could explain the large scatter in offset of the [Si ] 1.4300 emission line.
So from the offset measurements of the CLs, we get constraints of the overall geometry of the CLR. From the FWHM measurements, we have indications that the emission comes from different regions within the CLR and from the offset analysis, we have indications that most of the ionization takes place along the bicone axis, and as a result of orientation effects, this causes the more highly ionized CL species to show a different offset than the other emission lines.
The structure of the CLR has been addressed in past IFU studies (Müller-Sánchez et al. 2011;Mazzalay et al. 2013;Rodríguez-Ardila et al. 2017a;May et al. 2018May et al. , 2020Rodríguez-Ardila & Fonseca-Faria 2020) and extended CL emission with > 400 eV could be observed. However, resolving the innermost parsec region is not yet possible with current instrumentation.
We would like to emphasize that the analysis presented has some limitations. Outflows in Seyfert galaxies are likely complex and more complicated than the conical outflow depicted in the simple sketch in Figure 19. Besides linear outflow kinematics, rotational kinematics are possible for the CLR and NLR (Müller-Sánchez et al. 2011). In addition, a previous study focusing on the NLR has found that some ionized outflows are hollow (Fischer et al. 2013). In the case of a hollow structure, the geometry of the different ionization cones would be more similar, and the picture presented in Figure 19 would not sufficiently explain our findings. However, it is not clear that the CLR follows the same geometry as the NLR, as the CLR is generally situated closer to the center of the AGN (Oliva 1997;Mazzalay et al. 2010).

Hidden Broad Lines
We detect broad emission lines in the NIR in a handful of sources that are optically classified as Seyfert 2 galaxies. These sources presumably consist of AGN where the line of sight is impacted by a moderate column density, and hence by extinction, such that the BLR is completely obscured yet they have a column density above log( H /cm −2 ) = 21.9, to place them in context. Previous studies (e.g. Garcet et al. 2007;Oh et al. 2015;Kamraj et al. 2019) have found sources that have high column density yet show optical broad lines.
As can be seen in Figure 14, for 6/59 (10%) of Seyfert 2 galaxies, broad Pa or Pa components are detected. This is consistent with the 9% fraction found by Lamperti et al. (2017). Furthermore, if we include Seyfert 1.9 galaxies, we detect broad components in 12/75 (16%) sources. This is lower than the 31% fraction found by Lamperti et al. (2017) and the 32% fraction found by Onori et al. (2017). The reason why we have such a low fraction is most likely the low number statistics, as only 12/75 (15%) in our sample are classified as Seyfert 1.8 or 1.9 galaxies. Are these sources challenging the unified model? Lamperti et al. (2017) found that the Seyfert 2 cases with broad NIR components occupy the bottom 11 th percentile of column densities ( log( H /cm −2 ) = 22.4). Focusing on Seyfert 1.9 and Seyfert 2 galaxies, we find a much broader range, extending up to log( H /cm −2 ) = 23.8 (median column density: log( H /cm −2 ) = 23.3). We find that for at least 10 % of Seyfert 2 galaxies, one can detect broad components, so-called hidden broad lines, in the NIR, which then can be used to estimate the mass of the central black hole. The reason why broad lines are detected in the NIR and not in the optical is mainly the decreased obscuration at longer wavelengths. Lamperti et al. (2017) found that sources with hidden broad lines are often merger systems, so the optical broad emission component is most likely obscured by the host galaxy's dust rather than by the nuclear torus. The [O ] to X-ray luminosity ratio is also found to be lower in merging BAT AGN systems (e.g., Koss et al. 2010 and most of the late-stage, close nuclear (<3 kpc) mergers are found in optical Seyfert 2 systems ) rather than in broad-line AGN consistent with this claim. Higher X-ray obscuration is also found to correlate with later merger stages (Koss et al. 2016a,b;Ricci et al. 2017a) and has been predicted by theoretical studies (e. g., Hopkins et al. 2006;Blecha et al. 2018).
If true, this indicates that hidden broad lines are not a confutation of the unification model, because the obscuration of the broad components is not due to the torus, but rather due to more extended host galaxy dust and gas. Indeed, one of the AGN counterparts with hidden broad lines we detect, 2MASX J042340. 80+04080.17, shows a spiral companion indicating a possible merger event (Gonçalves et al. 1999). The second example, ESO 383-18, shows dust winds and Compton-thin dust lanes, which could cause the optical broad lines to be obscured (Ricci et al. 2010). NGC 4941 is a Seyfert 2 galaxy and is marked as a galaxy without hidden broad lines in Yu & Hwang (2005), even though we detect a broad Pa emission component. This galaxy shows no signs of large-scale interactions.
Sources with hidden broad lines are also interesting to investigate in terms of potential differences between the optical and NIR broad-line properties, and how they affect the estimation of the mass of the black hole. In the case of Seyfert 1.9 galaxies, broad H can be attenuated by dust leading to a low black hole mass. Figure 15 presents the value of the mass estimation from velocity dispersion measurements and the broad Paschen line method. No structural offset can be seen and points are spread out equally on both sides of the 1:1 relation line and hence the Paschen lines appear to yield reliable estimates. We compare the mass estimates from the various hydrogen recombination lines H , Pa and Pa in Table 5 and Figure 15. The mass estimation from H emission is lower by approximately 1 dex. A similar result is found when comparing the broad H estimates for black hole mass based on velocity dispersion measurements.
To study whether the cause of the bias when using the broad H to estimate the mass is dust extinction, we compare the ratio between the FWHMs of broad Pa and H with the hydrogen column density and Balmer decrement ( Figure 16). Based on the limited number of cases, a trend can be seen between the FWHM ratio and the hydrogen column density, potentially indicating that the Paschen lines are broader than H for higher column densities. This lets us conclude that obscuration indeed causes the H line to be attenuated as the Paschen lines are less affected by reddening. Is the obscuration indeed due to dust? To understand any potential trend, however, a larger sample of Seyfert 1.9 galaxies would be necessary. The reason for this is most likely that the Balmer decrement is more complex: Pottasch (1960) noted that the Balmer line optical depth can also lead to larger Balmer decrements when the gas is optically thick, H is scattered, and H is absorbed and reemitted as Pa or H . As a consequence, H emission gets stronger and the Balmer decrement increases. This means that the use of the Balmer decrement as an indicator of reddening due to dust is potentially not valid.
To conclude our analysis on the use of the Paschen lines, we find that they provide reliable estimates of the black hole mass for Seyfert 1.9 and 2 galaxies assuming that velocity dispersion measurements of the black hole mass are robust. The use of broad H for mass estimation is already established for Seyfert 1 -1.8 galaxies (Mejía-Restrepo et al. 2016). For Seyfert 1.9 galaxies, the use of H for the mass estimation is shown to be less robust, as there is clearly an offset when comparing with measurements based on Paschen lines and stellar velocity dispersions. To fully quantify the bias when using H in the case of Seyfert 1.9 galaxies, we need more sources. Our analysis relies on sources for which the column density is determined in order to find cases with hidden broad lines (see Figure 14). With upcoming data releases from the BASS project, we will have more Seyfert 1.9 Our complete X-shooter sample includes 168 sources. However, for some sources, we do not have column density or Pa measurements, which reduces our sample size from 168 total to 68 Seyfert 1-1.9 galaxies. B .
galaxies to work with.

Outlook for Studies Using JWST
In this study, we provide the largest NIR spectroscopic census and legacy database for nearby AGN using the large collecting area of the VLT. The AGN luminosities of our sample ( bol ∼ 10 43 −5×10 45 erg s −1 ) are similar to the luminosities of AGN at the epoch of the peak of black hole growth at ∼ 1−2 (e.g., Aird et al. 2015). Our spectra thus provide a useful high-resolution, high-S/N template for higher-redshift AGN (∼ = 1−2). With the advent of JWST, unprecedentedly deep CL surveys will be possible. On board the satellite is the Near-infrared Spectrograph (NIRSPEC; Dorner et al. 2016), which is an NIR multi-object dispersive spectrograph. It operates in the 1−5 m regime and can simultaneously observe more than 100 slits. This large spectroscopic sample will have immense legacy value for NIRSPEC/JWST in the full 1-5 m range (z=1-3, ∼0.3-2 m rest frame). While the spectral resolution is slightly lower than that of X-shooter spectra ( ∼ 1000), it is still sufficient for resolving CLs (e.g. see Lamperti et al. 2017). An exposure time of 10 5 s is expected to yield an S/N = 3 sensitivity of ∼ 2 × 10 −19 erg cm −2 s −1 at 2 m. This is 1000 times more sensitive than our [Si ] 1.9640 observations (∼ 2 × 10 −16 erg cm −2 s −1 ). Translating this to 14-195 keV X-ray flux using the line ratio we find between [Si ] 1.9640 and the X-ray flux, the limit corresponds to a flux of ∼ 5 × 10 −15 erg cm −2 s −1 . This is ∼1000 more sensitive than the sensitivity limit of the 105 month deep Swift-BAT survey (the 105 month survey reaches >50% completeness at that sensitivity; Oh et al. 2018). Consequently, with NIRSPEC/JWST, it is potentially possible to observe highly obscured (Compton-thick) AGN missed by X-ray surveys as well as much fainter sources. As X-ray confusion can be a problem for low-luminosity AGN, it may be possible to detect them in the NIR with the CLs discussed in this paper, besides other NIR (high-ionization) lines (Satyapal et al. 2021). We note, however, that with greater sensitivity AGN CL emission may be difficult to distinguish from other emission mechanisms such as shocks (Rich et al. 2011).

SUMMARY AND CONCLUSIONS
In this work, we analyze 168 NIR spectra of nearby ( < 0.6) hard X-ray selected AGN from BASS. First, we look at high-ionization lines in the NIR spectrum of these nearby AGN: • We find CLs in more cases than found by previous studies. We find that 49/109 (45%) Seyfert 2 and 35/58 (60%) Seyfert 1 -1.9 galaxies show at least one NIR high-ionization line.
• • Studying the sources with hidden broad lines case by case, we find indications of galaxy-scale interactions and obscuration from extended dust lanes. The lack of broad optical emission line components can be explained by obscuration due to dust or gas in the environment of the host galaxy rather than by obscuration by the nuclear torus.
• NIR hidden broad lines can be used to estimate the black hole mass. Mass estimations using the FWHM of Pa and Pa are in accordance with estimations from velocity dispersion measurements. On the other hand, the H width underestimates the mass in Seyfert 1.9 galaxies.
This study provides a benchmark investigation of the use of CL emission as a tracer of AGN activity using the largest assembled NIR rest-frame sample to date. With next-generation NIR instruments, particularly JWST, deeper and more sensitive observation will be possible. As such, it will be possible to observe highly obscured (Compton-thick) AGN missed by X-ray surveys and much fainter sources, expanding our understanding of the AGN population.  Figure A.1 we show the full observed NIR X-shooter arm observations for a set of 12 sources from our sample. B. EMISSION LINE TABLE   In Table B.1 the emission lines in the NIR range of 0.9 -2 m are shown based on the fitting routine. The lines are sorted by increasing wavelength and spectral region. In Table B.1 the CLs that are part of this study are listed, including their IP and critical density values. C. MEASURED DATA Table C.1 includes the flux measurements from the spectral fits done in this work. The table is available in its entirety in the online journal. Table C.1 includes the flux values of the Pa spectral range; the other spectral ranges can be found in the online journal.

D. SPECTRAL FITS
As an example, Figure D.1 shows the different spectral regions where the emission lines have been fitted. The spectrum taken as an example is of the source 2MASX J214805.31-535941.3 (BAT ID 1604). For each object, we show the spectra (in black), the components fitted for the emission lines (in blue), the overall best-fit model (in red), and the residuals (below).
In Figure D.2 source ESO 103-035 (BAT ID 988) is shown, where the fitting routine is more difficult to apply due to irregular line shapes and heavy telluric absorption. In the Pa spectral range, a spline fit is applied to estimate the continuum level.            Comparison of the two continuum fitting methods. The black spectrum in the upper part shows the raw spectrum. The spline fit is indicated in blue and the fourth-order polynomial fit in brown. The fourth-order polynomial fit corrected spectrum is shown in the center (in brown) and the spline fit corrected spectrum is shown at the bottom (in blue). The red shaded regions are excluded from the continuum fit due to intrinsic emission.