<EMAIL_ADDRESS>
NASA Hubble Fellowship Program, Einstein Fellow]

# Collisionless accretion onto black holes: dynamics and flares

Alisa Galishnikova  [  Department of Astrophysical Sciences, Princeton
University, 4 Ivy Lane, Princeton, NJ 08544, USA  Alexander Philippov
Department of Physics, University of Maryland, College Park, MD 20742, USA
Eliot Quataert  Department of Astrophysical Sciences, Princeton University, 4
Ivy Lane, Princeton, NJ 08544, USA  Fabio Bacchini  Centre for mathematical
Plasma Astrophysics, Department of Mathematics, KU Leuven, Celestijnenlaan
200B, B-3001 Leuven, Belgium  Royal Belgian Institute for Space Aeronomy,
Solar-Terrestrial Centre of Excellence, Ringlaan 3, 1180 Uccle, Belgium  Kyle
Parfrey  School of Mathematics, Trinity College Dublin, Dublin 2, Ireland
Bart Ripperda  [  Department of Astrophysical Sciences, Princeton University,
4 Ivy Lane, Princeton, NJ 08544, USA  School of Natural Sciences, Institute
for Advanced Study, 1 Einstein Drive, Princeton, NJ 08540, USA  Center for
Computational Astrophysics, Flatiron Institute, 162 Fifth Avenue, New York, NY
10010, USA

###### Abstract

We study the accretion of collisionless plasma onto a rotating black hole from
first principles using axisymmetric general-relativistic particle-in-cell
simulations. We carry out a side-by-side comparison of these results to
analogous general-relativistic magnetohydrodynamic simulations. Although there
are many similarities in the overall flow dynamics, three key differences
between the kinetic and fluid simulations are identified. Magnetic
reconnection is more efficient, and rapidly accelerates a nonthermal particle
population, in our kinetic approach. In addition, the plasma in the kinetic
simulations develops significant departures from thermal equilibrium,
including pressure anisotropy that excites kinetic-scale instabilities, and a
large field-aligned heat flux near the horizon that approaches the free-
streaming value. We discuss the implications of our results for modeling
event-horizon scale observations of Sgr A* and M87 by GRAVITY and the Event
Horizon Telescope.

Introduction.—The recent high-resolution images of synchrotron emission around
the central black holes (BH) in M87 and the Milky Way (Sgr A*) captured by the
Event Horizon Telescope (EHT) reveal asymmetric ring-like structures around
the event horizon EHT1 ; EHT2022 . The radiation is produced by relativistic
plasma on event-horizon scales. General-relativistic magnetohydrodynamic
(GRMHD) simulations are a conventional tool for modeling accretion onto BHs
EHTcode . In conjunction with GR radiative transfer, one can predict aspects
of the observed radiation, including spatially resolved images, from these
numerical models GRRTcomparison . This theoretical framework allows for a
direct comparison of GRMHD simulations and observations. However, the
accreting plasma in these systems is collisionless, which makes the
simplifying assumptions of GRMHD formally inapplicable. Theoretical models
thus require a kinetic approach, which describes collisionless plasmas from
first-principles. In this Letter we present global GR kinetic simulations of
BH accretion and determine the ways in which they differ from conventional
fluid models.

Supermassive BHs show emission across the electromagnetic spectrum. Besides a
relatively constant background emission, Sgr A* also exhibits episodic bright
flares in the near-infrared and X-rays (e.g., Dodds-Eden_2009 ; Xray ; GRAVITY
). The observed power-law emission implies a presence of accelerated particles
(electrons and possibly positrons) near the BH. Studying the generation of
non-thermal particles is not possible within GRMHD fluid models and requires a
kinetic approach. Additionally, GRMHD does not accurately capture reconnection
of magnetic field lines, which is conjectured to be responsible for particle
energization and flares Markoff_2005 ; BroderickLoeb_2006 ; Dodds-Eden_2009 ;
Dexter2020 ; Ripperda2020 ; Ripperda2022 . Specifically, the rate of
reconnection, which regulates the energization efficiency and can be
responsible for the duration of flares Bransgrove2021 , is known to be
substantially faster in collisionless plasma (e.g., Bhattacharjee2009 ;
SironiSpitkovsky2014 ). Ideal GRMHD also assumes an isotropic Maxwellian
plasma distribution function, while collisionless plasmas easily develop
pressure anisotropy along and across the magnetic field direction Kunz2014 ;
Riquelme2015 , which leads to the development of plasma instabilities. The
saturation of these instabilities regulates the thermodynamic state of the
plasma, potentially affecting the global accretion dynamics and its
observational properties.

Methods. In order to study accretion of collisionless plasmas onto BHs, we
perform GR kinetic simulations using the particle-in-cell (PIC) code ZELTRON
which solves Maxwell’s equations and the equations of motion for individual
macroparticles in the 3+1 formalism GRPIC . We also study the same problem
with identical initial conditions using the GRMHD code Athena++ White_2016 ,
which allows for a side-by-side comparison of the two approaches. Both
approaches utilize horizon-penetrating Kerr-Schild coordinates. We measure
distance in units of the gravitational radius, $r_{g}=GM/c^{2}$, where $G$ is
the gravitational constant, $M$ is the mass of the BH, and $c$ is the speed of
light; time is measured in light crossing times of the gravitational radius,
$r_{g}/c$.

In GRPIC, we resolve all the microphysical plasma scales and respect the
correct hierarchy of scales, i.e. all plasma scales are significantly smaller
than $r_{g}$. Due to their high computational cost, our simulations are
limited to two-dimensional, axisymmetric accretion onto a BH in the $r-\theta$
plane, which is aligned with the BH spin, $a=0.95$. For simplicity, and
motivated by the relevance to the accretion flow onto Sgr A* Quataert2003 , we
start with a zero-angular-momentum spherically symmetric distribution of
stationary plasma. Previous work shows that this accretion problem behaves
similarly in many respects to one incorporating rotating initial conditions
Ressler_2020 . We set an initially constant density, pressure, and magnetic
field aligned with the spin axis throughout the box, and add randomly
distributed magnetic loops Supplemental  to mimic magnetized turbulence
expected in real accretion disks.

We initialize a thermal plasma with $k_{B}T_{inj}\approx 0.02m_{i}c^{2}$,
which corresponds to a Bondi radius $r_{B}=2GM/c_{s}^{2}\approx 50r_{g}$,
where $k_{B}$ is the Boltzmann constant and $c_{s}$ is the sound speed. We
focus on two initial values of the plasma-$\beta$ parameter,
$\beta_{0}=P/P_{B_{0}}=4$ or $10$, where $P$ is the gas pressure and
$P_{B_{0}}$ is the magnetic pressure of the initial vertical magnetic field
$B_{0}$. Our GRPIC simulations use mass ratios of $m_{i}/m_{e}=1$ and $3$; the
thermal Larmor radius of ions is set to be
$\rho_{L}=v_{th}m_{i}c/eB_{0}=0.018r_{g}$. Below we show results of
simulations with mass ratio $m_{i}/m_{e}=1$, which we are able to run until
$1000r_{g}/c$. The dynamics of ions in simulations with $m_{i}/m_{e}=3$ until
$100r_{g}/c$ is similar compared to the case with $m_{i}/m_{e}=1$ Supplemental
.

Our simulation domain extends from the inner boundary, located below the event
horizon, to the outer boundary at $100r_{g}$. We employ constant boundary
conditions at the outer boundary in GRMHD, and use absorbing boundary
conditions supplemented with injection of fresh plasma in GRPIC Supplemental .
In highly magnetized regions plasma density can get depleted; we therefore
apply a ceiling value for the magnetization parameter $\sigma=B^{2}/(4\pi
n(m_{e}+m_{i})c^{2})\approx 30$ in both approaches. In GRPIC simulations, we
add electron-positron pairs when $\sigma$ is above this threshold, thus
mimicking a pair cascade expected in these regions BeskinIstominParev1992 ;
FFPairProduction1998 ; Crinquand2020 ; ChenYuan2020 .

Figure 1: Comparison of fluid and kinetic simulations of accretion onto a
rotating BH. The color shows $n/n_{0}$, which corresponds to the plasma number
density in GRMHD (left) and ion and positron number density in GRPIC (right)
in a region close to the BH ($20r_{g}$). The region inside the BH event
horizon, $r_{h}=r_{g}(1+\sqrt{1-a^{2}})$, is shown by a black circle. Thin
black lines represent the magnetic field lines, a thick black line outlines
the ergosphere. The same quantities are shown at a time of $160r_{g}/c$ (a)
and $300r_{g}/c$ (b). Several current sheets form and reconnect (a), leading
to a flaring state (b).

Results.—In Fig. 1 we show a side-by-side comparison of the evolution of the
number density $n$ of the accreting plasma in GRMHD (plasma number density,
left) and GRPIC (ion and positron number density, right) simulations with
$\beta_{0}=4$ at a time of $160r_{g}/c$ (a) and $300r_{g}/c$ (b); $n_{0}$
corresponds to the initial value. Both panels show a zoom into the inner
$20r_{g}$. Initially, the accreting plasma is free-falling onto the BH within
the Bondi radius, dragging the magnetic field lines towards the event horizon
(a). Inflow streams with a similar structure in both GRPIC and GRMHD are
formed: a thin inflow is formed just below the equator, and a squeezed large
loop is accreting above the equator. As the accretion proceeds, the BH’s
rotation and magnetic flux on the event horizon, which becomes dynamically
important, lead to the launching of magnetically dominated outflows. The
accretion stalls when the magnetic field becomes too strong, leading to
thinning of the inflow streams into current sheets, onset of reconnection (b)
Ripperda2020 ; Ripperda2022 , and the evacuation of the accretion flow in the
equatorial plane (eruption). Eventually, the magnetic loop above the equator
evacuates as an outflowing density bubble. Since the rate of collisionless
reconnection is faster by a factor of a few, compared to MHD, the two
numerical results ultimately diverge.

Figure 2: Evolution of averaged quantities in GRPIC (solid lines) and GRMHD
(dotted). (a) Time evolution of the accretion rate in units of the Bondi
accretion rate, $\dot{M}/\dot{M}_{B}$; (b) magnetic flux on the horizon
normalized by its initial value $\Phi/\Phi_{0}$. The gray dashed lines give
the exponential fit for the decay rate of $\Phi/\Phi_{0}$, highlighting the
difference of the reconnection rate. (c) Mean profiles of number density,
$\langle n/n_{0}\rangle$; (d) temperature, $\langle T\rangle$, and (e)
plasma-$\beta=\langle P\rangle/\langle P_{B}\rangle$, as a function of radius
$r/r_{g}$, averaged over $t=[100-200]r_{g}/c$. The runs initialized with the
turbulent magnetic field are shown by a darker color (loops), with a vertical
uniform magnetic field — a lighter color (no loops).

The reconnection physics manifests itself in the time evolution of the
accretion rate, $\dot{M}=-\int_{\theta}\int_{\phi}\sqrt{-g}\rho u^{r}d\theta
d\phi$, and magnetic flux on the horizon,
$\Phi=0.5\int_{\theta}\int_{\phi}\sqrt{-g}\lvert B^{r}\rvert d\theta d\phi$,
shown in Fig. 2a-b, which we normalize by the Bondi accretion rate,
$\dot{M}_{B}$, and initial value of the magnetic flux on the horizon,
$\Phi_{0}$. Here, $g$ is the metric determinant, and $u^{\mu}$ is the fluid
4-velocity. Initially the infalling plasma in both approaches causes an
increase in $\dot{M}/\dot{M}_{B}$ and $\Phi/\Phi_{0}$ at a similar rate,
reaching saturation in the magnetically arrested state Tchekhovskoy2011 .
GRMHD shows two $\dot{M}$ maxima followed by post-accretion eruption events
associated with the onset of the decline of $\Phi$ at $\approx
200,600r_{g}/c$. The GRPIC simulation, however, shows one $\dot{M}$ maximum
followed by an eruption event at $\approx 380r_{g}/c$ and another accretion
period starting at $\approx 800r_{g}/c$. Therefore, even though reconnection
is more efficient, the variability — the frequency of the eruption events —
might be smaller in the kinetic approach, as the accretion stalls due to its
regulation by the efficient large-scale reconnection. Consequently, both
$\dot{M}/\dot{M}_{B}$ and $\Phi/\Phi_{0}$ saturate at smaller values over a
longer time period in GRPIC.

A comparison of the two approaches due to the reconnection physics alone is
demonstrated by simulations with an initially vertical uniform magnetic field
(no loops, in Fig. 2a-b). Here, since kinetic reconnection is more efficient,
GRPIC shows a steeper exponential decline in $\Phi/\Phi_{0}$ compared to GRMHD
Bransgrove2021 .

We show radial profiles (integrated over $\theta$ during the first accretion
event, $100-200r_{g}/c$) of density $\langle n/n_{0}\rangle$ (c), temperature
$\langle T\rangle$ (d), and $\beta=\langle P\rangle/\langle P_{B}\rangle$ (e)
Supplemental . We find a striking similarity of the number density profiles,
while the temperature and $\beta$-profiles show a significant difference
between the two approaches due to the non-ideal physics described next.

Figure 3: Anisotropy and heat flux in GRPIC simulation initialized with
$\beta_{0}=10$. Each row corresponds to a different moment in time:
$40r_{g}/c$ (a-c), and $100r_{g}/c$ (d-f). First column: number density (left)
and magnetic field (right) fluctuations. Second column: temperature anisotropy
$T_{\parallel}/T_{\perp}$ (left) and the ratio of non-ideal and ideal matter
stress-energy tensor components in the tetrad frame
$\mathcal{T}^{(0)(3)}/\mathcal{T}^{(3)(3)}$, which represents the ratio of
parallel heat flux $q_{\parallel}$ and parallel pressure $P_{\parallel}$. The
black circle and black lines represent the event horizon, ergosphere, and
magnetic field lines as in Fig.1. Third column: probability density as a
function of $\beta_{\parallel}=P_{\parallel}/P_{B}$ and temperature anisotropy
$T_{\perp}/T_{\parallel}$. White dashed lines correspond to the boundaries of
the growth rate of mirror (top) and firehose (bottom) instabilities exceeding
10% of the ion cyclotron frequency calculated using NHDS . A mirror
instability develops and saturates at $40r_{g}/c$, for which we show a zoom
into the structure of the instability (a).

To quantify another key feature of the kinetic approach, departure of the
accreting plasma from thermal equilibrium, we calculate the matter stress-
energy tensor in GRPIC, $\mathcal{T}^{\mu\nu}_{\rm matter}$, which can be used
to derive the pressure tensor, $P^{\mu\nu}$, and the heat flux, $q^{\mu}$. We
project these quantities onto a tetrad, $e^{\mu}_{(\nu)}$, where
$e^{\mu}_{(0)}$ is directed along the (Eckart) fluid velocity and
$e^{\mu}_{(3)}$ is along the magnetic field in the fluid frame. The pressure
tensor in this tetrad frame is diagonal,
$diag(P_{\perp},P_{\perp},P_{\parallel})$, where parallel and perpendicular
components are measured with respect to the magnetic field direction in the
fluid frame.

Ion quantities from a GRPIC simulation initialized with $\beta_{0}=10$ are
shown in Fig. 3, where the two rows correspond to times of $40r_{g}/c$ and
$100r_{g}/c$. The first column shows number density $n/\langle n\rangle$
(left) and magnetic field $\sqrt{B^{2}/\langle B^{2}\rangle}$ (right)
fluctuations, second – temperature anisotropy $T_{\parallel}/T_{\perp}$ (left)
and ratio of the two stress-energy tensor components
$\mathcal{T}^{(0)(3)}/\mathcal{T}^{(3)(3)}$ in the tetrad frame (right), where
$\mathcal{T}^{(0)(3)}$ corresponds to parallel heat flux $q_{\parallel}$
(absent in ideal GRMHD). The last column shows an ion probability density
plotted as a function of $\beta_{\parallel}=P_{\parallel}/P_{B}$ and
$T_{\perp}/T_{\parallel}$. The polar inflow of the plasma leads to the build
up of magnetic field, and an associated increase in $P_{\perp}\propto B$. The
deviation from thermal equilibrium with $T_{\perp}>T_{\parallel}$ leads to the
excitation of small-scale plasma density and magnetic field fluctuations (a),
where we also show a zoom into a small region. This is a kinetic-scale mirror
instability which develops when plasma crosses a $\beta$-dependent temperature
anisotropy threshold, as shown in (c). The saturated strength of the magnetic
field fluctuations, $|1-\sqrt{B^{2}/\langle B^{2}\rangle}|$, is of order of a
few tens of percent, consistent with local simulations Kunz2014 ; Riquelme2015
; ICI . At earlier times in the simulation, transiently, we observe the
development of the electromagnetic firehose instability in the equatorial
region, where $T_{\parallel}>T_{\perp}$ Supplemental .

At early times (b), the effective collisions due to particle scattering by the
kinetic-scale fluctuations leads to a suppression of the heat flux,
$\mathcal{T}^{(0)(3)}\approx 0.1\mathcal{T}^{(3)(3)}$. We find that
$q_{\parallel}$ is also $\approx 0.1$ of the value corresponding to the free
streaming of particles along magnetic field lines Supplemental . As the
accretion proceeds and the value of $\beta$ near the event horizon drops below
$1$, the inflow is no longer accompanied by significant density or magnetic
field fluctuations (d), which is consistent with the plasma being pushed away
from the pressure-anisotropy instability boundaries (f). The absence of
scattering on micro-scale magnetic field fluctuations leads to larger values
of the non-ideal components of the stress-energy tensor,
$\mathcal{T}^{(0)(3)}/\mathcal{T}^{(3)(3)}\approx 1$, outside of the current
sheets (e), and $q_{\parallel}$ also approaches the free-streaming value.
These non-ideal effects contribute to substantial differences between GRPIC
and GRMHD temperature profiles (Fig. 2d).

Figure 4: Panel (a) shows the presence (yes or no) of energetic electrons with
$\gamma\geq 8$ throughout the GRPIC simulation at $300r_{g}/c$ (shown in Fig.
1b, right). To avoid low density regions near the axis, we show regions above
the density threshold, $n>n_{0}/2$. Panel (b) shows particle spectra measured
inside the outflowing bubble (blue) and in the current sheet (purple), which
are outlined in (a) by the respective colors.

To understand where particles are accelerated, we highlight regions with
highly energetic particles with a Lorentz factor $\gamma\geq 8$ in GRPIC at
$300r_{g}/c$ in Fig. 4 (the same time snapshot as Fig. 1b, right). These
particles are predominantly located around the current sheets and the
outflowing dense bubble. We measure particle spectra (b) in the regions
outlined by corresponding colored wedges in (a). The spectral slope of
$\approx-3$ (current sheet) is consistent with particle acceleration in
relativistic magnetic reconnection for the measured magnetization parameter
$\sigma\approx 5$ in the upstream SironiSpitkovsky2014 ; Guo2014 ; Werner2016
. As $\sigma$ increases due to evacuation of plasma in the jet region, we find
a harder slope of $\approx-2$, consistent with higher $\sigma\gtrsim 10$
around the current sheet. Positively charged particles are accelerated more
efficiently (Fig. 4b) because their acceleration by the electric fields inside
the current sheet is aligned with the direction of the outflow motion above it
Cerutti2015 .

Discussion.—A fully kinetic approach is crucial for understanding the dynamics
of plasmas accreting onto supermassive BHs such as Sgr A* and M87*. Using
global GRPIC simulations of accretion onto a rotating BH, we highlight three
significant differences relative to matched GRMHD simulations: (1) differences
in the physics of magnetic reconnection can lead to less frequent eruption
episodes in GRPIC; (2) GRPIC includes pressure anisotropy with respect to the
magnetic field and the associated kinetic instabilities; (3) in GRPIC a large
field-aligned heat flux near the horizon is important in regulating the plasma
temperature. Our kinetic approach allows for self-consistent modeling of
particle acceleration during flaring episodes powered by magnetic
reconnection, which opens up a unique opportunity for comparing theory with
observed radiation spectra and light curves. Studying the relative heating and
acceleration of ions and electrons will require a more realistic ion-to-
electron mass ratio, which we currently lack in our simulations. In
conjunction with general-relativistic radiative transfer, extension of our
simulations to larger mass ratio and 3D will allow us to compare spatially
resolved images, polarization maps and lightcurves constructed from GRPIC
simulations to GRAVITY and EHT data. Future GRPIC simulations will rigorously
measure the non-ideal corrections to the GRMHD stress tensor for a range of
plasma conditions, which can then be included in GRMHD simulations
Chandra_2015 ; Foucart2017  to improve their realism.

###### Acknowledgements.

This work was supported by NASA grant 80NSSC22K1054 and NSF grant PHY-2231698.
EQ and AG were supported in part by a Simons Investigator grant from the
Simons Foundation. Computing resources were provided and supported by
Princeton Institute for Computational Science and Engineering; and by the VSC
(Flemish Supercomputer Center), funded by the Research Foundation Flanders
(FWO) and the Flemish Government – department EWI. This research is part of
the Frontera computing project at the Texas Advanced Computing Center (LRAC-
AST21006). Frontera is made possible by NSF award OAC-1818253. F.B.
acknowledges support from the FED-tWIN programme (profile Prf-2020-004,
project “ENERGY”) issued by BELSPO. Support for this work was provided by NASA
through the NASA Hubble Fellowship grant HST-HF2-51518.001-A awarded by the
Space Telescope Science Institute, which is operated by the Association of
Universities for Research in Astronomy, Incorporated, under NASA contract
NAS5-26555. This research was facilitated by Multimessenger Plasma Physics
Center (MPPC), NSF grant PHY-2206607.

Supplemental Materials

In this Supplemental Material we provide additional details about our fluid
GRMHD and kinetic GRPIC simulations, as well as the calculation of the non-
ideal stress-energy tensor in GRPIC.

##  I Setup

We study accretion of plasma onto a spinning black hole with GRMHD and GRPIC
simulations using boundary conditions described below. Both approaches utilize
Kerr-Schild coordinates, which allows us to set the inner boundary inside the
event horizon. At the outer boundary of the GRMHD simulations, we set constant
boundary conditions which fix the fluid variables to their initial conditions.
We mimic this in GRPIC simulations by keeping the plasma number density at a
constant value of $n_{e}=0.9n_{0}$ at $r/r_{g}\in[95-96]$, close to the outer
boundary, by injecting particles. We also introduce an absorption layer for
the electromagnetic fields at the outer boundary, $r/r_{g}=100$, to allow an
outflow of the escaping electromagnetic waves (Cerutti2015, ).

MHD does not permit vacuum, hence we utilize radius-dependent density and
pressure floors in GRMHD simulations, $\rho\geq
10^{-6}\rho_{0}(r/r_{g})^{-3/2}$ and $P\geq 3.33\times
10^{-9}P_{0}(r/r_{g})^{-5/2}$, as well as additional limits for the plasma
magnetization parameter, $\sigma=2P_{B}/\rho\leq 30$, and $\beta=P/P_{B}\geq
10^{-3}$. In GRPIC, we utilize a pair-injection scheme close to the BH,
$r<15r_{h}$, by using a magnetization parameter threshold, $\sigma\approx 30$,
above which we introduce electron-positron pairs into the system. Due to these
floors, as the system evolves in time and $\beta$ decreases close to the BH, a
substantial region becomes affected by the floors. Thus, the comparison of
density and temperature profiles at late times has to be taken with caution.
In more realistic three-dimensional case, accretion is not completely halted
during the eruptions (Tchekhovskoy2011, ), and the effect of floors is less
significant.

In both approaches, we initialize plasma at rest (all spatial components of
$u^{\mu}$ are set to $0$) with a constant density ($\rho_{0}$ in GRMHD,
$n_{0}$ in GRPIC) and temperature. The magnetic field is initialized via an
out-of-plane vector potential, $A_{\phi}=A_{0}+A_{\rm loops}$, represented by
a sum of a uniform background vertical field along the $\hat{z}$-axis and
additional randomly distributed magnetic loops:

$\displaystyle A_{0}$ $\displaystyle=\frac{1}{2}B_{0}r^{2}\sin^{2}\theta,$
(1a) $\displaystyle A_{n}$
$\displaystyle=r(L_{n}-r_{c,n})\exp{\Big{(}1-\frac{L_{n}^{2}}{(L_{n}-r_{c,n}^{2})}\Big{)}},$
(1b) $\displaystyle A_{\rm loops}$ $\displaystyle=\begin{cases}0,\text{if}\
r_{c}<L_{n}\\\ \sum\limits_{n}^{N}kB_{0}A_{n},\ \text{otherwise},\\\
\end{cases}$ (1c)

Here, $B_{0}$ is the magnetic field strength in the absence of the loops,
which corresponds to the initial plasma-$\beta_{0}$. Here $x_{c,n}$,
$y_{c,n}$, and $L_{n}$ set the position of the $n$-th loop center and its
size, $r_{c,n}=\sqrt{(x-x_{c,n})^{2}+(y-y_{c,n})^{2}}$ is the distance from
the loop’s center, $N=1000$ is the total number of loops. The normalization
coefficient $k$ is chosen such that magnetic energy that is contained in the
background vertical magnetic field is equal to the magnetic energy that is
contained in the turbulent loops. Additionally, an exponential hole in
$A_{\phi}$ and initial density is set by multiplying by an additional factor
of $\exp{(5(1-6/r))}$ close to the BH at $r<6r_{g}$, same as in Ressler_2020 .
The loops are randomly distributed across the simulation box with random sizes
ranging from $r_{g}$ to $20r_{g}$, avoiding the spin axis, $x_{c}>20r_{g}$.

Figure 5: Comparison of fluid and kinetic simulations of plasma accreting onto
a rotating BH. The color shows $n/n_{0}$, which corresponds to the plasma
number density in GRMHD (left) and ion and positron number density in GRPIC
(right) in a full simulation box (up to $100r_{g}$). Thin black lines
represent magnetic field lines. The same simulations are shown at a time of
$0r_{g}/c$ (a), $160r_{g}/c$ (b), and $300r_{g}/c$ (c).

Our kinetic simulations are initialized with $20$ particles per cell ($10$
electrons and $10$ ions), which are thermally distributed with the thermal
velocity, $v_{th}$, corresponding to the Bondi radius, $r_{B}\approx 50r_{g}$.
In the absence of loops, the ion Larmor radius corresponds to
$v_{th}m_{i}c/eB_{0}\approx 0.018r_{g}$. Our grid is logarithmically
distributed in $r$, thus, the value of the electron skin depth
$d_{e}=c/\sqrt{4\pi n_{0}e^{2}/m_{e}}$ at the outer cell sets the resolution
of our simulations.

We utilize horizon-penetrating Kerr-Schild coordinates, logarithmic in $r$.
Our GRMHD simulations have resolution of $N_{r}\times N_{\theta}=4096\times
4096$ cells. GRPIC simulations with $\beta_{0}=4$ have a resolution of
$N_{r}\times N_{\theta}=7840\times 7936$ cells, while those with
$\beta_{0}=10$ have $N_{r}\times N_{\theta}=11720\times 11520$ cells,
representing the largest-to-date global kinetic simulations. Figure 5a shows
the equivalent initial conditions ($t=0r_{g}/c$) for the GRMHD (left) and
GRPIC (right) approaches with $\beta_{0}=4$. Color represents density,
initially uniform with an exponential hole in the center, black lines
represent magnetic field lines. We also show the time evolution of the setup
at $t=160r_{g}/c$ (b) and $t=300r_{g}/c$ (c) over the full box (compare to
Fig.1 in the main text which shows the inner $20r_{g}$ of the computational
domain).

##  II Calculation of the matter stress-energy tensor, pressure anisotropy and
heat flux

To quantify the departure of the accreting plasma from thermal equilibrium, we
calculate the matter stress-energy tensor in the kinetic simulations for every
species as $\mathcal{T}^{\mu\nu}_{{\rm
matter}}\equiv\int\frac{d^{3}p}{\sqrt{-g}p^{t}}p^{\mu}p^{\nu}f$, where $f$ is
the distribution function, and $p^{\mu}$ is the particle contravariant
momentum. The stress-energy tensor includes components that are present in
ideal GRMHD and some that are not. To distinguish between these, we find the
(Eckart) fluid velocity $u^{\mu}\equiv\frac{N^{\mu}}{\sqrt{N_{\nu}N^{\nu}}}$,
where $N^{\mu}\equiv\int\frac{d^{3}U}{\sqrt{-g}}(U^{\mu}/U^{t})f$ and
$U^{\mu}$ is the particle contravariant velocity. The pressure tensor is
derived as the spatial part of the stress-energy tensor,
$P^{\mu\nu}\equiv\Delta^{\mu}_{\alpha}\Delta^{\nu}_{\beta}\mathcal{T}^{\alpha\beta}_{\rm
matter}$, where $\Delta^{\mu\nu}\equiv g^{\mu\nu}+u^{\mu}u^{\nu}$ is the
projection tensor. In a non-ideal plasma, the time-space components of the
stress-energy tensor set the heat flux
$q^{\mu}\equiv-\Delta^{\mu}_{\alpha}u_{\beta}\mathcal{T}^{\alpha\beta}_{\rm
matter}$. The magnetic field in the fluid frame is expressed as
$b^{\mu}=\frac{1}{2}\epsilon^{\mu\nu\kappa\lambda}u_{\nu}F_{\lambda\kappa}$,
where $\epsilon$ is the Levi-Civita pseudo-tensor, and $F$ is the
electromagnetic Maxwell tensor. The pressure tensor in a tetrad frame with
$e^{\mu}_{(0)}=u^{\mu}$ and $e^{\mu}_{(3)}=\hat{e}_{b^{\mu}}$, is diagonal,
with the perpendicular $P_{\perp}$ and parallel $P_{\parallel}$ fluid pressure
components: $diag(P_{\perp},P_{\perp},P_{\parallel})$. The pressure can then
be used together with the number density $n$ to calculate the effective
temperature: $k_{B}T_{\parallel}=P_{\parallel}/n$,
$k_{B}T_{\perp}=P_{\parallel}/n$, and $T=(T_{\parallel}+2T_{\perp})/3$.

Figure 6: Anisotropy and heat flux for ions in GRPIC initialized with
$\beta_{0}=10$ at $t=15r_{g}/c$. (a) Number density (left) and magnetic field
(right) fluctuations. (b) Temperature anisotropy $T_{\parallel}/T_{\perp}$
(left) and the ratio of non-ideal and ideal matter stress-energy tensor
components in the tetrad frame $\mathcal{T}^{(0)(3)}/\mathcal{T}^{(3)(3)}$,
which represents the ratio of parallel heat flux $q_{\parallel}$ and parallel
pressure $P_{\parallel}$. The inside of the event horizon is shown by a black
circle in the center of the image, the thick black line outlines the
ergosphere, and thin black lines represent the magnetic field lines. (c) Ion
probability density as a function of $\beta_{\parallel}=P_{\parallel}/P_{B}$
and temperature anisotropy $T_{\perp}/T_{\parallel}$. White dashed lines
correspond to thresholds for mirror (top) and firehose (bottom) instabilities.
A mirror instability develops close to the BH; a firehose instability develops
in the equatorial region. Figure 7: Parallel heat flux $q_{\parallel}$ as a
fraction of the heat flux corresponding to particles free streaming along the
magnetic field lines $P_{\parallel}v_{{\rm th}}$, for ions in GRPIC
initialized with $\beta_{0}=10$. The Figure is a zoom-in of the region close
to the BH ($r<5r_{g}$). On the left a time of $40r_{g}/c$ is shown, on the
right – $100r_{g}/c$.

Figure 6 shows ion quantities from a GRPIC simulation initialized with
$\beta_{0}=10$ but at an early time of $t=15r_{g}/c$, compare to Fig.3 in the
main text which shows times of $t=40r_{g}/c$ and $t=100r_{g}/c$. We show
number density $n/\langle n\rangle$ (a, left) and magnetic field
$\sqrt{B^{2}/\langle B^{2}\rangle}$ (a, right) fluctuations, temperature
anisotropy $T_{\parallel}/T_{\perp}$ (b, left), ratio of the two stress-energy
tensor components $\mathcal{T}^{(0)(3)}/\mathcal{T}^{(3)(3)}$ in the tetrad
frame (b, right), and ion probability density as a function of
$\beta_{\parallel}=P_{\parallel}/P_{B}$ and $T_{\perp}/T_{\parallel}$. White
dashed lines in Fig. 6c correspond to the boundaries of the growth rate of
mirror (top line) and firehose (bottom line) instabilities exceeding 10% of
the ion cyclotron frequency calculated using NHDS  for our case,
$m_{i}/m_{e}=1$. Build-up of the magnetic field near the BH due to the polar
inflow of the plasma leads to an associated increase in $P_{\perp}\propto B$.
Thus, most of the plasma near the BH is dominated by the
$T_{\parallel}/T_{\perp}<1$ region (b, left), and the mirror instability
dominates. Because of our limited separation between macroscopic and
microscopic scales, the plasma transiently significantly overshoots the mirror
stability boundary Riquelme2015 , as seen in Figure 6c. This effect is
negligible for realistic accretion flows. Another anisotropy-driven
electromagnetic instability develops at the equator (a, right), where
$T_{\parallel}/T_{\perp}>1$ (b, left). Here the magnetic field fluctuations
are not accompanied by plasma density fluctuations, which is a signature of a
firehose instability.

Figure 8: Anisotropy in GRPIC simulation with $m_{i}/m_{e}=3$ initialized with
$\beta_{0}=10$ at $t=40r_{g}/c$. (a) Temperature anisotropy for ions
$T_{\parallel,i}/T_{\perp,i}$ (left) and electrons
$T_{\parallel,e}/T_{\perp,e}$ (right). Electron (b) and ion (c) probability
density as a function of $\beta_{\parallel}=P_{\parallel}/P_{B}$ and
temperature anisotropy $T_{\perp}/T_{\parallel}$. White dashed lines
correspond to thresholds for mirror (top) and firehose (bottom) instabilities.

Figure 7 shows the ratio of the heat flux, $q_{\parallel}$, to the value
corresponding to the free streaming of particles along magnetic field lines,
$P_{\parallel}v_{{\rm th}}$, for ions in our GRPIC simulations with
$\beta_{0}=10$, at two different times, $t=40r_{g}/c$ and $t=100r_{g}/c$. At
the earlier time, the kinetic-scale mirror instability dominates in the region
close to the BH, suppressing the heat flux, $q_{\parallel}\sim
0.1P_{\parallel}v_{{\rm th}}$. At the later time, in the absence of effective
scattering due to kinetic-scale instabilities, the heat flux approaches the
free-streaming value in low-$\beta$ regions. Thus, we find that the heat flux
is a significant fraction of the free-streaming value, even higher than what
was found in previous non-ideal GRMHD simulations (Foucart2017, ).

Figure 8 demonstrates anisotropies measured at a time of $t=40r_{g}/c$ in a
GRPIC simulation with $m_{i}/m_{e}=3$ initialized with $\beta_{0}=10$. We show
temperature anisotropy for electrons $T_{\parallel,e}/T_{\perp,e}$ (a, left)
and ions $T_{\parallel,i}/T_{\perp,i}$ (a, right), electron (b) and ion (c)
probability density as a function of $\beta_{\parallel}=P_{\parallel}/P_{B}$
and $T_{\perp}/T_{\parallel}$. White dashed lines in (b) and (c) correspond to
the boundaries of the growth rate of mirror (top line) and firehose (bottom
line) instabilities exceeding 10% of the ion cyclotron frequency calculated
using NHDS  for the case of $m_{i}/m_{e}=3$. Pressure anisotropy of ions in
this simulation is similar to the one in the simulation with $m_{i}/m_{e}=1$
(Figure 3 in the main text). The probability density of ions (c) extends up to
the mirror instability threshold. Scattering of electrons on the ion-induced
fluctuations leads to their lower degree of the anisotropy. At larger values
of the mass ratio, we expect excitation of electron-scale waves, which may
lead to qualitatively different saturated values of the electron anisotropy.

## References

  * [1] Event Horizon Telescope Collaboration 2019, Astrophys. J. Lett., 875, L1
  * [2] Event Horizon Telescope Collaboration 2022, Astrophys. J. Lett., 930, L12
  * [3] Porth, O., Chatterjee, K., Narayan, R., et al. 2019, Astrophys. J. Supp. Ser., 243, 2
  * [4] Gold, R., Broderick, A. E., Younsi, Z., et al. 2020, Astrophys. J., 897, 2
  * [5] GRAVITY Collaboration 2018, A&A, 618, L10
  * [6] Dodds-Eden, K., Porquet, D., Trap, G., et al. 2009, The American Astronomical Society, 698, 1
  * [7] Yusef-Zadeh, F., Bushouse, H., Wardle, M., et al. 2009, Astrophys. J., 706, 348
  * [8] Markoff, S. 2005, Astrophys. J., 618, L103
  * [9] Broderick, A. E. & Loeb, A. 2006, MNRAS, 367, 3
  * [10] Dexter, J., Tchekhovskoy, A., Jiménez-Rosales, A., et al. 2020, MNRAS, 497, 4
  * [11] Ripperda, B. and Bacchini, F. & Philippov, A. A. 2020, Astrophys. J., 900, 2
  * [12] Ripperda, B., Liska, M., Chatterjee, K., et al. 2022, Astrophys. J. Lett., 924, 2
  * [13] Bransgrove, A., Ripperda, B. & Philippov, A. 2021, Phys. Rev. Lett., 127, 055101
  * [14] Sironi, L. & Spitkovsky, A. 2014, Astrophys. J. Lett., 783, 1
  * [15] Bhattacharjee, A., Huang, Y., Yang, H. & Rogers, B. 2009, Phys. Plasmas, 16, 11
  * [16] Kunz, M. W., Schekochihin, A. A. & Stone J. M. 2014, 112, 20
  * [17] Riquelme, M. A., Quataert, E., & Verscharen, D. 2015, Astrophys. J., 800, 27
  * [18] Parfrey, K., Philippov, A. & Cerutti, B. 2019, Phys. Rev. Lett., 122, 035101
  * [19] White, C. J., Stone, J. M. & Gammie, C.F. 2016, The American Astronomical Society, 225, 2
  * [20] Quataert, E. 2003, Astronomische Nachrichten Supplement, 324, 1
  * [21] Ressler, S. M., White, C. J., Quataert, E. & Stone, J.M. 2020, The American Astronomical Society, 896, 1
  * [22] See Supplemental Material for additional discussion of the numerical setup of GRPIC and GRMHD simulations, measurement of the pressure anisotropy in the GRPIC simulation with $m_{i}/m_{e}=3$ and details about the calculation of the plasma stress-energy tensor, pressure tensor and heat flux.
  * [23] Beskin, V. S., Istomin, Y. N. & Parev, V. I. 1992, Soviet Astronomy, 36, 642
  * [24] Hirotani, K. & Okamoto, I. 1998, Astrophys. J., 497, 2
  * [25] Crinquand, B., Cerutti, B., Philippov, et al. 2020, Phys. Rev. Lett., 124, 145101
  * [26] Chen,A.Y. & Yuan,Y. 2020, Astrophys. J., 895, 2
  * [27] Tchekhovskoy, A., Narayan, R. & McKinney, J. C. 2011, MNRAS, 418, 1
  * [28] Verscharen,D. & Chandran, B.D.G. 2018, Research Notes of the AAS, 2, 13
  * [29] In principle, pair plasma at moderate $\beta$ is unstable to ion-cyclotron instability at lower values of the pressure anisotropy. However, we observe mirror instability to dominate in the saturated phase, in agreement with local simulations [17].
  * [30] Guo, F., Li, H., Daughton, W., et al. 2014, Phys. Rev. Lett., 113, 155005
  * [31] Werner, G. R., Uzdensky, D. A., Cerutti, B., et al. 2016, Astrophys. J. Lett., 816, L8
  * [32] Cerutti, B., Philippov, A., Parfrey, K. & Spitkovsky, A. 2015, MNRAS, 448, 1
  * [33] Chandra, M., Gammie, C. F., Foucart, F., Quataert, E. 2015, The American Astronomical Society, 810, 2
  * [34] Foucart, F., Chandra, M., Gammie, C. F, et al. 2017, MNRAS, 470, 2


# A Suite of Fairness Datasets for Tabular Classification

Martin Hirzel IBM Research, USA  Michael Feffer Carnegie Mellon University,
Pittsburgh, PA, USA

###### Abstract

There have been many papers with algorithms for improving fairness of machine-
learning classifiers for tabular data. Unfortunately, most use only very few
datasets for their experimental evaluation. We introduce a suite of functions
for fetching 20 fairness datasets and providing associated fairness metadata.
Hopefully, these will lead to more rigorous experimental evaluations in future
fairness-aware machine learning research.

Table 1: Static information about the datasets. Column ‘origin’ specifies from
where the data is downloaded. Columns ‘#rows’ and ‘#cols’ give the shape of X.
Columns ‘any cat.’ and ‘any mis.’ indicate whether X has any categorical
columns and any missing values, respectively. Column ‘#labels’ shows the
number of unique values in y and ‘target name’ shows the name of y. Columns
‘favorable labels’ and ‘protected attributes’ are part of the fairness
metadata.
*For details see https://github.com/IBM/lale/blob/master/examples/demo_fairness_datasets.ipynb.

name origin #rows #cols any any #la- target favorable protected attributes
cat. mis. bels name labels (first) (second) ricci OpenML 118 5 yes no 2
promotion Promotion race tae OpenML 151 5 no no 3 class_attribute 3 whether_…*
heart_disease OpenML 303 13 no no 2 target 1 age student_math OpenML 395 32
yes no 2 g3_ge_10 1 sex age student_por OpenML 649 32 yes no 2 g3_ge_10 1 sex
age creditg OpenML 1,000 20 yes no 2 class good personal…* age titanic OpenML
1,309 13 yes yes 2 survived 1 sex us_crime OpenML 1,994 102 yes no 2
crimegt70pct 0 blackgt6pct compas_violent ProPublica 4,020 51 yes yes 2
two_year_recid 0 sex race nlsy OpenML 4,908 15 yes no 2 income96gt17 1 age
gender compas ProPublica 6,172 51 yes yes 2 two_year_recid 0 sex race
speeddating OpenML 8,378 122 yes yes 2 match 1 samerace importa…* nursery
OpenML 12,960 8 yes no 5 class spec_prior parents meps19 AHRQ 16,578 1,825 yes
no 2 UTILIZATION 1 RACE meps21 AHRQ 17,052 1,936 yes no 2 UTILIZATION 1 RACE
meps20 AHRQ 18,849 1,825 yes no 2 UTILIZATION 1 RACE law_school OpenML 20,800
11 yes no 2 ugpagt3 TRUE race1 default_credit OpenML 30,000 24 no no 2
default.pay…* 0 sex bank OpenML 45,211 16 yes no 2 class 1 age adult OpenML
48,842 14 yes yes 2 class >50K race sex

##  1 Introduction

Many people share the goal of making artificial intelligence fairer to those
affected by it. There is extensive debate about which fairness interventions
are appropriate and effective to achieve this goal. This debate should be
informed, at least in part, by rigorous experimental evaluation. Rigorous
experiments can help stakeholders make more informed choices among existing
fairness interventions, as well as help researchers invent better ones.
Unfortunately, most papers about fairness interventions evaluate them on at
most a handful of datasets. This is because historically, it was hard to find
and fetch datasets relevant to fairness, as well as associate them with
fairness metadata, such as favorable labels or protected attributes.

###  1.1 Related Work

OpenML [8] provides thousands of datasets ready for machine learning
experiments, but does not identify which of them are relevant to fairness and
does not provide fairness metadata. AIF360 [2] provides functions for fetching
8 fairness datasets along with metadata, but requires using a special class or
a multi-level pandas index. Quy et al. [7] describe 15 fairness datasets, but
do not provide code for fetching them, do not provide machine-readable
metadata for them, and some of their datasets are difficult to obtain. We
applaud OpenML, AIF360, and Quy et al. for getting most of the way towards a
suite of fairness datasets and build upon their work to take the last missing
step.

###  1.2 Contribution

This paper describes a suite of 20 Python functions to fetch 20 datasets along
with fairness metadata (see Table 1). It focuses on tabular data with
classification targets, which is the most well-studied setting. (Other
settings also have merit but are beyond the scope of this paper.) To make
these functions easy to use, they simply return data in pandas format [6]
along with fairness metadata in JSON format.

##  2 Functions for Fetching Datasets

⬇

1{

2 "favorable_labels": ["good"],

3 "protected_attributes": [

4 {

5 "feature": "personal_status",

6 "reference_group": ["male div/sep", "male mar/wid", "male single"]

7 },

8 {

9 "feature": "age",

10 "reference_group": [[26, 1000]]

11 }

12 ]

13}

Figure 1: Example fairness metadata for the creditg dataset. The list
"favorable_labels" contains values of y that indicate a favorable outcome. The
list "protected_attributes" contains sub-objects indicating the name "feature"
of the attribute of X, along with a list "reference_group" specifying which
values or ranges of that attribute indicate membership in the privileged
group.

Our suite of dataset fetching functions grew over time in an effort to gather
fairness datasets that are available for easy download with reasonable terms
and usage restrictions. Most of them are used in the literature, and where
possible, the fairness metadata returned by our functions emulates prior work.
Each of the 20 functions does three things: first download the data, second
minimally process the data, and third provide fairness metadata to go along
with the data.

#### Download the data

Our library distributes only functions for fetching data, but the data itself
is not part of the library. When you consume the data, it is your
responsibility to honor its licenses and terms of use. As shown in Table 1, 15
of the datasets are from OpenML [8], which uses a CC-BY license as explained
here: https://www.openml.org/terms. Our functions download these 15 OpenML
datasets. On the other hand, for the remaining 5 datasets from AHRQ and
ProPublica, our functions do not download the data, but instead print
instructions for downloading them manually. Once downloaded, our functions
return them from local disk. The data use information for the MEPS datasets
from AHRQ is at https://meps.ahrq.gov/data_stats/data_use.jsp, and the
ProPublica data is at https://github.com/propublica/compas-analysis. Our
functions are themselves part of the open-source Lale library [1], distributed
under an Apache license.

#### Minimally process the data

Our functions perform only limited preprocessing, because preprocessing
impacts fairness and can be difficult to invert. That said, some preprocessing
already happened at source before downloading, beyond our control. Where the
prediction target is not yet categorical, our functions discretize it. Where
necessary, our functions drop the feature column from which the discretized
target was derived. There are some other cases where our functions drop
additional feature columns because they are not useful. Finally, some feature
columns lack a meaningful name and our functions rename them, e.g. from "v1"
to "age". See the code for details.

#### Provide fairness metadata

Each of the functions returns a JSON object with fairness metadata. Figure 1
shows an example. The fairness metadata comprises a list of favorable labels
(i.e., favorable values in y) and a list of protected attributes (i.e., column
names in X). For each protected attribute, it gives a list of either ranges or
values that indicate membership in the privileged group. In practice, features
and labels relevant to fairness considerations are subject to interpretation
and should be determined through careful consultation with stakeholders.
Hence, we opted for a simple format that is easy to change.

#### How to use the functions

First install the Lale library [1] by doing pip install lale. Then you can
call the functions as illustrated by the following two lines of Python code
for the creditg dataset:

⬇

1import lale.lib.aif360

2X, y, fairness_info = lale.lib.aif360.fetch_creditg_df()

After this code, X and y contain the features and labels of the data,
represented as a pandas dataframe and series [6], and fairness_info contains
the metadata, represented as a JSON object as illustrated in Figure 1. At this
point, you can use your favorite library to split and preprocess the data,
make predictions, evaluate metrics, and perhaps mitigate bias. A popular
choice for many of these tasks would be the sklearn library [3]. While our
dataset fetching functions are part of the Lale library [1], you do not need
to use Lale to process their results. On the other hand, Lale contains
additional code that uses the metadata, including bias mitigators and fairness
metrics.

Figure 2: Example pipeline for the speeddating dataset. (The preprocessing
shown here is not part of the fetching function, but instead, is done on its
result.) Referring to Table 1, since ‘any cat.’ is true, there are separate
sub-pipelines with Lale Project operators for categorical and numerical
columns. Since ‘any mis.’ is true, the sklearn SimpleImputer operators fill in
missing values. And since ‘#cols’ is greater than 32, an sklearn SelectKBest
operator reduces it to 32.

##  3 Characterization of the Datasets

We already saw some static information characterizing the datasets in Table 1.
This section provides additional information based on dynamic experiments. Our
experiments feed the result from the dataset fetching function into a pipeline
of operators from Lale [1], sklearn [3], and XGBoost [4]. Exactly which
operators are included in the pipeline depends on whether the dataset has
categorical values and missing values and whether it has a large number of
columns. Figure 2 shows an example pipeline for a dataset where all three of
these are true; pipelines for other datasets are simpler. When you use the
datasets, it is up to you what type of pipeline to use, which may or may not
resemble our example.

Figure 3: Metrics characterizing the datasets. Subplot ‘data_ci’ shows the
class imbalance of the data based on binarizing y using the metadata. Subplot
‘data_di’ shows the symmetric disparate impact of the data. Subplots ‘xgb_di’,
‘xgb_eo’, and ‘xgb_ba’ show the symmetric disparate impact, equal opportunity
difference, and balanced accuracy of predictions from XGBoost.

Figure 3 shows the results. Subplot ‘data_ci’ shows class imbalance; it is the
ratio between the number of unfavorable and favorable outcomes, and higher
values mean the data is more balanced. Subplot ‘data_di’ shows the symmetric
disparate impact [5]; it is the ratio of the favorable rates for the
unprivileged and privileged groups. Higher disparate impact values mean the
data is more fair, with values under 0.8 usually considered unfair [5]. The
remaining three subplots show averages from 5-fold cross-validation
experiments with a popular and well-performing classifier, XGBoost [4], with
error bars showing standard deviations. Subplot ‘xgb_di’ shows that while bias
in the data does not always exactly equal bias in predictions of a classifier
trained on the data, the trends are similar across the 20 datasets. Subplot
‘xgb_eo’ shows equal opportunity difference, which is the difference of true
positive rates between the unprivileged and privileged groups, with zero
indicating perfect fairness. Subplot ‘xgb_ba’ shows balanced accuracy, which
is the average recall for the all classes, where higher values are better and
the best value is 1. Despite using 5-fold cross-validation, the classifier
overfit a couple of datasets with 100% balanced accuracy.

##  4 Conclusion

We hope our functions for fetching fairness datasets are useful and we welcome
contributions to their open-source code. Ideally, future papers with
experimental evaluations of fairness interventions will use at least 20, if
not more, datasets.

## References

  * Baudart et al., [2021] Baudart, G., Hirzel, M., Kate, K., Ram, P., Shinnar, A., and Tsay, J. (2021).  Pipeline combinators for gradual AutoML.  In Advances in Neural Information Processing Systems (NeurIPS).  https://proceedings.neurips.cc/paper/2021/file/a3b36cb25e2e0b93b5f334ffb4e4064e-Paper.pdf.
  * Bellamy et al., [2018] Bellamy, R. K. E., Dey, K., Hind, M., Hoffman, S. C., Houde, S., Kannan, K., Lohia, P., Martino, J., Mehta, S., Mojsilovic, A., Nagar, S., Ramamurthy, K. N., Richards, J., Saha, D., Sattigeri, P., Singh, M., Varshney, K. R., and Zhang, Y. (2018).  AI Fairness 360: An extensible toolkit for detecting, understanding, and mitigating unwanted algorithmic bias.  https://arxiv.org/abs/1810.01943.
  * Buitinck et al., [2013] Buitinck, L., Louppe, G., Blondel, M., Pedregosa, F., Mueller, A., Grisel, O., Niculae, V., Prettenhofer, P., Gramfort, A., Grobler, J., Layton, R., VanderPlas, J., Joly, A., Holt, B., and Varoquaux, G. (2013).  API design for machine learning software: Experiences from the scikit-learn project.  https://arxiv.org/abs/1309.0238.
  * Chen and Guestrin, [2016] Chen, T. and Guestrin, C. (2016).  XGBoost: A scalable tree boosting system.  In Conference on Knowledge Discovery and Data Mining (KDD), pages 785–794.  http://doi.acm.org/10.1145/2939672.2939785.
  * Feldman et al., [2015] Feldman, M., Friedler, S. A., Moeller, J., Scheidegger, C., and Venkatasubramanian, S. (2015).  Certifying and removing disparate impact.  In Conference on Knowledge Discovery and Data Mining (KDD), pages 259–268.  https://doi.org/10.1145/2783258.2783311.
  * McKinney, [2011] McKinney, W. (2011).  pandas: a foundational Python library for data analysis and statistics.  Workshop on Python for High Performance and Scientific Computing (PyHPC), pages 1–9.  https://www.dlr.de/sc/Portaldata/15/Resources/dokumente/pyhpc2011/submissions/pyhpc2011_submission_9.pdf.
  * Quy et al., [2022] Quy, T. L., Roy, A., Iosifidis, V., Zhang, W., and Ntoutsi, E. (2022).  A survey on datasets for fairness-aware machine learning.  https://arxiv.org/abs/2110.00530.
  * Vanschoren et al., [2014] Vanschoren, J., van Rijn, J. N., Bischl, B., and Torgo, L. (2014).  OpenML: Networked science in machine learning.  SIGKDD Explorations Newsletter, 15(2):49–60.  http://doi.acm.org/10.1145/2641190.2641198.