© The Authors 2023

1. Introduction

Active galactic nuclei (AGN) are powered by accreting supermassive black holes (SMBHs), surrounded by a torus of obscuring material. According to the unification theory (Urry & Padovani 1995), this torus is uniform and obscures certain lines of sight, preventing us from observing the broad line region (BLR, composed of gas clouds closely orbiting the black hole) from certain lines of sight. However, more recent studies based on infrared (IR) spectral energy distributions (SEDs) favor a scenario in which this torus is clumpy or patchy, rather than uniform (e.g., Nenkova et al. 2002; Ramos Almeida et al. 2014). This has been further confirmed by direct observations of changes in the line-of-sight (l.o.s.) obscuration (N_H, los) in the X-ray spectra of nearby AGN (e.g., Risaliti et al. 2002).

Obscuration variability in X-rays has been detected in a large range of timescales, from ≲1 day (e.g., Elvis et al. 2004; Risaliti et al. 2009) to years (e.g., Markowitz et al. 2014). Similarly, a large range of obscuring density variations have been observed: from small variations of Δ(N_H, los)∼10²² cm⁻² (e.g., Laha et al. 2020) to the so-called changing-look AGN, which transition between Compton-thin (N_H, los < 10²⁴ cm⁻²) and Compton-thick (N_H, los > 10²⁴ cm⁻²) states (e.g., Risaliti et al. 2005; Bianchi et al. 2009; Rivers et al. 2015).

Despite the multiple works that detect a Δ(N_H, los) between two different observations of the same source, very few have observations covering a complete eclipsing event (e.g., Maiolino et al. 2010; Markowitz et al. 2014). This is because observing the ingress and egress of single clouds across the line of sight may require daily observations over years. In fact, the most extensive statistical study of N_H, los variability to date is the result of frequent monitoring of 55 sources (Markowitz et al. 2014), spanning a total of 230 years of equivalent observing time with the Rossi X-ray Timing Explorer (RXTE, Jahoda et al. 2006). And it resulted in the detection of variability in only five Seyfert 1 (Sy1) and three Seyfert 2 (Sy2) galaxies, with a total of eight and four eclipsing events respectively. This study has been used to calibrate the most recent X-ray emission models based on clumpy tori (e.g., Buchner et al. 2019).

While it is clear that further studies such as the one mentioned are not possible with the current X-ray telescopes, due to time constraints of pointed observations, studies including large samples of sources with sporadic observations can still be particularly helpful in understanding the torus structure. The Δ(N_H, los) between two different observations, separated by a given Δt, has been used to set upper limits to cloud sizes and/or their distances to the SMBH (e.g., Risaliti et al. 2002, 2005; Pizzetti et al. 2022; Marchesi et al. 2022).

Recently, Laha et al. (2020) studied the variability of 20 Sy2s and found that only seven out of the 20 sources showed changes in N_H, los over timescales from months to years. A particularly interesting source also showed an increase in N_H, los over a period of 3.5 yr, and then remained seemingly constant for ∼11 yr. Laha et al. (2020) further argued that obscured AGN in which N_H, los variability is not present, or is only present on yearly timescales, are difficult to reconcile with a simple clumpy torus scenario. Similarly, Hernández-García et al. (2015) found long-term N_H-variability in only 11 out of 25 Sy2 galaxies, and no short-term variability in ten sources analyzed. The presence of a two-phase medium (e.g., Siebenmorgen et al. 2015), or important contributions of larger-scale structures in the galaxy (e.g., gas lines or filaments) have been suggested as possible alternatives to obscuration in such cases.

Even now, the number of well-studied sources in the literature still remains small. In particular, very few works exist dedicated to analyzing larger samples of AGN with multiepoch X-ray observations. Even in such cases, they tend to use phenomenological models (e.g., Markowitz et al. 2014; Laha et al. 2020), which do not allow for a comparison between the N_H, los variability and general torus properties.

Recently, a variety of self-consistent physical torus models aiming to better fit the reflection component of AGN X-ray spectra have been developed. Some are based on a uniform torus assumption, such as MYTorus (Murphy & Yaqoob 2009) or borus02 (Baloković et al. 2018), and have been widely tested. Others, while more recent and perhaps not as robustly tested, include the option of a clumpy or patchy torus, such as UXCLUMPY (Buchner et al. 2019) and XCLUMPY (Tanimoto et al. 2019). All of these models, both uniform and patchy, take advantage of the high-energy coverage of telescopes such as the Nuclear Spectroscopic Telescope Array (hereafter NuSTAR, Harrison et al. 2013) to accurately model the reprocessed emission of the torus (i.e., the reflection component). Through this process, quantities such as the torus covering factor, the inclination angle, and the average torus column density can be estimated.

In this work, we aim to analyze a sample of 12 likely variable AGN that have multiple X-ray observations, covering timescales of weeks to decades. These have been selected from a parent sample of BAT-detected, Compton-thin AGN at low (z < 0.1) redshift, which have archival NuSTAR observations. We used three different physical torus models, with the objective of comparing our results on N_H, los variability to various torus properties.

The sample selection and data reduction processes are discussed in Sects. 2 and 3. In Sect. 4 we discuss the physical torus models used to model the spectra of the sources, and the various torus properties that can be derived from each of them. In Sect. 5 we discuss the methods we used to classify a source as N_H, los-variable, or non-variable. And finally, our results and a discussion on those are provided in Sects. 6 and 7, respectively. We add a conclusion in Sect. 8. Further information, such as tables listing fit parameters, images of the spectra, and comments on individual sources can be found in Appendices A, B, and C, respectively.

2. Sample selection

The sample in this work was selected from Zhao et al. (2021), a work performing a broadband X-ray spectral analysis of an unbiased sample of 93 heavily obscured AGN (with line-of-sight column density 23 ≤ log(N_H)≤24; i.e. Compton-thin AGN) in the nearby Universe, for which high-quality archival NuSTAR data were available. This sample, derived from the Swift-BAT catalog (Burst Alert Telescope, observing in the 15–150 KeV range, Oh et al. 2018) is the largest NuSTAR dataset analyzed to date. Zhao et al. (2021) estimated torus geometry and N_H, los for the whole sample by jointly fitting a NuSTAR observation and a non-simultaneous soft X-ray observation, from either XMM-Newton, Chandra, or Swift-XRT.

It is an ideal starting sample, first of all because a BAT detection already guarantees that the sources are X-ray bright and are typically at low redshift (z < 0.12). Secondly, all sources analyzed already have one NuSTAR observation, which is essential in breaking the degeneracy between reflection and line-of-sight components, allowing us to constrain torus parameters. On top of that, it is a sample of Compton-thin AGN. These are obscured enough to let the reflection component shine through, allowing us to study the torus geometry, while being unobscured enough to allow us to constrain N_H, los with low uncertainty (compared to e.g., Compton-thick AGN).

Through a preliminary study performed in their analysis of the sample, Zhao et al. (2021) found that at least 31¹ of the sources presented variability (either in N_H, los or flux). This was determined via simultaneous X-ray spectral analysis of two archival observations for each source: one NuSTAR observation, and one soft X-ray observation (XMM-Newton, Chandra, or Swift-XRT; in this order of priority). In the analysis, both intrinsic flux and l.o.s. N_H-variability were considered whenever required to obtain a good fit to the data. Flux variability can often be confused with N_H, los variability when the data quality is low; therefore, we considered all these sources possible candidates to perform an in-depth study of N_H, los variability.

Out of the mentioned 31 sources, only 18 had additional archival data to that analyzed by Zhao et al. (2021)². Out of those, NGC 7479 was analyzed and published as a pilot project (Pizzetti et al. 2022), and Mrk 477 is currently the subject of a monitoring campaign (Torres-Albà et al., in prep.). ESO 201-IG004 is part of a double system, which is not clearly resolved in the NuSTAR data, and was therefore removed from our sample, given the sensitivity required of the proposed analysis. 4C+73.08 was also removed as the XMM-Newton observations (additional to the one used by Zhao et al. 2021) were corrupted by flares. NGC 7582 and NGC 6300 both have a large number of observations, and have been studied in depth in previous works (e.g., Rivers et al. 2015; Jana et al. 2020, respectively) regarding N_H, los variability. Both sources require a more careful comparison with previous results, which is beyond the scope of this work. In order to complete a self-consistent analysis of the whole sample, we present their in-depth analysis in future works (Torres-Albà et al., in prep., Sengupta et al., in prep.).

This leaves us with 12 sources, with a total of 54 observations. These are listed in Table 1.

3. Data reduction

