© The Authors 2023

1. Introduction

Active galactic nuclei (AGN) emit over the entire electromagnetic spectrum and are powered by accretion onto a supermassive black hole (Rees 1984). The central engine emits UV photons which interact with energetic electrons in the so-called corona (Nandra & Pounds 1994), thus producing the X-ray emission (e.g., Haardt & Maraschi 1993). This emission interacts with the circumnuclear dust and gas, producing the obscuration observed in the spectra. The dust emission can be observed in the infrared energy range, as a result of the thermalization of the UV photons by the dust. In the optical, the nuclear emission will be obscured, removing the continuum and the broad components in the emission lines (Osterbrock 1981). On the other hand, the gas will absorb and scatter the X-ray continuum producing the X-ray absorption, most noticeably at energies below 10 keV in the X-ray spectrum (Brightman & Nandra 2011). We note that in general, the absorption in the optical and in the X-rays occurs concurrentl (Mainieri et al. 2007; Malizia et al. 2012; Merloni et al. 2014; Davies et al. 2015; Koss et al. 2017).

Obscuration gives evidence of material in the line of sight (LOS), which could be associated with the torus. Gas that is not in the LOS of the observer can also imprint some features on the X-ray spectrum. Between 10 keV and up to hundreds of keV, there is a reflection hump created by X-rays being reflected at the accretion disk (Fabian 2006) or more distant material, such as the torus (Brightman & Nandra 2011). Furthermore, the most robust emission line seen in the X-rays, the Fe Kα emission line (e.g., Fabian 2006), can be related to circumnuclear material, being broad and exhibiting relativistic effects due to its creation close to the supermassive black hole (SMBH), as well as narrow, presumably originating from more distant material. These reflection features are therefore a useful tool to study the configuration of the accretion disk and the torus. To better understand the properties of the reflector, many models have been developed, such as BORUS (Baloković et al. 2018), where the reprocessing medium is assumed to be a sphere with a conical cut-off at both poles, approximating a torus with variable covering factor; cTORUS (Liu & Li 2014), which is similar to BORUS but clumpy and with the half-opening angle of the torus fixed at 60 degrees; MYTORUS (Murphy & Yaqoob 2009) that proposes a toroidal geometry where the covering fraction is fixed to 0.5, RXTorus (Paltani & Ricci 2017) a model that assumes absorption and reflection from a torus with a varying ratio of the minor to major axis; and XILLVER (García et al. 2013) which calculates the reflected spectrum from the surface of an X-ray illuminated, ionized accretion disk by solving the equations of radiative transfer, energy balance, and ionization equilibrium in a Compton-thick and plane parallel medium.

It is not clear how the reflecting structure is formed, but clues can be gathered from the relation between reflection strength and the nuclear accretion rate. From the observational point of view, the torus in the infrared (IR) becomes weaker in the low-luminosity regime (i.e., for low accretion rates – below 10⁻³, González-Martín et al. 2017). Furthermore, in the X-ray, it has been seen that the Compton thin absorption (N_H <  1.5 × 10²⁴ cm⁻²) is less frequent in objects with low accretion rates: the fraction of Compton-thin obscured sources (10²² < N_H < 10²⁴ cm⁻²) decreases in the low-luminosity regime (Ricci et al. 2017b; Osorio-Clavijo et al. 2022), while the fraction of Compton-thick sources apparently remains constant. Both Compton thick and thin absorbers can produce reflection features, with different shapes and strengths. In addition, Osorio-Clavijo et al. (2022) found that in a sample of 81 AGN, ∼13% of the objects are lacking reflection signatures and the remaining galaxies (with a detected reflection component) are expected to be highly obscured. In the following, we attempt to measure the global distribution of gas around the nucleus, whether in the line of sight or otherwise, through their contribution to the reflection. In particular, we aim to establish whether the changes in the gas configuration become flattened or overall optically thinner as the accretion rate goes down.

Additionally, by modeling the X-ray reflection we are able to study the continuum emission, estimating the coronal parameters: the power law (Γ) and the high energy cut-off (E_cut). It has been shown that the slope of the power law depends on the accretion rate with changes at intermediate accretion rates (L_Bol/L_Edd = λ_Edd ∼ 10⁻³), pointing to a change in the accretion mechanism, for example between a corona on a thin disk to an advection-dominated accretion flow (ADAF, Narayan et al. 1994). The relationship toward low accretion rates is usually seen with a lot of scattering, which can be intrinsic or due to observational uncertainties (Shemmer et al. 2006; Gu & Cao 2009; Younes et al. 2011; Yang et al. 2015; She et al. 2018). Our second objective is to re-evaluate this relationship in the low accretion-rate range, through detailed modeling of the reflection and broadband X-ray data, using observations from XMM-Newton+NuSTAR+Swift).

This paper is organized as follows. In Sect. 2, we present details of the observations and sample. The data reduction is reported in Sect. 3. The methodology followed during this work is shown in Sect. 4. All the results are reported in Sect. 5. The implications of our X-ray spectral analysis are discussed in Sect. 6. Finally, a summary of our findings is presented in Sect. 7.

2. Sample and data

Hard X-rays (E ≥ 10 keV) are not significantly affected by obscuration, at least up to N_H ∼ 10²⁴ cm⁻² (Ricci et al. 2015), which allows us to obtain a highly complete AGN sample. Our work focuses on LLAGN selected through their hard-band X-ray emission as identified in the Swift/BAT 70-month catalogue (Baumgartner et al. 2013) on board the Neil Gehrels Swift Observatory (Gehrels et al. 2004). In particular, BAT operates in the 14–195 keV energy band, and the BAT AGN Spectroscopic Survey (BASS) provides high-quality multi-wavelength data for the BAT AGN, including black hole mass measurements (Koss et al. 2017) and X-ray spectroscopy modeling (Ricci et al. 2017a). The first data release (DR1) of the BASS project (Koss et al. 2017) includes 642 of Swift/BAT AGN and the second release of optical spectroscopy (BASS/DR2) will also soon be publicly available (Koss et al. 2022; Oh et al. 2022).

Our sample of galaxies was selected from the BASS/DR2 with accretion rates log(λ_Edd)≤ − 3.0 obtaining in total a sample of 24 AGN. We used the HEASARC¹ archive to search simultaneous and not simultaneous NuSTAR and XMM-Newton data public until August 2020. This analysis has provided data with both telescopes for 16 sources. We include the proprietary data of the galaxy NGC 5033 (PI: Diaz Y.; log(λ_Edd)= − 4.0), an AGN also contained in the BASS/DR2.

Our final sample of LLAGN contains 17 objects, 11 of which are classified as Seyfert 2 (i.e., only narrow lines are visible in the optical spectrum), and six are classified as Seyfert 1.9 (a broad component is visible in Hα but not in Hβ) in the BASS/DR2. Table 1 shows the general properties of our sample. Notes for the individual galaxy are given in Appendix A and Table D.1 shows the log of the observations.

3. Data reduction

The data reduction process is described below. Details on the observations can be found in Table D.1.

3.1. XMM-Newton data

This satellite has two X-ray instruments, a grating spectrometer, and the European Photon Imaging Camera (EPIC). The EPIC instrument has three detectors, two MOS (Turner et al. 2001) and one PN CCDs (Strüder et al. 2001), we only used the observations from the EPIC-PN because of its higher throughput (Strüder et al. 2001) and because of inclusion of the EPIC-MOS spectra resulted in too much statistical weight attributed to the low energy range data points, as compared to the NuSTAR and Swift/BAT data. We processed the Observation Data Files (ODFs) from the European Photon Imaging Camera (EPIC) PN detector using the Science Analysis System (SAS version 17.0.0). We followed standard procedures to obtain calibrated and concatenated event lists, filter them for periods of high background flaring activity, and extract light curves and spectra. Source events were extracted using a circular region of 49 arcsec centered on the target, and background events were extracted from a circular region of 98 arcsec on the same chip far from the source. We verified the photon pile-up is negligible in the filtered event list with the XMMSAS task EPATPLOT. We generated response matrix files (RMFs) and ancillary response files (ARFs) and rebinned the spectra in order to include a minimum of 25 counts in each background-subtracted spectral channel and to not oversample the intrinsic energy resolution by a factor greater than 3.

3.2. NuSTAR data

Nuclear Spectroscopic Telescope Array (NuSTAR) was successfully launched in 2012 June (Harrison et al. 2013). It has two identical co-aligned telescopes, each consisting of an independent set of X-ray mirrors and a focal-plane detector, referred to as focal plane modules A and B (FPMA and FPMB), which operate in the energy range of 3–79 keV. The data reduction was performed with NUSTARDAS V1.6.0, available in the NuSTAR data analysis software. The event data files were calibrated with the NUPIPELINE task using the response files from the Calibration Database CALDB v.20180409 and HEASOFT version 6.25. With the NUPRODUCTS script, we generated both the source and background spectra, plus the ARF and RMF files. For both focal plane modules (FPMA, FPMB), we used a circular extraction region of radius 49 arcsec centered on the position of the source. The background selection was made, taking a region free of sources of twice the radius of the target and located in the same detector quadrant. The spectral channels were grouped with the FTOOLS task GRPPHA to have a minimum of 20 counts per spectral bin in the 3.0–79.0 keV energy range.

3.3. Swift data

The Neil Gehrels Swift Observatory was launched on November 20, 2004. It carries three instruments to enable the most detailed observations: the Swift Burst Alert Telescope (BAT; Barthelmy et al. 2005; bandpass: 15–350 keV), the X-ray Telescope (XRT; Burrows et al. 2005; bandpass: 0.3–10 keV), and the UV/Optical Telescope (UVOT 170–650 nm). In this work, we focus on the Swift/BAT and Swift/XRT instruments.

3.3.1. Swift/BAT

We retrieved the binned and calibrated spectra, together with the response matrices for our targets, from the Swift/BAT 105-month All-sky Hard X-Ray Catalog reported in Oh et al. (2018). The observations were taken with the Burst Alert Telescope (BAT) on board the Swift observatory. This survey has a sensitivity of 8.4 × 10⁻¹² erg s⁻¹ cm⁻² in the 14–195 keV bands over 90% of the sky, with eight-channel spectra averaged over the 105-month duration of the survey. The complete analysis pipeline is described in the Swift/BAT 22 All-sky Hard X-Ray Survey (Tueller et al. 2010).

3.3.2. Swift/XRT

For three sources (NGC 7674, ESO 253-G003, and IGRJ 11366-6002) there are no simultaneous XMM-Newton and NuSTAR observations. We explored the Swift/XRT archived and found simultaneous observations for the three of them. The data reduction of the Swift/XRT in the Photon Counting mode was performed by following standard routines described by the UK Swift Science Data Centre (UKSSDC) and using the software in HEASoft version 6.30.1. Calibrated event files were produced using the routine XRTPIPELINE, accounting for bad pixels and effects of vignetting, and exposure maps were also created. Source and background spectra were extracted from circular regions with 25 arcsec and 50 arcsec radius. The XRTMKARF task was used to create the corresponding ancillary response files. The response matrix files were obtained from the HEASARC CALibration DataBase. The spectra were grouped to have a minimum of 20 counts per bin using the GRPPHA task.

Following the data reduction, we find that both NGC 7674 and ESO 253-G003 present very low counts, preventing us from doing a proper spectral fit. Thus, these are not used in the analysis.

4. Methodology

The analysis of the data comprises two steps: (1) a combination of XMM-Newton and NuSTAR observations and (2) a homogeneous spectral fitting of the sample. All the spectra have been fitted using xspec version 12.10.0 (Arnaud 1996) and all the errors reported throughout the paper correspond to 90% of confidence level.

4.1. Combination of the XMM-Newton and NuSTAR observations

In this work, we use NuSTAR observations, covering an energy range from 3 to 79 keV, that are vital to studies of the Compton hump as a key signature of the reflection. Additionally, we have XMM-Newton data, which provide the best combination of sensitivity, bandpass, and spectral resolution at energies ranging from 0.5–10.0 keV. Objects with simultaneous observations with XMM-Newton and NuSTAR were fitted with all model parameters tied between the different spectra, except for a free cross-normalization factor. Objects with non-simultaneous observations (denoted with the symbol * in our work) were tested for spectral variability between the observation epochs. In order to detect the spectral variability, we simultaneously fitted the XMM-Newton + NuSTAR spectra in the overlapping 3.0–10.0 keV range for each object with a power-law model under neutral absorption. In cases where the spectrum was not well-fitted with this model, we added a Gaussian component centered at 6.4 keV and studied the improvement of the fit.

