© The Authors 2023

1. Introduction

Obscuration events in active galactic nuclei (AGN) have started being detected more readily in recent years. The first major event was seen in NGC 5548 (Kaastra et al. 2014) where the XMM-Newton soft X-ray flux in 2013 was 25 times weaker compared to the 2002 Chandra data. Other AGN that have undergone obscuring events include NGC 3783 (Mehdipour et al. 2017; Kriss et al. 2019), Mrk 335 (Longinotti et al. 2013, 2019; Parker et al. 2019), NGC 985 (Ebrero et al. 2016), NGC 3227 (Turner et al. 2018), and ESO 033-G002 (Walton et al. 2021). In addition, Rossi X-ray Timing Explorer (RXTE) monitoring by Markowitz et al. (2014) detected 12 new events in eight objects, significantly increasing the number observed in AGN. It is possible that many AGN have undergone obscuration at some point in their lifetime and these events are more common than first thought (Kaastra et al. 2018a). For example, multiple obscuration events were seen in NGC 3783 between 1993 and 2016 (Kaastra et al. 2018a), including the December 2016 event lasting 32 days (Mehdipour et al. 2017), while NGC 3227 has undergone many occultations lasting days, weeks and months (Lamer et al. 2003; Markowitz et al. 2014; Beuchert et al. 2015; Turner et al. 2018). Comparatively, some AGN, such as Seyfert 2 or extreme Compton-thick AGN, are persistently obscured due to their geometry or a viewing angle that prevents our line of sight (LOS) from reaching the innermost regions near the supermassive black hole (SMBH). Obscuring clouds have different properties compared to the warm absorbers (WAs; e.g. Kaastra et al. 2002; Grafton-Waters et al. 2020), ultra-fast outflows (UFOs; e.g. Risaliti et al. 2005; Tombesi et al. 2013; Nardini et al. 2015), and the emission line regions (e.g. Kinkhabwala et al. 2002; Grafton-Waters et al. 2021), such as column densities (N_H, WA = 10²⁴ − 10²⁶ m⁻²) and outflow velocities (v_out, WA= few 100 − 1000 km s⁻¹; v_out, UFO = 0.05 − 0.4c).

Since the obscuring event in NGC 5548 (Kaastra et al. 2014), we started a Swift (Gehrels et al. 2004) monitoring programme to trigger joint Target of Opportunity (ToO) observations of Seyfert 1 AGN with XMM-Newton (Jansen et al. 2001), the Hubble Space Telescope (HST; Green et al. 2012), and NuSTAR (Harrison et al. 2013). NGC 3783 was triggered in December 2016 as the Swift hardness ratio was large enough to signify an obscuration event, caused by a high-velocity wind crossing our LOS (Mehdipour et al. 2017). The column densities of the two obscuring components in NGC 5548 and NGC 3783 were N_H = 10²⁶ − 10²⁷ m⁻² (Kaastra et al. 2014; Mehdipour et al. 2017). These obscuring clouds were likely to be located within the broad line region (BLR), originating from the accretion disk (Kaastra et al. 2014; Mehdipour et al. 2017; Kriss et al. 2019).

The type 1.5 AGN, NGC 3227 (z  =  0.003859; de Vaucouleurs et al. 1991), has shown strong variability in both X-ray and UV bands, particularly during periods of high flux (Markowitz et al. 2009; Arévalo & Markowitz 2014), with a ‘steeper-when-brighter’ trend consistent with other AGN (Lobban et al. 2020). The intrinsic X-ray continuum in NGC 3227 has been found to vary across many timescales, from days to months (e.g. George et al. 1998; Gondoin et al. 2003; Uttley & McHardy 2005; Lobban et al. 2020). However, obscuration events can cause additional contributions to the observed variability in sources such as NGC 3227 (e.g. Beuchert et al. 2015). This is similar to the explanations for the variability in NGC 4151, which are both driven by the intrinsic continuum and absorbing material within the LOS (e.g. Keck et al. 2015; Beuchert et al. 2017).

