# Optimal Power Flow in Four-Wire Distribution Networks: Formulation and
Benchmarking

Sander Claeys
Geert Deconinck  KU Leuven / EnergyVille
Leuven / Genk, Belgium
{sander.claeys<EMAIL_ADDRESS>Frederik Geth  CSIRO Energy
Newcastle, Australia
<EMAIL_ADDRESS>
###### Abstract

In recent years, several applications have been proposed in the context of
distribution networks. Many of these can be formulated as an optimal power
flow problem, a mathematical optimization program which includes a model of
the steady-state physics of the electricity network. If the network loading is
balanced and the lines are transposed, the network model can be simplified to
a single-phase equivalent model. However, these assumptions do not apply to
low-voltage distribution networks, so the network model should model the
effects of phase unbalance correctly. In many parts of the world, the low-
voltage distribution network has four conductors, i.e. three phases and a
neutral. This paper develops OPF formulations for such networks, including
transformers, shunts and voltage-dependent loads, in two variable spaces, i.e.
current-voltage and power-voltage, and compares them for robustness and
scalability. A case study across 128 low-voltage networks also quantifies the
modelling error introduced by Kron reductions and its impact on the solve
time. This work highlights the advantages of formulations in current-voltage
variables over power-voltage, for four-wire networks.

###### Index Terms:

nonlinear programming, power distribution, unbalanced optimal power flow

Submitted to the 22nd Power Systems Computation Conference (PSCC 2022).

##  I Introduction

In low-voltage (LV) distribution networks, single-phase loading and generation
can lead to unequal voltage magnitudes between phases. In networks with
neutral conductors (e.g. European-style four-wire networks), the neutral
current can lead to neutral-point shift, which can worsen the voltage
unbalance [1]. Barbato et al. [2] report on distribution network measurements
obtained through an extensive monitoring infrastructure. At all times, the
neutral current was above 10 % of the peak phase current, and at several times
it exceeded some of the phase currents in magnitude. These effects are largely
overlooked in the state-of-the-art on power flow optimization in distribution
networks. Therefore, this paper compares optimization models which are well-
suited for networks where modelling the neutral wire explicitly is paramount.

In the analysis of three-phase transmission networks, the network is often
modeled by its positive-sequence a.k.a. ‘single-phase’ approximation. In
distribution networks however, lines are not transposed and line currents are
unbalanced; the positive-sequence approximation is a poor substitute under
these circumstances [3]. The series impedance of a power line or cable is
modeled by a square matrix of self- and mutual impedances, for each and
between all conductors. A four-wire segment has a neutral conductor in
addition to three phase conductors. If the neutral is perfectly grounded
($0\,\Omega$) at both ends, a ‘Kron-reduction’ will transform the 4$\times$4
impedance matrix to an equivalent 3$\times$3 matrix. This simplifies the
analysis, and is often done even though no such grounding exists in reality
[4]. When not all nodes are actually grounded or have a nonzero grounding
resistance, a Kron-reduced model will under-estimate the voltage deviations
[5]. The Kron-reduced $3\times 3$ matrix finally is the basis for deriving the
sequence components. The whole process is summarized in Fig. 1. There are many
alternative metrics to capture unbalance levels; Girigoudar et al. discuss
these in-depth, and compare them analytically [6].

Figure 1: The approximation step applied to the series impedance. The Kron-
reduction reduces the original 4$\times$4 to a 3$\times$3 matrix. The single-
phase equivalent model represents the line with the complex number
${\color[rgb]{0,0,1}z}^{+}$.

###  I-A State-of-the-art of Unbalanced OPF

Ramos de Araujo et al. [7] develop an n-wire nonlinear UBOPF model in the
current-voltage variable space that is capable of representing neutral wires.
They use a NLP interior point method to solve the optimization problem.
Dall’Anese et al. propose a bus-injection model semidefinite relaxation (BIM-
SDP) of UBOPF [8], and Gan and Low introduced a branch-flow model semidefinite
relaxation (BFM-SDP) [9]. Furthermore, there is no well-developed UBOPF
benchmark library; instead, authors modify PF benchmarks ad-hoc.

Kron-reduced models assume that the return current of connected devices goes
into ground at the bus. In four-wire networks however, the return current is
injected in the neutral and - if applicable - the neutral grounding shunt. A
nonzero current in the neutral induces voltage in other phases. There are few
studies which use four-wire network data in the context of OPF problems. To
the best of the authors’ knowledge, Su et al. [10] were the first to do this.
They optimize the setpoints of PV inverters with respect to a weighted
combination of network losses, generation cost, voltage deviations and
unbalance. Recently, Gastalver-Rubio et al. [11] did a similar study. They
minimize network losses, which is also argued to be a proxy for the unbalance.
The UBOPF studies by Ramos de Araujo et al. [7] is arguably the most in-depth,
including analysis of neutral voltage rise.

In addition to exact formulations, few authors propose convex relaxations for
four-wire networks. Liu et al. extend previous work on BIM-SDP relaxations of
UBOPF to four-wire networks [12]. Additionally, they introduce a technique to
reduce the number of variables and to recover a feasible solution when the
relaxation is inexact. The solutions are not compared with solutions obtained
with an exact, non-linear formulation. Usman et al. also proposed a four-wire
extension to BIM-SDP relaxations, and further included ZIP loads [13]. The
case study examines five different test feeders, including IEEE13. However,
Usman et al. state that a wide range of solvers, including Ipopt, failed to
solve their implementation of the exact, non-linear formulation, and therefore
they cannot report results for the optimality gap. As part of the case
studies, they also examine relaxations of Kron-reduced network models, and
compare the objective and tightness with the corresponding four-wire model. In
later work, they add delta-connected loads [14] and extensions of BFM-SDP
relaxations [15]. Finally, Liu et al. [16] also apply four-wire models in the
context of state estimation, and compare rectangular and polar voltage
variables.

We conclude that there is lots of ambiguous information in the literature
surrounding the scalability of exact formulations of unbalanced OPF through
nonlinear programming methods. To make sure the results are sufficiently
generalizable, we must prioritize formulations which are expressive enough to
support the key features used in the IEEE unbalanced power flow test cases:

  * •

transformers: lossy, multi-conductor, multi-winding;

  * •

lines: multi-conductor with line-charging;

  * •

loads: multi-phase, ZIP voltage dependency;

  * •

shunts: multi-conductor admittance;

  * •

generators: multi-phase, PQ;

  * •

bounds: voltage unbalance factor, neutral shift, negative sequence voltage,
phase-to-phase voltage, phase-to-neutral voltage.

###  I-B Scope and contributions

Therefore this paper describes in detail two exact formulations to the UBOPF
problem, and explores computation speed and reliability. Special care was
taken to show how everything can be formulated as a _purely quadratic_
optimization problem. The first formulation is a quadratic one in rectangular
current and voltage variables (similar to [7]); the second is the equivalent
power-voltage rectangular form. These formulations include many modeling
features, validated against OpenDSS on several popular PF benchmarks with an
accuracy of at least 1E-7.

The main contributions of this paper are 1) a methodology to convert Kron-
reduced network data to four-wire data, 2) a set of numerical experiments
comparing balanced, Kron-reduced and four-wire OPF models across 128 LV test
feeders, in terms of constraint violation and scalability, and 3) the
observation that power balance constraints are a relaxation of Kirchhoff’s
current law, potentially leading to inexact solutions and numerical stability
issues four-wire OPF formulations using power balance constraints.

###  I-C Notation and preliminaries

The symbol for current is $I$; voltage $U$; power $S/P/Q$. Complex variables
get a blue typeface, real variables a red one; also, vectors and matrix
variables get a bold typeface. Superscript ‘r’ indicates the real part and ‘i’
the imaginary part of an otherwise complex-valued symbol.

Figure 2 illustrates a bus $i$ that connects to a line $l$, a single-phase
load $d$, a single-phase generator $g$ and a shunt $s$. We furthermore use $t$
as the transformer index. The buses $i$ have terminals
$k\in\mathcal{P}_{i}\subseteq\\{a,b,c,n\\}$. The component conductors are
indexed over $p\in\mathcal{U}\subseteq\\{x,y,z,w\\}$, e.g. for a single-phase
branch $l$ we have $\mathcal{U}_{l}=\\{x,w\\}$. The bus voltage vector
${\color[rgb]{0,0,1}\mathbf{U}}_{i}$ stacks the terminal voltages
${\color[rgb]{0,0,1}U}_{ip}$,
${\color[rgb]{0,0,1}\mathbf{U}}_{i}=[{\color[rgb]{0,0,1}U}_{ip}]_{p\in\mathcal{P}_{i}}$.
The branch current vector ${\color[rgb]{0,0,1}\mathbf{I}}_{lij}$ stacks the
conductor currents
${\color[rgb]{0,0,1}\mathbf{I}}_{lij}=[{\color[rgb]{0,0,1}I}_{lij,u}]_{u\in\mathcal{U}_{l}}$.
We therefore define an mapping for each component, from its conductors to (a
subset of) the bus terminals. For instance for the branch with
$\mathcal{U}_{l}=\\{x,w\\}$ we have the from-side conductor-terminal mapping
$\mathcal{X}_{lij}=\\{(x,b),(w,n)\\}\subset\mathcal{U}_{l}\times\mathcal{P}_{i}$.
Note that this data model allows for arbitrary connections without relying on
the order of the conductors.

Figure 2: Illustration of a four-wire line which connects to a single-phase
generator $g$, load $d$ and shunt $s$.

Where-ever more concise, we will make abstraction of delta or wye connections
for devices. Figure 3 shows that both can be interpreted as equivalent
compositions of single-phase components. For transformers, which are two-port
components, this can be done for each port separately, leading to Yy, Dd, Yd
and Dy transformers.

(a) wye-connected

(b) delta-connected

Figure 3: Three-phase components can be equivalently composed of three single-
phase components.

##  II Current-voltage Formulation

The derivation uses primarily complex current (I) and voltage (V) variables,
both split into rectangular (R) coordinates; hence, it is referred to as the
IVR formulation.

(a) multi-conductor line

(b) single-phase transformer

Figure 4: Two-port (excluding ground) components connect two buses $i$ and
$j$.

###  II-A Multi-conductor line

A multi-conductor line $l$ consists of $n$ conductors, which connect $n$
terminals at bus $i$ to $n$ terminals at bus $j$. A $\pi$-model of the line is
parameterized by three complex matrices,
${\color[rgb]{0,0,1}\mathbf{Z}}^{\text{s}}_{l}$,
${\color[rgb]{0,0,1}\mathbf{Y}}^{\text{fr}}_{l}$ and
${\color[rgb]{0,0,1}\mathbf{Y}}^{\text{to}}_{l}$. Figure 4a illustrates the
$\pi$-model, governed by

$\displaystyle{\color[rgb]{0,0,1}\mathbf{U}}_{i}-{\color[rgb]{0,0,1}\mathbf{U}}_{j}$
$\displaystyle={\color[rgb]{0,0,1}\mathbf{Z}}^{\text{s}}_{l}{\color[rgb]{0,0,1}\mathbf{I}}^{\text{s}}_{lij},$
(1a)
$\displaystyle{\color[rgb]{0,0,1}\mathbf{I}}_{lij}={\color[rgb]{0,0,1}\mathbf{Y}}^{\text{fr}}_{l}{\color[rgb]{0,0,1}\mathbf{U}}_{i}+{\color[rgb]{0,0,1}\mathbf{I}}^{\text{s}}_{lij}$
$\displaystyle,\quad{\color[rgb]{0,0,1}\mathbf{I}}_{lji}={\color[rgb]{0,0,1}\mathbf{Y}}^{\text{to}}_{l}{\color[rgb]{0,0,1}\mathbf{U}}_{j}-{\color[rgb]{0,0,1}\mathbf{I}}^{\text{s}}_{lij}.$
(1b)

###  II-B Single-phase transformer

Figure 4b illustrates an ideal, two-winding, single-phase transformer $t$,
with turns ratio ${\color[rgb]{1,0,0}n}_{t}$. Note that
${\color[rgb]{1,0,0}n}_{t}$ is a real and therefore does not model any phase
shifting. We observe the conductor-terminal maps
$\mathcal{X}_{lij}=\\{(x,a),(y,b)\\}$, $\mathcal{X}_{lji}=\\{(x,a),(y,b)\\}$.
The transformer’s feasible set is therefore,

$\displaystyle{\color[rgb]{0,0,1}I}_{tijx}+{\color[rgb]{0,0,1}I}_{tijy}$
$\displaystyle=0,$ (2a)
$\displaystyle{\color[rgb]{0,0,1}I}_{tjix}+{\color[rgb]{0,0,1}I}_{tjiy}$
$\displaystyle=0,$ (2b)
$\displaystyle{\color[rgb]{1,0,0}n}_{t}{\color[rgb]{0,0,1}I}_{tijx}+{\color[rgb]{0,0,1}I}_{tjix}$
$\displaystyle=0,$ (2c)
$\displaystyle{\color[rgb]{0,0,1}U}_{ia}-{\color[rgb]{0,0,1}U}_{ib}$
$\displaystyle={\color[rgb]{1,0,0}n}_{t}\left({\color[rgb]{0,0,1}U}_{ja}-{\color[rgb]{0,0,1}U}_{jb}\right).$
(2d)

