# Pseudogap effects in the strongly correlated regime of the two-dimensional
Fermi gas

S. Ramachandran1 S. Jensen2 Y. Alhassid1 1Center for Theoretical Physics,
Sloane Physics Laboratory, Yale University, New Haven, Connecticut 06520, USA
2Department of Physics, University of Illinois at Urbana-Champaign, Urbana,
Illinois 61801, USA

###### Abstract

The two-species Fermi gas with attractive short-range interactions in two
spatial dimensions provides a paradigmatic system for the understanding of
strongly correlated Fermi superfluids in two dimensions. It is known to
exhibit a BEC-BCS crossover as a function of $\ln(k_{F}a)$, where $a$ is the
scattering length, and to undergo a Berezinskii-Kosterlitz-Thouless superfluid
transition below a critical temperature $T_{c}$. However, the extent of a
pseudogap regime in the strongly correlated regime of $\ln(k_{F}a)\sim 1$, in
which pairing correlations persist above $T_{c}$, remains largely unexplored
with controlled theoretical methods. Here we use finite-temperature auxiliary-
field quantum Monte Carlo (AFMC) methods on discrete lattices in the canonical
ensemble formalism to calculate thermodynamical observables in the strongly
correlated regime. We extrapolate to continuous time and the continuum limit
to eliminate systematic errors and present results for particle numbers
ranging from $N=42$ to $N=162$. We estimate $T_{c}$ by a finite-size scaling
analysis, and observe clear pseudogap signatures above $T_{c}$ and below a
temperature $T^{*}$ in both the spin susceptibility and free-energy gap. We
also present results for the contact, a fundamental thermodynamic property of
quantum many-body systems with short-range interactions.

_Introduction_.— Cold atomic Fermi gases are of great interest in diverse
areas of physics in part because they provide a well-defined paradigm of
strongly correlated Fermi superfluids. They have been the subject of intensive
experimental and theoretical studies. Of particular interest is the two-
species uniform Fermi gas with attractive short-range interactions, whose
strength can be controlled experimentally through a Feshbach resonance.

The interaction strength is characterized by the two-particle $s$-wave
scattering length $a$. In three spatial dimensions (3D), this system makes a
crossover at low temperatures from a Bose-Einstein condensate (BEC) regime of
weakly interacting dimers for $(k_{F}a)^{-1}\to\infty$ to a Bardeen-Cooper-
Schrieffer (BCS) regime at $(k_{F}a)^{-1}\to-\infty$ ($k_{F}$ is the Fermi
wavenumber). The physics of the interacting Fermi gas in two spatial
dimensions (2D) differs qualitatively from that in 3D. Unlike the 3D case,
where a two-particle bound state is formed when the interaction is
sufficiently attractive, there is a bound state for arbitrarily weak
attractive interactions Randeria  _et al._ (1989); Brodsky  _et al._ (2006).
There is still a BEC to BCS crossover Bertaina and Giorgini (2011) as a
function of the scattering parameter $\eta=\ln(k_{F}a)$ Levinsen and Parish
(2015), with the BEC and BCS limits corresponding to $\eta\rightarrow-\infty$
and $\eta\rightarrow\infty$, respectively.

While both the 2D and 3D systems undergo a superfluid transition below a
critical temperature $T_{c}$, in 2D this phase transition does not have a non-
vanishing condensate fraction with off-diagonal long-range order as in the 3D
case, but instead exhibits a quasi-long-range order with algebraic decay of
correlations in the superfluid regime. This 2D superfluid transition is known
as a Berezinskii-Kosterlitz-Thouless (BKT) transition Berezinskii (1971,
1972); Kosterlitz and Thouless (1973).

The 2D BEC-BCS crossover and BKT transition have been intensively studied both
experimentally Fröhlich  _et al._ (2011); Vogt  _et al._ (2012); Murthy  _et
al._ (2015); Ries  _et al._ (2015); Boettcher  _et al._ (2016); Fenech  _et
al._ (2016); Toniolo  _et al._ (2017); Luciuk  _et al._ (2017); Hueck  _et
al._ (2018); Murthy  _et al._ (2018) and theoretically Watanabe  _et al._
(2013); Matsumoto and Ohashi (2014); Bauer  _et al._ (2014); Anderson and Drut
(2015); Marsiglio  _et al._ (2015); Shi  _et al._ (2015); Galea  _et al._
(2016); Vitali  _et al._ (2017); Madeira  _et al._ (2017); Schonenberg  _et
al._ (2017); Mulkerin  _et al._ (2018); Hu  _et al._ (2019); Wu  _et al._
(2020); Pascucci and Salasnich (2020); Zhao  _et al._ (2020); Zielinski  _et
al._ (2020); Mulkerin  _et al._ (2020a, b); Wang  _et al._ (2020); He  _et
al._ (2022). The strongly correlated regime, $\eta=\ln(k_{F}a)\sim 1$, of the
2D BEC-BCS crossover presents a particular theoretical challenge. In this
regime there is no controlled analytic approach, and controlled computational
methods provide the most reliable results. At zero temperature, the diffusion
Monte Carlo studies of Ref. Bertaina and Giorgini (2011) calculated the
energy, the pairing gap, and the contact, and Refs. Shi  _et al._ (2015);
Galea  _et al._ (2016); Vitali  _et al._ (2017); Zielinski  _et al._ (2020)
addressed pairing correlations. At finite temperature, lattice auxiliary-field
Monte Carlo (AFMC) methods were used in Ref. Anderson and Drut (2015) to
calculate the pressure, compressibility, and contact of the 2D crossover but
without a continuum extrapolation. The BKT critical temperature was recently
calculated in Ref. He  _et al._ (2022) as a function of $\eta$ using AFMC
methods on large lattices in the grand-canonical ensemble with a continuum
extrapolation.

An open problem is the extent of a pseudogap regime, in which signatures of
pairing correlations survive above the critical temperature $T_{c}$ for
supefluidity. Measurements of the spectral function of a harmonically trapped
gas Feld  _et al._ (2011) indicate a pairing gap at temperatures well above
the critical temperature in the strongly correlated regime. More recent radio-
frequency spectroscopy experiments in the normal phase revealed the presence
of an energy gap in the spectrum, which, at $\eta\sim 1$, far exceeds the two-
body binding energy Murthy  _et al._ (2018). In Ref. Bauer  _et al._ (2014),
the pseudogap regime was studied by calculating the single-particle spectral
function using the Luttinger-Ward self-consistent field theory approach. A
pronounced depression in the single-particle density of states was found in
the normal phase of the strongly correlated regime, whereas the same method
applied in 3D to the unitary Fermi gas showed a substantially reduced
pseudogap signature Haussmann  _et al._ (2009); Zwerger (2016); Pini  _et al._
(2019). While this self-consistent $T$-matrix method compares remarkably well
to AFMC for the unitary Fermi gas Wlazłowski  _et al._ (2013); Jensen  _et
al._ (2020a); Richie-Halford  _et al._ (2020); Rammelmüller  _et al._ (2021),
it is nevertheless an uncontrolled approximation and its accuracy in 2D has
not been established.

In this work, we calculate thermodynamic properties of the 2D Fermi gas across
the superfluid transition in its strongly correlated regime and, in
particular, explore pseudogap effects above $T_{c}$. We employ lattice
auxiliary-field quantum Monte Carlo (AFMC) methods in the canonical ensemble
Alhassid (2017); Jensen  _et al._ (2019, 2020a) and extrapolate to continuous
time and the continuum limit, thus eliminating any systematic errors and
yielding results that are accurate up to statistical errors for a given
particle number $N$. This approach has the advantage of being a well-
controlled computational method.

We estimate the critical temperature of the BKT transition using a finite-size
scaling analysis of the largest eigenvalue of the two-body density matrix. We
identify two pseudogap signatures above the critical temperature: (i) the
suppression of the spin susceptibility above $T_{c}$ and below a temperature
$T^{*}$ (also called the spin gap), and (ii) the increase in a model-
independent free-energy gap with decreasing temperature within the spin gap
regime. The calculation of the free-energy gap requires the use of canonical-
ensemble AFMC and particle-reprojection Alhassid  _et al._ (1999); Jensen  _et
al._ methods. We find that the pseudogap regime is broad at a coupling of
$\eta=1$ and its width decreases with the coupling $\eta$ as the system
approaches its BCS regime. We also calculate Tan’s contact, a fundamental
quantity for quantum many-body systems with short-range interactions that
describes the short-distance pair correlations of opposite spin particles Tan
(2008); Werner and Castin (2012). We find that the contact increases as the
temperature decreases within the pseudogap regime.

_Finite-temperature canonical ensemble AFMC_.— We briefly describe the
auxiliary-field quantum Monte Carlo (AFMC) method we use in this work; for
recent reviews see Refs. Gubernatis  _et al._ (2017); Alhassid (2017). We
consider a system of $N$ spin $1/2$ fermions interacting with a contact
interaction of strength $V_{0}$ on a finite area with periodic boundary
conditions. We discretize our system on a lattice of size $N_{L}^{2}$ with
$N_{L}$ points in each direction and a lattice spacing $\delta x$. The
corresponding lattice Hamiltonian is then

$\hat{H}=\sum_{\mathbf{k},s_{z}}\epsilon_{\mathbf{k}}\hat{a}^{\dagger}_{\mathbf{k},s_{z}}\hat{a}_{\mathbf{k},s_{z}}+g\sum_{\mathbf{x}_{i}}\hat{n}_{\mathbf{x}_{i},\uparrow}\hat{n}_{\mathbf{x}_{i},\downarrow}\;,$
(1)

where $\hat{a}^{\dagger}_{\bf{k},\sigma}$ and $\hat{a}_{\bf{k},\sigma}$ are
creation and annihilation operators for fermions with momentum ${\bf k}$ and
spin projection $s_{z}$ and
$\epsilon_{\mathbf{k}}=\frac{\hbar^{2}\mathbf{k}^{2}}{2m}$ for a quadratic
dispersion. Here $g=V_{0}/(\delta x)^{2}$ is the coupling strength, chosen to
reproduce the physical scattering length on a lattice, and
$\hat{n}_{\mathbf{x}_{i},s_{i}}$ is the density operator of the fermions at
position $\mathbf{x}_{i}$ and spin projection $s_{z}$. The first sum on the
r.h.s. of Eq. (1) is taken over the complete first Brillouin zone, and the
second sum is taken over all lattice sites $\mathbf{x}_{i}$.

We divide the inverse temperature $\beta$ into $N_{\tau}$ time slices of
length $\Delta\beta=\beta/N_{\tau}$, and apply a symmetric Trotter-Suzuki
decomposition to the imaginary-time propagator
$e^{-\beta\hat{H}}=\left(e^{-\frac{\Delta\beta}{2}\hat{H}_{0}}e^{-\Delta\beta\hat{V}}e^{-\frac{\Delta\beta}{2}\hat{H}_{0}}\right)^{N_{\tau}}+\mathcal{O}[\left(\Delta\beta\right)^{2}]$,
where $\hat{H}_{0}$ is the kinetic energy and $\hat{V}$ is the interaction
term in the Hamiltonian (1). Rewriting
$\hat{V}=g\sum_{\mathbf{x}_{i}}(\hat{n}_{\mathbf{x}_{i}}^{2}-\hat{n}_{\mathbf{x}_{i}})/2$,
we use a Hubbard-Stratonovich transformation to decouple the two-body
interaction introducing auxiliary fields
$\\{\mathbf{\sigma}_{\mathbf{x}_{i}}(\tau_{n})\\}$ at each lattice site
$\mathbf{x}_{i}$ and time slice $\tau_{n}$. For $g<0$, the thermal propagator
can then be written in the form

$e^{-\beta\hat{H}}=\int
D[\mathbf{\sigma}]G_{\mathbf{\sigma}}\hat{U}_{\sigma}+\mathcal{O}\left[(\Delta\beta)^{2}\right]$
(2)