Two (Markowitz et al. 2009, 2014) and three (Crenshaw et al. 2001; Beuchert et al. 2015 ionised absorbers have been found in NGC 3227 during different epochs. The variable absorption event in NGC 3227 during 2008, analysed by Beuchert et al. (2015), was caused by a mildly ionised component (log ξ ∼ 1.2 − 1.4 erg cm s⁻¹) with increasing covering fraction (C_f; from 0.7 to 0.9) and column density (N_H; from 5 to 18 × 10²⁶ m⁻²), possibly located within the BLR or inner torus (Beuchert et al. 2015). Furthermore, an XMM-Newton and NuSTAR campaign in 2016 found a rapid variability event lasting a day (Turner et al. 2018). To explain the observed spectral hardening and a depth change of the unresolved transition array (UTA; at ∼16 − 17 Å) in the RGS data, a component covering 60% of the X-ray continuum was required in the modelling (with N_H  =  5  ×  10²⁶ m⁻²; log ξ  =  2; v_out ∼ −800 km s⁻¹). This obscurer was also associated with clouds within the BLR.

During the Swift Cycle 16 ToO monitoring programme, an obscuration event in NGC 3227 was captured at the end of 2019, triggering observations with XMM-Newton, NuSTAR and HST. Strong spectral hardening and intrinsic X-ray variability was found between the two 2019 observations taken three weeks apart, along with strong soft X-ray absorption compared to the 2016 observations (see Figs. 1 and 2 in Mehdipour et al. 2021). In the first paper of this series (Mehdipour et al. 2021, hereafter Paper I), we determined the broadband continuum, giving rise to the spectral energy distribution (SED) of NGC 3227. Analysis of the archival 2006 and 2016 XMM-Newton observations were carried out by Wang et al. (2022, hereafter Paper II), where we used the broadband SED from Paper I to model the WA wind. Four WA components were found with ionisation parameters¹ ranging from log ξ = −1 to 3, and outflow velocities in the range v_out = −100 to −1300 km s⁻¹, located between the BLR and narrow line region (NLR). A missed obscuring event in 2006 was found for the first time in Paper II, while the event in 2016 (Turner et al. 2018) was re-evaluated. The 2006 obscurer was explained with one component, but the 2016 obscurer required two: log ξ = 1 − 1.9 and N_H = 10²⁶ m⁻² for both epochs, and log ξ = 2.8 and N_H = 10²⁷ m⁻² was needed for 2016 only. The obscurers were located within the BLR. In Mao et al. (2022, hereafter Paper III), we studied the 2019 obscuring event in detail, focusing on the time-averaged spectra through modelling the RGS, EPIC-PN, and NuSTAR data. Different models were applied to each 2019 observation, with one or two photoionisation components required to account for the obscurer. The estimated location for the obscurer was comparable to the BLR and inner torus. In this paper, we study the variability of the obscurer and continuum in NGC 3227, focusing on the EPIC-PN data only.

The outline of this paper is as follows. We describe how we reduced the data for analysis in Sect. 2, while in Sect. 3 we investigate the light curves of each observation and explain how we split each observation into equal length timing bins. The spectral fitting and modelling process is explained in detail in Sect. 4, before we present the results and discuss our findings in Sect. 5. Finally, we give our concluding remarks in Sect. 6.

2. Data reduction

XMM-Newton observed NGC 3227 twice during an X-ray obscured state at the end of 2019. The two observations (Obs 1 and Obs 2, hereafter) lasted 103 and 50 ks, respectively, with the details shown in the observation log (Table 1). Here we report on the analysis of the European Photon Imaging Camera pn-CCD (EPIC-PN; Strüder et al. 2001) spectrum, taken in small-window mode, modelled between 0.35–10 keV. The EPIC-PN data were reduced using the EPPROC command in the SAS (v18.0.0) pipeline, where we removed background flaring with counts greater than 0.4 counts s⁻¹ in 10–12 keV. We extracted the source spectrum using a circle of radius 40 arcseconds and the background spectra were extracted with two circular regions of 30 arcseconds in radius from a source free region on the same CCD; both with quality photons only (FLAG == 0) and using the #XMMEA_EP filter. We included both single and double events (PATTERN <  = 4). Finally, the instrumental response matrices were generated with the rmfgen and arfgen commands. The data were binned using the Kaastra & Bleeker (2016) optimal binning method and we fitted the spectra using SPEX (v3.05.00; Kaastra et al. 1996, 2018b). We use the Cash statistic (C-statistic, hereafter; Cash 1979; Kaastra 2017) for statistical significance, with all errors in this paper quoted at 1σ confidence.

3. Light curves and timing bins

3.1. Light curves

The three different energy band light curves for each observation are displayed in Fig. 1: 0.3–10 keV (red; full X-ray range), 0.3–2 keV (green; soft X-ray band), and 2–10 keV (blue; hard X-ray band). They were extracted using the EPICLCCORR command with SAS, and chosen to cover large enough energy ranges with sufficient counts within the full 0.3–10 keV band. Obs 1 (left side of Fig. 1) shows large variability over the course of the observation, whereas Obs 2 (right side) has far fewer counts and is less variable. To determine how the soft and hard X-ray energy bands varied with respect to each other, we took the ratios between the soft or hard X-ray light curves and the full 0.3–10 keV range (shown within the bottom two panels in Fig. 1). The only strong feature in Fig. 1 (bottom two panels) is during Obs 1 where, between 30 and 60 ks, the 0.3–2 keV/0.3–10 kev ratio (green line) decreases, but the 2–10 keV/0.3–10 kev (blue) ratio increases. No changes are seen in Obs 2.

Fig. 1. Light curves for Obs 1 (left) and Obs 2 (right). Top panels: Three different energy bands: 0.3–10 keV (red), 0.3–2 keV (green), and 2–10 keV (blue). Also displayed are the timing bins (TB) from which we extracted EPIC-PN spectra for our analysis (see Sect. 3.2 for details). The light curves were binned at 100 s. The errors on the energy bin fluxes are shown by the shaded regions. Bottom panels: Ratio light curves between the 0.3–2 keV (green) and 2–10 keV (blue) energy bands in the top panel divided by the full 0.3–10 keV light curve (red line). The Obs 1 ratios are on the left hand side and the Obs 2 ratios are on the right.

We also obtained the hardness ratio (HR) for each observation, which is defined as (H–S)/(H + S), where H and S are the 2–10 keV hard band and 0.3–2 keV soft band count rates, respectively. For Obs 2 (right side of Fig. 2), the ratio stays constant over the observing window. However, for Obs 1 (left side of Fig. 2), there is a strong spectral hardening between 30 and 60 ks, corresponding to the change in soft (0.3–2 keV) and hard (2–10 keV) flux ratios over this period (see left panel of Fig. 1). This is either caused by variations in the underlying continuum, or the possible variability of the obscurer, which is what we are investigating in this paper. The rest of the observation shows that NGC 3227 is in a relatively softer state (negative HR), and is consistent with the average value (red dashed line).

Fig. 2. Hardness ratio between the 2–10 keV and 0.3–2 keV energy bands, defined as (H–S)/(H+S), for Obs 1 (left) and Obs 2 (right), binned at 100 s. The errors on the hardness ratios are shown by the shaded regions. The red dashed lines show the average hardness ratio (HRavg) across each observation. For Obs 1 HRavg = −0.115, and for Obs 2 HRavg = −0.033.

3.2. Timing bins

To measure any spectral changes over the duration of the observations, we split each observation into equal length timing bins (TBs) and obtained the full 0.35–10 keV spectrum from each. Obs 1 was split up into 10 TBs of length 10.3 ks, while Obs 2 was split into 5 bins with a duration of 10 ks (see Fig. 1). The top panel of Fig. 3 exhibits a strong correlation between the averaged 0.3–2 and 2–10 keV count rates of each TB in the two observations. The bottom panel of Fig. 3 presents the average HR for each TB as a function of the averaged 0.3–10 keV count rates. The average HR is lower for each TB in Obs 1 compared to the TBs in Obs 2, but the average count rates are larger. This agrees with the overall flux change in the light curves for each observation, whereby both panels in Fig. 3 imply that NGC 3227 shows a softer-when-brighter correlation (similar to, e.g. Lobban et al. 2020; Paper II). This relation is typically associated with intrinsic variability in the continuum, either from steepening of the power law or an increased strength in the soft excess. The increased hardening in Obs 1 from TB1–TB2 to TB3–TB5, before decreasing again for TB6–TB10, is related to the Figs. 1 and 2 where the hard X-ray flux ratio (2–10 keV/0.3–10 keV) and the HR increased between 30 and 60 ks.

Fig. 3. Top: Relation between the averaged 0.3–2 keV (soft) and 2–10 keV (hard) count rates for each timing bin. Obs 1 (hollow circles) and Obs 2 (diamonds) are compared with the time-averaged values shown in red. Bottom: Hardness ratio (H–S/H+S) as a function of the average 0.3–10 keV count rates. NGC 3227 shows a softer-when-brighter trend on timescales of weeks as Obs 1 has a larger flux but lower hardness ratio compared to Obs 2.

We also obtained the 0.35–10 keV time-averaged spectrum from each observation. Figure 4 shows each TB spectrum for Obs 1 (left) and Obs 2 (right), with the time-averaged spectra shown in red as a comparison. A clear gap can be seen in Obs 1 between the spectra of TB1–TB5 and TB7–TB10, with TB6 somewhere in between, relating to the flux increase half way through the observation (see Fig. 1).

Fig. 4. EPIC-PN spectra for Obs 1 (left) and Obs 2 (right). The time-averaged spectra are shown in red and the TB numbers are presented in the legends by their respective colours.

We initially modelled the time-averaged spectrum from each observation, corresponding to the base model, before fitting each TB spectrum. This allowed us to study the obscurer and continuum parameters in order to explain the origin and nature of the X-ray variability within NGC 3227 during obscuration in 2019.

4. X-ray spectral fitting

4.1. Setting up the model

In order to analyse the spectral variations over the course of the two observations, we required the models from the previous work on this campaign. The intrinsic (unobscured) SED (Paper I) consisted of a power law (POW), a neutral reflection component that takes into account the neutral FeK[[INLINE51]] line at 6.4 keV (REFL; Magdziarz & Zdziarski 1995), a Comptonisation component (COMT; Titarchuk 1994) for the soft excess, and a disk blackbody component (DBB) for the optical to UV band. Here, we modelled the REFL component without an accretion disk profile, similar to the previous series papers, and do not assume any blurred ionised reflection that could partly explain the soft X-ray excess. Although this could be a plausible alternative explanation for the origin of the soft excess, the COMT component is sufficient for our purposes in the modelling of a broadband SED for NGC 3227 (see e.g. Paper I, in addition to Mehdipour et al. 2011, 2017, 2018; Petrucci et al. 2013). We are therefore neglecting ionised disk reflected emission.

In our modelling, the intrinsic continuum was absorbed by the obscurer, so most of the SED parameters were fixed to the values from Paper I (for example, the seed photon and plasma temperatures, T_seed = 10.2 eV and T_c = 60 eV, respectively, in the COMT component). However, the normalisations of the DBB, POW and COMT components (N_dbb, N_pow, N_comt, respectively), the reflection scaling parameter² (s), and photon index (Γ_pow), were measured for the two 2019 observations during the obscured state in Paper III; this is the obscured SED.

The initial model examined in this paper was model M1 from Paper III (see Tables 2 and 3 for Obs 1 and Obs 2, respectively), where the REFL normalisation (N_refl) and photon index (Γ_refl) values were coupled to the POW normalisation (N_pow) and photon index (Γ_pow), for the two 2019 observations. The other models studied in Paper III were also fitted and are compared in Appendix A: for M2 the N_refl is fixed to the 2016 N_pow value from Paper I, and M3 is the same as M1, but models the obscurer with two PION components. For the power law component, we applied upper (309 keV; Turner et al. 2018) and lower energy (equal to the disk energy, in this case 10.2 eV from Paper I) cutoffs. The high energy cutoff was also applied to the reflection component and fixed at 309 keV. The fitted continuum parameters in this paper were N_pow, Γ_pow, N_comt, s. In addition, we modelled the neutral gas in our Galaxy using the HOT model in SPEX; we fixed the column density and temperature at $N_{\mathrm{H}}^{\mathrm{HOT}}=2.07\times {10}^{24}$ m⁻² (Murphy et al. 1996) and T_Gal = 0.5 eV, respectively.

We also included the four WA components from Paper II, modelled with PION (Mehdipour et al. 2016a), but accounted for their less ionised nature due to the obscurer, which was also accounted for in the previous papers of this series. The de-ionised nature of the WA means only the ionisation parameters of the WAs are lowered as they receive less photons from the central source due to the presence of the obscurer. The other parameters (e.g. N_H, v_out) were assumed to be constant over time. This was also seen in NGC 5548 (Kaastra et al. 2014) and NGC 3783 (Mehdipour et al. 2017). One PION component was used to model the obscurer in both observations with the initial parameters from Table 2 for Obs 1 and Obs 2 (taken from Tables 2 and 3 in Paper III, respectively). The ionisation parameter of the obscurer in Obs 2 is significantly larger compared to the value in Obs 1 which is indicative of there being two completely separate obscuring components that are observed in Obs 1 and Obs 2, respectively, as stated in Paper III. The fitted obscurer parameters in this paper were the column density (N_H), ionisation parameter (ξ), and covering fraction (C_f). The flux changes between the 2016 and 2019 observations (see Fig. 2 in Paper I) is a result of the obscurer. The variations we observe are unlikely to be due to the known WA components (Paper II) as the ionising continuum is absorbed by the obscurer before it reaches them.

Finally, we accounted for the emission lines, produced by the warm emitter, with a single PION component in both observations. However, the ionising continuum was different to the one used for the obscurer and WA, so the intrinsic 2016 SED (same components as the obscured SED but with different parameter values) was used to ionise the warm emitter (see Paper III for details of the modelling and Table 2 in Paper I for the values). The warm emitter values were fixed to those from Table 2 (see also Table 2 in Paper III).

4.2. Fitting the spectra

We folded the model described above to the time-averaged 0.35–10 keV EPIC-PN spectrum for each observation. The initial C-statistics were C = 2425 (for 123 degrees of freedom, d.o.f hereafter) and C = 820 (for 113 d.o.f) for Obs 1 and Obs 2, respectively. However, we found that there were strong residuals between the data and model at softer energies in both observations (see the initial fit in Fig. 5). This was most likely due to the cross calibration issue between the reflection grating spectrometer (RGS; den Herder et al. 2001) and EPIC-PN instruments. Paper III modelled the RGS for the soft X-rays and EPIC-PN data for the hard X-rays, while here we model the full X-ray band using EPIC-PN data only. Therefore, to correct for the cross calibration issue we forced the EPIC-PN data to agree with RGS. We do it this way because RGS is able to resolve the O-edge, while EPIC-PN cannot. The process of how we corrected for this cross calibration effect is explained in Appendix B. The corrected C-statistics once we accounted for any instrumental effects (Appendix B) were C = 400 (123 d.o.f) and C = 142 (113 d.o.f) for Obs 1 and Obs2, respectively (see the corrected fit in Fig. 5).

Fig. 5. Initial and KNAK corrected models (see Appendix B) for Obs 1 (left) and Obs 2 (right). The initial models were from Paper III (green and pink), while the corrected models are in blue and purple, respectively. The C-statistic values are displayed, along with their respective degrees of freedom in the parentheses. The bottom panels show the residuals between the data and model. Appendix B explains the steps required to correct the models, shown in Fig. B.1.

For the time-averaged spectrum of each observation, we then fitted N_pow, Γ, N_comt, s of the 2019 continuum (obscured SED), and the column density (N_H), ionisation parameter (ξ), and covering fraction (C_f) of the obscurer component. Fitting these parameters improved the best fit to C = 221 (116 d.o.f) in Obs 1, and C = 103 (106 d.o.f) for Obs 2. We note that the best fit in Obs 1 is not formally acceptable (C/d.o.f = 1.91), while in Obs 2 the fit is very good (C/d.o.f = 0.97). However, for Obs 1 there are no clearly structured residuals left despite the formally unacceptable fit. One explanation is that the obscurer and continuum parameters for TB1-TB5 are so different from TB6-TB10 that a single model can’t explain the average spectrum. Furthermore, this is similar to the analysis of NGC 1068, where Grafton-Waters et al. (2021) found that the 2014 observation, owing to having more than 2.5 times the exposure time compared to 2000 epoch, had a significantly worse fit. The final best fit models of each observation are displayed in Fig. 6. Now that we obtained a good fit to the EPIC-PN data only, we used this as the baseline model to fit the spectra of each TB. The best fit parameter results are displayed in Table 3.

Fig. 6. Final best fit models fitted to the time-averaged spectra: Obs 1 grey crosses and blue model; Obs 2 black crosses and purple model. Also shown are the C-statistics and degrees of freedom. The best fit parameter values are displayed in Table 3. The bottom panel shows the residuals between the data and model for Obs 1 (blue) and Obs 2 (purple).

5. Results and discussion

5.1. Timing bins comparison

Here we compare the parameter results from each TB for the two observations. Table 3 displays the best fit parameter values for each TB in Obs 1 and Obs 2, along with the time-averaged spectrum results (All). Figure 7 displays how these parameters changed over the course of each observation. For each panel in Fig. 7, the time-averaged values are shown by the red dashed and solid horizontal lines for Obs 1 and Obs 2, respectively.

Fig. 7. Model parameter comparisons for each timing bin (TB) in Obs 1 (hollow circles) and Obs 2 (diamonds). The red horizontal dashed and solid lines represent the time-averaged parameter values for Obs 1 and Obs 2, respectively, and the shaded areas show the uncertainties. The TB colours are consistent with their respective spectra in Fig. 4. The model parameter is shown on the y-axis of each panel. The light purple lines (dashed for Obs 1 and solid for Obs 2) connect the TB values for each parameter in the two observations.

The top two panels in Fig. 7 display how the power law normalisation (N_pow; left) and photon index (Γ_pow; right) vary over time. For Obs 1, N_pow increases during the observation, jumping in flux between TB5 and TB7, similar to the count rate increase in the light curve of Fig. 1. For Obs 2, N_pow stays constant between TBs, within the uncertainties. The photon index in Obs 2 displays strong changes, although similar to N_pow, they are consistent with each other within the uncertainties. In Obs 1, Γ_pow does not vary throughout the observation, except for TB1 and TB2. For both observations, Γ_pow is consistent with a hard power law, similar to previous observations (Γ = 1.5 − 1.8; e.g. Gondoin et al. 2003; Cappi et al. 2006; Rivers et al. 2011; Arévalo & Markowitz 2014), although this parameter is overall slightly lower in Obs 1.

The top-middle panels in Fig. 7 show the Comptonisation normalisation (N_comt; left) and reflection scaling parameter (s; right) changes. For Obs 1, N_comt has the same, but more significant, trend as N_pow. For Obs 2, similarly to N_pow, N_comt is consistent between TBs. Although the reflection scaling parameter (s) of REFL in both observations appears to correlate well with Γ_pow (top right panel of Fig. 7), both parameter values are consistent within uncertainties, suggesting no overall correlation. For Obs 1, the reflection scaling of TB2 is larger compared to the other TBs, but from TB3 to TB10, s decreases over time, showing an anti-correlation between the flux (Fig. 1) and the amount of reflected continuum. In other words, we are observing more intrinsic flux during the latter half of Obs 1.

The bottom-middle panels in Fig. 7 present the column density (N_H; left) and covering fraction (f_cov; right) of the obscurer component. In Obs 2, it appears that N_H fluctuates significantly between TBs. However, with the exception of TB3, the N_H values in each TB are constant within the uncertainties of each other and the time-averaged value. Therefore, the given changes in N_H are unlikely to occur on timescales between bins, and we suggest that these fluctuations are not real. For Obs 1, on the other hand, N_H decreases through each TB, which follows the inverse trend compared to N_pow and N_comt, meaning at the same time there are changes in the obscurer and in the source continuum. We run further tests to explain this inverse trend in Sect. 5.2, and discuss the possible causes of this observed variability in Sect. 5.3 to exclude possible degeneracy between the model parameters. For the obscurer covering fraction, overall, f_cov stays constant within the uncertainties throughout each observation. However, f_cov is systematically higher in TB1–TB5 compared to TB6–TB10. This suggests that the covering fraction is not the parameter driving the observed variability. In Obs 2, TB3 gives apparent outliers for both N_H and f_cov compared to the other TBs. However, we are unable to constrain ξ (bottom panel of Fig. 7) in TB3, which, in addition to the lower flux of Obs 2, could explain this drop in obscurer parameter values compared to the other TBs.

Between Obs 1 and Obs 2 (∼3 weeks), the time-averaged results (Table 3) indicate that f_cov and N_H increased from 0.55 to 0.66, and from 3.87 × 10²⁶ m⁻² to 9.93 × 10²⁶ m⁻², respectively. Conversely, the observed N_pow and N_comp time-averaged parameters decreased from 2.60 to 0.72 × 10⁵⁰ ph s⁻¹ keV⁻¹ and from 1.89 to 0.16 × 10⁵³ ph s⁻¹ keV⁻¹, respectively. Similar to the obscurer results in this paper, Beuchert et al. (2015) found that the observed variability of NGC 3227 in 2008 could also be produced by an increase in covering fraction (0.7 to 0.9) and column density (5 to 18 × 10²⁶ m⁻²) in an absorption component, in addition to changes caused by the intrinsic continuum. Furthermore, in NGC 5548, the short-term variability (days to weeks) was driven by changes in both the column density and covering fraction of the obscurer (Cappi et al. 2016), whereas for the long-term variability (weeks to months), the covering fraction of the obscurer changed over time (Mehdipour et al. 2016b, 2022). Here we discover that over the course of Obs 1 (103 ks) the column density of the obscurer varied significantly in relation to the changes in the observed X-ray flux, despite ξ being free, while f_cov stayed constant. However, three weeks later, during Obs 2, both the covering fraction and column density increased, consistent with the strong X-ray flux decrease observed between the two observations. Paper III attribute this to two separate obscuring clouds occulting the LOS X-ray continuum due to large differences in the ionisation parameters. On the other hand, the continuum normalisations are significantly smaller in Obs 2, compared to Obs 1, and therefore it is not clear what produces the strong flux decrease.

Finally, the bottom panel of Fig. 7 displays the ionisation parameter (ξ) of the obscurer for each TB. For Obs 2, the TB values are constant with the time-averaged value, except for TB3 which is an upper limit. For Obs 1, on the other hand, it is very difficult to constrain ξ, even for the time-averaged value. For TB1–TB3, the uncertainties are very large, and for TB4–TB6 only upper limits can be found. The obscurer ionisation parameter can only be constrained for TB7 to TB10 when the flux is at its highest. Given the size of these uncertainties we cannot conclude one way or the other whether the ionisation state of the obscurer responds to changes in the ionising continuum (N_pow and N_comt). On the other hand, these large ξ errors could indicate an inhomogeneous obscurer that is made up of multiple components, each with a different N_H, ξ, or number density (n_e), such that any global change would be hard to identify with a single model component.

5.2. Parameter comparison

After fitting the data for each TB, we examined whether there were any correlations between the parameters. In doing so, we find that there are significant trends between N_pow and N_comt, as well as with N_H. These are displayed in the left side of Fig. 8. For Obs 1, N_pow and N_comt relate well with Fig. 1, while an inverse trend was seen with N_H (Fig. 7). In each panel of Fig. 8, the Pearson rank and p-values are shown to display how well correlated the parameters are. In all panels, the p-values are significantly lower than the quoted value and hence they are shown as p < 0.001. The absolute values of r are greater than 0.8 (both in the positive and negative directions) implying that these parameters in Fig. 8 are correlated or anti-correlated very well with each other.

Fig. 8. Correlations between model parameters. Left: The normalisations of the power law (Npow) and Comptonisation component (Ncomt) are correlated over the two observations (top); Npow against the column density of the obscurer (NH) (middle); and Ncomt against NH (bottom). Right: Correlations between Npow, Ncomt and NH against with the 0.3–10 keV count rate for top, middle and bottom, respectively. In each panel, the Pearson rank (r) and p-value is shown to display the significance of each correlation.

The right side of Fig. 8 shows the relation between N_pow, N_comt, and N_H with the 0.3–10 keV count rates, taken from the bottom panel of Fig. 3. Both normalisations are positively correlated with the 0.3–10 keV count rates (top right and middle right panels of Fig. 8), while there is an inverse correlation between N_H and the count rate (bottom right panel in Fig. 8). This possibly suggests that the observed variability is related to changes in the column density seen in Obs 1 (see Sect. 5.3). For Obs 2, the count rates are consistently low for each model parameter, which is a result of high N_H and f_cov values, and a low continuum, in Obs 2; no strong increase over the observation is seen.

We further analysed our modelling to test the statistical significance of thawing a parameter from its time-averaged value in the spectral fitting. To do this, we fixed each of the four parameters of interest – N_H, N_pow, N_comt, and f_cov – to their time-averaged values (Table 3) in turn, while all other model parameters were free, and we fitted this adjusted model to the ten TB spectra. We repeated this for each TB in Obs 1, before measuring the ratios between the C-statistic values and the d.o.f and summing them (S_fixed = ∑C_i/d.o.f._i, where d.o.f. = 117 and i is the TB number). The summed ratios for the four parameters are displayed in Table 4. As a comparison, we also sum the S_free = C_i/d.o.f._i ratios (where d.o.f. = 116 here) from the best fit TB results in Table 3, which is S_free = 11.87. We then calculate the difference between these two summations (ΔS = S_fixed − S_free), which we use to test the statistical significance each parameter has on the model; the differences are also shown in Table 4.

Next, we ran some F-tests to examine whether the variance (σ) calculated from the σ_fixed = C_i/d.o.f._i ratios in each TB when the parameters of interest are fixed, come from the same distribution as the best fit σ_free = C_i/d.o.f._i variance when all parameters in the model are free (see Table 3). The ratio between the variances give rise to the F-value, given as F = σ_fixed/σ_free. The F-value is used, along with the calculated p-value, to determine the statistical probability that the two variances come from the same underlying distribution, when compared to the difference in the summed ratios ΔS for each fixed parameter. If the F-value is larger than ΔS, then it implies that the parameter in question is statistically significant in the best fit model. The results for the F-values and corresponding p-values are displayed in Table 4.

From the results in Table 4, the F-values for N_comt and N_H are larger than the difference in the summed ratios (ΔS), for the respective parameter. In addition, N_comt and N_H have sufficiently small p-values (less than 0.05) which implies that it is unlikely for the C-statistic variances (σ_fixed) of these two parameters to come from the same distribution as the best fit C-statistic variance (σ_free). These results suggest that N_comt and N_H have the most statistically significant impact on the best model from Sect. 4.2 and Table 3. We also find that the F-value for f_cov is larger than the ΔS measure. However, the p-value for f_cov is greater than 0.05 (Table 4), which means we cannot conclude whether this parameter has a significant impact on the best fit with certainty, as this result indicates that it is likely for the C-statistic variance from f_cov to lie on the same distribution as the best fit model variance. Finally, the F-value for N_pow is smaller than the ΔS value and therefore means it does not have a significant impact on the best fit model in Obs 1, compared to N_comt and N_H.

To explore the anti-correlation between the obscurer and continuum (left side of Fig. 8), we create some contour plots of N_H against N_pow and N_comt (Fig. 9) for the time-averaged spectrum in Obs 1, and for the spectra of TB2, TB6, and TB9 as they correspond to a low flux bin, a change in flux bin, and a high flux bin, respectively (see Fig. 1). The far left panels in Fig. 9 show that the time-averaged parameters (N_pow top, and N_comt bottom) are relatively well constrained within the 68 % confidence level (black solid lines). This suggests that there are no significant degeneracies between these three parameters. For TB2 (middle left panels) and TB6 (middle right panels), the contours are less smooth and round compared to the time-averaged plots. The contours for TB9 (far right panels), however, display a similarity with the time-averaged contours. A possible explanation for this is that TB2 and TB6 have lower fluxes compared to TB9, and therefore a lower signal-to-noise ratio (S/N).

Fig. 9. Contour plots of NH (1026 m−1) against Npow (1050 ph s−1 keV−1; top panels) and NH (1026 m−1) against Ncomt (1053 ph s−1 keV−1; bottom panels). From left to right, the panels show the contours for the time-averaged model, TB2, TB6, and TB9, respectively, in Obs 1. The levels of confidence presented are at 68 % (ΔC = 2.30; black solid line), 90 % (ΔC = 4.61; red dashed line), and 95 % (ΔC = 6.17; blue dotted line). The black crosses show the parameter values (but not the errors) from the best fit models in Table 3.

There are slight differences between the widths of the 68% confidence contours (black lines) of N_pow (0.48 to 0.55 × 10⁵⁰ ph s⁻¹ keV⁻¹) and N_comt (0.40 to 0.48 × 10⁵³ ph s⁻¹ keV⁻¹) for TB2, TB6, and TB9, in Fig. 9. On the other hand, there is a large difference between the widths of the confidence intervals of N_H (0.4 to 1.0  ×  10²⁶ m⁻²). For the time-averaged contours, the intervals covered by the 68% confidence contours are ΔN_pow = 0.16 × 10⁵⁰ ph s⁻¹ keV⁻¹, ΔN_comt = 0.10 × 10⁵³ ph s⁻¹ keV⁻¹ and ΔN_H = 0.20 × 10²⁶ m⁻². Therefore, from Fig. 9 there does not appear to be any strong degeneracy between these model parameters, both from the time-averaged and TB spectra, and the apparent anti-correlation between the continuum and obscurer column density seen in Figs. 7 and 8 are due to a physical relation discussed in Sect. 5.3.

Finally, we carry out Monte Carlo spectral simulations to test whether there are any degeneracies between N_H and N_pow or N_comt. We do this by simulating spectra in SPEX using the best fit models (Table 3) fitted to the time-averaged spectrum of Obs 1 and the spectra of TB2, TB6, and TB9. We then refit the best fit models (all of the parameters in Table 3) to the simulated spectra (time-averaged, TB2, TB6, or TB9), repeating this process 50 times each. The results of these simulations are shown in Fig. 10, where we plot N_H against N_pow in the left side of each panel and N_H against N_comt in the right side. In addition, we also present the Pearson rank (r) and p-value in each panel, demonstrating the linear correlation coefficient and corresponding probabilities for the parameters.

Fig. 10. Monte Carlo spectral simulation results for the time-averaged model (top panel; red) and the models from TB2 (top-middle; blue), TB6 (bottom-middle; brown), and TB9 (bottom panel; cyan). The left and right sides of each panel display the results of NH against Npow and NH and Ncomt, respectively, after fitting the best fit models to the simulated spectra 50 times each. Also displayed in each panel are the Pearson rank (r) and p-values to demonstrate the linear correlation coefficient and corresponding probabilities for the parameters, used to determine parameter degeneracies.

Overall, there does not appear to be any correlation between N_H and N_comt for any TB or time-averaged simulations (right panels in Fig. 10), suggesting that there is no degeneracy between these two parameters (the Pearson rank values are between r = −0.23 and 0.20 and the p-values are larger than 0.05, suggesting the correlations are not significant). This is a similar result for N_H against N_pow in TB2, where there is no significant correlation between these parameters (left panels in Fig. 10). However, for the Obs 1 time-averaged, TB6, and TB9 simulated results, there appears to be some correlation between N_H and N_pow, where the Pearson rank values are r = 0.673, 0.523, and 0.755, respectively, and p < 1 × 10⁻⁵ (displayed as < 0.001 in Fig. 10 for visual purposes). These correlations in the left panels of Fig. 10 imply that there could be some degeneracy between N_H and N_pow. The results in Fig. 10 also suggests that if the obscurer was to change with the continuum, it would vary with N_pow (not N_comt). On the other hand, we conclude that there is no degeneracy between N_H and N_comt.

5.3. Cause of the observed variability

From Figs. 7 and 8, and the discussion in Sects. 5.1 and 5.2, there is a clear inverse correlation between the continuum (N_pow and N_comt) and the column density of the obscurer (N_H) in Obs 1, but no clear correlation between ξ and the continuum. We therefore carried out further tests to see whether the observed changes could be caused by the continuum, ionisation parameter, or the column density.

To do this, we firstly fixed the obscurer parameters in all of the Obs 1 TBs to the time-averaged values of N_H = 3.87 × 10²⁶ m⁻², log ξ = 0.44, and f_cov = 0.55 (see Table 3), and fitted the normalisations of POW and COMT, in addition to Γ_pow and s. The results for N_pow and N_comt can be seen in the left side of Fig. 11. We then fixed the continuum parameters to their time-averaged values in each TB, N_pow = 2.60 × 10⁵⁰ ph s⁻¹ keV⁻¹, N_comt = 1.89 × 10⁵³ ph s⁻¹ keV⁻¹, Γ_pow = 1.68 and s = 0.64, and fitted N_H, ξ and f_cov. The results can be seen in the right side of Fig. 11 for N_H and f_cov.

Fig. 11. Further tests to determine the cause of the observed variability. Left: Fixing the obscurer parameters to their time-averaged values in Obs 1 (Table 3), and fitting Npow (red circles) and Ncomt (blue triangles) in each TB. Right: Fixing the continuum parameters to their time-averaged values and fitting NH (magenta diamonds) and fcov (orange triangles). In both scenarios the changes in the parameters suggest that they could contribute to the observed variability. See Sect. 5.3 for details. The uncertainties on Npow (left) and NH (right) are relatively small that they are within the data points. The photon index and reflection scaling parameters from the continuum show the same trends as Fig. 7 when the obscurer parameters are fixed. The ionisation parameter of the obscurer does not change when the continuum parameters are fixed.

From the left side of Fig. 11, with the obscurer parameters fixed, both N_pow and N_comt increase over the course of Obs 1. This strong change suggests that the most of the variability is driven by continuum changes. The photon index and reflection scaling parameters show the same trend as their respective panels in Fig. 7. On the other hand, from the right side of Fig. 11, when the continuum parameters are fixed, both N_H and f_cov respond to changes in the spectrum independently to the continuum, which could also suggest that the obscurer parameters drive the variability. Nevertheless, given the large error bars on f_cov (Fig. 11) it seems that the covering fraction is consistent throughout Obs 1 and is therefore not changing, similar to Fig. 7.

Figure 11, therefore, is not very conclusive in determining the cause of the observed variability in Obs 1. There are two scenarios here: (1) the continuum varies and as a consequence, so does the obscurer; (2) the obscurer varies independently of the continuum changes. For scenario (1), the ionisation parameter ξ does not appear to vary with the continuum in Obs 1 (Fig. 7, which is supported further by the test in Fig. 11). However, from the left side of Fig. 11, it is clear that the continuum itself can drive the variability we observe, regardless of the obscurer. It is possible that ξ does change, but given the errors in Fig. 7, we cannot conclude one way or the other with full certainty; either possibility is plausible from these results. However, a likely explanation as to why the errors are so large is because we are fitting a multi-phased cloud, made up of different N_H, ξ and n_e regions, with a single model component. Alternatively for scenario (2), the obscurer varies independently of the continuum, which can be seen in the right side of Fig. 11, where N_H and f_cov have to change in order to obtain a good fit when the continuum parameters are fixed. However, it is not clear what causes the changes in the column density of the obscurer over the course of Obs 1.

One way in which the column density of the obscurer can change is by moving transversely across our LOS to the X-ray source. If the obscurer was clumpy, for example, then by crossing our LOS it would produce variations in N_H. The other two possible scenarios are that the obscurer expands, causing the column density to decrease, for an increase in covering fraction (but this requires a large expansion velocity of order 1120 km s⁻¹ – see Appendix D for the derivation), or alternatively, the material condenses to become visible dust (very unlikely as the ionising flux increases over the course of Obs 1). As the latter two explanations are less probable, we determine whether the obscurer could cross our LOS, and thus cause changes in N_H.

The location of the obscurer in Obs 1 was estimated to be between 0.012 and 0.055 pc (Paper III), assuming a spherical cloud geometry. Using ${\mathit{v}}_{\mathrm{cross}}=\sqrt{\frac{G{M}_{\mathit{BH}}}{r}}$, where G is the gravitational constant, M_BH = 5.96  ×  10⁶ M_⊙ is the black hole mass, and r is the distance of the obscurer from the SMBH, the crossing velocity of the obscurer, assuming it follows a Keplerian orbit (Paper III), is between 680 and 1470 km s⁻¹, depending on which distance is used. For this range in velocities, the obscurer would travel a distance of 6.80 × 10¹⁰ − 1.51 × 10¹¹ m (8 − 17 R_g) during Obs 1 (∼103 ks). Assuming that the radius of the X-ray corona is between 7 and 15 R_g (e.g. Chainakun et al. 2019; a fiducial value of 10 R_g was adopted in NGC 3783 by Mehdipour et al. 2017), then the diameter of the X-ray source would be within the range of 1.24 − 2.65  ×  10¹¹ m. Here, there is clearly an overlap in the size of the corona and the crossing distance travelled by the obscurer during Obs 1. Therefore, changes in the obscurer column density for Obs 1 (Fig. 7) could be explained with the obscurer moving transversely across our LOS towards the X-ray corona.

In comparison, the length of time between observations was 1810 ks (three weeks), such that the distance travelled by the obscurer, assuming v_cross = 680 − 1470 km s⁻¹, is 1.23 − 2.65 × 10¹² m (140 − 300 R_g). This is larger than the size of the ionising source, suggesting that the increase in N_H between Obs 1 and Obs 2 could be a result of the same obscurer moving transversely into our LOS for the two observations. This would explain the large drop in flux in the spectra and light curves between Obs 1 and Obs 2. However, the more likely argument to explain the difference in obscuration and continuum properties is that we are viewing two completely different components in the two observations. The main evidence is due to the ionisation parameter in Obs 2 being larger compared to Obs 1, even though the flux decreases greatly, as discussed in Paper III.

Therefore, from these arguments and different tests (in Sects. 5.2 and 5.3), we conclude that the observed variability in Obs 1 is likely to be driven by the continuum, reflecting changes in the intrinsic brightness of the central engine in NGC 3227. Some explanations for the continuum driving the observed changes could include: (1) changes in the accretion rate of NGC 3227; (2) fluctuations in the heating and cooling rates in the corona (e.g. Ballantyne & Xiang 2020, studied how a warm corona responds to the accretion disk, in addition to the heating and cooling effects); (3) variations in the corona size or location from the black hole (this was an argument to explain the observed changes in NGC 4151; Couto et al. 2016). However, we cannot rule out scenario 2), whereby the possible changes could be caused by N_H as a result of the obscurer moving transversely over the X-ray source within our LOS during the observation time; the crossing distance is comparable to the X-ray corona size.