A multi-winding, multi-phase, lossy transformer can be equivalently
represented by a composition of lines, shunts and ideal two-winding
transformers [17]. Combining these elementary building blocks, this
formulation can represent all transformers covered in [17].

(a) load

(b) generator

(c) shunt

Figure 5: One-port components connect to a single bus.

###  II-C Load

####  II-C1 Three-phase

A three-phase, wye-connected load connects to bus $i$ with the conductor-
terminal map $\\{(x,a),(y,b),(z,c),(w,n)\\}$. The load current and power
variables therefore satisfy,

$\displaystyle 0$
$\displaystyle={\color[rgb]{0,0,1}I}_{dx}+{\color[rgb]{0,0,1}I}_{dy}+{\color[rgb]{0,0,1}I}_{dz}+{\color[rgb]{0,0,1}I}_{dw},$
(3a) $\displaystyle{\color[rgb]{0,0,1}S}_{dx}$
$\displaystyle=({\color[rgb]{0,0,1}U}_{ia}-{\color[rgb]{0,0,1}U}_{in})\left({\color[rgb]{0,0,1}I}_{dx}\right)^{*},$
(3b) $\displaystyle{\color[rgb]{0,0,1}S}_{dy}$
$\displaystyle=({\color[rgb]{0,0,1}U}_{ib}-{\color[rgb]{0,0,1}U}_{in})\left({\color[rgb]{0,0,1}I}_{dy}\right)^{*},$
(3c) $\displaystyle{\color[rgb]{0,0,1}S}_{dz}$
$\displaystyle=({\color[rgb]{0,0,1}U}_{ic}-{\color[rgb]{0,0,1}U}_{in})\left({\color[rgb]{0,0,1}I}_{dz}\right)^{*},$
(3d)

where e.g. ${\color[rgb]{0,0,1}S}_{dx}$ is the power drawn between terminals
$a$ and $n$. The power flow for conductor $x$ is defined as
${\color[rgb]{0,0,1}S}^{\text{b}}_{dx}={\color[rgb]{0,0,1}U}_{ia}\left({\color[rgb]{0,0,1}I}_{dx}\right)^{*}$,
and similarly for the others. To define voltage-sensitive loads, we consider
three-phase loads as the composition of single-phase loads.

####  II-C2 Single-phase

Figure 5a illustrates a single-phase load $d$, which draws an amount of power
${\color[rgb]{0,0,1}S}_{d}$ between the terminals $a$ and $b$ at bus $i$, with
map $\mathcal{X}_{di}=\\{(x,a),(y,b)\\}$,

$\displaystyle{\color[rgb]{0,0,1}I}_{dx}+{\color[rgb]{0,0,1}I}_{dy}$
$\displaystyle=0,$ (4a) $\displaystyle{\color[rgb]{0,0,1}S}_{d}$
$\displaystyle=\left({\color[rgb]{0,0,1}U}_{ia}-{\color[rgb]{0,0,1}U}_{ib}\right)\left({\color[rgb]{0,0,1}I}_{dx}\right)^{*}.$
(4b)

Loads are further characterized by how the drawn power depends on the voltage
magnitude applied to them,
${\color[rgb]{0,0,1}S}_{d}=f(\left|{\color[rgb]{0,0,1}U}_{ia}-{\color[rgb]{0,0,1}U}_{ib}\right|)$.
Below, several common models are discussed. We note that ZIP loads can be
represented as a composition of a constant impedance, constant current and a
constant power load.

##### Constant power

The power is constant and independent of the applied voltage.

$\displaystyle{\color[rgb]{1,0,0}P}_{d}$
$\displaystyle={\color[rgb]{1,0,0}P}^{\text{nom}}_{d},$
$\displaystyle{\color[rgb]{1,0,0}Q}_{d}$
$\displaystyle={\color[rgb]{1,0,0}Q}^{\text{nom}}_{d}.$ (5)

##### Constant impedance

The power is proportional to the square of the voltage magnitude,

$\displaystyle{\color[rgb]{1,0,0}P}_{d}$
$\displaystyle=\frac{{\color[rgb]{1,0,0}P}^{\text{nom}}_{d}}{\left({\color[rgb]{1,0,0}U}^{\text{nom}}_{d}\right)^{2}}\left|{\color[rgb]{0,0,1}U}_{iab}\right|^{2},$
$\displaystyle{\color[rgb]{1,0,0}Q}_{d}$
$\displaystyle=\frac{{\color[rgb]{1,0,0}Q}^{\text{nom}}_{d}}{\left({\color[rgb]{1,0,0}U}^{\text{nom}}_{d}\right)^{2}}\left|{\color[rgb]{0,0,1}U}_{iab}\right|^{2}.$
(6)

where we introduced the shorthand
${\color[rgb]{0,0,1}U}_{iab}={\color[rgb]{0,0,1}U}_{ia}-{\color[rgb]{0,0,1}U}_{ib}$.
This seems like it is quadratic, but can be represented linearly in the
variables ${\color[rgb]{0,0,1}U}_{iab}$ and ${\color[rgb]{0,0,1}I}_{dx}$. If
${\color[rgb]{0,0,1}U}_{iab}\neq 0$,

$\displaystyle\Leftrightarrow\;$
$\displaystyle{\color[rgb]{0,0,1}I}_{dx}=\left({\color[rgb]{1,0,0}g}_{d}+j{\color[rgb]{1,0,0}b}_{d}\right){\color[rgb]{0,0,1}U}_{iab},$
(7a) $\displaystyle\hskip 10.00002pt{\color[rgb]{1,0,0}g}_{d}=\hskip
6.99997pt{\color[rgb]{1,0,0}P}^{\text{nom}}_{d}$
$\displaystyle/\left({\color[rgb]{1,0,0}U}^{\text{nom}}_{d}\right)^{2},\hskip
10.00002pt{\color[rgb]{1,0,0}b}_{d}=-{\color[rgb]{1,0,0}Q}^{\text{nom}}_{d}/\left({\color[rgb]{1,0,0}U}^{\text{nom}}_{d}\right)^{2}.$
(7b)

##### Constant current

The power is proportional to the voltage magnitude,

$\displaystyle{\color[rgb]{1,0,0}P}_{d}$
$\displaystyle=\frac{{\color[rgb]{1,0,0}P}^{\text{nom}}_{d}}{{\color[rgb]{1,0,0}U}^{\text{nom}}_{d}}\left|{\color[rgb]{0,0,1}U}_{iab}\right|,$
$\displaystyle{\color[rgb]{1,0,0}Q}_{d}$
$\displaystyle=\frac{{\color[rgb]{1,0,0}Q}^{\text{nom}}_{d}}{{\color[rgb]{1,0,0}U}^{\text{nom}}_{d}}\left|{\color[rgb]{0,0,1}U}_{iab}\right|.$
(8)

This can be represented using an auxiliary variable
${\color[rgb]{1,0,0}U}^{\text{m,sqr}}_{iab}$ defined as,

$\displaystyle{\color[rgb]{1,0,0}U}^{\text{m,sqr}}_{iab}=\left|{\color[rgb]{0,0,1}U}_{ia}-{\color[rgb]{0,0,1}U}_{ib}\right|^{2}=\left({\color[rgb]{1,0,0}U}^{\text{r}}_{ia}-{\color[rgb]{1,0,0}U}^{\text{r}}_{ib}\right)^{2}+\left({\color[rgb]{1,0,0}U}^{\text{i}}_{ia}-{\color[rgb]{1,0,0}U}^{\text{i}}_{ib}\right)^{2},$
(9)

and two quadratic equalities and two linear inequalities,

$\displaystyle\left({\color[rgb]{1,0,0}P}_{d}\right)^{2}$
$\displaystyle=\left(\frac{{\color[rgb]{1,0,0}P}^{\text{nom}}_{d}}{{\color[rgb]{1,0,0}U}^{\text{nom}}_{d}}\right)^{2}{\color[rgb]{1,0,0}U}^{\text{m,sqr}}_{iab},$
$\displaystyle\text{sign}\left({\color[rgb]{1,0,0}P}^{\text{nom}}_{d}\right){\color[rgb]{1,0,0}P}_{d}\geq
0,$ (10a) $\displaystyle\left({\color[rgb]{1,0,0}Q}_{d}\right)^{2}$
$\displaystyle=\left(\frac{{\color[rgb]{1,0,0}Q}^{\text{nom}}_{d}}{{\color[rgb]{1,0,0}U}^{\text{nom}}_{d}}\right)^{2}{\color[rgb]{1,0,0}U}^{\text{m,sqr}}_{iab},$
$\displaystyle\text{sign}\left({\color[rgb]{1,0,0}Q}^{\text{nom}}_{d}\right){\color[rgb]{1,0,0}Q}_{d}\geq
0.$ (10b)

###  II-D Single-phase generator

Figure 5b shows a single-phase generator $g$, with map
$\mathcal{X}_{gi}=\\{(x,a),(y,b)\\}$,

$\displaystyle{\color[rgb]{0,0,1}I}_{gx}+{\color[rgb]{0,0,1}I}_{gy}$
$\displaystyle=0,$ (11a) $\displaystyle{\color[rgb]{0,0,1}S}_{g}$
$\displaystyle=\left({\color[rgb]{0,0,1}U}_{ia}-{\color[rgb]{0,0,1}U}_{ib}\right)\left({\color[rgb]{0,0,1}I}_{gx}\right)^{*}.$
(11b)

The power
${\color[rgb]{0,0,1}S}_{g}={\color[rgb]{1,0,0}P}_{g}+j{\color[rgb]{1,0,0}Q}_{g}$
it injects, is bounded by

$\displaystyle{\color[rgb]{1,0,0}P}^{\text{min}}_{g}\leq{\color[rgb]{1,0,0}P}_{g}\leq{\color[rgb]{1,0,0}P}^{\text{max}}_{g},\quad$
$\displaystyle{\color[rgb]{1,0,0}Q}^{\text{min}}_{g}\leq{\color[rgb]{1,0,0}Q}_{g}\leq{\color[rgb]{1,0,0}Q}^{\text{max}}_{g},$
(12a)
$\displaystyle\left({\color[rgb]{1,0,0}P}_{g}\right)^{2}+\left({\color[rgb]{1,0,0}Q}_{g}\right)^{2}$
$\displaystyle\leq\left({\color[rgb]{1,0,0}S}^{\text{max}}_{g}\right)^{2}.$
(12b)

###  II-E Shunt

Figure 5c shows a multi-conductor shunt $s$, defined by,

$\displaystyle{\color[rgb]{0,0,1}\mathbf{I}}_{s}$
$\displaystyle={\color[rgb]{0,0,1}\mathbf{Y}}_{s}{\color[rgb]{0,0,1}\mathbf{U}}_{i}.$
(13)

Lossy grounding of the neutral conductor can be modelled as a single-conductor
shunt.

###  II-F Kirchhoff’s current law at buses

Let $\mathcal{L}_{ik}$ denote the set of all line conductors $lijp$ connected
to terminal $(i,k)$ at bus $i$; the same goes for $\mathcal{T}_{ik}$
(transformers), $\mathcal{D}_{ik}$ (loads), $\mathcal{S}_{ik}$ (shunts) and
$\mathcal{G}_{ik}$ (generators). The expression for Kirchhoff’s current law,
requiring that the sum of all incoming currents equals the sum of all outgoing
currents at each terminal in the network, at each terminal $(i,k)$,

$\displaystyle\sum_{gp\in\mathcal{G}_{ik}}\\!\\!{\color[rgb]{0,0,1}I}_{gc}=\\!\\!\\!\sum_{lijp\in\mathcal{N}_{ik}}\\!\\!{\color[rgb]{0,0,1}I}_{lij,c}+\\!\\!\sum_{tijp\in\mathcal{T}_{ik}}\\!\\!{\color[rgb]{0,0,1}I}_{tij,c}+\\!\\!\sum_{dp\in\mathcal{D}_{ik}}\\!\\!{\color[rgb]{0,0,1}I}_{dc}+\\!\\!\sum_{sp\in\mathcal{S}_{ik}}\\!\\!{\color[rgb]{0,0,1}I}_{sc}.$
(14)

For the other terminals at the bus, we define similar sets and expressions.
Figure 2 illustrates this with an example, where a bus connects to a line $l$,
a single-phase load $d$, a single-phase generator $g$ and a shunt $s$.
Applying (14) at each of the four terminals,

