# Heart rate estimation in intense exercise videos

###### Abstract

Estimating heart rate from video allows non-contact health monitoring with
applications in patient care, human interaction, and sports. Existing work can
robustly measure heart rate under some degree of motion by face tracking.
However, this is not always possible in unconstrained settings, as the face
might be occluded or even outside the camera. Here, we present IntensePhysio:
a challenging video heart rate estimation dataset with realistic face
occlusions, severe subject motion, and ample heart rate variation. To ensure
heart rate variation in a realistic setting we record each subject for around
1-2 hours. The subject is exercising (at a moderate to high intensity) on a
cycling ergometer with an attached video camera and is given no instructions
regarding positioning or movement. We have 11 subjects, and approximately 20
total hours of video. We show that the existing remote photo-plethysmography
methods have difficulty in estimating heart rate in this setting. In addition,
we present IBIS-CNN, a new baseline using spatio-temporal superpixels, which
improves on existing models by eliminating the need for a visible face/face
tracking. We will make the code and data publically available
soon.111https://github.com/ynapolean/IBIS-CNN

Index Terms—  Heart rate estimation, challenging new dataset, spatio-temporal
superpixels.

##  1 Introduction

Remote photo-plethysmography (rPPG) [1], allows non-contact heart rate
estimation. This facilitates applications where contact sensors are difficult
such as infant health monitoring, or athlete monitoring in large scale events
like marathons, where access to individual sensors is unavailable. Initial
methods of rPPG estimation [2, 3] needed the subject to sit motionless in
front of a camera. Other approaches [4, 5, 6, 7] typically used face detectors
and Fourier analysis based techniques. Some methods additionally use skin
tracking/segmentation to estimate heart rate [8, 9]. Other approaches cast
heart rate estimation from video as a blind source separation problem [10,
11]. Most existing rPPG methods including recent deep learning based
approaches [12, 13] typically operate within a ‘constrained’ setting: minor
pose and camera angle variations, subjects are cooperative, are not occluded
and do not perform high-speed movements. Enforcing such constraints is not
possible in real-life situations. For instance in a sports scenario, it is
hard to control lighting conditions and rapid pose changes.

Fig. 1: Methods tracking skin pixels and their color change for heart rate
estimation can perform poorly when the subject’s face is not visible. Here we
introduce the IntensePhysio dataset, which includes high heart rate
variability and is recorded in a relatively ‘unconstrained’ setting, making it
challenging for existing methods. We propose the IBIS-CNN method which relies
on spatio-temporal grouping of pixels instead of face tracking, which improves
on baseline models.

A recent analysis [14] shows that large, high speed motion, tracking based
methods can fail. There are methods that have been proposed to handle partial
facial occlusion [15][16], but these methods rely on extracting information
from other parts of the face that are visible and do not address facial
occlusion (as depicted in Fig. 1). To develop a method for heart rate
estimation from video in a practical setting (like sports, for example the
last image from the left in the second row of Fig. 2), this must be addressed.

Early datasets, MAHNOB [17] and COHFACE [18] feature almost no subject motion.
The PURE dataset features minimal controlled head motion [19]. An important
step towards more realistic scenarios is the seminal ECG-Fitness dataset [20]
where subjects perform fitness-related motions under varying lighting
conditions, demonstrating that previous visual heart rate estimation methods
performed poorly in this scenario. This ECG-Fitness dataset, however, assumes
that faces are still strictly visible, allowing for solutions limited to face
tracking.

In this work, we introduce IntensePhysio: a new, publicly available,
challenging dataset with ample variance in the heart rate, fast motions,
severe appearance changes. It presents a completely unconstrained setting
wherein subjects are allowed exacerbated motion, causing their face to be
sometimes occluded, or they can go out of frame and exhibit greatly varying
face angles. In addition to the new dataset, we analyze existing rPPG
approaches and provide an intentionally simple baseline method called IBIS-CNN
that replaces face tracking with spatio-temporal superpixels. Bobbia[21]
introduced an algorithm called Iterative Boundaries implicit Identification
for superpixel Segmentation (IBIS) [21] for real-time computation of temporal
superpixels for rPPG. We propose IBIS-CNN as a novel application of IBIS for
heart rate estimation in a deep learning setting. IBIS-CNN performs on-par to
existing methods on existing datasets, but significantly outperforms existing
work on more realistic scenarios as exemplified in our new IntensePhysio
dataset. The contributions of this paper are: (i) The new, public,
IntensePhysio dataset with facial occlusions, and high speed motion (ii) A
simple baseline method using superpixels which does not rely on the face to
extract heart rate. (iii) We demonstrate the difficulty of the dataset for
existing rPPG work and (iv) show that our simple baseline performs equally
well on existing datasets, yet significantly better on our more realistic
dataset.

##  2 Method

###  2.1 The IntensePhysio dataset

Fig. 2: Samples from our IntensePhysio rPPG dataset. Note the major pose
changes, specular reflections on the background and on the subject’s skin and
instances where the face region is occluded or not fully in frame.

We collected video recordings and heart rate as the ground truth. The videos
are recorded with a Go-Pro Hero 7 black camera at a resolution of $1,920\times
1,080$ pixels at $\sim$ 60 frames per second. The ground truth heart rate (in
beats per minute) is measured every second, using the Polar Vantage-M heart
rate sensor (chest strap).

For illumination we use LED lights (47W), natural light is also allowed to
enter through a small window. Recordings occur during different times of the
day, resulting in variation of lighting conditions. There are a total of 11
subjects (3 female and 8 male) and all participants are working out on an
ergometer and are given no instructions except to do a workout. This dataset
features the largest heart rate range and variation: from 51 bpm up to 186 bpm
with a standard deviation of 25 bpm. IntensePhysio dataset statistics are
presented in Table 1. Subjects exercised at a moderate to high intensity. They
did not receive any camera-related instructions, such as to face the camera or
be in the frame. This means that the subject’s motion is ’unconstrained’ and
consequently the subject’s face could be occluded or might not be visible at
all, as illustrated in Fig. 2.

Dataset Statistics
---
Mean heart rate (ground truth) | 129
Std. deviation heart rate (ground truth) | 25
Max. heart rate (ground truth) | 186
Min. heart rate (ground truth) | 51
No. of subjects | 11
No. of videos | 15
Avg. video runtime | 1hr. 12min.
Frames per second | 59.94
Table 1: Statistics of IntensePhysio. The dataset features a wide range of
heart rates, from resting to intense workout. Heart rate values are in beats
per minute.

###  2.2 Simple rPPG baseline: IBIS-CNN

To accompany the IntensePhysio dataset, we implement a simple baseline method
for rPPG estimation using temporal superpixels and a convolutional neural
network. The approach is shown in Fig. 3.

Fig. 3: An overview of the IBIS-CNN model. A clip is split into $M$ windows of
10 seconds each. For each frame, $K$ superpixels are extracted and per
superpixel, the average YUV color is stacked over time, where $T$ is the
number of frames and generate a spatio-temporal map. The predicted heart rate
per map from the CNN are averaged to estimate the HR per video.

The model is trained to predict the scalar valued instantaneous heart rate in
beats per minute (bpm) for an input video. We break an input video sequence
into non-overlapping clips of one minute each. For each clip we generate $K$
temporal superpixels using IBIS [21] where $K$ refers to a user-defined value
for the number of superpixels, which we set to 300.

The input images are converted to the $CIE\;l^{*}a^{*}b^{*}$ color space. Then
pixels are grouped iteratively based on their closest ‘seeds’, which are the
initial values for superpixel centers. Pixel grouping is done according to
both the chromatic similarity $D_{lab}$ and spatial proximity $D_{spatial}$ of
the pixel with the associated seed. The metric used to quantify this distance,
$D_{total}$ is given as

$\displaystyle D_{total}$
$\displaystyle=D_{lab}+\theta*D_{spatial},\;\mbox{where}$ $\displaystyle
D_{lab}$
$\displaystyle=\left\|(l,a,b)_{i^{th}pixel}-(l,a,b)_{k^{th}seed}\right\|$
$\displaystyle D_{spatial}$
$\displaystyle=\left\|(x,y)_{i^{th}pixel}-(x,y)_{k^{th}seed}\right\|$

for the $i^{th}$ pixel and the $k^{th}$ seed. $\left\|.\right\|$ is a
Euclidean distance and $\theta$ is defined by $\theta=1/c^{2}$ where $c$ is a
user-defined compacity parameter. The seeds value is propagated temporally, to
generate an output of average RGB color value per temporal superpixel per
frame.

Let $C_{p}(t)$ denote the 3-dimensional per-channel average YUV color signal
of a superpixel $p$ in frame $t$. The 3 average YUV values and grouping all
$K$ superpixels gives a 2d matrix of size $K\times 3$. Stacking these over a
temporal window of $T$ frames, we obtain a $3\times K\times T$ tensor, which
can be seen as a $K\times T$ color image which we call a spatial-temporal map.
This map is 10 seconds long, while each clip lasts around 1 minute. We slide a
temporal window over the clip to obtain M maps per clip using a window of 10
seconds (stride = 0.5s).

##  3 Experiments

We compare against two state of the art rPPG methods in HR-CNN [20] and
RhythmNet [12]. Both these methods involve pre-processing the input video
frames using face-detection (with an additional alignment step in RhythmNet),
cropping and re-sizing using specific toolboxes. We use train-test splits as
specified by the authors. For our dataset, we make use of a train-test split
featuring 8 subjects in the train set and 3 in the test set. The number of
superpixels is $K=300$.

Fig. 4: Predictions of the IBIS-CNN model compared with he ground truth and
HR-CNN on subjects from the test set. The HR-CNN model outputs a constant
value, while our predictions exhibit learning behavior, following the ground
truth.

| Model | PURE | ECG-Fitness | IntensePhysio
---|---|---|---|---
RMSE | HR-CNN | 11.00 | 19.15 | 25.27
RhythmNet | 19.67 | 20.47 | 32.67
IBIS-CNN (ours) | 11.99 | 17.03 | 22.01
MAE | HR-CNN | 8.72 | 14.48 | 22.19
RhythmNet | 17.46 | 16.82 | 28.36
IBIS-CNN (ours) | 9.39 | 13.75 | 16.53

Table 2: RMSE and MAE results for the proposed IBIS-CNN and two baseline
methods in HR-CNN and RhythmNet on the test sets of PURE, ECG-fitness and
IntensePhysio (ours).

Due to lack of publicly available code for RhythmNet [12], we re-implemented
it as per the authors description in their paper. The HR-CNN [20] model and
code is publicly available. We evaluate with mean absolute error (MAE) and
root mean square error (RMSE) calculated over the number of spatio-temporal
maps. Results in Table 2 show that our method achieves similar performance to
others on PURE and ECG-Fitness, validating our baseline. Yet, our baseline
significantly outperforms others on our new IntensePhysio dataset.

All models show poor performance on our dataset. To establish the challenging
nature of the IntensePhysio dataset for methods that rely on face tracking, we
run an off the shelf face detector (Dlib [22]) and find that _on average we
are able to detect faces in $\sim$28$\%$ of frames per video in our dataset.On
the other hand we observed a successful face detection in $\sim$74$\%$ frames
per video in ECG-Fitness_. Since the face region (particularly the cheeks and
forehead) contain the most information on heart rate [13][20], this leads to a
poor performance when models rely on face detection and tracking.