One final possible solution to explain the trend between the continuum and obscurer column density could be established if the obscurer is a wind. In NGC 1365, Connolly et al. (2014) found that the absorbing column density decreased for an increase in luminosity, with no change in the ionisation state of the material; there is a similar trend here and can be seen in Figs. 7 and 8 for NGC 3227. In NGC 1365, this trend was interpreted through a disk wind model, where if the X-ray luminosity increased, then the launching radius of the wind became larger, causing the amount of material in the LOS to decrease (Connolly et al. 2014). Therefore, if we view the X-ray source through the inner edge of the obscurer within NGC 3227, and the continuum flux increases, then the obscurer would move outwards from within our LOS, causing the NH to decrease. A schematic can be seen in Fig. 5 of Connolly et al. (2015). However, this solution is speculative in NGC 3227, because there is no evidence that the obscurer is a wind due to the lack of UV absorption lines (Paper I; Paper III), while in NGC 5548 (Kaastra et al. 2014) and NGC 3783 (Mehdipour et al. 2017; Kriss et al. 2019) the outflow velocities of the obscurers were measured with UV absorption features. In addition, we note that a leading model for X-ray continuum variability across a range of timescales is the model of propagating fluctuations in the local mass accretion rate (e.g. Lyubarskii 1997; Arévalo & Uttley 2006; Ingram & van der Klis 2013), that could also explain the changes we observe in NGC 3227 during the obscuration event at the end of 2019.