The data retrieved for both NuSTAR Focal Plane Modules (FPMA and FPMB; Harrison et al. 2013) were processed using the NuSTAR Data Analysis Software (NUSTARDAS) v1.8.0. The event data files were calibrated running the nupipeline task using the response file from the Calibration Database (CALDB) v. 20200612. With the nuproducts script, we generated both the source and background spectra, and the ancillary and response matrix files. For both focal planes, we used a circular source extraction region with a 50″ diameter centered on the target source. For the background, we used an annular extraction region (inner radius 100″, outer radius 160″) surrounding the source, excluding any resolved sources. The NuSTAR spectra have then been grouped with at least 20 counts per bin.

We reduced the XMM-Newton data using the SAS v18.0.0 after cleaning for flaring periods, adopting standard procedures. The source spectra were extracted from a 30″ circular region, while the background spectra were obtained from a circle that has a radius 45″ located near the source (avoiding contamination by nearby objects). All spectra were binned with at least 15 counts per bin.

The Chandra data were reduced using CIAO v4.12 (Fruscione et al. 2006). The source spectra were extracted from a 5″ circular region centered around the source, while the background spectra were obtained using an annulus (inner radius 6″, outer radius 15″) surrounding the source, excluding any resolved sources. All spectra were binned with at least 15 counts per bin.

All spectrum extracting regions have sizes and characteristics as specified above unless otherwise stated in the source comments in Appendix C. Likewise, any exceptions on the mentioned minimum counts per bin (which ensure good usage of χ² statistics) are mentioned in the same appendix.

We fit our spectra using the XSPEC software (Arnaud 1996, in HEASOFT version 6.26.1), taking into account the Galactic absorption measured by Kalberla et al. (2005). We used Anders & Grevesse (1989) cosmic abundances, fixed to the solar value, and the Verner et al. (1996) photoelectric absorption cross-section. The luminosity distances are computed assuming a cosmology with H₀ = 70 km s⁻¹ Mpc⁻¹, and Ω/italic>_Λ = 0.73. We used χ² as the fitting statistic unless otherwise mentioned.

4. X-ray spectral analysis

All sources were fit using a physically motivated torus model, with the addition of a soft component, generally of thermal origin. Three torus models, responsible for the reflection of the AGN emission in the spectra, were used (and are described below): MYTorus (Murphy & Yaqoob 2009), borus02 (Baloković et al. 2018) and UXCLUMPY (Buchner et al. 2019). To account for the soft excess present in most galaxies, we used the thermal emission model apec (Smith et al. 2001). In multiple occasions, sources required the use of two apec components to accurately describe the soft excess. This has been shown to reproduce the complex thermal emission in star-forming galaxies (Torres-Albà et al. 2018)³.

X-ray data for each source were fit simultaneously. That is, parameters that are not expected to change in the considered timescales (of up to ∼20 yr) were linked between different observations, and thus kept a constant value. As shown in previous works, this strategy can significantly reduce the error of the common parameters (e.g., Marchesi et al. 2022). Parameters kept constant include the intrinsic photon index of the AGN (i.e. Γ) and torus geometry parameters (see individual torus model sections for details). Any caveats and/or implications of this approach are discussed in Sect. 7.

The model used is

(1)

where C accounts for intrinsic flux variability and/or cross-calibration effects between different observations; and phabs is a photoelectric model that accounts for the Galactic absorption in the direction of the source (Kalberla et al. 2005). We note that, for the purposes of this paper, we considered N_H, los free to vary at all epochs. However, this is not the case for C. In order to minimize the number of free parameters in the models⁴, we did not consider intrinsic flux variability between two observations (A and B) when: 1) χ² did not improve significantly when adding the additional free parameter (which we ensured via f-test); 2) C_A and C_B were compatible with each other within errors at 1σ; and 3) forcing C_A = C_B did not result in a source that was N_H, los variable to become non-N_H, los variable (and vice-versa).

The Soft Model can take the two following forms:

(2)

(3)

and in which kT₂ > kT₁. As mentioned above, this is a first approximation to a multiphase medium, in which the material closer to the nucleus of the galaxy is hotter, as well as more obscured (Torres-Albà et al. 2018).

The AGN Model accounts for both line of sight and reflection components, as well as a scattered component. The latter characterizes the intrinsic powerlaw emission of the AGN that either leaks through the torus without interacting with it, or interacts with the material via elastic collisions. This component is set equal to the intrinsic powerlaw, multiplied by a constant, F_s, that represent the fraction of scattered emission (typically on the order of few percent, or less).

All sources were fit in the range from 0.6 keV⁵ to 25–55 keV, with the higher energy limit depending on the point in which NuSTAR data was overtaken by the background. For every source, all models have been consistently applied to the same energy range. Results of the X-ray spectral analysis of each source can be found in Sect. 4 and Appendix A. The obtained spectra along with the simultaneous borus02 best-fit can be found, for all sources, in Appendix B. Comments on the specific fitting details of each source can be found in Appendix C.

4.1. MYTorus

The MYTorus model (Murphy & Yaqoob 2009) assumes a uniform, neutral (cold) torus with half-opening angle fixed to 60°, containing a uniform X-ray source. It is decomposed into three different components: an absorbed line-of sight emission, a reflected continuum, and a fluorescent line emission. These components are linked to each other via the same power-law normalization and torus parameters (i.e. torus absorbing column density, N_H, and inclination angle θ_i). The inclination angle is measured from the axis of the torus, so that θ_i = 0° represents a face-on AGN, and θ_i = 90° an edge-on one.

Both the reflected continuum and line emission can be weighted via multiplicative constants, A_S and A_L, respectively. When left free to vary, these can account for differences in the fixed torus geometry (i.e. metallicity or torus half-opening angle) and time delays between direct, scattered and fluorescent line photons.

We use MYTorus in “decoupled configuration” (Yaqoob 2012), so as to better represent the emission from a clumpy torus. Generally, a better description of the data is possible when decoupling the line-of-sight emission from the reflection component (e.g., Marchesi et al. 2019; Torres-Albà et al. 2021). That is, the N_H associated to absorption, N_H, los, and the N_H associated to reflection, N_H, av, are not fixed to the same value. This allows for the flexibility of having a particularly dense line of sight in a (still uniform) Compton-thin torus, or vice versa.

In this configuration, the line of sight inclination angle is frozen to θ_i = 90°. In order to better represent scattering, two reflection and line components are included. One set with θ_i = 90° (forward scattering), weighted with A_S, L90; and one set with θ_i = 0° (backward scattering), weighted with A_S, L0. In this configuration θ_i is no longer a variable. We note however that the ratio between forward to backward scattering (i.e. A_S, L90/A_S, L0), can give a qualitative idea of the relative orientation of the AGN, as it indicates the predominant direction reflection comes from.

In the particular case of fitting multiple observations together, we considered that N_H, av does not vary with time, and neither do the constants A_S and A_L. All of these parameters are representative of properties of the overall torus, which is assumed to not vary in the considered timescales. However, N_H, los can change as the torus rotates and our line of sight pierces a different material. Therefore, each individual observation is associated to a different N_H, los.

In XSPEC this model configuration is as follows,

(4)

We fixed A_S, 90 = A_L, 90 and A_S, 0 = A_L, 0, as is standard.

4.2. BORUS02

borus02 (Baloković et al. 2018) is also a uniform torus model, but with a more flexible geometry: the opening angle is not fixed, and can be changed via the covering factor, C_F, parameter (C_F ∈ [0.1, 1]). The model consists of a reflection component, which accounts for both the continuum and lines. Therefore, an absorbed line-of-sight component must be added.

We also used this model in a decoupled configuration, with N_H, los and N_H, av set to vary independently. In this case, however, θ_i (with θ_i ∈ [18 − 87]) can still be fit in a decoupled configuration. borus02 also includes a high-energy cutoff (which we froze at ∼300 keV, consistent with the results of Baloković et al. (2020), on the local obscured AGN population) and iron abundance (which we froze at 1) as free parameters. We were not able to constrain these two parameters with the data available.

When considering our variability analysis, we again allowed N_H, los to vary between different observations, but forced all torus parameters (N_H, av, C_F, θ_i) to remain constant. In XSPEC this model configuration is as follows,

(5)

where zphabs and cabs are the photoelectric absorption and Compton scattering, respectively, applied to the line-of-sight component.

4.3. UXCLUMPY

UXCLUMPY is a clumpy torus model, which uses the Nenkova et al. (2008) formalism to describe the distribution and properties of clouds. Possible torus geometries are further narrowed down using known column density distributions (Aird et al. 2015; Buchner et al. 2015; Ricci et al. 2015), as well as by reproducing observed frequencies of eclipsing events (Markowitz et al. 2014).

Clouds are set in a Gaussian distribution of width σ (with σ ∈ [6 − 90]) away from the equatorial plane. This distribution is viewed from a given inclination angle, θ_i (with θ_i ∈ [0° − 90°]).

The model consists of one single component, which includes both reflection and line of sight in a self-consistent way, allowing for a high-energy cutoff, which we again freeze at E_cut = 300 keV. Although this model has the advantage of providing a clumpy distribution of material, it does not provide an estimate of the average column density of the torus, N_H, av, which can be compared to the that provided by MYTorus and borus02. Therefore, N_H, los is the sole column density provided by the model.