where
$D[\mathbf{\sigma}]=\prod_{n=1}^{N_{\tau}}\prod_{\mathbf{x}_{i}}\left(\sqrt{\Delta\beta|g|/2\pi}\right)\textup{d}\sigma_{\mathbf{x}_{i}}(\tau_{n})$
is the integration measure,
$G_{\mathbf{\sigma}}=e^{-\Delta\beta|g|\sum_{i,n}\sigma^{2}_{\mathbf{x}_{i}}(\tau_{n})/2}$
is a Gaussian weight, and
$\hat{U}_{\sigma}=\prod_{n}e^{-\Delta\beta\hat{H}_{0}/2}e^{-\Delta\beta\hat{h}_{\sigma}(\tau_{n})}e^{-\Delta\beta\hat{H}_{0}/2}$
is the time-ordered one-body thermal propagator for a given auxiliary-field
configuration, where
$\hat{h}_{\sigma}(\tau_{n})=g\sum_{\mathbf{x}_{i}}\mathbf{\sigma}_{\mathbf{x}_{i}}(\tau_{n})\hat{n}_{\mathbf{x}_{i}}-g\hat{N}/2$
is the corresponding one-body Hamiltonian at any given time slice $\tau_{n}$.
Each integral over a continuous auxiliary field
$\mathbf{\sigma}_{\mathbf{x}_{i}}(\tau_{n})$ is discretized by a three-point
Gaussian quadrature Dean  _et al._ (1993).

In canonical-ensemble AFMC we are interested in the calculation of thermal
observables at fixed particle numbers $N_{\uparrow},N_{\downarrow}$

$\\!\\!\\!\langle\hat{O}\rangle_{N_{\uparrow},N_{\downarrow}}=\frac{\operatorname{Tr}_{N_{\uparrow},N_{\downarrow}}(\hat{O}e^{-\beta\hat{H}})}{\operatorname{Tr}_{N_{\uparrow},N_{\downarrow}}(e^{-\beta\hat{H}})}=\frac{\int
D[\mathbf{\sigma}]\langle\hat{O}\rangle_{\sigma}G_{\mathbf{\sigma}}\operatorname{Tr}_{N_{\uparrow},N_{\downarrow}}\hat{U}_{\sigma}}{\int
D[\mathbf{\sigma}]G_{\mathbf{\sigma}}\operatorname{Tr}_{N_{\uparrow},N_{\downarrow}}\hat{U}_{\sigma}}$
(3)

where
$\langle\hat{O}\rangle_{\sigma}=\frac{\operatorname{Tr}_{N_{\uparrow},N_{\downarrow}}(\hat{O}\hat{U}_{\sigma})}{\operatorname{Tr}_{N_{\uparrow},N_{\downarrow}}\hat{U}_{\sigma}}$
is the expectation value of $\hat{O}$ for a given field configuration. The
traces at fixed particle numbers are calculated by using projection operators
$\hat{P}_{N_{s_{z}}}$, i.e.,
$\operatorname{Tr}_{N_{\uparrow},N_{\downarrow}}X=\operatorname{Tr}(\hat{P}_{N_{\uparrow}}\hat{P}_{N_{\downarrow}}\hat{X})$.
For a finite number $N_{s}$ of single-particle states (for given spin
projection $s_{z}$), the particle-number projection can be represented exactly
by a discrete finite Fourier sum Ormand  _et al._ (1994)

$\hat{P}_{N_{s_{z}}}=\frac{e^{-\beta\mu
N_{s_{z}}}}{N_{s}}\sum_{m=1}^{N_{s}}e^{-i\varphi_{m}N_{s_{z}}}e^{\left(\beta\mu+i\varphi_{m}\right)\hat{N}_{s_{z}}}\;.$
(4)

where $\varphi_{m}=2\pi m/N_{s}$ are quadrature points, and $\mu$ is a real
chemical potential introduced to stabilize the Fourier sum. Using (4), the
fixed particle-number traces on the r.h.s. side of Eq. (3) reduce to the
calculation of the unrestricted traces involving the propagator
$e^{(\beta\mu+i\phi_{m})\hat{N}}\hat{U}_{\sigma}$. Since the latter is a one-
body propagator, these traces can be expressed in terms of the single-particle
representation matrix $\mathbf{U}_{\sigma}$ of the many-particle propagator
$\hat{U}_{\sigma}$. For example

$\operatorname{Tr}\left[e^{(\beta\mu+i\varphi_{m})\hat{N}}\hat{U}_{\sigma}\right]=\det\left[1+e^{(\beta\mu+i\varphi_{m})}\mathbf{U}_{\sigma}\right]\;.$
(5)

Figure 1: Spin susceptibility $\chi$ (top row) and free energy gap
$\Delta_{F}$ (bottom row) vs. temperature $T/T_{F}$ for couplings of
$\eta=1.0$ (left column), $\eta=1.3$ (middle column), and $\eta=1.8$ (right
column) for particle numbers $N=74,90$ and $114$. The black dashed line with
the gray band is $T_{c}$, and the dotted-dashed line is $T^{*}$. $\chi$ is
measured in units of the free gas zero-temperature susceptibility $\chi_{0}$
and $\Delta_{F}$ is measured in units of the free Fermi gas energy $E_{F}$.
The insets in the bottom row show the free energy gap on a logarithmic scale
(see text).

In canonical-ensemble AFMC, we sample auxiliary-field configurations according
to the positive-definite distribution
$G_{\sigma}|\operatorname{Tr}_{N_{\uparrow},N_{\downarrow}}\hat{U}_{\sigma}|$
and use them to estimate the thermal canonical expectation values of
observables following Eq. (3).

We take the continuous time limit $\Delta\beta\rightarrow 0$ by performing
calculations for several values of $\Delta\beta$, and extrapolation in
$(\Delta\beta)^{2}$. We then take the continuum limit by performing
calculations for a range of lattices containing as many as $29^{2}$ points,
and performing a linear extrapolation in the filling factor
$\nu=N/N_{L}^{2}\rightarrow 0$ (see the Supplemental Material for detailed
fits). The continuum limit requires large lattice calculations, which were
made possible through the use of two algorithms: (i) a stable diagonalization
method that reduces the $N_{s}^{4}$ scaling of the stabilization of canonical-
ensemble AFMC to $N_{s}^{3}$ Gilbreth and Alhassid (2015), and (ii) a
controlled truncation of the single-particle model space that further reduces
the $N_{s}^{3}$ scaling to essentially $N_{s}N^{2}$ for most parts of the
algorithm Gilbreth  _et al._ (2021), offering a dramatic improvement for $N\ll
N_{s}$. A similar method was implemented in Ref. He  _et al._ (2019).

_Results_.— We present canonical-ensemble AFMC results for three different
couplings in the strongly correlated regime of the 2D BEC-BCS crossover:
$\eta=\ln(k_{F}a)=1.0,1.3$, and $1.8$. For each of these couplings, we carried
out calculations for $N=42,74,90,114$, and $162$ particles.

We estimate the critical temperature $T_{c}$ of the transition to
superfluidity for a given $\eta$ using the phenomenological finite-size
scaling approach of Refs. Nightingale (1982); dos Santos and Sneddon (1981)
which is suitable for the BKT universality class. In this approach we scale
the largest eigenvalue $\lambda_{\rm max}$ of the two-body density matrix
$\langle
a^{\dagger}_{\mathbf{k}_{1},\uparrow}a^{\dagger}_{\mathbf{k}_{2},\downarrow}a_{\mathbf{k}_{3},\downarrow}a_{\mathbf{k}_{4},\uparrow}\rangle$.
For details, see the Supplemental Material sup .

We find $T_{c}=0.11(2)\,T_{F}$ at $\eta=1.0$, $T_{c}=0.09(2)\,T_{F}$ at
$\eta=1.3$, and $T_{c}=0.06(2)\,T_{F}$ at $\eta=1.8$. These estimates agree
within error bars with the results of Ref. He  _et al._ (2022) obtained by a
different method. They are lower than the experimental estimates of $T_{c}$ in
Ref. Ries  _et al._ (2015).

(i) Spin susceptibility: We calculate the spin susceptibility $\chi$ using

$\chi=\frac{\beta}{V}\langle(N_{\uparrow}-N_{\downarrow})^{2}\rangle\;.$ (6)

In the calculation of $\chi$, we project only on the total number of particles
$N=N_{\uparrow}+N_{\downarrow}$, as opposed to the two-species projection used
for other observables.

The spin susceptibility is suppressed by pairing correlations. Our AFMC
results for $\chi$ in units of $\chi_{0}$, the zero-temperature susceptibility
of the non-interacting Fermi gas, are shown in the top row of Fig. 1 as a
function of $T/T_{F}$. We observe the suppression of $\chi$ as the temperature
is decreased below $T_{c}$. However, we also see moderate suppression of
$\chi$ above $T_{c}$ but below a temperature scale $T^{*}$, known as a spin
gap. We define $T^{*}$ to be the temperature at which $\chi/\chi_{0}$ reaches
95% of its maximal value. This spin gap regime becomes narrower with
increasing $\eta$. We estimate $T^{*}=0.31\,T_{F},0.22\,T_{F}$ and
$0.12\,T_{F}$ for $\eta=1,1.3$ and $1.8$, respectively.

(ii) Free energy gap: We define the free energy gap by

$\Delta_{F}=[2F\left(N/2-1,N/2\right)-F\left(N/2-1,N/2-1\right)\\\
-F\left(N/2,N/2\right)]/2\;,$ (7)

where $F(N_{\uparrow},N_{\downarrow})$ is the free energy for the system at
$N_{\uparrow}$ spin-up particles and $N_{\downarrow}$ spin-down particles. To
calculate $\Delta_{F}$, we rewrite it as

$\mathbb{\bigtriangleup}_{F}=-k_{B}T\ln\left(\frac{Z_{N/2,N/2+1}}{Z_{N/2,N/2}}\frac{Z_{N/2,N/2+1}}{Z_{N/2+1,N/2+1}}\right)\;,$
(8)

where $Z_{N_{\uparrow},N_{\downarrow}}$ is the partition function for
$N_{\uparrow}$ spin-up particles and $N_{\downarrow}$ spin-down particles.
Each partition function ratio can then be calculated using particle-number
reprojection Alhassid  _et al._ (1999); Jensen  _et al._ , e.g.,

$\frac{Z_{N/2,N/2+1}}{Z_{N/2,N/2}}=\left\langle\frac{{\rm
Tr}_{N/2,N/2+1}U_{\sigma}}{{\rm Tr}_{N/2,N/2}U_{\sigma}}\right\rangle_{W}\;,$
(9)

where we have introduced the notation $\langle
X_{\sigma}\rangle_{W}\equiv{\int D[\sigma]W_{\sigma}X_{\sigma}/\int
D[\sigma]W_{\sigma}}$ with $W_{\sigma}=G_{\sigma}{\rm Tr}_{N/2,N/2}U_{\sigma}$
a positive definite weight function for $N_{\uparrow}=N_{\downarrow}=N/2$.

Figure 2: The contact $C$ (in units of $Nk_{F}^{2}$) vs. temperature for
couplings of $\eta=\ln(k_{F}a)=1.0,1.3$, and $1.8$ for particle numbers
$N=74,90$ and $114$. The dashed black line with the gray band represents
$T_{c}$, and the dotted-dashed line is $T^{*}$.

We show our results for $\Delta_{F}$ as a function of $T/T_{F}$ in the bottom
row of Fig. 1. Using a linear extrapolation at low temperatures, we determine
the zero-temperature energy staggering pairing gap Carlson  _et al._ (2003);
Gezerlis and Carlson (2008) $\Delta$ from the finite-temperature $\Delta_{F}$.
For additional details, see the Supplemental Material sup . A direct AFMC
calculation of $\Delta$ is challenging due to a sign problem introduced by the
spin imbalance.

We observe that the free energy gap $\Delta_{F}$ is suppressed above $T_{c}$.
To see more closely the behavior of $\Delta_{F}$ above $T_{c}$, we show
$\Delta_{F}$ in the insets using a logarithmic scale. In the spin gap regime
between $T_{c}$ and $T^{*}$ we observe an increase in $\Delta_{F}$ with
decreasing temperature due to pairing correlations. We identify this behavior
of $\Delta_{F}$ as a pseudogap signature that correlates with the suppression
of the spin susceptibility $\chi$.

(iii) Contact: The contact $C$ describes the short-range correlations between
particles of opposite spin and is defined by

$\int d^{2}R\
g^{(2)}_{\uparrow\downarrow}\left(\mathbf{R}+\frac{\mathbf{r}}{2},\mathbf{R}-\frac{\mathbf{r}}{2}\right)\underset{r\to
0}{\sim}\frac{C}{(2\pi)^{2}}\ln^{2}r\;,$ (10)