$\displaystyle{\color[rgb]{0,0,1}I}_{gx}$
$\displaystyle={\color[rgb]{0,0,1}I}_{lijx},$ $\displaystyle 0$
$\displaystyle={\color[rgb]{0,0,1}I}_{lijy}+{\color[rgb]{0,0,1}I}_{dx},$ (15a)
$\displaystyle 0$ $\displaystyle={\color[rgb]{0,0,1}I}_{lijz},$
$\displaystyle{\color[rgb]{0,0,1}I}_{gy}$
$\displaystyle={\color[rgb]{0,0,1}I}_{lijw}+{\color[rgb]{0,0,1}I}_{dy}+{\color[rgb]{0,0,1}I}_{sx}.$
(15b)

###  II-G Technical envelopes

##### Current ratings

The magnitude of each of the current variables
${\color[rgb]{0,0,1}I}={\color[rgb]{1,0,0}I}^{\text{r}}+j{\color[rgb]{1,0,0}I}^{\text{i}}$
introduced before, can be enforced with a second-order-cone (SOC) constraint,

$\displaystyle({\color[rgb]{1,0,0}I}^{\text{r}})^{2}+({\color[rgb]{1,0,0}I}^{\text{i}})^{2}\leq({\color[rgb]{1,0,0}I}^{\text{max}})^{2}.$
(16)

##### Voltage magnitude bounds

Some voltage magnitude bounds apply to the difference of two voltage phasors
(phase-to-phase
${\color[rgb]{1,0,0}U}^{\text{min}}_{iab},{\color[rgb]{1,0,0}U}^{\text{max}}_{iab}$
and phase-to-neutral bounds
${\color[rgb]{1,0,0}U}^{\text{min}}_{ian},{\color[rgb]{1,0,0}U}^{\text{max}}_{ian}$),
whilst others apply directly to a single voltage phasor (neutral shift bound,
${\color[rgb]{1,0,0}U}^{\text{min}}_{in},{\color[rgb]{1,0,0}U}^{\text{max}}_{in}$).
Both can be expressed as quadratic inequalities,

$\displaystyle({\color[rgb]{1,0,0}U}^{\text{min}}_{ia})^{2}$
$\displaystyle\leq$ $\displaystyle({\color[rgb]{1,0,0}U}^{\text{r}}_{ia})^{2}$
$\displaystyle+$ $\displaystyle({\color[rgb]{1,0,0}U}^{\text{i}}_{ia})^{2}$
$\displaystyle\leq$
$\displaystyle({\color[rgb]{1,0,0}U}^{\text{max}}_{ia})^{2},$ (17)
$\displaystyle({\color[rgb]{1,0,0}U}^{\text{min}}_{iab})^{2}$
$\displaystyle\leq$
$\displaystyle({\color[rgb]{1,0,0}U}^{\text{r}}_{iab})^{2}$ $\displaystyle+$
$\displaystyle({\color[rgb]{1,0,0}U}^{\text{i}}_{iab})^{2}$
$\displaystyle\leq$
$\displaystyle({\color[rgb]{1,0,0}U}^{\text{max}}_{iab})^{2},$ (18)

where the upper bound is representable as a SOC constraint, whilst the lower
bound is non-convex.

##### Unbalance bounds

Take a bus $i$, with phase terminals $a$, $b$ and $c$, and neutral $n$. We
assume $b$ lags $a$, $c$ lags $b$. The positive sequence
${\color[rgb]{0,0,1}U}^{\text{p}}_{i}$, negative sequence
${\color[rgb]{0,0,1}U}^{\text{n}}_{i}$ and zero sequence
${\color[rgb]{0,0,1}U}^{\text{z}}_{i}$ phasors are then defined,

$\displaystyle{\color[rgb]{0,0,1}U}^{\text{z}}_{i}={\color[rgb]{1,0,0}U}^{\text{z,r}}_{i}+j{\color[rgb]{1,0,0}U}^{\text{z,i}}_{i}=$
$\displaystyle\tfrac{1}{3}$ $\displaystyle{\color[rgb]{0,0,1}U}_{ian}$
$\displaystyle+$ $\displaystyle\tfrac{1}{3}$
$\displaystyle{\color[rgb]{0,0,1}U}_{ibn}$ $\displaystyle+$
$\displaystyle\tfrac{1}{3}$ $\displaystyle{\color[rgb]{0,0,1}U}_{icn},$ (19a)
$\displaystyle{\color[rgb]{0,0,1}U}^{\text{p}}_{i}={\color[rgb]{1,0,0}U}^{\text{p,r}}_{i}+j{\color[rgb]{1,0,0}U}^{\text{p,i}}_{i}=$
$\displaystyle\tfrac{1}{3}$ $\displaystyle{\color[rgb]{0,0,1}U}_{ian}$
$\displaystyle+$ $\displaystyle\tfrac{\alpha}{3}$
$\displaystyle{\color[rgb]{0,0,1}U}_{ibn}$ $\displaystyle+$
$\displaystyle\tfrac{\alpha^{2}}{3}$
$\displaystyle{\color[rgb]{0,0,1}U}_{icn},$ (19b)
$\displaystyle{\color[rgb]{0,0,1}U}^{\text{n}}_{i}={\color[rgb]{1,0,0}U}^{\text{n,r}}_{i}+j{\color[rgb]{1,0,0}U}^{\text{n,i}}_{i}=$
$\displaystyle\tfrac{1}{3}$ $\displaystyle{\color[rgb]{0,0,1}U}_{ian}$
$\displaystyle+$ $\displaystyle\tfrac{\alpha^{2}}{3}$
$\displaystyle{\color[rgb]{0,0,1}U}_{ibn}$ $\displaystyle+$
$\displaystyle\tfrac{\alpha}{3}$ $\displaystyle{\color[rgb]{0,0,1}U}_{icn},$
(19c)

where $\alpha=\angle 120^{\circ}=-0.5+j\tfrac{\sqrt{3}}{2}$. The voltage
unbalance factor,
$\text{VUF}=|{\color[rgb]{0,0,1}U}^{\text{n}}_{i}|/|{\color[rgb]{0,0,1}U}^{\text{p}}_{i}|$,
is a common measure for the amount of ‘unbalance’. Some grid codes specify
$\text{VUF}_{i}\leq\text{VUF}^{\text{\,max}}_{\,i}$,

$\displaystyle({\color[rgb]{1,0,0}U}^{\text{n,r}}_{i})^{2}+({\color[rgb]{1,0,0}U}^{\text{n,i}}_{i})^{2}$
$\displaystyle\leq(\text{VUF}^{\text{\,max}}_{\,i})^{2}\left(({\color[rgb]{1,0,0}U}^{\text{p,r}}_{i})^{2}+({\color[rgb]{1,0,0}U}^{\text{p,i}}_{i})^{2}\right),$
(20)

whilst others limit the negative sequence,
$|{\color[rgb]{0,0,1}U}^{\text{n}}_{i}|\leq{\color[rgb]{1,0,0}U}^{\text{n,max}}_{i}$,

$\displaystyle({\color[rgb]{1,0,0}U}^{\text{n,r}}_{i})^{2}+({\color[rgb]{1,0,0}U}^{\text{n,i}}_{i})^{2}$
$\displaystyle\leq({\color[rgb]{1,0,0}U}^{\text{n,max}}_{i})^{2}.$ (21)

##### Power bounds

Where desired, bounds on the product of any ${\color[rgb]{0,0,1}U}$ and
${\color[rgb]{0,0,1}I}$ can be imposed quadratically by introducing auxiliary
variables ${\color[rgb]{1,0,0}P}$ and ${\color[rgb]{1,0,0}Q}$, and imposing
the constraints

$\displaystyle{\color[rgb]{1,0,0}P}$
$\displaystyle={\color[rgb]{1,0,0}U}^{\text{r}}{\color[rgb]{1,0,0}I}^{\text{r}}+{\color[rgb]{1,0,0}U}^{\text{i}}{\color[rgb]{1,0,0}I}^{\text{i}},$
$\displaystyle{\color[rgb]{1,0,0}Q}$
$\displaystyle=-{\color[rgb]{1,0,0}U}^{\text{r}}{\color[rgb]{1,0,0}I}^{\text{i}}+{\color[rgb]{1,0,0}U}^{\text{i}}{\color[rgb]{1,0,0}I}^{\text{r}}.$
(22)

Analogous to (12), bounds ${\color[rgb]{1,0,0}P}^{\text{min}}$,
${\color[rgb]{1,0,0}P}^{\text{max}}$, ${\color[rgb]{1,0,0}Q}^{\text{min}}$,
${\color[rgb]{1,0,0}Q}^{\text{max}}$ and ${\color[rgb]{1,0,0}S}^{\text{max}}$
can then be imposed.

##  III Power-Voltage Formulation

In this section we present derivations that primarily use _complex power_ and
voltage variables, again both split into rectangular (R) coordinates; hence,
it is referred to as the ACR form. Because IVR and ACR share rectangular
voltage variables, all expressions that just depend on voltage can be used as-
is. For brevity, we don’t repeat such constraints in this section, but list
the shared constraints later when defining the feasible sets (§VI-A).

###  III-A Kirchhoff’s current law at buses

We reformulate Kirchhoff’s current law, $\sum{\color[rgb]{0,0,1}I}=0$, by
multiplying its conjugate on both sides by the nodal voltage
${\color[rgb]{0,0,1}U}$, obtaining
$\sum{\color[rgb]{0,0,1}U}{\color[rgb]{0,0,1}I}^{*}=0$. Introducing variables
${\color[rgb]{0,0,1}S}={\color[rgb]{0,0,1}U}({\color[rgb]{0,0,1}I})^{*}$ and
substituting these, then leads to $\sum{\color[rgb]{0,0,1}S}=0$, the power
balance form. The ‘flow variable’ interfacing the components, now is power
instead of current:

$\displaystyle\sum_{gp\in\mathcal{G}_{ik}}\\!\\!{\color[rgb]{0,0,1}S}_{gc}=\\!\\!\\!\sum_{lijp\in\mathcal{N}_{ik}}\\!\\!{\color[rgb]{0,0,1}S}_{lij,c}+\\!\\!\sum_{tijp\in\mathcal{T}_{ik}}\\!\\!{\color[rgb]{0,0,1}S}_{tij,c}+\\!\\!\sum_{dp\in\mathcal{D}_{ik}}\\!\\!{\color[rgb]{0,0,1}S}_{dc}+\\!\\!\sum_{sp\in\mathcal{S}_{ik}}\\!\\!{\color[rgb]{0,0,1}S}_{sc}.$
(23)

###  III-B Multi-conductor line

Lines and shunts allow a compact representation in terms of power flow
variables. The from and to-side line flow variables,
${\color[rgb]{0,0,1}\mathbf{S}}_{lij}$ and
${\color[rgb]{0,0,1}\mathbf{S}}_{lji}$, are linked to the voltage,

$\displaystyle{\color[rgb]{0,0,1}\mathbf{S}}_{lij}$
$\displaystyle={\color[rgb]{0,0,1}\mathbf{U}}_{i}\odot\left({\color[rgb]{0,0,1}\mathbf{Y}}^{\text{fr}}_{l}{\color[rgb]{0,0,1}\mathbf{U}}_{i}+{\color[rgb]{0,0,1}\mathbf{Y}}^{\text{s}}_{l}({\color[rgb]{0,0,1}\mathbf{U}}_{i}-{\color[rgb]{0,0,1}\mathbf{U}}_{j})\right)^{*},$
(24) $\displaystyle{\color[rgb]{0,0,1}\mathbf{S}}_{lji}$
$\displaystyle={\color[rgb]{0,0,1}\mathbf{U}}_{j}\odot\left({\color[rgb]{0,0,1}\mathbf{Y}}^{\text{to}}_{l}{\color[rgb]{0,0,1}\mathbf{U}}_{j}+{\color[rgb]{0,0,1}\mathbf{Y}}^{\text{s}}_{l}({\color[rgb]{0,0,1}\mathbf{U}}_{j}-{\color[rgb]{0,0,1}\mathbf{U}}_{i})\right)^{*},$
(25)

where
${\color[rgb]{0,0,1}\mathbf{Y}}^{\text{s}}_{l}=({\color[rgb]{0,0,1}\mathbf{Z}}^{\text{s}}_{l})^{-1}$,
and $\odot$ denotes the element-wise product.

###  III-C Shunt

The shunt power flow variable, ${\color[rgb]{0,0,1}\mathbf{S}}_{s}$, is linked
to the voltage by

$\displaystyle{\color[rgb]{0,0,1}\mathbf{S}}_{s}={\color[rgb]{0,0,1}\mathbf{U}}_{i}\odot({\color[rgb]{0,0,1}\mathbf{Y}}_{s}{\color[rgb]{0,0,1}\mathbf{U}}_{i})^{*}.$
(26)

###  III-D Loads, generators and transformers

We redefine the load model with power flow variables by reformulating (3) with
${\color[rgb]{0,0,1}S}^{\text{b}}_{dx}={\color[rgb]{0,0,1}U}_{ia}\left({\color[rgb]{0,0,1}I}_{dx}\right)^{*}$
etc.,