In addition to the cloud distribution, UXCLUMPY offers the possibility of adding an inner ‘thick reflector’ ring of material, which was shown to be needed to fit sources with strong reflection (Buchner et al. 2019; Pizzetti et al. 2022). This material has a covering factor, C_F (with C_F ∈ [0 − 0.6]). Sources with C_F = 0 do not require this additional inner reflector.

When considering our variability analysis, we again allowed N_H, los to vary between different observations, but forced all torus parameters (C_F, θ_i, σ) to remain constant. In XSPEC this model configuration is as follows,

(6)

where uxclumpy-scattered is the scattered emission that leaks through the torus. UXCLUMPY however provides a more realistic version than a simple powerlaw, which includes the emission that leaks after being reflected.

5. Variability estimates

The main objective of this work is to measure the variability in obscuring column density, or N_H, los, for the proposed sample of sources. As such, a method to determine whether sources are variable is needed. Here, we propose two estimators of source variability. A detailed explanation on the interpretation of these comparisons for each source can be found in Appendix C.

5.1. Reduced χ² comparison

The parameters of the best-fit models to the data are reported in Table 2, and Tables A.1 through A.11. The reduced χ² () of the best-fit is reported for all three models used.

As a further test for the need to introduce variability in the models, we present a comparison with for the best fit under three different assumptions:

– There is no variability, either in intrinsic flux or N_H, los, at any epoch ( No Var).

– There is no intrinsic flux variability at any epoch, but N_H, los variability is allowed at all epochs ( No C Var.).

– There is no N_H, los variability at any epoch, but intrinsic flux variability is allowed at all epochs ( No N_H Var.).

A χ² distribution approximates a Gaussian for large values of N (number degrees of freedom), with a variance . can then be used to compare different models to select the one that best fits the data. The of the “true” model, the one with the true parameter values, is a Gaussian distributed around the mean value of 1 with standard deviation σ (see e.g., Andrae et al. 2010). A tension can then be defined between the proposed model and the data, as .

We considered that a model fit a source significantly better than another when the former had a T < 3σ, and the latter yielded T > 5σ (see e.g., Andrae et al. 2010). We used this system to classify sources as N_H, los-variable, by comparing the best-fit T with the no-N_H, los-variability T. When both models yielded T < 3σ we interpreted that N_H, los-variability is not required to fit the data, and thus classified the source as non-variable. Disagreement between the different torus models used resulted in classifying the source as ‘Undetermined’.

An exception to this rule was made for NGC 4388. No model fit the data with T < 3σ (see discussion in Appendix C), but the difference in significance between the best-fit (which includes N_H, los variability) and the non-variability scenarios is of 30 − 40σ. Therefore, we considered that including N_H, los variability results in a significant improvement to the fit, and thus we classified this source as N_H, los-variable.

We note that for two sources in our sample, NGC 612 and 4C+29.30, the fitting statistic used is a mix of C-stat and χ² (due to one or more of the spectra having very few cts/bin. See Sect. 6, and individual source comments in Appendix C). In such cases, we use T = |1 − Stat_red|/σ. However, given how this distribution does not necessarily approximate a Gaussian, the interpretation of T in such cases is not straightforward. We opt to still provide this value as a reference.

5.2. P-value

We took the derived best-fit values of N_H, los for all epochs (as depicted in Fig. 2) and estimated the probability that they all result from the same ‘true’ value. Here the null-hypothesis is that no N_H, los variability was found among different observations of the source. That is, the probability that the source is not N_H, los-variable. We did this via a χ² computation, that we later converted into a p-value (probability of the hypothesis: the source is not N_H, los variable). The χ² is generally computed as follows:

(7)

Fig. 1. borus02 fit to the data for NGC 612. Color code is as explained in Appendix B.

Fig. 2. NH, los as a function of time (data points) for MYTorus, borus02 and UXCLUMPY. Dashed horizontal lines and shaded areas correspond to the best-fit values of NH, av, and their error, respectively, for MYTorus and borus02. This quantity is considered constant with time.

Fig. 2. continued.

However, in our particular scenario, the errors of the N_H, los determinations are asymmetric (i.e. not Gaussian). In order to calculate the equivalent to Eq. (7) one needs to know (or, in its default, assume) the probability distribution of the error around the best-fit value. We followed the formalism detailed in Barlow (2003) and opted to assume a simple scenario to describe this function: two straight lines which meet at the central value. In such a case, in order to evaluate the χ² one needs only to assume as the error δ either σ⁺ or σ⁻, as appropriate.

From the obtained χ² we obtained the probability (p-value) of the null-hypothesis.

– We classified a source as N_H, los-variable if p-value < 0.01 for all three models used (MYTorus, borus02,UXCLUMPY).

– We classifyied a source as not N_H, los-variable if p-value > 0.01 for all three models used.

– We classified a source as ‘Undetermined’ if p-value was above the given threshold for at least one model, and below it for the others.

6. Results

In this section we present results on the analysis of all sources. Figure 1 shows an example of the reduced data and simultaneous borus02 best-fit for one of the sources, NGC 612. Table 2 is an example of the tabulated best-fit parameters for NGC 612. The table lists, for each of the three models used, the best-fit statistics (reduced χ² and χ²/d.o.f., i.e. degrees of freedom; or a mix of χ² and C-stat for sources with at least one spectra binned with < 15 cts/bin⁶) in the first block. It also includes the tension, T, between the data and the obtained best-fit model, derived as described in Sect. 5.1.

The second block shows parameters related to the soft emission. The third block shows the parameters corresponding to the AGN emission models. The fourth and fifth blocks refer to source variability, either of N_H or intrinsic flux (C, the cross-normalization constant), respectively. The sixth block shows the estimates of the intrinsic luminosity of the sources in two bands, 2–10 keV and 10–40 keV, for the first Chandra observation. To obtain the luminosity at any other epoch, one needs to multiply this number by the corresponding cross-normalization constant.

The final blocks show the best fit statistics that could be achieved when considering: a) No variability at all between observations; b) No intrinsic flux variability between observations; c) No obscuring column density variability between observations. For each of these scenarios, the tension between the data and the best-fit models is also computed, as described in Sect. 5.1. Finally, we computed the probability of the source being not variable in N_H, los (p-value), as described in Sect. 5.2.

Tables containing the best-fit results for the rest of the sample can be found in Appendix A. Table 3 contains a summary of the results of applying the variability determination methods described in Sect. 5 to all sources, for all three models used.

We classified a source as N_H, los-variable or as not N_H, los-variable if at least five out of six classifications (accounting for both variability estimation methods, applied on the N_H, los determinations from all three used models) agreed on the classification. If two or more determinations disagreed for any source, we classified it as ‘Undetermined’. This is the case for only two sources within the sample: NGC 612, for which borus02 resulted in variability according to both determinations; and IC 4518 A, for which the p-value and the determinations disagreed for both MYTorus and borus02. Further commentary on these disagreements can be found in Appendix C.

Following the method described above, out of the 12 sources analyzed in this work, five are not N_H, los-variable, five are N_H, los-variable, and two remain undetermined. It is worth noting that all sources required at least one type of variability (either N_H, los or intrinsic flux) in order to explain the data, as expected from our sample selection. This can be appreciated when comparing the best-fit to the no-variability in the tables presented in Appendix A.

Figure 2 shows the N_H, los variability as a function of time for all the sources analyzed, considering all three physical torus models: MYTorus, borus02 and UXCLUMPY. The dashed horizontal lines represent the best fit values for N_H, av obtained with MYTorus and borus02. The shaded areas correspond to the uncertainties associated to those values. All values of N_H depicted can be found in Table 2, and Tables A.1–A.11.

7. Discussion

Using the comparison between in the no-variability scenario and the best-fit scenario, it is easy to see that all sources in the sample required some form of variability in order to fit the data. About 42% of the sample (five out of 12) presented N_H, los variability for certain; a number that could be as high as ∼58% if all our ‘Undetermined’ cases turned out to be N_H, los variable. For five sources in the sample we can confidently say no N_H, los variability is present between the given observations.

When analyzing the results, however, one must take into account the following two factors: 1) The sample was intentionally biased toward variable sources, meaning that we expected to detect more N_H, los variability than in a blind survey. 2) The fact that we did not detect N_H, los variability for any given source does not mean it has never varied in N_H, los.

For the two ‘Undetermined’ sources, we were not able to claim whether flux variability or N_H, los variability was needed to fit the source, but we could claim that at least one of them is required. This showcases the difficulty in disentangling the two types of variability in X-ray datasets, even when dealing with nearby, bright AGN. In particular, this behavior was amplified when fitting NuSTAR data: for both 3C 445 and NGC 7319 the clumpy model UXCLUMPY favored higher flux variability and smaller N_H, los variability between other observations and the NuSTAR one, while the opposite was true for borus02 and MYTorus, the homogeneous models. It is likely that simultaneous NuSTAR and XMM-Newton observations would allow to properly disentangle the two scenarios.