5.4. Excess variance

Finally, we further explore the variability of each observation using the variance and fractional root mean square (rms) variability amplitude (F_var; Vaughan et al. 2003; Ponti et al. 2004), similar to Lobban et al. (2020) for NGC 3227 and Mehdipour et al. (2016b) for NGC 5548. We calculate F_var following the steps from Vaughan et al. (2003), where the difference between the source variance (S²) and the measured errors ($\overline{{\sigma }_{\mathrm{err}}^2}$) give rise to the excess variance ${\sigma }_{\text{XS}}^2$. Here, $S^2=\sum _{i=1}^N\frac{{(x_i-\overline{x})}^2}{(N-1)}$ and $\overline{{\sigma }_{\mathrm{err}}^2}=\frac{1}{N}\sum _{i=1}^N{\sigma }_{\mathrm{err},\mathrm{i}}^2$, where x_i is the measured value (flux), $\overline{x}$ is the mean, N is the number of measurements, and ${\sigma }_{\text{err},\text{i}}^2$ is the uncertainty, or mean square error, on each measurement. The fractional rms variability amplitude is related to the excess variance via

$$\begin{array}{c}\hfill F_{\mathrm{var}}=\sqrt{\frac{{\sigma }_{\mathrm{XS}}^2}{{\overline{x}}^2}}=\sqrt{\frac{S^2-\overline{{\sigma }_{\mathrm{err}}^2}}{{\overline{x}}^2}}.\end{array}$$(1)