$\displaystyle 0$ $\displaystyle=\begin{aligned}
&{\color[rgb]{0,0,1}U}_{ib}{\color[rgb]{0,0,1}U}_{ic}{\color[rgb]{0,0,1}U}_{in}{\color[rgb]{0,0,1}S}^{\text{b}}_{dx}+{\color[rgb]{0,0,1}U}_{ia}{\color[rgb]{0,0,1}U}_{ic}{\color[rgb]{0,0,1}U}_{in}{\color[rgb]{0,0,1}S}^{\text{b}}_{dy}\\\
+\;&{\color[rgb]{0,0,1}U}_{ia}{\color[rgb]{0,0,1}U}_{ib}{\color[rgb]{0,0,1}U}_{in}{\color[rgb]{0,0,1}S}^{\text{b}}_{dz}+{\color[rgb]{0,0,1}U}_{ia}{\color[rgb]{0,0,1}U}_{ib}{\color[rgb]{0,0,1}U}_{ic}{\color[rgb]{0,0,1}S}^{\text{b}}_{dw},\end{aligned}$
(27a) $\displaystyle{\color[rgb]{0,0,1}S}_{dx}{\color[rgb]{0,0,1}U}_{ia}$
$\displaystyle={\color[rgb]{0,0,1}S}^{\text{b}}_{dx}({\color[rgb]{0,0,1}U}_{ia}-{\color[rgb]{0,0,1}U}_{in}),$
(27b) $\displaystyle{\color[rgb]{0,0,1}S}_{dy}{\color[rgb]{0,0,1}U}_{ia}$
$\displaystyle={\color[rgb]{0,0,1}S}^{\text{b}}_{dy}({\color[rgb]{0,0,1}U}_{ia}-{\color[rgb]{0,0,1}U}_{in}),$
(27c) $\displaystyle{\color[rgb]{0,0,1}S}_{dz}{\color[rgb]{0,0,1}U}_{ia}$
$\displaystyle={\color[rgb]{0,0,1}S}^{\text{b}}_{dz}({\color[rgb]{0,0,1}U}_{ia}-{\color[rgb]{0,0,1}U}_{in}).$
(27d)

Note that (27a) requires the introduction of several auxiliary variables to
reformulate it as a set of quadratic constraints (i.e. repeated substitution
$z\leftarrow xy$). This approach, despite including a fourth degree
polynomial, is popular for Kron-reduced forms, because when
${\color[rgb]{0,0,1}U}_{in}=0$ and terminal $n$ is grounded, (27) reduces to
${\color[rgb]{0,0,1}S}^{\text{b}}_{dx}={\color[rgb]{0,0,1}S}_{dx}$,
${\color[rgb]{0,0,1}S}^{\text{b}}_{dy}={\color[rgb]{0,0,1}S}_{dy}$ and
${\color[rgb]{0,0,1}S}^{\text{b}}_{dz}={\color[rgb]{0,0,1}S}_{dz}$.

Generators and transformers are similarly reformulated in terms of power flow
variables.

###  III-E Technical envelopes

##### Current ratings

The magnitude of each of the current variables
${\color[rgb]{0,0,1}I}={\color[rgb]{1,0,0}I}^{\text{r}}+j{\color[rgb]{1,0,0}I}^{\text{i}}$
introduced before, can be constrained with a second-order-cone (SOC)
constraint in the SU variable space,

$\displaystyle({\color[rgb]{1,0,0}P})^{2}+({\color[rgb]{1,0,0}Q})^{2}\leq({\color[rgb]{1,0,0}I}^{\text{max}})^{2}\left(({\color[rgb]{1,0,0}U}^{\text{r}})^{2}+({\color[rgb]{1,0,0}U}^{\text{i}})^{2}\right).$
(28)

##  IV Assembling four-wire OPF data sets

To obtain a realistic amount of unbalance, we use the dataset released by
Rigoni et al., i.e. a set of 128 LV, four-wire distribution networks,
representative of feeders in the UK [18]. However, only the positive and zero
sequence components of the lines are included, whilst for the upcoming
numerical experiments, full 4$\times$4 impedance data is paramount. We now
present a methodology to obtain the missing 4$\times$4 data, and to derive OPF
problems with interesting dispatchable resources and envelopes.

###  IV-A Obtaining four-wire data

The original dataset only includes positive and zero sequence components.
Figure 1 shows how these are calculated based on the 4$\times$4 impedance
matrices.

Carson’s equations provide a method to calculate the 4$\times$4 matrices,
based on the spacing and properties of the conductors and the ground [3]. A
linecode specifies the series impedance of a line per unit of length. The
following procedure assigns to each line an alternative linecode, which
contains four-wire impedance data:

  1. 1.

Compile a library of candidate linecodes $c\in\mathcal{C}$, of which the
4$\times$4 impedance matrices are known.

  2. 2.

For each of these, calculate
${\color[rgb]{0,0,1}z}^{+}_{c}={\color[rgb]{1,0,0}r}^{+}_{c}+j{\color[rgb]{1,0,0}x}^{+}_{c}$.

  3. 3.

For each linecode $o$ in the original dataset, replace it with the linecode
that is closest to it in
(${\color[rgb]{1,0,0}r}^{+}$,${\color[rgb]{1,0,0}x}^{+}$) space,

$\displaystyle
o\leftarrow\underset{c\in\mathcal{C}}{\operatorname{argmin}}\begin{vmatrix}r^{+}_{o}-r^{+}_{c}\\\
x^{+}_{o}-x^{+}_{c}\end{vmatrix}_{2}.$ (29)

Urquhart et al. studied in depth the assumptions and approximations typically
made in simulating low-voltage distribution networks. As part of this study,
the full impedance data for several underground linecodes common in the UK,
were calculated [19]. Linecodes corresponding to conductors with cross section
areas of $16$ (x2), $25$, $35$, $70$, $95$, $120$, $185$ and $300\,mm^{2}$,
are used to populate the library of candidate linecodes. Figure 6 shows
graphically how each linecode in the original dataset was mapped to one of
these 9 candidate linecodes.

Figure 6: Each original linecode (small markers) is replaced by a candidate
linecode (big markers). Note that the x-axis has a logarithmic scale, which
makes the clustering seem a bit counter-intuitive in this figure. One original
linecode outlier, ($r^{+}_{o}=0.162$, $x^{+}_{o}=0.050$), is not displayed in
this plot.

Once assigned, the candidate linecodes are used to re-derive the Kron-reduced
and single-phase equivalent models, as shown in Figure 7. The neutral in the
four-wire models, is grounded only at the substation (lossless). Also, the
original data contains sections with many short lines series; this is a
consequence of tracing the network topology from GIS data. When these short
lines have zero shunt admittance, and no bounds are imposed on the lines or
buses connecting them, the network can be simplified. Segments of short lines
in series can be replaced by a single long line, leading to fewer buses but
identical results. The network sizes reported in this paper are the ones after
aforementioned reduction.

###  IV-B Setting up interesting OPF cases

The data set released by Rigoni et al. [18] contains 100 time series for
various devices (households, electric vehicles, etc.), with a horizon of 24 h
and a resolution of 5 min. We apply the load and photovoltaic time series with
a power factor of 0.95. This case study will compare dispatch results for a
single snapshot. We choose timestep $t=207$ such that the phase-to-neutral
voltage magnitude is within the upper bound. The upstream grid is modelled as
a voltage source: the voltage phasors are fixed at the reference bus, and an
unbounded generator models the import/export.

Next, for each network separately, we add a set of dispatchable, distributed
generators (DGs) which are allocated and sized based on the following
procedure.

  * •

Add a single-phase DG at every 4th load. Each DG has a generation cost which
is ten times smaller than the generation cost of the reference bus generator,
modelling the grid import. Intuitively, this will lead the DG output to
increase until some operational bound becomes active.

  * •

Each DG $g$ has the same capacity ${\color[rgb]{1,0,0}P}^{\text{max}}_{g}$
(and furthermore
${\color[rgb]{1,0,0}P}^{\text{min}}_{g}={\color[rgb]{1,0,0}Q}^{\text{min}}_{g}=0$,
${\color[rgb]{1,0,0}Q}^{\text{max}}_{g}={\color[rgb]{1,0,0}P}^{\text{max}}_{g}$),
which is sized such that under the balanced decision model, the DGs can
congest the network.

The objective is to minimize the total DG dispatch cost. The main operational
constraints are the generator bounds, and upper and lower bounds on the phase-
to-neutral voltage magnitude at buses which have loads or generators
connected.

##  V Implementation and Validation

The formulations presented focus on simplicity of the derivation rather than
efficiency in terms of number of optimization variables. Therefore this
subsection discusses implementation details that were left out before.

##### Complex to real

Most optimization modelling languages do not support complex parameters and
variables, and require the user to enter the problem in real variables. Most
expressions are straight-forward generalisations of those presented in [20].

##### Perfect grounding

When a terminal $(i,a)$ is perfectly grounded, the associated voltage
variables can be replaced by constants,
${\color[rgb]{1,0,0}U}^{\text{i}}_{ia}\leftarrow 0$,
${\color[rgb]{1,0,0}U}^{\text{r}}_{ia}\leftarrow 0$. Instead of introducing a
variable for the grounding current, which would act as a slack variable in
KCL, one can drop the KCL constraint for $(i,a)$ and recover the grounding
current in post-processing.

##### Unnecessary conductor or terminal variables

We do not assume the network is 4-wire throughout and/or all buses have four 4
terminals. Instead of eliminating those (free) variables, e.g. by fixing them
to 0 numerically, we chose not to define the variables to begin with. Because
of this, extensions to n-conductor forms and/or multiple neutrals should be
straight-forward and numerically efficient.

##### Efficient three-phase composition

A three-phase, wye-connected load only requires 4 current variables. However,
when it is composed of three equivalent single-phase loads, 6 current
variables are used. In the implementation, we did the proper substitutions to
eliminate these extra 2 current variables, and more generally as well for
other multi-phase components.

##### Solver

All OPF instances were solved with the primal-dual interior point solver Ipopt
[21] (version 3.13.2, with MUMPS) using the automatic differentiation as
implemented in JuMP/MOI [22].

##### Initialization

Proper initialization of the voltage variables was crucial to achieve reliable
convergence with Ipopt. The voltage variables were initialized at the value
they would take under no-load conditions (and ignoring any shunts in the
system), i.e. we assume a balanced phasor at the reference bus and propagate
that across transformers to include the vector group angle shift, while
keeping track of the terminal-conductor maps.

##### Validation

The proposed formulation can also be used to solve PF problems as feasibility
problems. This approach was used to validate the optimization models against
the simulation models included in OpenDSS. Table I shows that across several
well-known benchmarks for unbalanced PF, the largest relative deviation for
any voltage phasor is 1.2E-8. These results are obtained with Ipopt (version
3.13.2, tolerance set to 1E-10), and compared against OpenDSS with tolerance
set to 1E-10.

TABLE I: The max. relative error on the voltage does not exceed 1.2E-8.

| IEEE13 | IEEE34 | IEEE123∗ | LVTestCase (t=1000)
---|---|---|---|---
$\delta^{\text{max}}$ | 2.8E-8 | 7.7E-8 | 1.2E-8 | 3.4E-8

∗: Bus 610 is compared phase-to-phase.

##  VI Numerical Experiments

This section discusses the constraint violation and solve time of approximate
balanced and Kron-reduced models, compared to explicit four-wire models.

###  VI-A Formulation variants and calculation time

The results are obtained with explicit variables for the total and series line
current. However, the total current variables can be replaced by affine
transformations of the series current to make the formulation more efficient.
For example, for the four-wire network model, this further reduces the solve
time on average by 40% (note that line current / power are unbounded). Bounds
are often useful to keep thee solver from diverging; for these cases, voltage
magnitude bounds alone were sufficient to do so. The feasible sets are defined
by §II and §III for the IVR and ACR formulations.

###  VI-B Methodology for comparing UBOPF forms

The goal of the numerical experiments is to compare optimization results
obtained with three different models: ‘four-wire’, ‘Kron-reduced’ and ‘single-
phase equivalent’. Figure 7 illustrates how each of these can be represented
by the proposed formulation.

(a) four-wire

(b) Kron-reduced

(c) single-phase equivalent

Figure 7: Illustration of different network models, from more to less
detailed.

Figure 8 shows how the accuracy of the simplified models (Kron-reduced and
single-phase equivalent) is evaluated. Once the dispatch results are obtained,
the ‘true’ voltage profile is then obtained by solving a four-wire PF
_simulation_ problem. The difference between these models typically increases
as the networks are more unbalanced.

Figure 8: The dispatch results obtained with a less detailed OPF model (Kron-
reduced or single-phase equivalent), are evaluated with a four-wire PF model.
The voltages can then be compared, quantifying the approximation error.

###  VI-C Current-voltage formulation

####  VI-C1 Constraint violation

Following the methodology in Figure 8, the approximate balanced and Kron-
reduced models are used to optimize the generator setpoints. These setpoints
are then evaluated by solving the four-wire PF. Figure 9 shows the highest
phase-to-neutral voltage magnitude that this would cause for each of the 128
network instances, for an imposed upper bound of 1.1 pu. Even using the Kron-
reduced model, leads to a constraint violation of 1-4%.