At first, all the parameters were tied between the spectra of the different epochs and instruments. When this model produced a satisfactory fit (${\chi }_{\nu }^2\le 1.2$), the source is considered to be non-variable and treated in the same way as the objects with simultaneous observations. Two objects are well-fitted with a tied normalization of the power law (LEDA 96373* and NGC 2273*). For the remaining objects, NGC 4258*, NGC 3147*, NGC 2110* and IC 4518A*, we found that the normalization of the power law varying between epochs resulted in a satisfactory fit. For these objects, in the subsequent fitting, the normalization of the power law was left free between the epochs, but the remaining parameters were tied. For one object, ESO 253-G003*, allowing the slope of the power law to vary freely improved the fit significantly according to the F-test. Given the spectral variability in this source, we had to leave most parameters untied between the epochs, and, therefore, the inclusion of the lower energy spectrum would not constrain the reflection model further. For this reason, in this source, we used the high-energy spectra only. Finally, one object (NGC 7674) could not be fitted well with a free slope and normalization, as this is a known changing-look AGN and its spectrum changed significantly in shape between observations (Bianchi et al. 2005), so we retained only the NuSTAR spectrum for the following analysis. The best model and the final configuration for each object are summarized in Table D.2.

4.2. Spectral analysis

The aim of our spectral analysis is to quantify how much reflecting material can possibly exist around the central engine. For this purpose, we allow for different types of reflectors and maximize the freedom of the fitted parameters of all required models, even if the fit results in unconstrained values for many of them. We tested the significance of the reflection component in the 3–195 keV band in all objects, returning significant detections in most of them. These tests are summarized in Appendix B.

For all spectral fits, we included a multiplicative constant normalization between FPMA, FPMB, EPIC-PN, and Swift/BAT to account for calibration uncertainties between the instruments. We started with a baseline model and added different components until a satisfactory fit was obtained. We have selected three broad components in order to parametrize the following three scenarios:

The first is the cut off power law model (cPL) obscured by neutral material, whereby a single power-law model, which corresponds to the primary emission of a non-thermal source. The column density, N_H, los, is added as a free parameter to take the absorption by matter along our LOS to the target into account. The free parameters in this model are the column density, N_H, los, the slope of the power law, Γ, the high energy cut-off, E_cut and the normalization.

The second comprises a series of reflection models (Refl): when the X-ray continuum is scattered by the surrounding gas, it can produce fluorescent emission lines (most notably Fe Kα 6.4 keV) and a broad hump-like continuum peaking around 10–30 keV.

Three further potential scenarios form the basis for the reflection models (1) A neutral reflector with a semi-infinite column density modeled with PEXMON (Nandra et al. 2007). (2) A smooth spherical distribution of neutral gas, with conical cavities along the polar direction, is modeled with BORUS (Baloković et al. 2018); (3) accretion disk reflection modeled with XILLVER (García et al. 2013).

The neutral reflector with a semi-infinite column density modeled with PEXMON assumes the existence of optically thick and cold material, distributed in a slab and covering a given fraction of the X-ray source. The PEXMON model includes fluorescence, adding some spectral features, such the emission lines FeKα and Fekβ, following Monte Carlo calculations from George & Fabian (1991). This model represents both the reflected and intrinsic emission defined with Γ and the high energy cut-off (E_cut), as well as the reflection fraction, R_f. The free parameters in this model are the reflection fraction, R_f (to account for the reflection component and the contribution from the intrinsic power-law continuum), the spectral index, Γ, the high energy cut-off, E_cut, the inclination and the normalization.

The smooth spherical distribution of neutral gas, with conical cavities along the polar direction, modeled with BORUS calculates the reprocessed continuum of photons that are propagated through a cold and static medium. BORUS is similar to the torus model BNtorus of Brightman & Nandra (2011), but it has additional free parameters (E_cut, A_Fe), additional chemical elements included, calculation extending to higher energies and the LOS component separated out. Furthermore, this model has a variable covering factor which is an advantage compared with other models, as MYTORUS (Murphy & Yaqoob 2009) that proposes a toroidal geometry where the covering fraction is fixed to 0.5. In this work, we used the geometry of a smooth spherical distribution of gas, with conical cavities along the polar directions (borus02). The column density and the inclination of the torus are free parameters in this model. borus02 includes fluorescent emission lines, according to fluorescent yields for Kα and Kβ lines from Krause & Oliver (1979), for all elements up to zinc (Z <  31). The reflected spectrum of this torus is calculated for a cut-off power law illuminating continuum, where E_cut, Γ and normalization are free parameters. We modeled the direct coronal emission separately with a cut-off power law under a neutral absorber as described above. We have set as free parameters the column densities along the LOS, N_H, los, the inclination, Cos(θ_incl), the covering factor, CF, the column density of the reflector, log(N_H, refl), the spectral index of the primary emission, Γ, the high energy cut off, E_cut, and the normalization of the reflector tied to the primary emission.

In the accretion disk reflection modeled with XILLVER, the coronal spectrum is a power law with an exponential cut-off described by the photon index, Γ, and the high energy cut-off, E_cut. Another important parameter is the ionization parameter, ξ, defined as the incident flux divided by the density of the disk. This parameter is described by log(ξ) ranging from 0 for a neutral disk to 4.7 erg cm⁻² s⁻¹ for a heavily ionized disk (see García et al. 2013, for a more detailed description). Other parameters in this model are the iron abundance, A_Fe relative to the solar value (assumed to be solar in this work), redshift, reflection fraction, R_f, and the inclination. Also, this model takes into account both the reflected continuum and the FeKα. The free parameters in this model are the spectral index, Γ, the high energy cut-off, E_cut, the ionization degree, log(ξ), the inclination, incl, the reflection fraction, R_f (to normalize the reflection component relative to the intrinsic power-law continuum) and the normalization.

Finally, when the combination of the models described above does not produce a good fit, we explore if the addition of spectral component(s) improves the fit. For the model based on the soft X-ray emission (SE), the following three spectral components are considered:

An absorbed scattered power law: an absorbed power-law PL to model the scattered emission that is deflected by ionized gas. The photon index, Γ, of the scattered component is tied to the primary power law. We set as free parameters the column density, N_H, ext and the normalization of the scattered component but restricted to be less than 5% of the main one.

Thermal emission: An optically-thin thermal component, modeled by MEKAL in xspec, to model the soft excess observed below 1 keV, and potentially due to either star formation processes and/or thermal emission from a hot interstellar medium. We kept the hydrogen column density, abundance, and switch at their default values (1, 1, and 1 respectively) and we let the temperature, ionization, and normalization free to vary.

An ionized absorber (ab): a warm absorber was modeled with zxipcf within xspec. This model uses a grid of xstar photoionized absorption models (calculated assuming a microturbulent velocity of 200 km s⁻¹) for the absorption, and it assumes an absorbent covering some fraction of the source, cf_W (Reeves et al. 2008). zxipcf has as free parameters the column density, N_H, W, the ionization state, log(ξ_W), the covering fraction, cf_W, and redshift. We set the covering fraction to cf_W = 1 to mimic an absorber covering all the sky. We let as a free parameter N_H, W and log(ξ_W).

We started our analysis by fitting a baseline model that is defined as MOD = Refl + cPL to the data. Then we added one SE emission or absorption component (we tested one by one: MOD + PL, MOD + MEKAL and MOD*ab) and explored whether the inclusion of these components improves the fit (using an F-test in the case of the scattered power law and MEKAL, and we evaluated the improvement of the fit using the values of the χ² and a visual inspection of the residuals in the case of the ionized absorber). If any of the improvements was significant, we selected the model that returned the lowest value of χ²/d.o.f.² and chose it as the new baseline model and the process of including and testing an additional SE component was repeated. When none of the additional SE components provided a significant improvement, the iteration stopped. Up to four iterations were necessary for each object and reflection model. The method is represented in Fig. 1. The process was repeated separately for each reflection model. Thus, we report up to three best-fitting models for each object.

Fig. 1. Schematic view of the methodology followed to fit the data. Note that the loop in blue iterate a maximum of four times. For a detailed explanation of the method, we refer to the main text.

The models that were selected to fit the data are represented in xspec as:

$$\begin{array}{c}\hfill C\times N_{\mathrm{H},\mathrm{Gal}}\times ab\times (N_{\mathrm{H},\mathrm{ext}}\times SE+N_{\mathrm{H},\mathrm{los}}\times \text{cPL}+N_{\mathrm{H},\mathrm{los}}\times \mathrm{Refl})\end{array}$$

where C represents the cross-calibration constant between different instruments, N_H, Gal is the Galactic absorption (phabs in xspec) predicted using N_H tool within FTOOLS (Dickey & Lockman 1990; Kalberla et al. 2005). “ab” is the ionized absorption component modeled with zxipcf - in cases where this component is used, otherwise it is equal to unity. Two absorbing column densities are used, which will be called here N_H, ext and N_H, los (zphabs in xspec). N_H, los is assumed to cover the nuclear components (power law and disk reflection)³ and N_H, ext covers the SE component⁴. Moreover, cPL is a cutoff power law (cutoffpl in xspec) representing the primary X-ray emission and “Refl” represents the different reflection models used.

We note that we imposed the following conditions to the resulting best-fit Γ >  0.5, N_H, Gal ≤ N_H, ext and N_H, los >  N_H, ext. In the case of NGC 1052, additional Gaussian lines were required at soft energies from a visual inspection, we included S XIV at 2.4 keV and Si XIII at 1.85 keV, also in agreement with (Osorio-Clavijo et al. 2020). They were added as a narrow Gaussian line with fixed centroid energy and a width fixed at 0.01 keV.

5. Results

We refer to the following sections and tables for details on the analysis: comparison with previous works and our results on individual objects can be found in Appendix A. The coronal parameters (i.e., Γ, E_cut, and χ²) are listed in Table D.3. The reflection parameters, namely, R_f and inclination for PEXMON; log(N_H, refl), CF and inclination for borus02; log(ξ), R_f and inclination for XILLVER are listed in Table D.4. In Table D.5 we show all the additional components required for the fit with each of the reflection models, that is, the column density of the neutral absorbers in the LOS to the extended and nuclear components and the temperature of the optically thin thermal emission components. Additional parameters, namely, the column density and ionization parameter of the ionized absorbers in the LOS and the normalization of the scattered power law can be seen in Table D.6. All cross-calibration constants are listed in Table D.7. The plots of the spectra with the best-fit models and their residuals can be found in Appendix E.

5.1. Models

In this work, we used three reflection models (PEXMON, borus02, and XILLVER) that were used to fit the spectrum of each of the sources in the sample, namely, each of the sources is fitted by three different models.

The simplest model used in our work (PEXMON) is a good representation of the data, however, we focus on models that can explore different reflector geometries as borus02 and XILLVER. To decide which model provides the best description of the observations, we estimate the “evidence ratio” using the Akaike information criterion (AIC) for both models. This evidence ratio allows us to compare if one model is better than another one, it is defined using as ϵ = W(AIC_torus)/W(AIC_disk) where W(AIC_torus) and W(AIC_disk) are the “Akaike weight” (see Emmanoulopoulos et al. 2016 for more details). The evidence ratio is a measure of the relative likelihood of the torus versus the disk model. The torus model is 200 times more likely than the disk model when ϵ ≤ 0.0067. The disk model is 200 times more likely than the torus model when ϵ ≥ 150. The evidence ratio is listed in Table 2.

For nine (53%) objects (NGC 4258,NGC 1052, NGC 2110, LEDA 96373, NGC 2992, M 51, HE 1136−2304, IC 451A, and NGC 5033) borus02 is preferred. Then in the following sections, we describe how we chose this model as the best representation of the data in these objects. On the other hand, two (12%) objects (NGC 2273 and NGC 7674) are well fitted with a disk (XILLVER) model and in six (35%) objects (NGC 3998, NGC 3718, ESO 253-G003, NGC 2655, NGC 3147, and IGRJ 11366-6002) both models fit similarly well the data.