The predictions of our IBIS-CNN model and HR-CNN along with the ground truth
are shown in Fig. 4. A few sample frames from the video are also presented
with the dotted lines connecting them to the minute of extraction. In these
frames, there is a lot of specular reflection on his skin (due to sweating),
the subject’s face is not visible. This happens particularly around the
$30^{th}$ minute, where the error is larger. The HR-CNN model outputs a
constant value, indicating that it might not be learning relevant rPPG
features. The IBIS-CNN model is able to predict heart rate reasonably well
despite occlusions, specular reflections etc. However, we also note that the
predictions from IBIS-CNN are not always correlated with the ground truth,
especially in the second plot of 4. This could be because of rapid motion and
would be interesting future direction of research. IBIS-CNN does not solve all
the challenges posed by the IntensePhysio dataset but is an alternative
baseline with a new approach.

Comparison of input representations. We compare two existing methods for
temporal superpixel generation - IBIS [21] and TS-PPM [23]. We find that the
performance of TS-PPM is lower than that of IBIS-CNN (as seen in Table 3
possibly because, in the IBIS method pixel membership is constrained to not
vary beyond a certain threshold to maintain temporally coherent RGB traces.

As a pre-processing step, the IBIS superpixel generation is more
computationally efficient than the TS-PPM. To generate the superpixel results,
IBIS processes $5.33$ frames/sec on average whereas TS-PPM was averages $0.56$
frames/sec. Thus, the IBIS method is more suited to the rPPG estimation task.
These results were obtained using PURE dataset.

| Model | Validation error for PURE
---|---|---
RMSE | TS-PPM [23] | 13.97
IBIS-CNN (ours) | 11.99
MAE | TS-PPM [23] | 10.81
IBIS-CNN (ours) | 9.39

Table 3: Comparison between temporal superpixel methods for heart rate
estimation.

##  4 Conclusions

We present IntensePhysio, a challenging new dataset for heart rate estimation.
The dataset features large subject motion with frequent face occlusions and
cases of facial region absent from the frame entirely. Through a comparative
study, we observe a considerable degradation in the performance of the
existing state of the art methods on this new dataset, especially methods
relying on face detection and tracking. This highlights IntensePhysio as a
challenging dataset for heart rate estimation, indicating that occlusion and
non-visible facial regions are key factors for this performance degradation.
Hence, we propose IBIS-CNN as a new baseline method for heart rate estimation
(using temporal superpixels) which significantly outperforms state of the art
methods on our challenging new dataset in addition to the existing ones.
However, the IBIS-CNN baseline predictions are not strongly correlated with
the actual ground truth always.This shows that there are further problems
posed by our dataset that require addressing. Also, there is a need for a
thorough hyperparameter tuning while generating these temporal superpixels
using IBIS, specific to our task. So, it would be of significance to
investigate the possibility of developing IBIS-CNN as an end to end learnable
pipeline.

## References

  * [1] Qi Zhan, Wenjin Wang, and Gerard de Haan,  “Analysis of cnn-based remote-ppg to understand limitations and sensitivities,”  Biomedical optics express, vol. 11, no. 3, pp. 1268–1283, 2020\.
  * [2] Sungjun Kwon, Hyunseok Kim, and Kwang Suk Park,  “Validation of heart rate extraction using video imaging on a built-in camera system of a smartphone,”  in 2012 Conference of the IEEE Engineering in Medicine and Biology Society. IEEE, 2012, pp. 2174–2177.
  * [3] Wim Verkruysse, Lars O Svaasand, and J Stuart Nelson,  “Remote plethysmographic imaging using ambient light.,”  Optics express, vol. 16, no. 26, pp. 21434–21445, 2008.
  * [4] E Jonathan and Martin Leahy,  “Investigating a smartphone imaging unit for photoplethysmography,”  Physiological measurement, vol. 31, no. 11, pp. N79, 2010.
  * [5] Amogh Gudi, Marian Bittner, Roelof Lochmans, and Jan Van Gemert,  “Efficient real-time camera based estimation of heart rate and its variability,”  in ICCV Workshops, 2019.
  * [6] Amogh Gudi, Marian Bittner, and Jan van Gemert,  “Real-time webcam heart-rate and variability estimation with clean ground truth for evaluation,”  Applied Sciences, vol. 10, no. 23, 2020.
  * [7] Huishi Zhu, Yuejin Zhao, and Liquan Dong,  “Non-contact detection of cardiac rate based on visible light imaging device,”  in Optics and Photonics for Information Processing VI. International Society for Optics and Photonics, 2012, vol. 8498, p. 849806.
  * [8] Chuanxiang Tang, Jiwu Lu, and Jie Liu,  “Non-contact heart rate monitoring by combining convolutional neural network skin detection and remote photoplethysmography via a low-cost camera,”  in CVPR Workshops, 2018.
  * [9] Zitong Yu, Wei Peng, Xiaobai Li, Xiaopeng Hong, and Guoying Zhao,  “Remote heart rate measurement from highly compressed facial videos: an end-to-end deep learning solution with video enhancement,”  in ICCV, 2019.
  * [10] Ming-Zher Poh, Daniel J McDuff, and Rosalind W Picard,  “Non-contact, automated cardiac pulse measurements using video imaging and blind source separation.,”  Optics express, vol. 18, no. 10, pp. 10762–10774, 2010.
  * [11] Magdalena Lewandowska, Jacek Ruminski, Tomasz Kocejko, and Jedrzej Nowak,  “Measuring pulse rate with a webcam—a non-contact method for evaluating cardiac activity,”  in 2011 federated conference on computer science and information systems (FedCSIS). IEEE, 2011, pp. 405–410.
  * [12] Xuesong Niu, Hu Han, Shiguang Shan, and Xilin Chen,  “VIPL-HR: A multi-modal database for pulse estimation from less-constrained face video,”  in ACCV, 2018.
  * [13] Weixuan Chen and Daniel McDuff,  “Deepphys: Video-based physiological measurement using convolutional attention networks,”  in ECCV, 2018.
  * [14] Changchen Zhao, Peiyi Mei, Shoushuai Xu, Yongqiang Li, and Yuanjing Feng,  “Performance evaluation of visual object detection and tracking algorithms used in remote photoplethysmography,”  in ICCV Workshops, 2019.
  * [15] Yuichiro Maki, Yusuke Monno, Masayuki Tanaka, and Masatoshi Okutomi,  “Remote heart rate estimation based on 3d facial landmarks,”  in EMBC, 2020.
  * [16] Kun Zheng, Kangyi Ci, Jinling Cui, Jiangping Kong, and Jing Zhou,  “Non-contact heart rate detection when face information is missing during online learning,”  Sensors, vol. 20, no. 24, 2020.
  * [17] Mohammad Soleymani, Jeroen Lichtenauer, Thierry Pun, and Maja Pantic,  “A multimodal database for affect recognition and implicit tagging,”  IEEE transactions on affective computing, vol. 3, no. 1, pp. 42–55, 2011.
  * [18] Guillaume Heusch, André Anjos, and Sébastien Marcel,  “A reproducible study on remote heart rate measurement,”  arXiv preprint arXiv:1709.00962, 2017.
  * [19] Ronny Stricker, Steffen Müller, and Horst-Michael Gross,  “Non-contact video-based pulse rate measurement on a mobile service robot,”  in Symposium on Robot and Human Interactive Communication, 2014\.
  * [20] Radim Špetlík, Vojtech Franc, and Jirí Matas,  “Visual heart rate estimation with convolutional neural network,”  in BMVC, 2018.
  * [21] Serge Bobbia, Duncan Luguern, Yannick Benezeth, Keisuke Nakamura, Randy Gomez, and Julien Dubois,  “Real-time temporal superpixels for unsupervised remote photoplethysmography,”  in CVPR Workshop, 2018.
  * [22] Davis E. King,  “Dlib-ml: A machine learning toolkit,”  Journal of Machine Learning Research, vol. 10, pp. 1755–1758, 2009\.
  * [23] Se-Ho Lee, Won-Dong Jang, and Chang-Su Kim,  “Temporal superpixels based on proximity-weighted patch matching,”  in ICCV, 2017.


# Joint Estimation of System Inertia and Load Relief

Julius Susanto Senior Member, IEEE, Alireza Fereidouni Member, IEEE and Dean
Sharafi Senior Member, IEEE System Management, Australian Energy Market
Operator, Perth, WA 6000, Australia

###### Abstract

Modern power systems with high share of renewable generation are at the risk
of rapid changes in frequency resulting from contingencies. The importance of
an accurate assessment of system inertia and load relief, as well as frequency
changes in the critical timescales is more pronounced in these power systems.
This knowledge serves as an insight which will guide the solutions required
for enhanced security and resilience. A novel procedure for simultaneously
assessing system inertia and load relief using a system frequency response
(SFR) model to fit measured disturbance data is presented. Results of applying
the proposed approach on generator contingencies in the South West
Interconnected System (SWIS) in Western Australia demonstrate the validity of
the method and indicate that it can overcome some of the limitations observed
in conventional inertia estimation methods based on the linearised swing
equation, such as the sliding window and polynomial fit methods.

###### Index Terms:

System inertia estimation, load relief, load damping

##  I Introduction

Many power systems across the world are experiencing a gradual decline in
synchronous inertia levels as synchronous generators are increasingly being
displaced by generation interfaced to the grid via power-electronic converters
[1]. As inertia levels decrease, the average frequency in the power system is
subject to more severe and rapid fluctuations in response to active power
disturbances, e.g. from generation or load contingencies. It is therefore
becoming more crucial that estimates of system inertia are reliable and
accurate.

Most system inertia estimation techniques, whether offline estimates such as
[2], [3], [4] and [5] or online estimates such as [6], [7] and [8], rely on a
linearised representation of the swing equation evaluated at the onset of a
disturbance ($t=0$):

$\overline{KE}=-\frac{f_{n}}{2}\times\frac{P_{cont}}{\frac{\overline{d\Delta
f}}{dt}\bigg{|}_{t=0}}$ (1)

where $\overline{KE}$ is the system inertia estimate (MW.s), $f_{n}$ is the
nominal frequency (e.g. $50$ Hz), $P_{cont}$ is the disturbance size (MW) and
$\frac{d\overline{\Delta f}}{dt}\big{|}_{t=0}$ is an estimate of the rate of
change of frequency (RoCoF) evaluated at the onset of the disturbance (Hz/s).

For this linearised representation to be valid, it is assumed that i) there
are no frequency deviations at the onset of the disturbance, i.e. $\Delta
f(t=0)=0$, ii) primary frequency response has yet to begin acting, and ii)
load relief is negligible in the early stages after a disturbance. To satisfy
these assumptions, inertia estimates of the form described in (1) are often
performed using historical generator contingency events, filtered to only
select for single instantaneous disturbances (e.g. protection trip) as opposed
to ramp downs (e.g. due to a generating unit boiler malfunction / failure) or
contingencies involving multiple generating units. With such sudden single
unit trips, the contingency size $P_{cont}$ can normally be estimated quite
accurately [5]. Alternatively, known active power perturbations by normal
operational events such as the reversal of power in an HVDC interconnector [9]
or an intentionally injected signal [10] can also be used.

However, the RoCoF is arguably the most sensitive parameter in (1). Ideally,
the RoCoF would be calculated based on the sample-by-sample time derivative of
a high-resolution frequency trace at the onset of the disturbance. However,
high-resolution frequency measurements from Phasor Measurement Units (PMUs) or
fault recorders contain white noise introduced by transducer, quantisation and
signal processing errors [11]. Moreover, as depicted in Fig. 1, the frequency
trace can contain transient [5] and/or oscillatory components [2]. Therefore,
calculating the time derivative using consecutive frequency samples would
likely lead to significant errors [3].

