Double $J/\psi$ production as a test of
parton correlations in double parton scattering
Sergey Koshkarev and Stefan Groote
Institute of Physics, University of Tartu, 51010 Tartu, Estonia

###### Abstract

Using the GS09 model we predict the possible impact of the parton correlation
on the $J/\psi$-pair production at the Spin Physics Detector at the Nuclotron-
based Ion Collider Facility. The double $J/\psi$ production and the effective
cross sections are calculated.

##  1 Motivation

With the advent of high luminosity accelerators, multi parton scattering (MPS)
physics at high-energy hadron colliders like the Tevatron at Fermilab and the
Large Hadron Collider (LHC) at CERN, in particular those observing double
parton scattering (DPS), has become one of the hottest topics of the modern
experimental particle physics [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14,
15, 16, 18, 19, 20, 21]. DPS results are usually interpreted in terms of the
so-called “pocket formula” [22]

$\sigma^{AB}_{\rm DPS}=\frac{m}{2}\frac{\sigma^{A}_{\rm SPS}\sigma^{B}_{\rm
SPS}}{{\sigma_{\rm eff}}},$ (1)

where the combinatorial factor is $m=1$ for identical final states $A=B$ and
$m=2$ for $A\neq B$. This formula is derived under the assumption of
independent parton scatterings, with

${\sigma_{\rm eff}}=\bigg{[}\int d^{2}b~{}T^{2}({\bf b})\bigg{]}^{-1},$

where $T({\bf b})$ is the overlap function that characterizes the transverse
area occupied by the interacting partons in the impact parameter space ${\bf
b}$. Such a definition leads to the assumption that the value of ${\sigma_{\rm
eff}}$ should be universal, independent of the final states or the phase
space. On the other hand, the DØ collaboration provided two different
measurement of DPS for $\gamma+3$ jet events, namely $\langle{\sigma_{\rm
eff}}\rangle=16\pm 0.3\pm 2.3$ in 2010 [5] and ${\sigma_{\rm eff}}^{\rm
incl}=12.7\pm 0.2\pm 1.3$ in 2014 [9]. It is easy to see that not only the
central values are different, but the errors bars do not overlap neither. In
contrast to many other results, the only difference between these two
measurements is the kinematic region (cf. Tab. 1). This definitely contradicts
the “universality” assumption for ${\sigma_{\rm eff}}$.

In Refs. [23, 24, 25] the importance of the effect of the evolution of the
double parton distribution function (dPDF) is emphasized. Most of the
phenomenological analysis of a possible impact of the dPDF effect are focused
on the LHC energies (cf. Refs. [26, 27] and references therein). Indeed,
standard wisdom tells us that DPS should be more preferable at high energies.
On the other hand, possible dPDF effect is expected at higher values of the
Bjorken-$x$ or with a significant gap between the $x$ values. However, even
for double W-boson production the typical value for the Bjorken-$x$ is
$\langle x\rangle\sim 0.01$. In contrast to that, preferable conditions can be
easily achieved at lower energy experiments, where the DPS contribution could
be relatively small but still expected to be far from zero [36, 37].

The Spin Physics Detector at the Nuclotron-based Ion Collider fAcility (NICA)
collider (JINR, Dubna) is a universal facility to investigate the spin
structure of the proton and deuteron and the other spin-related phenomena with
polarized proton and deuteron beams at a collision energy up to
$27\mbox{\,GeV}$ and a luminosity up to $10^{32}$ ${\rm cm^{-2}~{}s^{-1}}$
[29]. It is easy to see that the production threshold of the pair of $J/\psi$
is already more than $20\%$ of the NICA energies, leading to the typical value
for the Bjorken-$x$ of $\langle x\rangle>0.1$. The relatively small cross
section can be compensated by the high luminosity.

In this paper we phenomenologically investigate the possible impact of the
evolution of the dPDF on the production cross section of double $J/\psi$
events and the measurements of ${\sigma_{\rm eff}}$ at NICA energies.

Table 1:  The effective cross section measured in $\gamma+3$ jets process in $p\bar{p}$ collisions at $\sqrt{s}=1.96\mbox{\,TeV}$ measured by the DØ experiment at the Tevatron. Year | Final State | $p_{T}^{\rm min}(\mbox{\,GeV}/c)$ |  $\eta(y)$ range |  ${\sigma_{\rm eff}}$ (mb)
---|---|---|---|---
| 2010
---
[5]
$\gamma+3$ jets |  | $60<p_{T}^{\gamma}<80$
---
$p_{T}^{\rm jet\,1}>25$
$p_{T}^{\rm jet\,2,3}>15$
| $|y^{\gamma}|<1.0$
---
$1.5<|y^{\gamma}|<2.5$
$|y^{\rm jet}|<3.0$
$\langle{\sigma_{\rm eff}}\rangle=16\pm 0.3\pm 2.3$
| 2014
---
[9]
|  $\gamma+3$ jets
---
$\gamma+c/b+2$ jets
| $p_{T}^{\gamma}>26$
---
$p_{T}^{\rm jet\,1}>15$
$15<p_{T}^{\rm jet\,2,3}<35$
| $|\eta^{\gamma}|<1.0$
---
$1.5<|\eta^{\gamma}|<2.5$
$|\eta^{\rm jet}|<2.5$
| ${\sigma_{\rm eff}}^{\rm incl}=12.7\pm 0.2\pm 1.3$
---
${\sigma_{\rm eff}}^{\rm HF}=14.6\pm 0.6\pm 3.2$

##  2 Phenomenological basement

If the final states $A$ and $B$ are produced in a DPS process independently,
Eq. (1) can be phenomenologically cast into the form

$\displaystyle\sigma^{AB}_{\rm DPS}$ $\displaystyle=$
$\displaystyle\frac{m}{2}\frac{1}{{\sigma_{\rm eff}}}\int dx_{1}\ldots dx_{4}\
~{}\ f(x_{1},Q_{A})f(x_{2},Q_{A})\hat{\sigma}_{A}(x_{1},x_{2})\times$
$\displaystyle
f(x_{3},Q_{B})f(x_{4},Q_{B})\hat{\sigma}_{B}(x_{3},x_{4})\theta(1-x_{1}-x_{3})\theta(1-x_{2}-x_{4}),$

where $f(x,Q)$ denotes the parton distribution function, $\hat{\sigma}$ is the
cross section at parton level, and $\theta$ is the Heaviside step
function.111We use odd indexes for the “first” proton and even indexes for the
“second” proton.

As we already mentioned, some models predict the correlation between partons
of a same hadron [23, 25]. In this case, Eq. (1) should be written in the form

$\sigma^{AB}_{\rm DPS}=\frac{m}{2}\frac{1}{{\sigma_{\rm eff}}}\int
dx_{1}\ldots
dx_{4}~{}D(x_{1},x_{3},Q_{A},Q_{B})D(x_{2},x_{4},Q_{A},Q_{B})\hat{\sigma}_{A}(x_{1},x_{2})\hat{\sigma}_{B}(x_{3},x_{4}),$

where e.g. the dPDF $D(x_{1},x_{3},Q_{A},Q_{B})\neq
f(x_{1},Q_{A})f(x_{3},Q_{B})\theta(1-x_{1}-x_{3})$ is not the product of the
two sPDFs. For the possible impact of the dPDFs we can write the ratio

$R_{\Delta}(x_{1},x_{2},x_{3},x_{4},Q_{A},Q_{B})=\frac{D(x_{1},x_{3},Q_{A},Q_{B})D(x_{2},x_{4},Q_{A},Q_{B})}{f(x_{1},Q_{A})f(x_{2},Q_{A})f(x_{3},Q_{B})f(x_{4},Q_{B})}.$
(2)

In single parton scattering, the evolution of the sPDF is described by the
DGLAP QCD evolution equations, formerly known as the Altarelli–Parisi
equations. As shown in Ref. [28], these renormalisation group equations can be
solved by imposing two sum rules, namely the momentum and the number sum rule.
While for single parton scattering the sums over momenta fractions and the
number of valence partons stay constant, the sum rules for double parton
scattering are more involved and, as it turns out, more restrictive.

Corresponding to the (s)DGLAP for single parton scattering, Gaunt and Stirling
postulate a dDGLAP which is based on a couple of principles explained in
detail in Ref. [23]. As for the sDGLAP, an increasing scale $t=\ln(Q^{2})$
allows for the splitting of partons (by, e.g., the emission of a gluon),
described by the splitting functions $P_{i\to j}(x)$ which denote a parton $i$
splitting to $j$ with a momentum fraction $x$. For double parton scattering,
the splitting functions $P_{i\to jk}(x)$ can be understood as the parton $i$
splitting to a parton $j$ with momentum fraction $x$, and a parton $k$ with
momentum fraction $1-x$. This latter is certainly true in leading order (LO)
scattering but has to be modified in next-to-leading order (NLO). For this
reason, in Ref. [23] Gaunt and Stirling deal with the LO scattering only,
postponing the NLO analysis to a future publication. In LO, the sum rules to
be postulated have to take into account two very basic correlations, namely

  1. 1.

having found a quark with a given flavour, the probability is smaller to find
a second quark with the same flavour, and

  2. 2.

having found a parton with momentum fraction $x$, the probability is smaller
to find a second parton with momentum faction $1-x$.

Based on three splitting diagrams for increasing and two diagrams for
decreasing population of partons, the dDGLAP equation is justified in detail
(cf. Eq. (2.1) and Fig. 3 in Ref. [23]). The new sum rules are postulated, and
a procedure is developed to solve the dDGLAP on a grid of momentum fractions
and scale, known as the implementation of the GS09 model.

##  3 $J/\psi$-pair production from DPS with the “pocket formula”

As a first step, it is interesting to estimate the DPS effect with the “pocket
formula” approximation (1). In order to provide such an estimate, we are going
to use experimentally measured $J/\psi$ production cross sections at similar
energies. Using the CERN proton beam at $400\mbox{\,GeV}$/c to produce charm
particles with incident on different nuclear targets, the NA3 experiment
provided data on the production of $J/\psi$ pairs on a platinum target with
the production cross sections of $27\pm 10\mbox{\,pb}$ per nucleon [30], and
the NA50 experiment measured single $J/\psi$ production
$Br(J/\psi\to\mu^{+}\mu^{-})\times\sigma(J/\psi)$ on Be, Al, Cu, Ag, W, and Pb
targets in the range between $4.717\pm 0.026\mbox{\,nb}$ and $3.715\pm
0.016\mbox{\,nb}$ per nucleon [31].

Combining the mean values of the single $J/\psi$ cross section for proton-W
and proton-Pb collisions, $\sigma(J/\psi)\approx 12.5\mbox{\,$\mu$b}$, and the
effective cross section $\langle{\sigma_{\rm eff}}\rangle=4.6\mbox{\,mb}$ for
double $J/\psi$ production (cf. Tab. 2), we can calculate the contribution of
DPS. However, $\sigma(J/\psi)$ as measured by the NA50 collaboration includes
$J/\psi$ from the feed-down effect. To consider this effect, we cast formula
(1) into the form

$\sigma_{DPS}(J/\psi J/\psi)=\frac{\sigma(J/\psi)^{2}}{{\sigma_{\rm
eff}}}\bigg{(}\frac{r^{2}_{J/\psi}+r^{2}_{\psi^{\prime}}+r^{2}_{\chi_{c}}}{2}+2\cdot(r_{J/\psi}r_{\psi^{\prime}}+r_{J/\psi}r_{\chi_{c}}+r_{\psi^{\prime}}r_{\chi_{c}})\bigg{)}$

(cf. also the discussion in Ref. [32]), where $r_{J/\psi}=0.62\pm 0.04$,
$r_{\psi^{\prime}}=0.08\pm 0.02$, and $r_{\chi_{c}}=0.30\pm 0.08$ are feed-
down fractions [33]. Performing this calculation, we obtain
$\sigma_{DPS}(J/\psi J/\psi)=2.6\mbox{\,pb}$. Utilizing the double $J/\psi$
production cross section measured by the NA3 experiment, the fraction of DPS
events can be estimate to be $f_{DP}\approx 9.6\%$.

