<figure><img src="image_1.png" /></figure>

HAL
open science

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

S Grafton-Waters, J Mao, M Mehdipour, G Branduardi-Raymont, M Page, J

Kaastra, Y Wang, C Pinto, G.A Kriss, D.J Walton, et al.

► To cite this version:

S Grafton-Waters, J Mao, M Mehdipour, G Branduardi-Raymont, M Page, et al. Transient obscu-ration event captured in NGC 3227 IV. Origin of the obscuring cloud variability. Astronomy and Astrophysics - A&A, 2023, 673, pp.A26. 10.1051/0004-6361/202243681. hal-04047654

HAL Id: hal-04047654

https://hal.science/hal-04047654

Submitted on 18 Aug 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Transient obscuration event captured in NGC 3227

## IV. Origin of the obscuring cloud variability

S. Grafton-Waters<fnref n="1" />, J. Mao<fnref n="2" />,<fnref n="3" />,<fnref n="4" />, M. Mehdipour<fnref n="5" />,<fnref n="4" />, G. Branduardi-Raymont<fnref n="1" />, M. Page<fnref n="1" />, J. Kaastra<fnref n="4" />,<fnref n="6" />, Y. Wang<fnref n="4" />,<fnref n="6" />,<fnref n="7" />,<fnref n="8" />, C. Pinto<fnref n="9" />, G. A. Kriss<fnref n="5" />, D. J. Walton<fnref n="10" />, P.-O. Petrucci<fnref n="11" />, G. Ponti<fnref n="12" />,<fnref n="13" />, B. De Marco<fnref n="14" />, S. Bianchi<fnref n="15" />, E. Behar<fnref n="16" />, and J. Ebrero<fnref n="17" />

Received 30 March 2022 / Accepted 6 March 2023

## ABSTRACT

*Context.* 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, which was observed with *XMM-Newton*, *NuSTAR*, and the *Hubble* Space Telescope.

*Aims.* We aim to determine the nature and origin of the observed spectral variability in the 2019 obscuration event.

*Methods.* We study the evolution of the obscurer by splitting the two *XMM-Newton* observations from 2019 into timing bins of length ~10 ks. We used the SPEX code to analyse the 0.35–10 keV EPIC-PN spectra of each timing bin.

*Results.* 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_{\text{pow}}$ and $N_{\text{comt}}$, respectively). The power-law 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 $N_H$ 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 the ionising continuum. The second observation (Obs 2) displays a significantly lower count rate due to the combination of a high $N_H$ and covering fraction of the obscurer, and a lower continuum flux.

*Conclusions.* 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.

**Key words.** X-rays: galaxies – techniques: spectroscopic – galaxies: active – galaxies: Seyfert – galaxies: individual: NGC 3227

## 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 under-

gone 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

<footnote n="1">¹ Mullard Space Science Laboratory, University College London, Holmbury St. Mary, Dorking, Surrey RH5 6NT, UK
e-mail: sam.waters.17@ucl.ac.uk</footnote>

<footnote n="2">² Department of Physics, Hiroshima University, 1-3-1 Kagamiyama, Higashi-Hiroshima, Hiroshima 739-8526, Japan</footnote>

<footnote n="3">³ Department of Physics, University of Strathclyde, Glasgow G4 0NG, UK</footnote>

<footnote n="4">⁴ SRON Netherlands Institute for Space Research, Niels Bohrweg 4, 2333 CA Leiden, The Netherlands</footnote>

<footnote n="5">⁵ Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA</footnote>

<footnote n="6">⁶ Leiden Observatory, Leiden University, PO Box 9513, 2300 RA Leiden, The Netherlands</footnote>

<footnote n="7">⁷ CAS Key Laboratory for Research in Galaxies and Cosmology, Department of Astronomy, University of Science and Technology of China, Hefei 230026, PR China</footnote>

<footnote n="8">⁸ School of Astronomy and Space Science, University of Science and Technology of China, Hefei 230026, PR China</footnote>

<footnote n="9">⁹ INAF-IASF Palermo, Via U. La Malfa 153, 90146 Palermo, Italy</footnote>

<footnote n="10">¹⁰ Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK</footnote>

<footnote n="11">¹¹ Univ. Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France</footnote>

<footnote n="12">¹² INAF – Osservatorio Astronomico di Brera, Via E. Bianchi 46, 23807 Merate (LC), Italy</footnote>

<footnote n="13">¹³ Max Planck Institute fur Extraterrestriche Physik, 85748 Garching, Germany</footnote>

<footnote n="14">¹⁴ Departament de Física, EEBE, Universitat Politècnica de Catalunya, Av. Eduard Maristany, 16, Barcelona 08019, Spain</footnote>

<footnote n="15">¹⁵ Dipartimento di Matematica e Fisica, Università degli Studi Roma Tre, via della Vasca Navale 84, 00146 Roma, Italy</footnote>

<footnote n="16">¹⁶ Department of Physics, Technion-Israel Institute of Technology, 32000 Haifa, Israel</footnote>