Since it is not possible to distinguish between a reflection dominated by a torus or a disk, we will treat them separately in the following sections. When referring to the torus case, we refer to all borus02 models for the whole sample (15 galaxies: indistinguishable and distinguishable cases). In the disk case, we refer to the XILLVER model, considering the indistinguishable cases and the cases where the disk is a good representation of the data (eight galaxies in total), since in the other cases a torus model is not a good representation of the data.

5.2. X-Ray continuum properties

The X-ray continuum of AGN is described by a power law with a high-energy cutoff. The free parameters of this component are the spectral index (Γ) and the high energy of the cut-off (E_cut). In Fig. 2, we show the histogram derived from the broadband spectral analysis with each reflection model. We note that the spectral index between both the reflected and intrinsic emission are tied. We find that the mean values (dashed vertical lines) of Γ for the sample using PEXMON, borus02 and XILLVER are consistent (1.73 ± 0.21, 1.72 ± 0.17 and 1.72 ± 0.20 respectively). We note that the simplest model used in our analysis is PEXMON, which shows values of the photon index consistent with more geometrical models.

Fig. 2. Comparison between the spectral index estimation, Γ, between the models. The dotted lines represent the mean values.

Considering the torus model (15 galaxies), we found median values of Γ = 1.76 and σ = 0.16, with values ranging between [1.40, 2.06]. Another important parameter that could be estimated with the NuSTAR data is the high energy cut-off (E_cut). This parameter can be considered as an indicator of the temperature of the X-ray corona. Consequently, its knowledge provides information about the dynamics of the corona and the physical processes occurring within it. Nevertheless, this parameter is poorly constrained. A lower (upper) limit of E_cut could be determined for eight (two) sources. The five AGN for which E_cut could be determined (NGC 3998, ESO 253-G003, NGC 2110, NGC 2992, and NGC 5033) have a mean value of E_cut = 193.28 keV, with a standard deviation of σ = 99.19 keV.

Furthermore, taking into account the disk model (eight galaxies), we found median values of Γ = 1.71 and σ = 0.23, with values ranging between [1.40, 2.06]. Regarding the high energy cut-off, we could find five (one) lower (upper) limits and for six objects we obtained a mean value of E_cut = 371.47.58 keV with a standard deviation of σ = 619.99 keV.

5.3. Soft band spectral fit

In the soft (0.3–10.0 keV) energy band, we added a thermal (MEKAL in xspec), scattered power-law absorption by ionized gas (also referred to as “warm absorption”), modeled with zxipcf in xspec, or a combination of these components to improve the spectral fit.

When considering the torus model (15 objects), six objects, NGC 3998, ESO 253-G003, NGC 3147, M 51, IGRJ 11366-6002, and NGC 5033, do not require an additional component to improve the fit. Two objects (NGC 3718 and HE 1136-2304) required a MEKAL component (with kT = 0.88${}_{0.67}^{\ast \ast }$ keV and kT = 0.59${}_{0.51}^{0.67}$ keV, respectively). Two objects are well-fitted with a combination of MEKAL, a power law, and warm absorber (NGC 1052 with MEKAL+PL and NGC IC 451A with MEKAL*ab). Composite models are needed for five galaxies (NGC 4258, NGC 2655, NGC 2110, LEDA 96373, and NGC 2992). In the cases where two MEKAL were required, the values of the temperatures are in the range of kT₁ = [0.58–0.62] keV, with a mean value of kT₁ = 0.60 keV and σ = 0.02 keV, and then kT₂ = [0.15–0.22] keV, with a mean value of KT₂ = 0.19 keV and σ = 0.03 keV. The mean value of the ionized absorber is N_H, W = 1.66 × 10²² cm⁻² and σ = 1.41 × 10²² cm⁻² The degree of ionization is in the range [−1.14, 4.30] with the mean log(ξ_W)=1.31 and σ = 1.99.

In relation to the disk model (eight objects), five galaxies do not require any component to improve the spectral fit: NGC 3998, ESO 253-G003, NGC 3147, IGRJ 11366-6002, and NGC 7674. One galaxy, NGC 3718, requires a MEKAL component to improve the fit. Two galaxies, NGC 2655 and NGC 2273, are well-fitted with a composite model, MEKAL*ab.

5.4. Line-of-sight column density

Absorption of X-rays by neutral material is the result of the combined effect of Compton scattering and photoelectric absorption. The Compton scattering and the photoelectric absorption were modeled using CABS and ZPHABS in xspec respectively. In ZPHABS, we fixed the redshift at the value of each source. The only free parameter is the column density, which is tied in all the fits (i.e., N_{H − ZPHABS} = N_{H − CABS} = N_{H − los}).

According to the torus model, we can classify six galaxies as unobscured (log(N_H, los)< 22), namely, NGC 3998, NGC 3147, NGC 2992, HE 1136–2304, IGRJ 11366-6002, and NGC 5033, with values in the range between log(N_H, los)=[20.0, 21.89] and then eight galaxies, NGC 3718, NGC 4258, ESO 253-G003, NGC 1052, NGC 2655, NGC 2110, IC 451A, and LEDA 96373, as obscured (22 <  log(N_H)< 24.18), with values aross the range of log(N_H, los)=[22.01, 24.09]. According to our spectral analysis, one galaxy (M 51) in our sample can be classified as Compton-thick (CT), using as a threshold N_H = 1.5 × 10²⁴ cm⁻², or log(N_H, los)=24.18. The mean values of the spectral index, column density in the LOS, and column density of the torus are reported in Table 3. All the parameters are consistent between the groups. We note that the values of the cross-calibration constant between the groups are consistent.

Regarding the disk model, two galaxies, NGC 3998 and IGRJ 11366-6002, can be classified as unobscured. Six galaxies, NGC 3718, ESO 253-G003, NGC 2655, NGC 3147, NGC 2273, and NGC 7674, as obscured with values log(N_H, los)=[22.03, 23.36]. The mean values of the spectral index, the column density in the LOS, the ionization degree of the accretion disk, and the reflection fraction are reported in Table 3 and showed values consistent between the categories. We note that according to a reflection dominated by an accretion disk, none of the galaxies in our sample can be classified as CT.

5.5. The reflection component

The reflection features observed in the hard X-ray spectra of AGN may be caused by neutral and distant material such as the torus or by the ionized material of the accretion disk.

For the case where the reflection is dominated by the torus, the mean value of our sample for the column density for this structure is log(N_H, refl)=23.69 and σ = 0.76 with values between [22.50, 25.40]. Four objects, NGC 1052, M 51, IGRJ 11366-6002, and IC451A, show a column density of the torus consistent with a Compton-thick structure. Another important parameter derived from the torus reflector model is the covering factor. We were only able to determine a lower (upper) limit for three (six) objects (with a lower limit for NGC 3998, NGC 2655, and IC451A, and an upper limit for NGC 3718, NGC 4258, ESO 253-G003, NGC 3147, NGC 2992, and M 51). This parameter was determined for six (40%) objects (NGC 1052, NGC 2110, LEDA 96373, HE 1136-2304, IGRJ 11366-6002, and NGC 5033) with a mean value of CF =0.59 and σ = 0.25. The half-opening angle of the polar cutouts, cos(θ_incl), is also measured with the torus model. However, we obtain only lower (upper) limits for seven (three) objects and properly constrained it for five sources.

For the disk-like reflection, we constrained the value of the ionization degree of the accretion disk in six galaxies, NGC 3718, ESO 253-G003, NGC 3147, NGC 2273, IGRJ 11366-6002 and NGC 7674, and we found median values of ionization degree of the disk of log(ξ)=2.44 and σ = 0.81. In one (one) objects, we only obtained an upper (lower) limit, namely: NGC 3998 with the lower limit and NGC 2655 with the upper limit. Regarding the reflection fraction, R_f, we performed a test by fixing the value of log(ξ) to 0 and comparing the reflection fraction obtained with this XILLVER configuration and PEXMON. For the sample, we found consistent values between both models. However, since our goal is to restrict the accretion disk features, we left the ionization degree as a free parameter. We obtain six lower limits, one upper limit (NGC 2655) and it is constrained in one case (NGC 7674). This model also allows us to estimate the inclination, and this parameter is constrained in five sources with mean value Incl=70.10 deg and σ = 27.23 deg and four upper limits.

5.6. Flux and luminosity

We computed the X-ray flux and luminosity in two energy bands: 2.0–10.0 keV and 10.0–79.0 keV using xspec. We note that the redshift of the sources was taken from NASA/IPAC Extragalactic Database (NED). The values can be seen in Table D.8. Taking into account a torus model, the mean value of the intrinsic luminosities in the borus02 case are log(L_{2.0 − 10.0})=41.73 and σ = 1.16 log(L_{10.0 − 79.0})=42.14 with σ = 1.29. In the XILLVER case, we found log(L_{2.0 − 10.0})=41.57 with σ = 1.06 and log(L_{10.0 − 79.0})=41.85 with σ = 1.24, and they are equivalent. The distribution of the intrinsic luminosity obtained in both cases can be seen in Fig. 3.

Fig. 3. Histogram of the intrinsic luminosity in the band 2.0–10.0 keV (left) and 10.0–79.0 (right) in the borus02 and XILLVER with number of objects in each group.

6. Discussion

We performed the X-ray spectral analysis of an AGN sample with accretion rates, log(L_Bol/L_Edd)≤ − 3 selected from the BASS/DR2 with available NuSTAR + XMM-Newton + Swift data. Models from a neutral reflector (PEXMON), reflection from an ionized accretion disk (XILLVER) and from the torus (borus02) have been used to fit the data. This sample is composed of 17. Our main results are summarized as follows: (1) In our sample, six (35%) objects are equally well-fitted with a disk or with a torus-like reflector. For nine (53%) galaxies, the torus reflection model is the best representation of the data. In two cases (12%), the disk model fits the data well; (2) When modeling the reflection with borus02, seven objects are well-fitted by a single neutrally-absorbed cutoff power law plus reflection (i.e., no components are required in the soft band). When modeling the reflection with XILLVER instead, five objects can be well modeled in the same way. The remaining objects require the addition of a MEKAL and/or scattered power law, an ionized absorber, or a combination of two or more of these components; (3) According to the torus model, six sources can be classified as unobscured (log(N_H)< 22), eight galaxies as obscured (22 <  log(N_H)< 24.18), and one object have a column density in the LOS consistent with a Compton thick source (log(N_H)> 24.18). According to the disk reflection, two (six) objects can be classified as unobscured (obscured). These classifications are consistent among the models, except in the case of NGC 3147 (unobscured according to the torus and obscured with the disk).

The high quality and broad spectral coverage available combining XMM-Newton+NuSTAR+Swift allowed us to put constraints on spectral parameters related to the accretion mechanism and reflection of LLAGN. Our analysis covers energies above 10.0 keV, where the reflection has an important role in the spectral fit, and considering this feature in the X-ray spectral analysis can affect the estimation of the coronal parameters (see Diaz et al. 2020). In the following, we discuss the physical interpretations of the results presented in this paper.

6.1. Determination of the L_Bol/L_Edd

The selection of the sample presented in this work was based on sources with log(L_Bol/L_Edd)≤ − 3 according to those values reported in BASS/DR2. However, because variability is one of the properties that characterize AGN, we estimate these accretion rates using the data analyzed here.

To estimate L_Bol/L_Edd, we followed the relation given in Eracleous et al. (2010), which uses the black hole mass and bolometric luminosities. According to Koss et al. (2017), the black hole masses available for the BASS sources were determined using different methods. For 14 of our sources, these were estimated using the velocity dispersion method, from the M_BH − σ_* relation by Kormendy & Ho (2013). Two galaxies have M_BH taken from the literature (NGC 3998 via the M-σ relation and NGC 4258 by a rotating H₂O maser disk), and for one source, it was estimated from the MgII emission line (ESO 253-G003). The uncertainties on these M_BH determinations are ∼0.3–0.4 dex (as explained in the BASS/DR1 paper – Koss et al. 2017). We conservatively assume that the typical uncertainty on M_BH is 0.4 dex.