Figure 1: Transients and oscillations in frequency measurements

To mitigate these issues, moving average windows or low-pass smoothing filters
have been proposed in the literature, as well as in practical applications.
For example, a sliding window has been adopted as a standard in Europe (e.g.
[3] and [5]), where the RoCoF estimate is calculated by applying a 500-ms
sliding window over the frequency trace and then finding the maximum value of
the sample-by-sample time derivative (since the maximum RoCoF theoretically
occurs at the onset of the contingency). Alternative techniques include
fitting a polynomial function to the raw frequency measurements [2] or simply
applying a 0.5 Hz low-pass filter to the frequency trace [4]. The common
thread is that the raw frequency measurements are smoothed out in some way to
prevent spurious RoCoF measurements.

It was observed in [12] that system inertia estimates based on this approach
are only accurate for events with relatively high RoCoFs, in which the
assumption that the initial frequency decline is retarded by system inertia
alone remains valid. For smaller disturbances or high inertia conditions where
the RoCoF is low, the influence of primary frequency response becomes more
prominent and may invalidate this assumption. Furthermore, the choice of
sliding window length (or cut-off frequency for the smoothing filter) will
also affect the quality of the estimate, with higher levels of smoothing
leading to larger estimation errors. It was also shown that inertia can be
materially overestimated when wind farms provide virtual inertia emulation or
fast frequency response [13].

Another factor that may confound system inertia estimates is the presence of
load relief (or damping), which is the natural sensitivity of loads in a
network to changes in system frequency. As frequency declines after a
generation contingency, so does power consumption from frequency-sensitive
loads such as induction motors [14]. Load relief acts to correct the active
power imbalance after a contingency and thus helps to arrest the decline in
frequency. As a result, system inertia estimates based on the linearised swing
equation after a disturbance may mistake the effects of load relief as
additional inertia.

In the literature, the system-level load relief factor (LRF), which is also
referred to as the load damping coefficient or load-frequency characteristic,
is typically estimated separately from system inertia. Earlier studies
attempted to directly measure load relief by conducting live tests such as the
intentional tripping of generators and transmission lines, e.g. in Great
Britain [15], Norway [16] and Ireland [17]. A less invasive approach is to
estimate the system-level LRF after a contingency event using measurements of
total system load from the Energy Management System (EMS). However, EMS
sampling resolution may not be high enough to accurately assess system load
changes after a contingency [18]. To mitigate EMS resolution issues, high-
speed fault recorder data can be used to estimate the system load change after
a contingency [12]. Finally, estimating LRF has also been investigated using a
bottom-up approach, e.g. using high-speed fault recorder data at substation
level [14] or at the level of individual consumer devices [19]. However, it is
remains unclear how to aggregate the individual LRF estimates into a
meaningful system-level LRF.

This paper takes a different approach to system inertia and load relief
estimation. Rather than using the linearised swing equation, the system
inertia and load relief factor are treated as unknown parameters in a low-
order system frequency response (SFR) model and the model is then fitted to
post-disturbance frequency and generator output measurements. This approach
exploits the SFR model’s ability to accurately predict the system frequency
trajectory after a disturbance, and thus use more measurement data as inputs
to the estimation model.

The main contribution of this paper is the introduction of a novel model
fitting approach to inertia estimation that uses a system frequency response
model to simultaneously estimate both system inertia and system-level load
relief. The proposed method is shown to be credible when compared with other
inertia estimation methods, while requiring less discretion in the choice of
parameters, e.g. sliding window length or polynomial order.

The remainder of this paper is organised as follows: In Section II, the system
inertia estimation method based on model-fitting is proposed, with the added
benefit that the method concurrently estimates the system-level load relief
factor. This is followed in Section III by several computer simulation case
studies comparing the performance of the proposed approach with other methods.
In Section IV, case studies are performed using data from a real power system
to illustrate some of the advantages of the proposed method over other
approaches. A discussion of the results in this paper is then presented in
Section V, and finally concluding remarks and avenues for future work are
described in Section VI.

##  II Proposed Methodology

Rather than using the linearised swing equation in (1) to estimate system
inertia, we propose to jointly estimate system inertia and load relief /
damping by fitting a low-order SFR model to measured frequency data.

At a high level, the proposed model-fitting estimation method has two steps as
follows:

  1. 1.

Compile generator active power measurement traces to form an aggregate system
primary frequency response (PFR) trace and contingency trace

  2. 2.

Estimate system inertia and load relief factor by fitting an SFR model to the
measured frequency trace

###  II-A Step 1: Measured system PFR and contingency traces

The measured system PFR trace is formed by aggregating the active power
measurement traces from individual generators and interruptible loads. In
order for the measured system PFR to be accurate (i.e. representative of the
actual system PFR), the following conditions should be met:

  * •

Active power measurements are high resolution and time-synchronised, e.g. from
fault recorders or PMUs with time synchronisation via GPS or IEEE Std 1588
Precision Time Protocol [20].

  * •

There are measurements at each major generating facility and interruptible
load in the system. Where this is not practical, there should be at least
measurements at each of the main PFR providers.

Note that with this approach, the influence of any PFR not explicitly included
in the measured system PFR trace will be considered as load relief. This could
include the influence of behind-the-meter embedded generation or small
distribution-connected generation that may not be visible by fault recorders
or PMUs.

Similarly, the contingency trace is formed by extracting the active power
measurement trace of the contingency (e.g. power flow from a generator or
radial network element).

###  II-B Step 2: Estimation of system inertia and LRF via model fitting

The linear ordinary differential equation (ODE) for a generic low-order SFR
model after a contingency can be formulated from the swing equation as
follows:

$\frac{d\Delta
f(t)}{dt}=\frac{f_{n}}{2KE}\left[p(t)-P_{cont}(t)-DP_{load}\Delta f(t)\right]$
(2)

where $\Delta f(t)$ is the frequency deviation at time $t$ (Hz), $f_{n}$ is
the nominal frequency (e.g. $50$ Hz), $KE$ is the post-contingent system
inertia (in MW.s), $P_{cont}(t)$ is the contingency at time $t$ (MW), $p(t)$
is the aggregate system PFR at time $t$ (MW), $D$ is the load relief factor (%
MW/Hz) and $P_{load}$ is the system load at the onset of the contingency (MW).

Equation (2) can be rewritten in discrete form with a fixed time step $\Delta
t$ as follows:

$\Delta f_{k}=\Delta
f_{k-1}+\frac{f_{n}}{2KE}\left[p_{k-1}-P_{cont,k-1}-DP_{load}\Delta
f_{k-1}\right]\Delta t$ (3)

In (3), the measured system PFR trace $p_{k}$ and contingency trace
$P_{cont,k}$ were compiled in Step 1 and are known values, as is the system
load pre-contingency $P_{load}$. This leaves the system inertia $KE$ and load
relief factor $D$ as the unknown parameters to be estimated.

There are many ways to solve this parameter estimation problem, for example
using non-linear least squares [21], Newton-Raphson algorithms [22] or with
Kalman filters [23]. In this paper, a recursive non-linear least squares
approach similar to the method described in [24] is used to minimise an error
function based on the Root Mean Squared Error (RMSE):

$\min\left\\{\sqrt{\frac{1}{N}\sum_{k=0}^{N}(\Delta
f_{k}-\Delta\hat{f_{k}})^{2}}\right\\}$ (4)

where $f$ is the measured frequency, $\hat{f}$ is the frequency calculated
from the SFR model and $N$ is the total number of samples.

##  III Computer Simulations

Computer simulations were performed on the IEEE 9-bus and New England IEEE
39-bus systems to compare the performance of the proposed model fitting
approach with two commonly used inertia estimation methods:

  1. 1.

Sliding window method [3]

  2. 2.

Inoue polynomial fit method [2]

The time-domain RMS simulations were conducted in DIgSILENT PowerFactory v2020
software and in each case, loads were modelled to be frequency-dependent with
a load relief factor of $D=4\%$.

###  III-A Case 1: IEEE 9-Bus System

In Case 1, the IEEE 9-bus system depicts a small islanded system with three
generators and 315 MW of load. The default IEEE 9-bus system was modified to
include a standard turbine-governor model (TGOV) for each generator. To
simulate a frequency disturbance, generator G3 was tripped from 85 MW.

A summary of the best and worst inertia and load relief estimates for each of
the different methods is shown in Table I. It should be noted that only a
single estimate is shown for the model fitting approach as there are no
adjustable parameters (such as window length or polynomial order).

TABLE I: Case 1: Summary of results Method | KE (MW.s) | D (%) | Remarks
---|---|---|---
Actual | 1,522 | 4.00 % |
Sliding window (best) | 1,524 | - | 60-ms window
Sliding window (worst) | 2,810 | - | 1,000-ms window
Inoue (best) | 1,547 | - | 20th order
Inoue (worst) | 1,730 | - | 14th order
Model fitting | 1,589 | 3.45 % |

System inertia estimates using the sliding window and Inoue polynomial fit
methods with varying window lengths and polynomial order are shown in Fig. 2a
and Fig. 2b respectively. The figures indicate that the accuracy of the
inertia estimate is sensitive to the choice of window length or polynomial
order. Fig. 2c shows the simulated (measured) frequency trace vs the fitted
frequency trace after applying the model fitting method.

(a) Sliding window method (variable window length)

(b) Inoue method (variable polynomial order)

(c) Model fitting method frequency fit

Figure 2: Case 1: IEEE 9-bus system

###  III-B Case 2: New England IEEE 39-Bus System

Relative to Case 1, the New England IEEE 39-bus system depicts a larger
interconnected power system comprising nine generators and 6,097 MW of load.
To simulate a frequency disturbance, generator G10 was tripped from 200 MW.

Like in Case 1, a summary of the best and worst inertia and load relief
estimates from each of the different methods is shown in Table II. The inertia
estimates with varying window length and polynomial order are shown in Fig. 3a
and Fig. 3b respectively, and the measured vs fitted frequency trace after
applying the model fitting method is shown in Fig. 3c.

TABLE II: Case 2: Summary of results Method | KE (MW.s) | D (%) | Remarks
---|---|---|---
Actual | 78,270 | 4.00 % |
Sliding window (best) | 78,161 | - | 510-ms window
Sliding window (worst) | 12,270 | - | 20-ms window
Inoue (best) | 66,443 | - | 12th order
Inoue (worst) | 266,772 | - | 10th order
Model fitting | 75,499 | 3.44 % |

An interesting feature of the simulation results in Case 2 is the presence of
oscillations in the frequency trace as generators in the system swing against
each other (see Fig. 3c). The oscillations cause the sliding window and
polynomial fit methods to overestimate the RoCoF at shorter window lengths (as
observed in [3]) and high polynomial orders (i.e. because of overfitting).

(a) Sliding window method (variable window length)

(b) Inoue method (variable polynomial order)

(c) Model fitting method frequency fit

Figure 3: Case 2: New England IEEE 39-bus system

###  III-C Case 3: IEEE 9-Bus System with Fast Frequency Response

In Case 3, the IEEE 9-bus system in Case 1 is augmented with a 30 MW BESS
providing FFR. This case is included to examine the effects of FFR on inertia
estimates and validate the findings in [13], i.e. material errors in inertia
estimates can arise when the sliding window method is applied to systems with
FFR.