<footnote n="17">¹⁷ Telespazio UK for the European Space Agency (ESA), European Space Astronomy Centre (ESAC), Camino Bajo del Castillo, s/n, 28692 Villanueva de la Cañada, Madrid, Spain</footnote>
AGN. It is possible that many AGN have undergone obscu-
ration 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^{24} - 10^{26} m^{-2}$) and outflow veloc-
ities ($v_{out,WA} = \text{few } 100 - 1000 \text{ 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 NuS-
TAR (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 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 compo-
nent ($\log \xi \sim 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 \times 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 last-
ing a day (Turner et al. 2018). To explain the observed spec-
tral hardening and a depth change of the unresolved transition
array (UTA; at ~16–17 Å) in the RGS data, a component cov-
ering 60% of the X-ray continuum was required in the mod-
elling (with $N_H = 5 \times 10^{26}$ m$^{-2}$; $\log \xi = 2$; $v_{out} \sim -800$ km s$^{-1}$).

<table><caption>**Table 1.** Observation log for NGC 3227 in 2019.</caption>
   <thead>
    <tr>
     <td>
      Obs
     </td>
     <td>
      Obs ID
     </td>
     <td>
      Start date
     </td>
     <td>
      Duration (ks)
     </td>
     <td>
      Avg. CR (count s
      <sup>
       -1
      </sup>
      )
     </td>
    </tr>
   </thead>
   <tbody>
    <tr>
     <td>
      1
     </td>
     <td>
      0844341301
     </td>
     <td>
      2019/11/15
     </td>
     <td>
      103
     </td>
     <td>
      8.01
     </td>
    </tr>
    <tr>
     <td>
      2
     </td>
     <td>
      0844341401
     </td>
     <td>
      2019/12/05
     </td>
     <td>
      50
     </td>
     <td>
      1.74
     </td>
    </tr>
   </tbody>
  </table>

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<fnref n="1" />
ranging from $\log \xi = -1$ to 3, and outflow velocities in the range
$v_{\text{out}} = -100$ to $-1300 \text{ 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 \xi = 1 - 1.9$ and $N_H = 10^{26} \text{ m}^{-2}$ for both epochs,
and $\log \xi = 2.8$ and $N_H = 10^{27} \text{ m}^{-2}$ 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 mod-
elling the RGS, EPIC-PN, and NuSTAR data. Different models
were applied to each 2019 observation, with one or two pho-
tionisation 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 inves-
tigate 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

<footnote n="1">¹ The ionisation parameter is defined as $\xi = \frac{L_{\text{ion}}}{n_{\text{H}}r^2}$, where $L_{\text{ion}}$ is the ionising luminosity (measured between 13.6 eV and 13.6 keV), $n_{\text{H}}$ is the hydrogen number density, and $r$ is the distance of the plasma from the black hole.</footnote>
<figure><img src="image_3.png" /><figcaption>**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.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**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 ($HR_{avg}$) across each observation. For Obs 1 $HR_{avg} = -0.115$, and for Obs 2 $HR_{avg} = -0.033$.</figcaption></figure>

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-NoEP 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 EPIC-LCORR 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 ratio (blue line) 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
<figure><img src="image_5.png" /><figcaption>**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.</figcaption></figure>

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

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

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

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 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 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_{\text{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_{\text{dbb}}, N_{\text{pow}}, N_{\text{comt}}$, respectively), the reflection scaling parameter<fnref n="2" /> (s), and photon index ($\Gamma_{\text{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_{\text{refl}}$) and photon index ($\Gamma_{\text{refl}}$) values were coupled to the POW normalisation ($N_{\text{pow}}$) and photon index ($\Gamma_{\text{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_{\text{refl}}$ is fixed to the 2016 $N_{\text{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_{\text{pow}}, \Gamma_{\text{pow}}, N_{\text{comt}}, s$. In addition, we modelled the neutral gas in our Galaxy using the HOT

<footnote n="2">² The reflection scaling parameter is the fraction of the observed power law that is reflected, given by $s = F_{\text{refl}}/F_{\text{pow}}$, where $F_{\text{refl}}$ and $F_{\text{pow}}$ are the fluxes of the reflected and power law continua, respectively. It is identical to reflection scaling factor of the PEXRAV model in XSPEC, and represents a multiplicative scaling of the reflection once $N_{\text{refl}}$ is coupled to $N_{\text{pow}}$.</footnote>
<figure><img src="image_6.png" /><figcaption>**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.</figcaption></figure>

<table><caption>**Table 2.** Parameter values for the obscurer and warm emitter, taken from Tables 2 and 3 from Paper III.</caption><thead><tr><th rowspan="2">Parameter</th><th colspan="2">Obscurer</th><th rowspan="2">Warm emitter<br>(Obs 1)</th></tr><tr><th>(Obs 1)</th><th>(Obs 2)</th></tr></thead><tbody><tr><td><i>N</i><sub>H</sub> (10<sup>26</sup> m<sup>-2</sup>)</td><td>4.00</td><td>10.50</td><td>0.96</td></tr><tr><td>log <i>ξ</i> (10<sup>-9</sup> Wm)</td><td>0.50</td><td>2.60</td><td>1.37</td></tr><tr><td>CF<sup>(*)</sup></td><td>0.61</td><td>0.72</td><td>0.01</td></tr></tbody></table>

**Notes.** The warm emitter parameters are fixed throughout the modelling, while the obscurer parameters are fitted (see Sect. 4). (*)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.

model in SPEX; we fixed the column density and temperature at $N_H^{HOT} = 2.07 \times 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).

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

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
<figure><img src="image_7.png" /></figure>

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

<figure><img src="image_8.png" /></figure>

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

this as the baseline model to fit the spectra of each TB. The best fit parameter results are displayed in Table 3.

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

The top two panels in Fig. 7 display how the power law normalisation ($N_{\text{pow}}$; left) and photon index ($\Gamma_{\text{pow}}$; right) vary over time. For Obs 1, $N_{\text{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_{\text{pow}}$ stays constant between TBs, within the uncertainties. The photon index in Obs 2 displays strong changes, although similar to $N_{\text{pow}}$, they are consistent with each other within the uncertainties. In Obs 1, $\Gamma_{\text{pow}}$ does not vary throughout the observation, except for TB1 and TB2. For both observations, $\Gamma_{\text{pow}}$ is consistent with a

hard power law, similar to previous observations ($\Gamma = 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_{\text{comt}}$; left) and reflection scaling parameter ($s$; right) changes. For Obs 1, $N_{\text{comt}}$ has the same, but more significant, trend as $N_{\text{pow}}$. For Obs 2, similarly to $N_{\text{pow}}$, $N_{\text{comt}}$ is consistent between TBs. Although the reflection scaling parameter ($s$) of REFL in both observations appears to correlate well with $\Gamma_{\text{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_{\text{H}}$; left) and covering fraction ($f_{\text{cov}}$; right) of the obscurer component. In Obs 2, it appears that $N_{\text{H}}$ fluctuates significantly between TBs. However, with the exception of TB3, the $N_{\text{H}}$ values in each TB are constant within the uncertainties of each other and the time-averaged value. Therefore, the given changes in $N_{\text{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_{\text{H}}$ decreases through each TB, which follows the inverse trend compared to $N_{\text{pow}}$ and $N_{\text{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_{\text{cov}}$ stays constant within the uncertainties throughout each observation. However, $f_{\text{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_{\text{H}}$ and $f_{\text{cov}}$ compared to the other TBs. However, we are unable to constrain $\xi$ (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_{\text{cov}}$ and $N_{\text{H}}$ increased from 0.55 to 0.66, and from $3.87 \times 10^{26} \text{ m}^{-2}$ to $9.93 \times 10^{26} \text{ m}^{-2}$, respectively. Conversely, the observed $N_{\text{pow}}$ and $N_{\text{comp}}$ time-averaged
<table><caption>**Table 3.** Best fit parameter results fitted to the 0.35–10 keV EPIC-PN time-averaged ('All') and TB spectra.</caption><thead><tr><th colspan="9">Obs 1 (d.o.f = 116)</th></tr><tr><th>Time<br>Bin</th><th>N<sub>pow</sub><br>(1)</th><th>Γ</th><th>N<sub>comt</sub><br>(2)</th><th>s</th><th>N<sub>H</sub><br>(3)</th><th>log ξ<br>(4)</th><th>f<sub>cov</sub></th><th>C-stat</th></tr></thead><tbody><tr><td>All</td><td>2.60 ± 0.05</td><td>1.68 ± 0.01</td><td>1.89 ± 0.03</td><td>0.64 ± 0.04</td><td>3.87 ± 0.07</td><td>0.44<sup>+0.90</sup><sub>-2.00</sub></td><td>0.55 ± 0.01</td><td>221</td></tr><tr><td>1</td><td>2.46 ± 0.25</td><td>1.83 ± 0.05</td><td>1.01 ± 0.20</td><td>1.10 ± 0.30</td><td>6.27 ± 0.40</td><td>0.48<sup>+1.70</sup><sub>-0.90</sub></td><td>0.60 ± 0.04</td><td>158</td></tr><tr><td>2</td><td>2.40 ± 0.20</td><td>1.85 ± 0.05</td><td>0.87 ± 0.10</td><td>1.93 ± 0.30</td><td>5.09 ± 0.30</td><td>0.03<sup>+0.60</sup><sub>-3.40</sub></td><td>0.58 ± 0.03</td><td>163</td></tr><tr><td>3</td><td>2.09 ± 0.20</td><td>1.73 ± 0.04</td><td>1.15 ± 0.10</td><td>1.20 ± 0.20</td><td>5.98 ± 0.40</td><td>0.08<sup>+2.30</sup><sub>-1.20</sub></td><td>0.59 ± 0.03</td><td>127</td></tr><tr><td>4</td><td>2.21 ± 0.30</td><td>1.66 ± 0.04</td><td>1.52 ± 0.10</td><td>0.77 ± 0.15</td><td>5.09 ± 0.95</td><td>&lt;2.37</td><td>0.60 ± 0.04</td><td>133</td></tr><tr><td>5</td><td>2.17 ± 0.20</td><td>1.65 ± 0.04</td><td>1.52 ± 0.10</td><td>0.61 ± 0.15</td><td>4.93 ± 0.30</td><td>&lt;0.38</td><td>0.62 ± 0.03</td><td>157</td></tr><tr><td>6</td><td>2.77 ± 0.20</td><td>1.68 ± 0.04</td><td>1.76 ± 0.20</td><td>0.63 ± 0.15</td><td>3.96 ± 0.50</td><td>&lt;1.62</td><td>0.54 ± 0.04</td><td>146</td></tr><tr><td>7</td><td>3.24 ± 0.20</td><td>1.69 ± 0.03</td><td>2.34 ± 0.20</td><td>0.48 ± 0.10</td><td>3.45 ± 0.20</td><td>1.39 ± 0.30</td><td>0.56 ± 0.02</td><td>101</td></tr><tr><td>8</td><td>3.60 ± 0.20</td><td>1.73 ± 0.03</td><td>2.82 ± 0.20</td><td>0.58 ± 0.10</td><td>3.04 ± 0.20</td><td>1.02 ± 0.25</td><td>0.59 ± 0.02</td><td>114</td></tr><tr><td>9</td><td>3.03 ± 0.20</td><td>1.64 ± 0.03</td><td>3.01 ± 0.20</td><td>0.37 ± 0.10</td><td>2.70 ± 0.20</td><td>0.73 ± 0.35</td><td>0.57 ± 0.02</td><td>142</td></tr><tr><td>10</td><td>3.21 ± 0.20</td><td>1.67 ± 0.03</td><td>2.56 ± 0.20</td><td>0.31 ± 0.10</td><td>2.93 ± 0.20</td><td>1.33 ± 0.25</td><td>0.57 ± 0.02</td><td>136</td></tr></tbody></table><table><thead><tr><th colspan="9">Obs 2 (d.o.f = 106)</th></tr><tr><th>Time<br>Bin</th><th>N<sub>pow</sub><br>(5)</th><th>Γ</th><th>N<sub>comt</sub><br>(6)</th><th>s</th><th>N<sub>H</sub><br>(3)</th><th>log ξ<br>(4)</th><th>f<sub>cov</sub></th><th>C-stat</th></tr></thead><tbody><tr><td>All</td><td>7.23 ± 0.40</td><td>1.73 ± 0.03</td><td>1.56 ± 0.30</td><td>1.66 ± 0.15</td><td>9.93 ± 1.00</td><td>2.72 ± 0.10</td><td>0.66 ± 0.02</td><td>102</td></tr><tr><td>1</td><td>6.80 ± 1.10</td><td>1.68 ± 0.10</td><td>1.52 ± 0.70</td><td>1.58 ± 0.30</td><td>10.72 ± 1.60</td><td>2.91 ± 0.15</td><td>0.63 ± 0.04</td><td>93</td></tr><tr><td>2</td><td>8.12 ± 1.10</td><td>1.84 ± 0.08</td><td>1.30 ± 0.70</td><td>2.19 ± 0.40</td><td>7.68 ± 1.50</td><td>2.50 ± 0.25</td><td>0.67 ± 0.05</td><td>103</td></tr><tr><td>3</td><td>6.02 ± 0.90</td><td>1.62 ± 0.08</td><td>2.46 ± 0.70</td><td>1.32 ± 0.30</td><td>5.09 ± 1.00</td><td>&lt;2.24</td><td>0.57 ± 0.06</td><td>112</td></tr><tr><td>4</td><td>7.47 ± 1.00</td><td>1.72 ± 0.07</td><td>1.37 ± 0.70</td><td>1.42 ± 0.30</td><td>12.24 ± 2.00</td><td>2.80 ± 0.15</td><td>0.68 ± 0.05</td><td>129</td></tr><tr><td>5</td><td>8.43 ± 1.10</td><td>1.90 ± 0.09</td><td>1.56 ± 1.00</td><td>2.41 ± 0.40</td><td>9.02 ± 1.50</td><td>2.46 ± 0.25</td><td>0.74 ± 0.04</td><td>107</td></tr></tbody></table>

**Notes.** Continuum parameters: power law normalisation and photon index ($N_{\text{pow}}$ and $\Gamma_{\text{pow}}$, respectively); Comptonisation normalisation ($N_{\text{comt}}$); reflection scaling (s). Obscurer parameters: column density ($N_{\text{H}}$); ionisation parameter ($\xi$); covering fraction ($f_{\text{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.

parameters decreased from 2.60 to $0.72 \times 10^{50}$ ph s$^{-1}$ keV$^{-1}$ and from 1.89 to $0.16 \times 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 \times 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, 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 $\xi$ being free, while $f_{\text{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 ($\xi$) 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 $\xi$, 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 parame-

ter 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_{\text{pow}}$ and $N_{\text{comt}}$). On the other hand, these large $\xi$ errors could indicate an inhomogeneous obscurer that is made up of multiple components, each with a different $N_{\text{H}}$, $\xi$, 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_{\text{pow}}$ and $N_{\text{comt}}$, as well as with $N_{\text{H}}$. These are displayed in the left side of Fig. 8. For Obs 1, $N_{\text{pow}}$ and $N_{\text{comt}}$ relate well with Fig. 1, while an inverse trend was seen with $N_{\text{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.

The right side of Fig. 8 shows the relation between $N_{\text{pow}}$, $N_{\text{comt}}$, and $N_{\text{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_{\text{H}}$ and the count rate (bottom right panel in Fig. 8).
<figure><img src="image_9.png" /><figcaption>**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.</figcaption></figure>

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 param-
eter, 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 signif-
icance of thawing a parameter from its time-averaged value in the
<figure><img src="image_10.png" /><figcaption>**Fig. 8.** Correlations between model parameters. Left: The normalisations of the power law ($N_{\text{pow}}$) and Comptonisation component ($N_{\text{comt}}$) are correlated over the two observations (top); $N_{\text{pow}}$ against the column density of the obscurer ($N_{\text{H}}$) (middle); and $N_{\text{comt}}$ against $N_{\text{H}}$ (bottom). Right: Correlations between $N_{\text{pow}}$, $N_{\text{comt}}$ and $N_{\text{H}}$ against 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.</figcaption></figure>

spectral fitting. To do this, we fixed each of the four parameters of interest – $N_{\text{H}}$, $N_{\text{pow}}$, $N_{\text{comt}}$, and $f_{\text{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_{\text{fixed}} = \sum C_i/\text{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_{\text{free}} = C_i/\text{d.o.f.}_i$ ratios (where d.o.f. = 116 here) from the best fit TB results in Table 3, which is $S_{\text{free}} = 11.87$. We then calculate the difference between these two summations ($\Delta S = S_{\text{fixed}} - S_{\text{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 ($\sigma$) calculated from the $\sigma_{\text{fixed}} = C_i/\text{d.o.f.}_i$ ratios in each TB when

the parameters of interest are fixed, come from the same distribution as the best fit $\sigma_{\text{free}} = C_i/\text{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 = \sigma_{\text{fixed}}/\sigma_{\text{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 $\Delta S$ for each fixed parameter. If the F-value is larger than $\Delta 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_{\text{comt}}$ and $N_{\text{H}}$ are larger than the difference in the summed ratios ($\Delta S$), for the respective parameter. In addition, $N_{\text{comt}}$ and $N_{\text{H}}$ have sufficiently small $p$-values (less than 0.05) which implies that it is unlikely for the C-statistic variances ($\sigma_{\text{fixed}}$) of these two
<table><caption>**Table 4.** Testing the significance of the parameters on the best fit model (see Sect. 5.2 for details).</caption><thead><tr><th>Parameter</th><th>S<sub>fixed</sub></th><th>ΔS</th><th>F-value</th><th>p-value</th></tr></thead><tbody><tr><td>N<sub>pow</sub></td><td>12.78</td><td>0.91</td><td>0.69</td><td>0.98</td></tr><tr><td>N<sub>comt</sub></td><td>14.72</td><td>2.85</td><td>4.81</td><td>3.33 × 10<sup>-16</sup></td></tr><tr><td>N<sub>H</sub></td><td>13.62</td><td>1.75</td><td>2.10</td><td>3.75 × 10<sup>-5</sup></td></tr><tr><td>f<sub>cov</sub></td><td>12.09</td><td>0.22</td><td>0.82</td><td>0.85</td></tr></tbody></table>

parameters to come from the same distribution as the best fit C-statistic variance ($\sigma_{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 $\Delta 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 $\Delta 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).

There are slight differences between the widths of the 68% confidence contours (black lines) of $N_{pow}$ (0.48 to $0.55 \times 10^{50}$ ph s$^{-1}$ keV$^{-1}$) and $N_{comt}$ (0.40 to $0.48 \times 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 \times 10^{26}$ m$^{-2}$). For the time-averaged contours, the intervals covered by the 68% confidence contours are $\Delta N_{pow} = 0.16 \times 10^{50}$ ph s$^{-1}$ keV$^{-1}$, $\Delta N_{comt} = 0.10 \times 10^{53}$ ph s$^{-1}$ keV$^{-1}$ and $\Delta N_H = 0.20 \times 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 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 \times 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}$.

### 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 $\xi$ 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 \times 10^{26}$ m$^{-2}$, $\log \xi = 0.44$, and $f_{cov} = 0.55$ (see Table 3), and fitted the normalisations of POW and COMT, in addition to $\Gamma_{pow}$ and $\xi$. 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 \times 10^{50}$ ph s$^{-1}$ keV$^{-1}$, $N_{comt} = 1.89 \times 10^{53}$ ph s$^{-1}$ keV$^{-1}$, $\Gamma_{pow} = 1.68$ and $\xi = 0.64$, and fitted $N_H$, $\xi$ 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 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 $\xi$ 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 $\xi$ does change, but given the errors in Fig. 7, we cannot conclude one way
<figure><img src="image_11.png" /><figcaption>**Fig. 9.** Contour plots of $N_H$ (10$^{26}$ m$^{-1}$) against $N_{\text{pow}}$ (10$^{50}$ ph s$^{-1}$ keV$^{-1}$; top panels) and $N_H$ (10$^{26}$ m$^{-1}$) against $N_{\text{comt}}$ (10$^{53}$ 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 % ($\Delta C = 2.30$; black solid line), 90 % ($\Delta C = 4.61$; red dashed line), and 95 % ($\Delta 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.</figcaption></figure>

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$, $\xi$ 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_{\text{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 clumsy, 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$^{-1}$ – 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 $v_{\text{cross}} = \sqrt{\frac{GM_{\text{BH}}}{r}}$, where G is the gravitational constant, $M_{\text{BH}} = 5.96 \times 10^6 M_{\odot}$ 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 \times 10^{10} - 1.51 \times 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 \times 10^{11}$ 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_{\text{cross}} = 680 - 1470$ km s$^{-1}$, is $1.23 - 2.65 \times 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
<figure><img src="image_12.png" /><figcaption>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 N<sub>H</sub> against N<sub>pow</sub> and N<sub>comt</sub>, 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.</figcaption></figure>

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 lumi-
nosity 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 observa-
tion using the variance and fractional root mean square (rms)
variability amplitude (F<sub>var</sub>; 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<sub>var</sub> fol-
lowing the steps from Vaughan et al. (2003), where the differ-
ence between the source variance (S<sup>2</sup>) and the measured errors
(σ<sup>2</sup><sub>err</sub>) give rise to the excess variance σ<sup>2</sup><sub>XS</sub>. Here, S<sup>2</sup> = Σ<sub>i=1</sub><sup>N</sup> (x<sub>i</sub> - x̄)<sup>2</sup> / (N-1)

and σ̄<sup>2</sup><sub>err</sub> = 1/N Σ<sub>i=1</sub><sup>N</sup> σ<sup>2</sup><sub>err,i</sub>, where x<sub>i</sub> is the measured value (flux),

x̄ is the mean, N is the number of measurements, and σ<sup>2</sup><sub>err,i</sub> is the uncertainty, or mean square error, on each measurement. The fractional rms variability amplitude is related to the excess variance via

$$
F_{\text{var}} = \sqrt{\frac{\sigma_{\text{XS}}^2}{\bar{x}^2}} = \sqrt{\frac{S^2 - \sigma_{\text{err}}^2}{\bar{x}^2}}. \quad (1)
$$

To compute the uncertainty on F<sub>var</sub>, 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 (σ<sup>2</sup><sub>XS</sub>) 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<sub>var</sub> 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<sub>var</sub> 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.
<figure><img src="image_13.png" /><figcaption>**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 $N_{\text{pow}}$ (red circles) and $N_{\text{comt}}$ (blue triangles) in each TB. Right: Fixing the continuum parameters to their time-averaged values and fitting $N_{\text{H}}$ (magenta diamonds) and $f_{\text{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_{\text{pow}}$ (left) and $N_{\text{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.</figcaption></figure>

To determine what causes the features in the $F_{\text{var}}$ curve of Obs 1 (blue line in the top panel of Fig. 12), we obtained the $F_{\text{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_{\text{var}}$ spectra, we binned each TB spectrum (Fig. 4) into 20 energy bins and obtained the $F_{\text{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_{\text{var}}$ spectrum from the EPIC-PN light curves (blue line in the middle panel). The $F_{\text{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_{\text{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_{\text{H}} = 0 \text{ m}^{-2}$. The uncertainties on the measured values for the TBs ($\sigma_{\text{err,i}}^2$) are the standard error on the mean for each energy bin (similar to Lobban et al. 2020).

From plotting the $F_{\text{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_{\text{comt}}$ and $N_{\text{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_{\text{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_{\text{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_{\text{var}}$ spectra for seven different timing bin lengths, $\Delta t = 0.1, 0.5, 1, 2, 10, 20, \text{ and } 50 \text{ ks}$. For the four smallest timing bins, we found that the $F_{\text{var}}$ values are similar and therefore average these values together. We convert these timing bins to frequency bands of $10^{-2}–5 \times 10^{-4} \text{ Hz}$ for the averaged values, and $10^{-4}, 5 \times 10^{-5}, \text{ and } 2 \times 10^{-5}$, 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_{\text{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_{\text{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_{\text{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.
<figure><img src="image_14.png" /><figcaption>**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 ($\sigma_{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 $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 \times 10^{-4}$ Hz ($\Delta t = 100, 500, 1000$ and 2000 s; blue), $10^{-4}$ Hz ($\Delta t = 10$ ks; magenta), $5 \times 10^{-5}$ Hz ($\Delta t = 20$ ks; red), and $2 \times 10^{-5}$ Hz ($\Delta 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.</figcaption></figure>

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 \times 10^{26}$ m$^{-2}$ to $9.93 \times 10^{26}$ m$^{-2}$, respectively. Comparatively, $N_{pow}$ and $N_{comp}$ decreased from $2.60 \times 10^{50}$ ph s$^{-1}$ keV$^{-1}$ at 1 keV, and $1.89 \times 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.

**Acknowledgements.** 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 obscured 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

* Arévalo, P., & Markowitz, A. 2014, ApJ, 783, 82

* Arévalo, P., & Uttley, P. 2006, MNRAS, 367, 801

* Ballantyne, D. R., & Xiang, X. 2020, MNRAS, 496, 4255

* Beuchert, T., Markowitz, A. G., Krauß, F., et al. 2015, A&A, 584, A82

* Beuchert, T., Markowitz, A. G., Dauser, T., et al. 2017, A&A, 603, A50

* Cappi, M., Panessa, F., Bassani, L., et al. 2006, A&A, 446, 459

* Cappi, M., De Marco, B., Ponti, G., et al. 2016, A&A, 592, A27

* Cash, W. 1979, ApJ, 228, 939

* Chainakun, P., Watcharangkool, A., Young, A. J., & Hancock, S. 2019, MNRAS, 487, 667

* Connolly, S. D., McHardy, I. M., & Dwelly, T. 2014, MNRAS, 440, 3503

* Connolly, S. D., McHardy, I. M., & Dwelly, T. 2015, ArXiv e-prints

[arXiv:1502.07501]

* Couto, J. D., Kraemer, S. B., Turner, T. J., & Crenshaw, D. M. 2016, ApJ, 833, 191

* Crenshaw, D. M., Kraemer, S. B., Bruhweiler, F. C., & Ruiz, J. R. 2001, ApJ, 555, 633

* den Herder, J. W., Brinkman, A. C., Kahn, S. M., et al. 2001, A&A, 365, L7

* Detmers, R. G., Kaastra, J. S., Costantini, E., et al. 2010, A&A, 516, A61

* de Vaucouleurs, G., de Vaucouleurs, A., Corwin, Herold G., J., et al. 1991, Third Reference Catalogue of Bright Galaxies (Springer, New York)

* Ebroiro, J., Kriss, G. A., Kaastra, J. S., & Ely, J. C. 2016, A&A, 586, A72

* Gandhi, P., Hönig, S. F., & Kishimoto, M. 2015, ApJ, 812, 113

* Gehrels, N., Chincarini, G., Giommi, P., et al. 2004, ApJ, 611, 1005

* George, I. M., Mushotzky, R., Turner, T. J., et al. 1998, ApJ, 509, 146
Gondoin, P., Orr, A., Lumb, D., & Siddiqui, H. 2003, A&A, 397, 883

Grafton-Waters, S., Branduardi-Raymont, G., Mehdipour, M., et al. 2020, A&A, 633, A62

Grafton-Waters, S., Branduardi-Raymont, G., Mehdipour, M., et al. 2021, A&A, 649, A162

Green, J. C., Froning, C. S., Osterman, S., et al. 2012, ApJ, 744, 60

Harrison, F. A., Craig, W. W., Christensen, F. E., et al. 2013, ApJ, 770, 103

Ingram, A., & van der Klis, M. 2013, MNRAS, 434, 1476

Jansen, F., Lumb, D., Altieri, B., et al. 2001, A&A, 365, L1

Kaastra, J. S. 2017, A&A, 605, A51

Kaastra, J. S., & Bleeker, J. A. M. 2016, A&A, 587, A151

Kaastra, J. S., Mewe, R., & Nieuwenhuijzen, H. 1996, in UV and X-ray Spectroscopy of Astrophysical and Laboratory Plasmas, 411

Kaastra, J. S., Steenbrugge, K. C., Raassen, A. J. J., et al. 2002, A&A, 386, 427

Kaastra, J. S., Kriss, G. A., Cappi, M., et al. 2014, Science, 345, 64

Kaastra, J. S., Mehdipour, M., Behar, E., et al. 2018a, A&A, 619, A112

Kaastra, J. S., Raassen, A. J. J., de Plaa, J., & Gu, L. 2018b, https://doi.org/10.5281/zenodo.2419563

Keck, M. L., Brenneman, L. W., Ballantyne, D. R., et al. 2015, ApJ, 806, 149

Kinkhabwala, A., Sako, M., Behar, E., et al. 2002, ApJ, 575, 732

Kriss, G. A., Mehdipour, M., Kaastra, J. S., et al. 2019, A&A, 621, A12

Lamer, G., Uttley, P., & McHardy, I. M. 2003, MNRAS, 342, L41

Lobban, A. P., Turner, T. J., Reeves, J. N., Braito, V., & Miller, L. 2020, MNRAS, 494, 5056

Longinotti, A. L., Krongold, Y., Kriss, G. A., et al. 2013, ApJ, 766, 104

Longinotti, A. L., Kriss, G., Krongold, Y., et al. 2019, ApJ, 875, 150

Lyubarskii, Y. E. 1997, MNRAS, 292, 679

Magdziarz, P., & Zdziarski, A. A. 1995, MNRAS, 273, 837

Mao, J., Kaastra, J. S., Mehdipour, M., et al. 2022, A&A, 665, A72

Markowitz, A., Reeves, J. N., George, I. M., et al. 2009, ApJ, 691, 922

Markowitz, A. G., Krumpe, M., & Nikutta, R. 2014, MNRAS, 439, 1403

Mehdipour, M., Branduardi-Raymont, G., Kaastra, J. S., et al. 2011, A&A, 534, A39

Mehdipour, M., Kaastra, J. S., & Kallman, T. 2016a, A&A, 596, A65

Mehdipour, M., Kaastra, J. S., Kriss, G. A., et al. 2016b, A&A, 588, A139

Mehdipour, M., Kaastra, J. S., Kriss, G. A., et al. 2017, A&A, 607, A28

Mehdipour, M., Kaastra, J. S., Costantini, E., et al. 2018, A&A, 615, A72

Mehdipour, M., Kriss, G. A., Kaastra, J. S., et al. 2021, A&A, 652, A150

Mehdipour, M., Kriss, G. A., Costantini, E., et al. 2022, ApJ, 934, L24

Murphy, E. M., Lockman, F. J., Laor, A., & Elvis, M. 1996, ApJS, 105, 369

Nardini, E., Reeves, J. N., Gofford, J., et al. 2015, Science, 347, 860

Parker, M. L., Longinotti, A. L., Schartel, N., et al. 2019, MNRAS, 490, 683

Petrucci, P. O., Paltani, S., Malzac, J., et al. 2013, A&A, 549, A73

Ponti, G., Cappi, M., Dadina, M., & Malaguti, G. 2004, A&A, 417, 451

Porquet, D., Kaastra, J. S., Page, K. L., et al. 2004, A&A, 413, 913

Risaliti, G., Bianchi, S., Matt, G., et al. 2005, ApJ, 630, L129

Rivers, E., Markowitz, A., & Rothschild, R. 2011, ApJS, 193, 3

Strüder, L., Briel, U., Dennerl, K., et al. 2001, A&A, 365, L18

Titarchuk, L. 1994, ApJ, 434, 570

Tombesi, F., Cappi, M., Reeves, J. N., et al. 2013, MNRAS, 430, 1102

Turner, T. J., Reeves, J. N., Braito, V., et al. 2018, MNRAS, 481, 2470

Uttley, P., & McHardy, I. M. 2005, MNRAS, 363, 586

Vaughan, S., Edelson, R., Warwick, R. S., & Uttley, P. 2003, MNRAS, 345, 1271

Walton, D. J., Baloković, M., Fabian, A. C., et al. 2021, MNRAS, 506, 1557

Wang, Y., Kaastra, J., Mehdipour, M., et al. 2022, A&A, 657, A77
## 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_{\text{pow}}$ and $N_{\text{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_{\text{cov},1}$ decreases over the duration of the observation, while $f_{\text{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 not show this plot.

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_{\text{pow}}$ and $N_{\text{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_{\text{H}}$ and $f_{\text{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.

<figure><img src="image_15.png" /><figcaption>**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 ($N_{H,1}$) and red for component 2 ($N_{H,2}$). Right: Covering fractions of the two obscurer components, orange for obscurer 1 ($f_{\text{cov},1}$) and magenta for obscurer 2 ($f_{\text{cov},2}$). The dashed lines in each panel show the time-averaged best fit parameter values, with their uncertainties displayed as the shaded areas.</figcaption></figure>
## 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_{PN} = D_{PN}/M_{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_{RGS} = D_{RGS}/M_{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 ($\lambda_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<fnref n="3" />). 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 ($\lambda_1$) was fixed to 0.1 Å and $T_1$ was coupled to $T_2$. For the opposite end, $\lambda_9$ was fixed at 100 Å and $T_9$ was coupled to $T_8$. 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 $\lambda_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.

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.

<figure><img src="image_16.png" /></figure>

**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 ($\lambda_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.

<footnote n="3">³ https://www.sron.nl/astrophysics-spex/manual</footnote>
**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).

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}$, $\Gamma$, 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 $\xi$ against the 0.3 - 10 keV count rates, and $N_{comt}$ against $\xi$. 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 $\xi$ 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.
<figure><img src="image_17.png" /><figcaption>**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.</figcaption></figure>
## 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

$$M = 4\pi r^2 m_p \mu N_H f_{cov}, \qquad (D.1)$$

where $m_p$ is the proton mass, $\mu$ 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\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}, \qquad (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 = \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_H}$ for TB1, where $n_H$ is the hydrogen number density. From Paper III (Table 4), $n_H = 0.1 - 2.7 \times 10^{15} \text{ m}^{-3}$ for a spherical cloud. Taking the average of $n_H$ gives an obscurer radius in TB1 of $r = 2.24 \times 10^{11} \text{ m}$, and therefore yields a radius after expansion in TB10 of $R = 3.36 \times 10^{11} \text{ m}$. So the expansion distance travelled during Obs 1 would therefore be $D = R - r = 1.12 \times 10^{11} \text{ m}$.

Finally, the expansion velocity can be estimated if we take the time to be $t = 103 \text{ ks}$, and as such $v_{exp} = D/t \sim 1120 \text{ km s}^{-1}$. Therefore, the obscurer would have to expand with a velocity of around $1120 \text{ 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.