where
$g^{(2)}_{\uparrow\downarrow}\left(\mathbf{r}_{\uparrow},\mathbf{r}_{\downarrow}\right)=\langle\hat{n}_{\uparrow}\left(\mathbf{r}_{\uparrow}\right)\hat{n}_{\downarrow}\left(\mathbf{r}_{\downarrow}\right)\rangle$
is the two-body correlation function with
$\hat{n}_{\uparrow}\left(\mathbf{r}_{\uparrow}\right)$ and
$\hat{n}_{\downarrow}\left(\mathbf{r}_{\downarrow}\right)$ the density of spin
up and spin down particles, respectively. The contact is of key interest due
to Tan’s relations Tan (2008); Werner and Castin (2012), which relate the
contact to several different properties of the interacting Fermi gas. In
particular, the contact characterizes the tail of the momentum distribution
$n_{\sigma}(\mathbf{k})\sim C/k^{4}$ in the limit $k\rightarrow\infty$. It can
also be obtained from the derivative of the thermal energy with respect to the
coupling parameter

$\frac{\hbar^{2}C}{2\pi m}=\frac{dE}{d\ln a}\;.$ (11)

In the lattice simulations we calculate $C$ from the thermal expectation value
$\langle V\rangle$ of the potential energy Jensen  _et al._ (2020b)

$C=\frac{m^{2}}{\hbar^{4}}V_{0}\langle\hat{V}\rangle\;.$ (12)

We show our results for the contact (in units of $Nk_{F}^{2}$) as a function
of temperature for the three coupling values in Fig. 2. The contact increases
rapidly with lowering $T$ in the vicinity of $T_{c}$. This increase occurs
within a narrower temperature range as the coupling parameter $\eta$
increases. In general, the contact decreases with increasing $\eta$. In the
pseudogap regime between $T_{c}$ and $T^{*}$, we observe a monotonic increase
in the contact with decreasing $T$.

_Conclusion and outlook_.— We explored thermal properties of the two-species
spin-balanced Fermi gas with short-range interactions in two spatial
dimensions in the strongly correlated regime $\eta=\ln(k_{F}a)\sim 1$. We used
canonical-ensemble AFMC methods on discrete lattices and extrapolated our
results to continuous time and the continuum limit. We estimated the critical
temperature $T_{c}$ for the superfluid transition and calculated several
thermodynamic observables across the superfluid phase transition for several
values of $\eta$. We observed pseudogap signatures in a temperature regime
$T_{c}<T<T^{*}$, including a suppression of the spin susceptibility and an
increase in the free energy gap with decreasing temperature. Our AFMC results
for the spin susceptibility and free-energy gap are the first controlled
calculation of observables that probe the pseudogap regime in two spatial
dimensions and thus provide an accurate benchmark for experiment. We also
calculated Tan’s contact and found it to increase monotonically with
decreasing temperature in the pseudogap regime. In future AFMC studies, it
will be interesting to calculate the spectral function of the 2D strongly
interacting Fermi gas across the BEC-BCS crossover as a dynamical probe of the
pseudogap regime.

_Acknowledgements_.— This work was supported in part by the U.S. DOE grants
No. DE-SC0019521 and No. DE-SC0020177. The calculations used resources of the
National Energy Research Scientific Computing Center (NERSC), a U.S.
Department of Energy Office of Science User Facility operated under Contract
No. DE-AC02-05CH11231. We thank the Yale Center for Research Computing for
guidance and use of the research computing infrastructure.

## References

  * Randeria  _et al._ (1989) M. Randeria, J.-M. Duan, and L.-Y. Shieh, Phys. Rev. Lett. 62, 981 (1989).
  * Brodsky  _et al._ (2006) I. V. Brodsky, M. Y. Kagan, A. V. Klaptsov, R. Combescot, and X. Leyronas, Phys. Rev. A 73, 032724 (2006).
  * Bertaina and Giorgini (2011) G. Bertaina and S. Giorgini, Phys. Rev. Lett. 106, 110403 (2011).
  * Levinsen and Parish (2015) J. Levinsen and M. M. Parish, “Strongly interacting two-dimensional fermi gases,” in  _Annual Review of Cold Atoms and Molecules_ (World Scientific, 2015) Chap. 1, pp. 1–75.
  * Berezinskii (1971) V. L. Berezinskii, Sov. Phys. JETP 32, 493 (1971).
  * Berezinskii (1972) V. L. Berezinskii, Sov. Phys. JETP 34, 610 (1972).
  * Kosterlitz and Thouless (1973) J. M. Kosterlitz and D. J. Thouless, Journal of Physics C: Solid State Physics 6, 1181 (1973).
  * Fröhlich  _et al._ (2011) B. Fröhlich, M. Feld, E. Vogt, M. Koschorreck, W. Zwerger, and M. Köhl, Phys. Rev. Lett. 106, 105301 (2011).
  * Vogt  _et al._ (2012) E. Vogt, M. Feld, B. Fröhlich, D. Pertot, M. Koschorreck, and M. Köhl, Phys. Rev. Lett. 108, 070404 (2012).
  * Murthy  _et al._ (2015) P. A. Murthy, I. Boettcher, L. Bayha, M. Holzmann, D. Kedar, M. Neidig, M. G. Ries, A. N. Wenz, G. Zürn, and S. Jochim, Phys. Rev. Lett. 115, 010401 (2015).
  * Ries  _et al._ (2015) M. G. Ries, A. N. Wenz, G. Zürn, L. Bayha, I. Boettcher, D. Kedar, P. A. Murthy, M. Neidig, T. Lompe, and S. Jochim, Phys. Rev. Lett. 114, 230401 (2015).
  * Boettcher  _et al._ (2016) I. Boettcher, L. Bayha, D. Kedar, P. A. Murthy, M. Neidig, M. G. Ries, A. N. Wenz, G. Zürn, S. Jochim, and T. Enss, Phys. Rev. Lett. 116, 045303 (2016).
  * Fenech  _et al._ (2016) K. Fenech, P. Dyke, T. Peppler, M. G. Lingham, S. Hoinka, H. Hu, and C. J. Vale, Phys. Rev. Lett. 116, 045302 (2016).
  * Toniolo  _et al._ (2017) U. Toniolo, B. C. Mulkerin, C. J. Vale, X.-J. Liu, and H. Hu, Phys. Rev. A 96, 041604 (2017).
  * Luciuk  _et al._ (2017) C. Luciuk, S. Smale, F. Böttcher, H. Sharum, B. A. Olsen, S. Trotzky, T. Enss, and J. H. Thywissen, Phys. Rev. Lett. 118, 130405 (2017).
  * Hueck  _et al._ (2018) K. Hueck, N. Luick, L. Sobirey, J. Siegl, T. Lompe, and H. Moritz, Phys. Rev. Lett. 120, 060402 (2018).
  * Murthy  _et al._ (2018) P. A. Murthy, M. Neidig, R. Klemt, L. Bayha, I. Boettcher, T. Enss, M. Holten, G. Zürn, P. M. Preiss, and S. Jochim, Science 359, 452 (2018).
  * Watanabe  _et al._ (2013) R. Watanabe, S. Tsuchiya, and Y. Ohashi, Journal of Low Temperature Physics 171, 341 (2013).
  * Matsumoto and Ohashi (2014) M. Matsumoto and Y. Ohashi, Journal of Physics: Conference Series 568, 012012 (2014).
  * Bauer  _et al._ (2014) M. Bauer, M. M. Parish, and T. Enss, Phys. Rev. Lett. 112, 135302 (2014).
  * Anderson and Drut (2015) E. R. Anderson and J. E. Drut, Phys. Rev. Lett. 115, 115301 (2015).
  * Marsiglio  _et al._ (2015) F. Marsiglio, P. Pieri, A. Perali, F. Palestini, and G. C. Strinati, Phys. Rev. B 91, 054509 (2015).
  * Shi  _et al._ (2015) H. Shi, S. Chiesa, and S. Zhang, Phys. Rev. A 92, 033603 (2015).
  * Galea  _et al._ (2016) A. Galea, H. Dawkins, S. Gandolfi, and A. Gezerlis, Phys. Rev. A 93, 023602 (2016).
  * Vitali  _et al._ (2017) E. Vitali, H. Shi, M. Qin, and S. Zhang, Phys. Rev. A 96, 061601 (2017).
  * Madeira  _et al._ (2017) L. Madeira, S. Gandolfi, and K. E. Schmidt, Phys. Rev. A 95, 053603 (2017).
  * Schonenberg  _et al._ (2017) L. M. Schonenberg, P. C. Verpoort, and G. J. Conduit, Phys. Rev. A 96, 023619 (2017).
  * Mulkerin  _et al._ (2018) B. C. Mulkerin, X.-J. Liu, and H. Hu, Phys. Rev. A 97, 053612 (2018).
  * Hu  _et al._ (2019) H. Hu, B. C. Mulkerin, U. Toniolo, L. He, and X.-J. Liu, Phys. Rev. Lett. 122, 070401 (2019).
  * Wu  _et al._ (2020) F. Wu, J. Hu, L. He, X.-J. Liu, and H. Hu, Phys. Rev. A 101, 043607 (2020).
  * Pascucci and Salasnich (2020) F. Pascucci and L. Salasnich, Phys. Rev. A 102, 013325 (2020).
  * Zhao  _et al._ (2020) H. Zhao, X. Gao, W. Liang, P. Zou, and F. Yuan, New Journal of Physics 22, 093012 (2020).
  * Zielinski  _et al._ (2020) T. Zielinski, B. Ross, and A. Gezerlis, Phys. Rev. A 101, 033601 (2020).
  * Mulkerin  _et al._ (2020a) B. C. Mulkerin, H. Hu, and X.-J. Liu, Phys. Rev. A 101, 013605 (2020a).
  * Mulkerin  _et al._ (2020b) B. C. Mulkerin, X.-J. Liu, and H. Hu, Phys. Rev. A 102, 013313 (2020b).
  * Wang  _et al._ (2020) X. Wang, Q. Chen, and K. Levin, New Journal of Physics 22, 063050 (2020).
  * He  _et al._ (2022) Y.-Y. He, H. Shi, and S. Zhang, Phys. Rev. Lett. 129, 076403 (2022).
  * Feld  _et al._ (2011) M. Feld, B. Fröhlich, E. Vogt, M. Koschorreck, and M. Köhl, Nature 480, 75 (2011).
  * Haussmann  _et al._ (2009) R. Haussmann, M. Punk, and W. Zwerger, Phys. Rev. A 80, 063612 (2009).
  * Zwerger (2016) W. Zwerger, “Strongly interacting fermi gases,” in  _Proceedings of the International School of Physics “Enrico Fermi” - Course 191 “Quantum Matter at Ultralow Temperatures”_ , edited by M. Inguscio, W. Ketterle, S. Stringari, and G. Roati (IOS Press, Amsterdam, SIF Bologna, 2016) pp. 63–141.
  * Pini  _et al._ (2019) M. Pini, P. Pieri, and G. C. Strinati, Phys. Rev. B 99, 094502 (2019).
  * Wlazłowski  _et al._ (2013) G. Wlazłowski, P. Magierski, J. E. Drut, A. Bulgac, and K. J. Roche, Phys. Rev. Lett. 110, 090401 (2013).
  * Jensen  _et al._ (2020a) S. Jensen, C. N. Gilbreth, and Y. Alhassid, Phys. Rev. Lett. 124, 090604 (2020a).
  * Richie-Halford  _et al._ (2020) A. Richie-Halford, J. E. Drut, and A. Bulgac, Phys. Rev. Lett. 125, 060403 (2020).
  * Rammelmüller  _et al._ (2021) L. Rammelmüller, Y. Hou, J. E. Drut, and J. Braun, Phys. Rev. A 103, 043330 (2021).
  * Alhassid (2017) Y. Alhassid, in  _Emergent Phenomena in Atomic Nuclei from Large-Scale Modeling: a Symmetry-Guided Perspective_ , edited by K. D. Launey (World Scientific, 2017).
  * Jensen  _et al._ (2019) S. Jensen, C. N. Gilbreth, and Y. Alhassid, European Journal of Physics: Special Topics 227, 2241 (2019).
  * Alhassid  _et al._ (1999) Y. Alhassid, S. Liu, and H. Nakada, Phys. Rev. Lett. 83, 4265 (1999).
  * (49) S. Jensen, C. N. Gilbreth, and Y. Alhassid, to be published .
  * Tan (2008) S. Tan, Annals of Physics 323, 2952 (2008).
  * Werner and Castin (2012) F. Werner and Y. Castin, Phys. Rev. A 86, 013626 (2012).
  * Gubernatis  _et al._ (2017) J. E. Gubernatis, N. Kawashima, and P. Werner,  _Quantum Monte Carlo Methods_ (Cambridge University Press, 2017).
  * Dean  _et al._ (1993) D. Dean, S. Koonin, G. Lang, W. Ormand, and P. Radha, Physics Letters B 317, 275 (1993).
  * Ormand  _et al._ (1994) W. E. Ormand, D. J. Dean, C. W. Johnson, G. H. Lang, and S. E. Koonin, Phys. Rev. C 49, 1422 (1994).
  * Gilbreth and Alhassid (2015) C. Gilbreth and Y. Alhassid, Computer Physics Communications 188, 1 (2015).
  * Gilbreth  _et al._ (2021) C. N. Gilbreth, S. Jensen, and Y. Alhassid, Computer Physics Communications 264 (2021).
  * He  _et al._ (2019) Y.-Y. He, H. Shi, and S. Zhang, Phys. Rev. Lett. 123, 136402 (2019).
  * Nightingale (1982) P. Nightingale, Journal of Applied Physics 53, 7927 (1982).
  * dos Santos and Sneddon (1981) R. R. dos Santos and L. Sneddon, Phys. Rev. B 23, 3541 (1981).
  * (60) See the Supplemental Material accompanying this article.
  * Carlson  _et al._ (2003) J. Carlson, S.-Y. Chang, V. R. Pandharipande, and K. E. Schmidt, Phys. Rev. Lett. 91, 050401 (2003).
  * Gezerlis and Carlson (2008) A. Gezerlis and J. Carlson, Phys. Rev. C 77, 032801 (2008).
  * Jensen  _et al._ (2020b) S. Jensen, C. N. Gilbreth, and Y. Alhassid, Phys. Rev. Lett. 125, 043402 (2020b).