To compute the uncertainty on F_var, we used Eq. (B2) from Vaughan et al. (2003).

The top panel of Fig. 12 shows the fractional rms variance spectrum for the full duration of Obs 1 (103 ks) and Obs 2 (50 ks), calculated using Eq. (1), with the EPIC-PN light curves binned at 1000 s for both observations. For Obs 1 (blue circles) the energies are binned logarithmically, such that the 20 data points are evenly spaced out. For Obs 2 (red diamonds in top panel of Fig. 12), we have used large energy bins, such that there is sufficient S/N to prevent negative excess variance (${\sigma }_{\text{XS}}^2$) values resulting from the much weaker signal in Obs 2, compared to Obs 1.

Fig. 12. Fractional root mean square (rms) variability amplitude (Fvar) spectra as a function of energy. Top panel: Obs 1 (blue circles) and Obs 2 (red diamonds), where the light curves were binned at 1000 s. For Obs 1, the energies are binned logarithmically such that the 20 data points are evenly spaced. For Obs 2, large energy bins are used to remove negative excess variance (${\sigma }_{\text{XS}}^2$) values caused by the significant drop in flux compared to Obs 1; the bin sizes are displayed by the x-axis error bars. The purple Fvar spectrum shows the results from averaging the two halves of Obs 1. See Sect. 5.4 for details. Middle panel: Fractional rms variability amplitude spectra for the best fit spectral model (red line) and the three main model components: POW (magenta line), COMT (purple line), and the obscurer (PION; green line). The blue line shows the spectrum from the light curves (see top panel). From this plot, it appears that the COMT component drives the fractional rms variability below 1 keV, while the obscurer drives the Fvar spectrum between 1 and 5 keV; above 5 keV POW dominates. Bottom panel: Fractional rms variability amplitude spectra comparing different timing bins for the EPIC-PN light curves from Obs 1. The frequencies are 10−2 – 5 × 10−4 Hz (Δt = 100, 500, 1000 and 2000 s; blue), 10−4 Hz (Δt = 10 ks; magenta), 5 × 10−5 Hz (Δt = 20 ks; red), and 2 × 10−5 Hz (Δt = 50 ks; green). The shape of each curve stays consistent between frequency bands, suggesting that the model components vary for all light curve timing bins.