Table 2:  The effective cross section values measured in double and triple $J/\psi$ production at DØ, ATLAS, and CMS experiments. Experiment |  $\sqrt{s}$,  TeV | Colliding Mode |  ${\sigma_{\rm eff}}$,  mb
---|---|---|---
DØ [10] | 1.96 | $p\bar{p}$ | $4.8\pm 0.5\text{(stat)}\pm 2.5\text{(syst)}$
ATLAS [17] | 8 | $pp$ | $6.3\pm 1.6\text{(stat)}\pm 1.0\text{(syst)}$
CMS [21] | 13 | $pp$ | $2.7^{+1.4}_{-1.0}\text{(exp)}^{+1.5}_{-1.0}\text{(th)}$

Figure 1:  Differential cross section distributions for double $J/\psi$
production, as predicted by single PDF (squared points) and double PDF
(circular points).

##  4 Possible effects of double PDF within the GS09 model

In order to distinguish between single and double PDF predictions, we use a
Pythia 8 Monte-Carlo simulation [34], where the sPDFs are calculated with
MSTW2008LO [35] and the dPDFs are calculated in the GS09 model employing the
ratio (2), where $R_{\Delta}$ was calculated for every single event.

The effect of double PDF for double $J/\psi$ production cross section at NICA,
differential in the normalized angular difference for the polar angle,
$\Delta\theta_{\pi}=(\theta(J/\psi_{1})-\theta(J/\psi_{2}))/\pi$, is presented
in Fig. 1. In case of the single PDF prediction, the first bin
$\Delta\theta_{\pi}<0.25$ contains more than $40\%$ of the statistics. In
contrast to that, double PDF predicts almost equal statistics for all bins,
except for the bin $\Delta\theta_{\pi}>0.75$, where almost no events are
found.

It is also interesting to calculate ${\sigma_{\rm eff}}$ in the style of a
“pocket formula” analysis. Our calculation gives $\langle
R_{\Delta}\rangle\approx 0.37$ that leads us to the effect of the dPDFs,
namely the amplification of ${\sigma_{\rm eff}}=\langle{\sigma_{\rm
eff}}\rangle/\langle R_{\Delta}\rangle\approx 12.4\mbox{\,mb}$. This value is
much higher than the value previously measured by DØ, ATLAS and CMS at low
Bjorken-$x$. Obviously, our result also reduces the DPS fraction to
$f_{DP}\approx 3.6\%$.

##  5 Discussion and Summary

Current theoretical and experimental knowledge about parton distribution
functions do not provide clear evidence for the parton correlation. To resolve
this issue, we propose to measure the double $J/\psi$ production cross section
at NICA energies, differential in the normalized angular difference for the
polar angle, $d\sigma_{DPS}(J/\psi J/\psi)/d\Delta\theta_{\pi}$. As it turns
out, the significant difference in distributions predicted in sPDF and dPDF
does not require a precision fit of data. therefore, in the analysis we can
rely on counting of events in different $\Delta\theta_{\pi}$ bins only.

In order to distinguish DPS from SPS, we proposed to use the (normalized)
angular difference for the azimuthal angle of radiation of the $J/\psi$ pair,
given by $\Delta\phi_{\pi}=(\phi(J/\psi_{1})-\phi(J/\psi_{2}))/\pi$. Having
taken into account the fact that $\Delta\phi_{\pi}$ has a peak near 1 for SPS
but a flat shape for DPS, we were able to exclude the region
$\Delta\phi_{\pi}\sim 1$ in order to maximize the DPS/SPS ratio.

### Acknowledgements

The authors would like to acknowledge J.R. Gaunt for his help in understanding
how to run the GS09 model code and providing code examples. We also would like
to thank D. Bandurin and G. Golovanov for very useful discussions about DPS,
and A. Gridin and A. Guskov for their comments about the Spin Physics Detector
and NICA. The research was supported by the European Regional Development Fund
under Grant No. TK133.

## References

  * [1] T. Åkesson et al. [Axial Field Spectrometer], Z. Phys. C 34 (1987), 163
  * [2] J. Alitti et al. [UA2], Phys. Lett. B 268 (1991), 145-154
  * [3] F. Abe et al. [CDF], Phys. Rev. D 47 (1993), 4857-4871
  * [4] F. Abe et al. [CDF], Phys. Rev. D 56 (1997), 3811-3832
  * [5] V. M. Abazov et al. [D0], Phys. Rev. D 81 (2010), 052012
  * [6] R. Aaij et al. [LHCb], JHEP 06 (2012), 141
  * [7] G. Aad et al. [ATLAS], New J. Phys. 15 (2013), 033038
  * [8] S. Chatrchyan et al. [CMS], JHEP 03 (2014), 032
  * [9] V. M. Abazov et al. [D0], Phys. Rev. D 89 (2014) no.7, 072006
  * [10] V. M. Abazov et al. [D0], Phys. Rev. D 90 (2014) no.11, 111101
  * [11] G. Aad et al. [ATLAS], Eur. Phys. J. C 75 (2015) no.5, 229
  * [12] R. Aaij et al. [LHCb], JHEP 07 (2016), 052
  * [13] V. M. Abazov et al. [D0], Phys. Rev. D 93 (2016) no.5, 052008
  * [14] V. M. Abazov et al. [D0], Phys. Rev. Lett. 116 (2016) no.8, 082002
  * [15] M. Aaboud et al. [ATLAS], JHEP 11 (2016), 110
  * [16] V. Khachatryan et al. [CMS], JHEP 05 (2017), 013
  * [17] M. Aaboud et al. [ATLAS], Eur. Phys. J. C 77 (2017) no.2, 76
  * [18] R. Aaij et al. [LHCb], JHEP 06 (2017), 047 [erratum: JHEP 10 (2017), 068]
  * [19] A. M. Sirunyan et al. [CMS], JHEP 02 (2018), 032
  * [20] M. Aaboud et al. [ATLAS], Phys. Lett. 790 (2019), 595
  * [21] S. Leontsinis [CMS], Rev. Mex. Fis. Suppl. 3 (2022) no.3, 0308046
  * [22] D. d’Enterria and A. Snigirev, Adv. Ser. Direct. High Energy Phys. 29 (2018), 159-187
  * [23] J. R. Gaunt and W. J. Stirling, JHEP 03 (2010), 005
  * [24] A. M. Snigirev, Phys. Rev. D 81 (2010), 065014
  * [25] M. Rinaldi, S. Scopetta, M. Traini and V. Vento, JHEP 12 (2014), 028
  * [26] F. A. Ceccopieri, M. Rinaldi and S. Scopetta, Phys. Rev. D 95 (2017) no.11, 114030
  * [27] J. R. Gaunt, C. H. Kom, A. Kulesza and W. J. Stirling, Eur. Phys. J. C 69 (2010), 53-65
  * [28] R. K. Ellis, W. J. Stirling and B. R. Webber, Camb. Monogr. Part. Phys. Nucl. Phys. Cosmol. 8 (1996), 1-435
  * [29] V. M. Abazov et al. [SPD proto], [arXiv:2102.00442 [hep-ex]].
  * [30] J. Badier et al. [NA3], Phys. Lett. B 158 (1985), 85
  * [31] B. Alessandro et al. [NA50], Eur. Phys. J. C 48 (2006), 329
  * [32] J. P. Lansberg and H. S. Shao, Phys. Lett. B 751 (2015), 479-486
  * [33] S. Digal, P. Petreczky and H. Satz, Phys. Rev. D 64 (2001), 094015
  * [34] T. Sjöstrand, S. Mrenna and P. Z. Skands, Comput. Phys. Commun. 178 (2008), 852-867
  * [35] A. D. Martin, W. J. Stirling, R. S. Thorne and G. Watt, Eur. Phys. J. C 63 (2009), 189-285
  * [36] J. P. Lansberg and H. S. Shao, Nucl. Phys. B 900 (2015), 273-294
  * [37] S. Koshkarev, [arXiv:1909.06195 [hep-ph]].


longtable captionstableheading

# Everything, Everywhere All in One Evaluation: Using Multiverse Analysis to
Evaluate the Influence of Model Design Decisions on Algorithmic Fairness

Jan Simson  Florian Pfisterer  Christoph Kern

###### Abstract

A vast number of systems across the world use algorithmic decision making
(ADM) to (partially) automate decisions that have previously been made by
humans. When designed well, these systems promise more objective decisions
while saving large amounts of resources and freeing up human time. However,
when ADM systems are not designed well, they can lead to unfair decisions
which discriminate against societal groups. The downstream effects of ADMs
critically depend on the decisions made during the systems’ design and
implementation, as biases in data can be mitigated or reinforced along the
modeling pipeline. Many of these design decisions are made implicitly, without
knowing exactly how they will influence the final system. It is therefore
important to make explicit the decisions made during the design of ADM systems
and understand how these decisions affect the fairness of the resulting
system.

To study this issue, we draw on insights from the field of psychology and
introduce the method of multiverse analysis for algorithmic fairness. In our
proposed method, we turn implicit design decisions into explicit ones and
demonstrate their fairness implications. By combining decisions, we create a
grid of all possible “universes” of decision combinations. For each of these
universes, we compute metrics of fairness and performance. Using the resulting
dataset, one can see how and which decisions impact fairness. We demonstrate
how multiverse analyses can be used to better understand variability and
robustness of algorithmic fairness using an exemplary case study of predicting
public health coverage of vulnerable populations for potential interventions.
Our results illustrate how decisions during the design of a machine learning
system can have surprising effects on its fairness and how to detect these
effects using multiverse analysis.

##  1 Introduction

Across the world, more and more decisions are being made with the support of
machine learning (ML) and algorithms; so called algorithmic decision making
(ADM). Examples of such systems can be found in finance for loan approvals
(Mukerjee et al. 2002), the labor market for hiring decisions or filtering
resumes (Faliagka, Ramantas, and Tzimas 2012) and the criminal justice system
to assess risks of recidivism (Angwin et al. 2016). While these systems are
very promising when designed well, raising hopes of more accurate and
objective decisions, their impact can be quite the opposite when designed
wrongly. There are many examples of ADM systems discriminating against people
(Mehrabi et al. 2021). One prominent example was the _robodebt_ system, where
the Australian government used an algorithm to detect potential social
security overpayments. Due to serious flaws in the design of the system it
often overestimated debts and put the burden on the accused to prove the
contrary (Henriques-Gomes 2023). Other examples include the Dutch childcare
benefits system using an ADM system that was much more likely to accuse
immigrants of having committed fraud (International 2021).

These fairness problems often occur because algorithms replicate biases in the
underlying training data. However, biases can be amplified throughout the
machine learning pipeline depending on how exactly data is processed and
turned into outputs (Kern et al. 2021; Rodolfa, Saleiro, and Ghani 2020).
Unfortunately, no single straightforward method exists to prevent biases in
the machine learning pipeline (Agrawal et al. 2021). Understanding how
modeling decisions interact with fairness is therefore a prerequisite for
effectively mitigating unintended outcomes in practice. A systematic mapping
of design decisions to fairness outcomes can critically guide the model
selection process as multiple models may achieve similar accuracy, but can
considerably differ in their fairness properties (Black, Raghavan, and Barocas
2022). As a result, preventing algorithms from introducing new or reinforcing
existing biases requires careful study and evaluation of the – often implicit
– decisions made while designing a machine learning system. To facilitate this
objective in a systematic and efficient way, we introduce the method of
multiverse analysis for algorithmic fairness. Multiverse analyses were
introduced to psychology with the intent to improve reproducibility and create
more robust research (Steegen et al. 2016). We adapt this methodology across
domains to work in the context of machine learning with a focus on evaluating
metrics of algorithmic fairness.

In the following, we present a generalizable approach of using multiverse
analysis to estimate the effect of decisions during the design of a machine
learning or ADM system on fairness outcomes. We demonstrate the feasibility of
this approach using a case study of predicting public health coverage in US
census data. We provide modular source code to allow streamlined adaptation of
the proposed method in other use cases and contexts.

###  1.1 Multiverse Analysis