Supplemental Material: Pseudogap effects in the strongly interacting regime of
the two-dimensional Fermi gas

We discuss a number of technical details including the continuum
extrapolations and finite-size scaling. We present additional results for the
zero-temperature pairing gap, the thermal energy, and the single-particle
momentum distribution.

##  .1 Extrapolations

The lattice AFMC calculations are carried out by dividing $\beta$ into $N_{t}$
time slices of finite length $\Delta\beta$ and for discrete lattices, and it
is necessary to take the continuous time limit $\Delta\beta\to 0$ and the
continuum limit $\nu\to 0$ (where $\nu$ is the filling factor). In the
following we demonstrate how these limits are calculated.

###  .1.1 Continuous time limit

There are two systematic errors in $\Delta\beta$. One arises from the
symmetric Trotter-Suzuki decomposition of $e^{-\beta\hat{H}}$ which is of
order $(\Delta\beta)^{2}$, and a second arises from using a three-point
Gaussian quadrature formula to evaluate the integral over each auxiliary
field, which is accurate to order $(\Delta\beta)^{2}$, We eliminate these
systematic errors by performing a linear extrapolation in $(\Delta\beta)^{2}$
to find the limit $\Delta\beta\rightarrow 0$. We demonstrate these
extrapolations for the contact in Fig. 1. We note that sufficiently small
values of $\Delta\beta$ are necessary in order to reach the quadratic regime.

For the the maximal eigenvalue of the two-body density matrix and the thermal
energy, we find only a weak dependence on $\Delta\beta$ for small values of
$\Delta\beta$, in which case we perform a flat extrapolation (i.e., take an
average value).

Figure 1: The contact $C$ (in units of $Nk_{F}^{2}$) versus the dimensionless
parameter $(E_{F}\Delta\beta)^{2}$ for a coupling parameter of
$\ln(k_{F}a)=1.3$ and $N=114$ particles. The dashed lines are linear
extrapolations in $(\Delta\beta)^{2}$ to find the limit
$\Delta\beta\rightarrow 0$. Only the points in red below a cutoff of
$(\Delta\beta E_{F})^{2}<0.002$ are included in the fit, and the extrapolated
results for $\Delta\beta=0$ are shown in black solid circles.  Figure 2: The
$\Delta\beta=0$ contact $C$ (in units of $Nk_{F}^{2}$) versus the filling
factor $\nu$ for $\ln(k_{F}a)=1$ (left column), $\ln(k_{F}a)=1.3$ (middle
column), and $\ln(k_{F}a)=1.8$ (right column) for various temperatures and
particle numbers. The dashed lines are linear extrapolations in $\nu$ to
determine the continuum limit at $\nu=0$.

###  .1.2 Continuum extrapolations

Figure 3: The contact $C$ versus $T/T_{F}$ for different finite lattice sizes
ranging from $13^{2}$ ($\nu=0.44$) to $23^{2}$ ($\nu=0.14$). The extrapolated
continuum results at $\nu\rightarrow 0$ are shown by the black solid circles.
We notice a strong dependence of the contact on the filling factor.  Figure 4:
$\lambda_{\rm max}$ (see text) scaled by $N^{-7/8}$ as a function of
$\beta/\beta_{F}$ at the three coupling values. We observe the merging of
curves with varying particle number in the critical region, marked by the gray
bands with a black dashed line at the critical temperature $T_{c}$. The gray
bands are used to estimate the errors in $T_{c}$.

The continuum limit is obtained by taking the limit
$\nu=N/N_{L}^{2}\rightarrow 0$, and is achieved by considering increasing
lattice sizes. Thermal observables behave linearly in $\nu$ for small values
of $\nu$ and we use linear extrapolations to extract the values of observables
in the limit $\nu\rightarrow 0$.

Examples of the linear extrapolations in $\nu$ are shown in Fig. 2 for the
contact. We find that this observable has a relatively strong dependence on
$\nu$. In Fig. 3, we show the contact (in units of $Nk_{F}^{2}$) for $N=74$
and a coupling of $\eta=1$ as a function of $T/T_{F}$ for several lattice
sizes $N_{L}^{2}$. The extrapolated contact in the limit $\nu\rightarrow 0$
(black solid circles) show significant differences with the finite filling
factor results even on a qualitative level. In particular, the contact in the
continuum limit above the critical temperature decreases with increasing
temperature. In contrast, the spin susceptibility (not shown) exhibits only a
weak dependence on $\nu$.

##  .2 Finite-size scaling

We estimate the critical temperature $T_{c}$ of the transition to
superfluidity through a finite-size scaling analysis of the largest eigenvalue
$\lambda_{\rm max}$ of the two-body density matrix $\langle
a^{\dagger}_{\mathbf{k}_{1},\uparrow}a^{\dagger}_{\mathbf{k}_{2},\downarrow}a_{\mathbf{k}_{3},\downarrow}a_{\mathbf{k}_{4},\uparrow}\rangle$.
We use the phenomenological renormalization group analysis of Refs.
Nightingale1982S ; Santos1981S . In the ordered phase $\lambda_{\rm max}$
scales as Moreo1991S

$\lambda_{\rm max}=L^{2-\eta}f(L/\xi)\;,$ (13)

where $L$ is the linear size of the system, $\xi$ is the correlation length,
and the exponent $\eta=1/4$ at $T=T_{c}$ Kosterlitz1973S . For a BKT
transition, the correlation length $\xi\rightarrow\infty$ for $T\leq T_{c}$,
and we expect the scaled curves $L^{-7/4}\lambda_{\rm max}$ across different
values of $L$ to merge for $T\rightarrow T_{c}$ Paiva2004S  (rather than
intersecting in the manner typical of the three-dimensional superfluid phase
transition Nightingale1982S ). At constant density, the particle number
$N\propto L^{2}$, and we scale $\lambda_{\rm max}$ according to
$L^{-7/4}\lambda_{\rm max}\sim N^{-7/8}\lambda_{\rm max}$ for systems with
different particle number $N$.

We demonstrate this finite-size scaling analysis in Fig. 4, where
$N^{-7/8}\lambda_{\rm max}$ is shown as a function of $\beta/\beta_{F}$ for
several values of $N$ at the three coupling values. We observe the merging of
the different curves in the critical regime and use it to estimate $T_{c}$ and
its associated error. We find $T_{c}=0.11(2)\,T_{F}$ at $\eta=1.0$,
$T_{c}=0.09(2)\,T_{F}$ at $\eta=1.3$, and $T_{c}=0.06(2)\,T_{F}$ at
$\eta=1.8$.

## Zero-temperature pairing gap

Figure 5: $\Delta-|\epsilon_{b}|/2$ (see text) as a function of the coupling
parameter $\eta=\ln(k_{F}a)$. Our results, based on extrapolations of the AFMC
free energy gaps to $T=0$, are shown by the green, red and blue symbols for
$N=74,90$ and $114$, respectively. The solid squares, x’s and diamonds are,
respectively, the diffusion Monte Carlo results of Refs. Bertaina2011S ;
Galea2016S ; Vitali2017S .

At $T=0$ the free energy gap coincides with the energy staggering pairing gap.
The statistical errors on the free-energy gap are relatively small and at low
temperatures we can use a linear extrapolation to determine the energy gap
$\Delta$ at $T=0$. In Fig. 5 we compare our results for
$\Delta-|\epsilon_{b}/2|$ (here $\epsilon_{b}=-4\hbar^{2}/(ma^{2}e^{2\gamma})$
is the two-particle binding energy with $\gamma\approx 0.577$ being the Euler-
Mascheroni constant) for several values of the particle number $N$ with the
diffusion Monte Carlo results of Ref. Bertaina2011S  (solid squares), Ref.
Galea2016S  (x’s) and Ref. Vitali2017S  (solid diamonds).

## Thermal energy

Figure 6:  Thermal energy versus $T/T_{F}$ for particle numbers ranging from
42 to 162 and for coupling parameters of $\ln(k_{F}a)=1.0$, $\ln(k_{F}a)=1.3$,
and $\ln(k_{F}a)=1.8$ We notice that the energy below $T_{c}$ increases with
increasing $\ln(k_{F}a)$, in agreement with Ref. Bertaina2011S . Figure 7:
Left column: single-particle momentum distribution $n(k)$ for various
couplings and lattice sizes at $T=0.07\ T_{\textup{F}}$ and $N=90$. Right
column: as in the left column but showing $n(k)k^{4}$ versus $k$. The black
dashed line is the continuum limit contact obtained from $\langle V\rangle$.

In Fig. 6 we show our results for the thermal energy, defined as the
expectation value of the Hamiltonian, as a function of $T/T_{F}$ for
$N=42,74,90,114$ and $162$. The thermal energy increases linearly with
temperature above the critical temperature, and depends weakly on temperature
below the critical temperature. Our results are consistent with the ground-
state energies calculated in Ref. Bertaina2011S .

## Single-particle momentum distribution

Our AFMC results for the single-particle momentum distribution, $n(k)$, are
shown in the top row of Fig. 7 for $N=90$ and $T=0.07\,T_{F}$. Of particular
interest is the tail of the momentum distribution. According to Tan’s
universality relations Tan2008S ; Werner2012S , $n(k)\rightarrow C/k^{4}$ at
large $k$, where $C$ is the contact. In the bottom row of Fig. 7, we compare
$n(k)k^{4}$ with the value of $C$ extracted from the average potential energy
(dashed lines) and find good agreement for the larger lattices. It is
difficult to take the continuum limit extrapolation of the momentum
distribution tail because of edge effects. Thus, calculating the contact from
$\langle V\rangle$ is a more effective method of obtaining the continuum
results for the contact.

## References

  * (1) P. Nightingale, Journal of Applied Physics 53, 7927 (1982).
  * (2) R. R. dos Santos and L. Sneddon, Phys. Rev. B 23, 3541 (1981).
  * (3) A. Moreo and D. J. Scalapino, Phys. Rev. Lett. 66, 946 (1991).
  * (4) J. M. Kosterlitz and D. J. Thouless, Journal of Physics C: Solid State Physics 6, 1181 (1973).
  * (5) T. Paiva, R. R. dos Santos, R. T. Scalettar, and P. J. H. Denteneer, Phys. Rev. B 69, 184501 (2004).
  * (6) G. Bertaina and S. Giorgini, Phys. Rev. Lett. 106, 110403 (2011).
  * (7) A. Galea, H. Dawkins, S. Gandolfi, and A. Gezerlis, Phys. Rev. A 93, 023602 (2016).
  * (8) E. Vitali, H. Shi, M. Qin, and S. Zhang, Phys. Rev. A 96, 061601 (2017).
  * (9) S. Tan, Annals of Physics 323, 2952 (2008).
  * (10) F. Werner and Y. Castin, Phys. Rev. A 86, 013626 (2012).