7.1. Disagreement between average torus N_H and l.o.s. N_H

One of the most obvious results of our analysis can be appreciated at first glance when looking at the plots in Fig. 2. For the majority of sources, there is a large difference between the column density in the line-of-sight (at all times) and the average column density of the torus.

If one assumes that the whole (or the majority) of the torus is responsible for both obscuration and reflection, one would expect that the time-averaged value of N_H, los (i.e. ⟨N_H, los⟩) would be similar to the value of N_H, av. This is because, as the torus rotates, our line-of-sight should intercept a variety of cloud densities, representative of the density of the torus.

To estimate the feasibility that we are probing a significant fraction of the torus, we made some simple calculations. We assumed Keplerian velocities, with black hole masses in the range M_SMBH = 10⁷ − 10⁸ M_⊙ (representative of the local Universe), distances in the range 1 − 10 pc (representative of the torus scales), and timescales in the 8 − 20 yr range (representative of our sample). Under these assumptions, we estimated the torus to have rotated between 0.003 − 0.3° within the timespan of our observations⁷. At the mentioned distances, this corresponds to a physical size of 6 × 10⁻⁴ − 6 × 10⁻³ pc.

The number of works that place constraints on torus cloud/clump size (hereafter r_c) is small. For reference, we list here a few determinations and/or commonly used values in the literature. Maiolino et al. (2010) placed the most direct lower limit on cloud size, based on their X-ray observations of a whole eclipsing event (i.e. from ingress to egress). They estimated the size of the cloud head (i.e. denser, spherical region) to be r_c > 10⁻⁷ pc, while the size of the following ‘cometary tail’ of less-dense material would be r_tail > 3 × 10⁻⁶ pc. However, one must take into account these estimates correspond to a cloud placed in the broad line region (BLR), which does not necessarily have the same size as clouds orbiting the SMBH at larger distances.

Infrared emission models of patchy or clumpy tori only require the clouds to be ‘small enough’ in order to reproduce the observed MIR SEDs (e.g., Nenkova et al. 2008). X-ray clumpy models based on the previous work assume cloud sizes on the order of r_c = 2 × 10⁻³ pc (Tanimoto et al. 2019), or θ_c = 0.1′−1°. All of these are larger than the region sizes we estimated. These, however, do not necessarily correspond to observed cloud sizes, but rather to modeling or computational requirements.

The region sizes we obtained from our estimates (6 × 10⁻⁴ − 6 × 10⁻³ pc) would not correspond to the size of a single cloud, given how multiple of our sources show variability at shorter timescales. However, in order to explain why we systematically see this N_H, los variability at a level incompatible to N_H, av, this would need to be the size of the underdense/overdense region.

While this is in principle not unfeasible, one needs to take into consideration the chances of systematically looking through overdense regions (as is the case of at least six out of 12 sources), while in only one (or two, depending on the model considered for NGC 3281) are observed through underdense ones. Furthermore, one should consider that the overdense regions are so by a factor 2−10 with respect to the torus average, while the underdense regions are so by orders of magnitude (see not only IC 4518 A and NGC 3281 in this work, but also NGC 7479 in Pizzetti et al. 2022).

A study of the actual feasibility of this geometry would require: 1) A dynamical model to generate and sustain these underdense/overdense regions within a torus; and 2) An analysis of the probability of systematically observing overdense regions in a sample of 12 sources. Both of these studies are beyond the scope of this paper.

In the sections below we explore other possibilities that could explain the observed disagreement, by assuming that the material responsible for obscuration (characterized by N_H, los and, hereafter, the obscurer) and the material responsible for reflection (characterized by N_H, av and, hereafter, the reflector) are not the same.

7.1.1. Inner reflector ring

The need for an additional, thick reflector, disentangled from the rest of the torus material, has been proposed in the past. As already mentioned above, Pizzetti et al. (2022) suggested this possibility to explain the N_H, los variability curve in NGC 7479. Furthermore, the only clumpy model used in this work, UXCLUMPY, required the addition of one such thick ring to reproduce the spectrum of sources with strong reflection (Buchner et al. 2019). In fact, both IC 4518 A and NGC 7479 require this inner ring component to model the spectrum when using UXCLUMPY, which is in agreement with the large column densities invoked by MYTorus and borus02.

This theory could explain the large differences in N_H between the two structures in the torus (of factors between 10 − 100) without the need to invoke a particularly underdense region of size up to ∼0.3° through which we observe the source. It has been suggested that such a ring could correspond to a launch site for a Compton-thick cloud wind (e.g., Krolik & Begelman 1988), an inner wall (e.g., Lightman & White 1988), the inner rim of a hot disk, as seen in proto-planetary disks (e.g., Dullemond & Monnier 2010), or a warped disk (e.g., Buchner et al. 2019, 2021, particularly suitable to explain the spectrum of Circinus).

However, we note that this UXCLUMPY inner ring component was also required to explain the spectrum of NGC 612 in this work. NGC 612 does not have a particularly large torus column density, as modeled by both MYTorus and borus02. An alternative solution for the source exists, with the caveat that it gives an unreasonable intrinsic luminosity estimate for the source (see Appendix C for details).

7.1.2. Multiple reflectors

The majority of sources in our sample have a thin reflector, rather than a thick one. This is of particular interest, given how even if one assumes a disentangled thinner reflector near the SMBH, one needs to explain why then the thicker cloud distribution does not reflect.

Figure 3 shows the overall X-ray spectrum in the 1 − 50 keV range resulting from an obscured l.o.s. (with N_H, los = 10²⁴ cm⁻², in red), a scattered component (with F_S = 10⁻², in green), a medium-thick reflector (with N_H, av = 10²⁴ cm⁻², in blue), and a thin reflector (with N_H, av = 10²³ cm⁻², in cyan).

Fig. 3. borus02 AGN X-ray spectrum resulting from an obscured l.o.s. (NH, los = 1024 cm−2, in red), a scattered component (FS = 10−2, in green), a medium-thick reflector (NH, av = 1024 cm−2, in blue), and a thin reflector (NH, av = 1023 cm−2, in cyan). We use Γ = 1.8, CF = 0.5 and cos(θObs) = 0.5.

As can be appreciated in the model, thin reflectors have more significant contributions in the 2–5 keV range, where the line-of-sight component (in the case of heavily obscured AGN) does not contribute. The medium-thick reflector, while also having a minor contribution in that range, has a shape more similar to that of the line-of-sight component. It is thus possible that when only one reflector is considered, the thin reflector is made necessary by the detected emission in the 2–5 keV range. However, the medium-thick reflector, if present, could could be more difficult to recognize given the degeneracies with the combined contribution of the line-of-sight component and the thin reflector.

While this possibility is brought forward when observing the spectra in Fig. 3, it must be thoroughly tested. We propose to do that in future works, using sources with good quality data, in which we may be able to disentangle the three components.

If such was the case, the idea of a two-phase medium (as propsoed by e.g., Siebenmorgen et al. 2015) could explain the observations: a thinner, inter-cloud medium could act as the thin reflector, while the cloud distribution itself would be the medium-thick reflector.

7.2. Torus geometry as a function of variability

Figure 4 shows a series of histograms, which showcase how certain torus properties depend on source variability. We computed the plots by averaging a given parameter for sources in each of the three variability categories defined (i.e. Variable, Not Variable, and Undetermined).

Fig. 4. Histograms containing the averaged best-fit properties of all sources in the sample, grouped by variability class. All models providing the plotted parameter are shown (MYTorus in blue, borus02 in orange, UXCLUMPY in red). Source properties are as follows: Top left, time average of all NH, los (i.e. average value of the obscurer column density) for each single source. Top right,NH, av (i.e. column density of the reflector) considered constant with time. Middle left, absolute value of the difference between the two properties plotted above. Middle right, cosine of the inclination angle, θObs. Bottom left, covering factor of the torus. Bottom right, dispersion of the torus cloud distribution.

Each of these categories contains a low number of sources (particularly, we only classify two sources as ‘Undetermined’, which results in large error bars), and thus we are unable to make strong claims about torus geometry differences for (N_H, los-) variable and non-variable sources. However, possible trends are seen in the plots in Fig. 4, which should be further explored with an increased sample size.

The top, left panel of the figure shows the histogram for the average value of N_H, los across time. Meaning, the average column density of the obscurer. We observe a tendency for N_H, los-variable sources to have thicker obscurers compared to their non-variable counterparts. However, the p-value of a t-test is between 0.20–0.25 (depending on the model), which is much larger than the p < 0.05 required to claim a significant difference.