Multiverse analyses were first introduced in psychology by Steegen et al.
(2016) in response to the reproducibility crisis affecting the field (OPEN
SCIENCE COLLABORATION 2015). The goal of this analysis type is to investigate
the invariance of results to researchers’ analysis decisions. Specifically,
when analyzing a dataset, researchers make many implicit and explicit choices
(Simmons, Nelson, and Simonsohn 2011), often without the option of confirming
whether a choice is correct or incorrect. This leads to many plausible
scenarios when analyzing data, as one traverses a _garden of forking paths_
(Gelman and Loken 2014), where each fork corresponds to a decision. The
multitude of these scenarios becomes especially evident when multiple
researchers analyze the same data, coming to staggeringly different results
(Breznau et al. 2022).

Multiverse analysis focuses on the pre-processing steps applied to a dataset:
Steps such as selecting the observations and predictor variables to include in
a dataset or scaling and binning their values. Based on the different
decisions made and paths taken when pre-processing a dataset, analysts will
end up with one of many possible datasets for the actual analysis. In a
multiverse analysis, the goal is to make this variation explicit by using the
complete grid of decisions and their options to generate all plausible
datasets. Using all potential datasets, a multiverse analysis re-runs the
analysis on each of them to receive the distribution of results instead of a
single result point (Figure 1, Steps 1 - 3). We adapt this methodology for the
machine learning context with a special focus on using it to generate insights
on metrics of algorithmic fairness.

Besides multiverse analyses, a highly related type of analysis emerged around
the same time in the specification curve analysis (Simonsohn 2018). Similarly
to a multiverse analysis, a specification curve analysis uses the complete
grid of possible decision combinations to estimate the variability of research
findings. This data is used to create a specification curve, a graph
displaying the distribution of the effect size or coefficient of interest
using a single curve. All decision combinations are ordered by the metric
creating a single, monotonically increasing curve. The actual decisions are
one-hot encoded and displayed as a binary rug below the curve, clearly
identifying all options that produced a certain value on the curve. Both
multiverse analysis and specification curve analysis focus on bringing the
diversity in pre-processing (and beyond) to light and we incorporated the
essential ideas from both approaches.

Figure 1: Steps to conduct a multiverse analysis for algorithmic fairness.
Steps 1 - 3 apply to multiverse analyses in general, whereas steps 4 - 5 are
unique to multiverse analyses for algorithmic fairness.

###  1.2 Multiverse Analysis for Algorithmic Fairness

In our proposed adaptation of multiverse analysis for algorithmic fairness,
one starts by compiling a list of all potentially relevant decisions that are
being made during the design of a particular system. We differentiate between
different kinds of decisions in this context: (1) decisions which are already
made explicitly with a consideration of their different options e.g. choice of
model and its hyperparameters, and (2) decisions which are made explicitly,
but without any consideration for alternatives e.g. log-transforming an income
column because it is common practice. In a multiverse analysis, the goal is to
turn both types of decisions into completely explicitly made decisions and
evaluate their impacts. There are also decisions which may initially not even
be considered as such e.g. modifying classification cutoffs post-hoc due to
external constraints. Conducting a multiverse analysis invites reflection on
the modeling pipeline such that implicit decisions may surface and are turned
into explicit ones. One of the key differences in the present analysis
compared to a classic multiverse analysis is that we will evaluate machine
learning systems, whereas classical multiverse analyses will typically
evaluate the outcomes of null-hypothesis-significance-tests (NHST) across
analysis choices. While many of the decision points apply to any machine
learning system (e.g., choice of algorithm, how to preprocess certain
variables, cross-validation splits), many of them are also domain-specific
(e.g., coding of certain variables, how to set classification thresholds, how
fairness is operationalized). We focus on decisions made during the pre-
processing of data, in line with the original approach of multiverse analysis
(Steegen et al. 2016). We extend this approach to incorporate decisions
relevant to algorithmic fairness, particularly with regard to protected
attributes and the translation of predictions into real-world actions or
interventions. Similarly to a classical multiverse analysis, we use the
resulting _garden of forking paths_ to generate a grid of all possible
universes of decisions combinations, the multiverse. For each of these
universes, we compute the resulting fairness metric of the machine learning
system and collect it as a data point. Based on the resulting dataset of
decision universes and corresponding fairness scores, we evaluate how
individual decisions influence the fairness metric and explore the most
important decisions in more detail (Figure 1).

####  1.2.1 Related Research

Existing work has described the effects of specific pre-processing or modeling
decisions in isolation, such as the influence of different imputation methods
(Caton, Malisetty, and Haas 2022), of the model architecture and
hyperparameters (Sukthanker et al. 2022) on fairness in different contexts.
Multiverse analyses have also been used to model the performance distribution
in hyperparameter-space (Bell et al. 2022), but not yet for analyzing
algorithmic fairness.

The field of hyperparameter-optimization (HPO) (Feurer and Hutter 2019; Bischl
et al. 2023), tries to optimize the process of tuning machine learning model
hyperparameters. This field typically focuses on optimizing algorithm
performance by employing efficient search strategies. These search strategies
allow achieving higher performance without requiring exploration of the
complete hyperparameter space. However, such search patterns usually focus on
finding the optimal configuration and usually yield non-i.i.d. optimization
traces. This makes them unsuitable for assessing the influence and robustness
of any particular decisions. While algorithmic fairness and the importance of
pre-processing are also explored (Perrone et al. 2021), they are typically
only examined on the sidelines. Here, we draw on insights and methodology from
the field of HPO, in particular the functional analysis of variance (FANOVA)
(Hooker 2007; Hutter, Hoos, and Leyton-Brown 2014) to allow a more
interpretable and efficient analysis of the results from the multiverse
analysis. Our focus, however, is on uncovering and systematically exploring
variation induced by the different decisions instead of finding the setting
that optimizes fairness metrics.

###  1.3 Case Study

We illustrate how multiverse analysis can enrich the machine learning fairness
toolkit using a case study of predicting public health insurance coverage.
Accurate and fair prediction of public health insurance coverage in the United
States is an important issue as access to healthcare is quite expensive in the
US, with the country spending almost $16\%$ of its gross domestic product per
capita on healthcare (Ortiz-Ospina and Roser 2017). Whether or not someone is
covered by health insurance has large effects on their health and financial
situation: People with insurance have better self-reported health, have more
preventative doctors appointments, improved depression outcomes, and fewer
personal bankruptcies (Sommers, Gawande, and Baicker 2017).

Given the complexity of US public healthcare system it is easy for people to
fall through its cracks, missing the chance for coverage due to not
understanding how it works or by not fulfilling all necessary requirements. To
combat this, one might want to set up either financial or informational
interventions targeting low-income individuals who are at risk of not being
covered by public health insurance. Determining who and who not to target with
such an intervention could be facilitated with an ADM system. However, given
the vulnerability of the target population and ethical implications of
distributing interventions, fair and well-calibrated predictions are of the
utmost importance in such a scenario.

We implement our case study using the ACSPublicCoverage dataset (Ding et al.
2021). We use this particular dataset as it is rich enough for us to implement
a wide range of design decisions and because many other well-established
datasets used in the fairness literature suffer from non-trivial quality
issues (Ding et al. 2021; Fabris et al. 2022; Bao et al. 2022): UCI Adult
(Kohavi and Becker 1996), the most popular dataset in the fairness literature
(Fabris et al. 2022), uses an arbitrary threshold of 50,000 USD to create a
binary task of income prediction. This threshold has been shown to greatly
influence the accuracy of predictions in certain groups, biasing measures of
algorithmic fairness and threatening external validity (Ding et al. 2021). The
ACSPublicCoverage dataset is one of the datasets which have been specifically
developed in response to the issues in UCI Adult.

Here, we operationalize having public insurance coverage as being covered by
either Medicare, Medicaid, Medical Assistance (or any kind of government-
assistance plan for those with low incomes or a disability) or Veterans
Affairs Health Care, following the official Guidance for Health Insurance Data
Users from the US Census Bureau (Bureau 2021). In line with the original task
setup by Ding et al. (2021), only individuals with an age below 65 years and a
yearly income of less than $30,000 are examined. Low-income households are
also more likely to rely on public health insurance (Keisler-Starkey and Bunch
2022).

As there are no clear guidelines on how to set up an ADM system within this
context (as would be the case in heavily regulated contexts such as credit
scoring) one faces a multitude of decisions when designing a solution for this
task, each of which can govern how bias is fed into the final system. A
multiverse analysis for algorithmic fairness requires developers to make these
design decisions explicit and shows their fairness implications in the present
context.

##  2 Methodology

###  2.1 Fairness Metric

While our proposed analysis works with multiple different fairness metrics, it
requires one to choose a primary metric for analysis. For the present case
study we used _equalized odds difference_ (Agarwal et al. 2018; Hardt et al.
2016) as the primary fairness metric, as it quantifies the degree to which a
system’s predictions are equally good across different groups defined by a
protected attribute. Equalized odds require both the _true positive rate_
(TPR) and the _false positive rate_ (FPR) of a system’s predictions to be
equal across all groups of the protected attribute. Values of the _equalized
odds difference_ can range from $0$ to $1$. A value of $0$ corresponds to a
perfectly fair model according to the metric, whereas a value of $1$
corresponds to a completely unfair model. We use the implementation from the
fairlearn package (Bird et al. 2020) to calculate the metric, where the
differences in both the _true positive rate_ and the _false positive rate_ are
calculated and the larger of the two is used as the metric. We consider _race_
as the protected attribute in our case study given the persisting racial
disparities in various domains, including health outcomes, in the US
(Obermeyer et al. 2019).

###  2.2 Decision Space

When conducting a multiverse analysis, the first step is the identification of
relevant and plausible decisions to be made. Based on the literature on data
science and machine learning workflows (Kuhn and Johnson 2020; Le Quy et al.
2022) we identified five distinct categories to structure and guide the
identification of decisions: Data Selection, Preprocessing, Modeling,
Evaluation and Post-Hoc decisions (Table 1).

For this case study, we considered $10$ distinct and orthogonal design
decisions. Each of these decisions has two to five unique choice options,
leading to a total of $N=122880$ combinations of decisions or universes. An
overview of the decisions and their respective options can be seen in Table 1,
and a detailed description of each is provided below. We consider decisions
roughly in the order they would be made during a typical analysis and sort
them under the list of typical decision categories. As there is a potentially
infinite list of possible decisions to consider, the present list is not
intended to be exhaustive, but rather to highlight the most common and
important categories of decisions one may typically encounter when designing a
machine learning or ADM system. We also deliberately set the focus on
decisions where alternative options are typically not considered or ones that
are not identified as decisions at all. When adapting the methodology to a new
system, this list can serve as an inspiration, however, one must also consider
the domain-specific decisions unique to each applied problem.

Table 1: Overview of the typical decision categories, the actual decisions examined in the case study and their respective options used to construct the multiverse.  | Decisions and Options Examined in Case Study
---|---
Category | Decision |  Options
Data Selection | Exclude Features |  (1) none; (2) race; (3) sex; (4) race-sex
| Exclude Subgroups |  (1) keep-all; (2) drop-smallest-1; (3) drop-smallest-2; (4) keep-largest-2; (5) drop-other
Preprocessing | Scale |  (1) do-not-scale; (2) scale
| Preprocess Age |  (1) none; (2) bins-10; (3) quantiles-3; (4) quantiles-4
| Preprocess Income |  (1) none; (2) bins-10000; (3) quantiles-3; (4) quantiles-4
| Encode Categorical |  (1) one-hot; (2) ordinal
Modeling | Model |  (1) logreg; (2) rf; (3) gbm; (4) elasticnet
Evaluation | Stratify Split |  (1) none; (2) target; (3) protected-attribute; (4) both
| Fairness Grouping |  (1) majority-minority; (2) separate
Post-Hoc | Cutoff |  (1) raw-0.5; (2) quantile-0.1; (3) quantile-0.25

####  2.2.1 Data Selection

#####  2.2.1.1 Excluding Variables as Predictors (Exclude Features)

Selecting features to train a model on presents a critical design decision. In
the ADM context, it can be required to exclude certain protected features
(such as sex/gender, race, ethnicity) as predictors due to legal constraints
when designing a machine-learning system. However, as prominently shown in
various studies this does not necessarily lead to increased fairness, as the
protected attribute is often correlated with other (“legitimate”) features
(Weerts 2021).