# Towards objective and systematic evaluation of bias in medical imaging AI

Emma A.M. Stanley  Raissa Souza  Anthony Winder  Vedant Gulve  Kimberly Amador
Matthias Wilms111Shared last authorship. Nils D. Forkert222Shared last
authorship.

###### Abstract

Background. Artificial intelligence (AI) models trained using medical images
for clinical tasks often exhibit bias in the form of disparities in
performance between subgroups. Since not all sources of biases in real-world
medical imaging data are easily identifiable, it is challenging to
comprehensively assess how those biases are encoded in models, and how capable
bias mitigation methods are at ameliorating performance disparities. In this
article, we introduce a novel analysis framework for systematically and
objectively investigating the impact of biases in medical images on AI models.

Methods. We developed and tested a framework for conducting controlled in
silico trials to assess bias in medical imaging AI using a tool for generating
synthetic magnetic resonance images with known disease effects and sources of
bias. The feasibility of this framework is showcased by using three
counterfactual bias scenarios to measure the impact of simulated bias effects
on a convolutional neural network (CNN) classifier and the efficacy of three
bias mitigation strategies.

Findings. The analysis revealed that the simulated biases resulted in expected
subgroup performance disparities when the CNN was trained on the synthetic
datasets. Moreover, reweighing was identified as the most successful bias
mitigation strategy for this setup, and we demonstrated how explainable AI
methods can aid in investigating the manifestation of bias in the model using
this framework.

Interpretation. Developing fair AI models is a considerable challenge given
that many and often unknown sources of biases can be present in medical
imaging datasets. In this work, we present a novel methodology to objectively
study the impact of biases and mitigation strategies on deep learning
pipelines, which can support the development of clinical AI that is robust and
responsible.

###### keywords:

Artificial intelligence , medical imaging , algorithmic bias , synthetic data

[label1]organization=Biomedical Engineering Graduate Program, University of
Calgary,city=Calgary, country=Canada [label2]organization=Department of
Radiology, University of Calgary,city=Calgary, country=Canada
[label3]organization=Hotchkiss Brain Institute, University of
Calgary,city=Calgary, country=Canada [label4]organization=Alberta Children’s
Hospital Research Institute, University of Calgary,city=Calgary,
country=Canada [label5]organization=Department of Electronics and Electrical
Communication Engineering, Indian Institute of Technology,city=Kharagpur,
country=India [label6]organization=Department of Pediatrics, University of
Calgary,city=Calgary, country=Canada [label7]organization=Department of
Community Health Sciences, University of Calgary,city=Calgary, country=Canada
[label8]organization=Department of Clinical Neuroscience, University of
Calgary,city=Calgary, country=Canada [label9]organization=Department of
Electrical and Software Engineering, University of Calgary,city=Calgary,
country=Canada

##  1 Introduction

Artificial intelligence (AI) systems are often praised for having the
potential to serve as an objective decision-maker in health care. However,
recent studies have shown that AI can produce systematic disparities in
performance outcomes between subgroups; such models are considered to be
biased or unfair. This is of particular concern in the domain of AI for
medical imaging, where the first systems have been approved for assisting
clinicians in making diagnostic and treatment decisions for their patients
(Yearley et al., 2023). In this context, models may overfit to majority
populations in the training data or learn spurious correlations that lead to
poor generalization ability. For instance, various studies have found that
deep learning models trained on clinical tasks achieve better performance on
racial subgroups that represent the majority in the training data (Piçarra and
Glocker, 2023; Puyol-Antón et al., 2021; Seyyed-Kalantari et al., 2021). Other
studies have shown that models trained on medical images can encode explicit
information about image acquisition systems or sociodemographic
characteristics (Burns et al., 2023; Gichoya et al., 2022; Souza et al.,
2023). This suggests that models could learn these “bias attributes” rather
than representations of the disease of interest itself, which may lead to so-
called shortcut learning (Banerjee et al., 2023). Given these concerning
findings, a considerable amount of recent research has focused on applying and
developing strategies to mitigate subgroup performance disparities in clinical
models. Noteworthy examples include mitigating racial bias in cardiac
segmentation (Puyol-Antón et al., 2021), removing confounders in HIV diagnosis
from brain scans (Adeli et al., 2021), pruning to debias chest
x-ray(Marcinkevics et al., 2022) and dermatological models (Wu et al., 2022),
and harmonizing information related to magnetic resonance (MR) scanners
encoded in models (Dinsdale et al., 2021).

Moreover, recent studies have started investigating how biases from medical
imaging data are encoded in models (Brown et al., 2023; Glocker et al., 2023;
Jones et al., 2023). An improved understanding about this holds great value
since it could facilitate the development of models that are inherently robust
to such biases. However, it is difficult to comprehensively and objectively
assess this phenomenon because medical imaging datasets contain numerous
biases or spurious correlations that deep learning models may use, but that we
as humans do not know about, cannot fully understand, or have no ability to
disentangle (Banerjee et al., 2023). The same is true for studies that
evaluate or propose bias mitigation methods – while it is valuable to know
that certain techniques are successful in reducing performance disparities
between subgroups on certain clinical tasks, it remains unknown if such
strategies are correcting for the true sources of bias in images, and if they
are also effective when applied to data from populations not represented in
the datasets used to evaluate them. Therefore, to objectively and
systematically study how biases manifest in deep learning models, and
confidently evaluate the robustness of bias mitigation strategies, it is
imperative to understand and control the biases in a medical imaging dataset
in the first place.

In this study, we introduce and utilize an analysis framework for conducting
such systematic investigations of the impact of bias on medical imaging AI. As
a key component, we use the recently proposed Simulated Bias in Artificial
Medical Images (SimBA) (Stanley et al., 2023) tool for generating structural
neuroimaging datasets with user-defined biases. We exploit the controlled and
customizable nature of this data-generating mechanism to generate
counterfactual subject-paired datasets with biases in different spatial
locations of the human brain, or without the presence of bias at all. Thus,
with identical model training schemata, we can conduct in silico trials
exploring bias in AI, in which we systematically study the impact of different
bias manifestations on the same deep learning pipeline. Moreover, we use these
synthetic datasets to assess how well existing bias mitigation techniques
compensate for the simulated biases by reducing subgroup performance
disparities.

##  2 Methods

###  2.1 Generation and study of dataset bias

As the basis of this work, we use and extend our recently proposed SimBA
(Stanley et al., 2023) tool to generate synthetic T1-weighted brain MR images.
Briefly described, morphological variations representing disease, bias, and
subject effects are sampled from principal component analysis (PCA)-based
generative models of non-linear deformations and are applied to a template
image that represents an average brain morphology. In this setup, disease and
bias effects are represented by spatially localized deformations, whereas
distinct subjects are represented by variation in global brain morphology.
This procedure of applying morphological effects to the template image is
repeated numerous times to generate a large dataset, where the user specifies
the images belonging to each target (disease) class and the images containing
bias effects. The whole process is user-defined such that the dataset contains
a specified number of images, number of target classes, and proportions of
images in each class with or without bias effects, which enables a controlled
setup for investigating the impact of bias on deep learning models.

In this work, we extend this basic methodology to pre-define stratified
distributions of subject and disease effect magnitude. Simply put, this
ensures that the subject and disease effect magnitude distributions for each
class and bias group are similar, which helps to prevent these distributions
from introducing an additional, unwanted source of bias in deep learning
models. This significant extension of SimBA is, therefore, crucial for the aim
of this work. Further details are depicted in Supplementary Figure 4.

The datasets generated for this study used the SRI24 atlas (Rohlfing et al.,
2010) with LPBA40 labels (Shattuck et al., 2008) as the brain template, with
dimensions of 173×211×155 voxels and an isotropic resolution of 1 mm3. The
morphological effects added to this template were derived from PCA models fit
to stationary velocity fields of T1-weighted MRIs (Log-Euclidean framework
(Arsigny et al., 2006)), which were obtained via non-linear diffeomorphic
registration of 50 images from the IXI (brain development.org, ) database to
the SRI24 atlas, as described in Stanley et al. (2023). Each full dataset
generated for the experiments in this work is comprised of 2002 three-
dimensional images of distinct simulated subjects, with two target classes
(“disease”, N=1000, and “non-disease”, N=1002), and two bias groups (“bias”,
containing an additional morphological effect, and “non-bias”, which does not
contain such an effect). 70% of the images in the disease class were defined
to belong to the bias group, while 30% of the images in the non-disease class
were defined to belong to the bias group. With this unequal representation of
the bias group in each class, it is expected that the deep learning model uses
the presence of the bias effect as a shortcut for predicting the disease
class, resulting in subgroup performance disparities.

Supplementary Table 2 provides the detailed dataset generation parameters. We
defined the disease region, which contains the effects that are meant to act
as the prediction target for the deep learning model, as the left insular
cortex. Next, we generated datasets for three bias scenarios: 1) No Bias
Effect, 2) “Near” Bias Effect (left putamen, adjacent to disease region), and
3) “Far” Bias Effect (right postcentral gyrus, opposite hemisphere as the
disease region and on a different axial slice). These three bias scenarios are
“paired” such that the images in each dataset contain the same subject effects
and disease effects – in other words, each dataset acts as a counterfactual,
where the same generated subjects have either no bias effect, or a bias effect
in the near or far regions, respectively. Thus, evaluation of a model with
identical weight initialization and deterministic training allows for direct
comparison of the impact of each bias scenario on the deep learning pipeline.
We applied this training schema for prediction of the target class from the
imaging data using a convolutional neural network (CNN) (implementation
details in A). CNNs automatically extract relevant features from images,
making them a commonly used model architecture in medical image analysis where
subtle details or patterns can be indicative of diseases or conditions.

###  2.2 Bias mitigation

To further showcase the utility of this framework, we quantitatively evaluated
the efficacy of three conceptually different and representative bias
mitigation strategies from the literature for each bias scenario. First, we
applied reweighing, a pre-processing strategy proposed by Calders et al.
(2009) that adds weights to each image during training to balance the
contribution of bias groups in each target class. Second, we applied the in-
processing adversarial bias “unlearning” strategy described in Dinsdale et al.
(2021) which encourages the feature encoder to unlearn information required to
predict the bias group label, while keeping disease classification performance
high. Finally, we fine-tuned bias group models, following Puyol-Antón et al.
(2021) by pretraining a model on the full dataset, and then training separate
models for each bias group until convergence. Implementation details of the
bias mitigation strategies are provided in A.

###  2.3 Evaluation

Model accuracy was used as a metric to determine efficacy in achieving the
target prediction task of disease classification, where the disease class is a
positive prediction, and the non-disease class is a negative prediction.
Subgroup performance disparities, or the difference in true positive rate
($\Delta$TPR) and false positive rate ($\Delta$FPR) between the bias group and
the non-bias group, were the primary outcome metrics for measuring bias
(Equations 1, 2):

$\Delta TPR=TPR_{bias\;group}-TPR_{non-bias\;group}$ (1) $\Delta
FPR=FPR_{bias\;group}-FPR_{non-bias\;group}$ (2)

Standard deviations across five model seeds are reported as error values.

Exploiting the benefits of the counterfactual setup, we measured performance
disparities relative to the baseline (No Bias) scenario (Equations 3, 4),
since there are minimal but non-zero performance disparities due to
differences in the underlying disease and subject effect distributions. Thus,
any measured relative performance disparity between bias groups can be
attributed purely to the presence of the morphological bias effects we wish to
study.

$Relative\;\Delta TPR=\Delta TPR-\Delta TPR_{No\;Bias}$ (3) $Relative\;\Delta
FPR=\Delta FPR-\Delta FPR_{No\;Bias}$ (4)

###  2.4 Explainability