The other key parameter is the bolometric luminosity and the best method for estimating it is through the integrated area under the spectral energy distribution (SED). However, observations from a variety of telescopes are necessary in order to build up a complete and detailed SED. The bolometric correction is another way to estimate it, which depends on the X-rays luminosity, and it is the one we will be using in this work. For instance, BASS/DR1 (Koss et al. 2017) focused on the bolometric correction derived by Vasudevan & Fabian (2007, 2009), where L_bol/L_{2 − 10keV} = 20 for L_Bol/L_Edd ≤ 0.4, and L_Bol/L_{2 − 10keV} = 70 for L_Bol/L_Edd ≥ 0.4.

As the bolometric luminosity is fundamental in the estimation of the accretion rate, we examined an alternative determination of L_Bol/L_Edd based on the available X-ray luminosity estimated using the data of our sample of AGN.

We also used the intrinsic luminosity in the 2–10 keV rest-frame energy range of L_{2.0 − 10.0keV}, derived from the best-fitting spectral models of the X-ray data. We note that in the case of indistinguishable cases, we used the values from the borus02 model. When using the XILLVER model, the results are the same. A comparison between our L_{2.0 − 10.0keV} calculation and the BASS/DR2 is presented in Fig. 4, showing differences in the luminosity, possibly related to the variability. The difference between the fluxes measured by BASS/DR2 integrated 70 months and the flux measured in the short exposures with NuSTAR that we use here can be quite large for highly variable sources such as NGC 2992 (Gilli et al. 2000; Shu et al. 2010; Hernández-García et al. 2017; Marinucci et al. 2018, 2020) and LEDA 96373 (Landi et al. 2009).

Fig. 4. Intrinsic luminosity in the 2.0–10.0 keV range from BASS/DR2 and our work. Dotted black line represents x = y.

To remain consistent between the state of each AGN when measuring Γ and other parameters, we recalculated L_Bol/L_Edd using the fluxes measured here and refined it by changing the bolometric correction as described below.

Furthermore, we use our L_{2.0 − 10.0keV} calculation in combination with the bolometric correction K(2.0–10.0 keV) from Duras et al. (2020), who used a sample of ∼1000 type 1 and type 2 AGN from five different AGN surveys for which they performed a SED -fitting. They reported a bolometric correction as a function of 2.0–10.0 keV X-ray luminosity. The resulting K(2.0–10.0 keV) are slightly smaller than those used previously (L_bol/L_{2 − 10keV} = 20), with a median value of K(2.0–10.0 keV)= 15.60 with a scatter of ∼0.37 dex (Duras et al. 2020). The values of bolometric luminosity and Eddington ratio are given in Table D.9. The errors in the bolometric luminosity correspond to the error propagation of M_BH (0.4 dex), K(2.0–10.0 keV) (0.37 dex), and L_{2 − 10keV} (estimated with xspec). In the following analysis, we use these L_Bol/L_Edd values to minimize the effects of source variability.

6.2. Accretion mechanism: Γ versus L_Bol/L_Edd relation

It has been suggested that the accretion mechanism in LLAGN (L_Bol/L_Edd <  10⁻³) is different from that in more powerful AGN (e.g., Seyferts) and similar to that of X-ray binaries (XRB) in their low or hard states (Yamaoka et al. 2005; Gu & Cao 2009; Younes et al. 2011; Xu 2011; Yuan & Narayan 2014; Hernández-García et al. 2016).

Some authors, following the relations obtained for XRB, have studied the accretion mechanism using the relation between the spectral index Γ and the accretion rate λ_Edd, finding a positive correlation between these quantities at high accretion rates and thereby suggesting a geometrically thin and optically thick disk, known as the standard model for accretion disks (Shakura & Sunyaev 1973; Koratkar & Blaes 1999). A negative correlation has also been found at low accretion rates, indicating radiatively inefficient accretion (e.g., Yuan et al. 2007). In this configuration, the accretion disk becomes truncated near the SMBH, with a geometrically thick and optically thin disk at lower radii and a thin disk at higher radii. However, these correlations show a large scatter (Shemmer et al. 2006; Gu & Cao 2009; Younes et al. 2011; Yang et al. 2015; Trakhtenbrot et al. 2017; She et al. 2018), with Γ values between [1,3] (Gu & Cao 2009; Younes et al. 2011) and [0.5, 3.5] (She et al. 2018). The high scatter in the spectral index estimate is still not understood – it could be due to the sensitivity of the measurements or to the intrinsic properties of the galaxies.

Thanks to the excellent statistics of NuSTAR in combination with XMM-Newton, we were able to better constrain the spectral index Γ in our low accretion rate sample. In Fig. 5, we show the relation between Γ and λ_Edd using the best fitting reflection model (borus02). We added data from Esparza-Arredondo et al. (2021), who studied the torus configuration of 36 AGN, using NuSTAR and Spitzer data and estimated the spectral parameters using the same reflection model used in this work (borus02). We applied the same bolometric correction to these data (see Sect. 6.1). In Fig. 5, the blue points correspond to this work and the light yellow stars represent the data points from Esparza-Arredondo et al. (2021).

Fig. 5. Correlation between the spectral index, Γ, from individual fits, vs. the Eddington ratio, log(λEdd)=log(LBol/LEdd), for our sample of galaxies of the best fit models. The dotted and dashed green line represent the relation given by Gu & Cao (2009), the orange dotted represents Younes et al. (2011), the magenta dashed line is the relation obtained by She et al. (2018), while the solid black line is the correlation obtained in this work. The purple dashed line corresponds to the relation found by Fanali et al. (2013). The blue points represented the binned data. The pink points and light yellow stars are the data point of the best fit model in this work and the ones obtained by Esparza-Arredondo et al. (2021).

In order to check whether a break exists in this relation, as in XRB, we sought to test different scenarios. We started by fitting a 1st-degree polynomial to the data, using the polyfit tool in Python to perform this analysis. However, previous works reported a break in the correlation (see, for example, She et al. 2018), where a 1st-degree polynomial with a negative slope is found on one side and a positive slope on the other. Then, we used piecewise-regression⁵ tool in Python (Pilgrim 2021), to test the existence of the breakpoint. This package fits breakpoint positions and linear models for the different fit segments, and it gives confidence intervals for all the model estimates (see Muggeo 2003 for more details). We found a breakpoint at log(λ_Edd, break)= − 2.39 with σ = 0.45, in agreement with what was previously obtained by She et al. (2018), who proposed a break at −2.5. To determine which model (with a breakpoint or not) best represents the data, we use the Bayesian information criterion (BIC), where the model with the lowest BIC value is considered the best. With no break, BIC is −113.2, and with a break, it is −115.5. This suggests that the model with a breakpoint at log(λ_Edd, break) represents well the data, that is, AGN seems to follow the same relation as XRB. We also explored, the possibility of more breakpoints, and it does not change the previous result.

In Fig. 5, we also plot the relations given by other authors for comparison. For high-luminosity AGN (log(λ_Edd)>  − 2.39), we compare with Fanali et al. (2013), who studied a sample of 71 type 1 AGN using XMM-Newton data (purple dashed line). In the low-luminosity branch (log(λ_Edd)<  − 2.39), we compare with Gu & Cao (2009), which used a sample of 55 LLAGN using Chandra and XMM-Newton data (green dashed line); She et al. (2018) used a sample of 314 AGN with Chandra (cyan dashed line); and Younes et al. (2011) used Chandra and XMM-Newton data from a sample of 13 LINER with accretion rates below −4.5 (black dashed line).

In this work, we have shown that the inclusion of XMM-Newton + NuSTAR data and reflection models in the spectral fit improves the estimation of the spectral index – as also reported in Hinkle & Mushotzky (2021) – which could improve the scatter compared to what was previously found by Gu & Cao (2009), Younes et al. (2011), She et al. (2018). For details on the improvement of the uncertainties in the spectral index estimation, see Appendix C. Indeed, in Fig. 5 can be seen that our results, when compared with previous studies, seem to agree with the correlations found by Gu & Cao (2009), She et al. (2018), and Younes et al. (2011), but the effect of the large scatter in previous studies can be appreciated. The same is true for the high-accretion branch, where the relation of Fanali et al. (2013) (at log(λ_Edd)>  − 2.39) fits well the data of Esparza-Arredondo et al. (2021).

To determine whether there is a relation between Γ and λ_Edd, we use the tool pymccorrelation in Python (Isobe et al. 1986; Curran 2014; Privon et al. 2020) to test the relationship between two variables. This tool is able to calculate Pearson’s r, Spearman’s ρ, and Kendall’s τ correlation coefficients. In this work, we use the Kendall τ correlation test, a non-parametric method for measuring the degree of association of two variables in censored data (upper and lower limits) and taking into account the uncertainties in the parameters (see Akritas & Siebert 1996; Isobe et al. 1986 for a detailed explanation of the calculus). A Kendall’s τ close to zero indicates that there is no trend, and if they are perfectly related, Kendall’s τ becomes 1.0 (or −1.0 for an anti-correlation). For the LLAGN, log(λ_Edd)<  − 2.39, Kendall’s correlation coefficient is τ = −0.27. However, possibly because of the small number of sources, the associated p-value is 0.06, so the correlation is not formally significant and confirmation would require a larger sample. We note that our result is also consistent with a flat correlation for these objects. In the high-luminosity branch (log(λ_Edd)>  − 2.39), we obtain τ = 0.39 and a corresponding p-value of check 0.02, consistent with a positive correlation for high accreting sources. Thus, it appears that our sample provides evidence of a Γ − λ_Edd relation that is consistent with previous studies, although at lower statistical significance. In any case, the change in correlation between these parameters at log(λ_Edd)∼ − 2.39 highlights the change in accretion physics between high- and low-luminosity AGN, consistent with previous studies (Shemmer et al. 2006; Younes et al. 2011 and references therein).

Despite the small number of sources in our sample, we study the anti-correlation of the sample presented here using the tool linregress in Python. Then, for the low-luminosity branch where log(λ_Edd)<  − 2.39:

$$\begin{array}{c}\hfill \mathrm{\Gamma }=(-0.128\pm 0.272)\times log({\lambda }_{\mathrm{Edd}})+(1.348\pm 0.075).\end{array}$$

Our work allowed us to identify the change in correlation between the spectral index and the accretion rate at log(λ_Edd)∼ − 2.39, which is highly suggestive of a change in accretion physics in AGN. We recall that a larger sample of sources combining XMM–Newton and NuSTAR data and fitting physical reflection models would be very useful to confirm this relation.

6.3. Reflection

An important feature in the spectra of AGN is the reflection that imprints its mark at X-ray energies. The shape of this reflection component is characterized by the FeKα emission line and the Compton hump, peaking at ∼30 keV (Pounds et al. 1990). The gas producing the X-ray reflection in AGN could be related to the accretion disk, a neutral reflector such as the torus, or a combination of both emissions. Because we cannot separate these scenarios, in the following, we analyze the following scenarios, in which each of the structures dominates the X-ray spectra.

We started our analysis by studying the torus-like reflector. Previous studies have suggested a relation between the torus properties and the accretion rate. For example, Müller-Sánchez et al. (2013) used observations of VLT/SINFONI AO-assisted integral-field spectroscopy of H2 1-0 S(1) emission of four LLAGN, NGC 1052, NGC 2911, NGC 3169 and, NGC 1097, and found that on scales of 50–150 pc, the spatial distribution and kinematics of the molecular gas are consistent with a rotating thin disk, where the ratio of rotation (V) to dispersion (σ) exceeds the unity. However, in the central 50 pc in their sample, the observations reveal a geometrically and optically thick structure of molecular gas (V/σ <  1 and N_H >  10²³ cm⁻²). This can be associated with the outer extent of any smaller-scale obscuring structure. In contrast to Seyfert galaxies, the molecular gas in LLAGN has V/σ <  1 over an area that is about nine times smaller and with column densities that are on average ∼three times smaller. They interpret these results as evidence for a gradual disappearance of the nuclear-obscuring structure, also in agreement with what was previously found by González-Martín et al. (2017) using a sample of 109 AGN using IRS/Spitzer observations.