We implement the following options for this decision in our case study: (1)
use all features as predictors (incl. protected ones), (2) exclude race, the
protected attribute in the case study, (3) exclude sex, a sensitive attribute
and (4) exclude both race and sex as protected / sensitive attributes.

#####  2.2.1.2 Excluding Subgroups of the Protected Attribute (Exclude
Subgroups)

When working with variables with an uneven distribution or very rare
categories one may focus only on the most common groups, dropping data for
smaller ones. This can be done to preserve the privacy of small groups or out
of convenience to allow for an easier model interpretation downstream.
However, the exclusion of subgroups of the population can potentially be
harmful, with discriminatory differences in downstream model predictions.
While we decided to include this practice as a decision in our analysis to (1)
raise awareness of the issue and (2) represent the effects of the practice in
our analysis, this should not be taken as an endorsement of this practice.

We try to capture the implications of this practice via the attribute race. We
therefore chose to include a decision of dropping certain groups from the
training data based on their prevalence. To accurately compute the fairness
metric, groups were _not_ dropped from the test data used for evaluation.

We include six options for this decision, with the fraction of discarded data
in brackets111Fractions of discarded training data are only reported for a
non-stratified train-test split, as there are only _very slight_ differences
in the fraction of discarded data based on stratification strategy.: (1) to
keep all groups ($0.00\%$), (2) to drop the smallest group ($0.01\%$), (3) to
drop the two smallest groups ($0.33\%$), (4) to keep the two largest groups
($27.45\%$) and (5) to drop the category “Some Other Race alone” specifically
($15.81\%$).

####  2.2.2 Preprocessing

#####  2.2.2.1 Scaling of Continuous Variables (Scale)

It is common to scale continuous variables during preprocessing, centering
them on a mean of $\mu=0$ and standard deviation of $\sigma=1$ (also referred
to as z-scaling). Scaling may be particularly advisable if kernel-based
learners are used as it typically leads to improved performance for such
models.

We include two options for this decision: (1) to keep continuous variables as
they are and (2) to scale continuous variables.

#####  2.2.2.2 Binning of Continuous Variables (Preprocess Age, Preprocess
Income)

Another common practice is binning continuous variables, i.e., turning
continuous variables into ordinal variables with discrete categories. The
reasons to do this are plentiful: To deal with outliers, to address privacy
concerns, or for a more tangible interpretation to name a few.

We provide two distinct and orthogonal decisions here on whether or how to bin
the variables $age$ and $income$. We include four options for the variable
$age$: (1) perform no binning, (2) bin into bins of size $10$, (3) bin into
three evenly sized quantiles, (4) bin into four evenly sized quantiles.
Likewise, we include four options for the variable $income$: (1) perform no
binning, (2) bin into bins of size $10,000$, (3) bin into three evenly sized
quantiles, (4) bin into four evenly sized quantiles.

#####  2.2.2.3 Encoding of Categorical Variables (Encode Categorical)

Another common pre-processing step includes transforming categorical variables
into a numerical format. When doing this one typically has two options: (1)
One-hot (or dummy) coding each variable with $K$ categories into $K$ (or
$K-1$) new binary variables or (2) ordinally encoding each variable by
assigning an integer value from $1$ to $K$ for each category. Ordinal encoding
is only applicable, however, for variables with a natural ordering.

For all ordinal variables (including continuous variables that have been
binned), we include both options. Any variables without a natural ordering are
always one-hot coded.

####  2.2.3 Modeling

#####  2.2.3.1 Model Type (Model)

A major choice when designing any statistical or machine learning system is
which model type one decides to use. While there is a large number of
potential models to explore here, we focused on the most commonly used ones in
the context of ADM in the literature. We note that hyperparameter selection
has shown to have an impact on fairness, but choose to focus on the simple
case, as HPO has already been studied elsewhere (Perrone et al. 2021).

We therefore support the following model types as options for this decision:
(1) logistic regression (Cox 1958), (2) random forest (Ho 1995), (3) gradient
boosted machine (Friedman 2001), and (4) elastic net (Zou and Hastie 2005).

####  2.2.4 Evaluation

#####  2.2.4.1 Stratification of Train-Test Split (Stratify Split)

Training and test sets are often created by simple random splitting of the
full dataset. It can be beneficial, however, to perform this split conditional
on certain groupings to ensure equal representation of all labels within both
the train and test sets.

We include four options for this decision: (1) to not stratify at all, using a
completely random split instead, (2) to stratify using the target variable
(_public coverage_), (3) to stratify using the protected attribute (_race_)
and (4) to stratify using a combination of both variables.

#####  2.2.4.2 Grouping of Protected Attribute (Fairness Grouping)

When working with a fairness metric, it is necessary to specify for which
groups of the protected attribute it is calculated. The present case study
uses _race_ as the protected attribute. For protected attributes with more
than two categories, however, multiple comparisons can be computed. Depending
on the application context one may, e.g., simplify these groups into the
largest group (_majority_) and all other groups (_minority_)222Majority group:
‘White alone’; Minority group: ‘Asian alone’, ‘Two or More Races’, ‘Some Other
Race alone’, ‘Black or African American alone’, ‘American Indian alone’,
‘Native Hawaiian and Other Pacific Islander alone’, ‘American Indian and
Alaska Native tribes specified; or American Indian or Alaska Native, not
specified and no other races’ and ‘Alaska Native alone’..

An important note regarding this decision is that it changes how the fairness
metric is calculated: with two groups, the difference between those two groups
is calculated, however, with more than two groups all possible differences
between group-pairs are calculated and the largest difference between them is
used (the standard procedure for this metric). Naturally, this has a strong
influence on the fairness metric. We therefore conduct many of our later
analyses separately for each of its options.

We include two options for this decision: (1) The fairness metric is computed
between the _majority_ group and _minority_ group and (2) the fairness metric
is computed as the maximum of the metric as computed between all groups of the
protected attribute (_race_)333This corresponds to the default behavior in the
fairlearn library..

####  2.2.5 Post-Hoc

#####  2.2.5.1 Cutoff for Final Classification (Cutoff)

At the end of the ML pipeline, the prediction models’ (risk) scores can be
used to classify new observations based on a pre-specified classification
threshold. By default a threshold of $0.5$ would be used with every score
equal or above classified as $1$ (_having coverage_) and everything below as
$0$ (_not having coverage_). Actual interventions, however, are often based on
the ranked list of scores such that (costly) interventions are targeted at the
top $X$ percent with the highest risk. With real-world scenarios often coming
with resource-bound restrictions, one may for example only be able to provide
an intervention for, say, $10\%$ or $25\%$ of the most in-need in the
population. These real-world restrictions are typically not taken into account
in fairness evaluations, despite having potentially devastating implications.

We therefore also consider different cutoff values for the final predictions
of the system. We support the following options for this decision: (1) use the
default raw cutoff value of $0.5$, (2) only treat the lowest $0.1$ quantile as
_not having coverage_ , (2) only treat the lowest $0.25$ quantile as _not
having coverage_.

###  2.3 Analysis

We examined the overall variation of the fairness metric and the degree to
which it is explained by the outlined decisions. To estimate the importance of
individual decisions and their interactions, we used a functional analysis of
variance (FANOVA) (Hooker 2007; Hutter, Hoos, and Leyton-Brown 2014). Based on
the results of the FANOVA, we examined the most important decisions and
interactions in detail.

As it can be costly to iterate over the complete multiverse grid, we also
examined the feasibility of running our analyses on smaller subsets of the
data and comparing results with analysis performed on the full dataset. The
implementation for running a FANOVA by Hutter, Hoos, and Leyton-Brown (2014)
used here has also been demonstrated to work well for assessing hyperparameter
importance with only a subset of data from the hyperparameter grid available.

###  2.4 Technology

Analyses were conducted using Python Version 3.8 (Van Rossum and Drake 2009)
and pipenv (P. M. Team 2017) for reproducibility. The Python package scikit-
learn (Pedregosa et al. 2011) was used for pre-processing and fitting of
models, pandas (team 2020) for loading and modification of data, folktables
(Ding et al. 2021) for retrieval of data, fairlearn (Bird et al. 2020) for
computation of fairness metrics, fANOVA (Hutter, Hoos, and Leyton-Brown 2014)
for calculation of variable importance and papermill (contributors 2017) for
parameterized computation of decision universes. This reproducible document
was generated using quarto (Allaire et al. 2022), R (R. C. Team 2022) Version
4.2, the R packages from the tidyverse (Wickham et al. 2019) and ggpubr
(Kassambara 2023) for generation of figures. The source code of the analyses
and this publication is available at https://github.com/reliable-ai/fairml-
multiverse.

##  3 Results

###  3.1 Distribution of Metrics

The multiverse analysis in our case study produced a total of $N=122880$
values of the fairness metric across all possible analysis choices. When
examining the distribution of the fairness metric across the multiverse of
decisions, the large variation of the fairness metric becomes apparent, with
values spanning the entire possible range of the metric from $0$ to $1$
(Figure 2). Overall performance of the resulting models was moderate with raw
accuracies between $0.419$ and $0.722$. Performance and the fairness metric
were almost uncorrelated with a Pearson correlation of $r=0.082$. Raw
performance varied largely based on the decision _Cutoff_ , with three large
clusters of similar performance (Figure 3). Within these three clusters of
almost equal performance there was a large variance of the fairness metric,
highlighting the opportunity to optimize algorithmic fairness without
sacrificing performance.

Figure 2: Using a different grouping of the protected attribute strongly
influences the fairness metric. Distribution of fairness metric (equalized
odds difference) split by grouping of the protected attribute with vertical
lines corresponding to mean values. Lower values on the fairness metric
indicate smaller _TPR_ and _FPR_ differences across groups. Figure 3:
Performance and fairness are largely unrelated with three plateaus of low
variance in performance, but high variance in fairness. Distribution of
overall performance (raw accuracy) and fairness metric (equalized odds
difference) across all multiverses. The marginal histogram shows performance
for different options of the _Cutoff_ decision. A marginal histogram of the
fairness metric can be seen in Figure 2.

The fairness metric is calculated differently depending on how the groupings
of the protected attribute are made: for the _majority-minority_ grouping only
a single comparison can be made, whereas the largest value from multiple
comparisons is used with the _separate_ grouping. Because of this, we examined
this decision first and in isolation. The implications of the two groupings on
the fairness metric are visible in Figure 2: The _majority-minority_ grouping
leads to a roughly normal distribution of low values on the fairness metric,
centered around a mean of $M=0.047$ ($SD=0.030$). The _separate_ grouping on
the other hand leads to much wider distribution of fairness values centered
around a mean value of $M=0.533$ ($SD=0.291$). Since the decision has such a
strong influence on the fairness metric and interacts with many of the other
decisions, we chose to conduct further analyses only within these two
groupings to differentiate results more clearly and allow for easier
interpretation as well as visualization. The following analyses are therefore
conducted once for the _majority-minority_ grouping and once for the
_separate_ grouping of the protected attribute; the sample size in each arm of
the following analyses is thus $n=61440$.

###  3.2 Importance of Decisions

Table 2: The 10 most important decisions or decision interactions and their
relative importance for both groupings of the protected attribute.

(a) Most important decisions for majority-minority grouping.  Effect Type | Decision / Interaction of Decisions | Importance | Std. Deviation
---|---|---|---
main | $\mathit{ExcludeSubgroups}$ | 0.254 | 0.001
main | $\mathit{Cutoff}$ | 0.211 | 0.001
2-way int. | $\mathit{ExcludeFeatures}\times\mathit{ExcludeSubgroups}$ | 0.110 | 0.001
2-way int. | $\mathit{Cutoff}\times\mathit{ExcludeSubgroups}$ | 0.082 | 0.000
3-way int. | $\mathit{Cutoff}\times\mathit{ExcludeFeatures}\times\mathit{ExcludeSubgroups}$ | 0.027 | 0.000
2-way int. | $\mathit{Model}\times\mathit{PreprocessIncome}$ | 0.020 | 0.000
main | $\mathit{Model}$ | 0.017 | 0.000
2-way int. | $\mathit{ExcludeSubgroups}\times\mathit{Model}$ | 0.016 | 0.000
main | $\mathit{ExcludeFeatures}$ | 0.014 | 0.000
3-way int. | $\mathit{Model}\times\mathit{PreprocessIncome}\times\mathit{Scale}$ | 0.012 | 0.000