When it comes to the average torus column density, N_H, av, this trend is not necessarily maintained. When considering the MYTorus results, we find overall thin reflectors for the whole sample, as already mentioned. However, the results are apparently different when considering borus02. We note that the error bar of the borus02 bar for Variable sources is particularly large, and that the high average value is largely due to the borus02 model yielding N_H, av > 10²⁵ cm⁻² for a single source (NGC 3281, but also IC 4518 A for the Undetermined sources data point).

This effect is similarly present in the center, left plot. In here, we show the absolute value of the difference between the N_H of the obscurer and that of the reflector. The large value and large error bar of borus02 are again due to the two sources mentioned above. However, MYTorus also suggests a larger difference between the absorber and the reflector for variable sources. Meaning, non-variable sources are more consistent with having homogeneous tori. However, the p-value of a t-test is 0.12 for MYTorus and 0.33 for borus02.

We see no significant difference between inclination angles for the two different source populations. This means the observed variability (or lack thereof) is not a result of relative orientation.

We again see no difference between the two samples when it comes to C_F, as determined by borus02. However, some difference is present when considering σ_T, as determined by UXCLUMPY. This is interesting, as both parameters are representative of the height of the material responsible for reflection. It is not obvious what could be the cause of such discrepancy, but it likely lays in the different shapes assumed for the reflector: for borus02, a homogeneous sphere with two conical cut-outs; for UXCLUMPY, a cloud distribution of different densities. UXCLUMPY thus already contains the ‘multiple reflector’ concept, and is perhaps more representative of the whole shape of the torus. If we assume, however, that borus02 only models the thin reflector, the actual C_F of the medium-thick material is left unknown. In any case, UXCLUMPY results suggest that N_H, los-variable sources have broader cloud distributions. However, the p-value of a t-test is 0.16, meaning this trend is also not significant enough with the current data-set.

Previous work by Marchesi et al. (2022) successfully used a small borus02C_F to select a variable source, NGC 1358. They argued that, in some cases, as small C_F can represent a patchy and broad cloud distribution, rather than a homogeneous and flat one. If the theory is correct, one should expect a difference in the average values for variable and non-variable sources. However, once again, the discrepancy may be due to our inability to model all reflectors in the source.

We observed no clear difference in average X-ray luminosity among the three different populations.

7.3. Δ(N_H, los)vs. Δ(t)

Figure 5 shows the change in N_H, los between any two observations of the same source, as a function of the time difference between said observations. We opt to show results of only one model, borus02, in order to make the plot more easily readable.

Fig. 5. Δ(NH, los) as a function of Δ(t). Top: borus02-obtained values of Δ(NH, los) between all observation pairs for each source, as a function of the time difference between said observations. Bottom: fractional difference in NH, los between all observation pairs for each single source, with respect to the minimum NH, los of the two, as a function of the time difference between said observations. The 25%, 50% and 75% quartiles of the distribution (Q1, Q2, Q3, respectively) are also shown as dashed horizontal lines.

As can be appreciated in the figure, while small changes in N_H, los can be observed at all given time differences between observations (Δ(t)∼1 − 5000 days), large changes in N_H, los (Δ(N_H, los) > 50 × 10²² cm⁻²) are only observed with large Δ(t) (> 100 d).

This is likely a consequence of the fact that individual clouds are not homogenous in N_H (as already shown for BLR clouds by e.g., Maiolino et al. 2010), but rather present a density gradient toward their centers. Performing calculations similar to those in Sect. 7.1, imposing that a Δ(t) > 100 d is needed for a significant change in N_H, los implies clouds are generally larger than r_c > 6  ×  10⁻⁶ − 2  ×  10⁻⁵ pc, depending on underlying assumptions (such as black hole mass and cloud distance to the black hole).

Considering that events with Δ(N_H, los) > 50 × 10²² cm⁻² are still rare for Δ(t) < 1000 d, one could further infer that the majority of clouds have minimum sizes r_c > 6 × 10⁻⁵ − 2 × 10⁻⁴ pc. The lower limits we derive are ∼2 − 60 times larger than the ones for the ‘cometary tails’ of BLR clouds obtained by Maiolino et al. (2010).

However, this estimate is highly dependent on the fact that the majority of timescales probed are at Δ(t) > 1000 d. A much larger sample than the one considered in this work is needed to fully populate the plot in Fig. 5 and derive more reliable constraints on typical torus cloud size.

Figure 5 also shows the fractional change of N_H, los between two different observations (i.e. normalized to the value of the lowest N_H, los in each pair). The tendency to higher variability with larger timescales is maintained. We also provide the three quartiles for the ΔN_H, los/N_H, los distribution for the whole sample, which are Q1 = 0.1, Q2 = 0.36, Q3 = 0.94. Meaning, the median percentual variability between any two random observations of the same source is ∼36% (with respect to the observation with lowest N_H, los). For a quarter of the observation pairs in the sample, the increase is above 100%.

7.4. Constant parameters and treatment of reflection

In order to fit the data across multiple observations, we assumed that the following parameters remain unchanged across time: Γ for all three models, N_H, av for MYTorus and borus02, θ_Obs and C_F for borus02 and UXCLUMPY, and σ_T for UXCLUMPY. The inclination angle of the torus with respect to the observer, θ_Obs is not a quantity that is expected to change with time. Similarly, due to the large scale of the torus (∼1 − 10 pc), its overall geometry is not expected to vary significantly in timescales of up to ∼20 yr. Therefore, all parameters associated to the reflection component (N_H, av,C_F,σ_T), can be considered constant across different observations.

A recent work on multiepoch observations of NGC 1358 performed by Marchesi et al. (2022) found that fitting the torus parameters individually at each epoch produced results that were compatible with those of the joint fit, but with much higher uncertainties. This is compatible with our assumption. We note that an equivalent test cannot easily be performed unless one possesses multiple sets of simultaneous XMM-Newton and NuSTAR observations, which is unlikely to be the case for any other source.

For a handful of sources in the literature, with extremely good data quality, further tests on the treatment of the reflection component may also be performed. One such example is NGC 4388 in this work, which is not well-fit under our assumptions. While large variations of torus geometry still seem unlikely, other assumptions are present in our treatment of reflection. One of them is the already-discussed assumption of one single reflector. As such, NGC 4388 is a good candidate for a future study including multiple reflectors. Another assumption lays in the relation between the normalization of the line-of-sight component and the reflection component. In the analysis of obscured AGN, the widely used assumption is that the two components have the same normalization (e.g., Baloković et al. 2018; Marchesi et al. 2019; Zhao et al. 2021; Torres-Albà et al. 2021; Esparza-Arredondo et al. 2021; Tanimoto et al. 2022). However, due to the nonsimultaneous origin of the intrinsic and the reflected emission, this is not necessarily the case. In sources with very large flux variability, it is possible that the normalization of the reflection component corresponds to a past flux level of the intrinsic emission. We will explore these possibilities for sources with good data quality in the future.

We also assumed that the photon index did not vary between different observations. While some works have suggested variability of Γ with strong luminosity variability in AGN (e.g., Connolly et al. 2016)⁸, we note that none of the sources for which we had multiple NuSTAR observations suggested a need for Γ variability. Furthermore, we did not observe extreme intrinsic luminosity variability for the sources in this sample⁹.

7.5. Agreement with previous results and model comparison

Our results show satisfactory agreement with those obtained by Zhao et al. (2021). However, for 4/12 sources we obtained N_H, av values that are incompatible with (and in 3 sources, much lower than) those of their work. This could be a result of introducing the 0.5 − 2 keV emission into the fit, which Zhao et al. (2021) did not do. If the hypothesis of the thin reflector is correct, this could result in a different sub-component disentanglement needed to explain the emission at around ∼2 keV. Alternatively, it could also mean that a larger number of observations is needed to break degeneracies between parameters, and obtain reliable values of N_H, av (i.e. not pinned at the model hard limit).

Within our sample, there is reasonable agreement within the three used models. The most notable differences are the following:

– As already mentioned, borus02 has a slight tendency to move to very large values of N_H, av, sometimes even pegged at the upper limit, in sources for which MYTorus suggests more moderate densities.

– UXCLUMPY may favor scenarios in which, instead of higher obscuration, a combination of lower obscuration and lower intrinsic flux is preferred. This is particularly true for NuSTAR data (see Fig. 2, sources 3C 445 and NGC 7319).

– The three models tend to give slightly different N_H, los results. While the agreement is still remarkable, and very often the values stay within errors, Fig. 4 (top, left) shows a systematic trend between the three models. MYTorus yields the highest N_H, los values, followed by borus02 and further followed by UXCLUMPY, with the lowest values. Interestingly, this is in disagreement with the results obtained by Saha et al. (2022, see their Fig. 13), who saw large agreement between MYTorus and borus02 while UXCLUMPY had a tendency to yield larger N_H, los values. Both our results and theirs, however, agree that these differences tend to remain small.

