# Developing Optimization-Based Inverter Models for Short Circuit Studies

Thabiso R. Mabote  School of Electrical Engineering
& Computer Science
Oregon State University
Corvallis, OR, USA
<EMAIL_ADDRESS>Jose E. Tabarez  Advanced Network Science Initiative
Los Alamos National Laboratory
Los Alamos, NM, USA
<EMAIL_ADDRESS>Arthur K. Barnes  Advanced Network Science Initiative
Los Alamos National Laboratory
Los Alamos, NM, USA
<EMAIL_ADDRESS>Adam Mate  Advanced Network Science Initiative
Los Alamos National Laboratory
Los Alamos, NM, USA
<EMAIL_ADDRESS>Russell W. Bent  Advanced Network Science Initiative
Los Alamos National Laboratory
Los Alamos, NM, USA
<EMAIL_ADDRESS>Eduardo Cotilla-Sanchez  School of Electrical Engineering
& Computer Science
Oregon State University
Corvallis, OR, USA
<EMAIL_ADDRESS>
###### Abstract

As inverter-based generation becomes more common in distribution networks, it
is important to create models for use in optimization-based problems that
accurately represent their non-linear behavior when saturated. This work
presents models for grid-following and grid-forming inverters, and
demonstrates their use in optimization-based fault studies. The developed
models are shown to provide results in line with experimental tests from
literature; however, grid-forming inverter models fail to provide feasible
solutions for certain fault types even when currents, voltages, and powers
from the models seem reasonable. This work lays down the foundation for the
development of relaxed models for both grid-following and grid-forming
inverter models for use in optimization problems.

###### Index Terms:

power system operation, microgrid, distribution network, protection,
optimization, protective relaying.

††publicationid: pubid: 979-8-3503-4743-2/23/$31.00 ©2023 IEEE

##  I Introduction

With the increasing adoption of inverter-interfaced generation in microgrids
to ensure resilient operation of critical energy infrastructures and meet
decarbonization goals, comes the challenge of ensuring protective
coordination. As the size and mixture of distributed energy resources
increases in microgrids, so comes the requirement to ensure protection
coordination during different modes of operation to minimize the size of
affected areas during faults to avoid interruption of critical loads. To
achieve this goal, it is important to accurately model the short-circuit
current contributions from inverter-interfaced generation for various types of
faults and modes of operation [1, 2].

Experimental tests have demonstrated that the control schema and pre-fault
state of inverters have a profound effect on the fault currents injected by
the inverters [3, 4, 5, 6, 7, 8]. Developing mathematical models to represent
the inverters in optimization problems, based on the experimental tests,
results in non-linear non-convex constraints leading to feasibility issues [9,
10]. These constraints limit the size and complexity of optimization problems
that the inverter models can be used in, especially in optimization-based
fault studies used to develop protection coordination. To overcome this issue,
this work investigates the use of Thevenin circuits to approximate continuous
mathematical models of grid-forming and grid-following inverters to determine
their short-circuit current contributions, and compares the performance of
such models to the experimental results found in literature.

The inverter models presented in this work are intended to be used in
PowerModelsProtection.jl (PMsP) [9], a structure established on optimization
for short-circuit analysis. PMsP is an extension to a collection of open-
source packages using the PowerModels.jl (PMs) [11] framework; PMs allows for
different optimization problems in power systems to be developed and solved.
Presently, PMsP only produces mixed linear/nonlinear non-convex short-circuit
formulations, and it is also dependent on PowerModelsDistribution.jl (PMsD)
[12] for creating the optimization-based formulation and mathematical models
used for distribution networks.

This work focuses on developing inverter models for use in optimization based
short-circuit studies, however, the models are expressed in a general form for
use in power flow (pf) and optimal power flow (opf) analysis of balanced and
unbalanced distribution networks. The goal here is to develop continuous
variable inverter models that lay the foundation for developing relaxed linear
convex approximations in future work.

The contributions include the following: Section II presents the optimization-
based short-circuit formulation implemented in PMsP. Section III presents the
developed grid-following and grid-forming inverter models. Section IV
introduces the case study system used to test and verify the operation of the
inverter models, and Section V discusses and interprets the observed results.
Last, Section VI summarizes the conclusions regarding the state and importance
of the work.

##  II Short-Circuit Formulation

The symbols below are used in the following constraint formulation:
$\mathcal{B}$ is a collection of buses in a system; $\mathcal{G}_{ref_{b}}$ is
a collection of reference voltage sources at bus $b$;
$\Phi_{b}\subseteq\\{A,B,C\\}$ is the collection of phases at bus $b$;
$\mathcal{G}_{b}$ is a collection of generators at bus $b$; $\mathcal{F}$ is
the fault in a system; $\mathcal{E}_{b}$ is a collection of lines linked to
bus $b$; $\mathcal{T}_{b}$ is a collection of transformers linked to bus $b$;
and $\mathcal{S}_{b}$ is a collection of shunts at bus $b$.

The optimization-based short-circuit current formulation is implemented using
rectangular coordinates inside of PMsP. To indicate the real and imaginary
parts the subscripts $r$ and $i$ are applied to the variables that are
complex. The short-circuit formulation presented here is based on [9] and
[10]. The formulation defines a set of constraints that are based on Ohm’s law
and Kirchhoff’s current law, which define the phase voltages at the buses due
to current flow through various components within the system. These
constraints resemble constraints used in power flow, except the load
constraints are neglected and constraints that define the faulted condition in
the system are added.

Reference Voltage Sources Constraints:
$\forall g\in\mathcal{G}_{ref_{n}},~{}\forall b\in{\cal
B},~{}\forall\phi\in\Phi_{b}$

$\displaystyle V^{\phi}_{g_{r}}=V^{\phi}_{g_{setp.}}\cdot
cos\left(\theta_{g_{setp.}}\right)$ (1a) $\displaystyle
V^{\phi}_{g_{i}}=V^{\phi}_{g_{setp.}}\cdot sin\left(\theta_{g_{setp.}}\right)$
(1b) $\displaystyle V^{\phi}_{b_{r}}=V^{\phi}_{g_{r}}-r^{\phi}_{g}\cdot
I^{\phi}_{g_{r}}+x^{\phi}_{g}\cdot I^{\phi}_{g_{i}}$ (1c) $\displaystyle
V^{\phi}_{b_{i}}=V^{\phi}_{g_{i}}-r^{\phi}_{g}\cdot
I^{\phi}_{g_{i}}-x^{\phi}_{g}\cdot I^{\phi}_{g_{r}}$ (1d)

Generator Power Constraints:
$\forall g\in{\cal G},~{}\forall b\in{\cal B},~{}\forall\phi\in\Phi_{b}$

$\displaystyle V^{\phi}_{g_{r}}\cdot I^{\phi}_{g_{r}}+V^{\phi}_{g_{i}}\cdot
I^{\phi}_{g_{i}}$ $\displaystyle=P^{\phi}_{g}$ (2a) $\displaystyle
V^{\phi}_{g_{i}}\cdot I^{\phi}_{g_{r}}-V^{\phi}_{g_{r}}\cdot I^{\phi}_{g_{i}}$
$\displaystyle=Q^{\phi}_{g}$ (2b) $\displaystyle\Delta P^{\phi}_{g}$
$\displaystyle\leq 5\%$ (2c) $\displaystyle\Delta Q^{\phi}_{g}$
$\displaystyle\leq 5\%$ (2d)
$\displaystyle-\frac{P^{\phi}_{g}}{pf}\cdot\sin(\textrm{acos}(pf))\leq
Q^{\phi}_{g}$
$\displaystyle\leq\frac{P^{\phi}_{g}}{pf}\cdot\sin(\textrm{acos}(pf))$ (2e)

Fault Current Constraints:
$\forall f_{b}\in\mathcal{F},~{}\forall b\in{\cal
B},~{}\forall\phi\in\Phi_{b}$

$\displaystyle I^{\phi}_{{f_{b}}_{r}}=\sum_{c\in\Phi}\left(G_{\phi,c}\cdot
V^{c}_{b_{r}}\right)$ (3a) $\displaystyle
I^{\phi}_{{f_{b}}_{i}}=\sum_{c\in\Phi}\left(G_{\phi,c}\cdot
V^{c}_{b_{i}}\right)$ (3b)

Kirchhoff’s Current Constraints to Unfaulted Buses:
$\forall b\in{\cal B},~{}\forall\phi\in\Phi_{b}$

$\displaystyle\sum_{(b,j)\in\mathcal{E}_{b}}I^{\phi}_{{(b,j)}_{r}}+\sum_{(b,j)\in\mathcal{T}_{b}}I^{\phi}_{{(b,j)}_{r}}=\sum_{g\in\mathcal{G}_{b}}I^{\phi}_{{g}_{r}}+\sum_{s\in\mathcal{S}_{b}}I^{\phi}_{{s}_{r}}$
(4a)
$\displaystyle\sum_{(b,j)\in\mathcal{E}_{b}}I^{\phi}_{{(b,j)}_{i}}+\sum_{(b,j)\in\mathcal{T}_{b}}I^{\phi}_{{(b,j)}_{i}}=\sum_{g\in\mathcal{G}_{b}}I^{\phi}_{{g}_{i}}+\sum_{s\in\mathcal{S}_{b}}I^{\phi}_{{s}_{i}}$
(4b)

Kirchhoff’s Current Constraints to Faulted Buses:
$\forall b\in{\cal B},~{}\forall\phi\in\Phi_{b}$

$\displaystyle\sum_{b,j\in\mathcal{E}_{b}}I^{\phi}_{{(b,j)}_{r}}+\sum_{b,j\in\mathcal{T}_{b}}I^{\phi}_{{(b,j)}_{r}}=\sum_{g\in\mathcal{G}_{b}}I^{\phi}_{{g}_{r}}+\sum_{s\in\mathcal{S}_{b}}I^{\phi}_{{s}_{r}}+I^{\phi}_{{f}_{r}}$
(5a)
$\displaystyle\sum_{b,j\in\mathcal{E}_{b}}I^{\phi}_{{(b,j)}_{i}}+\sum_{b,j\in\mathcal{T}_{b}}I^{\phi}_{{(b,j)}_{i}}=\sum_{g\in\mathcal{G}_{b}}I^{\phi}_{{g}_{i}}+\sum_{s\in\mathcal{S}_{b}}I^{\phi}_{{s}_{i}}+I^{\phi}_{{f}_{i}}$
(5b)

Voltage Drop Constraints:
$\forall(b,j)\in{\cal E},~{}\forall b\in{\cal B},~{}\forall\phi\in\Phi_{b}$