Later, Ricci et al. (2017b) found that the probability of a source being obscured in the X-rays (covering factor of gas) depends primarily on the Eddington ratio instead of on absolute luminosity. They propose that the radiation pressure on dusty gas is responsible for regulating the distribution of obscuring material around the central black hole. At high accretion rates, radiation pressure expels the obscuring material in the form of outflows (Fabian 2006). However, this work was made for the LOS column density, which is different from the torus column density (N_{H − LOS} ≠ N_{H − refl}). Here, we go on to analyze, the relation between the column density of the torus-like reflector and the Eddington ratio. We plot this relation in Fig. 6 and the values obtained in this work are presented in Table D.4. The pink circles and light yellow stars are the data points of the best fit model (borus02 in the indistinguishable cases) in this work and the ones obtained by Esparza-Arredondo et al. (2021), respectively. The blue points represent the binned data points for a bin size equal to 0.5 dex in λ_Edd.

Fig. 6. Relation between the column density of the torus-like reflector (in log) vs. the Eddington ratio, λEdd = Lbol/LEdd, for the sample of this work. The pink points and light yellow stars are the data point of the best fit model of this work (borus02) and the one obtained by Esparza-Arredondo et al. (2021). The blue points represented the binned data point for a bin size equal to 0.5. The black solid line represents the best fit and the light blue zone the 3σ confidence level.

Using Kendall’s tau correlation coefficient, we found a correlation coefficient of τ = 0.22 and p-value of 0.04 for a torus-like reflector log(N_H, refl) and λ_Edd, suggestive of a correlation – but confirmation is required using a larger sample. As the parameters seem to be positively correlated, we perform a linear regression of the data using polyfit in python, and we found the following relation:

$$\begin{array}{c}\hfill log(N_{\mathrm{H},\mathrm{refl}})=(0.126\pm 0.303)\times {\lambda }_{\mathrm{Edd}}+(24.166\pm 0.102).\end{array}$$

Therefore, we find that our data is consistent with the scenario where lower accretion rate objects have, on average, lower column density material in their surroundings. However, due to the size of the sample, our correlation shows a high dispersion, then it is also consistent with being flat at 3σ. We note that our torus fits allow for a free covering factor, so the lower column densities are not a consequence of a fixed covering factor in the model and a geometrically thinner reflector in lower accretion rate objects. For the LLAGN (log(λ_Edd)<  − 2.39) we obtain a mean value of the torus column density log(N_H, refl cm⁻²)=23.76, with σ = 0.74 and in the high-luminosity regime, log(N_H, refl cm⁻²)=24.09 with a standard deviation of σ = 0.56. Consequently, our result is in line with the results in the infrared, which suggests a gradual disappearance of the torus (Müller-Sánchez et al. 2013; González-Martín et al. 2017), and in agreement with the scenario proposed by Ricci et al. (2022), where it is expected that LLAGN has lower column density. These authors proposed a model in which AGN move in the obscuration–accretion rate plane during their life cycle. The growth of AGN begins with an unobscured AGN accreting at log(λ_Edd)≤ − 4. Then, an accretion event then takes place, in which the SMBH is fueled, and as a result the accretion rate, column density, and covering fraction all increase. As a consequence, obscured AGN are preferentially observed. When the Eddington limit for dusty gas is reached, the covering factor and the column density will decrease, leading to an unobscured AGN being typically observed. As the remaining fuel is depleted, the SMBH goes back into a quiescent phase (see Ricci et al. 2022 for more details). Even with small statistics, our results can be interpreted within the framework of this evolutionary model, in which radiation pressure regulates their evolution.

Then, we compare the column density of the reflector and the column density in the LOS. Zhao et al. (2021), using all AGN in the 100-month Palermo Swift/BAT catalog with LOS column density between 10²³ and 10²⁴ cm⁻² with available NuSTAR data shows that the average torus column density is similar for both Compton thin and CT-AGN, independent of the observing angle, with log(N_{H − refl} cm⁻²)∼24.15. In Fig. 7, we compare the column density of the torus and the absorption in the LOS of our work. The black dotted line represents the mean value of log(N_{H − refl} cm⁻²), previously found by Zhao et al. (2021), and the green zone the interval of log(N_{H − LOS}) of their work. We note that our data points and their fit are in agreement in the interval of log(L_{H − LOS}) of their work, namely, for the moderately obscured sources in our sample. The majority of galaxies with log(L_{H − LOS})< 23.0 in our sample are clearly below the value previously obtained, with a mean value of log(N_H−refl cm⁻²) ∼ 23.36 and σ = 0.59.

Fig. 7. Relation between the column density of the torus-like reflector (in log) vs. column density in the LOS (in log). The red points represent the data point of this work. The black dashed line corresponds to the value obtained by Zhao et al. (2021) and the black dotted represents x = y. The green zone is the interval of the column density in the LOS analyzed in their work.

In order to explore any correlation between these parameters, we calculated the Kendall τ correlation coefficient, and we found τ = −0.17 and p-value of 0.51, suggestive of a negative correlation but compatible with a lack of correlation as well between them. Therefore, more data points are necessary to establish any correlation between these parameters. The majority of the objects show a larger log(N_{H − refl}) than log(L_{H − LOS}), suggesting that the torus is not seen through its densest part, which is consistent with what was reported by Zhao et al. (2021).

Regarding the covering factor of the torus-like reflector, we obtained a mean value of CF=0.64 and σ = 0.26. We note that this parameter could be constrained for seven sources. For another seven sources, the best-fit value and an upper limit could be placed, and for an additional three, the best-fit value and lower limit could be placed. In addition, we analyze the correlation between this parameter and the accretion rate by the Kendall τ correlation test. We find a correlation coefficient τ = −0.06 and p-value of 0.73, suggesting that these parameters are not correlated.

Considering a disk-like reflector, we could constrain R_f only for one object (NGC 7674), while for the others we only obtain lower limits. Also, this model allows us to study the ionization degree of the disk, however, this parameter is also unconstrained, with only six galaxies having well-constrained values and one upper limit and five lower limits.

The results presented in this paper suggest that the distribution of the gas in the torus in AGN is a dynamic and very complex structure, showing changes in the physical properties of the torus linked to the luminosity of the AGN, in agreement with what was previously found in the literature in the X-rays and the infrared. Certainly, combining XMM-Newton + NuSTAR is key to exploring the structure and distribution of the reflector and constraining its physical and geometrical parameters, especially in the low-luminosity range.

7. Conclusions

In this work, we study the reflection of LLAGN by analyzing the broadband X-ray spectra of a BASS/DR2 sample with log(λ_Edd)<  − 3 (17 objects) using XMM-Newton+NuSTAR+Swift observations and characterizing the reflection features using the borus02 model to represent torus reflection and XILLVER to model accretion disk emission. The goal was to investigate the accretion mechanism by the relation between the spectral index and the accretion rate, as well as constraining the properties of the potential reflector. The main results are summarized below:

• 1.

  All objects in our sample are well-fitted with a torus-like reflector. Of these, eight objects are equally well-fitted with a torus and a disk (they are indistinguishable from a statistical point of view and visual inspection). These eight objects have consistent values for the spectral index Γ and luminosities when modeled with a torus or a disk reflector.

• 2.

  In our sample, we can classify six objects as unobscured (log(N_H)< 22), nine galaxies as obscured (22 <  log(N_H)< 24.18), and two as Compton-thick (using as a threshold N_H = 1.5 × 10²⁴ cm⁻², according to the torus model). According to the disk case, all the galaxies can be classified as Compton-thin.

• 3.

  Combining XMM-Newton+ NuSTAR and considering the reflection component in the spectral fitting, the uncertainties on the spectral index and the scatter in the relation between this parameter and the accretion rate are reduced when compared to previous works over similar ranges in accretion rate.

• 4.

  Our work is consistent with the negative slope found in previous works at log(λ_Edd)≤ − 3, and also consistent with the change in the Γ − log(λ_Edd) relation at log(λ_Edd)∼ − 3, where in sources with high accretion rates, the slope is known to be positive.

• 5.

  We found a tentative correlation between the torus properties (column density) and the accretion rate, suggesting that the torus has a decreasing column density with decreasing accretion rate. Consequently, AGN at log(λ_Edd)<  − 3 has a lower torus column density compared with more luminous AGN. This column density is derived from reflection as opposed to absorption in the LOS, so it is representative of the global column density of gas around the X-ray corona.

• 6.

  All AGN in our sample with a column density in the LOS, log(N_{H − LOS})< 23.0 have a torus with a column density higher than their log(N_{H − LOS}) and thereby the torus could be observed through an underdense region.

In the future, new X-ray missions such as HEX-P (Madsen et al. 2019)⁶ facilities will detect a large sample of LLAGN, which could help us further constrain the evolution of the AGN reflection and the accretion physics behind SMBH.

Acknowledgments

We thank the referee for the valuable comments that improved the manuscript. D.Y. acknowledges financial support from the Doctorate Fellowship program FIB-UV of the Universidad de Valparaíso and the Max Planck Society by a Max Planck partner group. LHG acknowledges funds by ANID – Millennium Science Initiative Program – ICN12_009 awarded to the Millennium Institute of Astrophysics (MAS). ELN acknowledges financial support from ANID Beca 21200718. CR acknowledges support from the Fondecyt Iniciacion grant 11190831 and ANID BASAL project FB210003. MB acknowledges support from the YCAA Prize Postdoctoral Fellowship. NOC acknowledges support from CONACyT. J.A.G. acknowledges support from NASA grant 80NSSC21K1567.

References

Appendix A: Notes and comparisons with previous results for individual objects

A.1. NGC 3998

This object is spectroscopically classified as a LINER (Heckman 1980; Ho et al. 1997), and no significant broad-line polarization (Barth et al. 1999). HST/WFPC2 optical images revealed an unobscured nucleus and the presence of a bright circumnuclear ionized gas disk (Pogge et al. 2000; Cazzoli et al. 2018).

In radio, it contains a nucleus (Hummel et al. 1984) displaying a weak jet-like northern structure (Filho et al. 2002) and according to Frank et al. (2016), the radio source in NGC 3998 shows two S-shaped lobes.

Pian et al. (2010) found that the spectra of NGC 3998 is better fitted with a simple absorbed power-law model, allowing the overall normalization of the model fit to each detector to be independent. Later, Kawamuro et al. (2016) found the same result using Suzaku and Swift/BAT and fitting the reflection with PEXRAV reflection model. They found that the reflection strength is very weak with an upper-limit of R_f <  0.10 and no significant iron Kα line emission is detected, either. They suggest that there is little surrounding matter around the nucleus.

Later, Younes et al. (2019), using XMM-newton and NuSTAR data found that it is an unabsorbed AGN and its spectrum in the energy range 0.5-60 keV using NuSTAR data is best fit with a power law and cutoff energy. They also studied the reflection hump and found only an upper limit for the reflection fraction, with value R_f <  0.3 at 3σ.

In the present work, the data were fitted with PEXMON and the result is consistent with a reflection fraction R_f <  0.15. Using borus02 we found that the data is consistent with a Compton thin torus (with log(N_H, refl)=22.50${}_{22.28}^{23.11}$) covering more than 20% of the sky. The inclination in this model is unconstrained. In case of XILLVER we found a model consistent with a disk with log ξ>3.99 and R_f>0.26. From a statistical point of view, borus02 and XILLVER models are indistinguishable. All the models are consistent with a column density in the LOS to classify this galaxy as a Compton-thin source with log(N_H, los)∼20.5.

A.2. NGC 3718

NGC 3718 was optically classified as a type 1.9 LINER (Ho et al. 1997). Cazzoli et al. (2018) showed that its broad Hα line is produced in the BLR rather than an outflow. The X-ray spectrum shows a small Fe Kα line, indicative of a reflection component (Younes et al. 2011; Hernández-García et al. 2014; Diaz et al. 2020).