As in the previous cases, a summary of the best and worst inertia and load
relief estimates from each of the different methods is shown in Table III. The
inertia estimates with varying window length and polynomial order are shown in
Fig. 4a and Fig. 4b respectively, and the measured vs fitted frequency trace
after applying the model fitting method is shown in Fig. 4c.

(a) Sliding window method (variable window length)

(b) Inoue method (variable polynomial order)

(c) Model fitting method frequency fit

Figure 4: Case 3: IEEE 9-bus system with FFR

The results indicate a deterioration in estimation accuracy from the sliding
window method from Case 1, particularly with increasing window length. The
Inoue polynomial fit method performs better, with a narrower range of
estimates, but the accuracy is sensitive to the choice of polynomial order and
it may be difficult to ascertain a priori whether the order selected is the
appropriate one.

TABLE III: Case 3: Summary of results Method | KE (MW.s) | D (%) | Remarks
---|---|---|---
Actual | 1,522 | 4.00 % |
Sliding window (best) | 1,785 | - | 30-ms window
Sliding window (worst) | 4,209 | - | 1,000-ms window
Inoue (best) | 1,522 | - | 20th order
Inoue (worst) | 1,812 | - | 13th order
Model fitting | 1,662 | 3.76 % |

##  IV Application to a Real Power System

The following case studies are based on actual generator contingency events
from the South West Interconnected System (SWIS) in Australia. The SWIS is a
medium-sized islanded power system serving approximately 1.2 million customers
in the southwest region of Western Australia. The SWIS covers a vast
geographic area of 261,000 square kilometres (larger than the land area of the
United Kingdom), but with a historical record peak demand of just 4.3 GW.

TABLE IV: Actual contingency case descriptions Case | Description
---|---
Case 4 | Generator prime-mover trip
Case 5 | Network fault causing a generator trip

###  IV-A Case 4: Generator turbine trip

On 17 March 2020 at 18:06pm, an open-cycle gas turbine generator tripped from
161 MW causing frequency to drop to 49.55 Hz (at the nadir). At the onset of
the contingency, the system load was 2,441 MW and the total post-contingent
generator inertia was 15,009 MW.s (calculated by summing the known generator
inertia values for all online machines).

The inertia estimates with varying window length and polynomial order are
shown in Fig. 5a and Fig. 5b respectively, and the measured vs fitted
frequency trace after applying the model fitting method is shown in Fig. 5c.

Unlike in the computer simulations, the actual system inertia is not known.
The post-contingent generator inertia is known, but this does not include load
inertia so is more representative of the floor for the system inertia
estimate. The summary of the results in Table V show the range (minimum and
maximum) of estimates from the sliding window and Inoue methods, as well as
the single point estimate from the model fitting method.

TABLE V: Case 4: Summary of results Method | KE (MW.s) | D (%) | Remarks
---|---|---|---
Generator inertia | 15,009 | - |
Sliding window (min) | 20,735 | - | 170-ms window
Sliding window (max) | 26,833 | - | 20-ms window
Inoue (min) | 18,482 | - | 4th order
Inoue (max) | 56,750 | - | 15th order
Model fitting | 19,610 | 2.85 % |

The sliding window method in Fig. 5a exhibits a similar shape to that in Cases
1 and 3, where beyond a threshold window length (in this case 170-ms), the
inertia estimate continues to rise steadily without settling. In Cases 1 and
3, the actual system inertia was close to the estimate at the threshold window
length, which in this case would be 20,735 MW.s.

At high polynomial orders, the Inoue method in Fig. 5b appears to settle at
around 26,000 MW.s. However, such an estimate would suggest a load inertia of
roughly 42% of system inertia, which is inconsistent with previous inertia
estimates in the SWIS, as well as estimates reported in the literature (such
as [25] and [26]), where load inertia is typically estimated to be 10 - 25% of
system inertia.

The model fitting results in Fig. 5c show good alignment between the simulated
and measured frequency traces ($RMSE=0.0095$) and the inertia estimate of
19,610 MW.s is on the lower end of the sliding window and and Inoue estimate
range. This suggests that the actual system inertia is likely in the 19,000 -
21,000 MW.s range.

(a) Sliding window method (variable window length)

(b) Inoue method (variable polynomial order)

(c) Model fitting method frequency fit

Figure 5: Case 4: Generator turbine trip

###  IV-B Case 5: Network fault causing a generator trip

On 22 October 2020 at 16:25pm, a 330 kV network fault caused the trip of a
combined-cycle gas power plant, resulting in an instantaneous loss of 231 MW
and frequency to fall to 49.46 Hz (at the nadir). At the onset of the
contingency, the system load was 1,991 MW and the total post-contingent
generator inertia was 11,612 MW.s. The summary of the results in Table VI show
the range (minimum and maximum) of estimates from the sliding window and Inoue
methods, as well as the single point estimate from the model fitting method.

TABLE VI: Case 5: Summary of results Method | KE (MW.s) | D (%) | Remarks
---|---|---|---
Generator inertia | 11,612 | - |
Sliding window (min) | 3,609 | - | 40-ms window
Sliding window (max) | 18,392 | - | 1,000-ms window
Inoue (min) | 6,567 | - | 30th order
Inoue (max) | 22,309 | - | 3rd order
Model fitting | 15,534 | 2.01 % |

The network fault causes a transient component in the frequency trace, which
results in overestimates of system inertia with the sliding window method when
the window length is small, as shown Fig. 6a for window lengths $<$200-ms.
This is broadly similar to how the sliding window method performs with the
oscillatory frequency trace observed in Case 2. In that previous case, the
sliding window method is more accurate beyond a ”knee” point window length,
where the effects of the oscillations are filtered out. In this case, the knee
point occurs at a window length of 210-ms and a system inertia estimate of
14,000 MW.s.

Similarly, the Inoue method will tend to over-fit the frequency trace (and
particularly the transient component) at high polynomial order. This can be
seen in Fig. 6b, where the system inertia estimate is below the known
generator inertia at polynomial orders $>$15\. The Inoue method does not
credibly provide any useful information regarding the correct range for system
inertia.

The model fitting results in Fig. 6c show reasonably good alignment between
the simulated and measured frequency traces ($RMSE=0.0146$) and yields an
inertia estimate of 15,534 MW.s. This estimate is broadly consistent with the
sliding window method using a window length near the knee point.

(a) Sliding window method (variable window length)

(b) Inoue method (variable polynomial order)

(c) Model fitting method frequency fit

Figure 6: Case 5: Network fault causing a generator trip

##  V Discussion

The results from the case studies indicate that system inertia estimation
using the proposed model-fitting approach is a credible alternative to
conventional approaches based on the linearised swing equation. Good alignment
between the measured frequency and the simulated frequency implied by the
fitted SFR model suggests that the parameters in the SFR model are able to
endogenously explain the frequency trajectory. The proposed approach also has
the practical advantage of providing a single point estimate without needing
to select parameters such as the sliding window length or polynomial order.

Interestingly, the case studies highlight that system inertia estimates using
the sliding window and Inoue polynomial fit methods are very sensitive to the
choice of these parameters. The computer simulations show that depending on
the scenario and choice of window length or polynomial order, system inertia
estimates can be either very accurate or have large errors, e.g. $\pm$20% to
50%. While there are broad patterns that are observed in the case studies, for
example, systematic under-estimates when using the Inoue method with high
polynomial orders in frequency traces containing transient or oscillatory
components, there are no specific rules or heuristics to guide parameter
selection. As the case studies in Section IV show, it is difficult to select
the most appropriate parameter in practice, particularly as measurement noise
present in frequency traces can also affect the inertia estimate.

The proposed model-fitting approach does not suffer from these issues and in
the computer simulations, the method yields system inertia estimates within 3%
to 9% of the true value.

##  VI Conclusion

This paper proposed a novel approach for simultaneously estimating system
inertia and the load relief factor. The proposed approach is based on fitting
an SFR model to measured contingency / disturbance data. Case studies
demonstrated the applicability and validity of applying the new approach on
computer simulations with the IEEE 9-bus and 39-bus systems, and on real
events in the SWIS in Western Australia. Furthermore, the proposed method is
not sensitive to the selection of algorithm parameters such as the sliding
window length or polynomial order.

An avenue for future work is to investigate heuristics for the appropriate
selection of sliding window length and polynomial order, potentially using the
proposed model-fitting approach to calibrate the range of estimates.

##  Appendix A Removal of inertial components from generator output traces

High-speed active power measurements of synchronous generators contain both
the inertial response and PFR from the generators. Since inertia is endogenous
to models based on the swing equation, the system PFR trace can only be used
if the inertial components are removed first (otherwise the effects of inertia
are double counted).

A simple washout filter is proposed to isolate inertial components from the
PFR using the measured frequency as an input and the known inertia constants
of the generators contributing to the system PFR as a parameter (see Fig. 7).
The inertial components are then removed from the measured system PFR trace.

Figure 7: Estimating inertial components of PFR with a washout filter

The selection of the washout time constant $T_{w}$ is a trade-off between
accuracy and noise. Ideally, a very small time constant is applied so that the
washout filter more closely approximates a true differentiator. However, small
time constants are more sensitive to noise in the frequency measurement (refer
to $T_{w}=25$-ms in Fig. 8). On the other hand, a large time constant will
reduce the accuracy of the estimate despite appearing less noisy (refer to
$T_{w}=200$-ms in Fig. 8). It was empirically found that a washout time
constant of 60 ms to 80 ms yielded an acceptable balance between accuracy and
noise (refer to $T_{w}=60$-ms in Fig. 8).

Figure 8: Inertia removal with different washout time constants

