Transient obscuration event captured in NGC 3227 IV. Origin of the obscuring cloud variability

Obscuration events in type I active galactic nuclei (AGN) have been detected more frequently in recent years. The strong flux decrease in the soft X-ray band between observations has been caused by clouds with large column densities transiting our line-of-sight (LOS) and covering the central AGN. Another event has been captured in NGC 3227 at the end of 2019. We aim to determine the nature of the observed spectral variability in 2019 obscuration event. We split the two XMM-Newton observations from 2019 into timing bins of length $\sim$ 10 ks. We used the SPEX code to analyse the 0.35-10 keV EPIC-PN spectra of each timing bin. In the first observation (Obs 1), there is a strong anti-correlation between the column density ($N_H$) of the obscurer and the continuum normalisations of the X-ray power-law and soft Comptonisation components ($N_{pow}$ and $N_{comt}$, respectively). The powerlaw continuum models the hard X-rays produced by the corona, and the Comptonisation component models the soft X-ray excess and emission from the accretion disk. Through further testing we conclude that the continuum is likely to drive the observed variability, but we cannot rule out a possible contribution from NH of the obscurer if it fully transverses across the ionising source within our LOS during the observation. The ionisation parameter ($\xi$) of the obscurer is not easily constrained, and therefore it is not clear whether it varies in response to changes in ionising continuum. The second observation (Obs 2) displays a significantly lower count rate due to the combination of a high NH and covering fraction of the obscurer, and a lower continuum flux. The observed variability seen during the obscuration event of NGC 3227 in 2019 is likely driven by the continuum, but the obscurer varies at the same time, making it difficult to distinguish between the two possibilities with full certainty.


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. 2013Parker 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 Article number, page 1 of 19 arXiv:2303.04717v2 [astro-ph.GA] 12 Mar 2023 and these events are more common than first thought . For example, multiple obscuration events were seen in NGC 3783 between 1993 and 2016 , 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 Comptonthick 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 24 − 10 26 m -2 ) and outflow velocities (v out,WA = few 100 − 1000 km s -1 ; 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 26 − 10 27 m -2 (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 Xray and UV bands, particularly during periods of high flux (Markowitz et al. 2009;Arévalo et al. 2014), with a 'steeperwhen-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 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 -1 ) with increasing covering fraction (C f ; from 0.7 to 0.9) and column density (N H ; from 5 to 18 × 10 26 m -2 ), 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 26 m -2 ; log ξ = 2; v out ∼ −800 km s -1 ). 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 1 ranging from log ξ = −1 to 3, and outflow velocities in the range v out = −100 to −1300 km s -1 , 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 26 m -2 for both epochs, and log ξ = 2.8 and N H = 10 27 m -2 was needed for 2016 only. The obscurers were located within the BLR. In Mao et al. (2022) (Paper III, hereafter), 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.

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;Struü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 -1 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 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. 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. 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.

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.
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).

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. We also obtained the 0.35 -10 keV time-averaged spectrum from each observation. Fig. 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).
We initially modelled the time-averaged spectrum from each observation, corresponding to the base model, before fitting each Article number, page 3 of 19 A&A proofs: manuscript no. NGC3227_Paper_Final_arXiv  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.
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.

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 Fe kα 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 Xray 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. 2011Mehdipour et al. , 2017Mehdipour et al. , 2018Petrucci 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 2 (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 compo-  Notes. * The covering fraction (CF) for the obscurer components is f cov and for the warm emitter is C cov = Ω/4π, where Ω is the solid angle subtended by the warm emitter with respect to the SMBH. nent, 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 HOT H = 2.07 × 10 24 m -2 (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).

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).
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 for-Article number, page 5 of 19 mally 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.

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.
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 et al. 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 26 m -2 to 9.93 × 10 26 m -2 , respectively. Conversely, the observed N pow and N comp time-averaged parameters decreased from 2.60 to 0.72×10 50 ph s -1 keV -1 and from 1.89 to 0.16 × 10 53 ph s -1 keV -1 , 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 26 m -2 ) 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. 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.

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 sig-   Table 3. The bottom panel shows the residuals between the data and model for Obs 1 (blue) and Obs 2 (purple).
nificantly 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.
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 /dof i , where dof = 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 /dof i ratios (where dof = 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 /dof i ratios in each TB when the parameters of interest are fixed, come from the same distribution as the best fit σ free = C i /dof 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 Article number, page 7 of 19 A&A proofs: manuscript no. NGC3227_Paper_Final_arXiv  . Right: Correlations between N pow , N comt and N H 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. 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 Figure 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).
There are slight differences between the widths of the 68% confidence contours (black lines) of N pow (0.48 to 0.55 ×10 50 ph s -1 keV -1 ) and N comt (0.40 to 0.48 ×10 53 ph s -1 keV -1 ) 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 26 m -2 ). For the time-averaged contours, the intervals covered by the 68% confidence contours are ∆N pow = 0.16 × 10 50 ph s -1 keV -1 , ∆N comt = 0.10 × 10 53 ph s -1 keV -1 and ∆N H = 0.20 × 10 26 m -2 . 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 Table 3: Best fit parameter results fitted to the 0.35 -10 keV EPIC-PN time-averaged ('All') and TB spectra. Continuum parameters: power law normalisation and photon index (N pow and Γ pow , respectively); Comptonisation normalisation (N comt ); reflection scaling (s). Obscurer parameters: column density (N H ); ionisation parameter (ξ); covering fraction ( f cov ). Obs 1 results are shown in the top half and the bottom half displays the best fit Obs 2 results. (1) 10 50 ph s -1 keV -1 at 1 keV; (2) 10 53 ph s -1 keV -1 at 1 keV; (3) 10 26 m -2 ; (4) 10 −9 Wm; (5) 10 49 ph s -1 keV -1 at 1 keV; (6) 10 52 ph s -1 keV -1 at 1 keV. 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.
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 −5 (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 .  Table 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 26 m -2 , 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 50 ph s -1 keV -1 , N comt = 1.89 × 10 53 ph s -1 keV -1 , Γ 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 .
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. Fig. 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 multiphased 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 -1see 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 v cross = GM BH r , where G is the gravitational constant, M BH = 5.96 × 10 6 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 -1 , depending on which distance is used. For this range in velocities, the obscurer would travel a distance of 6.80 × 10 10 − 1.51 × 10 11 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 11 m. Here, there is clearly an overlap Article number, page 11 of 19 A&A proofs: manuscript no. NGC3227_Paper_Final_arXiv Fig. 10: Monte Carlo spectral simulation results for the timeaveraged model (top panel; red) and the models from TB2 (topmiddle; blue), TB6 (bottom-middle; brown), and TB9 (bottom panel; cyan). The left and right sides of each panel display the results of N H against N pow and N H and N comt , 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 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 -1 , is 1.23 − 2.65 × 10 12 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 Figure 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 et al. 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.  (Table 3), and fitting N pow (red circles) and N comt (blue triangles) in each TB. Right: Fixing the continuum parameters to their time-averaged values and fitting N H (magenta diamonds) and f cov (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 N pow (left) and N H (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.

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 2 ) and the measured errors (σ 2 err ) give rise to the excess variance σ 2 XS . Here, σ 2 err,i , where x i is the measured value (flux), x is the mean, N is the number of measurements, and σ 2 err,i is the uncertainty, or mean square error, on each measurement. The fractional rms variability amplitude is related to the excess variance via To compute the uncertainty on F var , we used Eq. B2 in Appendix B 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 (σ 2 XS ) values resulting from the much weaker signal in Obs 2, compared to Obs 1.
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 -2 . The uncertainties on the measured values for the TBs (σ 2 err,i ) 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 Cstatistic 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 Article number, page 13 of 19 A&A proofs: manuscript no. NGC3227_Paper_Final_arXiv 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 et al. 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 −2 −5×10 −4 Hz for the averaged values, and 10 −4 , 5×10 −5 , and 2×10 −5 , for 10, 20 and 50 ks, respectively (Arévalo et al. 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 et al. 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).

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 26 m -2 to 9.93 × 10 26 m -2 , respectively. Comparatively, N pow and Fig. 12: Fractional root mean square (rms) variability amplitude (F var ) 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 (σ 2 XS ) 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 F var 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 F var 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. N comp decreased from 2.60 to 0.73 × 10 50 ph s -1 keV -1 at 1 keV, and 1.89 to 0.16 × 10 53 ph s -1 keV -1 , 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.  (N H,2 ). Right: Covering fractions of the two obscurer components, orange for obscurer 1 ( f cov,1 ) and magenta for obscurer 2 ( f cov,2 ). The dashed lines in each panel show the time-averaged best fit parameter values, with their uncertainties displayed as the shaded areas.
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. 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 spheri-cal, it will have a mass given by 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 M = 4πr 2 m p µN H,1 f cov,1 = 4πR 2 m p µN H,10 f cov,10 , (D.2) where r and R are the radii of the obscurer in TB1 and TB10, respectively. A quick rearrangement gives R 2 = Ar 2 where A = 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 = √ 2.25 r. Here, we set r = N H,1 n H for TB1, where n H is the hydrogen number density. From Paper III (Table 4), n H = 0.1 − 2.7 × 10 15 m -3 for a spherical cloud. Taking the average of n H gives an obscurer radius in TB1 of r = 2.24 × 10 11 m, and therefore yields a radius after expansion in TB10 of R = 3.36 × 10 11 m. So the expansion distance travelled during Obs 1 would therefore be D = R − r = 1.12 × 10 11 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 -1 . Therefore, the obscurer would have to expand with a velocity of around 1120 km s -1 , which is in the middle of the estimated crossing velocity (v cross ) of the obscurer (between 680 and 1470 km s -1 ), 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.
Article number, page 17 of 19 A&A proofs: manuscript no. NGC3227_Paper_Final_arXiv 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 T N ) 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.
Article number, page 18 of 19 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.