(b) Most important decisions for separate grouping.  Effect Type | Decision / Interaction of Decisions | Importance | Std. Deviation
---|---|---|---
main | $\mathit{StratifySplit}$ | 0.375 | 0.001
2-way int. | $\mathit{Cutoff}\times\mathit{StratifySplit}$ | 0.313 | 0.001
main | $\mathit{Cutoff}$ | 0.082 | 0.000
4-way int. | $\mathit{Cutoff}\times\mathit{ExcludeFeatures}\times\mathit{Model}\times\mathit{StratifySplit}$ | 0.007 | 0.000
3-way int. | $\mathit{Cutoff}\times\mathit{Model}\times\mathit{StratifySplit}$ | 0.007 | 0.000
3-way int. | $\mathit{Cutoff}\times\mathit{Model}\times\mathit{PreprocessIncome}$ | 0.007 | 0.000
2-way int. | $\mathit{Model}\times\mathit{PreprocessIncome}$ | 0.006 | 0.000
2-way int. | $\mathit{ExcludeFeatures}\times\mathit{Model}$ | 0.006 | 0.000
3-way int. | $\mathit{Model}\times\mathit{PreprocessIncome}\times\mathit{Scale}$ | 0.006 | 0.000
2-way int. | $\mathit{Cutoff}\times\mathit{PreprocessIncome}$ | 0.005 | 0.000

We conducted two FANOVAs (Hooker 2007) as described in Hutter, Hoos, and
Leyton-Brown (2014) to assess the importance of decisions on the fairness
metric. This analysis decomposes the overall variance of the fairness metric
into the fractions which are explained by each decision. These variance
decompositions are used to assess the relative importance of decisions.
Moreover, the FANOVA also allows computing explained variance for interactions
of decisions. This is highly useful, as the overall interaction space between
decisions is quite large with $511$ possible (interaction and main) effects.

Using the resulting importance values from the FANOVA, one can see which
decisions are associated with a high variation in fairness scores, whether it
be by themselves or in conjunction with others. This allows assessing the most
consequential decisions on a one-by-one case. Table 2 contains a ranked list
of the most important decisions and decision interactions in our case study
alongside their respective importance. When using a _majority-minorit_ y
grouping the most important decision is whether and which subgroups one
chooses to exclude during model training, whereas the most important decisions
for the _separate_ grouping is how the stratification of the train-test split
is performed. The cutoff value used for the final predictions is important in
both cases and often has effects in conjunction with other decisions as well.
Specifically the interaction of the chosen cutoff value with the
stratification strategy is highly important when using the _separate_
grouping, accounting for more than $30\%$ of the variance in the fairness
metric.

#####  3.2.0.1 Examining Individual Decisions

We analyzed the three most important decisions or decision-interactions per
grouping approach to further illustrate the methodology and how one would
explore the results of the analysis. The results also highlight why one should
investigate the decisions in a detailed manner and not just pick the most-fair
and highest-performing universe’s model.

For the _majority-minority_ grouping, _Exclude Subgroups_ was the most
influential decision. Examining the decision leads to surprising results: One
might expect that the exclusion of any subgroups in the training data will
reduce fairness (as measured by the equalized odds difference), however,
results show that dropping the “Other” group from the training data actually
lead to slightly fairer models on average (Figure 4 A). Whether one should
actually drop the “Other” group from the training data should still be
evaluated very carefully, however. The least fair models were produced when
retaining only the largest two groups of the protected attribute and dropping
all others. The second most important decision was which _Cutoff_ value to
choose when making the final predictions. Here, the default value of $0.5$
typically lead to higher values of the fairness metric (Figure 4 B). It is
important, to be aware of the importance of this decision, as one might be
forced by practical circumstances to use a different cutoff after deployment
in which case model fairness would need to be re-evaluated. Last, examining
the interplay between _Exclude Subgroups_ and _Exclude Features_ in Figure 4
(C), illustrates how the effect of keeping only the two largest groups of race
or dropping the “Other” group is greatly amplified when race is included as a
predictive feature.

For the _separate_ grouping, the effects are slightly more indistinct as
computing the maximum fairness metric between all comparisons leads to more
volatile results. The decisions _Stratify Split_ , _Cutoff_ and their
interaction account for all three of the most important decisions with this
grouping. When examining the decision separately, it can be seen how
stratifying by the target variable leads to noticeably less fair models
(Figure 4 D, most important) and how the raw cutoff value of $0.5$ is suddenly
not leading to the most fair models anymore (Figure 4 F, third most
important). The effects of both variables become most clear, however, when
examining their interaction, which was identified as explaining almost as much
variance as the most important decision. While using a cutoff value
corresponding to the top $10\%$ quantile leads to the least fair model when
stratifying by the target variable it surprisingly leads to the fairest models
when using any other stratification strategy (Figure 4 E, second most
important).

Figure 4: Detailed analysis of most important decisions. Visualization of the
fairness metric depending on the three most important decision / decision
combinations and their respective options for _majority-minority_ (A-C) and
_separate_ (D-F) groupings of the protected attribute.

###  3.3 Scaling Multiverse Analysis for Algorithmic Fairness

To assess the feasibility of running the multiverse analysis on a smaller
subset of the grid, we also conducted the FANOVAs on different subsamples of
the collected _multiverse_ dataset. Specifically, we ran the analysis on
random subsets of $1\%$, $5\%$, $10\%$ and $20\%$ of the data and calculated
the correlation of variance decomposition or importance values with the FANOVA
estimated on the full multiverse dataset. The estimates of variance
decomposition are highly skewed, with a few highly important decisions and a
very larger number of very low-importance decisions. We therefore calculated
both, the Pearson correlation which is more sensitive to correlations of the
more important decisions and the Spearman rank-correlation which is also
sensitive to decisions with low importance estimates. To assess the
consistency of this approach we computed the FANOVA on each subsample 50 times
and calculated the correlation with the results from the full _multiverse_
dataset every time.

When calculating the Pearson correlation, the resulting mean correlation
coefficient ranged from $\bar{r}_{1\%}=0.989$ ($SD=0.004$) at $1\%$ to
$\bar{r}_{20\%}\geq 0.999$ ($SD=0$) at $20\%$ for the _majority-minority_
grouping and $\bar{r}_{1\%}=0.997$ ($SD=0.002$) at $1\%$ to
$\bar{r}_{20\%}\geq 0.999$ ($SD=0$) at $20\%$ for the _separate_ grouping.
Spearman rank-correlations were also high, but lower than the Pearson
correlation coefficients and more inconsistent (Figure 5), which indicates
that using sparse data to estimate the importance of decisions works well for
the most important decisions and less-so to identify nuances between the less-
important decisions. The resulting Spearman rank-correlation mean coefficients
ranged from $\bar{\rho}_{1\%}=0.517$ ($SD=0.032$) at $1\%$ to
$\bar{\rho}_{20\%}=0.860$ ($SD=0.012$) at $20\%$ for the _majority-minority_
grouping and $\bar{\rho}_{1\%}=0.530$ ($SD=0.03$) at $1\%$ to
$\bar{\rho}_{20\%}=0.937$ ($SD=0.007$) at $20\%$ for the _separate_ grouping.

Figure 5: Conducting the analysis with smaller subsets of the complete
multiverse leads to similar results. Correlations of variance decomposition /
importance estimates between full dataset and random subsets of different
sizes. Random subsets were drawn 50 times with points corresponding to mean
correlations and lines to +/- 1 standard deviation. Pearson correlation
coefficients are consistently higher than Spearman correlation coefficients,
indicating better estimation of high-importance decisions.

##  4 Discussion

We demonstrate how multiverse analysis for algorithmic fairness provides a
useful new method for evaluating the robustness of machine learning and ADM
systems with respect to decisions along the modeling pipeline and their
implications for algorithmic fairness. Our method provides a promising new
methodology which can empower analysts to better understand how their
decisions affect a system’s fairness and which decisions matter. We highlight
the importance of making decisions during model design explicitly rather than
implicitly.

By applying this new methodology in a use case of predicting public health
care coverage, we demonstrate the feasibility of this approach and show which
decisions from our list of decisions affect fairness the most: Unsurprisingly
the grouping of the protected attribute used for calculating the fairness
metric has considerable effect on algorithmic fairness. More interestingly,
this grouping does also affect the influence of almost all other decisions.
Besides the grouping, we showed that the cutoff value used for making final
decisions has a significant effect on the fairness metric, a decision often
implemented post-hoc after model deployment without any consideration of
fairness. We also observe that the exclusion of certain subgroups of the
protected attribute during training affects fairness downstream, especially
when the protected attribute is kept as a feature. Surprisingly, we also saw
that the stratification strategy used for the train-test split had strong
effects on the fairness metric.

When interpreting the results from a multiverse analysis for algorithmic
fairness, one should evaluate results with care and strictly avoid merely
selecting the combination of decisions with the best fairness metric. Results
should be seen as an indication of how susceptible the fairness of the model
is to design decisions and which decisions warrant closer examination. Results
from the analysis can also be used to guide the search of new options for the
most important decisions. Final choices regarding the design of the system
should be made using a combination of empirical results from the multiverse
analysis and practical as well as ethical considerations within the context of
the use case. This can often mean that decisions may be made that do not
correspond to the optimal value in the fairness metric for a decision. The
main goal of a multiverse analysis for algorithmic fairness is to facilitate
making educated and explicit decisions. We recommend including complete
results from the analysis alongside the final system.

As we explored only a single use-case, we do not make any generalizable claims
regarding the importance of any particular decisions, beyond the fact that
these decisions _can_ matter and are worth investigating. Another limitation
of this case-study is that we only examined 10 distinct decisions, with many
plausible alternative decisions which could’ve been examined in their place or
additionally. As there is an infinite space of decisions one may consider, we
decided to draw the line at these 10 decisions for illustrative purposes. A
successful adoption of multiverse analysis for algorithmic fairness in
different use cases and reporting of results could help identify a more
exhaustive list of the most important decisions across contexts. Potential
concerns regarding the computational cost of conducting a multiverse analysis
for algorithmic fairness are valid, but can be addressed as we demonstrate
that estimates of importance are robustly detected for important decisions
even when exploring only $1\%$ of the full _multiverse_.

We encourage the use of the method during the design of future machine
learning or ADM systems and provide an overview of the most important areas of
decisions to guide analysts when adapting multiverse analysis for algorithmic
fairness in their own context. We further provide a non-exhaustive list of
exemplary decisions to serve as inspiration to identify potentially relevant
decisions and a modular implementation that makes adoption to different use-
cases easy. We posit that results from a multiverse analysis for algorithmic
fairness can critically inform discussions between developers and stakeholders
and advise joint reflections on the ultimate design of ADM systems.

By successfully adapting multiverse analysis across disciplines we also
highlight the feasibility of adapting methods and techniques between
disciplines and hope that this will inspire further cross-discipline
pollination of ideas and methodologies between research fields.

##  5 References

reAgarwal, Alekh, Alina Beygelzimer, Miroslav Dudík, John Langford, and Hanna
Wallach. 2018. “A Reductions Approach to Fair Classification.”

preAgrawal, Ashrya, Florian Pfisterer, Bernd Bischl, Francois Buet-Golfouse,
Srijan Sood, Jiahao Chen, Sameena Shah, and Sebastian Vollmer. 2021.
“Debiasing Classifiers: Is Reality at Variance with Expectation?”
https://doi.org/10.48550/arXiv.2011.02407.

preAllaire, J. J., Charles Teague, Carlos Scheidegger, Yihui Xie, and
Christophe Dervieux. 2022. _Quarto_. https://doi.org/10.5281/zenodo.5960048.

preAngwin, Julia, Jeff Larson, Surya Mattu, and Lauren Kirchner. 2016.
“Machine Bias.” _ProPublica_ , May, 254264.
https://www.propublica.org/article/machine-bias-risk-assessments-in-criminal-
sentencing.