## References

  * [1] B. Hartmann, I. Vokony, and I. Táczi, “Effects of decreasing synchronous inertia on power system dynamics—overview of recent experiences and marketisation of services,” _International Transactions on Electrical Energy Systems_ , vol. 29, no. 12, p. e12128, 2019, e12128 ITEES-18-1219.R1. [Online]. Available: https://onlinelibrary.wiley.com/doi/abs/10.1002/2050-7038.12128
  * [2] T. Inoue, H. Taniguchi, Y. Ikeguchi, and K. Yoshida, “Estimation of power system inertia constant and capacity of spinning-reserve support generators using measured frequency transients,” _IEEE Transactions on Power Systems_ , vol. 12, no. 1, pp. 136–143, 1997.
  * [3] ENTSO-E, “Frequency measurement requirements and usage,” https://docstore.entsoe.eu/documents/soc%20documents/regional_groups_continental_europe/2018/tf_freq_meas_v7.pdf, European Network of Transmission System Operators for Electricity (ENTSO-E), Tech. Rep., 2018.
  * [4] D. P. Chassin, Z. Huang, M. K. Donnelly, C. Hassler, E. Ramirez, and C. Ray, “Estimation of wecc system inertia using observed frequency transients,” _IEEE Transactions on Power Systems_ , vol. 20, no. 2, pp. 1190–1192, May 2005.
  * [5] P. M. Ashton, C. S. Saunders, G. A. Taylor, A. M. Carter, and M. E. Bradley, “Inertia estimation of the gb power system using synchrophasor measurements,” _IEEE Transactions on Power Systems_ , vol. 30, no. 2, pp. 701–709, March 2015.
  * [6] P. Wall and V. Terzija, “Simultaneous estimation of the time of disturbance and inertia in power systems,” _IEEE Transactions on Power Delivery_ , vol. 29, no. 4, pp. 2018–2031, 2014.
  * [7] G. Cai, B. Wang, D. Yang, Z. Sun, and L. Wang, “Inertia estimation based on observed electromechanical oscillation response for power systems,” _IEEE Transactions on Power Systems_ , vol. 34, no. 6, pp. 4291–4299, 2019\.
  * [8] D. Wilson, J. Yu, N. Al-Ashwal, B. Heimisson, and V. Terzija, “Measuring effective area inertia to determine fast-acting frequency response requirements,” _International Journal of Electrical Power & Energy Systems_, vol. 113, pp. 1 – 8, 2019. [Online]. Available: http://www.sciencedirect.com/science/article/pii/S0142061519304892
  * [9] R. J. Best, P. V. Brogan, and D. J. Morrow, “Power system inertia estimation using hvdc power perturbations,” _IEEE Transactions on Power Systems_ , vol. 36, no. 3, pp. 1890–1899, 2021.
  * [10] N. Hosaka, B. Berry, and S. Miyazaki, “The world’s first small power modulation injection approach for inertia estimation and demonstration in the island grid,” in _2019 8th International Conference on Renewable Energy Research and Applications (ICRERA)_ , 2019, pp. 722–726.
  * [11] M. Brown, M. Biswal, S. Brahma, S. J. Ranade, and H. Cao, “Characterizing and quantifying noise in pmu data,” in _2016 IEEE Power and Energy Society General Meeting (PESGM)_ , 2016, pp. 1–5.
  * [12] J. Susanto, F. Shahnia, D. Sharafi, and L. Kwek, “Estimating a power system’s load relief factor using the high-resolution data of fault recorders,” in _Proc. 9th International Conference on Power and Energy Systems (ICPES 2019)_ , Dec 2019.
  * [13] A. Fernandez-Guillamon, A. Vigueras-Rodriguez, and A. Molina-Garcia, “Analysis of power system inertia estimation in high wind power plant integration scenarios,” _IET Renewable Power Generation_ , vol. 13, no. 15, pp. 2807–2816, 2019. [Online]. Available: https://ietresearch.onlinelibrary.wiley.com/doi/abs/10.1049/iet-rpg.2019.0220
  * [14] N. Masood, R. Yan, and T. K. Saha, “Investigation of load frequency relief from field measurements and its impact on contingency reserve evaluation,” _IEEE Transactions on Power Systems_ , vol. 33, no. 1, pp. 567–577, 2018.
  * [15] M. Davies, F. Moran, and J. Bird, “Power/frequency characteristics of the british grid system,” _Proceedings of the IEE - Part A: Power Engineering_ , vol. 106, pp. 154–162(8), April 1959. [Online]. Available: https://digital-library.theiet.org/content/journals/10.1049/pi-a.1959.0065
  * [16] G. Berg, “System and load behaviour following loss of generation. experimental results and evaluation,” _Proceedings of the Institution of Electrical Engineers_ , vol. 119, pp. 1483–1486(3), October 1972. [Online]. Available: https://digital-library.theiet.org/content/journals/10.1049/piee.1972.0293
  * [17] J. W. O’Sullivan and M. J. O’Malley, “Identification and validation of dynamic global load model parameters for use in power system frequency simulations,” _IEEE Transactions on Power Systems_ , vol. 11, no. 2, pp. 851–857, 1996.
  * [18] R. Pearmine, Y. Song, T. Williams, and A. Chebbo, “Identification of a load–frequency characteristic for allocation of spinning reserves on the british electricity grid,” _IEE Proceedings - Generation, Transmission and Distribution_ , vol. 153, pp. 633–638(5), November 2006. [Online]. Available: https://digital-library.theiet.org/content/journals/10.1049/ip-gtd_20050404
  * [19] H. Omara, “A methodology for determining the load frequency sensitivity,” Ph.D. dissertation, University of Manchester, 2012.
  * [20] “IEEE Standard for a Precision Clock Synchronization Protocol for Networked Measurement and Control Systems,” _IEEE Std 1588-2019 (Revision ofIEEE Std 1588-2008)_ , pp. 1–499, 2020.
  * [21] M. Burth, G. Verghese, and M. Velez-Reyes, “Subset selection for improved parameter estimation in on-line identification of a synchronous generator,” _IEEE Transactions on Power Systems_ , vol. 14, no. 1, pp. 218–225, 1999\.
  * [22] M. Karrari and O. Malik, “Identification of physical parameters of a synchronous generator from online measurements,” _IEEE Transactions on Energy Conversion_ , vol. 19, no. 2, pp. 407–415, 2004.
  * [23] M. A. M. Ariff, B. C. Pal, and A. K. Singh, “Estimating dynamic model parameters for adaptive protection and control in power system,” _IEEE Transactions on Power Systems_ , vol. 30, no. 2, pp. 829–839, 2015.
  * [24] P. Pourbeik, “Automated parameter derivation for power plant models from system disturbance data,” in _2009 IEEE Power Energy Society General Meeting_ , 2009, pp. 1–10.
  * [25] M. R. B. Tavakoli, M. Power, L. Ruttledge, and D. Flynn, “Load inertia estimation using white and grey-box estimators for power systems with high wind penetration,” _IFAC Proceedings Volumes_ , vol. 45, no. 21, pp. 399–404, 2012, 8th Power Plant and Power System Control Symposium. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S1474667016320031
  * [26] Y. Bian, H. Wyman-Pain, F. Li, R. Bhakar, S. Mishra, and N. P. Padhy, “Demand side contributions for system inertia in the gb power system,” _IEEE Transactions on Power Systems_ , vol. 33, no. 4, pp. 3521–3530, 2018.


alignment phase is significantly prolonged as the SMBH mass grows more slowly,
leading to a smaller truncation radius. This leads to a reduced precession
angle in the second solid-body precession phase. For the secondary SMBH, on
the other hand, the smaller truncation radius means that the BH-disc system
remains in the Bardeen-Petterson configuration until $t\sim 40$ Myr and then
transfers to the solid body precession regime with a significantly reduced
precession angle.

Finally, we also consider the modified set-up, where the precession rates are
reduced by 95 per cent in the truncated disc mode, uniadios_slowprec. As with
the bin_q1i0 set-up, this means that the SMBH versors are lagging behind the
total versor alignment due to external inflows throughout, so that the
precession angle increases, leading to a larger final offset angle for the
primary SMBH. The secondary SMBH experiences a similar albeit much smaller
effect as it spends significantly less time in the solid body precession
phase.

####  4.3.5 Summary and implications of angular momentum evolution

Figures 10 and 13 clearly show that the accretion disc subgrid model has a
significant impact on the spin alignment and precession process. In the
truncated disc regime, precession is much more prevalent with the SMBH spin
versors covering a large solid angle. This has crucial implications for
distributed feedback, in the particular in the context of spin-driven jets
(see Section 5.2.2 for a more detailed discussion of implications for AGN
feedback). Also note that VLA and Merlin observations have uncovered
precessing jets with inferred precession periods of $1-10$ Myr (Krause et al.,
2019). Whilst this could be attributed to geodetic precession, it would also
be consistent with solid-body precession in the truncated disc regime within
our framework.

Furthermore, it is interesting to note that despite the alignment torques
between the accretion disc and SMBH being extremely weak in the ADAF/ADIOS
regime, the SMBH may still realign with the disc if the external inflow
timescale matches the precession timescale. Furthermore, Bardeen-Petterson
alignment may be accelerated if the solid body precession reorientates the
SMBH versor towards the mini discs. On the other hand, if the precession
process is slowed down significantly due to additional torques in the
truncated disc state (Bollimpalli et al., 2023), then the spins will not align
within the inspiral timescale, with crucial implications for recoil
predictions, influencing merger rates and the gravitational wave background.

##  5 Discussion

We have developed a unified accretion disc model suitable for galaxy formation
simulations that can resolve the outer edge of the accretion disc. Our work
builds on the thin disc model prescription presented in Fiacconi et al. (2018)
and extends this scheme by harnessing the results from high-resolution
GR(R)MHD simulations of accretion discs. Our model spans a range of radiative
regimes from the optically thin, geometrically thick discs to the optically
thick, geometrically thin discs via a truncated disc configuration. We also
track the SMBH spin evolution with a disc-state-dependent Lense-Thirring
precession model. Our predictions for the SMBH spin alignment and precession
as well as SMBH luminosities are highly sensitive to the disc state with
important implications for future multi-messenger observational campaigns
which we outline in Section 5.1. We also discuss the caveats of our model and
important future extensions in Section 5.2.

###  5.1 Implications for observations

Parsec-scale SMBH binaries have proved extremely challenging to observe, with
only a handful of detections so far (Rodriguez et al., 2006; Bansal et al.,
2017; Kharb et al., 2017), as these are in the regime were extremely high
angular resolution is necessary to resolve the binary, however, the SMBHs are
not close enough to display significant variability on observable timescales.
This makes radio interferometers the tool of choice for searching for these
binaries and future facilities such as SKA or ngVLA are expected to
significantly increase sample sizes (e.g. Burke-Spolaor et al., 2018), though
see a cautionary note from D’Orazio & Charisi (2023).

To estimate radio luminosities, it is crucial to have reliable estimates for
the mass accretion rates, spin orientations and solid angles (e.g. Mangiagli
et al., 2022). With our unified accretion disc model, we can estimate all of
the relevant quantities and make predictions for the jet luminosities
depending on the disc state which reveals a much more promising avenue for
detecting SMBH binaries accreting at intermediate luminosities than the thin
disc modelling (see Figure 9).

According to analytical calculations, our binary system should have an
inspiral timescale of $\lesssim 10^{3}$Myr (see Figure 15), so this represents
a precursor to the gravitational-wave-emitting binaries. Nevertheless it is
crucial to make predictions for these precursor luminosities, in particular
for X-ray confirmations (Piro et al., 2023).

Furthermore, accurate AGN modelling will be crucial for dual AGN (and close
pairs) searches with JWST where recent tantalizing observations have revealed
hints at a large fraction of candidate dual AGN (Maiolino et al., 2023b; Übler
et al., 2023b), highlighting the importance of imaging spectroscopy in
tracking SMBH mergers throughout cosmic time.

Our spin evolution predictions also play an important role in estimating the
gravitational wave background signal as they determine the recoil velocity of
the SMBH remnant – in particular large and misaligned SMBH spins may lead to
recoil velocities of $\sim 5000\ \mathrm{km\,s^{-1}}$ (Lousto & Zlochower,
2011, 2013; Lousto & Healy, 2019) which would have a significant impact on
SMBH occupation fractions and merger rates.

###  5.2 Caveats and future work

####  5.2.1 Accretion modelling

We have extended the thin disc model from Fiacconi et al. (2018) into the
heavily sub-Eddington regime based on the ADAF/ADIOS models. However, we note
that there are several other disc models that have been proposed for the
radiatively inefficient accretion regime, including convection-dominated
accretion flows (CDAFs; Narayan et al., 2000; Quataert & Gruzinov, 2000) or
more recently simple convective flows (Xu, 2023) and strongly magnetised Bondi
models (Cho et al., 2023). Li et al. (2013) carried out an extensive
simulation suite, capturing the accretion flows from Bondi scale to the small
scale and demonstrated how the inflow and outflow patterns will change as a
function of cooling. Similarly, they show how strong outflows suppress the
Bondi rate, however, they also demonstrate a strong dependency on cooling.
Note that in principle, we could incorporate these accretion flow scenarios
within our framework by varying the $s$ parameter from the ADIOS wind loss
model, which effectively sets the flow density profile. However, this is
beyond the scope of this paper and we plan to investigate the impact of $s$
parameter choices more systematically in the future.