Modeling the reflection with PEXMON we found R_f <  0.24. This low reflection signature can be produced by a torus-like neutral reflector (borus02), if it is Compton thin ($log(N_{H,refl})=22.{78}_{22.38}^{23.98}$), with an unconstrained covering factor (CF< 0.99${}_{0.05}^{1.00}$). Fitting a disk-like reflector instead, using XILLVER requires high levels of ionization with $log\xi =3.{12}_{1.35}^{3.61}$ and a R_f >  0.86. A (MEKAL) component was necessary to improve the fit at low energies in all the models. The XILLVER and borus02 are indistinguishable. All the models are consistent with a Compton-thin absorber in the LOS to the primary coronal emission with log(N_H, los)∼22.

A.3. NGC 4258*

NGC 4258* is a SABbc spiral galaxy, spectroscopically classified as a 1.9 Seyfert (Ho et al. 1997) and as LINER (Balmaverde & Capetti 2014). This source has a highly obscured central X-ray source and is well known for its anomalous arms, discovered on the basis of Hα imaging (Wilson et al. 2001) and for its nuclear H₂O megamaser, which traces a dense, edge-on disk on sub-parcsec scales (Watson & Wallin 1994; Greenhill et al. 1995). The nucleus of NGC 4258* also contains a relativistic radio jet (Doi et al. 2013).

Reynolds et al. (2009) combined Suzaku, XMM-Newton and Swift, and detected robust flux variability of the 6.4 keV iron line and suggested a model in which the line originates from the surface of a warped accretion disk. Also, during their Suzaku observation, they detected high amplitude intraday variability, with fluctuations on timescales as short as 5 ksec. Herrnstein et al. (2005) found that the absorption may well arise in the outer layers of the warped geometrically thin accretion disk, further reducing the need for any cold structure other than the accretion disk itself.

NuSTAR and XMM–Newton observations of this object are not simultaneous, so we explored possible spectral variability. We found variations in the normalization of the power law, so we permitted the normalizations of the power law to be free. In the next step, we fit the models. Using PEXMON we found R_f=0.78${}_{0.20}^{3.744}$. Using borus02 we found log(N_H, refl)=23.29${}_{22.77}^{24.39}$, so borderline Compton-thin/thick torus. The torus covers less than 81% of the sky. Fitting with XILLVER we found a disk with log ξ <  1.77 and R_f>3.80. In order to improve the fit in the soft energy band, it was necessary to add three components, modeled with two MEKAL plus a scattered power law in all the models. We note that in case of XILLVER, was necessary to use a cutoff power law in the scattered component as in this component, dominate above ∼70 keV. The borus02 model is the best representation of the data. All the models are consistent with a Compton-thin absorber in the LOS to the primary coronal emission with log(N_H, los)∼22.94.

A.4. ESO 253-G003

This source is spectroscopically catalogued as a Seyfert 2 in Véron-Cetty & Véron (2006). The XMM-Newton and NuSTAR observations of this object are separated by 3 days, so we fit the overlapping 3.0–10.0 keV range of both instruments with a simple absorbed power-law model and checked whether they were consistent. This test showed there was spectral variability between both epochs, with a significant change in the slope Γ.

Due to its complexity, the XMM-Newton data was excluded for our analysis. We performed a simultaneous fit of the NuSTAR and Swift/BAT data with PEXMON and significant reflection was detected with R_f>1.59. Modeling this reflection with borus02 instead of PEXMON results in a marginally Compton thin torus reflector (log(N_H, refl)=24.04${}_{23.93}^{24.30}$) covering more than 73% of the sky. Fitting the reflector using XILLVER instead of borus02 results in a highly ionized disk with $log\xi =2.{68}_{1.34}^{2.84}$ and R_f >  0.85. The PEXMON model shows the lowest ${\chi }_{\nu }^2=0.89$ although all models show similar quality fits (${\chi }_{\nu }^2$=0.905 for borus02 and ${\chi }_{\nu }^2$=0.903 for XILLVER). All the models are. indistinguishable, and they are consistent with a Compton-thin absorber in the LOS to the primary coronal emission with log(N_H, los)∼23.

A.5. NGC 1052

NGC 1052 is the brightest elliptical galaxy in the Cetus I group. It was optically classified as a LINER by Heckman (1980), then it was classified as LINER type 1.9 (Ho et al. 1997). Cazzoli et al. (2018) proposed that this object presents signs of outflowing winds by using optical 2D spectra.

This source is a radio-loud galaxy (Maoz 2007). At 1.4 GHz Very Large Array (VLA) image of NGC 1052 shows a core-dominated radio structure, with only about 15% of the flux density in extended emission: there are two lobes spanning 3 kpc, possibly with hot spots (Wrobel 1984).

At X-rays, Osorio-Clavijo et al. (2020) presented an extensive study of NGC 1052 using observations from Chandra, XMM-Newton, NuSTAR, and Suzaku. They reported variability in the nucleus and found variations both in the intrinsic continuum flux, photon index, and in the obscuration along the LOS. The reflection component is a steady emission both in flux and shape, fully consistent with reflection in a distant structure, perhaps the torus. They argue that NGC 1052 is in the regime of Compton-thin sources, consistent with the fact that the flux of the reflection component is not dominant in the hard band. In addition, Baloković et al. (2021) using NuSTAR, XMM-Newton, Suzaku and BeppoSAX observations and fitting the borus02 model, found a covering factor of ∼80-100% and a column density in the range (1-2)×10²³ cm⁻² that is well-matched with LOS column density.

In this work, we fit the data with PEXMON, and it is consistent with R_f=1.26${}_{0.79}^{3.05}$. Using borus02 model it has (log(N_H, refl)=24.35${}_{24.19}^{24.57}$), so a borderline Compton-thin or thick torus, and CF=0.50${}_{0.27}^{0.54}$. Using XILLVER we found log ξ=1.88${}_{1.78}^{2.00}$ and R_f>9.18. In case of borus02 it was necessary to add a MEKAL and a scattered power law to improve the fit. In case of XILLVER, it was necessary to add two components to improve the fit (an absorber and a MEKAL) and in case of PEXMON three components were necessary (MEKAL, power law and an absorber). From a visual inspection, we had to add two Gaussian emission lines centered at: S XIV at 2.4 keV and Si XIII at 1.85 keV. The borus02 model is the best representation of the data. All the models are consistent with a Compton-thin absorber in the LOS to the primary coronal emission with log(N_H, los)∼23.2.

A.6. NGC 2655

NGC 2655 is the brightest member of a group NBGG 12-10 in the Nearby Galaxy Groups catalog of Tully (1988). It is classified as Seyfert 2 in Swift/BAT 70 month catalog (Baumgartner et al. 2013) and as a Type-2 LINER according to the optical classification done by Ho et al. (1997), Véron-Cetty & Véron (2010).

In the X-rays, (González-Martín et al. 2009) found that the 0.2–10 keV XMM-Newton spectrum is well modeled by a MEKAL plus 2 power-law components. Kawamuro et al. (2016) fit 2 APEC models to the Suzaku spectrum of this galaxy and obtained a LOS $log(N_H)=2.{61}_{-0.17}^{+0.27}\times {10}^{23}$ cm⁻² and $\mathrm{\Gamma }=1.{77}_{-0.07}^{+0.19}$, which are consistent with earlier fits to ASCA data (Terashima & Wilson 2002). Kawamuro et al. (2016) also detected a Fe Kα line at 6.4 keV, which was not detected in the previous ASCA and XMM-Newton observations Terashima et al. (2002), González-Martín et al. (2009).

We fit the data with PEXMON placing an upper limit on R_f <  0.52. Modeling the reflector with a torus instead, using borus02, we found a column density of the reflector of log(N_H, refl)=23.24${}_{22.48}^{25.19}$ and CF>0.11. Modeling the reflection with an ionized disk instead of a torus, using XILLVER results in a highly ionized disk with a large reflection fraction log ξ <  3.55 and R_f=9.84${}_{\ast \ast }^{10}$. In the case of borus02 it was necessary to add two MEKAL and scattered power law to improve the fit. In case of XILLVER it was necessary to add a MEKAL and an absorber (modeled with zxipcf) components to improve the fit at low energies. In the case of PEXMON an additional MEKAL component was required. The three reflection models showed similar quality fits with ${\chi }_{\nu }^2=0.97$ for PEXMON, ${\chi }_{\nu }^2=1.01$ for both borus02 and XILLVER, they are indistinguishable.

All the models are consistent with a heavy but Compton-thin absorber in the LOS to the primary coronal emission with log(N_H, los)∼23.4. We note that if a torus is producing the reflection, then its average column density is lower than the absorber column density in the LOS.

A.7. NGC 3147*

NGC 3147* is an isolated Seyfert 2 galaxy (Ho et al. 1997) of Hubble morphological type SA(rs)bc. It was suggested as a true type 2 Seyfert candidate (Bianchi et al. 2012), but later confirmed as an AGN accreting at low rate (Bianchi et al. 2019). In the radio, it shows a point-like structure (Ulvestad & Ho 2001; Krips et al. 2006).

In the X-rays, Bianchi et al. (2017) using NuSTAR data, shows that spectrum of NGC 3147* can be simply modeled by a power-law with a standard Γ∼1.7 and an iron emission line. These spectral properties, together with significant variability on timescales as short as weeks, strongly support a LOS free of absorption for this source. They suggested that NuSTAR data adds further evidence in favor of an X-ray spectrum completely unaffected by absorption, confirming NGC 3147 as one of the best cases of true Type 2 Seyfert galaxies, intrinsically characterized by the absence of a BLR.

NuSTAR and XMM-Newton observations of these objects are not simultaneous, then we explored possible spectral variability. We found variations in the normalization of the power law, so this parameter was left free. We perform a simultaneous fit with this configuration of the data with PEXMON and found R_f=3.99${}_{0.97}^{6.41}$. We also fitted the data using borus02 model and found that the data is consistent with a Compton thin torus with density log(N_H, refl)=23.27${}_{22.98}^{23.94}$ covering less than 89% of the sky. Using the XILLVER model, the fit was consistent with log ξ=0.79${}_{0.06}^{1.15}$ and an R_f>0.91. In case of PEXMON, it was necessary to add one component to improve the fit (a warm absorber modeled with zxipcf). All the models are indistinguishable and consistent with a Compton-thin absorber in the LOS to the primary coronal emission with log(N_H, los)< 20.9.

A.8. NGC 2110*

NGC 2110* is a nearby S0 galaxy, (de Vaucouleurs 1991) Seyfert 2 AGN. NGC 2110* shows a prominent Fe Kα line accompanied by variable intrinsic emission (Hayashi et al. 1996).

Marinucci et al. (2014) reported the X-ray spectral analysis of NGC 2110* observed by NuSTAR in 2012, when the source was at the highest flux level ever observed, and in 2013, when the source had more typical flux levels. They found an upper limit on R_f < 0.14, confirming results from past high-energy BeppoSAX and Suzaku observations (Malaguti et al. 1999; Reeves et al. 2006; Rivers et al. 2014). Using MyTorus (Murphy & Yaqoob 2009) to model the reflection, they found a CF of 0.5 with equatorial N_H = 2.0 ± 1.1 × 10²³ cm⁻².

We re-examined the 2013 NuSTAR observations together with the 2003 XMM-Newton observation in the overlapping 3–10 keV band, finding a consistent slope and variable normalization. We therefore fit jointly the broad band spectra with untied normalization for the cutoff power law component between instruments. Modelling the reflection with PEXMON gives R_f=1.01${}_{0.91}^{1.12}$. Modeling the reflection with a torus instead, using borus02, we found a Compton thin reflector with log(N_H, refl)=23.49${}_{23.03}^{23.76}$ covering 13% of the sky, consistent with Baloković et al. (2018). Modeling the reflection with a disk instead of the torus, with XILLVER, we obtain log ξ=1.99${}_{1.82}^{2.01}$ and R_f>7.80. In the case of XILLVER it was necessary to add two components to improve the fit at low energies: a MEKAL and zxipcf. In case of borus02 three components were necessary to improve the fit: a scattered power law, MEKAL, and a zxipcf, and in the other case, four components were necessary: 2 zxipcf + 2*MEKAL. The torus reflector produces the best fit.

All the models are consistent with a Compton-thin absorber in the LOS to the primary coronal emission with log(N_H, los)∼22.01.