preBao, Michelle, Angela Zhou, Samantha Zottola, Brian Brubach, Sarah
Desmarais, Aaron Horowitz, Kristian Lum, and Suresh Venkatasubramanian. 2022\.
“It’s COMPASlicated: The Messy Relationship Between RAI Datasets and
Algorithmic Fairness Benchmarks.” https://doi.org/10.48550/arXiv.2106.05498.

preBell, Samuel J., Onno P. Kampman, Jesse Dodge, and Neil D. Lawrence. 2022\.
“Modeling the Machine Learning Multiverse.”
https://doi.org/https://doi.org/10.48550/arXiv.2206.05985.

preBird, Sarah, Miroslav Dudík, Richard Edgar, Brandon Horn, Roman Lutz,
Vanessa Milan, Mehrnoosh Sameki, Hanna Wallach, Kathleen Walker, and Allovus
Design. 2020. “Fairlearn: A Toolkit for Assessing and Improving Fairness in
AI.”

preBischl, Bernd, Martin Binder, Michel Lang, Tobias Pielok, Jakob Richter,
Stefan Coors, Janek Thomas, et al. 2023. “Hyperparameter Optimization:
Foundations, Algorithms, Best Practices, and Open Challenges.” _WIREs Data
Mining and Knowledge Discovery_ 13 (2). https://doi.org/10.1002/widm.1484.

preBlack, Emily, Manish Raghavan, and Solon Barocas. 2022. “Model
Multiplicity: Opportunities, Concerns, and Solutions.”

preBreznau, Nate, Eike Mark Rinke, Alexander Wuttke, Hung H. V. Nguyen, Muna
Adem, Jule Adriaans, Amalia Alvarez-Benjumea, et al. 2022. “Observing Many
Researchers Using the Same Data and Hypothesis Reveals a Hidden Universe of
Uncertainty.” _Proceedings of the National Academy of Sciences_ 119 (44):
e2203150119. https://doi.org/10.1073/pnas.2203150119.

preBureau, US Census. 2021. “ACS Health Insurance Coverage Recoding
Programming Code.” https://www.census.gov/topics/health/health-
insurance/guidance/programming-code/acs-recoding.html.

preCaton, Simon, Saiteja Malisetty, and Christian Haas. 2022. “Impact of
Imputation Strategies on Fairness in Machine Learning.” _Journal of Artificial
Intelligence Research_ 74 (September). https://doi.org/10.1613/jair.1.13197.

precontributors, nteract. 2017. _Papermill: Parametrize and Run Jupyter and
Nteract Notebooks_. https://github.com/nteract/papermill.

preCox, David R. 1958. “The Regression Analysis of Binary Sequences.” _Journal
of the Royal Statistical Society: Series B (Methodological)_ 20 (2): 215232.

preDing, Frances, Moritz Hardt, John Miller, and Ludwig Schmidt. 2021.
“Retiring Adult: New Datasets for Fair Machine Learning,” 13.

preFabris, Alessandro, Stefano Messina, Gianmaria Silvello, and Gian Antonio
Susto. 2022. “Algorithmic Fairness Datasets: The Story so Far.” _Data Mining
and Knowledge Discovery_ , September.
https://doi.org/10.1007/s10618-022-00854-z.

preFaliagka, Evanthia, Kostas Ramantas, and Giannis Tzimas. 2012. “Application
of Machine Learning Algorithms to an Online Recruitment System.”

preFeurer, Matthias, and Frank Hutter. 2019. “Hyperparameter Optimization.”
In, edited by Frank Hutter, Lars Kotthoff, and Joaquin Vanschoren, 3–33. The
Springer Series on Challenges in Machine Learning. Cham: Springer
International Publishing. https://doi.org/10.1007/978-3-030-05318-5_1.

preFriedman, Jerome H. 2001. “Greedy Function Approximation: A Gradient
Boosting Machine.” _The Annals of Statistics_ 29 (5): 1189–1232.
https://doi.org/10.1214/aos/1013203451.

preGelman, Andrew, and Eric Loken. 2014. “The Statistical Crisis in Science.”
_American Scientist_ 102 (6): 460.

preHardt, Moritz, Eric Price, Eric Price, and Nati Srebro. 2016. “Equality of
Opportunity in Supervised Learning.”

preHenriques-Gomes, Luke. 2023. “Robodebt: Five Years of Lies, Mistakes and
Failures That Caused a $1.8bn Scandal.” _The Guardian_ , March.
https://www.theguardian.com/australia-news/2023/mar/11/robodebt-five-years-of-
lies-mistakes-and-failures-that-caused-a-18bn-scandal.

preHo, Tin Kam. 1995. “Random Decision Forests.” In, 1:278282. IEEE.

preHooker, Giles. 2007. “Generalized Functional ANOVA Diagnostics for High-
Dimensional Functions of Dependent Variables.” _Journal of Computational and
Graphical Statistics_ 16 (3): 709–32. https://www.jstor.org/stable/27594267.

preHutter, Frank, Holger Hoos, and Kevin Leyton-Brown. 2014. “International
Conference on Machine Learning.” In, 754–62. PMLR.
https://proceedings.mlr.press/v32/hutter14.html.

preInternational, Amnesty. 2021. “Xenophobic Machines.”
https://www.amnesty.org/en/wp-
content/uploads/2021/10/EUR3546862021ENGLISH.pdf.

preKassambara, Alboukadel. 2023. _Ggpubr: ’Ggplot2’ Based Publication Ready
Plots_. https://CRAN.R-project.org/package=ggpubr.

preKeisler-Starkey, Katherine, and Lisa N Bunch. 2022. “Health Insurance
Coverage in the United States: 2021 - Appendix Table C3.”
https://www.census.gov/content/dam/Census/library/publications/2022/demo/p60-278.pdf.

preKern, Christoph, Ruben L. Bach, Hannah Mautner, and Frauke Kreuter. 2021\.
“Fairness in Algorithmic Profiling: A German Case Study.”
https://doi.org/10.48550/arXiv.2108.04134.

preKohavi, Ronny, and Barry Becker. 1996. “Adult Data Set.” _UCI Machine
Learning Repository_ 5: 2093.

preKuhn, Max, and Kjell Johnson. 2020. _Feature engineering and selection: a
practical approach for predictive models_. Chapman & Hall/CRC data science
series. Boca Raton London New York: CRC Press, Taylor & Francis Group.
www.feat.engineering.

preLe Quy, Tai, Arjun Roy, Vasileios Iosifidis, Wenbin Zhang, and Eirini
Ntoutsi. 2022. “A Survey on Datasets for Fairness-Aware Machine Learning.”
_WIREs Data Mining and Knowledge Discovery_ 12 (3): e1452.
https://doi.org/10.1002/widm.1452.

preMehrabi, Ninareh, Fred Morstatter, Nripsuta Saxena, Kristina Lerman, and
Aram Galstyan. 2021. “A Survey on Bias and Fairness in Machine Learning.” _ACM
Computing Surveys_ 54 (6): 115:1115:35. https://doi.org/10.1145/3457607.

preMukerjee, Amitabha, Rita Biswas, Kalyanmoy Deb, and Amrit P. Mathur. 2002\.
“Multi–objective Evolutionary Algorithms for the Risk–return Trade–off in Bank
Loan Management.” _International Transactions in Operational Research_ 9 (5):
583–97. https://doi.org/10.1111/1475-3995.00375.

preObermeyer, Ziad, Brian Powers, Christine Vogeli, and Sendhil Mullainathan.
2019. “Dissecting Racial Bias in an Algorithm Used to Manage the Health of
Populations.” _Science_ 366 (6464): 447–53.
https://doi.org/10.1126/science.aax2342.

preOPEN SCIENCE COLLABORATION. 2015. “Estimating the Reproducibility of
Psychological Science.” _Science_ 349 (6251): aac4716.
https://doi.org/10.1126/science.aac4716.

preOrtiz-Ospina, Esteban, and Max Roser. 2017. “Healthcare Spending.” _Our
World in Data_ , June. https://ourworldindata.org/financing-healthcare.

prePedregosa, F., G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel,
M. Blondel, et al. 2011. “Scikit-Learn: Machine Learning in Python.” _Journal
of Machine Learning Research_ 12: 28252830.

prePerrone, Valerio, Michele Donini, Muhammad Bilal Zafar, Robin Schmucker,
Krishnaram Kenthapadi, and Cédric Archambeau. 2021. “AIES ’21: AAAI/ACM
Conference on AI, Ethics, and Society.” In, 854–63. Virtual Event USA: ACM.
https://doi.org/10.1145/3461702.3462629.

preRodolfa, Kit T., Pedro Saleiro, and Rayid Ghani. 2020. “Bias and Fairness.”
In, 2nd ed. Chapman; Hall/CRC.

preSimmons, Joseph P., Leif D. Nelson, and Uri Simonsohn. 2011. “False-
Positive Psychology: Undisclosed Flexibility in Data Collection and Analysis
Allows Presenting Anything as Significant.” _Psychological Science_ 22 (11):
1359–66. https://doi.org/10.1177/0956797611417632.

preSimonsohn, Uri. 2018. “Two Lines: A Valid Alternative to the Invalid
Testing of U-Shaped Relationships With Quadratic Regressions.” _Advances in
Methods and Practices in Psychological Science_ 1 (4): 538–55.
https://doi.org/10.1177/2515245918805755.

preSommers, Benjamin D., Atul A. Gawande, and Katherine Baicker. 2017. “Health
Insurance Coverage and Health — What the Recent Evidence Tells Us.” _New
England Journal of Medicine_ 377 (6): 586–93.
https://doi.org/10.1056/NEJMsb1706645.

preSteegen, Sara, Francis Tuerlinckx, Andrew Gelman, and Wolf Vanpaemel.
2016\. “Increasing Transparency Through a Multiverse Analysis.” _Perspectives
on Psychological Science_ 11 (5): 702–12.
https://doi.org/10.1177/1745691616658637.

preSukthanker, Rhea, Samuel Dooley, John P. Dickerson, Colin White, Frank
Hutter, and Micah Goldblum. 2022. “On the Importance of Architectures and
Hyperparameters for Fairness in Face Recognition.”
https://doi.org/10.48550/arXiv.2210.09943.

preTeam, Pipenv Maintainer. 2017. _Pipenv: Python Development Workflow for
Humans._ https://github.com/pypa/pipenv.

preTeam, R Core. 2022. _R: A Language and Environment for Statistical
Computing_. Vienna, Austria: R Foundation for Statistical Computing.
https://www.R-project.org/.

preteam, The pandas development. 2020. _Pandas-Dev/Pandas: Pandas_. Zenodo.
https://doi.org/10.5281/zenodo.3509134.

preVan Rossum, Guido, and Fred L. Drake. 2009. _Python 3 Reference Manual_.
Scotts Valley, CA: CreateSpace.

preWeerts, Hilde J. P. 2021. “An Introduction to Algorithmic Fairness.”
https://doi.org/10.48550/arXiv.2105.05595.

preWickham, Hadley, Mara Averick, Jennifer Bryan, Winston Chang, Lucy McGowan,
Romain François, Garrett Grolemund, et al. 2019. “Welcome to the Tidyverse.”
https://joss.theoj.org.

preZou, Hui, and Trevor Hastie. 2005. “Regularization and Variable Selection
via the Elastic Net.” _Journal of the Royal Statistical Society. Series B
(Statistical Methodology)_ 67 (2): 301–20.
https://www.jstor.org/stable/3647580.

p


# Deforming ice with drops