$\displaystyle V^{\phi}_{{b}_{r}}=V^{\phi}_{{b}_{r}}-r^{\phi}_{(b,j)}\cdot
I^{\phi}_{{(b,j)}_{r}}+x^{\phi}_{(b,j)}\cdot I^{\phi}_{{(b,j)}_{i}}$ (6a)
$\displaystyle V^{\phi}_{{b}_{i}}=V^{\phi}_{{b}_{i}}-r^{\phi}_{(b,j)}\cdot
I^{\phi}_{{(b,j)}_{i}}-x^{\phi}_{(b,j)}\cdot I^{\phi}_{{(b,j)}_{r}}$ (6b)
$\displaystyle
V^{\phi}_{b_{min}}\leq(V^{\phi}_{b_{r}})^{2}+(V^{\phi}_{b_{i}})^{2}\leq
V^{\phi}_{b_{max}}$ (6c)
$\displaystyle(I^{\phi}_{{(b,j)}_{r}})^{2}+(I^{\phi}_{{(b,j)}_{i}})^{2}\leq
I^{\phi}_{(b,j)_{thermal}}$ (6d)

Transformer Constraints:
$\forall(b,j)\in{\cal T},~{}\forall b\in{\cal B},~{}\forall\phi\in\Phi_{b}$

$\displaystyle W_{(b,j)}\cdot V^{\phi}_{{b}_{r}}=\eta_{(b,j)}\cdot
V^{\phi}_{{j}_{r}}$ (7a) $\displaystyle W_{(b,j)}\cdot
V^{\phi}_{{b}_{i}}=\eta_{(b,j)}\cdot V^{\phi}_{{j}_{i}}$ (7b) $\displaystyle
I^{\phi}_{{b}_{i}}=\eta_{(b,j)}I^{\phi}_{{j}_{i}}$ (7c)
$\displaystyle(I^{\phi}_{{(b,j)}_{r}})^{2}+(I^{\phi}_{{(b,j)}_{i}})^{2}\leq
I^{\phi}_{(b,j)_{thermal}}$ (7d)

The fault currents, as shown in Constraints (3a)-(3b), are derived using a
fault admittance matrix, which denotes a fault impedance between phase and
ground. Line-to-line and 3-phase to ground faults cannot be modeled directly
due to the presence of the impedances between the phases and ground,
consequently an equivalent impedance admittance matrix is derived using a
star-mesh transformation [9]. Along with fault admittance constraints and
Constraints (4) and (5), which couples the nodes of the network using
Kirchhoff’s current law, allows for the formulation not to require an explicit
objective function. Constraints (1a)-(1d) represents the reference voltage
sources; normally the distribution substation is the main voltage source in
grid-connected mode of operation, whereas it is modeled as a voltage source
behind an impedance during a fault, and this impedance is obtained from the
single-phase and 3-phase short-circuit powers (Constraint (1b)). Synchronous
generation within the network is defined by Constraints (2a)-(2e). The last
set of constraints, Constraints (7a)-(7d), define the voltage and current
changes across lines and transformers within the network.

##  III Inverter Models

The models presented in [9] define inverters with complex mathematical models
that result in infeasible solutions under certain fault conditions. This can
be observed in Fig. 1 and Fig. 2: for low fault impedances located away from
the inverter (Fig. 1), the optimization solver is able to find a feasible
solution for both grid-forming and grid-following inverters; however, as the
fault was moved closer to the inverters (Fig. 2), the grid-following inverter
model resulted in infeasible solutions causing large deviations in current.

Figure 1: Fault Current (A) against Fault Resistance (R) for pv_fl_current
(grid-following) model and pv_fm_current (grid-forming) model. Figure 2: Fault
Current (A) against Fault Resistance (R) for pv_fl_current (grid-following)
model and pv_fm_current (grid-forming) model.

This work aims to develop simplified inverter models based on Thevenin
approximations for more general use in fault studies. Thevenin approximations
are widely used to model generation sources in fault studies; the difficulty
is to incorporate the inverter’s non-linear current response within them. When
inverters are not saturated, they respond similarly to the Thevenin models.
Incorporating the current saturation into the mathematical constraints of the
optimization problem is difficult when trying to use continuous decision
variables.

The developed inverter models in this work are compared to experimental
results to evaluate their performance on different fault types. It must be
noted that the models presented below do not directly model the grid
interfacing transformers; the constraints for the inverters are defined such
that they help compensate for the transformer models not being directly
modeled.

###  III-A Grid-Following Inverter Model

The grid-following inverter model was designed to act as a constant-power
source up to the point where the inverter’s current saturates, then become a
constant-current source on the saturated phase.

$\forall g\in{\cal G}_{gfli},~{}\forall b\in{\cal
B},~{}\forall\phi\in\Phi_{b}$

$\displaystyle{I^{\phi}_{g_{r}}}^{2}+{I^{\phi}_{g_{i}}}^{2}\leq{I^{\phi}_{g_{max}}}^{2}$
(8a) $\displaystyle V^{\phi}_{g_{r}}\cdot
I^{\phi}_{g_{r}}+V^{\phi}_{g_{i}}\cdot I^{\phi}_{g_{i}}=P^{\phi}_{g}-z^{\phi}$
(8b) $\displaystyle
M\cdot(1-\frac{1}{{I^{\phi}_{g_{max}}}^{2}}({I^{\phi}_{g_{r}}}^{2}+{I^{\phi}_{g_{i}}}^{2}))\cdot
z^{\phi}\geq 0$ (8c) $\displaystyle V^{\phi}_{g_{r}}\cdot
I^{\phi}_{g_{r}}+V^{\phi}_{g_{i}}\cdot I^{\phi}_{g_{i}}\geq 0$ (8d)
$\displaystyle V^{\phi}_{g_{i}}\cdot I^{\phi}_{g_{r}}-V^{\phi}_{g_{r}}\cdot
I^{\phi}_{g_{i}}\leq Q^{\phi}_{g_{max}}$ (8e) $\displaystyle
V^{\phi}_{g_{i}}\cdot I^{\phi}_{g_{r}}-V^{\phi}_{g_{r}}\cdot
I^{\phi}_{g_{i}}\geq Q^{\phi}_{g_{min}}$ (8f) $\displaystyle 0\leq
z^{\phi}\leq 1$ (8g)

Constraint (8a) ensures that the inverter-sourced current is below the maximum
possible sourced current $I^{\phi}_{g_{max}}$ that the inverter can supply,
which is a parameter defined by the User or set by PMsP. Constraint (8b)
defines the active power sourced from the inverter as a function of bus
voltage and inverter current; $P^{\phi}_{g}$ is the set-point for the active
power per phase of the inverter. Under non-faulted or light faulted cases the
inverter will supply the defined power, but as a fault drives down the
voltage, the current will increase to meet the power requirement; once the
current reaches its limit, Constraint (8c) activates the continuous variable
Z, which in Constraint (8b) allows the inverter to drop its active power
output based on the voltage collapse due to a fault. The last set of
constraints, Constraints (8d)-(8f), ensure that active power is supplied by
the inverter and never absorbed, and that the inverter operates within a fixed
power factor range by setting the limits on the reactive power supplied or
absorbed.

###  III-B Grid-Forming Inverter Model

Two Thevenin approximation models were developed and tested to simulate a
grid-forming inverter. The first is a simplified model of a voltage source
behind an adjustable impedance. The second is a voltage source that can adjust
its voltage magnitude on each phase.

Simplified Grid-Forming Inverter

$\forall g\in{\cal G}_{gfmi},~{}\forall b\in{\cal
B},~{}\forall\phi\in\Phi_{b}$

$\displaystyle{I^{\phi}_{g_{r}}}^{2}+{I^{\phi}_{g_{i}}}^{2}\leq{I^{\phi}_{g_{max}}}^{2}$
(9a) $\displaystyle
M\cdot(1-\frac{1}{{I^{\phi}_{g_{max}}}^{2}}({I^{\phi}_{g_{r}}}^{2}+{I^{\phi}_{g_{i}}}^{2}))\cdot
r^{\phi}\geq 0$ (9b) $\displaystyle
V^{\phi}_{g_{r}}=V^{\phi}_{g_{r0}}-r^{\phi}\cdot I^{\phi}_{g_{r}}$ (9c)
$\displaystyle V^{\phi}_{g_{i}}=V^{\phi}_{g_{i0}}-r^{\phi}\cdot
I^{\phi}_{g_{i}}$ (9d) $\displaystyle 0\leq
r^{\phi}\leq\frac{1}{I^{\phi}_{g_{max}}}$ (9e)

Constraint (9a) ensures that the inverter sourced current is below the maximum
possible sourced current $I^{\phi}_{g_{max}}$. Constraint (9b) makes the
resistor an activation variable that is zero when the current is below
$I^{\phi}_{g_{max}}$ and is non-zero when the inverter is saturated in order
to allow the inverter’s bus voltage to collapse. Constraints (9c) and (9d)
define the Thevenin approximation, using a resistor as the Thevenin impedance;
the value of the resistor is determined by the solver as a value between 0 and
a maximum value based on the maximum current supplied by the inverter
(Constraint (9e)).

Complex Grid-Forming Inverter