The comparison between Obs 1 and Obs 2 does not take into account red-noise variability, as the duration of Obs 2 is half of Obs 1, and hence would infer lower variability for the Obs 2. Therefore, to test this, we split Obs 1 into two halves and found F_var in each energy bin (same 20 bins in Fig. 12), averaging the values between the two halves. This can be seen by the purple F_var spectrum in the top panel of Fig. 12. The result produces a fractional rms variability spectrum that is slightly above Obs 2 when judged on the same timescales (although the difference is less than implied by the original comparison), and the spectrum is more constant across the X-ray energy range compared to the blue line.

To determine what causes the features in the F_var curve of Obs 1 (blue line in the top panel of Fig. 12), we obtained the F_var values for the best fit model and the POW, COMT, and obscurer components (similar to Mehdipour et al. 2016b; Cappi et al. 2016, in their analysis of NGC 5548). This is displayed in the middle panel of Fig. 12. To obtain these F_var spectra, we binned each TB spectrum (Fig. 4) into 20 energy bins and obtained the F_var values for each energy bin from all ten TB model spectra. We did this for the best fit model fitted to each TB, which is the red line in Fig. 12. The shape of the fractional rms variability amplitude spectrum for the model follows a similar shape as the F_var spectrum from the EPIC-PN light curves (blue line in the middle panel). The F_var spectra for the three main components in the best fit models are also displayed in the middle panel – POW (magenta line), COMT (purple line), and the obscurer PION (green line). In order to obtain the F_var values for the obscurer, we took the difference between the absorbed best fit model (red line in Fig. 12) and the model when N_H = 0 m⁻². The uncertainties on the measured values for the TBs (${\sigma }_{\text{err},\text{i}}^2$) are the standard error on the mean for each energy bin (similar to Lobban et al. 2020).

From plotting the F_var spectra of the model components, we can see that the COMT component accounts for the variability at energies below 1 keV, while the obscurer causes changes between 1 and 5 keV; POW dominates above 5 keV. This could give additional evidence that the obscurer contributes significantly to the observed variability, and therefore this scenario should not be ruled out compared to the continuum, as discussed in Sect. 5.3. Moreover, Fig. 12 further confirms the results from the C-statistic summation test in Sect. 5.2, which shows that N_comt and N_H are statistically significant parameters in the best fit model. In the middle panel of Fig. 12, the curves for the COMT component (purple line) and the obscurer (green line) follow the best fit model F_var spectrum (red line), suggesting that they likely contribute significantly to the fractional rms variability values within the specific energy ranges. Finally, there is a drop in F_var at 6 keV which likely corresponds to the Fe K emission energy band, from material far from the SMBH, and therefore varies less compared to the underlying continuum. This was also found in the 2016 observations by Lobban et al. (2020).

In addition to obtaining the fractional rms variability amplitude from the spectral model components (middle panel of Fig. 12), we also explored how different timing bins for the EPIC-PN light curves affected the fractional rms variability amplitudes in Obs 1 (similar to, e.g. Arévalo & Markowitz 2014; Lobban et al. 2020, for NGC 3227). For each energy bin of the EPIC-PN light curves from Obs 1 we constructed F_var spectra for seven different timing bin lengths, Δt= 0.1, 0.5, 1, 2, 10, 20, and 50 ks. For the four smallest timing bins, we found that the F_var values are similar and therefore average these values together. We convert these timing bins to frequency bands of 10⁻² − 5 × 10⁻⁴ Hz for the averaged values, and 10⁻⁴, 5 × 10⁻⁵, and 2 × 10⁻⁵, for 10, 20 and 50 ks, respectively (Arévalo & Markowitz 2014; Lobban et al. 2020). The result is displayed in the bottom panel of Fig. 12. For all frequencies, the energy bins correlate with each other. In other words, if there is an increase in F_var for a particular energy range, such as between 1 and 5 keV (caused by the obscurer; middle panel of Fig. 12), all frequency fractional variability amplitudes follow this increase. This result suggests that all model components that vary, vary on all timescales (Arévalo & Markowitz 2014). As the frequency increases, the shape of the curve stays the same but the whole curve rises in F_var, as shown in the bottom panel of Fig. 12. This result is different from the analysis by Lobban et al. (2020) on NGC 3227 in 2016, whereby they found that as the frequency increased, the steepness of the curve decreased, until the highest frequency band displayed little change in F_var over the full energy range (top panel in their Fig. 7).

6. Summary and conclusions

NGC 3227 has undergone many occultation events (Lamer et al. 2003; Markowitz et al. 2014; Beuchert et al. 2015; Turner et al. 2018) lasting a variety of time scales. At the end of 2019, we observed this AGN to undergo obscuration again, detected by our Swift monitoring programme, where follow-up observations from XMM-Newton, NuSTAR and HST were obtained. In Paper I, we studied the unobscured and obscured SEDs from 2016 and 2019, respectively, in order to probe the properties and origin of the obscuring material in the sequential papers in the series. Paper II analysed the WAs and past occultation periods using the broadband continuum model, before the 2019 observations were studied in great detail in Paper III.

In this paper, using the initial model parameter values from Paper III, we investigated the possible origins of variability during the two observations of NGC 3227 taken at the end of 2019. In Obs 1, we found an anti-correlation between N_H and the continuum normalisations, N_pow and N_comt (Figs. 7 and 8). We conclude that the observed variation in Obs 1 is likely to be driven by the continuum. However, we cannot rule out the possible changes could be caused by N_H if the obscurer moved transversely over the X-ray source within our LOS during the observation time.

For Obs 2, the count rate was significantly lower compared to Obs 1, due to the combination of a high column density and covering fraction (f_cov) in the obscurer and strong flux decrease in the continuum. The lack of strong variability suggests that the continuum driving mechanism is not as significant compared to Obs 1. Between Obs 1 and Obs 2, the time-averaged results for f_cov and N_H increased from 0.55 to 0.66, and from 3.87 × 10²⁶ m⁻² to 9.93 × 10²⁶ m⁻², respectively. Comparatively, N_pow and N_comp decreased from 2.60 to 0.73 × 10⁵⁰ ph s⁻¹ keV⁻¹at 1 keV, and 1.89 to 0.16 × 10⁵³ ph s⁻¹ keV⁻¹, respectively.

In conclusion, we argue that the continuum is likely to drive the strong variability seen in NGC 3227 during Obs 1, but we cannot rule out the possible effects from N_H of the obscurer. There is no apparent change in the ionisation parameter of the obscurer in Obs 1, but given the large uncertainties, we cannot confirm this with certainty. A likely explanation here is that the obscurer is multi-phased, which we fit with only one component. Obs 2, on the other had, shows little change over the course of the observation, and the large increase in N_H and f_cov, causing the significant drop in X-ray flux, is a result of a different component in our LOS compared to the one in Obs 1.

Acknowledgments

This work is based on observations obtained with XMM-Newton, an ESA science mission with instruments and contributions directly funded by ESA Member States and NASA. We would like to thank the referee for their constructive and helpful feedback, and thank you to C. Done for suggesting a relevant paper to help explain the observed relationships between obscurer and continuum. S.G.W. acknowledges the support of a PhD studentship awarded by the UK Science & Technology Facilities Council (STFC). G.P. acknowledges funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme. B.D.M. acknowledges support from Ramón y Cajal Fellowship RYC2018-025950-I (grant agreement No 865637). S.B. acknowledges financial support from ASI under grants ASIINAFI/037/12/0 and n. 2017-14-H.O and from PRIN MIUR project ‘Black Hole winds and the Baryon Life Cycle of Galaxies: the stone-guest at the galaxy evolution supper’, contract no. 2017PH3WAT. P.O.P. acknowledges financial support from the French High Energy National Programme (PNHE) and from the French space agency (CNES). D.J.W. acknowledges support from STFC in the form of an Ernest Rutherford Fellowship (ST/N004027/1).

References

Appendix A: Alternative models

In Paper III, we examined three different models applied to the November 2019 observation, and four models to the December 2019 observation, of NGC 3227 (see Tables 2 and 3, and Sect. 3, in Paper III for details). In the modelling carried out in this paper (Sects. 4 and 5), we assumed that the neutral reflector was located close enough to the SMBH such that changes in the reflected continuum in 2019 were caused by changes in the intrinsic power law (Gandhi et al. 2015). However, in addition to the modelling in Sect. 4, we also fitted the M2 model from Table 2 in Paper III, coupling the REFL component normalisation and photon index in 2019 to the 2016 power law values in both observations. This model produced a better, albeit not significantly different, C-statistic compared to M1, but the parameter values were consistent between the models; a similar result was found in Paper III. As a result of this fitting, we found that these changes on the REFL component do not affect the properties of the obscurer. Therefore, no actions to update the current model were taken, and the results from Sects. 4 and 5 explain the properties of the obscurer and variations in the observed continuum for Obs 1 and Obs 2.

The third model applied to the 2019 observations in Paper III (M3 in Table 2 and M4 in Table 3) was the two obscurer model, fitted here with two PION components. Similarly to the modelling in Sect. 4 of this paper, the column densities, ionisation parameters and covering fractions of the obscurer components were fitted, along with the continuum parameters. After fitting the time-averaged spectra with the two obscurer model, we found large differences with the parameters in the current model from Fig. 7 and Table 3.