Whilst this work has focused on the radiatively-inefficient regime, in the
future it will also be important to further consider the disc modelling in the
radiatively efficient regime. Several recent studies have questioned and
improved the standard thin $\alpha$-disc model (e.g. Jiang et al., 2019) and
this will require careful consideration to improve the predictions of our
model in the quasar regime.

Misaligned discs experience radius-dependent Lense-Thirring torques. How these
differential torques are communicated depends on the accretion disc regime. In
our unified accretion disc model, we evolve the SMBH spin according to the
Bardeen-Petterson alignment (Bardeen & Petterson, 1975) in the thin disc state
and following the solid body precession model from Fragile et al. (2007) in
the thick disc state. We also include the possibility of slowed down solid
body precession in the truncated disc regime (Bollimpalli et al., 2023).
However, GR(R)MHD simulations paint an even more complex picture, including
disc tearing (e.g. Liska et al., 2023b) which could suppress Bardeen-Petterson
alignment (e.g. Steinle & Gerosa, 2023) and induce phase lags between the
inner hot flow and outer thin disc (e.g. Liska et al., 2023a). Furthermore, if
the thick disc is in the magnetically arrested state, strong magneto-spin
alignment could significantly suppress solid body precession (e.g. McKinney et
al., 2013; Chatterjee et al., 2023).

Finally, another crucial future extension to our model will be to incorporate
the super-Eddington regime. This accretion disc regime may be particularly
important in the early Universe, where JWST is starting to uncover a
population of overmassive AGN (compared to the local scaling relations) which
could stem from sustained super-Eddington accretion episodes (e.g. Goulding et
al., 2023; Maiolino et al., 2023b). Note that many of the key physical
properties of our ADIOS modelling in the sub-Eddington regime will also apply
to the super-Eddington regime. The latter is advection-dominated due to the
long radiative diffusion times (rather than long cooling times, see Yuan &
Narayan, 2014), so that the radiative efficiency drops in the super-Eddington
regime and accretion (as well as jets) cause the SMBH to spin down (also see
Lowell et al., 2023; Ricarte et al., 2023).

####  5.2.2 Feedback modelling

Within our current framework we only model gas accretion through a unified
accretion disc model and the obvious next step is to incorporate AGN feedback
mechanisms into our model. The most straightforward extension here would be to
incorporate ADIOS winds with the mass loading directly informed by the wind
mass loss in our model. The momentum and energy loading of the wind could then
be inferred from the kinetic wind efficiencies obtained from GRMHD simulations
(e.g. Tchekhovskoy et al., 2012; Sądowski et al., 2013; Sądowski & Gaspari,
2017). We note that the thin disc component is also expected to drive winds,
in particular at high Eddington fractions, and including these winds would be
an important extension for future work (also see Sala et al., 2021).

Moreover, AGN jets are an ubiquitous phenomenon and are likely powered by the
SMBH spin energy extraction via the Blandford-Znajek (BZ) mechanism (Blandford
& Znajek, 1977). Talbot et al. (2021) have coupled spin-driven BZ jets self-
consistently with the thin disc model from Fiacconi et al. (2018) and
considered their evolution in Seyfert galaxies (Talbot et al., 2022) and a
galaxy merger (Talbot et al., 2023). In future work, it would be timely to
extend this scheme to include spin-driven jets across different accretion disc
types (also see Huško et al., 2022). In particular, the solid body precession
model could have a crucial impact here leading to more ‘distributed’ jet
feedback. Note that the precession periods in the solid body regime are highly
uncertain (see Figure 4), and there are hints at Lense-Thirring precession in
AGN on the scale of $\sim 10$ yr (Cui et al., 2023; von Fellenberg et al.,
2023). In terms of AGN subgrid models, this would imply that for timescales
relevant for galaxy formation the AGN feedback would be effectively
distributed over a cone, which may significantly increase its impact.

####  5.2.3 Porting the unified accretion disc model to cosmological
simulations

In addition to improving the subgrid SMBH modelling, it will also be crucial
to move towards more realistic environments. We found that more chaotic
accretion, such as in the misaligned binary case, can have a significant
impact on luminosities and spin evolution so a natural next step would be to
test this framework within galaxy merger simulations (also see Liao et al.,
2023; Talbot et al., 2023). Furthermore, extending our model to cosmological
zoom-in simulations will allow us to track the mass and angular momentum flows
from the cosmological environment down to the scale of the accretion disc. For
the zoom-in simulations, the super-Lagrangian refinement technique will permit
us to reach the required resolution (also see Curtis & Sijacki, 2016; Bourne
et al., 2019; Bourne & Sijacki, 2021; Anglés-Alcázar et al., 2021; Hopkins et
al., 2023a, b). For large cosmological volumes, applying the super-Lagrangian
refinement technique to the central regions of tens of thousands of galaxies,
however, is computationally infeasible. One avenue to circumvent these issues
could involve using machine-learning techniques such as generative adversarial
networks or learned coarse models (e.g. Li et al., 2021; Stachenfeld et al.,
2021), though this would require appropriate training data sets covering the
vast range of relevant scales and different environments.

##  6 Conclusions

We have developed a novel unified accretion disc model designed to be embedded
in galaxy formation simulations. Drawing inspiration from high-resolution
GR(R)MHD simulations, our model tracks mass flows through three distinct disc
states and evolves the SMBH spin based on disc-state-dependent Lense-Thirring
precession models. This innovative accretion disc framework has been
integrated into the moving mesh code arepo, and our work includes a
comprehensive suite of validation simulations. These simulations encompass
idealized scenarios featuring single SMBHs and binary SMBHs embedded in
gaseous discs, with variations in mass ratio, inclination angle, and subgrid
accretion disc models. Our key findings are summarized below.

  * •

The ADIOS model predicts substantial wind loss in the thick disc regime,
leading to a significant suppression of gas mass accretion rates onto SMBHs by
an order of magnitude or more. Additionally, the decrease in radiative
efficiency within the advection-dominated accretion regime results in
luminosity estimates that are decreased by another order of magnitude.
Consequently, predictions for the luminosities of electromagnetic counterparts
of SMBH binaries exhibit significant sensitivity to the disc state.

  * •

Hence, it is imperative to incorporate these realistic accretion disc models
when forecasting outcomes for future missions both for single and binary
SMBHs. Specifically, our unified accretion disc model reveals a drastic
reduction in the redshift range for accretion disc luminosities compared to
predictions based on the thin $\alpha$-disc model, regardless of whether winds
are accounted for. This has important implications for SMBH binary searches in
the X-ray, optical and infrared. Furthermore, the predicted spectrum undergoes
significant shifts, with the luminosity being dominated by the hard inner hot
flow.

  * •

For the radio jet luminosities, on the other hand, the situation is reversed
with the unified accretion disc model predictions much more favourable than
the thin disc predictions, since the geometrically thick ADIOS flow leads to
significantly more efficient jet production. We predict that future VLBI
searches with SKA and the ngVLA should significantly increase the known
samples of parsec-scale binaries.

  * •

The evolution of the spin magnitude depends only mildly on the disc state,
primarily due to relatively low external inflow rates in all cases considered
here. We note however, that depending on the disc state adopted we predict for
the same large-scale gas inflow, the SMBH spin may either increase or
decrease. Clearly, in case of higher and more variable accretion rates,
significant differences in the spin magnitude evolution are expected.

  * •

Importantly, the SMBH spin direction’s evolution is markedly influenced by the
assumed disc model, ranging from slow alignment in the Bardeen-Petterson
configuration to fast solid-body precession in the truncated disc regime. This
has important implications for the impact of jet-mode feedback on galaxies and
the circumgalactic medium, as more rapid precession would result in
distributing the jet energy over a larger volume.

  * •

For binary systems, the variations in spin evolution could significantly
impact recoil velocities, consequently influencing BH merger rates and the
gravitational wave background.

This model provides a crucial advancement in AGN accretion modelling, and
there are numerous promising avenues for extensions, including spatially
resolved feedback. In the near future, next generation cosmological
simulations, incorporating this type of subgrid physics, will be able to give
us unprecedented insights into the complex interplay of SMBHs and their host
galaxies, and provide the scientific community with much needed robust
predictions for the era of multi-messenger science.

## Acknowledgements

The authors would like to thank Alexander Tchekhovskoy, Beverly Lowell,
Matthew Liska, Jenny Greene, Shy Genel, Angelo Ricarte, Doug Rennehan and Sam
Turner for useful discussions during the development of this work. The authors
are also grateful to Lucy Reading-Ikkanda at the Simons Foundation for
creating the illustration of our model in Figure 1. S.K. is supported by a
Flatiron Research Fellowship at the Flatiron Institute, a division of the
Simons Foundation and a Junior Research Fellowship from St Catharine’s
College, Cambridge. S.K, M.A.B., and D.S. acknowledge support by European
Research Council Starting Grant 638707 “Black holes and their host galaxies:
coevolution across cosmic time”. M.A.B and D.S. additionally acknowledge
support from the Science and Technology Facilities Council (STFC). K.P. is
supported by the Simons-NSBP Scholars Program. This work was supported by the
Simons Collaboration on “Learning the Universe”. The computations in this work
were, in part, run at facilities supported by the Scientific Computing Core at
the Flatiron Institute, a division of the Simons Foundation. This work used
the DiRAC Data Intensive service (CSD3) at the University of Cambridge,
managed by the University of Cambridge University Information Services on
behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). The DiRAC component
of CSD3 at Cambridge was funded by BEIS, UKRI and STFC capital funding and
STFC operations grants. This work used the DiRAC Memory Intensive service
(Cosma8) at Durham University, managed by the Institute for Computational
Cosmology on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). The
DiRAC service at Durham was funded by BEIS, UKRI and STFC capital funding,
Durham University and STFC operations grants. DiRAC is part of the UKRI
Digital Research Infrastructure.

## Data Availability

The data underlying this article will be shared on reasonable request to the
corresponding author.