8. Conclusions

In this work we have analyzed multiepoch X-ray data for a sample of 12 local Compton-thin AGN, selected from the work of Zhao et al. (2021). We have derived the amount of obscuring column density in our line-of-sight (N_H, los) for each source, for each epoch available. We have also obtained values of the average torus column density, N_H, av, covering factor, C_F, inclination angle, θ_Obs, and cloud dispersion, σ_T, among others. In this section we summarize our main conclusions:

– At least 42% (five out of 12) sources in the sample present N_H, los variability (through the available observations). All sources required some form of variability, either in flux, in N_H, los, or both. This is expected, given how the sample was selected to target variable sources.

– For the sources in this sample, the median variation in N_H, los for any two observations of the same source is of ∼36% (with respect to the lowest N_H, los value in the pair).

– The majority of sources show strong disagreement between the time-average of N_H, los (or average density of the obscurer) and N_H, av (average density of the reflector). This behavior is particularly strong in N_H, los-variable sources. The difference between the two oscillates between a factor of ∼2 − 100.

– Based on the previous point, if the reflector and the obscurer are the same (and representative of the density of the torus), we must be observing the torus through overdense/underdense regions. We estimate those to have angular sizes between 0.003 − 0.3° (i.e. 6 × 10⁻⁴ − 6 × 10⁻³ pc). These regions would have to contain a number of clouds of different densities to explain the observed N_H, los variability at shorter timescales. Furthermore, it is unclear how statistically feasible it is that we observe six out of 12 sources through underdense regions, while observing only one (or two) through an overdense one. It is equally unclear if such structures are dynamically feasible.

– We provide alternative explanations to the disagreement between N_H, los and N_H, av. These imply the possibility that the material responsible for reflection and the material responsible for obscuration are not the same. We suggest the possible presence of an inner, thicker ring for sources with N_H, av > N_H, los. We suggest the possibility of a two-phase medium (or the presence of multiple reflectors) for sources with N_H, los > N_H, av.

– We observe a tendency for N_H, los-variable sources to have, on average, larger obscuring density (i.e. N_H, los) and broader cloud distributions than their non-variable counterparts. These trends however are not significant up to a 95% confidence level, and thus a larger number of sources is needed to confirm (or disprove) the claims.

– We observe no difference between inclination angle or torus covering factors for variable and non-variable sources.

– We observe small changes in Δ(N_H, los) at all timescales, but we only observe large changes (Δ(N_H, los) > 50 × 10²² cm⁻²) at large timescales (> 100 d). This suggests clouds are extended, with a density profile increasing toward their centers. While this is not unexpected, we use these numbers to place rough constraints on minimum cloud sizes. We obtain that, even in the most rapid variability scenarios, r_c > 6 × 10⁻⁶ − 2 × 10⁻⁵ pc for smaller clouds. And, for the majority of cases, r_c > 6 × 10⁻⁵ − 2 × 10⁻⁴ pc. However, we note that these estimates are highly dependent on availability of observations spanning smaller timescales.

– We observe a tendency for UXCLUMPY to result in systematically lower N_H, los values than MYTorus and borus02. This is in disagreement with behavior observed in previous works.

Future work will extend this analysis to include the following: 12 more sources, for which new observations have been taken since 2019 (Pizzetti et al., in prep.); NGC 6300 (Sengupta et al., in prep.), Mrk 477 and NGC 7582 (Torres-Albà et al., in prep.) and NGC 4507 (Cox et al., in prep.). This will result in the completion of the ∼30 source sample of variable sources selected from Zhao et al. (2021). We will further expand the sample by selecting potential N_H, los-variable galaxies by applying the newly developed method of Cox et al. (2023).

Acknowledgments

N.T.A., M.A., R.S., A.P. and I.C. acknowledge funding from NASA under contracts 80NSSC19K0531, 80NSSC20K0045 and, 80NSSC20K834. S.M. acknowledges funding from the INAF “Progetti di Ricerca di Rilevante Interesse Nazionale” (PRIN), Bando 2019 (project: “Piercing through the clouds: a multiwavelength study of obscured accretion in nearby supermassive black holes”). The scientific results reported in this article are based on observations made by the X-ray observatories NuSTAR and XMM-Newton, and has made use of the NASA/IPAC Extragalactic Database (NED), which is operated by the Jet Propulsion Laboratory, California Institute of Technology under contract with NASA. We acknowledge the use of the software packages XMM-SAS and HEASoft.

References

Appendix A: X-ray fitting results

This Appendix is a compilation of tables showing the best-fit results for all sources analyzed in this work (except for NGC 612, which can be found in Table 2, in the main text).

Appendix B: Source spectra

In this section we present the best fit borus02 models to the multiepoch spectra of all sources in the sample, shown in Figs. B.1 and B.2. We opt to show the borus02 fits over those of the other models, since MYTorus has a reflection component divided into four different individual subcomponents, which makes the spectra much more difficult to interpret. UXCLUMPY, on the other hand, does not show a distinction between l.o.s. and reflection components, therefore providing less information in the spectral decomposition. The spectra shown in Figs. B.1 and B.2 should be read as follows:

• All observations for a single source are shown together, each one in a different color. Meaning, all detectors in the same telescope are colored the same in each individual observation (i.e. MOS1, MOS2, PN for XMM-Newton, and FPMA, FPMB for NuSTAR).

• Soft band observations (XMM-Newton and Chandra) are colored chronologically, as listed in Tables 2 and A.1-A.11. The color order is as follows, from first to last observation: Black, red, green, blue, cyan, magenta.

• Hard band observations (i.e. NuSTAR) are colored, also chronologically, but separated from the soft-band observations. This is done to avoid confusion between different bands. From first to second, the colors are gray and orange.

• For each individual observation, we plot the overall best-fit model as a solid line, the l.o.s. component as a dashed line, the reflection as a dotted line, the scattering as a dot-dash line, and the soft emission component (single or double mekal and any added lines) as a dash-dot-dot-dot. We note that 3C 452 has a jet component instead of a soft component + scattering, and we use a dash-dot-dot-dot (equivalent to the soft emission component) to represent it.

Fig. B.1. From left to right, top to bottom: borus02 fits to the data for NGC 788, NGC 833, NGC 835, 3C 105, 4C+29.30, NGC 3281. Color code is as explained in Appendix B.

Fig. B.2. From left to right, top to bottom: borus02 fits to the data for NGC 4388, IC 4518 A, 3C 445, NGC 7319 and 3C 452. Color code is as explained in Appendix B.

Appendix C: Comments on individual sources

In this section we provide a detailed explanation about specific analysis and fitting details for each source, that may deviate (or need clarification) from the methods described in sections 3 and 6. We also comment on the fitting results for each specific source, add comments on model comparison if discrepancies are present, and compare the obtained fitting parameters to those obtained by Zhao et al. (2021), from which this sample is selected, and who used borus02 on only two observations per source.

C.1. NGC 612

Data reduction/fitting: C-statistic was used to fit Chandra observations 1 and 2, given how the data quality forced us to bin them with 3 and 5 cts/bin, respectively. Table 2 thus refers to Stat. (total statistic, a mix of χ² and C-statistic) instead of χ². apec was applied to model solely the XMM-Newton data, as the Chandra data did not show any excess (again, probably due to the lower quality data).

Analysis of results: All models fit this source well, and our results are compatible with those derived by Zhao et al. (2021). The best-fit values for the torus parameters are in good agreement, within errors, for all models. However, that is not the case when it comes to the variability determination. While all models required some form of variability (T> 10σ for the non-variability scenario), MYTorus is not able to discern between a pure N_H, los variability scenario and a pure flux variability with enough significance. borus02, on the other hand, clearly favors an N_H, los-variable scenario¹⁰. And finally, UXCLUMPY favors a scenario in which the spectral variability is predominantly caused by intrinsic flux changes, rather than absorption. Another significant difference between model determinations is the needed addition of a thicker, inner reflection ring by UXCLUMPY. Generally, when this component is needed, it is an indicator of strong reflection, which both MYTorus and borus02 reproduce with a high N_H, av (see e.g., Pizzetti et al. 2022, and IC 4518 A in this work). This source is an exception to the trend, making the case for the existence of a strong reflection component unclear. We note that an alternative solution exists for this source, without the additional inner ring component, but with an unreasonable normalization value (i.e. resulting in an intrinsic luminosity estimation two orders of magnitude above that of the other two models). Given the two different variability determinations between MYTorus and borus02, and UXCLUMPY, we classified this source as ‘Undetermined’.

C.2. NGC 788

Data reduction/fitting: Three Gaussian lines (zgauss in xspec) were added to model the source soft emission. The reduced χ² showed significant improvement for all models, justifying this decision (1.24 to 1.13 for MYTorus, 1.27 to 1.13 for borus02 and 1.29 to 1.17 for UXCLUMPY).