As a result of this, therefore, we fitted the two obscurer time-averaged best fit models to each TB spectrum in Obs 1 and Obs 2. For Obs 2, as expected, there is little change in the parameters between TBs and they are consistent with the time-averaged values, similar to Fig. 7. For Obs 1, on the other hand, there are changes in parameter values between TBs, such as N_pow and N_comt, which increase during the observation, similar to the model in Sect. 5 and Fig. 7. However, the difference between the one and two obscurer models is clearer when comparing the obscurer parameters. For both obscurers in the two obscurer model, the column densities stay fairly consistent between TBs, and compared to the time-average values, as displayed in the left side of Fig. A.1. However, it is the covering fractions that appear to change in Obs 1: for obscurer 1, f_cov, 1 decreases over the duration of the observation, while f_cov, 2 increases. This comparison is displayed in the right panel of Fig. A.1. These results suggest that the covering fraction could be related to the continuum changes, rather than the column density; the opposite of what we initially discussed in Sect. 5.3. On the other hand, when the TB errors are taken into account, the covering fraction values overlap, suggesting that the covering fraction correlations are not significant throughout Obs 1. Similarly with the model in Sects. 4 and 5, the ionisation parameters of both obscurer components do not show a clear trend with over the Obs 1 and were difficult to constrain. As a result, we do no show this plot.

Fig. A.1. Comparing the parameters of the two obscurer model; details are discussed in Appendix A. Left: Column densities of the two obscurer components, blue for obscurer component 1 (NH, 1) and red for component 2 (NH, 2). Right: Covering fractions of the two obscurer components, orange for obscurer 1 (fcov, 1) and magenta for obscurer 2 (fcov, 2). The dashed lines in each panel show the time-averaged best fit parameter values, with their uncertainties displayed as the shaded areas.

Further investigations were then implemented, similar to Fig. 11 in Sect. 5.3 - we fixed either the obscurer or continuum parameters to their time-averaged values and fitted the other. From this test, we find that both the column densities and covering fractions in both obscurers decrease, however the first obscurer component changes are more significant. These parameters show an anti-correlation with the continuum parameters N_pow and N_comt, similar to the result found in Figs. 7 and 11; albeit this time the covering fractions also vary. This test implies that both obscurer parameters can vary in accordance with the continuum, and are therefore valid to explain the change in the obscurer properties. Therefore, from this argument, we conclude that the observed variability in Obs 1 is likely caused by the intrinsic properties of the continuum, but changes caused by the obscurer cannot be ruled out as we see N_H and f_cov vary independently of the continuum. This is the same conclusion we found for the one obscurer component model in Sect. 6.

In summary, the model fitted to the EPIC-PN spectra in Obs 1 and Obs 2 (Sects. 4 and 5) explains the properties of the continuum and obscurer of NGC 3227 in 2019. There is no significant consequence on the obscurer parameters when coupling either the 2016 or 2019 power law to the reflection component, and similar trends are seen between the one and two obscurer models. Therefore, we argue that the model in this paper is the correct one, compared to the alternative models discussed in this Appendix.

Appendix B: KNAK corrections

This section explains how we corrected for the cross-calibration issue between the RGS and EPIC-PN instruments, presented in Sect. 4.2. We first obtained the ratio between the data (D_PN) and the initial model (M_PN) from Paper III for the EPIC-PN spectra ($R_{\mathrm{PN}}=\frac{D_{\mathrm{PN}}}{M_{\mathrm{PN}}}$; from Fig. 5), between 6 and 37 Å (as this corresponds to the RGS band, shown in Fig. B.1). Here we assumed that the model from Paper III fitted the RGS data very well, such that the ratio between the data (D_RGS) and model (M_RGS) was approximately one ($R_{\mathrm{RGS}}=\frac{D_{\mathrm{RGS}}}{M_{\mathrm{RGS}}}\sim 1$). Next, we applied a KNAK component to the model (similar to Porquet et al. 2004; Detmers et al. 2010), which is a piecewise broken power law. KNAK splits the model into continuous segments, multiplying the model by a factor T_N, set to R_PN taken from Fig. B.1 (top panel), for a given wavelength (λ_N), where N are the segment number (the maximum is nine per KNAK component). KNAK only works for a wavelength grid, hence Fig. B.1 is in units of Angstroms rather than in keV (see the SPEX manual for more details³). The KNAK model describes the instrumental effects and not the physics of NGC 3227, therefore when applying the KNAK component we kept the initial model parameters fixed. In order to be consistent outside the selected range, the lowest wavelength (λ₁) was fixed to 0.1 Å and T₁ was coupled to T₂. For the opposite end, λ₉ was fixed at 100 Å and T₉ was coupled to T₈. For the remaining seven segments, we chose R_PN and its respective wavelength by eye from Fig. B.1. Once the KNAK parameters (T_N and λ_N) were fitted, we kept them fixed throughout the rest of the EPIC-PN modelling (Sect. 4.2). We repeated the above for Obs 2, as shown in the right side of Fig. B.1.

Fig. B.1. Steps taken to correct for the poor initial model fit caused by the RGS and EPIC-PN cross calibration issue (see Appendix B for details). Top: The ratio between the data and model was used in the KNAK model (λN and TN) for the correction. Bottom: The EPIC-PN spectrum between 6 - 36 Å with the two models on top. The KNAK corrected model was fitted to the data here. In Obs 1 (left) green is the initial model from Paper III and blue is the corrected model; in Obs 2 (right) pink and purple are the initial and corrected models, respectively. The KNAK component only works for a wavelength grid, and therefore the units of this figure are in Angstroms rather than keV.

The top panels of Fig. B.1 shows the ratio between the data and the model for Obs 1 (left) and Obs 2 (right). The ratios were taken from Fig. 5, where the initial model was taken directly from Paper III, and the corrected model was achieved after the KNAK corrections. The bottom panels in Fig. B.1 display the EPIC-PN spectrum from each observation in the RGS wavelength range. Again, the initial and corrected models are shown in the bottom panels of Fig. B.1, where the C-statistic and d.o.f (in parentheses) are shown. After making the corrections in each observation, the KNAK parameters were fixed for the remainder of the spectral fitting (Sect. 4.2).

In both observations (Fig. B.1), there are significant differences between the data and the model at longer wavelengths. However, for Obs 1 (left side), there are strong residuals between 7 and 12 Å. The difference between the data and model in the initial fits for each observation, although taken less than a month a part, could be due to how well the initial model from Paper III was constrained, or whether there are some unresolved features present in the EPIC-PN data that were fitted with RGS. Furthermore, as Obs 1 has a larger flux compared to Obs 2, any significant differences between the RGS and EPIC-PN spectra will be more obvious, and hence why there are more residuals in Obs 1 compared to Obs 2. This is why we carried out the KNAK corrections.

Appendix C: Additional plots

Figure C.1 displays additional correlations between the fitted parameters in the models for Obs 1 and Obs 2, further to the parameter comparisons in Fig. 8 (Sect. 5.2).

Fig. C.1. Comparing the relations between parameters in the model, in addition to the parameters discussed in Sect. 5.2 and Fig. 8. Ionisation parameter (top left) and reflection scaling parameter (top right) against the 0.3 - 10 keV count rates. The remaining panels display the correlations between the different parameters in the model. In each panel, the Pearson rank (r) and p-value is shown to display the significance of each correlation.

In addition to Fig. 8, we found some, albeit less conclusive, trends between parameters which we present in Fig. C.1 in Appendix C. These include N_pow, N_comt, Γ, and f_cov against s, N_H against f_cov, and s compared to the 0.3 - 10 keV count rates. Furthermore, Fig. C.1 displays ξ against the 0.3 - 10 keV count rates, and N_comt against ξ. The former comparison shows no strong correlation, as discussed in Sect. 5.1, and the latter comparison does not show a reliable trend, with some outliers. In each panel of Fig. C, the Pearson rank and p-values are shown to display how well correlated the parameters are, and the significance of each trend. In the majority of the panels in Fig. C, the correlations are not that significant as p = 0.001 − 0.02, however some panels in Fig. C.1 have values that are lower. In some panels, the absolute values of r values are greater than 0.7 (both in the positive and negative directions) implying that these parameters are somewhat well correlated or anti-correlated with each other, however less significantly compared to the parameters in Fig. 8. The panels comparing log ξ with either the count rates or N_comt values have r = −0.5 to −0.6, suggesting that the ionisation state does not correlate or change with the continuum; similar to the bottom panel in Fig. 7.

Appendix D: Expansion velocity derivation

Following on from Sect. 5.3, one way in which the column density of the obscurer could change is if the cloud expanded, causing N_H to decrease. Assuming that the obscurer cloud is spherical, it will have a mass given by

$$\begin{array}{c}\hfill M=4\pi r^2{m}_p\mu N_{\mathrm{H}}{f}_{\mathit{cov}},\end{array}$$(D.1)

where m_p is the proton mass, μ is the mean atomic mass per proton (∼1.4 for solar values), r is the radius of the cloud, and N_H and f_cov are the column density and covering fraction of the obscurer from the modelling in Sect. 4. If the obscuring cloud does expand during Obs 1, and we assume the mass of the obscurer stays constant, we can equate the masses from TB1 and TB10 as follows

$$\begin{array}{c}\hfill M=4\pi r^2{m}_p\mu N_{H,1}{f}_{cov,1}=4\pi R^2{m}_p\mu N_{H,10}{f}_{cov,10},\end{array}$$(D.2)

where r and R are the radii of the obscurer in TB1 and TB10, respectively.