---
network instance
(a) balanced

---
network instance
(b) Kron-reduced

Figure 9: Evaluation of the highest constraint violation (phase-to-neutral
limited to 1.10 pu) caused by approximate network models, for all 128
instances (ordered by violation amount).

Figure 10 shows the highest VUF and neutral voltage magnitude for each
instance when using the four-wire model and only imposing upper and lower
bounds on the phase-to-neutral voltage magnitudes. Figure 10 also shows the
resulting VUF and neutral shift when these are explicitly constrained to 2 %
and 5 V respectively.

---
network instance
(a) balanced

---
network instance
(b) Kron-reduced

Figure 10: When only bounding the phase-to-neutral voltages, the maximum VUF
and neutral shift reach high values, which can be limited with explicit
constraints (ordered by increasing value).

####  VI-C2 Solve time comparison

Figure 11 shows how the solve time increases with the network size, for each
model type. The solve time increases approximately linearly for all model
types, up to networks with about 500 buses. The four-wire model takes about 40
% longer to solve than a Kron-reduced model. This corresponds approximately to
a $\times$4/3 increase in the number of variables, as the four-wire model
typically uses 4 variables for every 3 that the Kron-reduced model introduces.

network size (number of buses)

Figure 11: Each marker shows the solve time for one of the 128 instances,
together with the network size of that instance. This was done for the four-
wire (circles), Kron-reduced (diamonds) and balanced (triangles) model. The
lines show a least-square linear regression for each model respectively. For
the balanced model, outliers at (246, 2.89s), (282, 7.06s) and (355, 6.63s)
were left out.

###  VI-D Power-voltage formulation

####  VI-D1 Solve time comparison

Figure 12 compares the solve time with the ACR to IVR, for each of the 128
network instances. Both were initialized with the same, flat, balanced voltage
profile. Unlike IVR, ACR did not even converge within 500 iterations for 35 %
of the instances, and takes several times longer to solve.

Figure 12: The power-voltage formulation takes a lot longer to solve for the
four-wire model across all instances; 44 instances did not even converge
within 500 iterations.

####  VI-D2 Constraint violation

The authors argue that this is a consequence of the equivalence between KCL
and the power balance, which breaks down when the nodal voltage goes to zero,

$\displaystyle\begin{aligned} &\hskip 15.00002pt\text{KCL}\\\
&\sum{\color[rgb]{0,0,1}I}=0\end{aligned}$ $\displaystyle\begin{aligned}
&\;\;\Rightarrow\\\ &\xLeftarrow[{\color[rgb]{0,0,1}U}\neq 0]{}\end{aligned}$
$\displaystyle\begin{aligned} &\text{power balance}\\\
&\sum{\color[rgb]{0,0,1}S}=0\end{aligned}.$ (30)

The power variables, across all components, are linked to the other variables
by constraints of the form
${\color[rgb]{0,0,1}S}={\color[rgb]{0,0,1}U}({\color[rgb]{0,0,1}I})^{*}$. When
${\color[rgb]{0,0,1}U}=0$, ${\color[rgb]{0,0,1}S}=0$ , even when
${\color[rgb]{0,0,1}I}\neq 0$. Therefore, $\sum{\color[rgb]{0,0,1}S}=0$ is
satisfied, whilst $\sum{\color[rgb]{0,0,1}I}=0$ is no longer enforced.
Clearly, all IVR-feasible solutions are feasible in ACR, but not the other way
around. Therefore, the ACR formulation is a relaxation of IVR, which allows a
non-physical grounding of any of the nodes. Figure 13 summarizes this
graphically.

(a)

(b)

(c)

(d)

Figure 13: Figure (a) shows the topology of a single-phase, 3-bus network as
an example. The ACR formulation for (a) also allows all solutions to the
networks (b-d), which are incorrect, non-physical solutions, introduced by the
lifting of KCL to power variables.

These additional, incorrect solutions can be excluded by putting a lower bound
on the voltage magnitude. In a four-wire network model, when a neutral
terminal is not grounded, the voltage will have a relatively small, non-zero
magnitude. Since it can be arbitrarily close to zero depending on the network
loading, it is very difficult to impose valid lower bounds a priori. The phase
voltage variables however, will rarely drop below e.g. 0.5 pu; often, such
bounds are even included in the problem specification. Because of this, this
problem does not appear when using ACR to solve Kron-reduced and balanced
network models.

The incorrect solutions introduced by ACR, are equivalent to grounding the
neutral at one or many buses. Since doing this typically improves the network
transfer capacity, it can be expected that this is also happens when
maximizing the DG output in the case study. Figure 14 shows that this is
exactly what happens. The neutral voltage variables are initialized at
$\epsilon+j0$ pu. The neutral voltage profile obtained with IVR, is only zero
at the secondary of the transformer, where the neutral is actually grounded.
Moving towards the load throughout the network, the neutral voltage steadily
increases. Independent of $\epsilon$, the objective value always converged to
-2.743.

ACR has many solutions, most of which correspond to solutions which are IVR
infeasible. The result is very sensitive to the value of $\epsilon$, leading
to different profiles for the neutral voltages with non-physical groundings at
intermediate buses. The objective value ranges from -2.780 to -2.757. Even
with an initialization as high as 0.1 pu, the problem persists. The authors
argue that this is an inherent disadvantage of lifting KCL to power variables
for four-wire networks, potentially leading to wrong solutions and unstable
solver runs.

Figure 14: This figure shows the neutral voltage magnitude along the path from
the substation to load 17, for the ENWL network instance (2,5), for solutions
obtained with the IVR and ACR formulation. The dots indicate buses where the
voltage magnitude is nearly zero ($<$1E-8). Figure 15: Summarizing the
properties of the lifted and approximated formulations.

####  VI-D3 Comparison of optimality

Fig. 15 summarizes how the various model approximations and formulations are
related. Reformulating KCL in terms of power flow variables introduces
additional, non-physical solutions if the voltage magnitude is not constrained
below. This is a relaxation, and allows non-physical solutions which can have
a lower objective value. Also, note that these non-physical solutions can be
interpreted as grounding any combination of neutral terminals. Therefore, the
ACR four-wire formulation also includes all solutions of the ACR Kron-reduced
formulation, i.e. when all neutral terminals are grounded.

TABLE II: Objective values for solutions obtained for a selection of ENWL test cases. Note that ‘=’ denotes that the ACR formulation yielded the same objective value as the IVR one. | four-wire | Kron-reduced | balanced
---|---|---|---
n-f | IVR | ACR | IVR | ACR | IVR | ACR
04-2 | -1.5384 | -1.6700 | -1.6589 | = | -2.6392 | =
23-3 | -0.7868 | -1.1020 | -1.1038 | = | -5.3645 | =
10-5 | -2.1162 | -2.6966 | -2.7441 | = | -5.3718 | =
11-4 | -6.3137 | -6.6470 | -6.7165 | = | -8.5624 | =
09-3 | -10.022 | -11.452 | -11.984 | = | -17.529 | =
18-5 | -4.1726 | -4.4614 | -4.6412 | = | -7.3296 | =
05-2 | -2.4577 | -2.5544 | -2.6838 | = | -3.6570 | =
11-1 | -5.8191 | -6.0016 | -6.4291 | = | -8.6007 | =
23-1 | -4.6631 | -4.7166 | -4.8993 | = | -6.0126 | =
12-2 | -3.9818 | -3.9793 | -4.2093 | = | -4.9350 | =
19-4 | -3.3815 | -3.4719 | -3.3767 | = | -5.4268 | =

Table II shows the objective value for a selection of ENWL test cases,
obtained for all formulations shown in Fig. a15. The objective value of the
ACR four-wire solution is usually in-between the IVR four-wire and ACR Kron-
reduced ones, but exceptions in both directions occur. Since ACR is a
relaxation of IVR for the four-wire model, its global optimum should have a
lower objective value. However, these optimization problems are non-linear and
are only solved to _local_ optimality by Ipopt. Fig. 14 already illustrated
how sensitive the obtained solution is to changes in the variable
initialization. It is evident that the solution is not globally optimal for
the entries where the objective value for the ACR four-wire solution is higher
than the ACR Kron-reduced one, because it is a relaxation.

##  VII Conclusions

This article developed two quadratic formulations for unbalanced OPF, one in
current-voltage variables and another one in power-voltage. Both forms include
transformers, lines, shunts, loads and generators. The optimization models are
agnostic of the number of phases, and can therefore model single-phase
equivalent, Kron-reduced and four-wire networks.

The proposed current-voltage formulation is especially well-suited for four-
wire networks. It is shown that using rectangular power-voltage variables
leads to problems, previously unobserved for positive-sequence and Kron-
reduced representations. When the neutral voltage is (close to) 0, the
equivalence between KCL in current and power breaks down.

A study across 128 low-voltage networks illustrated the effect of the network
model detail on the set points obtained for the distributed generators. The
Kron-reduced model tends to overestimate the network capacity, and caused the
phase-to-neutral voltage magnitudes to exceed their upper bound (limited 10 %
rise) by about 1-4 %. The solve time scales linearly with the number of
variables, and the four-wire model takes therefore about 30 % longer to solve.

## Acknowledgement

S. Claeys carried out this research whilst holding a doctoral (PhD) strategic
basic research grant (1S82518N) of Research Foundation - Flanders (FWO).
Andrew Urquhart and Murray Thomson at Loughborough University kindly provided
the detailed line impedance data used in the case study. We also thank David
Fobes and Carleton Coffrin for their assistance in integrating the proposed
formulation into the open-source software package PowerModelsDistribution.jl
[23].

## References

  * [1] V. Klonari, B. Meersman, D. Bozalakov, T. L. Vandoorn, L. Vandevelde, J. Lobry, and F. Vallée, “A probabilistic framework for evaluating voltage unbalance mitigation by photovoltaic inverters,” _Sust. Energy Grids Networks_ , vol. 8, pp. 1–11, Dec. 2016.
  * [2] A. Barbato, A. Dedè, D. Della Giustina, G. Massa, A. Angioni, G. Lipari, F. Ponci, and S. Repo, “Lessons learnt from real-time monitoring of the low voltage distribution network,” _Sust. Energy Grids Networks_ , vol. 15, pp. 76–85, Sep. 2018.
  * [3] W. Kersting, “The computation of neutral and dirt currents and power losses,” in _2004 IEEE PES Power Syst. Conf. Expo._ , pp. 213–218 vol.1.
  * [4] R. Ciric, A. Feltrin, and L. Ochoa, “Power flow in four-wire distribution networks-general approach,” _IEEE Trans. Power Syst._ , vol. 18, no. 4, pp. 1283–1290, Nov. 2003.
  * [5] A. J. Urquhart, “Accuracy of low voltage electricity distribution network modelling,” Ph.D., Loughborough University, 2016. [Online]. Available: https://hdl.handle.net/2134/21799
  * [6] K. Girigoudar, D. K. Molzahn, and L. A. Roald, “Analytical and empirical comparisons of voltage unbalance definitions,” in _2015 Workshop Eng. Appl. - Int. Congr. Eng._ , p. 6.
  * [7] L. R. D. Araujo, D. R. R. Penido, and F. D. A. Vieira, “A multiphase optimal power flow algorithm for unbalanced distribution systems,” _Int. J. Elec. Power Energy Syst._ , vol. 53, no. 1, pp. 632–642, 2013.
  * [8] E. Dall’Anese, G. B. Giannakis, and B. F. Wollenberg, “Optimization of unbalanced power distribution networks via semidefinite relaxation,” in _2012 North American Power Symp._ , pp. 1–6.
  * [9] L. Gan and S. H. Low, “Convex relaxations and linear approximation for optimal power flow in multiphase radial networks,” in _2014 Power Syst. Comp. Conf._ , pp. 1–9.
  * [10] X. Su, M. A. Masoum, and P. Wolfs, “Optimal PV inverter reactive power control and real power curtailment to improve performance of unbalanced four-wire LV distribution networks,” _IEEE Trans. Sustain. Energy_ , vol. 5, no. 3, pp. 967–977, Jul. 2014.
  * [11] A. Gastalver-Rubio, E. Romero-Ramos, and J. M. Maza-Ortega, “Improving the performance of low voltage networks by an optimized unbalance operation of three-phase distributed generators,” _IEEE Access_ , vol. 7, pp. 177 504–177 516, 2019.
  * [12] Y. Liu, J. Li, and L. Wu, “ACOPF for three-phase four-conductor distribution systems: semidefinite programming based relaxation with variable reduction and feasible solution recovery,” _IET Gener. Transm. Dis._ , vol. 13, no. 2, pp. 266–276, 2019.
  * [13] M. Usman, A. Cervi, M. Coppo, F. Bignucolo, and R. Turri, “Centralized OPF in unbalanced multi-phase neutral equipped distribution networks hosting ZIP loads,” _IEEE Access_ , vol. 7, pp. 177 890–177 908, 2019.
  * [14] ——, “Bus injection relaxation based OPF in multi-phase neutral equipped distribution networks embedding wye- and delta-connected loads and generators,” _Int. J. Elec. Power_ , vol. 114, p. 105394, Jan. 2020.
  * [15] ——, “Cheap conic OPF models for low-voltage active distribution networks,” _IEEE Access_ , vol. 8, pp. 99 691–99 708, 2020.
  * [16] Y. Liu, J. Li, and L. Wu, “State estimation of three-phase four-conductor distribution systems With real-time data From selective smart meters,” _IEEE Trans. Power Syst._ , vol. 34, no. 4, pp. 2632–2643, Jul. 2019.
  * [17] S. Claeys, G. Deconinck, and F. Geth, “Decomposition of n-winding transformers for unbalanced optimal power flow,” _IET Gener. Transm. Distrib._ , 2020.
  * [18] V. Rigoni, L. F. Ochoa, G. Chicco, A. Navarro-Espinosa, and T. Gozel, “Representative residential LV feeders: a case study for the north west of england,” _IEEE Trans. Power Syst._ , vol. 31, no. 1, pp. 348–360, Jan. 2016\.
  * [19] A. J. Urquhart and M. Thomson, “Series impedance of distribution cables with sector-shaped conductors,” _IET Gener. Transm. Dis._ , vol. 9, no. 16, pp. 2679–2685, 2015.
  * [20] F. Geth, S. Claeys, and G. Deconinck, “Current-voltage formulation of the unbalanced optimal power flow problem,” in _2020 Workshop on Model. and Simul. of Cyber-Physical Energy Syst._
  * [21] A. Wächter and L. T. Biegler, “On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming,” _Math. Program._ , vol. 106, no. 1, pp. 25–57, Mar. 2006.
  * [22] B. Legat, O. Dowson, J. D. Garcia, and M. Lubin, “Mathoptinterface: a data structure for mathematical optimization problems,” _Informs J. Computing_ , vol. 34, pp. 671–1304, C2, 2 2022.
  * [23] D. M. Fobes, S. Claeys, F. Geth, and C. Coffrin, “Powermodelsdistribution.jl: An open-source framework for exploring distribution power flow formulations,” _Electric Power Systems Research_ , vol. 189, p. 106664, 2020\. [Online]. Available: http://www.sciencedirect.com/science/article/pii/S0378779620304673