Duco van Buuren  Physics of Fluids group, Max Planck Center Twente for Complex
Fluid Dynamics, Department of Science and Technology, Mesa+ Institute and J.
M. Burgers Center for Fluid Dynamics, University of Twente, P.O. Box 217, 7500
AE Enschede, The Netherlands  Pallav Kant<EMAIL_ADDRESS>School of
Engineering, University of Manchester, M13 9PL, United Kingdom  Jochem G.
Meijer  Physics of Fluids group, Max Planck Center Twente for Complex Fluid
Dynamics, Department of Science and Technology, Mesa+ Institute and J. M.
Burgers Center for Fluid Dynamics, University of Twente, P.O. Box 217, 7500 AE
Enschede, The Netherlands  Christian Diddens  Physics of Fluids group, Max
Planck Center Twente for Complex Fluid Dynamics, Department of Science and
Technology, Mesa+ Institute and J. M. Burgers Center for Fluid Dynamics,
University of Twente, P.O. Box 217, 7500 AE Enschede, The Netherlands  Detlef
Lohse<EMAIL_ADDRESS>Physics of Fluids group, Max Planck Center Twente
for Complex Fluid Dynamics, Department of Science and Technology, Mesa+
Institute and J. M. Burgers Center for Fluid Dynamics, University of Twente,
P.O. Box 217, 7500 AE Enschede, The Netherlands  Max Planck Institute for
Dynamics and Self-Organization, Am Fassberg 17, 37077 Göttingen, Germany

###### Abstract

A uniform solidification front undergoes non-trivial deformations when
encountering an insoluble dispersed particle in a melt. For solid particles,
the overall deformation characteristics are primarily dictated by heat
transfer between the particle and the surrounding, remaining unaffected by the
rate of approach of the solidification front. In this Letter we show that,
conversely, when interacting with a droplet or a bubble, the deformation
behavior exhibits entirely different and unexpected behavior. It arises from
an interfacial dynamics which is specific to particles with free interfaces,
namely thermal Marangoni forces. Our study employs a combination of
experiments, theory, and numerical simulations to investigate the interaction
between the droplet and the freezing front and unveils its surprising
behavior. In particular, we quantitatively understand the dependence of the
front deformation $\Delta$ on the front propagation velocity, which, for large
front velocities, can even revert from attraction ($\Delta<0$) to repulsion
($\Delta>0$).

During the solidification of an aqueous suspension, the dispersed objects
(particles, drops, bubbles, etc.) either become an integral part of the
solidified bulk or get rejected by it (Shangguan  _et al._ , 1992; Rempel and
Worster, 2001; Park  _et al._ , 2006; Dedovets  _et al._ , 2018; Tyagi  _et
al._ , 2020; Meijer  _et al._ , 2023a). The prevailing outcome therefore
controls the micro-structure of the solidified material and thus its
functional properties. Crucially, the conditions dictating the engulfment or
rejection of a dispersed particle rely on a delicate interaction between the
particle/drop/bubble and the moving solidification front, where the relevant
physical processes act across various length scales. At length scales much
smaller than the size of the dispersed particle, this interaction is
controlled by a combination of repulsive van der Waals interactions between
the dispersed particle and the freezing interface and fluid flow within thin
premelted films enveloping the particle Wettlaufer and Worster (2006); Rempel
and Worster (2001); Meijer  _et al._ (2023b, a). In contrast, at the particle-
scale, it is the heat exchange between the particle and the surrounding
liquid, the interfacial energy of the solidification interface, and the
density change caused by the liquid–solid phase transition that controls the
dynamics. Previous investigations have revealed that, as a solidification
front approaches a dispersed particle, it experiences deformation in a non-
trivial manner Shangguan  _et al._ (1992); Tyagi  _et al._ (2020). It is
largely accepted that this deformation of the solidification front is mainly
controlled by the mismatch in thermal conductivity between that of the
particle $k_{p}$ and the surrounding liquid $k_{l}$ which, in turn, dictates
the heat transfer between the two Shangguan  _et al._ (1992). Most
importantly, in case of solid particles, the rate of approach of the
solidification front has no effect on the deformation of the solidification
front (see Fig. 1(a)) (Shangguan  _et al._ , 1992; Tyagi  _et al._ , 2020).

Figure 1: Deformation of the ice front by a solid particle. (a) Shape of the
solidification front approaching a ceramic and polystyrene (PS) particle at
different advancing velocities $V$. Dependent on the thermal conductivity
mismatch, the front is either deflected away (left, $k_{p}/k_{l}=4$) or
towards (right, $k_{p}/k_{l}=0.3$) the particle. The direction of the
deflection is independent on $V$. The scale-bars are
$50\text{\,}\mathrm{\SIUnitSymbolMicro m}$. (b) Schematic of the origin of the
front deflection for particles more (left) or less (right) conductive than the
surrounding liquid when an uniform temperature gradient is applied. The
solidification front corresponds to the isotherm (red lines) that matches the
melting temperature.

In contrast, our Letter introduces a remarkable phenomena occurring when an
immiscible oil droplet immersed in water encounters an advancing
solidification front: We demonstrate that, while the deformation of the
freezing front at low advancing velocities conforms to the expected behaviour
for a given thermal conductivity mismatch, i.e., deflection towards the drop
since $k_{p}/k_{l}<1$ (see Fig. 1(b)), an unexpected reversal of the ice front
deflection occurs at faster freezing rates. We employ a comprehensive
approach, integrating controlled experiments, theoretical insights, and
numerical simulations to elucidate the origin of the observed reversal. Our
analysis reveals that this reversal is attributed to a self-induced flow along
the droplet interface, driven by temperature-induced gradients in surface
tension - a phenomenon known as the thermal Marangoni effect. Such a fluid
flow transports warmer water toward the colder ice front, leading to localized
increase of the temperature field and thus delayed freezing. Since the shape
of the interface plays a crucial role in determining the critical engulfment
velocity, beyond which particles become trapped in the solidifying material
Park  _et al._ (2006); Meijer  _et al._ (2023a), our findings contribute novel
insights also to this aspect. Consequently, this study holds significance not
only in terms of fundamental understanding but also for practical
applications. Note that the deformation of the freezing interface is rather
distinct - both in underlying physical origin and in phenomenology - from the
deformation of the oil droplet during engulfment into the solidifying bulk,
which we presented in our earlier work Ref. Meijer  _et al._ (2023b).

The experimental setup used to investigate the deformation of the advancing
solidification front by dispersed droplets is similar to the one used in our
previous investigations Meijer  _et al._ (2023b, a, c). For brevity, here, we
highlight only the important details for the setup and refer to Ref. Meijer
_et al._ (2023b) for detailed descriptions. We investigate the interaction
between an advancing water-ice freezing front and silicone oil droplets in a
horizontal Hele-Shaw cell of thickness 200 $\mu$m, filled with an oil-in-water
emulsion. The emulsion is prepared using a microfluidic setup, as explained in
Ref. Van Elburg  _et al._ (2021). To ensure the stability of the emulsion, a
surfactant TWEEN80 (Sigma-Aldrich, Germany) is added to the water (0.5 vol%
initially). Since surfactants are known to affect freezing of complex liquids
Dedovets  _et al._ (2018); Kao  _et al._ (2009); Tyagi  _et al._ (2022), the
emulsion sample is diluted after fabrication to achieve low surfactant
concentration (0.01 vol%). As a consequence, any such solutal effects are
irrelevant for our experiments. In each experiment, the Hele-Shaw cell is
translated at constant velocity $V$ through a fixed temperature gradient
$G=\nabla T$ using a linear actuator (Physik Instrumente, M-230.25) such that
the position of the moving solidification front remains fixed in the lab frame
of reference. The temperature gradient $G\sim 10\,\mathrm{K\,cm^{-1}}$ is
applied over the entire length of the Hele-Shaw cell by placing it over a
copper block with cold and warm ends, maintained at fixed temperatures using
Peltier elements. The applied thermal gradient sets the ice growth rate and
also thus the approach velocity $V$ of a droplet to the front. Features of
interest along the moving freezing front are recorded in top-view through a
camera connected to a long working distance lens (Thorlabs, MVL12X12Z), and
diffused back lighting.

Figure 2: Deformation of the ice front by a droplet. (a) Shape of the
solidification front as it approaches a silicone oil droplet
($5\text{\,}\mathrm{c}\mathrm{S}\mathrm{t}$) of size
$R\approx$90\text{\,}\mathrm{\SIUnitSymbolMicro m}$$ at three different
advancing velocities ($V=$0.4\text{\,}\mathrm{\SIUnitSymbolMicro
m}\text{\,}{\mathrm{s}}^{-1}$,$0.9\text{\,}\mathrm{\SIUnitSymbolMicro
m}\text{\,}{\mathrm{s}}^{-1}$,$1.6\text{\,}\mathrm{\SIUnitSymbolMicro
m}\text{\,}{\mathrm{s}}^{-1}$$, see Suppl. Movies 1-3). Although the mismatch
in thermal conductivity remains unaltered ($k_{p}/k_{l}=0.2$), a reversal of
the solidification front deflection is observed as the front approaches
faster. The scale-bars are $100\text{\,}\mathrm{\SIUnitSymbolMicro m}$. (b)
Furthest deflection $\Delta$ of the front, with respect to the planar far
field, as function of time for the three different advancing velocities
(different symbols) corresponding to (a). (c) Deflection of the front when the
drop touches the ice (right panels in (a)) as a function of the advancing
velocity $V$. When increasing the viscosity of the drop from low (blue
circles, $5\text{\,}\mathrm{c}\mathrm{S}\mathrm{t}$) to high (red squares,
$50\text{\,}\mathrm{c}\mathrm{S}\mathrm{t}$), the reversal of the
solidification front deflection at larger velocities is less pronounced.

The sequences of images in Fig. 2(a) illustrate distinct deformation behaviors
of a planar solidification front as it encounters a dispersed oil droplet
($R\sim$90\text{\,}\mathrm{\SIUnitSymbolMicro m}$$) at different velocities.
At a low advancing velocity ($V=$0.4\text{\,}\mathrm{\SIUnitSymbolMicro
m}\text{\,}{\mathrm{s}}^{-1}$$), as shown in Fig. 2(a), the initially planar
solidification interface gradually deforms towards the droplet with its
approach, creating an impression that the droplet is exerting local pull on
the solidification front. This deformation behavior, consistent with previous
studies, is expected in the presence of a particle with lower thermal
conductivity than its surroundings, due to the thermal conductivity of the oil
droplet being less than that of water. In this scenario, the heat flow chooses
the path of least resistance, thereby here avoiding the particle itself. This,
in turn, causes the otherwise uniform isotherms to deflect as illustrated in
Fig. 1(b). Correspondingly, the solidification front simply follows the
isotherm, which matches the melting temperature of water (Shangguan  _et al._
, 1992; Rempel and Worster, 2001; Park  _et al._ , 2006; Meijer  _et al._ ,
2023a). Correspondingly, for small Péclet numbers
$\mathrm{Pe}=V\,R/\kappa_{l}\ll 1$, with $\kappa_{l}$ the thermal diffusivity
of the liquid, the distortion of the isotherms surrounding a particle placed
in a uniform thermal gradient can be calculated by solving the heat equation
$\nabla^{2}T_{B}=0$ in all phases. The analytical solution of the heat
equation with appropriate heat flux boundary conditions at the particle
surface and the simplifying assumption that both the solidified bulk and the
melt have the same thermal conductivity yield Chernov  _et al._ (1977);
Shangguan  _et al._ (1992)

$T_{B}(\mathbf{x})=T_{\mathrm{ref}}+\mathbf{G}\cdot(\mathbf{x}-\mathbf{x}_{0})+k_{e}R^{3}\,\frac{\mathbf{G}\cdot(\mathbf{x}-\mathbf{x}_{0})}{|\mathbf{x}-\mathbf{x}_{0}|^{3}},$
(1)

where $k_{e}=(k_{l}-k_{p})/(2k_{l}+k_{p})$, with $-1<k_{e}<1/2$, is the
dimensionless parameter that quantifies the mismatch in the thermal
conductivities and $T_{\mathrm{ref}}$ the reference temperature at the center
of the drop, located at $\mathbf{x}_{0}$. When the particle and liquid have
different thermal conductivities, $k_{e}$ is either negative ($k_{p}/k_{l}>1$)
or positive ($k_{p}/k_{l}<1$) and the interface deflects away from or towards
the particle, respectively (see Fig. 1). Note that despite the simplified
approach of this model, it produces result that closely match with
experimental observations Shangguan  _et al._ (1992).