SmoothGrad (Smilkov et al., 2017) saliency maps were used to explore the
extent to which the bias regions influenced the model in predicting the target
class. Final saliency maps were generated from an average of 10 images from
each bias group that were correctly classified as belonging to either the
disease or non-disease class. Quantitative explainability is reported using
weighted saliency scores (Stanley et al., 2022), which represents the relative
intensity of salient voxels in a given brain region across the whole saliency
map.

###  2.5 Statistics

We performed a two-way repeated measures ANOVA on $\Delta$TPR and $\Delta$FPR
with bias scenario (No Bias, Near Bias, Far Bias) and bias mitigation method
(naïve, reweighing, unlearning, group models) as the factors. Following
significant effects, we performed Shapiro-Wilk tests to check for normality
and two-tailed comparison tests (paired t-tests if normality test was passed,
otherwise Wilcoxon matched-pairs signed rank test) with Bonferroni correction
to compare between specific bias scenarios and mitigation methods
(Supplementary Tables 4 and 5). The corrected significance level was
$\alpha$=0.005. Analysis was performed in GraphPad Prism 10.0.2.

##  3 Results

###  3.1 Naïve model

All models trained using data from the three bias scenarios demonstrate high
overall classification accuracy, with 84·81±0·62%, 87·68±0·20%, and
86·77±0·71% for the No Bias, Near Bias, and Far Bias datasets, respectively.
The presence of the bias effect resulted in significant performance
disparities compared to the No Bias scenario, with a relative $\Delta$TPR of
22·12±2·72% (p$<$0·0001) and $\Delta$FPR of 20·32±2·81% (p$<$0·0001) in the
Near Bias scenario, and a relative $\Delta$TPR of 17·35±2·78% (p=0·0002) and
relative $\Delta$FPR of 11·88±6·22% (p=0·0008) in the Far Bias scenario. These
results indicate that the models likely used the bias effect as a shortcut for
predicting the (positive) disease class, since the disease class contains a
higher proportion of the bias group compared to the non-disease class.
However, the models may have used this shortcut learning to a lesser extent
when the bias region is farther from the disease region, as evidenced by the
lower relative $\Delta$FPR and $\Delta$TPR in the Far Bias scenario compared
to the Near Bias scenario.

Table 1: Performance of each bias mitigation strategy for each bias scenario. Errors are reported as standard deviation across five model seeds. Naïve=No bias mitigation, RW=reweighing, UL=unlearning, GM=group models, TPR=True positive rate, FPR=False positive rate. | Aggregate | Bias Group | Non-Bias Group
---|---|---|---
| Naïve | RW | UL | GM | Naïve | RW | UL | GM | Naïve | RW | UL | GM
| Accuracy
---
(%)
|  |  |  |  |  |  |  |  |  |  |
No Bias |  | 84.81
---
± 0.62
| 85.01
---
± 0.39
| 79.28
---
± 2.76
- |  | 85.24
---
± 0.61
| 85.08
---
± 0.64
| 79.84
---
± 6.66
| 88.23
---
± 0.34
| 84.38
---
± 1.18
| 84.94
---
± 0.44
| 78.73
---
± 10.47
| 87.01
---
± 0.78
Near Bias |  | 87.68
---
± 0.20
| 84.33
---
± 0.39
| 86.65
---
± 0.72
- |  | 88.21
---
± 0.39
| 84.68
---
± 0.29
| 87.34
---
± 1.05
| 88.06
---
± 0.88
| 87.15
---
± 0.75
| 83.98
---
± 0.77
| 85.98
---
± 1.07
| 87.41
---
± 0.83
Far Bias |  | 86.77
---
± 0.71
| 84.89
---
± 0.52
| 85.29
---
± 0.99
- |  | 87.10
---
± 0.99
| 85.32
---
± 0.36
| 85.48
---
± 2.61
| 87.90
---
± 0.40
| 86.45
---
± 0.69
| 84.46
---
± 0.75
| 85.10
---
± 2.59
| 87.17
---
± 0.71
| TPR
---
(%)
|  |  |  |  |  |  |  |  |  |  |
No Bias |  | 85.73
---
± 1.01
| 85.00
---
± 1.67
| 80.24
---
± 17.02
- |  | 85.44
---
± 0.91
| 84.78
---
± 1.65
| 79.71
---
± 20.57
| 93.44
---
± 1.54
| 86.47
---
± 1.61
| 85.59
---
± 1.92
| 81.18
---
± 17.96
| 70.59
---
± 6.15
Near Bias |  | 88.31
---
± 0.74
| 84.44
---
± 1.44
| 89.60
---
± 4.34
- |  | 93.89
---
± 0.79
| 84.33
---
± 1.44
| 92.67
---
± 3.69
| 92.67
---
± 2.56
| 73.53
---
± 2.08
| 84.71
---
± 1.68
| 81.4
---
± 6.29
| 72.06
---
± 1.80
Far Bias |  | 85.08
---
± 2.53
| 85.89
---
± 0.29
| 85.24
---
± 6.48
- |  | 89.56
---
± 3.10
| 85.44
---
± 0.25
| 86.22
---
± 5.91
| 92.67
---
± 0.72
| 73.24
---
± 1.92
| 87.06
---
± 0.66
| 82.65
---
± 8.54
| 72.94
---
± 3.83
| FPR
---
(%)
|  |  |  |  |  |  |  |  |  |  |
No Bias |  | 16.10
---
± 1.28
| 14.98
---
± 1.40
| 21.67
---
± 20.70
- |  | 15.28
---
± 0.81
| 14.12
---
± 3.05
| 20.29
---
± 20.57
| 25.29
---
± 2.87
| 16.39
---
± 1.64
| 15.30
---
± 0.86
| 22.19
---
± 20.76
| 6.89
---
± 1.26
Near Bias |  | 12.95
---
± 0.95
| 15.78
---
± 1.53
| 16.25
---
± 3.73
- |  | 26.84
---
± 2.21
| 14.41
---
± 3.35
| 26.76
---
± 6.10
| 24.12
---
± 5.46
| 7.79
---
± 1.63
| 16.28
---
± 1.41
| 12.35
---
± 3.01
| 6.89
---
± 1.13
Far Bias |  | 11.55
---
± 1.67
| 16.10
---
± 1.21
| 14.66
---
± 6.13
- |  | 19.41
---
± 4.81
| 15.00
---
± 1.61
| 16.47
---
± 6.44
| 24.71
---
± 1.23
| 8.63
---
± 1.18
| 16.50
---
± 1.18
| 13.99
---
± 6.14
| 7.54
---
± 1.25

Saliency maps highlighted both the disease and the corresponding bias region
in each of the bias scenarios (Figure 1). Quantitatively, weighted saliency
scores demonstrated a higher intensity saliency within the regions affected by
the bias effect for both the Near and Far Bias scenarios relative to the No
Bias baseline (Figure 3).

Figure 1: Average SmoothGrad saliency maps of correctly classified “subjects”
for the bias group in the disease class, for each bias scenario in the naïve
models.

###  3.2 Bias mitigated models

Each of the bias mitigated models demonstrated high overall classification
accuracy, comparable to that of the naïve model (Table 1). It was also
observed that all bias mitigation strategies reduce TPR and FPR disparities
between the bias groups, albeit to different degrees (Figure 2).

The reweighing strategy almost perfectly mitigated any measured performance
disparity between the bias groups, to respective values of relative
$\Delta$TPR and relative $\Delta$FPR of 0·44±1·80% and –0·69±1·45% for the
Near Bias scenario, and –0·80±0·64% and –0·32±2·13% for the Far Bias scenario.

The unlearning strategy also reduced the performance disparities between
groups, although to a lesser degree. More precisely, the relative $\Delta$TPR
was reduced to 12·22±3·04% for the Near Bias scenario and 4·60±4·20% for the
Far Bias scenario. The relative $\Delta$FPR was reduced to 15·51±4·87% and
3·58±3·01% for the Near and Far Bias scenarios, respectively. Notably,
unlearning reduced performance disparities to a greater degree for the Far
Bias scenario compared to the Near Bias scenario. For all bias scenarios, the
ability of the model to predict the bias label was reduced to chance level
(Supplementary Table 6). Performing unlearning with the dataset for the No
Bias scenario resulted in lower and more variable model performance across all
metrics (Table 1).

Figure 2: Relative true (left) and false (right) positive rate disparities
between bias groups relative to the No Bias baseline for each bias mitigation
strategy and bias scenario across five model seeds. Markers indicate
significant differences between the naïve model performance ($\Delta$TPR,
$\Delta$FPR) and each bias mitigation strategy (two-tailed paired t-test with
Bonferroni correction, ***p$<$0·001).

Training separate models for each bias group also reduced the average
$\Delta$TPR and $\Delta$FPR to near-zero, relative to the group models trained
on the No Bias dataset. The relative $\Delta$TPR and $\Delta$FPR were
–2·25±6·50% and –1·47±6·45% for the Near Bias scenario, and –3·13±6·62% and
–1·54±1·59% for the Far Bias scenario. However, it was observed that the group
model for the No Bias baseline demonstrated significant subgroup performance
disparities compared to the No Bias naïve model ($\Delta$TPR=22·86±5·52%,
p=0·0007 and $\Delta$FPR=18·70±2·42%, p=0·0001). Although the No Bias dataset
did not have bias effects added, each bias group model was trained on a
different number of images from each disease class (Supplementary Table 3),
resulting in a typical class imbalance problem. Thus, while the addition of
the near and far bias effects did not lead to any further performance
disparities into the group models, we instead observed disparities implicitly
introduced due to the different levels of bias group representation in each
target class (Supplementary Figure 5).

Weighted saliency score results were in line with the efficacy of each bias
mitigation strategy (Figure 3). For reweighing and group models, performance
disparities were reduced to near-zero, and the corresponding saliency scores
were similar to the No Bias baseline. The higher magnitude of performance
disparities that remained after unlearning was reflected in the higher bias
region saliency scores for this strategy.

Figure 3: Weighted saliency scores calculated from average saliency maps over
five model seeds for the bias group in the disease class for the brain regions
most affected by “near” bias effect (left) and “far” bias effect (right). The
horizontal line represents the weighted saliency score average ± standard
deviation in the corresponding region for the No Bias naïve model.

##  4 Discussion

Overall, the results of this study highlight the utility and benefits of using
the proposed analysis framework for performing a systematic and controlled
evaluation of bias in medical imaging AI. We generated counterfactual
neuroimaging datasets for three bias scenarios using an extended version of
the SimBA tool and simulated in silico trials of each dataset with identical
and deterministic deep learning pipelines. Since the exact bias effects that
we aim to study in this setup are known, and it is also known how models
perform on the counterfactual dataset without bias, the naïve models provide a
strong baseline for objectively interpreting results of subsequent analyses.
Thus, in the evaluation of each bias mitigation strategy, we can better
understand the situations these mitigation strategies are effective in, or the
reasons they may fail.

The results showed that reweighing completely mitigated performance
disparities caused by the single source of bias we introduced. Thus, several
real-world evaluations of reweighing (Larrazabal et al., 2020; Meissen et al.,
2023; Puyol-Antón et al., 2021) may have been successful as a result of being
applied to subgroups associated with the predominant sources of bias in the
data. On the other hand, systematically introducing more complexity (e.g.,
additional sources of bias) to these synthetic datasets may help to improve
our understanding why reweighing fails to alleviate performance disparities in
various other real-world datasets (Ioannou et al., 2022; Meissen et al., 2023;
Weng et al., 2023; Zhang et al., 2022). For the unlearning strategy, it was
found that the model’s ability to predict the bias label was reduced to chance
level, but this did not effectively mitigate subgroup performance disparities.
This may be due to the bias features and the features necessary for predicting
the target class being too similar to be fully disentangled. Understanding why
this occurs likely requires deeper investigation into the model’s internal
representations of disease and bias effects, which would be a straightforward
task with these synthetic datasets since there are only two user-defined
effects that the model can use for the prediction task. Furthermore, given the
simulated counterfactual No Bias scenario, it can be seen that unlearning
applied to this dataset even harms the model, resulting in lower performance
on average and higher variability across model seeds. This finding is
especially important considering that potential bias attributes in real-world
datasets, such as scanner models or self-reported demographics, could be
inaccurately labeled. Further work could be performed using this framework to
analyze how mislabeling of bias attributes affects the result of this and
other mitigation strategies. Moreover, the No Bias counterfactual datasets
made it clear that performance disparities in the group models were not caused
by the bias effects, but rather by the target class imbalance within each bias
group. If this approach was used to de-bias a classification task on real
data, the cause of these performance disparities may have been unclear since
it could have been related to other subgroup attributes unaccounted for.
Finally, we observed that the saliency scores align with performance of the
model for each bias mitigation strategy. With knowledge of the exact regions
that should be relevant for model decision-making, and of how a particular
explainability method manifests in the absence of the bias effect, this
framework facilitates a highly objective interpretation of explainable AI in
medical imaging. Further work could use this setup for studying the extent to
which different explainability methods highlight regions associated with
predictive targets or bias effects, or developing new methods that are better
at identifying salient regions of medical images.