# Quenching star formation with low-luminosity AGN winds

Ivan Almeida1,2, Rodrigo Nemmen1,3 and Rogemar A. Riffel4
1Instituto de Astronomia, Geofísica e Ciências Atmosféricas, Universidade de
São Paulo, São Paulo, SP 05508-090, Brazil
2 School of Mathematics, Statistics and Physics, Newcastle University, NE1
7RU, UK
3Kavli Institute for Particle Astrophysics and Cosmology, Stanford University,
Stanford, CA 94305, USA
4Departamento de Física, Centro de Ciências Naturais e Exatas, Universidade
Federal de Santa Maria, 97105-900 Santa Maria, RS, Brazil E-mail:
<EMAIL_ADDRESS>
(Accepted 2023 August 30. Received 2023 August 2; in original form 2023
February 28)

###### Abstract

We present a simple model for low-luminosity active galactic nucleus (LLAGN)
feedback through winds produced by a hot accretion flow. The wind carries
considerable energy and deposits it on the host galaxy at kiloparsec scales
and beyond, heating the galactic gas thereby quenching star formation. Our
model predicts that the typical LLAGN can quench more than $10\%$ of star
formation in its host galaxy. We find that long-lived LLAGN winds from
supermassive black holes (SMBH) with masses $\geq 10^{8}M_{\odot}$ and mass
accretion rates $\dot{M}>10^{-3}\dot{M}_{\rm Edd}\ (0.002M_{\odot}/yr)$ an
prevent gas collapse and significantly quench galactic star formation compared
to a scenario without AGN, if the wind persists over 1 Myr. For sustained wind
production over timescales of 10 Myr or longer, SMBHs with $10^{8}M_{\odot}$
or larger masses have important feedback effects with
$\dot{M}>10^{-4}\dot{M}_{\rm Edd}\ (0.0002M_{\odot}/yr)$.

###### keywords:

black hole physics – galaxies: active – galaxies: evolution – accretion,
accretion discs

††pubyear: 2023††pagerange: Quenching star formation with low-luminosity AGN
winds–References

##  1 Introduction

Once an early-type galaxy forms, that does not mean it will remain quiescent
forever and ever. Early-type galaxies have abundant gas (e.g. Binette et al.
1994) and should also accrete fresh amounts of it. If all this gas cooled and
led to star formation, the global stellar mass density should currently be
larger than observations by a factor of a few (Benson et al., 2003).
Furthermore, the number of galaxies in the red sequence is steadily growing
since the peak epoch of quasars and starbursts (e.g. Bell et al. 2004; Bundy
et al. 2006). This implies that galaxies are still transitioning to
quiescence. Taken together, these are evidence for an unceasing feedback
process which suppresses star formation in red sequence galaxies and keeps it
quenched.

In this work, we explore the possibility that the feedback mechanism keeping
these galaxies quiescent is due to winds from accreting supermassive black
holes (SMBH) hosted in low-luminosity active galactic nuclei (LLAGN). This
idea is quite promising because most SMBH activity in the nearby universe is
happening in LLAGNs (e.g. Ho 2008). These SMBHs are weakly accreting via
radiatively inefficient accretion flows (RIAF; Yuan & Narayan 2014). RIAFs are
prone to producing profuse winds (e.g. Yuan et al. 2015; Almeida & Nemmen
2020; Yang et al. 2021). In addition, there is increasing evidence of a new
class of early-type galaxies hosting galaxy-scale LLAGN winds from spatially
resolved spectroscopy (Cheung et al., 2016; Roy et al., 2021; Sánchez et al.,
2021) and radio observations (Roy et al., 2018).

Given the potential importance of AGN winds in quenching star formation at
late times, here we perform an analytical study of LLAGN winds as a feedback
mechanism. We build a simplified model of RIAF winds based on the latest
results from numerical simulations and analyze how the presence of an LLAGN
could impact the gas and stellar content of a galaxy.

RIAF winds are very hot, subrelativistic and non-collimated. They carry
considerable energy, with powers up to 1% of the rest mass energy
$\dot{M}c^{2}$ associated with accretion (Almeida & Nemmen, 2020). The kinetic
and thermal energy of the ejected wind must be deposited in the environment,
and its most plausible fate is depositing its energy in the interstellar
medium. By exploring the properties of these winds and their impact on the
host galaxy, we tackle the following questions: Are LLAGN powerful enough to
quench star-formation in an early-type galaxy? Can LLAGN winds keep a red-and-
dead galaxy quiescent?

This paper is structured as follows. In section 2, we present the details of
the model. In section 3 we present the results, which include the predicted
relation between LLAGN power and star-formation quenching. We compare our
results to the literature in section 4. Finally, section 5 presents a summary
and some perspectives.

##  2 Model

In order to quantify the effect of LLAGN feedback, we approximated a galaxy as
an isothermal sphere of dark matter with a fixed fraction of gas. The wind
itself is an expanding sphere. In the following subsections, we describe our
model in more details.

###  2.1 Galaxy

We followed Silk & Rees (1998) and modelled the galaxy as an isothermal sphere
characterized by a velocity dispersion $\sigma$. Stars dominate the total mass
of the galaxy’s central region, and only a small fraction is gaseous
corresponding to a fraction $f_{g}\approx 0.05-0.1$ of the total mass. The gas
density profile is described as

$\rho(R)=\frac{f_{g}\sigma^{2}}{2\pi GR^{2}}.$ (1)

The total gas mass enclosed in a radius $R$ is

$\displaystyle M_{\rm gas}(R)$ $\displaystyle=\int_{0}^{R}4\pi
r^{2}\rho(r)dr=\frac{2f_{g}\sigma^{2}R}{G}$ $\displaystyle=9.6\times
10^{9}f_{g}\left(\frac{\sigma}{200\,{\rm
km/s}}\right)^{2}\left(\frac{R}{1\,{\rm kpc}}\right)M_{\odot}$ (2)

and is in the form of atomic hydrogen. The gravitational binding energy
$E_{\rm gal}$ is

$E_{\rm gal}(R)=\frac{3GM_{\rm total}M_{\rm
gas}}{5R}=\frac{6M_{G}\sigma^{2}}{5}.$ (3)

Adopting $f_{g}=0.05$ and replacing equation (2) in (3) gives

$E_{\rm gal}(R)=4.5\times 10^{56}\left(\frac{\sigma}{200\,{\rm
km/s}}\right)^{4}\left(\frac{R}{1\,{\rm kpc}}\right)\text{ erg }.$ (4)

The system is isothermal with a temperature of $T_{\rm Gal}=1.5\times
10^{6}\sigma_{200}^{2}\text{ K }$ where $\sigma_{200}\equiv\sigma/200\,{\rm
km/s}$.

###  2.2 LLAGN Energy Output

The LLAGN is able to inject a $\Delta E$ amount of energy into the galaxy via
winds given by $\Delta E=L_{w}\Delta t$. Being $L_{w}$ the wind power and
$\Delta t$ is the LLAGN lifetime. We parameterise the wind power as a fraction
of the Eddington luminosity, $L_{w}=\eta L_{\rm Edd}$. Following Almeida &
Nemmen (2020), the wind power is $\sim 0.1-1$ per cent of the rest-mass energy
$\dot{M}c^{2}$ accreted by the SMBH. Given that for a LLAGN we expect
$\dot{M}\lesssim 10^{-3}\dot{M}_{Edd}$ and $L_{\rm Edd}\equiv 0.1\dot{M}_{\rm
Edd}c^{2}$, we have $\eta\lesssim 10^{-4}$. Thus, in our calculations we
assume $\eta=10^{-4}$ and thereby

$\Delta E=4\times 10^{56}\left(\frac{\eta}{10^{-4}}\right)\left(\frac{M_{\rm
BH}}{10^{9}\,{\rm M}_{\odot}}\right)\left(\frac{\Delta t}{1\,{\rm
Myr}}\right)\ \text{erg}.$ (5)

With these considerations, the impact of the AGN on the host galaxy increases
trivially with its lifetime and decreases with the distance from the SMBH, as
can be seen by taking the ratio of the LLAGN energy output with the galactic
gravitational binding energy,

$f_{\rm AGN}\equiv\frac{\Delta E}{E_{\rm gal}}=0.24\left(\frac{\Delta
t}{1\,{\rm Myr}}\right)\left(\frac{R}{1\,{\rm
kpc}}\right)^{-1}\left(\frac{M_{\rm BH}}{10^{9}\,{\rm
M}_{\odot}}\right)^{0.22},$ (6)

where we have used the $M-\sigma$ relation from McConnell et al. (2011),

$\left(\frac{M_{\rm BH}}{10^{8}M_{\odot}}\right)=1.95\left(\frac{\sigma}{200\
km\ s^{-1}}\right)^{5.12}.$ (7)

As we will see, the LLAGN energy output can be comparable to the galactic
gravitational binding energy.

###  2.3 Star-formation

Star formation usually occurs in giant molecular clouds (GMC), massive
reservoirs of cold gas prone to star formation.

In our model, we assume that the entirety of the wind kinetic power couples to
GMCs and is converted to thermal energy. This approximation amounts to $f_{\rm
AGN}$ translating directly into the fractional temperature increase caused by
AGN feedback.

We describe the protostellar core mass function as

$\frac{dN}{dlnM}=N_{0}\Big{(}\frac{M}{M_{0}}\Big{)}^{-\xi}\text{, }(M\lesssim
M_{0}).$ (8)

following Rosolowsky (2005); Dib et al. (2008). Equation (8) gives the
distribution of protostellar cores inside GMCs as a function of mass and
sizes. We considered in our model dense clouds with $M_{0}\lesssim
100M_{\odot}$ and $0.3\leq\xi\leq 2.7$ (Drapatz & Zinnecker, 1984; Rosolowsky,
2005; Dib et al., 2008; Mok et al., 2020; Rosolowsky et al., 2021).

Cores able to generate stars are those with masses exceeding the Jeans mass

$M_{J}=20M\odot\Big{(}\frac{T}{10\ {\rm K}}\Big{)}^{1.5}\Big{(}\frac{n}{100\
{\rm cm}^{-3}}\Big{)}^{-0.5}.$ (9)

Assuming, in our model, a constant external pressure around the cloud and
$n=100\ {\rm cm}^{-3}$, equation (9) simplifies to $M_{J}=20M_{\odot}(T/10\
{\rm K})^{2}$.

##  3 Results

###  3.1 Energetics

Figure 1 illustrates the characteristic values of $f_{\rm AGN}$ for a range of
AGN timescales and distances. The figure indicates that an LLAGN can inject a
significant amount of energy into the inner 10 kpcs of the host galaxy. The
effect is more prominent in galaxies with more massive SMBHs. For instance, a
galaxy hosting a $10^{8}M_{\odot}$ SMBH can undergo a $10\%$ temperature
increase in the innermost 2 kpc in one million years; a $10^{9}M_{\odot}$ SMBH
active over 2 Myr with achieve a heating fraction higher than 50%. Moreover,
if the LLAGN is active for 5 Myr or longer, the galactic heating within 5 kpc
will be energetically relevant regardless of the mass.