Analysis of results: The models and the data show a more significant tension than for the majority of sources in this sample, at around the 3σ level. For this source we present two borus02 configurations that can explain the data with the same goodness of fit. The two configurations can be described as a low-N_H, av scenario and a high-N_H, av one. The former is statistically preferred by MYTorus, which cannot reproduce the latter without forcing N_H, av to stay at a very high value. UXCLUMPY, while not directly comparable (it does not provide a value for N_H, av), results in values of N_H, los that are more similar to those of the high-N_H, av borus02 option. Given how the first configuration is practically identical to the MYTorus results, we opt to show the second borus02 configuration (the high-N_H, av scenario) in all plots regarding the source. The degeneracy between the reflection and line-of-sight component modeling results in different estimates for N_H, los for each model, although the upward trend of N_H, los vs time is maintained (see Fig. 2).

The analysis of Zhao et al. (2021) favored the high-N_H, av scenario, and preferred pure flux variability over the pure N_H, los variability depicted here. However, as shown by our comparisons, either option can explain the data at a similar level for all models. UXCLUMPY is the only model that, when considering the p-value determination, flags this source as variable. This is likely due to the smaller errors and slightly larger differences between N_H, los values at different epochs, compared to the MYTorus and borus02 results. However, given how the comparison doesn’t show a significant preference for N_H, los variability over intrinsic flux variability, we classified this source as ‘Non-variable in N_H, los’.

C.3. NGC 833

Data reduction/fitting: NGC 833 is part of a closely interacting system with NGC 835 (separation ∼1′). The second Chandra observation (Obs. ID: 10394) considered for this merging system does not include NGC 833, but rather only NGC 835. We opted to add this observation to the table (with blank data) to avoid confusion with the epochs shown for NGC 835. Similarly, in the XMM-Newton observation we used data from only the MOS modules, as NGC 833 falls on a prominent CCD line on the PN observation. The NuSTAR extraction region was limited to 40″ to avoid contamination from NGC 835. For the same reason, the background was extracted from a circular region (instead of the usual annulus) of radius 60″. Nearby source NGC 838, a starburst galaxy at ∼3.5′ from NGC 835, shows no NuSTAR emission, and therefore is not contaminating the spectrum. The Chandra spectrum was also extracted from a circular region (15″) radius, to avoid contamination.

Analysis of results: This source is well-fit by all models. The torus parameters are highly unconstrained, likely due to a very subdominant reflection component (see e.g., Torres-Albà et al. 2021). The comparison shows, for all models, that N_H, los variability is unnecessary to explain the data. Likewise, the p-value of all N_H, los being the same is large enough that one cannot rule out the hypothesis. Thus, we classified this source as ‘Non-variable in N_H, los’.

C.4. NGC 835

Data reduction/fitting: The NuSTAR extraction region was limited to 40″ to avoid contamination from NGC 833. For the same reason, the background was extracted from a circular region (instead of the usual annulus) of radius 60″. Nearby source NGC 838, a starburst galaxy at ∼3.5′ from NGC 835, shows no NuSTAR emission, and therefore is not contaminating the spectrum. The Chandra data was taken using a larger-than-usual 8″ circular region to include all the soft emission (this source is a known Luminous Infrared Galaxy, or LIRG), for easier comparison to the XMM-Newton data. Again, the background was extracted from a circular region (15″) radius, to avoid contamination. To fit the soft emission in this source we tried both adding Gaussian lines, or adding a second apec component (justified by this source being in a merging system, as well as a known LIRG, see Torres-Albà et al. 2018, for details). Adding two lines improved the χ² over adding a second apec, and the apec addition resulted in inverted temperatures (i.e. the ‘cooler’ gas was more obscured that the ‘hotter’ gas, which is physically implausible). We thus opted to use the Gaussian lines.

Analysis of results: The data is well-fit by all models, which are in reasonable agreement. However, the best-fit values for cos(θ_Obs) derived with borus02 and UXCLUMPY are incompatible. The former favors an edge-on configuration, while the latter favors an almost face-on one. Our results are compatible with those of Zhao et al. (2021), whose analysis also favors an edge-on scenario. The N_H, los determinations are also in perfect agreement with those of (González-Martín et al. 2016), for the Chandra data. All models agree that this source shows significant N_H, los variability. We classified this source as ‘variable in N_H, los’.

C.5. 3C 105

Data reduction/fitting: No issues to report.

Analysis of results: The data is well-fit by all models, which are in good agreement. Our results are also consistent with those of Zhao et al. (2021). Introducing N_H, los variability is not necessary to explain the data, and the p-value is also > 0.01 for all models. We thus classified this source as ‘Non-variable in N_H, los’.

C.6. 4C+29.30

Data reduction/fitting: The Chandra data shows a complex morphology in the soft band, including a jet further out from the nucleus (see e.g., Siemiginowska et al. 2012). The usual 5″-radius source region was used, but the background was extracted from a nearby 10″-radius circle, rather than an annulus, in order to avoid contamination. Furthermore, Chandra observation 1 has low quality, which forced us to use 5 cts/bin, and fit with C-statistic. The table shows therefore total Stat. instead of χ². The XMM-Newton emission was extracted as usual (avoiding the jet emission), but the larger region (needed to include the XMM-Newton PSF) resulted in including a larger fraction of hot gas. An additional constant was used to weight the normalization of apec, but both Chandra and XMM-Newton data were compatible with having the same exact kT. A second XMM-Newton observation exists (Obs. ID: 0504120201) which was not used, as it fell on the same day as the used XMM-Newton observation (Obs. ID: 0504120101) and was much shorter (see e.g., Sobolewska et al. 2012). All emission at > 2 keV originates in the nucleus, therefore the NuSTAR data is not affected by the jet presence.

Even though the cross-normalization constants are compatible with 1 within errors, forcing them all to stay equal to 1 resulted in meaningful shifts in N_H, los. Therefore, we opted to leave the necessary ones free to vary in this case.

Analysis of results: The data is well-fit by all models, which are in good agreement. We note that Chandra observations 2−5 took place within ∼1 week, which likely explains the lack of flux/N_H, los variability among those observations. While it is clear from the comparison that the data requires some form of variability (T > 20σ), neither intrinsic flux nor N_H, los variability is preferred over the other. The one exception to this is perhaps borus02, which shows a tension of > 3σ between model and data when no N_H, los variability is allowed. This is likely due to the high obscuration the model predicts for the XMM-Newton observation. In any case, the tension is not significant enough, and we classified this source as ‘Non-variable in N_H, los’.

C.7. NGC 3281

Data reduction/fitting: The Chandra data was extracted using a circle of radius 10″ (background region, annulus 11−20″) to include all the extended emission (thus, making the comparison with the XMM-Newton data easier). An additional NuSTAR observation exists that was not public at the moment this analysis took place.

Analysis of results: The data is well-fit by all models, although they are not in strong agreement: MYTorus favors a low-N_H, av scenario, while borus02 favors a high-N_H, av one. Both models are able to find an equivalent scenario to the best fit of the other, although with worse statistics ( = 1.14 for a MYTorus configuration with high N_H, av, and = 1.09 for a borus02 one with low N_H, av). Our borus02 best-fit is consistent with the results of Zhao et al. (2021).

The models show significant disagreement in the best-fit values of N_H, los, probably arising from different disentanglements of the degeneracy with Γ and N_H, av. borus02 and UXCLUMPY show the best agreement, although the NuSTAR observation is significantly more obscured in the UXCLUMPY best fit. MYTorus, on the other hand, generally prefers higher obscuration. However, the NuSTAR observation is compatible with the borus02 determination. Overall, this results in UXCLUMPY painting a much more variable picture of the source. In any case, all models agree that the source is indeed ‘N_H, los variable’, and we thus classified it as such.

C.8. NGC 4388

Data reduction/fitting:Chandra observations with Obs. ID 9276, 9277 and 2983 were not considered because they used HETG/LETG grating. This galaxy has a prominent extended emission, likely the result of star formation. We used a 12″-radius region (background annulus at radii 20−30″) to include it all in the analysis of Chandra data. The brightness and closeness of this galaxy results in great data quality, and therefore more substructure is appreciated in the soft emission. We used a two-apec model to describe it.

We note that the third XMM-Newton and the second NuSTAR observations took place simultaneously. Another XMM-Newton observation (Obs. ID: 0110930301) was not included, as it was completely affected by flares.

Analysis of results: The best-fit of all models to the data shows significant tension (T∼20σ). This may be a result of the large number of counts available for this source, compared to that of the rest of the sample. It may be that our model is too simple to adequately fit it. However, no obvious problem is seen in the fit residuals that may point toward any specific issues. This source is likely a good candidate to implement a more complex treatment of the reflection component, such as the scenarios mentioned in Sect. 7.