Constraint (10a) ensures that the inverter sourced current is below the
maximum possible sourced current $I^{\phi}_{g_{max}}$. Constraints (10c) and
(10d) define the inverter voltage applied to a bus: $V^{\phi}_{o_{g}}$ defines
the voltage magnitude of the inverter pre-fault, which has to be equal to the
bus voltage when the current is not saturated; when the current saturates, the
activation variable $z^{\phi}$ allows the inverter’s voltage magnitude to
adjust (Constraint (10b). Constraints (10e)-(10g) enforce balanced phase
angles across the inverter’s voltages. The last set of constraints,
Constraints (10h)-(10j), ensure that the inverter is operating within its
power limits.

$\forall g\in{\cal G}_{gfmi},~{}\forall b\in{\cal
B},~{}\forall\phi\in\Phi_{b}$

$\displaystyle{I^{\phi}_{g_{r}}}^{2}+{I^{\phi}_{g_{i}}}^{2}\leq{I^{\phi}_{g_{max}}}^{2}$
(10a) $\displaystyle
M\cdot(1-\frac{1}{{I^{\phi}_{g_{max}}}^{2}}({I^{\phi}_{g_{r}}}^{2}+{I^{\phi}_{g_{i}}}^{2}))\cdot
z^{\phi}\geq 0$ (10b)
$\displaystyle{V^{\phi}_{g_{r}}}^{2}+{V^{\phi}_{g_{i}}}^{2}\leq{V^{\phi}_{g_{0}}}^{2}\cdot(1+z^{\phi})$
(10c)
$\displaystyle{V^{\phi}_{g_{r}}}^{2}+{V^{\phi}_{g_{i}}}^{2}\geq{V^{\phi}_{g_{0}}}^{2}\cdot(1-z^{\phi})$
(10d) $\displaystyle V^{\phi}_{g_{r}}\cdot
V^{\phi}_{g_{i}0}-V^{\phi}_{g_{r}}\cdot V^{\phi}_{g_{z}^{\phi}}=0$ (10e)
$\displaystyle V^{\phi}_{g_{r}}\cdot V^{\phi}_{g_{r}0}\geq 0$ (10f)
$\displaystyle V^{\phi}_{g_{i}}\cdot V^{\phi}_{g_{i}0}\geq 0$ (10g)
$\displaystyle\sum_{\phi\in\Phi}\left(V^{\phi}_{g_{r}}\cdot
I^{\phi}_{g_{r}}+V^{\phi}_{g_{i}}\cdot I^{\phi}_{g_{i}}\right)=P$ (10h)
$\displaystyle\sum_{\phi\in\Phi}\left(V^{\phi}_{g_{i}}\cdot
I^{\phi}_{g_{r}}-V^{\phi}_{g_{r}}\cdot I^{\phi}_{g_{i}}\right)=Q$ (10i)
$\displaystyle P+Q\leq S_{max}$ (10j) $\displaystyle 0\leq z^{\phi}\leq 1$
(10k)

##  IV Case Study System

The case study system used to analyze and verify the operation of the
developed inverter models is a 4-bus synthetic microgrid, with photovoltaic
solar inverter-interfaced distributed generation. Fig. 3 presents the one-line
diagram of the system; it can be accessed in the PMsP GitHub as
case3_balanced_pv.dss[13].

The grid-following and grid-forming inverters were tested independently. Each
was installed on the PV Bus, and line-to-line, line-to-ground, and 3-phase-to-
ground faults were applied to the Load Bus, which is electrically close to the
inverters. The grid-following inverter was set to have a rated power of 25
kVA, while the grid-forming inverter to 100 kVA; the reason for the difference
is to ensure that the grid-forming inverter’s current was not constantly
saturated, trying to energize the circuit and feed the fault, since in the
tests the line between the substation (Source Bus) and the Primary Bus was
opened to create an island.

Figure 3: One-line diagram of the 4-bus case study system.

##  V Results and Discussion

The fault results of the developed inverter models were compared to the
experimental results presented in [3, 4, 5, 6, 7, 8]. To evaluate these new
inverter models, sweeps of fault impedance were performed to produce plots of
fault currents, voltages, and active and reactive powers versus the fault
impedances.

###  V-A Grid-Following Inverter Model

The performance of the grid-following inverter model under line-to-line,
single line-to-ground, and 3-phase faults was evaluated; the results are
presented in Figs. 4-6.

Figure 4: Grid-Following Inverter – Fault Currents [$A$] For Line-Line Fault.

Figure 5: Grid-Following Inverter – Fault Currents [$A$] For 3-Phase Fault.

Figure 6: Grid-Following Inverter – Fault Currents [$A$] For Line-Ground
Fault.

It can be observed that the inverter current saturates when the voltage
collapses, resulting in a drop in the supplied active power from the inverter;
the unfaulted phases, however, maintain the active power output. The
experimental results presented in [5] showed that most of the sequence
currents injected by grid-following inverters were positive sequence currents,
which is in agreement with the sequence currents supplied by the developed
grid-following inverter model.

Fig. 7 demonstrates that the solver returned a feasible solution for each
inverter over the faulted impedance sweep.

Figure 7: Feasibility of the Grid-Following Inverter Under Faults.

###  V-B Grid-Forming Inverter Model

The performance of the grid-forming inverter models under line-to-line, single
line-to-ground, and 3-phase faults were evaluated; the results are presented
in Fig. 8–13.

Figure 8: Simple Grid-Forming – Fault Currents [$A$] For Line-Line Fault.

Figure 9: Simple Grid-Forming – Fault Currents [$A$] For 3-Phase Fault.

The voltage-current relationship for the grid-forming inverters during and
outside saturation match the relationship presented in [7]. The experimental
results presented in [6] showed the sequence currents supplied by the grid-
forming inverters for various fault types, but provided little to no
information about the injected powers during faulted conditions. Based on the
sequence current results, the developed simple and complex grid-forming
inverter models appear to provide comparable sequence current injections.
However, it can be observed that unlike the developed grid-following inverter
model, the grid-forming models fail to be feasible for certain fault types.

Figure 10: Simple Grid-Forming – Fault Currents [$A$] For Line-Ground Fault.

The simple grid-forming inverter model failed to provide feasible solutions
for all fault impedances in the line-line fault (Fig. 8). This likely occurred
because the source was too stiff, as in [6] it was observed that the angles of
the faulted phases jumped by 40 degrees. This model was incapable of allowing
the voltage source to adjust enough to handle the fault, which can be seen in
the jumps in voltage in Fig. 8.

Figure 11: Complex Grid-Forming Fault Currents [$A$] For Line-Line Fault.

Figure 12: Complex Grid-Forming Fault Currents [$A$] For 3-Phase Fault.

Figure 13: Complex Grid-Forming Fault Currents [$A$] For Line-Ground Fault.

The complex grid-forming inverter model failed to provide feasible solutions
for all fault impedances in the 3-phase fault (Fig. 12). From the results, the
voltages and currents do not seem to exhibit any harsh jumps and do not seem
to indicate an obvious issue resulting from the voltage or current
constraints. It appears that the power constraints are the source of
infeasibility, which can be seen by comparing Fig. 9 and Fig. 12.

It must be noted that the voltages in the line-to-line fault exhibited similar
jumps as was seen in the simple grid-forming inverter model, but unlike that
model, the complex model was able to allow the voltage magnitudes to decrease
more to allow the solution to be feasible. The jumps in the voltage magnitudes
can be attributed to the voltage magnitude constraint, which is applied to
both faulted phases, not establishing a preference on which phase should be
lower than the other during the fault, and hence both solutions are feasible.

##  VI Conclusions

The introduced continuous variable inverter models for balanced and unbalanced
distribution networks are improved versions of the preliminary models used in
PMsP for pf, opf, and fault studies. The developed grid-following inverter
model presented in this work is able to provide feasible and reasonable
results for a variety of common fault types. On the other hand, the developed
simple and complex grid-forming inverter models still require future work to
improve behavior: to correct the infeasible results for line-to-line and
3-phase faults observed in the simple and the complex models, respectively.
From the simulation results, the modeling of the inverter as a Thevenin
equivalent source did not allow the voltage angles to adjust enough to allow
for a feasible solution. Usually, Thevenin sources are capable of sourcing a
large amount of current to allow for the proper voltages at the terminals.
This is not possible with the grid-forming inverter models because of current
saturation especially at low fault impedances. The complex model is capable of
adjusting the terminal voltage of the inverter to compensate for current
saturation, but has a limit on the amount of power that can be injected into
the system which caused infeasibility on the transition between saturated and
non-saturated current injection in the 3-phase fault. The power limit was not
an issue with the simple model because there was no constraint on power
injection with a Thevenin source. Future work will look on combining the
simple and complex model to overcome the limitation of the two models.

This work has laid down the foundation for developing relaxed models for both
grid-following and grid-forming inverter models for use in optimization
problems. The use of continuous activation variables (to indicate inverter
current saturation) over the use of mixed-integer formulation will allow for
the inverter models to be integrated into optimization problems on larger,
more complex distribution networks.

##  VII Acknowledgement

This work was performed with the support of the U.S. Department of Energy
(DOE) Office of Electricity (OE) Microgrid Research and Development (MRD)
Research Program under program manager Dan Ton. We gratefully acknowledge
Dan’s support of this work. The research work conducted at Los Alamos National
Laboratory is done under the auspices of the National Nuclear Security
Administration (NNSA) of the U.S. Department of Energy (OE) under Contract No.
89233218CNA000001. Approved for public release; distribution is unlimited –
LA-UR-22-32454.

## References

  * [1] J. A. Martinez et al.  Impact of Distributed Generation on Distribution Protection and Power Quality.  In Proc. of the 2009 IEEE Power & Energy Society General Meeting, pages 1–6, Jul. 2009.
  * [2] M. Meskin et al.  Impact of Distributed Generation on the Protection Systems of Distribution Networks.  IET Generation, Transmission & Distribution, 14(24):5944–5960, Nov. 2020.
  * [3] A. Mishra et al.  Fault Current Characterisation of Single Phase Inverter Systems.  In Proc. of the 2017 IEEE Power & Energy Society General Meeting, pages 1–5, 2017.
  * [4] D. Duckwitz et al.  Experimental Short-Circuit Testing of Grid-Forming Inverters in Microgrid and Interconnected Mode.  In Proc. of the 2018 Conference on Sustainable Energy Supply and Energy Storage Systems, pages 1–6, 2018.
  * [5] G. Kou et al.  Fault Characteristics of Distributed Solar Generation.  IEEE Transactions on Power Delivery, 35(2):1062–1064, Apr. 2020\.
  * [6] N. S. Gurule et al.  Experimental Evaluation of Grid-Forming Inverters Under Unbalanced and Fault Conditions.  In Proc. of the 46th Annual Conference of the IEEE Industrial Electronics Society, pages 4057–4062, 2020.
  * [7] N. S. Gurule et al.  Grid-forming Inverter Experimental Testing of Fault Current Contributions.  In Proc. of the 2019 IEEE 46th Photovoltaic Specialists Conference, pages 3150–3155, 2019.
  * [8] J. Hernandez-Alvidrez et al.  PV-Inverter Dynamic Model Validation and Comparison Under Fault Scenarios Using a Power Hardware-in-the-Loop Testbed.  In Proc. of the 2018 IEEE 7th World Conference on Photovoltaic Energy Conversion, pages 1412–1417, 2018.
  * [9] A. K. Barnes et al.  Optimization-Based Formulations for Short-Circuit Studies with Inverter-Interfaced Generation in PowerModelsProtection.jl.  MDPI Energies, 14(8):1–27, Apr. 2021.
  * [10] J. E. Tabarez et al.  Fault Current-Constrained Optimal Power Flow on Unbalanced Distribution Networks.  In Proc. of the 2022 IEEE ISGT ASIA Conference, pages 1–5, Nov. 2022.
  * [11] C. Coffrin et al.  PowerModels.jl: An Open-Source Framework for Exploring Power Flow Formulations.  In Proc. of the 2018 Power Systems Computation Conference, pages 1–8, Jun. 2018.
  * [12] D. M. Fobes et al.  PowerModelsDistribution.jl: An Open-Source Framework for Exploring Distribution Power Flow Formulations.  Electric Power Systems Research, 189(106664):1–7, 2020.
  * [13] PowerModelsProtection.jl/test/data/dist at master · lanl-ansi/PowerModelsProtection.jl.


# Constraining Energy-dependent emissivity profiles of AGN inflows

Dominic I. Ashton,1 & Matthew J. Middleton
1Department of Physics & Astronomy, University of Southampton, Southampton,
SO17 1BJ, UK

E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

The emissivity of the accretion flow is a key parameter affecting the shape of
both the energy and variability power spectrum of AGN. We explore the energy-
dependence of the power spectrum for five AGN, across the XMM-Newton bandpass,
and across the 0.01-1 mHz frequency range, finding a ubiquitous flattening of
the power spectrum towards higher energies. We develop a framework to explore
this behaviour and thereby extract the energy dependence of the emissivity
assuming a simple disc-like geometry for the inflow. We find that the
emissivity ranges from $\mathchoice{\mathrel{\raise 1.07639pt\hbox{\hbox
to0.0pt{\hbox{$\displaystyle\propto$}\hss}\lower
4.03563pt\hbox{$\displaystyle\sim$}}}}{\mathrel{\raise 1.07639pt\hbox{\hbox
to0.0pt{\hbox{$\textstyle\propto$}\hss}\lower
4.03563pt\hbox{$\textstyle\sim$}}}}{\mathrel{\raise 0.75346pt\hbox{\hbox
to0.0pt{\hbox{$\scriptstyle\propto$}\hss}\lower
2.82494pt\hbox{$\scriptstyle\sim$}}}}{\mathrel{\raise 0.5382pt\hbox{\hbox
to0.0pt{\hbox{$\scriptscriptstyle\propto$}\hss}\lower
2.0178pt\hbox{$\scriptscriptstyle\sim$}}}}~{}R^{-2}$ at energies around the
soft excess and increases to $\mathchoice{\mathrel{\raise 1.07639pt\hbox{\hbox
to0.0pt{\hbox{$\displaystyle\propto$}\hss}\lower
4.03563pt\hbox{$\displaystyle\sim$}}}}{\mathrel{\raise 1.07639pt\hbox{\hbox
to0.0pt{\hbox{$\textstyle\propto$}\hss}\lower
4.03563pt\hbox{$\textstyle\sim$}}}}{\mathrel{\raise 0.75346pt\hbox{\hbox
to0.0pt{\hbox{$\scriptstyle\propto$}\hss}\lower
2.82494pt\hbox{$\scriptstyle\sim$}}}}{\mathrel{\raise 0.5382pt\hbox{\hbox
to0.0pt{\hbox{$\scriptscriptstyle\propto$}\hss}\lower
2.0178pt\hbox{$\scriptscriptstyle\sim$}}}}~{}R^{-4}$ or steeper above 4-6 keV.
We describe the changing emissivity index with a linear function in energy,
finding the best-fitting slopes to vary between AGN. We attempt to correlate
the slope of the linear function against key AGN parameters but, as yet, the
sample size is too small to confirm hints of a correlation with Eddington
ratio.

###### keywords:

methods: data analysis – methods: statistical – galaxies: active – galaxies:
Seyfert – X-rays: galaxies

††pubyear: 2022††pagerange: Constraining Energy-dependent emissivity profiles
of AGN inflows–Constraining Energy-dependent emissivity profiles of AGN
inflows

##  1 Introduction

Studying the variability of AGN on short ($<$ Ms) timescales provides a means
by which to probe the otherwise unresolvable geometry of the inner accretion
flow (Fabian et al. 2009; Wilkins & Fabian 2013; Kara et al. 2019; Alston et
al. 2020), search for quasi-periodic oscillations (QPOs) resulting from plasma
dynamics or instabilities (Gierliński et al. 2008; Middleton & Done 2010;
Alston et al. 2014), and explore scaling relations for the SMBH mass (e.g.
Papadakis 2004; McHardy et al. 2006). Where data quality permits, the energy-
dependence of the variability can provide even greater insight via, e.g.
energy-lag spectra (Kara et al. 2019), deeper QPO searches (Ashton & Middleton
2021), rms and covariance spectra (Middleton et al. 2009; Middleton et al.
2011; Uttley et al. 2014).

It has been established that the shape of the power spectrum’s broad-band
noise is formed as a result of propagation of mass accretion rate fluctuations
(Lyubarskii 1997; Churazov et al. 2001; Ingram & Done 2011; Turner & Reynolds
2021), and is energy-dependent due to the sampling of different spatial
scales, differing emissivities (Arévalo & Uttley 2006), and the different
processes dominating various energy bands (e.g. the soft excess versus bands
in which the intrinsic coronal emission dominates). Studies of the energy-
dependent shape of the power spectrum have been rather limited, partly due to
the requirement for high signal-to-noise data and long observations (to probe
down to low frequencies and well-constrain the shape of the power spectrum).
When studies have been performed, they are usually restricted to a single AGN
of interest, e.g. NGC 7469 (in three energy bands: Nandra & Papadakis 2001),
MCG-6-30-15 (in three energy bands: Vaughan et al. 2003a; McHardy et al.
2005), NGC 4051 (in two energy bands: McHardy et al. 2004, and in three energy
bands using concatenated observations totalling $\sim 570$ ks: Vaughan et al.
2011), and most recently in RE J1034+396 (in three energy bands: Jin et al.
2020). In all of these previously published cases, the observed trend is that
the index of the power-law describing the power spectrum becomes flatter with
increasing energy (i.e. there is relatively more variability power at higher
frequencies at higher energies).

In this paper we explore in greater detail, the energy-dependence of the power
spectral shape for a sample of AGN. Ignoring reflection, we successfully model
the changing power spectral index with energy as a result of energy-dependent
emissivities and a low pass filter (Arévalo & Uttley 2006). The paper is
structured as follows: in Section 2 we discuss the AGN sample and the data
reduction methods. In Section 3 we describe the modelling of the power
spectra. In Section 4 we present our results, showing the energy-dependent
nature of these power spectra, across multiple AGN. In Section 5 we present
our theoretical framework and in Section 6, we discuss the limitations of our
approach and highlight the wider relevance and impact of our results.

##  2 Data analysis

###  2.1 AGN Sample and Data reduction

The data and extraction methods we have used are presented in Ashton &
Middleton (2021) in which we consider a sample of 38 bright AGN and a total of
200 XMM-Newton observations. The combination of XMM-Newton with bright AGN
results in high signal-to-noise data and typically long (up to 130 ks)
exposures. The source sample is drawn from the Palomar Green QSO (PGQSO)
catalogue as used by Crummy et al. (2006) and Middleton et al. (2007), with a
number of additional, bright and well-observed AGN: 1H 0707-495, MS
22549-3712, NGC 4151, PHL 1092 and IRAS 13224-3809. The total sample consists
of 22 NLS1s and 16 Type-1 Seyferts, as defined by Veron-Cetty & Veron (2010).

We considered all observations in the XMM-Newton public HEASARC
archive111https://heasarc.gsfc.nasa.gov/ up to September 2019, and follow the
same data extraction method described in full in Ashton & Middleton (2021),
which we do not replicate here for the sake of brevity. Of importance is that
we choose to use only the PN (rather than PN and MOS), and follow standard
procedure by selecting only single/double pixel events (with pattern $\leq$
4), before using the standard XMM-Newton data reduction pipeline
222https://www.cosmos.esa.int/web/xmm-newton/sas-thread-timing in SAS v17.0.0.
We use evselect to extract source and background light curves (from 40"
regions) with a binsize of 100s and a sliding energy window across the nominal
$0.3-10$ keV range of the instrument. This window is broken into 50 non-
overlapping energy bins, allowing for 1225 different energy combinations to be
studied (e.g. $0.3-0.5$ keV, $0.3-0.7$ keV etc.). We then proceed to use
epiclccorr to subtract the background and apply the necessary corrections to
each of the light curves.

Times corresponding to flares in the high-energy ($10-12$ keV), full-field
background were excised (see Ashton & Middleton 2021 for details). Such flare
subtraction can introduce gaps into the lightcurves, which, if left
unmodified, can affect the shape of the power spectrum. Standard approaches
include linear interpolation (González-Martín & Vaughan 2012; Alston et al.
2015) but in our analysis, we take the most conservative approach and consider
only the longest continuous segment between flares. Whilst this has the
potential to restrict the low frequency end of the power spectrum, we find
that many of the longest continuous segments in our sample allow us to reach
down to $\sim$0.01 mHz.

##  3 AGN Power Spectra

Figure 1: Example energy-resolved periodograms of IRAS 13224-3809
(OBSID:0780561601), fitted using a maximum likelihood plc model, with best-
fitting $\beta$ values illustrated. The decreasing number of frequency bins
reflects the increasing levels of white noise with increasing energy.

We calculate the periodogram for each energy-resolved lightcurve in
fractional-rms units (Vaughan et al. 2003b) – examples of which are shown in
Figure 1 for the case of an observation of IRAS 13224-3809. We analyse the
periodogram rather than the PSD in order to maximise the number of Fourier
frequencies and access the lowest available frequencies in the lightcurve. To
ensure sufficient data for modelling, we restrict our analysis to those
periodograms containing $\geq 13$ frequency bins at frequencies below the
point where white noise starts to dominate (Ashton & Middleton 2021).

Similar to González-Martín & Vaughan (2012), we fit each periodogram using a
maximum likelihood statistic and determine the preferred description by either
a power-law + constant (plc) or broken power-law + constant model (bknplc).
The plc takes the form of:

$P(\nu)=N_{0}\nu^{\hskip 1.84949pt\beta}+C$ (1)

which consists of free parameters: $\beta$, the spectral index, and $N_{0}$,
the normalisation, as well as $C$, a constant to account for Poisson noise
(which we set to the predicted level of the white noise following the method
described in Vaughan et al. 2003b). The bknplc model is:

$P_{1}(\nu)=N_{1}\nu^{\hskip 1.84949pt\beta_{1}}+C\hskip
14.22636pt\textnormal{for }\nu\leq\nu_{\rm b}$ $P_{2}(\nu)=N_{2}\nu^{\hskip
1.84949pt\beta_{2}}+C\hskip 14.22636pt\textnormal{for }\nu>\nu_{\rm b}$ (2)

At frequencies at, or below the break ($\nu\leq\nu_{\rm b}$), we assume a
spectral index $\beta_{1}$ and normalisation $N_{1}$, and at frequencies above
the break ($\nu>\nu_{\rm b}$), we assume a spectral index $\beta_{2}$. The
normalisation above the break, $N_{2}$, is calculated directly from the
previous parameters. A Bayesian Information Criterion (BIC) test (Liddle 2007)
is applied to determine the statistically preferred model. As reported in
Ashton & Middleton (2021), $\sim 10\%$ of the total sample of power spectra
which match our data-quality criteria are preferentially modelled with a
broken power-law, a proportion in general agreement with González-Martín &
Vaughan (2012). Indeed, this 10% includes the re-identification of break
features in observations first reported in González-Martín & Vaughan (2012).
To avoid the complicating effects of modelling a break in the periodogram,
energy bands where a break is identified are excluded from our subsequent
analysis.

In Ashton & Middleton (2021), $7$ AGN in our sample of $38$ were found to show
statistically significant QPO-like features in their power spectra, when the
broadband noise was fitted with a (statistically preferred) plc model. We
therefore take the additional precaution of identifying the Fourier bin
corresponding to the most prominent outlier in power above the maximum
likelihood plc model, removing it, and refitting using a bootstrap-with-
replacement method. This approach minimises the bias from an individual,
outlier bin on the shape of the power spectrum.

###  3.1 Energy-Dependent Power Spectra

We explore the energy-dependent shape of the best-fitting plc model for each
AGN in turn. As we lack the necessary data quality to examine the energy
dependence of the power spectrum at the highest resolution of our sliding
energy window, we obtain our maximum likelihood plc model parameters in five
larger energy bins: $E_{1}=0.3-1.0$ keV, $E_{2}=1.0-2.0$ keV, $E_{3}=2.0-4.0$
keV, $E_{4}=4.0-6.0$ keV and $E_{5}=6.0-10.0$ keV. For each AGN, each
observation, and in each energy bin, we obtain the mean of the power spectral
index, $\bar{\beta}$ and its 1$\sigma$ error (we are averaging over a
sufficient number of observations in each case to provide quasi-Gaussian
errors).

We search across all $38$ AGN in our sample, and report those cases where we
have constrained values for $\bar{\beta}$ across a minimum of four energy
bins. Due to diminishing data quality at high energies, we are unable to
constrain values for $\bar{\beta}$ as we approach 10 keV for the majority of
our sample. It is unsurprising that, of the five AGN for which we can
constrain $\bar{\beta}$ values across a minimum of four energy bins (IRAS
13224-3809, 1H 0707-495, MRK 766, NGC 4051 and ARK 564), all are bright and
amongst the most well-studied AGN by XMM-Newton.

##  4 Results

In Figure 2, we show how $\bar{\beta}$ changes with energy across the five AGN
mentioned above. For each AGN, we see that the energy dependence of
$\bar{\beta}$ follows an approximately linear trend, with less negative
$\bar{\beta}$ values, i.e. a flattening of the power spectrum, as we approach
higher energies. We proceeded to model the change in index with energy,
$\frac{d\bar{\beta}}{dE}$, for each AGN, with a simple linear model using a
widely available Markov-chain Monte Carlo (MCMC) routine (emcee). The high
data quality of IRAS 13224-3809 and NGC 4051 allows us to constrain the mean
index in the highest energy bin, but for ARK 564, MRK 766 and 1H 0707-495, we
are limited to a value from only a single observation in $E_{5}$. We highlight
these single points in Figure 2, but do not use them to constrain our linear
model. For the MCMC fitting, we assume uniform priors, and implement a routine
with $100$ walkers over 10,000 iterations, with a discarded burn-in time of
the first $2000$ steps. We return the median of each parameter’s posterior
distribution, along with the percentiles required to highlight $68\%$ credible
intervals about the MCMC fit, as shown in Figure 2. The best-fitting
$\frac{d\bar{\beta}}{dE}$ values returned by this method are provided in Table
LABEL:tab:LinearModels, and the total number of observations contributing to
each energy bin for each AGN is indicated in Figure 2.

As can be seen from Figure 2 and Table LABEL:tab:LinearModels, IRAS
13224-3809, ARK 564 and MRK 766 all have similar values of
$\frac{d\bar{\beta}}{dE}$, between $0.15$ keV-1 and $0.16$ keV-1 which differ
compared to NGC 4051 and 1H 0707-495, which have $\frac{d\bar{\beta}}{dE}$ of
$0.08~{}\pm 0.01$ keV-1 and $0.22~{}\pm 0.02$ keV-1 respectively.

Although we have limited confidence in the values for $\bar{\beta}$ above 6
keV for most of our AGN, it is interesting to note that there is a tendency
for those values to lie below the linear model, i.e. at steeper indexes than a
linear model would predict. It is unclear if this is real but may be confirmed
in future using NuSTAR. We also note that in all cases, the very softest index
lies below the linear relation, indicating that the true description of the
energy-dependence likely requires a steeper rise before flattening to higher
energies.

Table 1: Summary of the best fitting $\frac{d\bar{\beta}}{dE}$ values shown in Figure 2. Columns indicate: (1) AGN, (2) The gradient of the best-fit linear model, with associated 68% credible interval. AGN | $\frac{d\bar{\beta}}{dE}$
---|---
IRAS 13224-3809 |  $0.16~{}\pm^{0.02}_{0.03}$ keV-1
1H 0707-495 |  $0.22~{}\pm^{0.02}_{0.02}$ keV-1
NGC 4051 |  $0.08~{}\pm^{0.01}_{0.01}$ keV-1
ARK 564 |  $0.16~{}\pm^{0.04}_{0.04}$ keV-1
MRK 766 |  $0.15~{}\pm^{0.04}_{0.05}$ keV-1
Figure 2: Energy-dependence of the mean power spectral index, $\bar{\beta}$,
when modelled with a PLC across the five AGN for which we have constraints in
a minimum of four energy bins. We illustrate the maximum likelihood fit and
68% credible intervals from the MCMC based on a linear model. We show the
number of observations that are averaged in order to calculate $\bar{\beta}$
in each energy bin, along the bottom of each panel. In those cases where we
are limited to a single observation in an energy bin, we plot the single
measurement of $\beta$ (indicated by the filled triangles) and do not use this
to contribute to the linear fit due to insufficient statistics.

###  4.1 Excluded AGN

We briefly mention notable AGN where constraints were not possible for $\geq$
four energy bins. We find that, in many cases, our AGN power spectra can only
be constrained at soft energies, typically $<2$ keV or across a small number
of observations in total. An example of this is the case of RE J1034+396,
which we find to only have sufficient data to meet our selection criteria
below $2$ keV (due to its spectrum being dominated by the soft excess, see
Puchnarewicz et al. 2001; Casebeer et al. 2006; Middleton et al. 2009).

Both PG 1211+143 and RE J1034+396 show the start of a positive linear trend
between index and energy (as in Figure 2) but there are only constraints based
on multiple observations in $E_{1}$ and $E_{2}$ for RE J1034+396, and, for PG
1211+143, we have constraints from only a single observation. PKS 0558-504 is
similar in appearing to show a weak positive linear trend across three bins
but we only have a single observation contributing to the $E_{3}$ bin and none
at higher energies for this AGN.

In the case of PG 0003+199, we observe a strong linear trend across four
energy bins ($E_{1}-E_{4}$), comparable to our AGN in Figure 2, but in this
case we are limited to a single observation throughout. Considering only this
one observation of PG 0003+199, we would infer $\frac{d\bar{\beta}}{dE}=0.35$,
which is even larger than in the case of 1H 0707-495. However, as with only a
single observation we cannot obtain reliable errors for
$\frac{d\bar{\beta}}{dE}$, we merely highlight this AGN as worthy of future
investigation.

Lastly we report two cases of statistically limited, apparent negative or non-
linear trends in PG 1244+026 and PG 1322+659. Both show a flattening in index
between $E_{1}$ and $E_{2}$, before a decrease in index in $E_{3}$ to values
which are even steeper than in the lower energy bins. In both cases however,
these trends are derived from only a single observation.

##  5 Theoretical Framework

The spectra of NLS1s tend to be a complex combination of components (Jin et
al. 2012; Done et al. 2013), including a soft excess created by partially
ionised reflection (Crummy et al. 2006, Done & Nayakshin 2007) and/or Compton
upscattering of seed disc photons in a cool corona (Middleton et al. 2009;
Done et al. 2012), a power-law from a hot corona and its reflection from the
disc, producing – in addition to a soft excess – a characteristic Fe Kα line
and Compton hump (George & Fabian 1991). Imprinted onto this continuum are
lines due to photo-ionised absorption by outflows (Tombesi et al. 2010;
Tombesi et al. 2011; Pinto et al. 2018). Each component in the spectrum can
vary on a range of timescales and are interlinked either through direct
propagation (which produces the canonical rms-flux relationship: Uttley &
McHardy 2001, and lag to hard energies at low frequencies: Walton et al.
2013), or via irradiation (which can yield a soft lag at high frequencies:
Zoghbi et al. 2012; Kara et al. 2016).

Any model which aims to describe the emergent variability must account for the
interplay between these components on various timescales. Here, for the sake
of simplicity we focus on only the propagation and reflection components,
which will be the most salient to our study.

###  5.1 Reflection

Reflection is a well-studied feature of AGN energy spectra (e.g. George &
Fabian 1991) which not only produces the well known Fe line and Compton hump,
but can also provide some or all of the soft-excess (e.g. Fabian et al. 2004).
Due to being irradiated by an intrinsically variable source of seed photons
(the corona), the reflected emission will also vary, with the details of this
depending on the energy- and frequency-dependent transfer function (which
describes the co-location of the system components in curved space-time).

The potential impact of reflection on the shape of the AGN power spectrum has
been investigated by Papadakis et al. (2016), Emmanoulopoulos et al. (2016)
and Chainakun (2019). Initial investigations focussed on single lamp-post
configurations (Papadakis et al. 2016), although no signatures of
reverberation could be located in the power spectrum (Emmanoulopoulos et al.
2016). Most recently, Chainakun (2019) explored how a vertically extended
(two-blob) corona can lead to predictable changes due to reverberation,
utilising ray tracing in full GR to obtain the response functions. In all
cases – and as expected – the more extended the corona (in height), the lower
the frequency at which the PSD is predicted to change. However, these changes
are expected to be small, even where the energy band is dominated by
reflection.

Although the effect of reverberation/reflection on the power spectrum has not
been calculated for a radially extended corona (e.g. where the corona is a
sandwich type geometry or is the hot innermost flow), the effect can be well
reasoned. In all models (regardless of geometry), and for any chosen energy
band, $E_{A}\rightarrow E_{B}$, the observed lightcurve is a summation
$s(E,t)=p(E,t)+r(E,t)$, where $p(E,t)$ represents emission from the
propagation (coronal) component, and $r(E,t)$ from reflection. Assuming that
$r(E,t)$ is linearly correlated with $p(E,t)$, we can relate them through an
impulse response function $\phi(E,t)$:

$r(E,t)=\phi(E,t)\circledast p(E,t)$ (3)

and its Fourier transform (and transfer function $\Phi(\nu)$):

$R(E,\nu)=\Phi(E,\nu)P(E,\nu)$ (4)

These components may be expressed in terms of their complex Fourier amplitudes
in the standard way:

$\displaystyle|P|^{2}=P^{*}P$ (5) $\displaystyle|R|^{2}=R^{*}R=(\Phi
P)^{*}(\Phi P)=|\Phi|^{2}|P|^{2}$ (6)

The final power spectrum we observe is then a weighted sum of these components
with distortion due to reflection encoded in the impulse response function,
with the lowest frequency effects occurring for the largest light travel times
between corona and reflector (regardless of whether the response is an
approximation or is ray-traced: Chainakun 2019). In a disc-like coronal
geometry (which we will be considering here), the light-travel time between
seed photons and reflector is small and so any distortion is similarly small
(certainly across the $0.01-1$ mHz frequency range we are investigating here).
Given the lack of evidence for reverberation/reflection induced changes to the
AGN power spectrum, we will consider only the effects of propagation
hereafter.

###  5.2 Propagation

The propagation of material through the accretion flow imprints variability
across a broad range of radii (and therefore timescales) by coupling to local
mass accretion rate fluctuations (see Lyubarskii 1997; Kotov et al. 2001).
Arévalo & Uttley (2006) explored the effect of propagation in acting as a low-
pass filter on the shape of the intrinsic power spectrum of the AGN, a process
which suppresses high frequency variability. This low-pass filter is described
as:

$\displaystyle\textrm{PSD}_{\rm filtered}(\nu)=\textrm{PSD}_{\rm
intrinsic}(\nu)\left(\frac{\int^{r_{\nu}}_{r_{\rm min}}\epsilon(r)2\pi
rdr}{\int^{\infty}_{r_{\rm min}}\epsilon(r)2\pi rdr}\right)^{2}$ (7)

where $\epsilon(r)$ is the radial emissivity profile of the flux. From Arévalo
& Uttley (2006), $\epsilon(r)$ is given by:

$\displaystyle\epsilon(r)=r^{-\gamma}\left(1-\sqrt{\frac{r_{\rm
min}}{r}}\right)$ (8)

with emissivity index $\gamma$, from a minimum radius $r_{{\rm min}}$ out to
arbitrary radius $r$. A commonly assumed value for the emissivity index is
$\gamma=3$ (i.e. that of a geometrically thin disc). The value of $\gamma$ may
vary broadly with energy, from a suggested minimum of $\gamma\sim 2$ for low
energies up to $\gamma>5$ at the highest energies in the XMM-Newton bandpass
(Arévalo & Uttley 2006). Assuming the input power spectrum takes the form of a
power-law we can therefore explore the impact of propagation on the changing
index of the filtered power spectrum.

Figure 3: Simulated power spectra for a set of pre-defined emissivity indices,
$\gamma$(E), black hole spins, $a_{*}$ and $\alpha(\frac{H}{R})^{2}$ values,
assuming propagation determines the shape of the broadband noise (distorted
from an intrinsic index of $-1$). We show the frequency bands of interest
between $0.01-1$ mHz.

###  5.3 Modelling propagation

In order to utilise a radial emissivity, as in the model of (Arévalo & Uttley
2006), we are explicitly assuming a model where the corona is embedded in the
disc as a radial structure (see e.g. the model of Done et al. 2012).

We assume the timescales we observe here are generated in the corona itself
(and note that longer timescale changes will be driven by the disc outside of
the corona) and are made visible by a changing electron temperature and
density as a result of local turbulence and coupling to mass accretion rate
propagations through the flow from $r_{\rm out}$ to $r_{\rm in}$. We note that
we can discount variations in the seed photon population as the driver for
changes in the power spectrum, as the UV emission from the disc is established
to be considerably less variable than the corona in NLS1s (Leighly 1999; Smith
& Vaughan 2007; Ai et al. 2013; Alston et al. 2013; Done & Jin 2016). We
assume that the variability generated locally at each radius ($r_{\nu}$) is at
the viscous frequency (see Churazov et al. 2001) such that:

$\displaystyle
r_{\nu}=\left[\frac{2\pi\nu}{\alpha}\left(\frac{H}{R}\right)^{-2}\right]^{-2/3}$
(9)

where the frequency is in units of $c/R_{g}$ (see e.g. Kato et al. 1998;
Arévalo & Uttley 2006). In the above, $\alpha$ and $\frac{H}{R}$ are the
dimensionless viscosity parameter of Shakura & Sunyaev (1973), and scale
height of the accretion disk respectively. We assume that our frequency range
of interest, $0.01-1$ mHz (i.e. the range over which which we can practically
fit to the data) corresponds to radii between the ISCO ($r_{\rm in}$) and some
radius within the true outer edge of the corona (i.e. $r_{\rm out}\leq r_{\rm
corona}$). The actual frequencies generated in our model therefore depend on
the SMBH spin and the combination $(\frac{H}{R})^{2}\alpha$ (and somewhat on
the SMBH mass – although here the range is small). Given the reported high
spin values for these bright AGN (Ogle et al. 2004; Fabian et al. 2013; Done &
Jin 2016; Kara et al. 2017; Buisson et al. 2018b), we expect the ISCO to sit
at $\sim 1.25R_{\rm g}$. For the corona at the ISCO to produce variability
above our upper frequency limit of 1 mHz requires
$(\frac{H}{R})^{2}\alpha\gtrsim 0.01$. We note that for the mass range
subtended by our AGN sample (from $10^{6.00}-10^{6.63}M_{\odot}$, see Table
LABEL:tab:Gamma_Parameterisation), should we instead assume zero spin ($r_{\rm
in}=6R_{g}$), the viscous frequency at $r_{\rm in}$ is lower and we have
strong curvature in our observed bandpass.

As neither $\frac{H}{R}$ nor $\alpha$ are well determined in practice, in the
following we simply assume $(\frac{H}{R})^{2}\alpha=0.3$ to be consistent with
Arévalo & Uttley (2006), and which is appropriate for large scale height
inflows (which is a possibility in the inner regions of these high Eddington
ratio AGN) and/or large values of $\alpha$. A frequency of 0.01 mHz then
corresponds to an outer radius of $\sim 100R_{\rm g}$ (for a $10^{6}M_{\odot}$
SMBH). It should be recognised that this is sensitive to $\alpha(H/R)^{2}$ and
instead, at a value of 0.1, the outer radius is then $\lesssim 50R_{\rm g}$.
We will discuss the possible implications of such large truncation radii in
the Discussion.

We assume the intrinsic (unfiltered) power spectrum is well described by a
power law with index, $\beta=-1$ which is appropriate for frequencies above
any low frequency break (analogous to those detected in black hole XRBs, see
e.g. McHardy et al. 2006) but below the high frequency break at the ISCO
(Papadakis et al. 2002; Arévalo & Uttley 2006; McHardy et al. 2006). Following
Arévalo & Uttley (2006), and as explained above, we set $r_{\rm in}$ to be the
position of the ISCO with $r_{\rm in}=1.25R_{g}$, whilst we set $r_{\rm
out}=100R_{g}$ (although we reiterate that this need not be the actual extent
of the corona). Following equation 7, we can obtain the filtered power
spectrum in a given energy range assuming a value for the emissivity index
$\gamma(E)$. To illustrate this process, Figure 3 shows a range of simulated,
filtered power spectra, all produced from an intrinsic slope of $\beta=-1$,
assuming values of $(\frac{H}{R})^{2}\alpha$ including
$(\frac{H}{R})^{2}\alpha=0.3$ (following Arévalo & Uttley 2006), and using
various parameter combinations for $\gamma(E)$ and $a_{*}$. In practice, we
limit the output (filtered) power spectrum to a frequency range of $0.01-1$
mHz, corresponding to the range of frequencies typically accessible in our
observations and fit a power-law to each simulated power spectrum and obtain a
new value of the power-law index, $\beta$ (corresponding to the effect of
propagation).

As can be seen from Figure 3, we are able to obtain substantial energy-
dependent changes in $\beta$, obtaining less negative (flatter) values at
higher energies, similar to observation (Figure 2). For a given $\gamma(E)$,
we can therefore produce expectation values for $\bar{\beta}$, or conversely
use our observed $\bar{\beta}$ to constrain $\gamma(E)$ in this simplified
picture. As the exact profile of $\gamma(E)$ is uncertain, we parameterise it
in both a linear form:

$\gamma(E)=\frac{d\gamma}{dE}E+b$ (10)

and power law form:

$\gamma(E)=\gamma_{0}E^{\alpha}$ (11)

and determine the values of parameters $(\frac{d\gamma}{dE},b)$ and
$(\gamma_{0},\alpha)$ which maximise the likelihood function between our
observed values of $\bar{\beta}_{\rm obs}$ and those theoretical values
extracted via the above process $\bar{\beta}_{\rm mod}$. This is equivalent to
minimising the negative log-likelihood, $S\equiv-2\ln{[\mathcal{L}]}$ (Vaughan
2005). Assuming emissivity indices $\hat{\gamma}(E_{\rm k})$ across energy
bands k$=1,2,...,n-1$, $S$ may be written as:

$S=2\sum_{k=1}^{n-1}\left\\{\ln[\bar{\beta}_{\rm mod,\rm
k}]+\frac{\bar{\beta}_{\rm obs,\rm k}}{\bar{\beta}_{\rm mod,\rm k}}\right\\}$
(12)

We obtain similarly satisfactory values of the negative log likelihood for
both models, with a breakdown of our results provided in Table
LABEL:tab:Gamma_Parameterisation. Given the small statistical improvement
afforded by the linear model (the fit statistic $S_{\rm lin}\lesssim S_{\rm
PL}$ in most cases), we proceed to use this model for $\gamma(E)$ hereafter.
We plot a direct comparison between the data and simulated results for the
linear model in Figure 4.

Figure 4: A comparison of mean power spectral index values, $\bar{\beta}$, as a function of energy between our model (red) and observed values (black), across the five AGN. $68\%$ credible intervals are shown as in Figure 2. We determine the best-fitting linear model, $\gamma(E)=\frac{d\gamma}{dE}E+b$, and select values for $\gamma(E)$ that maximise the likelihood by comparing the theoretical to observed values of $\bar{\beta}$. Table 2: AGN properties and results from fitting for $\gamma(E)$, with maximum likelihood determined parameter values $(\frac{d\gamma}{dE},b)$ and $(\gamma_{0},\alpha)$, and associated negative likelihood values $S_{\rm lin}$ and $S_{\rm PL}$, for linear and power law models respectively, for each AGN in our study meeting our data criteria. Black hole masses, $M_{\rm BH}$ and bolometric luminosities, $L_{\rm bol}/L_{\rm Edd}$ (calculated using the mass estimate), as well as spin, $a_{*}$, and inclination angle $i$, are sourced from the literature: (a) Alston et al. (2019), (b) Buisson et al. (2018a); (c) Fabian et al. (2013); (d) Done & Jin (2016); (e) Dauser et al. (2012); (f) Denney et al. (2010), (g) Woo & Urry (2002); (h) Ogle et al. (2004); (i) Zhang & Wang (2006); (j) Vasudevan & Fabian (2007); (k) Kara et al. (2017); (l) Wang, T. & Lu, Y. (2001); (m) Vasudevan et al. (2010); (n) Buisson et al. (2018b). AGN | $\log_{10}(\frac{{}{M_{\rm BH}}}{M_{\odot}})$ | $L_{\rm bol}/L_{\rm Edd}$ | $a_{*}$ | $i(\rm deg)$ | $\frac{d\gamma}{dE}$ | $b$ | $S_{\rm lin}$ | $\alpha$ | $\gamma_{0}$ | $S_{\rm PL}$
---|---|---|---|---|---|---|---|---|---|---
IRAS 13224-3809 |  $6.0$ (a) | 2.96 (b) |  $0.989\pm 0.001$ (c) |  $62.3\pm 1.3$ (c) | 0.29 | 1.83 | 13.41 | 0.22 | 2.23 | 13.45
1H 0707-495 |  $6.3$ (d) | 1.18 (d) |  $0.998^{+0.000}_{-0.002}$ (e) |  $52.0^{+1.7}_{-1.8}$ (e) | 0.53 | 1.73 | 11.21 | 0.32 | 2.37 | 11.22
NGC 4051 |  $6.2$ (f) | 0.19 (g) |  $0.973^{+0.025}_{-0.000}$ (h) |  $48.0^{+4.0}_{-3.0}$ (h) | 0.11 | 1.99 | 15.54 | 0.11 | 2.16 | 15.55
ARK 564 |  $6.4$ (i) | 0.59 (j) |  $0.998^{+0.000}_{-0.000}$ (k) |  $22.0^{+7.0}_{-2.0}$ (k) | 0.35 | 2.11 | 11.57 | 0.29 | 2.36 | 11.57
MRK 766 |  $6.6$ (l) | 0.59 (m) |  $0.92^{+0.078}_{-0.000}$ (n) |  $46.0^{+1.0}_{-2.0}$ (n) | 0.36 | 2.31 | 12.25 | 0.22 | 2.75 | 12.24
Figure 5: The maximum-likelihood determined values for the emissivity index
function, $\gamma(E)$, calculated from the $\bar{\beta}$ values observed
across each AGN (as in Figure 4), each with simulated $\sim 68\%\ (1\sigma)$
confidence intervals. All calculations assume a linear parameterisation,
$\gamma(E)=\frac{d\gamma}{dE}E+b$.

As can be seen from Figure 4, our simple linear model for $\gamma(E)$ provides
a reasonable description for the changes in power-law index across all five
AGN, and performs extremely well in the case of ARK 564. The model is not
capable of providing an ideal description in each case however, with clear
systematic offsets visible in the cases of IRAS 13224-3809, 1H 0707-495 and
NGC 4051.

The differing nature of the best-fitting emissivity profiles for each AGN is
visible from Figure 5. The uncertainty on the $\gamma(E)$ values is evaluated
through a Monte-Carlo routine, whereby sets of $\beta$ values (equal in number
to the number of observations per energy band), are drawn for each AGN across
all energy bands (see Figure 2) which yields an individual MC value for
$\bar{\beta}(E)$, with this routine assuming Gaussian errors on the observed
values of $\bar{\beta}(E)$, exactly as in Figure 2. This process is repeated
$1000$ times and for each set of $\bar{\beta}(E)$ values, for each AGN, a
maximum likelihood set of $\gamma(E)$ values is determined which best describe
it. Figure 5 shows both the maximum likelihood set of $\gamma(E)$ values
determined from the observed power spectra and these simulated $68\%$
confidence intervals, where notably, the emissivity appears considerably
flatter with energy in NGC 4051 relative to IRAS 13224-3809 and 1H 0707-495 in
particular. We speculate as to the origin of such differences in the
Discussion.

##  6 Discussion

Modelling our sample of AGN indicates that, in all cases where the data allow
us to explore the changing shape of the power spectrum, there is a positive
correlation between power spectral index between $0.01-1$ mHz (where the
periodogram can be well described as a single power-law) and energy (across
the $0.3-10$ keV range), such that the index flattens with increasing energy.
Assuming a simple model of propagating fluctuations – which is observationally
supported by the presence of hard lags at low frequencies in AGN (Walton et
al. 2013), the ubiquitous linear rms-flux relationship (Uttley & McHardy
2001), and log-normal flux distribution (Gaskell 2004) – we have been able to
obtain corresponding constraints on the emissivity as a function of energy,
(discounting reflection – see Emmanoulopoulos et al. 2016).

Our model of an in-flow corona (e.g. Done et al. 2012) allows us to assume
that the inflow acts a low pass filter, weighted by the radial emissivity
profile. We have found that our model requires the emissivity to be around
$R^{-2}$ for energies around the soft excess before increasing to a maximum of
$\sim R^{-5}$. It is perhaps intriguing to note that radiation pressure or
advection dominated flows are expected to have flatter emissivity profiles
(with $\gamma<3$: e.g. Watarai et al. 2000) which may be consistent with a
high accretion rate and such a component producing some of the soft excess
(Middleton et al. 2009; Done et al. 2012). Alternatively, it may be that the
corona has a larger scale-height (perhaps due to radiation pressure), which
pushes the variability to higher frequencies (as the viscous timescale goes as
$(\frac{H}{R})^{-2}$), flattening the index to a fixed value above a certain
energy. Indeed, we can see from Figure 2 that the index appears to saturate
and possibly begin to inflect; whilst we presently lack the statistics to
constrain this inflection, with future observations using instruments with
higher energy coverage e.g. NuSTAR, we will be able to investigate this
further. Finally, whilst our analysis assumes an input power-law index of
$-1$, other values for the shape of the intrinsic noise are excluded by virtue
of producing incompatible results (section 5.3).

We find our best-fitting function for $\gamma(E)$ to vary across our AGN
sample (Figure 5). It is plausible that such differences could be due to
geometrical effects (e.g. different scale-heights) or some bulk property of
the AGN, notably accretion rate, mass and spin. We obtain the spins,
inclinations, black hole masses and Eddington scaled accretion rates using
$L_{\rm bol}$ from the literature – see Table
LABEL:tab:Gamma_Parameterisation. The spin and inclination values are derived
from use of reflection models which assume a lamp-post geometry and as such
are not consistent with our assumed geometry. By plotting these spin and
inclination values against our derived emissivity gradients we are instead
exploring whether an underlying degeneracy might exist. In addition,
statistical errors on each parameter are very likely to be underestimates
(e.g. see Reynolds 2014; Middleton et al. 2016); these values are therefore
used with due caution.

Figure 6: Key macro-properties for our AGN sample – spin, inclination,
Eddington ratio and BH mass plotted against $\frac{d\gamma}{dE}$ (see Table
LABEL:tab:Gamma_Parameterisation).

In Figure 6 we plot $\frac{d\gamma}{dE}$ versus the AGN parameters of interest
and find a lack of a correlation with spin, inclination or BH mass but
possible hints of a positive correlation with Eddington ratio. In examining
the Eddington ratio, there appears to be slightly stronger evidence of a
correlation, with the only outlier being IRAS 13224-3809 with a highly
uncertain value for $L_{\rm bol}/L_{\rm Edd}$ (see Buisson et al. 2018a). We
note that the Eddington ratio is also highly uncertain for 1H 0707-495 (Done &
Jin 2016); clearly both a considerably larger AGN sample is required as well
as a more accurate handle on the key macro-properties before any claim of a
correlation can be made.

One of the most important assumptions in our analysis is the geometry which
dictates how our variability is filtered and modified by an assumed radial
emissivity (c.f. Arévalo & Uttley 2006). Our assumption here of a continuous
disc-like inflow is in contrast to models where the corona sits above the
black hole and is probably the base of a jet (Wilkins & Gallo 2015). Such
lamp-post geometries lead to increased emissivities in the inner portions of
the disc – and therefore the reflected emission – due to light-bending
(Wilkins & Fabian 2011; Gonzalez et al. 2017). Caution is needed in assuming a
disc-like coronal geometry; whilst an extent of $\sim$ 50Rg would be broadly
consistent with micro-lensing studies (Mosquera et al. 2013), the magnitude of
hard lags may imply considerably smaller coronae. In addition, a radially
extended corona of the type we have utilised here would clearly be at odds
with relativistic broadening of Fe lines (Fabian et al. 2009) unless for some
reason the disc reappeared at radii closer to the ISCO.

Importantly, the variability of the corona (regardless of geometry) is subject
to other heating/cooling effects besides propagation which are beyond the
scope of this paper. Such a model may well explain the changing shape of the
power spectrum with energy without the need to invoke a changing emissivity.
Finally, in future it will be interesting to compare the changing shape of the
power-spectrum with energy in our AGN sample, to that seen in XRBs (in
particular the very high state).

##  7 Conclusion

The well-evidenced propagation model for accretion (Lyubarskii 1997; Kotov et
al. 2001; Ingram & Done 2011; Turner & Reynolds 2021) makes a clear prediction
that – where we assume a simple disc-like geometry for the corona – the shape
of the power spectrum we observe should be modified by the area of the inflow
producing the variability and the radial emissivity profile as a function of
energy (Arévalo & Uttley 2006). Modelling the energy-dependence of the power
spectrum should therefore provide interesting insights into the changing
radial emissivity under the conditions of the assumed geometry.

Following from the analysis of $38$ bright AGN in Ashton & Middleton (2021),
we model the energy-resolved periodograms using XMM-Newton data across a
frequency range of 0.01-1 mHz (at frequencies typically lower than the viscous
timescale at the ISCO for the high spins and low masses of the AGN in our
sample). We apply strict data criteria (both in rejecting observations with an
insufficient number of Fourier bins at frequencies below the white noise
cutoff, and sources for which there were too few observations to provide error
estimates), leaving us with five AGN: IRAS 13224-3809, 1H 0707-495, NGC 4051,
ARK 564 and MRK 766. In all five cases, we observe that the power-law index –
which well describes the power spectrum in this frequency range (and is
preferred over a broken power-law according to maximum likelihood fitting) –
is steep at low energies and becomes increasingly flat towards higher
energies. To-date this has only been reported for a very small number of AGN
(with only one such AGN within our sample – NGC 4051, see McHardy et al. 2004;
Vaughan et al. 2011) and across only a small number of energy bins (Nandra &
Papadakis 2001; Vaughan et al. 2003a; McHardy et al. 2005; Jin et al. 2020);
our analysis therefore represents a much broader investigation into the
energy-dependence of the shape of AGN power spectra.

Adopting the framework of Arévalo & Uttley (2006), we model the power-spectra
of our AGN sample in five energy bands, assuming energy-dependent changes are
driven by propagation only, and extract the emissivity index as a function of
energy, parameterising it as a linear function. This modelling is remarkably
successful in most cases, with the index starting at $\gamma<3$ at low
energies (which may well be consistent with an ADAF providing some of the soft
excess emission, e.g. Done et al. 2012) and rising to $\gamma>4$ above 4-6
keV.

The best-fitting emissivity function varies between sources; we attempt to
correlate $\gamma(E)$ with some of the key macro-properties of the AGN (mass,
spin, inclination and Eddington ratio), finding either a lack of a correlation
or only a tentative hint (with Eddington ratio). In future we will increase
the number of sources for correlation studies and extend the energy range by
utilising data from NuSTAR. It will also be extremely interesting to compare
these results for a simplified disc-like geometry to a more complex geometry
incorporating a more realistic Compton scattering scenario, and to power
spectra of XRBs in a variety of spectral states.

## Acknowledgements

This research was possible through an STFC studentship. We thank Phil Uttley
and Peter Boorman for helpful discussion, as well as the anonymous referee for
their constructive comments. We also reference the use of python libraries
which were essential in this analysis: Numpy (van der Walt et al. 2011), Scipy
(Virtanen et al. 2020), pandas (Wes McKinney 2010), Matplotlib (Hunter 2007)
and emcee (Foreman-Mackey et al. 2013), in addition to XSPEC (Arnaud 1996).

## Data Availability

The data underlying this article are freely accessible in the public HEASARC
archives (https://heasarc.gsfc.nasa.gov/).

## References

  * Ai et al. (2013) Ai Y. L., Yuan W., Zhou H., Wang T. G., Dong X. B., Wang J. G., Lu H. L., 2013, AJ, 145, 90
  * Alston et al. (2013) Alston W. N., Vaughan S., Uttley P., 2013, MNRAS, 429, 75
  * Alston et al. (2014) Alston W. N., Markevičiūtė J., Kara E., Fabian A. C., Middleton M., 2014, MNRAS, 445, L16
  * Alston et al. (2015) Alston W. N., Parker M. L., Markevičiūtė J., Fabian A. C., Middleton M., Lohfink A., Kara E., Pinto C., 2015, MNRAS, 449, 467
  * Alston et al. (2019) Alston W. N., et al., 2019, MNRAS, 482, 2088
  * Alston et al. (2020) Alston W. N., et al., 2020, Nature Astronomy, p. 2
  * Arévalo & Uttley (2006) Arévalo P., Uttley P., 2006, MNRAS, 367, 801
  * Arnaud (1996) Arnaud K. A., 1996, in Jacoby G. H., Barnes J., eds, Astronomical Society of the Pacific Conference Series Vol. 101, Astronomical Data Analysis Software and Systems V. p. 17
  * Ashton & Middleton (2021) Ashton D. I., Middleton M. J., 2021, MNRAS, 501, 5478
  * Buisson et al. (2018a) Buisson D. J. K., et al., 2018a, MNRAS, 475, 2306
  * Buisson et al. (2018b) Buisson D. J. K., et al., 2018b, MNRAS, 480, 3689
  * Casebeer et al. (2006) Casebeer D. A., Leighly K. M., Baron E., 2006, ApJ, 637, 157
  * Chainakun (2019) Chainakun P., 2019, ApJ, 878, 20
  * Churazov et al. (2001) Churazov E., Gilfanov M., Revnivtsev M., 2001, MNRAS, 321, 759
  * Crummy et al. (2006) Crummy J., Fabian A. C., Gallo L., Ross R. R., 2006, MNRAS, 365, 1067
  * Dauser et al. (2012) Dauser T., et al., 2012, MNRAS, 422, 1914
  * Denney et al. (2010) Denney K. D., et al., 2010, ApJ, 721, 715
  * Done & Jin (2016) Done C., Jin C., 2016, MNRAS, 460, 1716
  * Done & Nayakshin (2007) Done C., Nayakshin S., 2007, MNRAS, 377, L59
  * Done et al. (2012) Done C., Davis S. W., Jin C., Blaes O., Ward M., 2012, MNRAS, 420, 1848
  * Done et al. (2013) Done C., Jin C., Middleton M., Ward M., 2013, MNRAS, 434, 1955
  * Emmanoulopoulos et al. (2016) Emmanoulopoulos D., Papadakis I. E., Epitropakis A., Pecháček T., Dovčiak M., McHardy I. M., 2016, MNRAS, 461, 1642
  * Fabian et al. (2004) Fabian A. C., Miniutti G., Gallo L., Boller T., Tanaka Y., Vaughan S., Ross R. R., 2004, MNRAS, 353, 1071
  * Fabian et al. (2009) Fabian A. C., et al., 2009, Nature, 459, 540
  * Fabian et al. (2013) Fabian A. C., et al., 2013, MNRAS, 429, 2917
  * Foreman-Mackey et al. (2013) Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, PASP, 125, 306
  * Gaskell (2004) Gaskell C. M., 2004, ApJ, 612, L21
  * George & Fabian (1991) George I. M., Fabian A. C., 1991, MNRAS, 249, 352
  * Gierliński et al. (2008) Gierliński M., Middleton M., Ward M., Done C., 2008, Nature, 455, 369
  * González-Martín & Vaughan (2012) González-Martín O., Vaughan S., 2012, A&A, 544, A80
  * Gonzalez et al. (2017) Gonzalez A. G., Wilkins D. R., Gallo L. C., 2017, MNRAS, 472, 1932
  * Hunter (2007) Hunter J. D., 2007, Computing in Science & Engineering, 9, 90
  * Ingram & Done (2011) Ingram A., Done C., 2011, MNRAS, 415, 2323
  * Jin et al. (2012) Jin C., Ward M., Done C., 2012, MNRAS, 425, 907
  * Jin et al. (2020) Jin C., Done C., Ward M., 2020, arXiv e-prints, p. arXiv:2005.05857
  * Kara et al. (2016) Kara E., Alston W. N., Fabian A. C., Cackett E. M., Uttley P., Reynolds C. S., Zoghbi A., 2016, MNRAS, 462, 511
  * Kara et al. (2017) Kara E., García J. A., Lohfink A., Fabian A. C., Reynolds C. S., Tombesi F., Wilkins D. R., 2017, Monthly Notices of the Royal Astronomical Society, 468, 3489
  * Kara et al. (2019) Kara E., et al., 2019, Nature, 565, 198
  * Kato et al. (1998) Kato S., Fukue J., Mineshige S., 1998, in Black-hole accretion disks. Edited by Shoji Kato.
  * Kotov et al. (2001) Kotov O., Churazov E., Gilfanov M., 2001, MNRAS, 327, 799
  * Leighly (1999) Leighly K. M., 1999, ApJS, 125, 297
  * Liddle (2007) Liddle A. R., 2007, Monthly Notices of the Royal Astronomical Society: Letters, 377, L74
  * Lyubarskii (1997) Lyubarskii Y. E., 1997, MNRAS, 292, 679
  * McHardy et al. (2004) McHardy I. M., Papadakis I. E., Uttley P., Page M. J., Mason K. O., 2004, MNRAS, 348, 783
  * McHardy et al. (2005) McHardy I. M., Gunn K. F., Uttley P., Goad M. R., 2005, MNRAS, 359, 1469
  * McHardy et al. (2006) McHardy I. M., Koerding E., Knigge C., Uttley P., Fender R. P., 2006, Nature, 444, 730
  * Middleton & Done (2010) Middleton M., Done C., 2010, MNRAS, 403, 9
  * Middleton et al. (2007) Middleton M., Done C., Gierliński M., 2007, MNRAS, 381, 1426
  * Middleton et al. (2009) Middleton M., Done C., Ward M., Gierliński M., Schurch N., 2009, MNRAS, 394, 250
  * Middleton et al. (2011) Middleton M., Uttley P., Done C., 2011, MNRAS, 417, 250
  * Middleton et al. (2016) Middleton M. J., Parker M. L., Reynolds C. S., Fabian A. C., Lohfink A. M., 2016, MNRAS, 457, 1568
  * Mosquera et al. (2013) Mosquera A. M., Kochanek C. S., Chen B., Dai X., Blackburne J. A., Chartas G., 2013, ApJ, 769, 53
  * Nandra & Papadakis (2001) Nandra K., Papadakis I. E., 2001, ApJ, 554, 710
  * Ogle et al. (2004) Ogle P. M., Mason K. O., Page M. J., Salvi N. J., Cordova F. A., McHardy I. M., Priedhorsky W. C., 2004, ApJ, 606, 151
  * Papadakis (2004) Papadakis I. E., 2004, MNRAS, 348, 207
  * Papadakis et al. (2002) Papadakis I. E., Brinkmann W., Negoro H., Gliozzi M., 2002, A&A, 382, L1
  * Papadakis et al. (2016) Papadakis I., Pecháček T., Dovčiak M., Epitropakis A., Emmanoulopoulos D., Karas V., 2016, A&A, 588, A13
  * Pinto et al. (2018) Pinto C., et al., 2018, MNRAS, 476, 1021
  * Puchnarewicz et al. (2001) Puchnarewicz E. M., Mason K. O., Siemiginowska A., Fruscione A., Comastri A., Fiore F., Cagnoni I., 2001, ApJ, 550, 644
  * Reynolds (2014) Reynolds C. S., 2014, Space Sci. Rev., 183, 277
  * Shakura & Sunyaev (1973) Shakura N. I., Sunyaev R. A., 1973, A&A, 24, 337
  * Smith & Vaughan (2007) Smith R., Vaughan S., 2007, MNRAS, 375, 1479
  * Tombesi et al. (2010) Tombesi F., Cappi M., Reeves J. N., Palumbo G. G. C., Yaqoob T., Braito V., Dadina M., 2010, A&A, 521, A57
  * Tombesi et al. (2011) Tombesi F., Cappi M., Reeves J. N., Palumbo G. G. C., Braito V., Dadina M., 2011, ApJ, 742, 44
  * Turner & Reynolds (2021) Turner S. G. D., Reynolds C. S., 2021, MNRAS,
  * Uttley & McHardy (2001) Uttley P., McHardy I. M., 2001, MNRAS, 323, L26
  * Uttley et al. (2014) Uttley P., Cackett E. M., Fabian A. C., Kara E., Wilkins D. R., 2014, A&ARv, 22, 72
  * Vasudevan & Fabian (2007) Vasudevan R. V., Fabian A. C., 2007, MNRAS, 381, 1235
  * Vasudevan et al. (2010) Vasudevan R. V., Fabian A. C., Gandhi P., Winter L. M., Mushotzky R. F., 2010, MNRAS, 402, 1081
  * Vaughan (2005) Vaughan S., 2005, A&A, 431, 391
  * Vaughan et al. (2003a) Vaughan S., Fabian A. C., Nandra K., 2003a, MNRAS, 339, 1237
  * Vaughan et al. (2003b) Vaughan S., Edelson R., Warwick R. S., Uttley P., 2003b, MNRAS, 345, 1271
  * Vaughan et al. (2011) Vaughan S., Uttley P., Pounds K., Nandra K., Strohmayer T. E., 2011, MNRAS, 413, 2489
  * Veron-Cetty & Veron (2010) Veron-Cetty M. P., Veron P., 2010, VizieR Online Data Catalog, p. VII/258
  * Virtanen et al. (2020) Virtanen P., et al., 2020, Nature Methods,
  * Walton et al. (2013) Walton D. J., et al., 2013, ApJ, 777, L23
  * Wang, T. & Lu, Y. (2001) Wang, T. Lu, Y. 2001, A&A, 377, 52
  * Watarai et al. (2000) Watarai K.-y., Fukue J., Takeuchi M., Mineshige S., 2000, PASJ, 52, 133
  * Wes McKinney (2010) Wes McKinney 2010, in Stéfan van der Walt Jarrod Millman eds, Proceedings of the 9th Python in Science Conference. pp 56 – 61, doi:10.25080/Majora-92bf1922-00a
  * Wilkins & Fabian (2011) Wilkins D. R., Fabian A. C., 2011, MNRAS, 414, 1269
  * Wilkins & Fabian (2013) Wilkins D. R., Fabian A. C., 2013, Monthly Notices of the Royal Astronomical Society, 430, 247
  * Wilkins & Gallo (2015) Wilkins D. R., Gallo L. C., 2015, MNRAS, 449, 129
  * Woo & Urry (2002) Woo J.-H., Urry C. M., 2002, ApJ, 579, 530
  * Zhang & Wang (2006) Zhang E.-P., Wang J.-M., 2006, ApJ, 653, 137
  * Zoghbi et al. (2012) Zoghbi A., Fabian A. C., Reynolds C. S., Cackett E. M., 2012, MNRAS, 422, 129
  * van der Walt et al. (2011) van der Walt S., Colbert S. C., Varoquaux G., 2011, Computing in Science Engineering, 13, 22