Figure 1: Energy injected by LLAGN winds scaled by the galactic binding energy
as a function of distance to the supermassive black hole, based on equation 6.
Different AGN durations and black hole masses are displayed in the different
panels, with the mass in solar masses.

###  3.2 How far does the wind reach?

Simulations suggest strong winds coming from RIAFs, with powers reaching up to
one percent of the rest mass associated with accreted gas (Almeida & Nemmen,
2020). These winds have thermal energies greater than the gravitational
binding energy, which means they have enough energy to escape the black hole’s
gravitational sphere of influence. Nevertheless, the spatial extent of these
winds remains an open question. We investigated the wind extension using two
different approaches. In the first one, we model the wind as expanding bubble
which cools via bremsstrahlung. In the second one, we consider a central
heating source and a heat transfer through the gas—here, the wind carries only
energy and not mass.

In the first scenario, we computed the distance travelled by the bubble front
over the cooling time, $R_{\rm wind}=vt_{\rm cool}$, where we assume $v=300\
{\rm km\ s}^{-1}$ (Cheung et al., 2016; Almeida & Nemmen, 2020) and that the
density follows $\rho_{\rm wind}\propto r^{\alpha}$. The resulting expression
is

$R_{\rm wind}=\Bigg{(}\frac{5.9\times
10^{-5}\eta^{-1}\sigma_{200}^{2}M^{\alpha}}{10^{7\alpha}\sqrt{1-\alpha}}\Bigg{)}^{\frac{1}{1+\alpha}}\
{\rm kpc}$ (10)

where we assume $\eta\sim 10^{-4}$, related to the efficiency of the wind
production. This is roughly

$R_{\rm wind}\gtrsim\begin{dcases}3\ {\rm kpc},\alpha<-0.1\\\ 100\ {\rm
kpc},\alpha<-0.3\\\ \end{dcases}$ (11)

We find that for $\alpha<0$, the wind can reach distances larger than ten kpc
which are beyond the visible size of most galaxies.

For the second case, we numerically solve the one-dimensional radial heat
transfer equation for a sphere made of hydrogen with a central heat point
source,

$\frac{1}{r^{2}}\partial_{r}(r^{2}\partial_{r}T)=\frac{\rho
c_{P}r^{2}}{\kappa}\partial_{t}T+Q_{\rm AGN}$ (12)

We modelled the AGN impact as a spherical boundary with constant temperature
and hotter than the medium. This can be translated as the boundary condition
in equation (13) and initial condition in equation (14). For practical
reasons, we assumed $r_{\rm AGN}=0$ since the AGN scales are too small
compared to the galaxy.

$\displaystyle T(r=r_{\rm AGN})\leq T_{AGN}\text{,}$ (13) $\displaystyle
T(t=0,r)=\begin{dcases}T_{\rm AGN},&r\leq r_{\rm AGN}\\\ T_{\rm gal},&r>r_{\rm
AGN}\end{dcases}.$ (14)

Solving equation (12) and assuming the characteristic values from Fabian et
al. (2005) (their equation 4), we found that the resulting temperature profile
follows $T(R)\propto R^{-1}$. This is the same radial dependence as in
equation (6). After about 5 Myr, even gas at kiloparsec scales will undergo a
20% temperature increase. For this model $R_{\rm wind}$ is the radius at which
$\lim_{r\to R_{\rm wind}}T(r)=T_{\rm gal}$. We find that typically $R_{\rm
wind}\gtrsim 1$ kpc.

Both models indicate that winds can get to the galactic outskirts, reaching
distances up to kpc. We stress that the multiscale physics of the ISM and its
interaction with hot winds is quite complex. We leave the numerical modeling
of these phenomena for a future work.

###  3.3 Star formation quenching

The number of protostellar cores able to collapse and form stars can be
calculated using equations 8 and 9 as

$\mathcal{N}(M\geq M_{J})=\int_{M_{J}}^{M_{0}}N(M)dM.$ (15)

We use $\mathcal{N}$ to quantify the impact of LLAGN feedback in quenching
star formation by computing it in two different ways: $\mathcal{N}_{0}$ is the
number of protostellar cores able to collapse into stars when the AGN effect
is not taken into account, whereas $\mathcal{N}_{\rm AGN}$ is the
corresponding quantity with the AGN turned on. In particular, we are
interested in comparing how much lower $\mathcal{N}_{\rm AGN}$ is compared to
$\mathcal{N}_{0}$ as a function of the main accreting BH parameters: the BH
and mass accretion rate. When estimating $\mathcal{N}_{0}$, we consider a
temperature $T_{\rm PC}\sim 10$K and corresponding Jeans mass is denoted by
$M_{J}$ (see equation (9)); for $\mathcal{N}_{\rm AGN}$, we adopt $T_{\rm
PC}^{\rm AGN}=(1+f_{\rm AGN})T_{\rm PC}$ as the AGN increase the average
temperature and the appropriate Jeans mass is $M_{J}^{\rm AGN}$. This implies
that $M_{J}<M_{J}^{\rm AGN}$. Protostellar cores with masses in the range
$M_{J}<m<M_{J}^{\rm AGN}$ will suffer gravitational collapse when the impact
of the AGN is not considered; they would not if the LLAGN is taken into
account.

We define the fraction of star formation quenched by the LLAGN—the quenching
fraction $Q$—as

$Q\equiv 1-\frac{\mathcal{N}_{\rm
AGN}}{\mathcal{N}_{0}}=1-\frac{1-(M_{J}/M_{0})^{1-\xi}(1+f_{\rm
AGN})^{2-2\xi}}{1-(M_{J}/M_{0})^{1-\xi}}.$ (16)

where $\xi$ is a power-law index and $M_{0}$ is the mass scale related to the
protostellar core mass distribution, (see equation (8)). The meaning of $Q$ is
the following: in the extreme case when $Q=1$, the entirety star formation is
aborted due to AGN feedback; on the other hand, when $Q=0$ there is no
quenching at all. Therefore, $Q$ and the star-formation rate are inversely
correlated.

We plot in figure 2 the relation between star formation quenching and the AGN
heating fraction $f_{\rm AGN}$, where we explore the dependence on the
parameter $\xi$ (equation (8)). As expected, quenching becomes more pronounced
as the amount of energy dumped by the LLAGN increase though this proceeds in a
nonlinear fashion.

Figure 2: The quenching fraction as a function of the average heating of the
region. As the temperature increases, the fraction of shut-down stellar
formation sites increases. The different lines represent the different
distribution possibilities for the protostellar cores (see equation (8)).

Figure 3 illustrates the dependence of quenching on the SMBH mass accretion
rate. Each shaded region with a different color corresponds to a given SMBH
mass, with the interior spanning all allowed $\xi$ values assuming $R=20$ kpc
(a typical galaxy size). The different panels explore the impact of the
duration of the LLAGN activity varying from 1 Myr (upper left panel) to 50 Myr
(bottom right panel). For illustration, let’s consider a SMBH accreting at the
$10^{-3}\dot{M}_{\rm Edd}$ level. If its mass is $10^{8}M_{\odot}$
($10^{9}M_{\odot}$) and the wind is produced for only 1 Myr, it can quench
less than one per cent ($5\%$) of star formation in the host galaxy; now, if
the LLAGN is active for 10 Myr, it can quench up $10\%$ ($30\%$); moreover, if
it is active for 50Myr, the quenched grows to $40\%$ ($60\%$).

Figure 3: The plot shows the quenching fraction inside a region of 20kpc as a
function of the LLAGN accretion rate. The increase in the accretion rate has a
significant effect on the gas. Each colour represents a different SMBH mass.
We can observe the importance of the system’s total mass; the quenching only
occurs for the most massive SMBHs. The three different panels refer to the
LLAGN activity time $\Delta t$, long-lived LLAGN have a much more substantial
impact on the gas temperature and subsequent quenching. The denoted regions
represent the different distributions of the protostellar cores (see equation
(8)); they are the region delimited by the lines shown in figure 2.

Figure 4 displays the SMBH activation function for effective AGN feedback, as
predicted in our calculations. This figure displays the family of accreting
SMBH parameters required to produce a ten per cent quenching of star
formation, i.e. the combination of mass accretion rates and masses that result
in $Q=0.1$. Figure 4 shows that a $10^{8}M_{\odot}$ or $10^{9}M_{\odot}$ SMBH
that experiences an accretion episode lasting 1 Myr with $\dot{M}>4\times
10^{-3}\dot{M}_{\rm Edd}$ will be able to abort more than $10\%$ of star
formation in its host galaxy. For an accretion episode lasting 10 Myr, a
$10^{8}M_{\odot}$ SMBH needs $\dot{M}>4\times 10^{-4}\dot{M}_{\rm Edd}$ to
significantly impact its host galaxy via winds; a $10^{9}M_{\odot}$ SMBH needs
$\dot{M}>3\times 10^{-4}\dot{M}_{\rm Edd}$.

Figure 4: Family of accreting SMBH parameters required to produce quenching of
star formation of at least ten per cent, as a function of BH mass. Black (top)
and red (bottom) circles correspond to LLAGN lifetimes of 1 Myr and 10 Myr,
respectively. This figure serve as a guide to evaluate whether a LLAGN
feedback is effective as a function of black hole mass and accretion rate.

Correspondingly, Figure 5 displays the wind power resulting in effective AGN
feedback with $Q\geq 0.1$. Similarly to the story told in Figure 4, a
$10^{8}M_{\odot}$ SMBH that produces a wind lasting 1 Myr with power larger
than $5\times 10^{39}{\rm erg\ s}^{-1}$ will be able to abort more than $10\%$
of star formation in its host galaxy. For winds lasting 10 Myr, a
$10^{8}M_{\odot}$ ($10^{9}M_{\odot}$) SMBH needs a wind power larger than
$4\times 10^{38}{\rm erg\ s}^{-1}$ ($2\times 10^{39}{\rm erg\ s}^{-1}$) for
effective quenching.

Figure 5: Wind power required to produce quenching of star formation of at
least ten per cent. Color scheme is the same of Figure 4. The star indicates
the values measured for the Akira LLAGN.

Overall, the LLAGN will only have an impact larger than ten per cent on the
host galaxy if it persists for durations longer than 10 Myr, regardless of the
SMBH mass. This timescale is one order of magnitude larger than the typical
quasar lifetime. Long LLAGN durations are needed in order to significantly
suppress star formation.

##  4 Discussion

Going back to the questions posed at the beginning of this work: Are LLAGN
powerful enough to quench star-formation in an early-type galaxy? Can LLAGN
winds keep a red-and-dead galaxy quiescent? With our simple models we find
that the answer to both questions is yes. The quenching intensity, however,
depends on the black hole mass, accretion rate and on the duration of the
accretion episode.

The accretion rate is a crucial parameter in our model. By converting the
Eddington units to solar masses per year ($M_{\odot}$/yr), we obtain:

$\dot{M}_{\rm Edd}=2.23\left(\frac{M_{\rm
SMBH}}{10^{8}M_{\odot}}\right)M_{\odot}\ \text{yr}^{-1}.$ (17)

For a black hole with a mass of $10^{8}M_{\odot}$ accreting at
$10^{-3}\dot{M}_{\rm Edd}$, the total mass accreted per year is
$0.002M_{\odot}/\text{yr}$. In our extreme case of a $10^{9}M_{\odot}$ black
hole accreting at $10^{-2}\dot{M}_{\rm Edd}$, the required amount of gas is
approximately $0.2M_{\odot}/\text{yr}$. The mass necessary to trigger our
winds is highly feasible, considering the dense nature of the central regions
of galaxies. Now, if we consider the total amount of gas needed to sustain
these two scenarios for a duration of 1 Myr, we would require $\sim
10^{3}M_{\odot}$ and $\sim 10^{5}M_{\odot}$, respectively. These values still
remain significantly below the mass scales of the supermassive black hole or
the stellar mass of the galactic centre.

Let’s consider now the particular case of the “Akira” galaxy. Cheung et al.
(2016) reported evidence for winds emerging from the LLAGN in Akira. The
authors dubbed this putative class of objects “red geysers” (e.g. Roy et al.
2018). Our work supports the notion that LLAGN winds can indeed be energetic
enough to originate the red geyser phenomenon. Cheung et al. (2016) find that
Akira hosts a $10^{8}M_{\odot}$ SMBH currently accreting with $\lambda\equiv
L/L_{\rm Edd}=4\times 10^{-4}$ and that the wind lasts at least 10Myr. This
value of $\lambda$ corresponds to $\dot{M}=3\times 10^{-3}\dot{M}_{\rm Edd}$,
for a typical RIAF radiative efficiency of $1\%$ (Xie & Yuan, 2012). Our model
predicts that LLAGN winds in Akira can reach quenching fractions of about
$30\%$ if those accretion rates are sustained over 10 Myr, and potentially
much more for longer times. Star formation in the so-called red geyser
galaxies can be significantly impacted by winds produced from underfed SMBHs.