A.9. LEDA 96373*

LEDA 96373 (also 2MASX J07262635-3554214 or IGR J07264-3553) was first observed at high energies during the all-sky hard X-ray IBIS survey (Krivonos et al. 2007) and reported in the Palermo Swift/BAT survey. This galaxy is classified in NED as a Seyfert 2 and as a Compton thick galaxy (Koss et al. 2016).

Landi et al. (2009) observed this galaxy with Swift XRT data and reports an excess emission below 2 keV and found that a double power law model is a good fit. They found an intrinsic column density of ∼7 × 10²² cm⁻² and a photon index Γ∼2.5. The source seems to vary by a factor of 2 within a timescale of few days, with an average of 2-10 keV flux of ∼4.2 × 10⁻¹³ erg cm⁻² s⁻¹.

Simultaneous NuSTAR and XMM-Newton data are not available for this source, so we explored possible variability. We found variations in the normalization of the power law, so this parameter was left free. Using PEXMON we found R_f>9.38. Using borus02 we found that the data is consistent with a Compton-thin torus (log(N_H, refl)=23.40${}_{23.33}^{23.69}$) covering 59% of the sky. Then we fitted XILLVER and found log ξ=2.08${}_{2.03}^{2.17}$ and R_f>9.37. In case of PEXMON and borus02, it was necessary to add three components to improve the fit (two MEKAL*zxipcf), and in case of XILLVER, (MEKAL + a scatter power law)*zxipcf components were necessary. The borus02 model is the best representation of the data. The fit with borus02 requires Compton-thick absorption in the LOS to the coronal emission with log(N_H, los)= 24.09${}_{-24.06}^{+24.17}$.

A.10. NGC 2992

NGC 2992 is a spiral galaxy classified as a Seyfert 1.9 in the optical. This Seyfert galaxy is a changing-look AGN that varies from type 2 to intermediate type sometimes accompanied by extreme X-ray activity and back, over the span of a few years, (Shu et al. 2010) which is attributed to intrinsic variations of the power-law flux and to changes in the absorption in the LOS (Parker et al. 2015; Hernández-García et al. 2017).

At higher X-ray energies, a combined study with INTEGRAL, Swift and BeppoSAX data published by Beckmann et al. (2007) showed that variations in the normalization of the power law were needed when using an absorbed broken power-law model to fit the data simultaneously. They found a constant Γ and flux variations by a factor of 11 in timescales of months to years.

We used only the XMM-Newton data that overlaps in time with the NuSTAR exposure and perform a simultaneous fit. Modeling the reflection with PEXMON results in R_f=0.37${}_{0.33}^{0.39}$. Using borus02 model for the reflection instead of PEXMON we found a Compton-thin torus reflector with log(N_H, refl)=23.17${}_{23.13}^{23.29}$ covering more than 75% of the sky. Then we fit XILLVER instead of borus02 and the fit was consistent with a disk with log ξ <  0.08 producing R_f>8.80. In the case of PEXMON and borus02, it was necessary to add a MEKAL, a scattered power law and zxicpf components; in the case of XILLVER, two MEKAL components were required. The borus02 model is the best representation of the data.

All the models are consistent with a Compton-thin absorber in the LOS to the primary coronal emission with log(N_H, los)∼21.9.

A.11. M 51

M 51 (also known as Messier 51a, M 51a, and NGC 5194) is a very nearby (7.1 Mpc; Takáts & Vinkó 2006) Compton-thick AGN. It hosts a Seyfert 2 nucleus (Ho et al. 1997; Dumas et al. 2011) and shows two radio lobes that are filled with hot X-ray gas (Terashima et al. 2000) and an outflow of ionized gas (Bradley et al. 2004). This is an interacting spiral galaxy and lies in the constellation Canes Venatici.

Brightman et al. (2018) reported a ULX located close enough to the nucleus of M 51 that NuSTAR could not resolve it. They modeled this emission with a cut-off power-law model, where $\mathrm{\Gamma }=-2.{21}_{-2.16}^{+2.26}$ and $E_{\mathit{cut}}=1.2_{-1.1}^{+1.4}$ keV. We included the emission from the ULX using the same model setting Γ and E_cut frozen to the parameters already obtained by them. Here, we perform a simultaneous fit of the XMM-Newton, NuSTAR and Swift/BAT data. Using PEXMON reflection model we find R_f >  7.86. Using borus02 we found a Compton thick torus reflector (log(N_H, refl)> 24.68) covering less than 57% of the sky. Then we fit XILLVER and the fit was consistent with an accretion disk with log ξ <  1.26 and with R_f >  8.96. In the cases of PEXMON and XILLVER, one additional component was needed (modeled with zxipcf). The torus model gives a good representation of the data. The fit with this model requires Compton-thick absorption in the LOS to the coronal emission with log(N_H, los)> 24.45 otherwise with the other models, a Compton-thin column density is a good representation of the data.

A.12. NGC 2273*

NGC 2273* is a Seyfert 2 galaxy (Risaliti et al. 1999) considered a Compton-thick due to the detection of a strong Fe Kα line with an equivalent width > 1 keV and a low ratio of L_X/L[OIII] Guainazzi et al. (2005). Observations with Suzaku suggested that the nucleus of NGC 2273 is obscured by a Compton-thick column of 1.5 × 10²⁴ cm⁻² (Awaki et al. 2009). A hidden nucleus is also suggested by the detection of polarized broad lines (Moran et al. 2000).

We explore the spectral variability in this galaxy between the 2003 XMM-Newton observation and the 2014 NuSTAR observation by fitting an absorbed power law plus 6.4 keV Gaussian model to the overlapping 3–10 keV spectra. we find that allowing only the power-law normalization to vary between instruments or epochs produced an acceptable fit with ${\chi }_{\nu }^2=1.08$ and with consistent values of the normalization, so we will fit these spectra with all the parameters tied. We perform a simultaneous fit of the XMM-Newton, NuSTAR and Swift/BAT data with PEXMON reflection model and the result is consistent with R_f>9.02. We modeled the reflection using borus02 instead of PEXMON finding a Compton thin torus ($log(N_{H,refl})=23.{75}_{23.56}^{24.19}$) with a covering factor lower than 59% of the sky. Modeling the reflection withXILLVER instead produced moderately ionized disk with $log\xi =2.{15}_{1.74}^{2.31}$ and R_f>8.80. In all the models, it was necessary to add two components (MEKAL and absorber). The disk-like reflector modeled with XILLVER is the best representation of the data. The fit with this model requires Compton-thick absorption in the LOS to the coronal emission with log(N_H, los)> 24.55.

A.13. HE 1136-2304

HE 1136-2304 changed its optical spectral classification from 1994 (Seyfert 2) to 2014 (Seyfert 1.5) and can be considered an optical changing look AGN (Zetzl et al. 2018). HE 1136-2304 has been detected as a variable X-ray source by the XMM-Newton slew survey in 2014 (Parker et al. 2016). The 0.2–2 keV flux increased by a factor of about 30 in comparison to the 1990 ROSAT All-Sky Survey flux (RASS; Voges et al. (2000)). However, no clear evidence of X-ray absorption variability has been seen.

Parker et al. (2016) found an absorbing column density on the X-ray spectrum around 10²¹ cm⁻² in addition to Galactic absorption. This spectrum shows a moderate soft excess and a narrow Fe line at 6.4 keV and a high energy cut-off at ∼ 100 keV.

Although the X-ray flux is highly variable in this source, a joint spectral fitting with XMM-Newton and NuSTAR could be made since the observations were simultaneous. Using PEXMON reflection model we found R_f=3.99${}_{2.85}^{5.10}$. We also fit the data using borus02 model instead of PEXMON, and we found that reflection features can be produced by a Compton thin torus with log(N_H, refl)=23.52${}_{23.38}^{23.64}$ covering 59% of the sky. Modeling the reflection with a disk instead, using XILLVER results in a highly ionized disk with log ξ=3.78${}_{3.77}^{3.79}$ and R_f>1.45. To improve the fit in the low energy range, it was necessary to add one MEKAL component in all the models. The borus02 model is the best representation of the data. All models required a Compton-thin absorber on the LOS to the primary coronal emission, with log(N_H, los)∼20.97.

A.14. IGRJ 11366-6002:

We did not find previous information about this galaxy in the literature.

Here, we perform the fit of the NuSTAR and Swift/BAT and Swift/XRT data. When using PEXMON reflection model, we find R_f=0.72${}_{0.49}^{1.57}$. We also fit the data using borus02 model, and we found that the data is consistent with a Compton thick reflector torus (log(N_H, refl)>24.47) covering ∼93% of the sky. Using XILLVER, the fit was consistent with log ξ=3.14${}_{3.03}^{3.32}$ and R_f>1.12. All the models are indistinguishable and they are consistent with a Compton-thin absorber in the LOS to the primary coronal emission with log(N_H, los)< 21.2.

A.15. IC 4518A

This galaxy was optically classified as a type 2 Seyfert galaxy (Zaw et al. 2009) and it is classified as a Compton-thin source (Bassani et al. 1999; de Rosa et al. 2008; Hernández-García et al. 2015).

de Rosa et al. (2008) presented a 0.2–200 keV broad-band study of this galaxy with INTEGRAL, XMM-Newton, Chandra and ASCA to investigate the continuum shape and the absorbing and reflecting medium properties. They fitted PEXRAV the reflection model and found that in this object the presence of the reflection component above 10 keV is statistically required by the data. However, they found the best fit value larger than 1 that they suggested that could be related with the geometry of the reflector, which should be more complex than that used in their work.

Hernández-García et al. (2015) studied XMM-Newton data of this galaxy and found variations in an eight-day period, that correspond to a flux variation of 40% (41%) in the soft (0.5 - 2 keV) (hard - 5 - 10.0 keV) energy band.

For the purposes of this work, NuSTAR and XMM-Newton data were not simultaneous, thus we explored possible spectral variability. The best fit resulted when the normalization of the power law is a free parameter. Using PEXMON we found R_f=3.41${}_{1.10}^{6.17}$. Using borus02 we find a Compton-thick reflector torus with log(N_H, refl)=24.59${}_{24.11}^{25.16}$, covering more than 39% of the sky. Using XILLVER we found a model consistent with log ξ <  1.73 and R_f>2.94. To improve the fit, three additional components are necessary (MEKAL + scatter power law)*zxipcf for PEXMON, MEKAL + zxipcf for borus02) and XILLVER). The borus02 model is the best representation of the data. In all the models, the absorption in the LOS is consistent with being Compton-thin (log(N_H, los)∼ 23.22).

A.16. NGC 7674*

NGC 7674* (Mrk 533) is a nearby luminous infrared galaxy (LIRG). This galaxy is the brightest member of the Hickson 96 interacting galaxy groups and it is a known as a Seyfert 2 galaxy with broad Hα and Hβ components in polarized light (Young et al. 1996). Later, Bianchi et al. (2005) classified it as a changing look AGN in the X-ray range, switching between Compton-thin and Compton thick absorption in the LOS.

In the X-rays, NGC 7674* was first reported to be a reflection-dominated AGN by Malaguti et al. (1998) from BeppoSAX X-ray observations carried out in 1996, with the direct (intrinsic) continuum being fully absorbed by a Compton-thick gas column. Bianchi et al. (2005) studied the XMM-Newton spectra of this galaxy, finding a reflection fraction R_f∼1.5 and arguing that the changes on the CT regimes in this galaxy is associated with material in the LOS. Later, Gandhi et al. (2017) presented NuSTAR, Suzaku and Swift reporting a flat X-ray spectrum, suggesting that it is obscured by Compton-thick gas. Based upon long-term flux dimming, previous work suggested the alternate possibility that the source is a recently switched-off AGN with the observed X-rays being the lagged echo from the torus. Their data show the source to be reflection-dominated in hard X-rays, but with a relatively weak neutral Fe Kα emission line and a strong Fe XXVI ionised line. Also, they construct an updated long term X-ray light curve of NGC 7674* and find that the observed 2–10 keV flux has remained constant for the past 20 years.