Since the above model description excludes any dynamic effects arising due to
heat transfer due to convection and viscous dissipation, it only follows that
for a given thermal conductivity mismatch between the particle and the
surrounding liquid, the deformation of the solidification front remains
unaffected by the rate of engulfment. While this holds for solid particles
(see Fig. 1(a)), unexpected deformation behaviors of the solidification front
are observed at higher solidification rates when encountering a droplet. We
find that as the solidification rate is increased to a higher value
$V=$0.9\text{\,}\mathrm{\SIUnitSymbolMicro m}\text{\,}{\mathrm{s}}^{-1}$$, the
solidification front remains planar as it approaches the dispersed droplet. A
subsequent further increase in the rate of solidification to
$V=$1.6\text{\,}\mathrm{\SIUnitSymbolMicro m}\text{\,}{\mathrm{s}}^{-1}$$ even
leads to the reversal of the deformation of the solidification front, giving
an impression as if the droplet is indenting the solid ice.

The physical mechanism responsible for the strange deformation behaviour of
the solidification front at faster freezing rates is intimately connected to
the interfacial dynamics along the droplet’s interface. Due to the temperature
dependence of the liquid’s surface tension, it is expected that a droplet in a
thermal gradient experiences differential surface tension along its interface.
Accordingly, in our experiments, the part of the droplet closest to the
solidification front experiences a higher surface tension than the rest of the
surface, given that $\mathrm{d}\sigma/\mathrm{d}T$ is negative for silicone
oil. This temperature-induced difference in surface tension along the droplet
interface induces an interfacial flow, drawing liquid from the warmer end to
the colder end, which is known as thermal Marangoni flow (see Fig. 3). This
fluid flow surrounding the droplet disturbs the local thermal equilibrium and
distorts the base temperature field given by Eq. (1).

Figure 3: Mechanism of the reversal of the ice front deflection. Schematic of
the thermal Marangoni advection induced by a gradient in surface tension over
the interface of the drop $\sigma(T)$, caused by the applied thermal gradient
$G$. Warmer water at the top is advected towards the colder ice front, locally
increasing the temperature field, leading to enhanced melting in front of the
droplet.

Figure 4: Isotherm perturbation due to thermal Marangoni advection. (a)
Isotherms of the perturbed temperature field around a drop with increasing
Marangoni number for $k_{p}/k_{l}=0.2$ according to our perturbation theory.
The red lines correspond to the isotherms where $T(\mathbf{x})=T_{m}$, i.e.,
the melting temperature of water. (b) Normalised self-induced velocity field
around the drop due the thermal Marangoni effect Young  _et al._ (1959). The
flow is strongest at the equator of the drop (red color). Arrows only indicate
the flow direction. (c)-(d) Numerical simulations of a solid-liquid interface
approaching an oil droplet experiencing a thermal gradient of (c)
$G=$1\text{\times}{10}^{2}\text{\,}\mathrm{K}\text{\,}{\mathrm{m}}^{-1}$$
(implying $Ma=9$) and (d)
$G=$1\text{\times}{10}^{4}\text{\,}\mathrm{K}\text{\,}{\mathrm{m}}^{-1}$$
(implying $Ma=0.07$), respectively. The results qualitatively agree with our
experimental observations (see Fig. 2(a)).

To account for this change in the temperature field due to the thermal
Marangoni flow around the droplet, we develop a simplified model. We assume
that the droplet is placed in an external thermal gradient without any
boundary effects. For such a case, the analytical description of the flow-
field generated by the thermal Marangoni effect was originally derived in
Refs. Rybczynski (1911); Young  _et al._ (1959). Following the same analysis
in our case yields a polar velocity of the fluid flow $u_{\theta}$ outside the
drop as (see Suppl. Mat. Sec. I):

$u_{\theta}(r,\theta)=-\frac{1}{2}\frac{1}{\mu+\mu^{\prime}}\left(\frac{R}{r}+\frac{R^{3}}{r^{3}}\right)\frac{\mathrm{d}\sigma}{\mathrm{d}T}\frac{k_{l}}{2k_{l}+k_{p}}GR\sin\theta,$
(2)

where $\mu$ and $\mu^{\prime}$ are the viscosity of water and oil,
respectively, and $\mathrm{d}\sigma/\mathrm{d}T<0$ quantifies the temperature
dependence of the surface tension. We evaluate Eq. (2) at the equator of the
drop and define the thermal Marangoni velocity as

$V_{M}\equiv
u_{\theta}|_{r=R,\theta=\pi/2}=-\frac{1}{\mu+\mu^{\prime}}\frac{\mathrm{d}\sigma}{\mathrm{d}T}\frac{k_{l}}{2k_{l}+k_{p}}GR.$
(3)

From Eq. (3) we see that the flow at the interface of the drop, characterised
by $V_{M}$, is proportional to the size of the drop (not varied in
experiments), the strength of the applied thermal gradient $G$, and inversely
proportional to $\mu+\mu^{\prime}$, the sum of the viscosities of water and
oil.

To assess the influence of the thermal-Marangoni flow on the local heat
transfer near the droplet, we solve the linearized energy equation. The base
state is assumed to be established by heat-transfer via purely conductive
effects and is given as $\mathbf{u}_{B}=0$ and $\nabla^{2}\mathbf{T}_{B}=0$,
where $\mathbf{T}_{B}$ is given by Eq. (1). This base state is superimposed by
the thermal field due to the thermal-Marangoni flow as
$T(\mathbf{x})=T_{B}(\mathbf{x})+T_{1}(\mathbf{x})$ and
$\mathbf{u}=\mathbf{u}_{B}+\mathbf{u}_{1}$, with
$\mathbf{u}_{1}=[u_{r},u_{\theta}]$ (see Suppl. Mat. Sec. I). Accordingly, the
linearized, non-dimensionalized energy equation then reads:

$Ma\,\tilde{\mathbf{u}}_{1}\cdot\nabla\tilde{T}_{B}=\nabla^{2}\tilde{T}_{1},$
(4)

with tildes denoting dimensionless quantities. Here, we define the Marangoni
number $Ma\equiv V_{M}R/\kappa$ as the ratio of advective over diffusive
thermal transport and it is the same as the Péclet number introduced earlier.
Typically, here we have $Ma\sim\mathcal{O}(1)$. All lengths are non-
dimensionalized by the droplet radius $R$, velocities by $V_{M}$, and
temperatures as
$\tilde{T}=T-T_{\mathrm{cold}}/(T_{\mathrm{m}}-T_{\mathrm{cold}})$, where
$T_{\mathrm{cold}}$ and $T_{\mathrm{m}}$ are temperatures at the cold end of
the experimental cell and the melting point of water, respectively. We solve
Eq. (4) numerically to determine $T_{1}(\mathbf{x})$ and hence $T(\mathbf{x})$
(see Fig. 4(a)). The Marangoni number $Ma$ is the natural control parameter
dictating the strength of the velocity field around the drop (see Fig. 4(b))
and hence the extent of changes in the local thermal field. Considering that
the freezing interface follows the isotherm matching the melting temperature
of water (red lines in Fig. 4(a)), we qualitatively capture the same
deformation behaviour as seen in the experiments (see Fig. 2(a)). In line with
Eq. (3), by varying the temperature gradient $G$, and hence the advancing
velocity of the front, or by changing the viscosity of the oil (see Fig.
2(c)), we vary the strength of the flow around the drop in our experiments,
and thus the extend of the ice deformation.

The simplified approach discussed earlier offers a qualitative depiction of
the unique deformation behaviors of the solidification front and its
relationship with the applied thermal gradient. However, its precision in
faithfully simulating the influence of the solidification front on the flow
field around the droplet, and thus the overall deformation, is limited.
Consequently, here we have also introduced a detailed numerical approach that
systematically captures all interactions (refer to Supplementary Material,
Section II). The results obtained from finite element simulations (Fig 4(c,d))
quantitatively agree with the experimental observations.

In summary, our integrated approach involving experiments, theoretical
analysis, and numerical simulations disentangles various effects occurring
when dispersed droplets in a molten medium interact with an advancing freezing
front. Our findings reveal that, contingent on the freezing rate, a droplet
can induce unexpected deformations in the freezing interface, in both
directions. This phenomenon is attributed to the interfacial flow generated by
the thermal Marangoni effect, leading to alterations in the temperature field
and, consequently, the shape of the solid-liquid interface in its proximity.
Our findings contribute to a better understanding of the mechanism that govern
particle rejection or engulfment during the solidification of a complex
medium, which is crucial to step towards designing advanced material with
functional properties Deville  _et al._ (2009).

## Acknowledgements

The authors thank Gert-Wim Bruggert and Martin Bos for the technical support.
The authors acknowledge the funding by Max Planck Center Twente and the Balzan
Foundation.

## References

  * Shangguan  _et al._ (1992) D. Shangguan, S. Ahuja, and D. Stefanescu, An analytical model for the interaction between an insoluble particle and an advancing solid/liquid interface, Metallurgical Transactions A 23, 669 (1992).
  * Rempel and Worster (2001) A. Rempel and M. Worster, Particle trapping at an advancing solidification front with interfacial-curvature effects, Journal of Crystal Growth 223, 420 (2001).
  * Park  _et al._ (2006) M. S. Park, A. A. Golovin, and S. H. Davis, The encapsulation of particles and bubbles by an advancing solidification front, Journal of Fluid Mechanics 560, 415 (2006).
  * Dedovets  _et al._ (2018) D. Dedovets, C. Monteux, and S. Deville, Five-dimensional imaging of freezing emulsions with solute effects, Science 360, 303 (2018).
  * Tyagi  _et al._ (2020) S. Tyagi, H. Huynh, C. Monteux, and S. Deville, Objects interacting with solidification fronts: Thermal and solute effects, Materialia 12, 100802 (2020).
  * Meijer  _et al._ (2023a) J. G. Meijer, V. Bertin, and D. Lohse, Frozen cheerios effect: Particle-particle interaction induced by an advancing solidification front, arXiv preprint arXiv:2311.09477 (2023a).
  * Wettlaufer and Worster (2006) J. Wettlaufer and M. G. Worster, Premelting dynamics, Annu. Rev. Fluid Mech. 38, 427 (2006).
  * Meijer  _et al._ (2023b) J. G. Meijer, P. Kant, D. Van Buuren, and D. Lohse, Thin-film-mediated deformation of droplet during cryopreservation, Physical Review Letters 130, 214002 (2023b).
  * Meijer  _et al._ (2023c) J. G. Meijer, P. Kant, and D. Lohse, Freezing-induced topological transition of double-emulsion, arXiv preprint arXiv:2312.04317 (2023c).
  * Van Elburg  _et al._ (2021) B. Van Elburg, G. Collado-Lara, G.-W. Bruggert, T. Segers, M. Versluis, and G. Lajoinie, Feedback-controlled microbubble generator producing one million monodisperse bubbles per second, Review of Scientific Instruments 92 (2021).
  * Kao  _et al._ (2009) J. C. Kao, A. A. Golovin, and S. H. Davis, Particle capture in binary solidification, Journal of Fluid Mechanics 625, 299 (2009).
  * Tyagi  _et al._ (2022) S. Tyagi, C. Monteux, and S. Deville, Solute effects on the dynamics and deformation of emulsion droplets during freezing, Soft Matter 18, 4178 (2022).
  * Chernov  _et al._ (1977) A. Chernov, D. Temkin, and A. Mel’nikova, The influence of the ther- mal conductivity of a macroparticle on its capture by a crystal growing from a melt, Sov. Phys. Crystallogry 22, 656 (1977).
  * Young  _et al._ (1959) N. Young, J. S. Goldstein, and M. Block, The motion of bubbles in a vertical temperature gradient, Journal of Fluid Mechanics 6, 350 (1959).
  * Rybczynski (1911) W. Rybczynski, Uber die fortschreitende bewegung einer flussigen kugel in einem zahen medium, Bull. Acad. Sci. Cracovie A 1, 40 (1911).
  * Deville  _et al._ (2009) S. Deville, E. Maire, G. Bernard-Granger, A. Lasalle, A. Bogner, C. Gauthier, J. Leloup, and C. Guizard, Metastable and unstable cellular solidification of colloidal suspensions, Nature Materials 8, 966 (2009).