## References

  * Abbott et al. (2020) Abbott R., et al., 2020, PRL, 125, 101102
  * Agazie et al. (2023a) Agazie G., et al., 2023a, ApJ, 951, L8
  * Agazie et al. (2023b) Agazie G., et al., 2023b, ApJ, 952, L37
  * Anglés-Alcázar et al. (2013) Anglés-Alcázar D., Özel F., Davé R., 2013, ApJ, 770, 5
  * Anglés-Alcázar et al. (2017) Anglés-Alcázar D., Davé R., Faucher-Giguère C.-A., Özel F., Hopkins P. F., 2017, MNRAS, 464, 2840
  * Anglés-Alcázar et al. (2021) Anglés-Alcázar D., et al., 2021, ApJ, 917, 53
  * Antoniadis et al. (2023) Antoniadis J., et al., 2023, A&A, 678, A50
  * Aravindan et al. (2023) Aravindan A., Liu W., Canalizo G., Veilleux S., Bohn T., Sexton R. O., Rupke D. S. N., U V., 2023, ApJ, 950, 33
  * Balbus (1991) Balbus S. A., 1991, ApJ, 372, 25
  * Balbus & Hawley (1998) Balbus S. A., Hawley J. F., 1998, Reviews of Modern Physics, 70, 1
  * Bambic et al. (2023) Bambic C. J., Russell H. R., Reynolds C. S., Fabian A. C., McNamara B. R., Nulsen P. E. J., 2023, MNRAS, 522, 4374
  * Bansal et al. (2017) Bansal K., Taylor G. B., Peck A. B., Zavala R. T., Romani R. W., 2017, ApJ, 843, 14
  * Barausse (2012) Barausse E., 2012, MNRAS, 423, 2533
  * Bardeen & Petterson (1975) Bardeen J. M., Petterson J. A., 1975, ApJL, 195, L65
  * Bardeen et al. (1972) Bardeen J. M., Press W. H., Teukolsky S. A., 1972, ApJ, 178, 347
  * Beckmann et al. (2019) Beckmann R. S., et al., 2019, A&A, 631, A60
  * Bennett & Sijacki (2020) Bennett J. S., Sijacki D., 2020, MNRAS, 499, 597
  * Berti & Volonteri (2008) Berti E., Volonteri M., 2008, ApJ, 684, 822
  * Blandford & Begelman (1999) Blandford R. D., Begelman M. C., 1999, MNRAS, 303, L1
  * Blandford & Payne (1982) Blandford R. D., Payne D. G., 1982, MNRAS, 199, 883
  * Blandford & Znajek (1977) Blandford R. D., Znajek R. L., 1977, MNRAS, 179, 433
  * Bogdanović et al. (2022) Bogdanović T., Miller M. C., Blecha L., 2022, LRR, 25, 3
  * Bollati et al. (2023a) Bollati F., Lupi A., Dotti M., Haardt F., 2023a, MNRAS, 520, 3696
  * Bollati et al. (2023b) Bollati F., Lupi A., Dotti M., Haardt F., 2023b, arXiv e-prints, arXiv.2311.07576
  * Bollimpalli et al. (2023) Bollimpalli D. A., Fragile P. C., Kluźniak W., 2023, MNRAS, 520, L79
  * Bondi (1952) Bondi H., 1952, MNRAS, 112, 195
  * Bondi & Hoyle (1944) Bondi H., Hoyle F., 1944, MNRAS, 104, 273
  * Bourne & Sijacki (2021) Bourne M. A., Sijacki D., 2021, MNRAS, 506, 488
  * Bourne et al. (2019) Bourne M. A., Sijacki D., Puchwein E., 2019, MNRAS, 490, 343
  * Bourne et al. (2023) Bourne M. A., Fiacconi D., Sijacki D., Piotrowska J. M., Koudmani S., 2023, arXiv e-prints, arXiv:2311.17144
  * Braun et al. (2019) Braun R., Bonaldi A., Bourke T., Keane E., Wagg J., 2019, arXiv e-prints, arXiv.1912.12699
  * Burke-Spolaor et al. (2018) Burke-Spolaor S., et al., 2018, arXiv e-prints, arXiv.1808.04368
  * Bustamante & Springel (2019) Bustamante S., Springel V., 2019, MNRAS, 490, 4133
  * Cenci et al. (2021) Cenci E., Sala L., Lupi A., Capelo P. R., Dotti M., 2021, MNRAS, 500, 3719
  * Chatterjee et al. (2023) Chatterjee K., Liska M., Tchekhovskoy A., Markoff S., 2023, arXiv e-prints, arXiv.2311.00432
  * Cho & Narayan (2022) Cho H., Narayan R., 2022, ApJ, 932, 97
  * Cho et al. (2023) Cho H., Prather B. S., Narayan R., Natarajan P., Su K.-Y., Ricarte A., Chatterjee K., 2023, arXiv e-prints, arXiv.2310.19135
  * Civano et al. (2023) Civano F., et al., 2023, arXiv e-prints, arXiv.2311.04832
  * Cui et al. (2023) Cui Y., et al., 2023, Nature, 621, 711
  * Curtis & Sijacki (2015) Curtis M., Sijacki D., 2015, MNRAS, 454, 3445
  * Curtis & Sijacki (2016) Curtis M., Sijacki D., 2016, MNRAS, 457, L34
  * D’Orazio & Charisi (2023) D’Orazio D. J., Charisi M., 2023, arXiv e-prints, arXiv.2310.16896
  * Davis et al. (2022) Davis F., et al., 2022, MNRAS, 511, 4109
  * De Rosa et al. (2019) De Rosa A., et al., 2019, New Astronomy Reviews, 86, 101525
  * Dittmann & Ryan (2022) Dittmann A. J., Ryan G., 2022, MNRAS, 513, 6158
  * Dittmann & Ryan (2023) Dittmann A. J., Ryan G., 2023, arXiv e-prints, arXiv.2310.07758
  * Dotti et al. (2013) Dotti M., Colpi M., Pallini S., Perego A., Volonteri M., 2013, ApJ, 762, 68
  * Dubois et al. (2014) Dubois Y., Volonteri M., Silk J., Devriendt J., Slyz A., 2014, MNRAS, 440, 2333
  * Dubois et al. (2021) Dubois Y., et al., 2021, A&A, 651, A109
  * Esin et al. (1997) Esin A. A., McClintock J. E., Narayan R., 1997, ApJ, 489, 865
  * Fanidakis et al. (2011) Fanidakis N., Baugh C. M., Benson A. J., Bower R. G., Cole S., Done C., Frenk C. S., 2011, MNRAS, 410, 53
  * Fiacconi et al. (2018) Fiacconi D., Sijacki D., Pringle J. E., 2018, MNRAS, 477, 3807
  * Fragile et al. (2007) Fragile P. C., Blaes O. M., Anninos P., Salmonson J. D., 2007, ApJ, 668, 417
  * Franchini et al. (2022) Franchini A., Lupi A., Sesana A., 2022, ApJ, 929, L13
  * Gaspari & Sądowski (2017) Gaspari M., Sądowski A., 2017, ApJ, 837, 149
  * Giustini & Proga (2019) Giustini M., Proga D., 2019, A&A, 630, A94
  * Goulding et al. (2023) Goulding A. D., et al., 2023, ApJ, 955, L24
  * Greene et al. (2020) Greene J. E., Strader J., Ho L. C., 2020, ARA&A, 58, 257
  * Guo et al. (2023) Guo M., Stone J. M., Kim C.-G., Quataert E., 2023, ApJ, 946, 26
  * Haidar et al. (2022) Haidar H., et al., 2022, MNRAS
  * Harikane et al. (2023) Harikane Y., et al., 2023, ApJ, 959, 39
  * Hogg & Reynolds (2018) Hogg J. D., Reynolds C. S., 2018, ApJ, 854, 6
  * Honma (1996) Honma F., 1996, PASJ, 48, 77
  * Hopkins & Quataert (2011) Hopkins P. F., Quataert E., 2011, MNRAS, 415, 1027
  * Hopkins et al. (2023a) Hopkins P. F., et al., 2023a, arXiv e-prints, arXiv.2309.13115
  * Hopkins et al. (2023b) Hopkins P. F., et al., 2023b, arXiv e-prints, arXiv.2310.04506
  * Hoyle & Lyttleton (1939) Hoyle F., Lyttleton R. A., 1939, Math. Proc. Camb. Phil. Soc., 35, 405
  * Huško et al. (2022) Huško F., Lacey C. G., Schaye J., Schaller M., Nobels F. S. J., 2022, MNRAS, 516, 3750
  * Inayoshi et al. (2019) Inayoshi K., Ichikawa K., Ostriker J. P., Kuiper R., 2019, MNRAS, 486, 5377
  * Ingram & Done (2011) Ingram A., Done C., 2011, MNRAS, 415, 2323
  * Ingram & Motta (2019) Ingram A. R., Motta S. E., 2019, New Astronomy Reviews, 85, 101524
  * Ingram et al. (2009) Ingram A., Done C., Fragile P. C., 2009, MNRAS, 397, L101
  * Jiang & Blaes (2020) Jiang Y.-F., Blaes O., 2020, ApJ, 900, 25
  * Jiang et al. (2019) Jiang Y.-F., Blaes O., Stone J. M., Davis S. W., 2019, ApJ, 885, 144
  * Kharb et al. (2017) Kharb P., Lal D. V., Merritt D., 2017, Nature Astronomy, 1, 727
  * King et al. (2005) King A. R., Lubow S. H., Ogilvie G. I., Pringle J. E., 2005, MNRAS, 363, 49
  * Kormendy & Ho (2013) Kormendy J., Ho L. C., 2013, ARA&A, 51, 511
  * Koudmani et al. (2021) Koudmani S., Henden N. A., Sijacki D., 2021, MNRAS, 503, 3568
  * Koudmani et al. (2022) Koudmani S., Sijacki D., Smith M. C., 2022, MNRAS, 516, 2112
  * Krause et al. (2019) Krause M. G. H., et al., 2019, MNRAS, 482, 240
  * Lagos et al. (2009) Lagos C. D. P., Padilla N. D., Cora S. A., 2009, MNRAS, 395, 625
  * Lagos et al. (2023) Lagos C. D. P., et al., 2023, arXiv e-prints, arXiv.2309.02310
  * Lalakos et al. (2022) Lalakos A., et al., 2022, ApJ, 936, L5
  * Lalakos et al. (2023) Lalakos A., Tchekhovskoy A., Bromberg O., Gottlieb O., Jacquemin-Ide J., Liska M., Zhang H., 2023, arXiv e-prints, 2310.11487
  * Latif et al. (2024) Latif M. A., Whalen D. J., Mezcua M., 2024, MNRAS, 527, L37
  * Lense & Thirring (1918) Lense J., Thirring H., 1918, Physikalische Zeitschrift, 19, 156
  * Li et al. (2013) Li J., Ostriker J., Sunyaev R., 2013, ApJ, 767, 105
  * Li et al. (2021) Li Y., Ni Y., Croft R. A. C., Di Matteo T., Bird S., Feng Y., 2021, Proc Natl Acad Sci USA, 118, e2022038118
  * Liao et al. (2023) Liao S., et al., 2023, MNRAS, 520, 4463
  * Liska et al. (2018) Liska M., Hesp C., Tchekhovskoy A., Ingram A., van der Klis M., Markoff S., 2018, MNRAS, 474, L81
  * Liska et al. (2019) Liska M., Tchekhovskoy A., Ingram A., van der Klis M., 2019, MNRAS, 487, 550
  * Liska et al. (2022) Liska M. T. P., Musoke G., Tchekhovskoy A., Porth O., Beloborodov A. M., 2022, ApJ, 935, L1
  * Liska et al. (2023a) Liska M., Hesp C., Tchekhovskoy A., Ingram A., van der Klis M., Markoff S. B., 2023a, New Astronomy, 101, 102012
  * Liska et al. (2023b) Liska M. T. P., Kaaz N., Musoke G., Tchekhovskoy A., Porth O., 2023b, ApJ, 944, L48
  * Liska et al. (2023c) Liska M. T. P., Kaaz N., Chatterjee K., Emami R., Musoke G., 2023c, arXiv e-prints, arXiv.2309.15926
  * Liu & Taam (2009) Liu B. F., Taam R. E., 2009, ApJ, 707, 233
  * Liu et al. (2020) Liu W., Veilleux S., Canalizo G., Rupke D. S. N., Manzano-King C. M., Bohn T., U V., 2020, ApJ, 905, 166
  * Lodato & Pringle (2007) Lodato G., Pringle J. E., 2007, MNRAS, 381, 1287
  * Lousto & Healy (2019) Lousto C. O., Healy J., 2019, Physical Review D, 100, 104039
  * Lousto & Zlochower (2011) Lousto C. O., Zlochower Y., 2011, Physical Review Letters, 107, 231102
  * Lousto & Zlochower (2013) Lousto C. O., Zlochower Y., 2013, Physical Review D, 87, 084027
  * Lowell et al. (2023) Lowell B., Jacquemin-Ide J., Tchekhovskoy A., Duncan A., 2023, arXiv e-prints, arXiv.2302.01351
  * Maiolino et al. (2023a) Maiolino R., et al., 2023a, arXiv e-prints, arXiv.2305.12492
  * Maiolino et al. (2023b) Maiolino R., et al., 2023b, arXiv e-prints, arXiv.2308.01230
  * Mangiagli et al. (2022) Mangiagli A., Caprini C., Volonteri M., Marsat S., Vergani S., Tamanini N., Inchauspé H., 2022, Physical Review D, 106, 103017
  * Manzano-King et al. (2019) Manzano-King C. M., Canalizo G., Sales L. V., 2019, ApJ, 884, 54
  * Martin et al. (2007) Martin R. G., Pringle J. E., Tout C. A., 2007, MNRAS, 381, 1617
  * Martin et al. (2019) Martin R. G., Nixon C. J., Pringle J. E., Livio M., 2019, New Astronomy, 70, 7
  * Massonneau et al. (2023) Massonneau W., Dubois Y., Volonteri M., Beckmann R. S., 2023, A&A, 669, A143
  * McKinney et al. (2012) McKinney J. C., Tchekhovskoy A., Blandford R. D., 2012, MNRAS, 423, 3083
  * McKinney et al. (2013) McKinney J. C., Tchekhovskoy A., Blandford R. D., 2013, Science, 339, 49
  * Meier (2001) Meier D. L., 2001, ApJ, 548, L9
  * Moesta et al. (2012) Moesta P., Alic D., Rezzolla L., Zanotti O., Palenzuela C., 2012, ApJL, 749, L32
  * Motta et al. (2018) Motta S. E., Franchini A., Lodato G., Mastroserio G., 2018, MNRAS, 473, 431
  * Narayan & McClintock (2008) Narayan R., McClintock J. E., 2008, New Astronomy Reviews, 51, 733
  * Narayan & Yi (1995) Narayan R., Yi I., 1995, ApJ, 452, 710
  * Narayan et al. (2000) Narayan R., Igumenshchev I. V., Abramowicz M. A., 2000, ApJ, 539, 798
  * Narayan et al. (2022) Narayan R., Chael A., Chatterjee K., Ricarte A., Curd B., 2022, MNRAS, 511, 3795
  * Nemmen et al. (2014) Nemmen R. S., Storchi-Bergmann T., Eracleous M., 2014, MNRAS, 438, 2804
  * Papaloizou & Pringle (1983) Papaloizou J. C. B., Pringle J. E., 1983, MNRAS, 202, 1181
  * Penny et al. (2018) Penny S. J., et al., 2018, MNRAS, 476, 979
  * Perego et al. (2009) Perego A., Dotti M., Colpi M., Volonteri M., 2009, MNRAS, 399, 2249
  * Pfeifle et al. (2023a) Pfeifle R., Weaver K., Cann J., Ogorzalek A., Bogdan A., Chakraborty P., 2023a, BAAS, 55, 110.27
  * Pfeifle et al. (2023b) Pfeifle R. W., et al., 2023b, arXiv e-prints, arXiv.2311.05154
  * Pillepich et al. (2018) Pillepich A., et al., 2018, MNRAS, 473, 4077
  * Piro et al. (2023) Piro L., et al., 2023, MNRAS, 521, 2577
  * Plotkin & Reines (2018) Plotkin R. M., Reines A. E., 2018, arXiv e-prints, arXiv.1810.06814
  * Power et al. (2011) Power C., Nayakshin S., King A., 2011, MNRAS, 412, 269
  * Pringle (1981) Pringle J. E., 1981, ARA&A, 19, 137
  * Quataert & Gruzinov (2000) Quataert E., Gruzinov A., 2000, ApJ, 539, 809
  * Reardon et al. (2023) Reardon D. J., et al., 2023, ApJL, 951, L6
  * Rennehan et al. (2023) Rennehan D., Babul A., Moa B., Davé R., 2023, arXiv e-prints, arXiv.2309.15898
  * Reynolds (2019) Reynolds C. S., 2019, Nature Astronomy, 3, 41
  * Reynolds et al. (2023) Reynolds C. S., et al., 2023, arXiv e-prints, arXiv.2311.00780
  * Ricarte et al. (2023) Ricarte A., Narayan R., Curd B., 2023, ApJ, 954, L22
  * Rodriguez et al. (2006) Rodriguez C., Taylor G. B., Zavala R. T., Peck A. B., Pollack L. K., Romani R. W., 2006, ApJ, 646, 49
  * Rosas-Guevara et al. (2015) Rosas-Guevara Y. M., et al., 2015, MNRAS, 454, 1038
  * Russell et al. (2015) Russell H. R., Fabian A. C., McNamara B. R., Broderick A. E., 2015, MNRAS, 451, 588
  * Ryan et al. (2017) Ryan B. R., Ressler S. M., Dolence J. C., Tchekhovskoy A., Gammie C., Quataert E., 2017, ApJL, 844, L24
  * Sala et al. (2021) Sala L., Cenci E., Capelo P. R., Lupi A., Dotti M., 2021, MNRAS, 500, 4788
  * Schaye et al. (2015) Schaye J., et al., 2015, MNRAS, 446, 521
  * Sesana et al. (2014) Sesana A., Barausse E., Dotti M., Rossi E. M., 2014, ApJ, 794, 104
  * Shakura & Sunyaev (1973) Shakura N. I., Sunyaev R. A., 1973, A&A, 24, 337
  * Shen et al. (2020) Shen X., Hopkins P. F., Faucher-Giguère C.-A., Alexander D. M., Richards G. T., Ross N. P., Hickox R. C., 2020, MNRAS, 495, 3252
  * Sijacki et al. (2015) Sijacki D., Vogelsberger M., Genel S., Springel V., Torrey P., Snyder G. F., Nelson D., Hernquist L., 2015, MNRAS, 452, 575
  * Siwek et al. (2023a) Siwek M., Weinberger R., Muñoz D. J., Hernquist L., 2023a, MNRAS, 518, 5059
  * Siwek et al. (2023b) Siwek M., Weinberger R., Hernquist L., 2023b, MNRAS, 522, 2707
  * Somerville & Davé (2015) Somerville R. S., Davé R., 2015, ARA&A, 53, 51
  * Stachenfeld et al. (2021) Stachenfeld K., et al., 2021, arXiv e-prints, arXiv.2112.15275
  * Steinle & Gerosa (2023) Steinle N., Gerosa D., 2023, MNRAS, 519, 5031
  * Sądowski & Gaspari (2017) Sądowski A., Gaspari M., 2017, MNRAS, 468, 1398
  * Sądowski et al. (2013) Sądowski A., Narayan R., Penna R., Zhu Y., 2013, MNRAS, 436, 3856
  * Sądowski et al. (2016) Sądowski A., Lasota J.-P., Abramowicz M. A., Narayan R., 2016, MNRAS, 456, 3915
  * Talbot et al. (2021) Talbot R. Y., Bourne M. A., Sijacki D., 2021, MNRAS, 504, 3619
  * Talbot et al. (2022) Talbot R. Y., Sijacki D., Bourne M. A., 2022, MNRAS, 514, 4535
  * Talbot et al. (2023) Talbot R. Y., Sijacki D., Bourne M. A., 2023, arXiv e-prints, arXiv.2306.07316
  * Tartėnas & Zubovas (2022) Tartėnas M., Zubovas K., 2022, MNRAS, 516, 2522
  * Tchekhovskoy et al. (2010) Tchekhovskoy A., Narayan R., McKinney J. C., 2010, ApJ, 711, 50
  * Tchekhovskoy et al. (2012) Tchekhovskoy A., McKinney J. C., Narayan R., 2012, JPhCS, 372, 012040
  * The International Pulsar Timing Array Collaboration et al. (2023) The International Pulsar Timing Array Collaboration et al., 2023, arXiv e-prints, arXiv.2309.00693
  * Trump et al. (2011) Trump J. R., et al., 2011, ApJ, 733, 60
  * Turner & Reynolds (2023) Turner S. G. D., Reynolds C. S., 2023, MNRAS, 525, 2287
  * Vasudevan & Fabian (2007) Vasudevan R. V., Fabian A. C., 2007, MNRAS, 381, 1235
  * Verbiest et al. (2016) Verbiest J. P. W., et al., 2016, MNRAS, 458, 1267
  * Vogelsberger et al. (2020) Vogelsberger M., Marinacci F., Torrey P., Puchwein E., 2020, Nat Rev Phys, 2, 42
  * Volonteri et al. (2005) Volonteri M., Madau P., Quataert E., Rees M. J., 2005, ApJ, 620, 69
  * Volonteri et al. (2013) Volonteri M., Sikora M., Lasota J. P., Merloni A., 2013, ApJ, 775, 94
  * Xie & Yuan (2012) Xie F.-G., Yuan F., 2012, MNRAS, 427, 1580
  * Xu (2023) Xu W., 2023, ApJ, 954, 180
  * Yu et al. (2011) Yu Z., Yuan F., Ho L. C., 2011, ApJ, 726, 87
  * Yuan & Narayan (2014) Yuan F., Narayan R., 2014, ARA&A, 52, 529
  * Yuan et al. (2012) Yuan F., Wu M., Bu D., 2012, ApJ, 761, 129
  * Yuan et al. (2018) Yuan F., Yoon D., Li Y.-P., Gan Z.-M., Ho L. C., Guo F., 2018, ApJ, 857, 121
  * von Fellenberg et al. (2023) von Fellenberg S. D., Janssen M., Davelaar J., Zajaček M., Britzen S., Falcke H., Körding E., Ros E., 2023, A&A, 672, L5
  * Übler et al. (2023a) Übler H., et al., 2023a, A&A, 677, A145
  * Übler et al. (2023b) Übler H., et al., 2023b, arXiv e-prints, arXiv.2312.03589