As this galaxy is a changing look galaxy and the XMM-Newton and NuSTAR observations are separated by several years, we do not attempt a joint fit and exclude the XMM-Newton observation for our analysis. Using PEXMON we found that the data is consistent with R_f>2.21. Using borus02 to model the reflection instead of PEXMON we found a torus with log(N_H, refl)>24.43 covering 70% of the sky. Modeling the reflection with a disk instead of a torus, using XILLVER the fit was consistent with log ξ=2.82${}_{2.77}^{3.02}$ and R_f=3.25${}_{2.26}^{6.07}$. The XILLVER model is the best representation of the data.

We find heavy but Compton thin absorption in the LOS to the primary coronal emission with log(N_H, los)∼23 for all the models used.

A.17. NGC 5033

NGC 5033 is a nearby spiral galaxy with a low-luminosity Seyfert 1.8 type nucleus (Véron-Cetty & Véron 2010) with a point-like central X-ray source (Terashima et al. 1999). This galaxy has been alternatively classified as a Seyfert of type 1.5 (Ho et al. 1997). This object has a giant neighbor within a short distance: NGC 5005, a Seyfert SAB(rs) galaxy.

In the radio, it is predominated by a compact core (Ho & Ulvestad 2001), and showing extended jet-like features to the East-West (Pérez-Torres & Alberdi 2007).

In the X-ray, Terashima et al. (1999) using ASCA observations found a point-like X-ray source in the 2-10 a keV band. Their X-ray light curve showed variability on a timescale of ∼ 10⁴ s with an amplitude of ∼20%.

In this work, we combined the NuSTAR observations and perform a simultaneous fit of the XMM-Newton, NuSTAR, and Swift/BAT data. Using PEXMON reflection model, the results is R_f>1.0. Using borus02, we found that the data are consistent with a Compton thin torus (log(N_H, refl)=23.49${}_{23.34}^{23.69}$) covering 77% of the sky. Then we fit XILLVER instead of borus02 and the fit was consistent with log ξ <  1.0 and R_f>4.88. To improve the fit, it was necessary to two MEKAL to XILLVER. The borus02 model is the best representation of the data. In all the models, the required LOS absorption is compatible with being Compton-thin (log(N_H, los)< 21.0).

Appendix B: Existence of the reflection component

The aim of our spectral analysis is to quantify how much reflecting material can exist around the central engine, assuming that there must be material around the SMBH. In this section, we fit the nuclear continuum with and without reflection to probe that this component is usually required by the data. We based our analysis on the statistical significance of thereflection component compared to the model without reflection (i.e. F-test< 10⁻³), we examine the existence of the reflection component.

We start the analysis by fitting the spectra to the following model: constant*phabs*zphabs*cutoffpl. This model represents a power-law component (cutoffpl - associated with the intrinsic continuum) absorbed by the material along the LOS to the observer (phabs). The cross-normalization constant and the galactic absorption are denoted by constant and zphabs respectively. The column density (N_H), photon index (Γ), high energy cut off (E_cut), and normalization of the power law are the free parameters for this model. In order to avoid including the soft emission in this section of the analysis, we omit the spectra below 3.0 keV.

Then, we used constant*phabs*zphabs*(cutoffpl+pexmon), to study the case considering the reflection. We will use the pexmon model to fit the data. This model assumes Compton reflection from neutral X-ray photons in an optically thick material with plane-parallel geometry. The photon index, high energy cut-off, reflection fraction (R_f), metal and iron abundances, inclination angle, and normalization are all free parameters in pexmon. We have set the photon index of the pexmon to be the same as that of the power law. The iron abundance was assumed to be solar. The high energy cut off and normalization are free parameters. To take into account the reflection component and exclude the intrinsic power-law continuum in pexmon model, we have set R_f to -1. Then, we fit the spectra with this model and test if reflection is required by the data using the F-test tool within Xspec.

We present the results of this analysis in Table B.1, the addition of a reflection component improves the spectral fit for 14 out of the 17 sources (82% of the sample). Statistically, the reflection component is not required in three galaxies in our sample. NGC 2655 which has the lowest number of counts, and NGC 3998 and NGC 3718 which have the lowest Eddington ratios in our sample. Due to their luminosity, the contribution of the reflection could be very small, and as a consequence we are not able to quantify the reflection component. Because most of the sources do require a reflection component, we will use reflection models to fit all the sources in order to put constraints on this component.

Appendix C: Effect of the inclusion of the NuSTAR data and reflection models in the spectral index estimation

In this work, we use NuSTAR+XMM-Newton+Swift observations of a sample of LLAGN from BASS /DR2. We fit all our data with reflection models such as borus02 to model a torus-like reflector and XILLVER to model a disk-like emission. This may affect the best-fit parameters we obtain compared to results in the literature that use data from other instruments and with other reflection models.

In previous studies, the Γ parameter has large error bars (see, for example, She et al. 2018). To investigate the improvement in the spectral index estimation, we will compare the error bars of the present work using PEXMON and NuSTAR+XMM-Newton+Swift observations, with those obtained by Ricci et al. (2017a) using PEXRAV of the same sample with Swift/XRT, Swift/ BAT, ASCA, Chandra, and Suzaku observations. In Fig. C.1, we compare these error bars, with the blue stars representing the lower limits and the red circles the upper limits of the parameter. The black dotted line represents x=y. We found that the errors in this work are smaller compared to Ricci et al. (2017a). Then, the uncertainties in the photon index improve by including NuSTAR data and/or using models such as borus02 or XILLVER.

Fig. C.1. Relation between the error between our work using NuSTAR data and PEXMON model and the error values obtained using PEXRAV and observation from Swift/XRT, Swift/BAT, ASCA, Chandra, and Suzaku by Ricci et al. (2017a). The red circles represent the upper limit of the error bar and the blue stars the lower limit. The dotted black line is x=y.

Moreover, the Γ-log(λ_Edd) relation shows a high scatter (especially in the case of LLAGN), which is still not understood - it could be due to the sensitivity of the measurements or the intrinsic diversity of the nuclei (see Gu & Cao 2009; Younes et al. 2011; Yang et al. 2015; She et al. 2018). In this work, we found a significant decrease in the scatter in this correlation. To understand this improvement, the natural question is whether this is related to the inclusion of NuSTAR data or could be an effect of using more physical models than borus02 to perform the spectral analysis. In the upper panel in Fig. C.2, the spectral index is plotted against the accretion rate previously determined by Ricci et al. (2017a) using PEXMON. In the middle, we plot the relation between the same sample and the values determined in this work with PEXMON including the NuSTAR data. We note that we have not indicated the error bars in the accretion rate for illustrative purposes only. There is an improvement in the scatter and error bars, with smaller uncertainties in the photon index than in the results without including the NuSTAR data. Thus, the inclusion of the NuSTAR data is crucial for determining one of the fundamental parameters describing the X-ray emission, such as the spectral index Γ. We also investigate the improvement of the estimate of the parameter by including the borus02 model. Comparing the lower panel in Fig. C.2 with the middle panel in the same figure, we see a significant improvement in the scatter of the relation and the uncertainties of the parameter have improved significantly.

We found that we get lower uncertainties with a more physical model such as borus02Ṫhe process of combining NuSTAR data with more physical models such as borus02 is key to improving the scatter previously observed in this relationship. With our small error bars and improved constraints on the parameters, we can continue the search for the expected correlations in more detail.

Fig. C.2. Relation between the spectral index, Γ, and the accretion rate, log(λEdd). The blue points (top panel) represent the values obtained by (Ricci et al. 2017a) using PEXRAV model from the DR1. The middle panel (red points) represents the data points obtained in this work using PEXMON reflection model. In the lower panel (green points) are the values obtained in this work using borus02 reflection model. The dotted black lines are the best fit in all the panels.

Appendix D: Tables

Appendix E: Spectral models

Fig. E.1. Spectral modeling of NGC 3998. Plots correspond to borus02 (left) and XILLVER (right).

Fig. E.2. Spectral modeling of NGC 3718, NGC 4258*, ESO 253-G003* and NGC 1052. Plots correspond to borus02 (left) and XILLVER (right).

Fig. E.3. Spectral modeling of NGC 2655, NGC 3147, NGC 2110 and LEDA 96373. Plots correspond to borus02 (left) and XILLVER (right).

Fig. E.4. Spectral modeling of NGC 2992, M 51, NGC 2273 and HE 1136-2304. Plots correspond to borus02 (left) and XILLVER (right).

Fig. E.5. Spectral modeling of IGRJ 11366-6002, IC4518, NGC 7674 and NGC 5033. Plots correspond to borus02 (left) and XILLVER (right).

All Tables

All Figures

	Fig. 1. Schematic view of the methodology followed to fit the data. Note that the loop in blue iterate a maximum of four times. For a detailed explanation of the method, we refer to the main text.
In the text

	Fig. 2. Comparison between the spectral index estimation, Γ, between the models. The dotted lines represent the mean values.
In the text

	Fig. 3. Histogram of the intrinsic luminosity in the band 2.0–10.0 keV (left) and 10.0–79.0 (right) in the borus02 and XILLVER with number of objects in each group.
In the text

	Fig. 4. Intrinsic luminosity in the 2.0–10.0 keV range from BASS/DR2 and our work. Dotted black line represents x = y.
In the text

Fig. 5. Correlation between the spectral index, Γ, from individual fits, vs. the Eddington ratio, log(λEdd)=log(LBol/LEdd), for our sample of galaxies of the best fit models. The dotted and dashed green line represent the relation given by Gu & Cao (2009), the orange dotted represents Younes et al. (2011), the magenta dashed line is the relation obtained by She et al. (2018), while the solid black line is the correlation obtained in this work. The purple dashed line corresponds to the relation found by Fanali et al. (2013). The blue points represented the binned data. The pink points and light yellow stars are the data point of the best fit model in this work and the ones obtained by Esparza-Arredondo et al. (2021).

In the text

Fig. 6. Relation between the column density of the torus-like reflector (in log) vs. the Eddington ratio, λEdd = Lbol/LEdd, for the sample of this work. The pink points and light yellow stars are the data point of the best fit model of this work (borus02) and the one obtained by Esparza-Arredondo et al. (2021). The blue points represented the binned data point for a bin size equal to 0.5. The black solid line represents the best fit and the light blue zone the 3σ confidence level.

In the text

	Fig. 7. Relation between the column density of the torus-like reflector (in log) vs. column density in the LOS (in log). The red points represent the data point of this work. The black dashed line corresponds to the value obtained by Zhao et al. (2021) and the black dotted represents x = y. The green zone is the interval of the column density in the LOS analyzed in their work.
In the text

	Fig. C.1. Relation between the error between our work using NuSTAR data and PEXMON model and the error values obtained using PEXRAV and observation from Swift/XRT, Swift/BAT, ASCA, Chandra, and Suzaku by Ricci et al. (2017a). The red circles represent the upper limit of the error bar and the blue stars the lower limit. The dotted black line is x=y.
In the text

Fig. C.2. Relation between the spectral index, Γ, and the accretion rate, log(λEdd). The blue points (top panel) represent the values obtained by (Ricci et al. 2017a) using PEXRAV model from the DR1. The middle panel (red points) represents the data points obtained in this work using PEXMON reflection model. In the lower panel (green points) are the values obtained in this work using borus02 reflection model. The dotted black lines are the best fit in all the panels.

In the text

	Fig. E.1. Spectral modeling of NGC 3998. Plots correspond to borus02 (left) and XILLVER (right).
In the text

	Fig. E.2. Spectral modeling of NGC 3718, NGC 4258*, ESO 253-G003* and NGC 1052. Plots correspond to borus02 (left) and XILLVER (right).
In the text

	Fig. E.3. Spectral modeling of NGC 2655, NGC 3147, NGC 2110 and LEDA 96373. Plots correspond to borus02 (left) and XILLVER (right).
In the text

	Fig. E.4. Spectral modeling of NGC 2992, M 51, NGC 2273 and HE 1136-2304. Plots correspond to borus02 (left) and XILLVER (right).
In the text

	Fig. E.5. Spectral modeling of IGRJ 11366-6002, IC4518, NGC 7674 and NGC 5033. Plots correspond to borus02 (left) and XILLVER (right).
In the text