A quick rearrangement gives R² = Ar² where $A=\frac{N_{H,1}{f}_{cov,1}}{N_{H,10}{f}_{cov,10}}$. Substituting in the column density and covering fraction values for TB1 and TB10 from Table 3, the result yields $R=\sqrt{2.25}r$.

Here, we set $r=\frac{N_{H,1}}{n_{\mathrm{H}}}$ for TB1, where n_H is the hydrogen number density. From Paper III (Table 4), n_H = 0.1 − 2.7 × 10¹⁵ m⁻³ for a spherical cloud. Taking the average of n_H gives an obscurer radius in TB1 of r = 2.24 × 10¹¹ m, and therefore yields a radius after expansion in TB10 of R = 3.36 × 10¹¹ m. So the expansion distance travelled during Obs 1 would therefore be D = R − r = 1.12 × 10¹¹ m.

Finally, the expansion velocity can be estimated if we take the time to be t = 103 ks, and as such v_exp = D/t ∼ 1120 km s⁻¹. Therefore, the obscurer would have to expand with a velocity of around 1120 km s⁻¹, which is in the middle of the estimated crossing velocity (v_cross) of the obscurer (between 680 and 1470 km s⁻¹), to explain the observed changes in the column density between TB1 and TB10. This expansion scenario to explain the observed decrease in N_H for Obs 1 is rather unlikely as the expansion speed is similar to the travelling velocity.

All Tables

All Figures

Fig. 1. Light curves for Obs 1 (left) and Obs 2 (right). Top panels: Three different energy bands: 0.3–10 keV (red), 0.3–2 keV (green), and 2–10 keV (blue). Also displayed are the timing bins (TB) from which we extracted EPIC-PN spectra for our analysis (see Sect. 3.2 for details). The light curves were binned at 100 s. The errors on the energy bin fluxes are shown by the shaded regions. Bottom panels: Ratio light curves between the 0.3–2 keV (green) and 2–10 keV (blue) energy bands in the top panel divided by the full 0.3–10 keV light curve (red line). The Obs 1 ratios are on the left hand side and the Obs 2 ratios are on the right.

In the text

	Fig. 2. Hardness ratio between the 2–10 keV and 0.3–2 keV energy bands, defined as (H–S)/(H+S), for Obs 1 (left) and Obs 2 (right), binned at 100 s. The errors on the hardness ratios are shown by the shaded regions. The red dashed lines show the average hardness ratio (HRavg) across each observation. For Obs 1 HRavg = −0.115, and for Obs 2 HRavg = −0.033.
In the text

Fig. 3. Top: Relation between the averaged 0.3–2 keV (soft) and 2–10 keV (hard) count rates for each timing bin. Obs 1 (hollow circles) and Obs 2 (diamonds) are compared with the time-averaged values shown in red. Bottom: Hardness ratio (H–S/H+S) as a function of the average 0.3–10 keV count rates. NGC 3227 shows a softer-when-brighter trend on timescales of weeks as Obs 1 has a larger flux but lower hardness ratio compared to Obs 2.

In the text

	Fig. 4. EPIC-PN spectra for Obs 1 (left) and Obs 2 (right). The time-averaged spectra are shown in red and the TB numbers are presented in the legends by their respective colours.
In the text

Fig. 5. Initial and KNAK corrected models (see Appendix B) for Obs 1 (left) and Obs 2 (right). The initial models were from Paper III (green and pink), while the corrected models are in blue and purple, respectively. The C-statistic values are displayed, along with their respective degrees of freedom in the parentheses. The bottom panels show the residuals between the data and model. Appendix B explains the steps required to correct the models, shown in Fig. B.1.

In the text

	Fig. 6. Final best fit models fitted to the time-averaged spectra: Obs 1 grey crosses and blue model; Obs 2 black crosses and purple model. Also shown are the C-statistics and degrees of freedom. The best fit parameter values are displayed in Table 3. The bottom panel shows the residuals between the data and model for Obs 1 (blue) and Obs 2 (purple).
In the text

Fig. 7. Model parameter comparisons for each timing bin (TB) in Obs 1 (hollow circles) and Obs 2 (diamonds). The red horizontal dashed and solid lines represent the time-averaged parameter values for Obs 1 and Obs 2, respectively, and the shaded areas show the uncertainties. The TB colours are consistent with their respective spectra in Fig. 4. The model parameter is shown on the y-axis of each panel. The light purple lines (dashed for Obs 1 and solid for Obs 2) connect the TB values for each parameter in the two observations.

In the text

Fig. 8. Correlations between model parameters. Left: The normalisations of the power law (Npow) and Comptonisation component (Ncomt) are correlated over the two observations (top); Npow against the column density of the obscurer (NH) (middle); and Ncomt against NH (bottom). Right: Correlations between Npow, Ncomt and NH against with the 0.3–10 keV count rate for top, middle and bottom, respectively. In each panel, the Pearson rank (r) and p-value is shown to display the significance of each correlation.

In the text

Fig. 9. Contour plots of NH (1026 m−1) against Npow (1050 ph s−1 keV−1; top panels) and NH (1026 m−1) against Ncomt (1053 ph s−1 keV−1; bottom panels). From left to right, the panels show the contours for the time-averaged model, TB2, TB6, and TB9, respectively, in Obs 1. The levels of confidence presented are at 68 % (ΔC = 2.30; black solid line), 90 % (ΔC = 4.61; red dashed line), and 95 % (ΔC = 6.17; blue dotted line). The black crosses show the parameter values (but not the errors) from the best fit models in Table 3.

In the text

Fig. 10. Monte Carlo spectral simulation results for the time-averaged model (top panel; red) and the models from TB2 (top-middle; blue), TB6 (bottom-middle; brown), and TB9 (bottom panel; cyan). The left and right sides of each panel display the results of NH against Npow and NH and Ncomt, respectively, after fitting the best fit models to the simulated spectra 50 times each. Also displayed in each panel are the Pearson rank (r) and p-values to demonstrate the linear correlation coefficient and corresponding probabilities for the parameters, used to determine parameter degeneracies.

In the text

Fig. 11. Further tests to determine the cause of the observed variability. Left: Fixing the obscurer parameters to their time-averaged values in Obs 1 (Table 3), and fitting Npow (red circles) and Ncomt (blue triangles) in each TB. Right: Fixing the continuum parameters to their time-averaged values and fitting NH (magenta diamonds) and fcov (orange triangles). In both scenarios the changes in the parameters suggest that they could contribute to the observed variability. See Sect. 5.3 for details. The uncertainties on Npow (left) and NH (right) are relatively small that they are within the data points. The photon index and reflection scaling parameters from the continuum show the same trends as Fig. 7 when the obscurer parameters are fixed. The ionisation parameter of the obscurer does not change when the continuum parameters are fixed.

In the text

Fig. 12. Fractional root mean square (rms) variability amplitude (Fvar) spectra as a function of energy. Top panel: Obs 1 (blue circles) and Obs 2 (red diamonds), where the light curves were binned at 1000 s. For Obs 1, the energies are binned logarithmically such that the 20 data points are evenly spaced. For Obs 2, large energy bins are used to remove negative excess variance (${\sigma }_{\text{XS}}^2$) values caused by the significant drop in flux compared to Obs 1; the bin sizes are displayed by the x-axis error bars. The purple Fvar spectrum shows the results from averaging the two halves of Obs 1. See Sect. 5.4 for details. Middle panel: Fractional rms variability amplitude spectra for the best fit spectral model (red line) and the three main model components: POW (magenta line), COMT (purple line), and the obscurer (PION; green line). The blue line shows the spectrum from the light curves (see top panel). From this plot, it appears that the COMT component drives the fractional rms variability below 1 keV, while the obscurer drives the Fvar spectrum between 1 and 5 keV; above 5 keV POW dominates. Bottom panel: Fractional rms variability amplitude spectra comparing different timing bins for the EPIC-PN light curves from Obs 1. The frequencies are 10−2 – 5 × 10−4 Hz (Δt = 100, 500, 1000 and 2000 s; blue), 10−4 Hz (Δt = 10 ks; magenta), 5 × 10−5 Hz (Δt = 20 ks; red), and 2 × 10−5 Hz (Δt = 50 ks; green). The shape of each curve stays consistent between frequency bands, suggesting that the model components vary for all light curve timing bins.

In the text

Fig. A.1. Comparing the parameters of the two obscurer model; details are discussed in Appendix A. Left: Column densities of the two obscurer components, blue for obscurer component 1 (NH, 1) and red for component 2 (NH, 2). Right: Covering fractions of the two obscurer components, orange for obscurer 1 (fcov, 1) and magenta for obscurer 2 (fcov, 2). The dashed lines in each panel show the time-averaged best fit parameter values, with their uncertainties displayed as the shaded areas.

In the text

Fig. B.1. Steps taken to correct for the poor initial model fit caused by the RGS and EPIC-PN cross calibration issue (see Appendix B for details). Top: The ratio between the data and model was used in the KNAK model (λN and TN) for the correction. Bottom: The EPIC-PN spectrum between 6 - 36 Å with the two models on top. The KNAK corrected model was fitted to the data here. In Obs 1 (left) green is the initial model from Paper III and blue is the corrected model; in Obs 2 (right) pink and purple are the initial and corrected models, respectively. The KNAK component only works for a wavelength grid, and therefore the units of this figure are in Angstroms rather than keV.

In the text

Fig. C.1. Comparing the relations between parameters in the model, in addition to the parameters discussed in Sect. 5.2 and Fig. 8. Ionisation parameter (top left) and reflection scaling parameter (top right) against the 0.3 - 10 keV count rates. The remaining panels display the correlations between the different parameters in the model. In each panel, the Pearson rank (r) and p-value is shown to display the significance of each correlation.

In the text