We explored two different assumptions on the radial expansion of the wind.
Both of them indicate that the kinetic and thermal energies can be carried
over kiloparsec scales way beyond the SMBH gravitational sphere of influence.

An important parameter in our results is the activity time of the LLAGN. If we
want to explain the quiescence of the local universe galaxies as the effect of
a steady and weak wind from a very faint AGN, this object must be active for a
very long time. In figure 3, we can see in the left panel that only SMBHs with
mass $M_{\rm SMBH}\gtrsim 10^{9}M_{\odot}$ and $\dot{m}\gtrsim 5\times
10^{-3}$ can noticeably impact the star-formation in $\Delta t_{\rm
AGN}=1$Myr. However, for a longer time as $\Delta t_{\rm AGN}=10$Myr, one
LLAGN with $\dot{m}\gtrsim 10^{-3}$ and masses $M_{\rm SMBH}\gtrsim
10^{8}M_{\odot}$ can turn off more than 50% the stellar formation sites. The
star formation can be severely suppressed if the galaxy inflow can sustain the
LLAGN accretion rate for a long enough time.

One limitation of our model is that we are unable to give more details on
specific types of stellar populations arising after quenching by the LLAGN
winds. Modeling the vast dynamical range and the nonlinear physics involved in
star formation is complex problem and outside the scope of this work – a
simulation of effect feedback for an elliptical galaxy treated in much more
detail can be seen in Yuan et al. (2018). One broad brush consequence of the
suppression of star formation is that there will be a smaller amount of heavy
elements being spewed out throughout the galaxy. Thus, galaxies under the
influence of LLAGN feedback will have smaller metallicities. At the same time,
and for the same reasons, we expect a smaller number of younger stars, so
LLAGN winds tend to redden the host galaxy.

Our model assumes a smooth wind that interacts with molecular clouds, heating
them up over Myr timescales. In a more realistic setting, outflows likely
strip gas clouds. The ensuing cloud mass decrease would further boosting the
quenching fraction to higher values than we reported in figure 3. This
possibility remains to be investigated in the future.

Another aspect worth discussing is the feedback impact on the black hole mass
supply itself, i.e. does the LLAGN feedback reduces significantly $\dot{M}$?
According to Bu & Yang (2019) at first feedback inhibits $\dot{M}$, which also
leads to a reduction in the wind strength; feedback then becomes less
important and accretion eventually is reestablished to its previous levels. Bu
& Yang (2019) find that this cycle occurs on timescales orders of magnitude
shorter than those we consider here. For instance, when the $\dot{M}$ valus
reported on Bu & Yang are averaged over timescales longer than $10^{5}$ years,
the accretion rate is well-behaved (cf. their Figure 2c). This supports our
assumption that $\dot{M}$ does not vary much over timescales longer than 1
Myr.

Furthermore, another argument can be made that $\dot{M}$ is not severely
impacted by the LLAGN feedback. Keeping the assumption of a spherically
symmetric, stationary accretion flow, we consider the outflow depositing
energy at the Bondi radius $r_{B}$. By solving the one-dimensional fluid
conservation equations with appropriate boundary conditions (cf. appendix in
Di Matteo et al. 2003), the Bondi radius is modified due to the outflow
heating as

$r_{B}=\frac{GM}{2c_{s}^{2}}-\frac{2}{3}\frac{Hr_{B}^{2}}{c_{s}^{3}},$ (18)

where the first term on the right side is the usual Bondi radius and the
second term represents the heating effect; $H$ is the heating rate in units of
erg s${}^{-1}\ {\rm g}^{-1}$ given by $H=L_{w}/M_{g}$ where $M_{g}\sim 4\pi
r_{B}^{3}\rho/3$; $\rho$ and $c_{s}$ are evaluated at $r_{B}$. Since
$L_{w}=\eta\dot{M}c^{2}$ and assuming $\dot{M}=\epsilon\dot{M}_{B}$, i.e. the
BH accretion rate is a small fraction of the Bondi rate per ADIOS scenarios,
we find to first order:

$r_{B}\approx\frac{GM}{2c_{s}^{2}}\left[1-2\eta\epsilon\left(\frac{c}{c_{s}}\right)^{2}\right],$
(19)

being $c$ is the speed of light. Assuming typical wind feedback values of
$\eta\sim 0.01$ (e.g. Almeida & Nemmen 2020), $\epsilon\sim 0.01$ (e.g. Yuan &
Narayan 2014) and $c/c_{s}\sim 10^{3}$, we find that $r_{B}$ could be reduced
by $\sim 10\%$ due to the LLAGN feedback. Given that $\dot{M}\propto
r_{B}^{2}$, this implies that the accretion rate is affected only at the per
cent level. This is not enough to impact our assumption of a constant
$\dot{M}$ in the feedback calculations of the paper.

The focus of our model is on the evolutionary processes of a single isolated
galaxy. The consideration of events such as mergers or interactions with other
galaxies, which can result in material transfer between them, falls outside
the scope of this study. We adopt a fixed gas mass fraction ranging from
$0.05$ to $0.1$ of the total system mass. While we recognise that even
isolated galaxies can experience gas infall from the surrounding environment
over time, the typical infall rate at redshift $z\sim 0$ is generally not
significantly larger than $1M_{\odot}/\text{yr}$ (Sancisi et al., 2008; Mollá
et al., 2016). This rate has a minimal impact on the assumed value for the gas
mass fraction within the timescales considered in our study. It is worth
noting that gas infall plays a more substantial role on much larger scales
than those investigated in this work.

##  5 Summary

The main conclusions of our investigation can be summarised as follows:

(i) Low-luminosity active galactic nuclei can have important feedback effects
in their host galaxies by quenching star formation. This occurs via winds
emerging from the hot accretion flow which are able to heat up protostellar
clouds and prevent them from gravitationally collapsing.

(ii) The relevance of star formation quenching by LLAGN feedback is a function
of the SMBH mass, mass accretion rate and the duration of the accretion
episodes. In general, quenching is only relevant for accretion lasting longer
than 1 Myr.

(iii) For an accretion episode lasting 1 Myr, a $10^{8}M_{\odot}$ or
$10^{9}M_{\odot}$ SMBH needs $\dot{M}\gtrsim 10^{-3}\dot{M}_{\rm Edd}$ to
abort more than $10\%$ of star formation.

(iv) For an accretion episode lasting 10 Myr, a $10^{8}M_{\odot}$ or
$10^{9}M_{\odot}$ SMBH needs $\dot{M}\gtrsim 10^{-4}\dot{M}_{\rm Edd}$ to
significantly impact its host galaxy via winds.

(v) LLAGN winds can reach kiloparsec scales, and beyond.

Our model is subject to the limitations of our assumptions, mainly: the
assumption of a spherical isothermal galaxy, steady state, lack of details on
the treatment of the interstellar medium and the wind physics. Despite these
idealizations, we hope that our calculations can offer insights on the galaxy-
SMBH coevolution.

In conclusion, our model demonstrates that feedback via winds from LLAGNs is
an important suppressor of star formation in red sequence galaxies. LLAGNs,
despite their low Eddington ratios, will keep a red-and-dead galaxy quiescent
at late times. Winds from underfed SMBHs offer a third mode of AGN feedback,
in addition to the quasar or radiative mode relevant at the peak of galaxy
mergers, and the radio or jet mode relevant for radio galaxies in galaxy
clusters.

## Acknowledgements

We acknowledge useful discussions with Raniere de Menezes, Paula R. T. Coelho,
Stephane V. Werner, Feng Yuan, Roger Blandford, Thaisa Storchi-Bergmann, and
Ciriaco Goddi. This work was supported by FAPESP (Fundação de Amparo à
Pesquisa do Estado de São Paulo) under grants 2017/01461-2, 2019/10054-7 and
2022/10460-8. We acknowledge funding from an United Kingdom Research and
Innovation grant (code: MR/V022830/1). RN acknowledges a Bolsa de
Produtividade from Conselho Nacional de Desenvolvimento Científico e
Tecnológico. RAR acknowledges support from CNPQ (400944/2023-5 &
404238/2021-1) and FAPERGS (21/2551-0002018-0).

We used Python (Oliphant, 2007; Millman & Aivazis, 2011) to produce all
scientific results and plots of this paper, including several packages such as
NumPy (Van Der Walt et al., 2011), SciPy (Virtanen et al., 2019), and
Matplotlib (Hunter, 2007).

## Data Availability

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

## References

  * Almeida & Nemmen (2020) Almeida I., Nemmen R., 2020, Monthly Notices of the Royal Astronomical Society, 492, 2553
  * Bell et al. (2004) Bell E. F., et al., 2004, ApJ, 608, 752
  * Benson et al. (2003) Benson A. J., Bower R. G., Frenk C. S., Lacey C. G., Baugh C. M., Cole S., 2003, ApJ, 599, 38
  * Binette et al. (1994) Binette L., Magris C. G., Stasińska G., Bruzual A. G., 1994, A&A, 292, 13
  * Bu & Yang (2019) Bu D.-F., Yang X.-H., 2019, The Astrophysical Journal, 871, 138
  * Bundy et al. (2006) Bundy K., et al., 2006, ApJ, 651, 120
  * Cheung et al. (2016) Cheung E., et al., 2016, Nature, 533, 504
  * Di Matteo et al. (2003) Di Matteo T., Allen S. W., Fabian A. C., Wilson A. S., Young A. J., 2003, The Astrophysical Journal, 582, 133
  * Dib et al. (2008) Dib S., Brandenburg A., Kim J., Gopinathan M., André P., 2008, The Astrophysical Journal, 678, L105
  * Drapatz & Zinnecker (1984) Drapatz S., Zinnecker H., 1984, Monthly Notices of the Royal Astronomical Society, 210, 11P
  * Fabian et al. (2005) Fabian A. C., Reynolds C., Taylor G., Dunn R., 2005, Monthly Notices of the Royal Astronomical Society, 363, 891
  * Ho (2008) Ho L. C., 2008, ARA&A, 46, 475
  * Hunter (2007) Hunter J. D., 2007, Computing in science & engineering, 9, 90
  * McConnell et al. (2011) McConnell N. J., Ma C.-P., Gebhardt K., Wright S. A., Murphy J. D., Lauer T. R., Graham J. R., Richstone D. O., 2011, Nature, 480, 215
  * Millman & Aivazis (2011) Millman K. J., Aivazis M., 2011, Computing in Science & Engineering, 13, 9
  * Mok et al. (2020) Mok A., Chandar R., Fall S. M., 2020, The Astrophysical Journal, 893, 135
  * Mollá et al. (2016) Mollá M., Díaz Á. I., Gibson B. K., Cavichia O., López-Sánchez Á.-R., 2016, Monthly Notices of the Royal Astronomical Society, 462, 1329
  * Oliphant (2007) Oliphant T. E., 2007, Computing in Science & Engineering, 9, 10
  * Rosolowsky (2005) Rosolowsky E., 2005, Publications of the Astronomical Society of the Pacific, 117, 1403
  * Rosolowsky et al. (2021) Rosolowsky E., et al., 2021, Monthly Notices of the Royal Astronomical Society, 502, 1218
  * Roy et al. (2018) Roy N., et al., 2018, ApJ, 869, 117
  * Roy et al. (2021) Roy N., et al., 2021, The Astrophysical Journal, 913, 33
  * Sánchez et al. (2021) Sánchez S., Walcher C., Lopez-Cobá C., Barrera-Ballesteros J., Mejía-Narváez A., Espinosa-Ponce C., Camps-Fariña A., 2021, Revista mexicana de astronomía y astrofísica, 57, 3
  * Sancisi et al. (2008) Sancisi R., Fraternali F., Oosterloo T., Van Der Hulst T., 2008, The Astronomy and Astrophysics Review, 15, 189
  * Silk & Rees (1998) Silk J., Rees M. J., 1998, arXiv preprint astro-ph/9801013
  * Van Der Walt et al. (2011) Van Der Walt S., Colbert S. C., Varoquaux G., 2011, Computing in Science & Engineering, 13, 22
  * Virtanen et al. (2019) Virtanen P., et al., 2019, arXiv e-prints, p. arXiv:1907.10121
  * Xie & Yuan (2012) Xie F.-G., Yuan F., 2012, MNRAS, 427, 1580
  * Yang et al. (2021) Yang H., Yuan F., Yuan Y.-F., White C. J., 2021, ApJ, 914, 131
  * Yuan & Narayan (2014) Yuan F., Narayan R., 2014, ARA&A, 52, 529
  * Yuan et al. (2015) Yuan F., Gan Z., Narayan R., Sadowski A., Bu D., Bai X.-N., 2015, ApJ, 804, 101
  * Yuan et al. (2018) Yuan F., Yoon D., Li Y.-P., Gan Z.-M., Ho L. C., Guo F., 2018, The Astrophysical Journal, 857, 121