##  Appendix A CND relaxation for single SMBH simulation

Figure 16: Gas surface density maps during the refinement relaxation (panels 1
– 4) and the production run with the subgrid disc model (panels 4 – 5). This
illustrates the creation and dissipation of overdensities and spiral arms
during the refinement relaxation phase.

Here we present more detailed plots illustrating the evolution of the cavity
for the single SMBH simulations during the refinement relaxation and once the
accretion disc subgrid model is switched on. Note that we deliberately do not
let the cavity fully fill with gas before activating our model as we are
mainly interested in the low-accretion-rate regime.

Figure 16 shows gas density projections of the cavity during the refinement
relaxation (panels 1 – 4) and the production run with the subgrid disc model
(panels 4 – 5). After the refinement model is activated, several instabilities
develop, amplified by cooling from the CND onto the central region (also see
Bennett & Sijacki, 2020; Anglés-Alcázar et al., 2021). This leads to the
formation of spiral arms and a dense core which then propagates outwards (see
third panel). By the time that the subgrid model is activated (fourth panel),
this overdensity has dissipated. However, during the simulation, the cavity
continues to fill with gas leading to a mild increase in inflow rates.

Figure 17: Gas surface density, cell sizes, sound speed and azimuthal velocity
as a function of radius during the refinement relaxation phase (green lines)
and the subgrid model simulation phase (purple lines).

Figure 17 shows this more quantitatively, plotting the gas density, cell
sizes, sound speed and azimuthal velocity as a function of radius during the
refinement relaxation phase (green lines) and the subgrid model simulation
phase (purple lines).