Although the SimBA tool enables the generation of increasingly complex dataset
bias scenarios, real-world medical imaging data can be even more complex than
the synthetic data we analyze in this controlled setup. This is because real
datasets inherently contain numerous confounding and interacting biases, many
of which may be imperceptible to humans. For instance, it may be necessary to
account for disease prevalence and population shift when evaluating
performance disparities in real data, since the underlying distribution of
other subgroup-associated factors may confound analysis. Glocker and
colleagues (Glocker et al., 2023) account for this by performing test set
resampling to balance the representation of subgroups for a more accurate
representation of algorithmic bias. Conversely, using the SimBA data
generation tool in our framework, these underlying distribution shifts can be
mitigated, as we do here with stratified sampling, or systematically simulated
and evaluated by generating datasets with shifts in the distributions of
subject or disease effects. Thus, this analysis framework provides a powerful
technique for investigating algorithmic encoding of bias and efficacy of bias
mitigation strategies on a range of dataset bias scenarios, in which the
biases could be represented by localized spatial deformations, sampled
distribution of global morphology, or intensity-based simulated artifacts, to
name a few. Therefore, models and methods can be studied and benchmarked on
fully controlled setups prior to making conclusions on real-world data. The
control and flexibility over the dataset composition provided by SimBA, and
the ability to rigorously test in silico trials of the same simulated
“subjects” with different bias manifestations as we present in this framework
facilitates a comprehensive and objective exploration of biases in medical
imaging AI, which was not previously possible.

This first feasibility study was limited by restricting the scope to
simulation of spatially localized morphological bias effects in structural
neuroimaging data and analysis of subsequent impacts on a single CNN model
architecture. However, using this structured analysis framework as a baseline,
we believe that there are numerous avenues that the research community hoping
to gain a better understanding of bias in medical imaging AI can pursue. For
instance, the potential relationship between the magnitude of performance
disparities and proximity of disease and bias regions warrants further
investigation into whether these results are consistent when using vision
transformer-like architectures, where the model learns a more complex
representation of global and local interactions than what is possible in a
standard CNN. Furthermore, this framework could provide particular utility in
the development of new bias mitigation strategies (especially unsupervised
ones), since the user can ensure that all the sources of bias in the data are
being addressed, which is not the case if such a strategy were to be evaluated
on real data. Aside from studying other organs, arbitrary biases, and disease
effects, clinical knowledge of well-defined pathologies (e.g., lesions) could
be integrated to more faithfully model real diseases. Additionally, as users
define the disease effects, predictive tasks are not limited to
classification, but could also be extended to regression or segmentation
scenarios.

In summary, this work introduced a structured framework for analyzing bias
manifestation and mitigation in medical imaging AI. The results show that this
framework has the potential to greatly improve understanding of bias in
medical imaging deep learning pipelines, which alongside studies investigating
the impacts of real-world biased data, will support the development of
clinical AI tools that are more robust, responsible, and fair.

##  5 Data Sharing

Datasets, code, and train/validation/test splits are available for download
from https://github.com/estanley16/SimBA.

##  6 Declaration of Interests

The authors have no competing interests to declare.

## References

  * Adeli et al. (2021) Adeli, E., Zhao, Q., Pfefferbaum, A., Sullivan, E.V., Fei-Fei, L., Niebles, J.C., Pohl, K.M., 2021.  Representation learning with statistical independence to mitigate bias, in: 2021 IEEE Winter Conference on Applications of Computer Vision (WACV), IEEE, Waikoloa, HI, USA. p. 2512–2522.  URL: https://ieeexplore.ieee.org/document/9423445/, doi:10.1109/WACV48630.2021.00256.
  * Arsigny et al. (2006) Arsigny, V., Commowick, O., Pennec, X., Ayache, N., 2006.  A log-euclidean framework for statistics on diffeomorphisms, in: Larsen, R., Nielsen, M., Sporring, J. (Eds.), Medical Image Computing and Computer-Assisted Intervention – MICCAI 2006, Springer, Berlin, Heidelberg. p. 924–931.  doi:10.1007/11866565_113.
  * Banerjee et al. (2023) Banerjee, I., Bhattacharjee, K., Burns, J.L., Trivedi, H., Purkayastha, S., Seyyed-Kalantari, L., Patel, B.N., Shiradkar, R., Gichoya, J., 2023.  “shortcuts” causing bias in radiology artificial intelligence: Causes, evaluation, and mitigation.  Journal of the American College of Radiology 20, 842–851.  doi:10.1016/j.jacr.2023.06.025.
  * Brown et al. (2023) Brown, A., Tomasev, N., Freyberg, J., Liu, Y., Karthikesalingam, A., Schrouff, J., 2023.  Detecting shortcut learning for fair medical ai using shortcut testing.  Nature Communications 14, 4314.  doi:10.1038/s41467-023-39902-7.
  * Burns et al. (2023) Burns, J.L., Zaiman, Z., Vanschaik, J., Luo, G., Peng, L., Price, B., Mathias, G., Mittal, V., Sagane, A., Tignanelli, C., Chakraborty, S., Gichoya, J.W., Purkayastha, S., 2023.  Ability of artificial intelligence to identify self-reported race in chest x-ray using pixel intensity counts.  Journal of Medical Imaging 10, 061106.  doi:10.1117/1.JMI.10.6.061106.
  * Calders et al. (2009) Calders, T., Kamiran, F., Pechenizkiy, M., 2009.  Building classifiers with independency constraints, in: 2009 IEEE International Conference on Data Mining Workshops, IEEE, Miami, FL, USA. p. 13–18.  URL: http://ieeexplore.ieee.org/document/5360534/, doi:10.1109/ICDMW.2009.83.
  * (7) brain development.org, .  Ixi dataset.  URL: https://brain-development.org/ixi-dataset/.
  * Dinsdale et al. (2021) Dinsdale, N.K., Jenkinson, M., Namburete, A.I.L., 2021.  Deep learning-based unlearning of dataset bias for mri harmonisation and confound removal.  NeuroImage 228, 117689.  doi:10.1016/j.neuroimage.2020.117689.
  * Gichoya et al. (2022) Gichoya, J.W., Banerjee, I., Bhimireddy, A.R., Burns, J.L., Celi, L.A., Chen, L.C., Correa, R., Dullerud, N., Ghassemi, M., Huang, S.C., Kuo, P.C., Lungren, M.P., Palmer, L.J., Price, B.J., Purkayastha, S., Pyrros, A.T., Oakden-Rayner, L., Okechukwu, C., Seyyed-Kalantari, L., Trivedi, H., Wang, R., Zaiman, Z., Zhang, H., 2022.  Ai recognition of patient race in medical imaging: a modelling study.  The Lancet Digital Health 4, e406–e414.  doi:10.1016/S2589-7500(22)00063-2.
  * Glocker et al. (2023) Glocker, B., Jones, C., Bernhardt, M., Winzeck, S., 2023.  Algorithmic encoding of protected characteristics in chest x-ray disease detection models.  eBioMedicine 89, 104467.  doi:10.1016/j.ebiom.2023.104467.
  * Ioannou et al. (2022) Ioannou, S., Chockler, H., Hammers, A., King, A.P., 2022.  A study of demographic bias in cnn-based brain mr segmentation, in: Abdulkadir, A., Bathula, D.R., Dvornek, N.C., Habes, M., Kia, S.M., Kumar, V., Wolfers, T. (Eds.), Machine Learning in Clinical Neuroimaging, Springer Nature Switzerland, Cham. p. 13–22.  doi:10.1007/978-3-031-17899-3_2.
  * Jones et al. (2023) Jones, C., Roschewitz, M., Glocker, B., 2023.  The role of subgroup separability in group-fair medical image classification, in: Greenspan, H., Madabhushi, A., Mousavi, P., Salcudean, S., Duncan, J., Syeda-Mahmood, T., Taylor, R. (Eds.), Medical Image Computing and Computer Assisted Intervention – MICCAI 2023, Springer Nature Switzerland, Cham. p. 179–188.  doi:10.1007/978-3-031-43898-1_18.
  * Larrazabal et al. (2020) Larrazabal, A.J., Nieto, N., Peterson, V., Milone, D.H., Ferrante, E., 2020.  Gender imbalance in medical imaging datasets produces biased classifiers for computer-aided diagnosis.  Proceedings of the National Academy of Sciences 117, 12592–12594.  doi:10.1073/pnas.1919012117.
  * Marcinkevics et al. (2022) Marcinkevics, R., Ozkan, E., Vogt, J.E., 2022.  Debiasing deep chest x-ray classifiers using intra- and post-processing methods, in: Proceedings of the 7th Machine Learning for Healthcare Conference, PMLR. p. 504–536.  URL: https://proceedings.mlr.press/v182/marcinkevics22a.html.
  * Meissen et al. (2023) Meissen, F., Breuer, S., Knolle, M., Buyx, A., Müller, R., Kaissis, G., Wiestler, B., Rückert, D., 2023.  (predictable) performance bias in unsupervised anomaly detection URL: http://arxiv.org/abs/2309.14198. arXiv:2309.14198 [cs, eess].
  * Piçarra and Glocker (2023) Piçarra, C., Glocker, B., 2023.  Analysing race and sex bias in brain age prediction, in: Wesarg, S., Puyol Antón, E., Baxter, J.S.H., Erdt, M., Drechsler, K., Oyarzun Laura, C., Freiman, M., Chen, Y., Rekik, I., Eagleson, R., Feragen, A., King, A.P., Cheplygina, V., Ganz-Benjaminsen, M., Ferrante, E., Glocker, B., Moyer, D., Petersen, E. (Eds.), Clinical Image-Based Procedures, Fairness of AI in Medical Imaging, and Ethical and Philosophical Issues in Medical Imaging, Springer Nature Switzerland, Cham. p. 194–204.  doi:10.1007/978-3-031-45249-9_19.
  * Puyol-Antón et al. (2021) Puyol-Antón, E., Ruijsink, B., Piechnik, S.K., Neubauer, S., Petersen, S.E., Razavi, R., King, A.P., 2021.  Fairness in cardiac mr image analysis: An investigation of bias due to data imbalance in deep learning based segmentation, in: de Bruijne, M., Cattin, P.C., Cotin, S., Padoy, N., Speidel, S., Zheng, Y., Essert, C. (Eds.), Medical Image Computing and Computer Assisted Intervention – MICCAI 2021, Springer International Publishing, Cham. p. 413–423.  doi:10.1007/978-3-030-87199-4_39.
  * Rohlfing et al. (2010) Rohlfing, T., Zahr, N.M., Sullivan, E.V., Pfefferbaum, A., 2010.  The sri24 multichannel atlas of normal adult human brain structure.  Human Brain Mapping 31, 798–819.  doi:10.1002/hbm.20906.
  * Seyyed-Kalantari et al. (2021) Seyyed-Kalantari, L., Zhang, H., McDermott, M.B.A., Chen, I.Y., Ghassemi, M., 2021.  Underdiagnosis bias of artificial intelligence algorithms applied to chest radiographs in under-served patient populations.  Nature Medicine 27, 2176–2182.  doi:10.1038/s41591-021-01595-0.
  * Shattuck et al. (2008) Shattuck, D.W., Mirza, M., Adisetiyo, V., Hojatkashani, C., Salamon, G., Narr, K.L., Poldrack, R.A., Bilder, R.M., Toga, A.W., 2008.  Construction of a 3d probabilistic atlas of human cortical structures.  NeuroImage 39, 1064–1080.  doi:10.1016/j.neuroimage.2007.09.031.
  * Smilkov et al. (2017) Smilkov, D., Thorat, N., Kim, B., Viégas, F., Wattenberg, M., 2017.  Smoothgrad: removing noise by adding noise.  arXiv:1706.03825 [cs, stat] URL: http://arxiv.org/abs/1706.03825. arXiv: 1706.03825.
  * Souza et al. (2023) Souza, R., Wilms, M., Camacho, M., Pike, G.B., Camicioli, R., Monchi, O., Forkert, N.D., 2023.  Image-encoded biological and non-biological variables may be used as shortcuts in deep learning models trained on multisite neuroimaging data.  Journal of the American Medical Informatics Association , ocad171doi:10.1093/jamia/ocad171.
  * Stanley et al. (2023) Stanley, E.A.M., Wilms, M., Forkert, N.D., 2023.  A flexible framework for simulating and evaluating biases in deep learning-based medical image analysis, in: Greenspan, H., Madabhushi, A., Mousavi, P., Salcudean, S., Duncan, J., Syeda-Mahmood, T., Taylor, R. (Eds.), Medical Image Computing and Computer Assisted Intervention – MICCAI 2023, Springer Nature Switzerland, Cham. p. 489–499.  doi:10.1007/978-3-031-43895-0_46.
  * Stanley et al. (2022) Stanley, E.A.M., Wilms, M., Mouches, P., Forkert, N.D., 2022.  Fairness-related performance and explainability effects in deep learning models for brain image analysis.  Journal of Medical Imaging 9, 061102.  doi:10.1117/1.JMI.9.6.061102.
  * Weng et al. (2023) Weng, N., Bigdeli, S., Petersen, E., Feragen, A., 2023.  Are sex-based physiological differences the cause of gender bias for chest x-ray diagnosis?, in: Wesarg, S., Puyol Antón, E., Baxter, J.S.H., Erdt, M., Drechsler, K., Oyarzun Laura, C., Freiman, M., Chen, Y., Rekik, I., Eagleson, R., Feragen, A., King, A.P., Cheplygina, V., Ganz-Benjaminsen, M., Ferrante, E., Glocker, B., Moyer, D., Petersen, E. (Eds.), Clinical Image-Based Procedures, Fairness of AI in Medical Imaging, and Ethical and Philosophical Issues in Medical Imaging, Springer Nature Switzerland, Cham. p. 142–152.  doi:10.1007/978-3-031-45249-9_14.
  * Wu et al. (2022) Wu, Y., Zeng, D., Xu, X., Shi, Y., Hu, J., 2022.  Fairprune: Achieving fairness through pruning for dermatological disease diagnosis, in: Wang, L., Dou, Q., Fletcher, P.T., Speidel, S., Li, S. (Eds.), Medical Image Computing and Computer Assisted Intervention – MICCAI 2022, Springer Nature Switzerland, Cham. p. 743–753.  doi:10.1007/978-3-031-16431-6_70.
  * Yearley et al. (2023) Yearley, A.G., Goedmakers, C.M.W., Panahi, A., Doucette, J., Rana, A., Ranganathan, K., Smith, T.R., 2023.  Fda-approved machine learning algorithms in neuroradiology: A systematic review of the current evidence for approval.  Artificial Intelligence in Medicine 143, 102607.  doi:10.1016/j.artmed.2023.102607.
  * Zhang et al. (2022) Zhang, H., Dullerud, N., Roth, K., Oakden-Rayner, L., Pfohl, S., Ghassemi, M., 2022.  Improving the fairness of chest x-ray classifiers, in: Proceedings of the Conference on Health, Inference, and Learning, PMLR. p. 204–233.  URL: https://proceedings.mlr.press/v174/zhang22a.html.