Despite the poorer fit, the models show remarkable agreement, particularly in the N_H, los determinations. The largest discrepancy is in the photon index obtained by UXCLUMPY, which is largely incompatible with those of borus02 and MYTorus. The values of θ_obs obtained via UXCLUMPY and borus02 are also incompatible, with UXCLUMPY favoring an edge-on scenario, while borus02 suggests a much more inclined viewing angle. Our borus02 results are mostly in agreement with those of Zhao et al. (2021), although they obtain much higher N_H, av, on the order of 10²⁴ cm⁻².

Even if the fit to the data might be improved by using more complex models, it is clear that allowing both intrinsic flux and N_H, los variability significantly improves the fit. Taking this into account, as well as the derived p-values, we classify this source as ‘N_H, los variable’.

C.9. IC 4518 A

Data reduction/fitting: For the second XMM-Newton observation, MOS2 was not used as it was corrupted.

Analysis of results: The data is well-fit by MYTorus and borus02, with UXCLUMPY showing poorer statistics. This may be a result of the strong reflection seemingly needed to fit the data. In fact, this is one of the only two sources in our sample that require the addition of an inner, CT reflection ring in UXCLUMPY. This component was introduced into the UXCLUMPY model precisely because of difficulty fitting sources with strong reflection with only a cloud distribution (see Buchner et al. 2019). MYTorus and borus02 also yield large values of N_H, av, which agrees with this interpretation. This scenario is remarkably similar to that described in Pizzetti et al. (2022).

The results obtained from our best fit are consistent with those of Zhao et al. (2021), although we obtained higher values for N_H, av (∼2 × 10²⁴ cm⁻² in the mentioned work).

While both the comparison for all models and the p-value obtained for UXCLUMPY suggest the need for N_H, los variability, the p-values for MYTorus and borus02 remain above the threshold. Therefore, we classified this source as ‘Undetermined’.

C.10. 3C 445

Data reduction/fitting: We used an extraction region of 7″ for Chandra, as some extended emission is present. The background was taken from an annulus, of radii 10−20″. The source spectra shows a prominent excess at around 2 keV that is best-fit with a second, very hot apec component. It is not obvious whether star-formation, or perhaps the presence of a jet, could result in such very hot gas. Torres-Albà et al. (2018) used the two-apec model to explain the soft emisson of a large sample of U/LIRGs, and obtained a T₂ distribution of median 0.97±0.18 keV, with a long tail extending up to 4.5 keV. However, this galaxy is not classified as a U/LIRG, nor does it show obvious morphological signs of a merger (that could explain the dense star formation required). The detection of radio emission points toward the presence of a jet, as does the slightly elongated Chandra morphology. However, it is not obvious if the jet presence could justify the addition of the second apec component, from a physical point of view. We still opted to use it in the model, given how it was required to explain the data.

Analysis of results: The data is well-fit by all models, although MYTorus requires unusually large reflection constants (A_s90 and A_s0). It also results in a larger Γ than the other models. Furthermore, MYTorus and borus02 are barely in agreement in their N_H, los determinations, while UXCLUMPY results in systematically lower values (incompatible with the other models in three out of five observations). The most remarkable difference is in the NuSTAR observation, in which UXCLUMPY models the observed flux with lower obscuration than the other models, and compensates this with a lower intrinsic flux value. Precisely because of this, UXCLUMPY is the only model that classifies the source as ‘N_H, los variable’, according to the p-value. However, the comparison shows that, even for UXCLUMPY, an alternative fit exists when imposing no N_H, los variability, with T < 3σ. Therefore, we opted to classify this source as ‘Non-variable in N_H, los’.

Our borus02 results are in good agreement with those of Zhao et al. (2021), with the exception of the N_H, av, for which they obtain a much higher value of 10²⁴ cm⁻².

C.11. NGC 7319

Data reduction/fitting: We used an annulus of radii 10−20″ to extract the Chandra background, in order to avoid a nearby source. Similarly, we used a circular source extraction region of only 15″ for XMM-Newton, to avoid both extreme soft excesses and CCD lines present around the source. No source was detected in XMM-Newton Obs. ID 0021140401, and therefore it was not used in this analysis.

A double-apec model was used to fit the soft emission of this galaxy, since it is part of a closely interacting system, which is known to increase star forming activity.

Analysis of results: The data is well-fit by all models. However, UXCLUMPY yields significantly different values for N_H, los for the NuSTAR observations. Similarly to the case of 3C 445, it models the NuSTAR observed flux by using both lower N_H, los and lower intrinsic flux values. This scenario is more similar to the best-fit Zhao et al. (2021) found for the source using borus02. They detected no significant N_H, los variability between Chandra and NuSTAR, while needing a much lower intrinsic flux for the NuSTAR observation. We recovered this solution with MYTorus and borus02, with worse statistics. Interestingly, the Zhao et al. (2021)/UXCLUMPY solution is statistically the best when not accounting for the soft X-ray emission (for borus02 and MYTorus). However, this solution always has N_H, av at the maximum value allowed by the models.

Despite the mentioned differences, all models agree that N_H, los variability is required to explain the data, although this effect is larger for MYTorus and borus02. We thus classified this source as ‘N_H, los variable’.

C.12. 3C 452

Data reduction/fitting: We used an annulus of radii 12−20″ to extract the Chandra background, in order to avoid a nearby source. 3C 452 also shows diffuse, soft, (very) extended emission, which is coming from a jet (see Isobe et al. 2002). Given how the extraction region used by Chandra is smaller than that of NuSTAR and XMM-Newton, when including the jet emission it is necessary to use different jet normalization (i.e., the variation of the parameter norm_jet does not imply that the jet is varying in flux). This source did not require any cross-normalization for the AGN emission, with the exception of that associated to the jet.

Analysis of results: The data is well-fit by all models, even if they are not in perfect agreement. borus02 yields a significantly smaller Γ value, and the determinations of θ_Obs by UXCLUMPY and borus02 are incompatible within errors. borus02 favors an edge-on scenario, while UXCLUMPY favors a face-on one, although with very large errors. Additionally, UXCLUMPY results in a much harder spectrum for the jet emission, when compared to MYTorus and borus02. Our results for borus02 are compatible with those obtained by Zhao et al. (2021), with the exception of θ_Obs, which in their case results in a face-on scenario. Additionally, Zhao et al. (2021) introduced some AGN flux variability, which in our case was modeled via changes in the normalization of jet flux.

All models agree that N_H, los variability is required to explain the data, but UXCLUMPY yields smaller values for all observations. We thus classified this source as ‘N_H, los variable’.

All Tables

All Figures

	Fig. 1. borus02 fit to the data for NGC 612. Color code is as explained in Appendix B.
In the text

	Fig. 2. NH, los as a function of time (data points) for MYTorus, borus02 and UXCLUMPY. Dashed horizontal lines and shaded areas correspond to the best-fit values of NH, av, and their error, respectively, for MYTorus and borus02. This quantity is considered constant with time.
In the text

	Fig. 2. continued.
In the text

	Fig. 3. borus02 AGN X-ray spectrum resulting from an obscured l.o.s. (NH, los = 1024 cm−2, in red), a scattered component (FS = 10−2, in green), a medium-thick reflector (NH, av = 1024 cm−2, in blue), and a thin reflector (NH, av = 1023 cm−2, in cyan). We use Γ = 1.8, CF = 0.5 and cos(θObs) = 0.5.
In the text

Fig. 4. Histograms containing the averaged best-fit properties of all sources in the sample, grouped by variability class. All models providing the plotted parameter are shown (MYTorus in blue, borus02 in orange, UXCLUMPY in red). Source properties are as follows: Top left, time average of all NH, los (i.e. average value of the obscurer column density) for each single source. Top right,NH, av (i.e. column density of the reflector) considered constant with time. Middle left, absolute value of the difference between the two properties plotted above. Middle right, cosine of the inclination angle, θObs. Bottom left, covering factor of the torus. Bottom right, dispersion of the torus cloud distribution.

In the text

Fig. 5. Δ(NH, los) as a function of Δ(t). Top: borus02-obtained values of Δ(NH, los) between all observation pairs for each source, as a function of the time difference between said observations. Bottom: fractional difference in NH, los between all observation pairs for each single source, with respect to the minimum NH, los of the two, as a function of the time difference between said observations. The 25%, 50% and 75% quartiles of the distribution (Q1, Q2, Q3, respectively) are also shown as dashed horizontal lines.

In the text

	Fig. B.1. From left to right, top to bottom: borus02 fits to the data for NGC 788, NGC 833, NGC 835, 3C 105, 4C+29.30, NGC 3281. Color code is as explained in Appendix B.
In the text

	Fig. B.2. From left to right, top to bottom: borus02 fits to the data for NGC 4388, IC 4518 A, 3C 445, NGC 7319 and 3C 452. Color code is as explained in Appendix B.
In the text