##  Appendix A Implementation Details

Model and Training. The convolutional neural network (CNN) we use for all
experiments contains five convolutional blocks (filters=32, 64, 128, 256, 512)
comprised of 3x3x3 convolution filters, batch normalization, sigmoid
activation, and 2x2x2 max pooling, followed by average pooling, 20% dropout, a
flatten layer, and a dense classification layer. A batch size of four and
learning rate of 1e-4 was used, and models were trained until convergence, by
validation loss early stopping patience=15. For each bias scenario (No Bias,
Near Bias, Far Bias), we trained the same model architecture using the same
five weight initialization seeds with deterministic GPU state, in which each
seed replicate uses the same dataset train, validation, and test splits of
50%/25%/25%. Each split was stratified by subject effects, disease effects,
bias group label, and disease class label. Thus, evaluation of each bias
scenario is a counterfactual of the same model and same dataset, where the
difference is in the absence of the bias effect, or the bias effect in
different locations. The models were implemented in Keras/Tensorflow v. 2.10
using a NVIDIA GeForce RTX 3090 GPU.

Reweighing. The weight W for samples in target class D and bias group B is
given by:

$W(B=b|D=d)=\dfrac{P_{exp}(b\land d)}{P_{exp}(b\land d)}$

Where:

$P_{exp}(b\land d)=\text{expected probability of an unbiased dataset.}$

And:

$P_{act}(b\land d)=\text{actual probability in dataset.}$

Thus, the weights for each class and bias group are:

$W(B=1|D=1)=\dfrac{d(\gamma+\alpha)}{\gamma}$
$W(B=1|D=0)=\dfrac{(1-d)(\gamma+\alpha)}{\alpha}$
$W(B=1=0|D=1)=\dfrac{d(2-\gamma-\alpha)}{1-\gamma}$
$W(B=0|D=0)=\dfrac{(1-d)(2-\gamma-\alpha)}{1-\alpha}$

Where:

$\alpha=\text{fraction of bias group in non-disease class.}$
$\gamma=\text{fraction of bias group in disease class.}$ $d=\text{fraction of
disease class in full dataset.}$

Bias Unlearning. The goal of the unlearning strategy is to modify the weights
of the feature encoder (i.e., the entire CNN except for the dense
classification layer) such that the model is no longer able to predict the
bias group while retaining the ability to predict the disease. The same
feature encoder backbone was used with two classification heads (i.e., the
final dense layer), one for the disease classification and one for the bias
group classification. First, the encoder and disease prediction head were
trained until convergence. Then, the encoder was frozen and the disease
prediction head was replaced with a bias prediction head, which was trained to
convergence. Subsequently, the unlearning process takes place, in which the
cross-entropy loss for predicting the disease is minimized, while a confusion
loss for predicting the bias group is maximized. This iterative training takes
place until the disease prediction accuracy has stopped improving and the bias
prediction accuracy has stopped decreasing. In our experiments, only five
epochs were necessary to unlearn the bias group. Bias prediction results are
reported in Supplementary Table 6.

Group Models. First, the model was pre-trained for 5 epochs on the full
dataset. Then, using the pre-trained model as a basis, separate models were
trained on both the bias and non-bias groups. Accuracy, true positive rates,
and false positive rates were computed for each group model individually, and
then compared to get the reported subgroup performance disparities.

##  Appendix B Supplementary Tables

Table 2: SimBA effect sampling parameters.

Subject Disease Class Non-Disease Class Bias Group Non-Bias Group Number of
Samples 2002 1000 1002 1003 999 Region Whole brain Left insular cortex Left
insular cortex Near: Left putamen Far: Right postcentral gyrus - Sampling
Distribution N(0,1)$\in$[-3.5,3.5] N(-1,1) N(1,1) N(2,1) - Number of Principal
Components Sampled From 10 1 1 1 - Number of Stratification Bins 10 10 10 - -

Table 3: Composition of dataset splits.

Train Validation Test Disease Non-Disease Disease Non-Disease Disease Non-
Disease Non-Bias 140 360 72 176 68 183 Bias 360 145 180 70 180 68

Table 4: True positive rate statistical tests. Bold font indicates groups in
comparison test. Two-tailed paired t-test unless indicated by *, then Wilcoxon
matched-pairs signed rank test. †significant after Bonferroni correction
($\alpha$=0.005).

Model Bias Scenario P-value Naïve No Bias – Near Bias <0.0001† Naïve No Bias –
Far Bias 0.0002† Naïve Near Bias – Far Bias 0.1037 Naïve – Reweighing Near
Bias 0.0001† Naïve – Unlearning Near Bias 0.0075 Naïve – Group Models Near
Bias 0.8534 Naïve – Reweighing* Far Bias 0.0625 Naïve – Unlearning Far Bias
0.0104 Naïve – Group Models* Far Bias 0.0625 Naïve – Group Models No Bias
0.0007†

Table 5: False positive rate statistical tests. Bold font indicates groups in
comparison test. Two-tailed paired t-test unless indicated by *, then Wilcoxon
matched-pairs signed rank test. †significant after Bonferroni correction
($\alpha$=0.005).

Model Bias Scenario P-value Naïve No Bias – Near Bias <0.0001† Naïve No Bias –
Far Bias 0.0008† Naïve Near Bias – Far Bias 0.0552 Naïve – Reweighing Near
Bias 0.0002† Naïve – Unlearning Near Bias 0.0525 Naïve – Group Models Near
Bias 0.6071 Naïve – Reweighing Far Bias 0.0058 Naïve – Unlearning* Far Bias
0.0625 Naïve – Group Models Far Bias 0.0289 Naïve – Group Models No Bias
0.0001†

Table 6: Results of bias group prediction after unlearning.

Bias Scenario Bias Prediction Accuracy (%) No Bias 49.94 ± 0.33 Near Bias
47.17 ± 4.23 Far Bias 46.73 ± 7.98

Table 7: Weighted saliency scores (%) for the Near Bias region across five
model seeds.

Disease Class Non-Disease Class Bias Scenario Bias Mitigation Strategy Non-
Bias Group Bias Group Non-Bias Group Bias Group No Bias Naïve 47.32 ± 7.23
47.18 ± 8.35 30.34 ± 7.55 30.55 ± 6.85 Reweighing 40.97 ± 4.37 39.88 ± 4.73
23.25 ± 5.97 23.48 ± 5.62 Unlearning 46.14 ± 7.66 44.98 ± 6.86 30.67 ± 6.84
29.52 ± 6.51 Group Models 46.91 ± 1.40 45.92 ± 6.64 30.76 ± 7.95 29.53 ± 8.58
Near Bias Naïve 77.09 ± 12.11 75.68 ± 14.42 45.78 ± 7.01 55.81 ± 16.84
Reweighing 41.76 ± 4.45 42.7 ± 3.89 19.22 ± 5.68 20.51 ± 5.52 Unlearning 59.64
± 13.28 54.97 ± 9.37 36.02 ± 12.65 37.34 ± 13.98 Group Models 42.15 ± 4.03
39.28 ± 6.95 27.1 ± 4.92 21.59 ± 4.70

Table 8: Weighted saliency scores (%) for the Far Bias region across five
model seeds.

Disease Class Non-Disease Class Bias Scenario Bias Mitigation Strategy Non-
Bias Group Bias Group Non-Bias Group Bias Group No Bias Naïve 1.69 ± 0.76 1.67
± 0.72 1.22 ± 0.42 1.14 ± 0.32 Reweighing 1.34 ± 0.70 1.31 ± 0.69 1.35 ± 0.67
1.34 ± 0.67 Unlearning 1.13 ± 0.52 1.12 ± 0.53 0.98 ± 0.25 0.95 ± 0.28 Group
Models 1.91 ± 0.99 1.5 ± 0.55 1.59 ± 1.00 1.41 ± 0.63 Near Bias Naïve 7.77 ±
3.35 7.85 ± 3.32 7.47 ± 3.56 7.39 ± 3.55 Reweighing 2.23 ± 0.79 2.09 ± 0.66
1.48 ± 0.70 1.48 ± 0.70 Unlearning 10.47 ± 6.04 10.41 ± 6.16 8.66 ± 4.25 8.65
± 4.46 Group Models 1.91 ± 0.60 3.12 ± 1.53 1.6 ± 0.49 2.64 ± 1.70

##  Appendix C Supplementary Figures

Figure 4: Kernel density estimate/scatter plots of subject and disease effect
magnitude sampled from principal component analysis models. The distribution
of subject effect magnitude within each target class, or the joint
distribution of subject and disease effect magnitudes within each bias group
can be identified by deep learning models and lead to performance disparities
if there are differences between subgroups. Therefore, it is important to
stratify such effects as shown in the figure to minimize this distribution-
induced bias. The sampling distributions shown in the figure are identical for
each bias scenario. Left: All samples, middle: non-disease class, right:
disease class.  Figure 5: True (left) and false (right) positive rate
disparities for the naïve models and group models. ***p$<$0.001.