# Feasible and Desirable Counterfactual Generation by Preserving Human Defined
Constraints

Homayun Afrabandpey
Nokia Technologies
Finland
<EMAIL_ADDRESS>
&Michael Spranger
Sony AI
Japan
<EMAIL_ADDRESS>
###### Abstract

We present a human-in-the-loop approach to generate counterfactual (CF)
explanations that preserve global and local feasibility constraints. Global
feasibility constraints refer to the causal constraints that are necessary for
generating actionable CF explanation. Assuming a domain expert with knowledge
on unary and binary causal constraints, our approach efficiently employs this
knowledge to generate CF explanation by rejecting gradient steps that violate
these constraints. Local feasibility constraints encode end-user’s constraints
for generating desirable CF explanation. We extract these constraints from the
end-user of the model and exploit them during CF generation via user-defined
distance metric. Through user studies, we demonstrate that incorporating
causal constraints during CF generation results in significantly better
explanations in terms of feasibility and desirability for participants.
Adopting local and global feasibility constraints simultaneously, although
improves user satisfaction, does not significantly improve desirability of the
participants compared to only incorporating global constraints.

##  1 Introduction

Complex Machine Learning (ML) models have been adopted in many real-world
decision-making tasks either to support humans or even substitute them.
Despite their superior performance, the black-box nature of these models
necessitates need for interpretability methods to explain their automated
decisions for individuals who are subject to these decisions. Among the large
body of literature on interpretable ML [2, 10], Counterfactual (CF)
explanations have shown promise for practitioners. A CF explanation contains
one or more CF instances. A CF instance is a perturbed version of the original
instance that flips the black–box model’s prediction. By comparing a CF
instance with the original instance, a human user receives hints on what
changes to the current situation would have resulted in an alternative
decision, i.e., “If $\bm{X}$ was $\bm{X}^{\prime}$, the outcome would have
been $y^{\prime}$ rather than $y$.”

Generating CF explanations that are useful in real-world is still challenging.
First, CF explanations generated by many existing works do not take into
account causal relationships among features. This results in CF instances that
are not actionable in the real-world. Take a loan application as example, a CF
explanation approach that does not adopt causal relationships among features
could suggests to change the present employment type from “newbie”111In German
credit dataset, this qualitative value is defined as a person who has less
than $1$ year experience in his/her current job. to “senior”222Between $4$ to
$7$ years of experience according to German credit dataset. while the age is
unchanged. Second, CF explanations are subjective and should be personalized,
while existing works do not take into account constraints from the end-users
of the ML models. Back to the loan application example, a user might find it
feasible to change the housing type, while another might instead prefer
changing, e.g., number of installments (duration in month).

To bridge these gaps, we propose a new approach to generate CF explanations
for any differentiable classifier via feasible perturbations. For this, we
extend [8] by formulating an objective function for generating CF instances
that takes into account two types of feasibility constraints:

  * •

Global feasibilities: unary and binary monotonic causal constraints extracted
from a domain expert,

  * •

Local feasibilities: constraints in the form of feature perturbation
difficulty values, given by the end-users.

The objective function is optimized using gradient descent and feasibility
constraints are satisfied during the optimization by rejecting gradient steps
that do not satisfy them. It is important to note that, here we differentiate
between end-user and domain expert. An end-user is the individual who is
subject to the decision of the ML model, e.g. a bank customer whose loan
application is rejected. A domain expert, on the other hand, knows the data
and the application. We believe domain experts are naturally able to give
feedback on causal relationship among (at least) some features, without being
constrained to know the exact functional relationship.

The same feasibility constraints were also considered in [7] for CF
generation. They propose a generative model based on an encoder-decoder
framework, where the encoder projects features into a latent space and the
decoder generates CF instances from the latent space. Their approach, however,
requires complete information about the structural causal model including the
causal graph and the structural equations. This assumption is highly
restrictive for applicability of the method in real-world applications. To
cope with this issue, [7] proposed a data driven approach to approximate unary
and binary monotonic causal constraints and adopt the approximated
relationships in the CF generation. For local feasibility constraints, they
considered implicit user preferences, i.e., given a pair of original instance
and CF instance, $\left(\bm{x},\bm{x}^{\prime}\right)$, the user outputs $1$
if CF instance is locally feasible and $0$ otherwise. However, since there is
no access to the $\left(\bm{x},\bm{x}^{\prime}\right)$ query pairs apriori,
they approximate the user by first asking user preferences on some
$\left(\bm{x},\bm{q}\right)$, where $\bm{q}$ are sample CF instances generated
by a CF generator without considering user preferences, and then learn a model
that generates scores for each pair that mimics user preferences.

Our approach is different from [7] in several aspects:

  * •

in [7], for approximating each binary constraint, the model learns $2$ extra
parameters. This hinders the scalability of the method. Furthermore, these
approximated binary constraints could be imprecise as they are learned from
the data, while in our approach we rely on domain experts to provide such
constraints which is more reliable,

  * •

local feasibility constraints are incorporated via implicit feedbacks that are
approximated using a function. These feedbacks are not directly related to the
final CF instances to be generated. This could result in undesirable CF
instances that do not satisfy user’s constraints. On the other hand, we adopt
explicit user feedbacks directly into the optimization function,

  * •

the type of the user feedback considered in [7] for local feasibility is
difficult to provide and restrictive. It is difficult to provide since the
user must compare the CF instance with the original instance to find out if
perturbations are locally feasible or not. It is restrictive because the
approach provides no tool for the user to state the level of local
infeasibilty. As an example, assume a CF instance is generated by perturbing
more than one feature of the original instance where all but one perturbation
satisfy user’s feasibility constraints. In our approach, user feedbacks are
"feature level" and they are not restricted to $\\{0,1\\}$,

  * •

last but not least, [7] did not test they approach in a real user study and it
is not evident from the paper how a real user could be adopted in-the-loop to
obtain desirable CF explanation.

To explore the effectiveness of our explanations, we design user studies where
users are asked to rank CF instances generated under different conditions.
Through these studies, we found that users tend to give significantly better
ranks to CF instances generated by considering global feasibility constraints
compared to the case where such constraints are not considered. Furthermore,
CF instances generated by adopting both local and global feasibility
constraints are better than those generated by only considering global
feasibility constraints. However, their difference is not statistically
significant.

In summary, we make the following contributions:

  * •

we propose a novel method to generate CF explanations that preserve both local
and global feasibility constraints extracted from end-users and domain
experts, respectively. This is obtained via an optimization task rather than
relying on heuristics [4, 11],

  * •

we design and conduct user studies to demonstrate the quality of generated CF
instances. Our studies confirm that CF instances generated by capturing causal
relationships are more favorable for end-users compared to those generated
without causal constraints. Adding local feasibility constraints to the CF
generation can further improve user satisfaction from CF instances.

##  2 Counterfactual Generation With Local Feasibility via User-Defined
Metrics

Assume we want to explain the undesirable prediction $y$, of a binary black-
box classifier $f$, for an instance $\bm{x}\in\mathbb{R}^{p}$. Throughout the
paper, we assume that the black–box model is differentiable and does not
change over time. In their seminal paper, Mothilal et. al [8] proposed the
following optimization function to generate a set of $K$ diverse CF instances
to explain the prediction $f\left(\bm{x}\right)$:

$\arg\min_{\bm{c}_{1},\ldots,\bm{c}_{K}}\frac{1}{K}\sum_{k=1}^{K}\mathcal{L}\left(f(\bm{c}_{k}),y^{\prime}\right)+\frac{\lambda_{1}}{K}\sum_{k=1}^{K}d\left(\bm{c}_{k},\bm{x}\right)-\lambda_{2}\mbox{d\\_div}\left(\bm{c}_{1},\ldots,\bm{c}_{K}\right),$
(1)

where $\mathcal{L}$ is the loss function that pushes the predictions of $f$
for CF instances toward the desirable prediction $y^{\prime}$,
$d\left(.,.\right)$ is a distance measure to keep the CF instances close to
$\bm{x}$, $\mbox{d\\_div}\left(.,.\right)$ is the diversity metric among CF
instances and $\lambda_{1}$ and $\lambda_{2}$ are regularizers that control
the relative importance of diversity among the CF instances and proximity of
CF instances to $\bm{x}$.

Of particular importance is the choice of distance measure $d\left(.,.\right)$
and diversity metric $\mbox{d\\_div}\left(.\right)$. These metrics are
subjective. When generating CF instances, a user might prefer to keep the
distance between some features of the CF instance and the original instance
zero, i.e., user does not tolerate any change to values of these features,
while other features are adjustable. The same argument is valid for the
diversity measure. To take into account these subjective constraints, we
assume each user provides her preferences for perturbing features during CF
generation via feature perturbation difficulty values. These values are
exploited in the CF generation process via a similarity metric to be used when
computing proximities.

###  2.1 Proximity

The closeness of a set of CF instances to their original instance is equal to
the mean of their negative distances. Following [12, 8], we define separate
distance metrics for continuous and categorical features.

Continuous Features. The distance between each continuous feature of a CF
instance to the corresponding feature in the original instance can be
formulated as the Mahalanobis distance with a user defined metric $\bm{A}$:

$d_{\bm{A}}^{\mbox{cont}}\left(\bm{c}_{k},\bm{x}\right)=\left(\bm{c}_{k}-\bm{x}\right)^{T}\bm{A}_{\bm{\gamma}}\left(\bm{c}_{k}-\bm{x}\right)=\sum_{j=1}^{p_{cont}}\gamma_{j}\left(\bm{c}_{kj}-\bm{x}_{j}\right)^{2}$
(2)

where $\bm{A}_{\bm{\gamma}}\in\mathbb{R}^{p_{cont}\times p_{cont}}$ is a
diagonal matrix with diagonal elements being encoded in the vector
$\bm{\gamma}$, $\gamma_{j}$ is the perturbation difficulty value of the $j$th
feature, and $p_{cont}$ is the total number of continuous features. Users are
restricted to assign only positive values to $\gamma_{j}$ to ensure positive
semi-definitiy of the Mahalanobis distance metric $\bm{A}_{\bm{\gamma}}$. Each
$\gamma_{j}$ indicates how much influential the feature is in determining the
distance between the CF instance and the original instance; a large value for
a feature implies that the feature is very difficult for the user to perturb
while decreasing the value toward $1$ means that it is increasingly easier to
change [1]. To cancel out the effect of different feature ranges, feature-wise
distances can be divided by the standard deviation of the feature.

It has been shown in [12] that Manhattan distance normalized by the median
absolute deviation of features has desirable properties. This metric can also
adopt $\\{\gamma_{j}\\}_{j=1}^{p_{cont}}$ as multipliers for feature-wise
distances.

Categorical Features. For categorical features, we use the overlap distance
with feature perturbation difficulty values as:

$d_{\bm{A}}^{\mbox{cat}}\left(\bm{c}_{k},\bm{x}\right)=\sum_{j}^{p_{cat}}\gamma_{j}I\left(\bm{c}_{kj}\neq\bm{x}_{j}\right),$
(3)

where $I(.)$ is the index function that returns $1$ if the condition inside it
is true and 0 otherwise, and $p_{cat}$ is the total number of categorical
features. Based on this metric, when $\gamma$ is large, a miss–match of the
values of two categorical features results in high cost and vice versa.

###  2.2 Diversity

Following [8], we use the determinant of the kernel matrix of CF instances:

$\mbox{d\\_div}=\left|\mathbf{K}\right|,$ (4)

where
$\mathbf{K}_{ij}=\frac{1}{1+d\left(\mathbf{c}_{i},\mathbf{c}_{j}\right)}$, and
$d\left(\mathbf{c}_{i},\mathbf{c}_{j}\right)$ is the distance between two
counterfactual instances $i$ and $j$ as defined in the previous subsection.

###  2.3 Loss Function

We use hinge loss defined as:

$\max\left(0,1-z\times\mbox{logit}\left(f\left(\bm{c}\right)\right)\right),$
(5)

where $z=-1$ for $y=0$ and $z=1$ for $y=1$ and $f(\bm{c})$ is the output of
the model before entering the softmax layer (unscaled output).

##  3 Modeling Global Feasibility Constraints

We consider the case where the structural causal model of the observed data is
unknown, however there is a domain expert who is able to provide unary and
binary causal constraints at least for some features. Unary constraints
stipulate whether a feature can increase or decrease. For example, a unary
constraint on “Age” restricts it to only increase since in real-world it is
not possible to decrease the age. Binary constraints capture the causal
relationship between two features. An example of a binary causal relationship
between “education” and “age” states that increasing “education” results in
increase of “age”. We focus on monotonic binary constraints where increasing
(decreasing) an upstream feature, causes an increase (decrease) in its child.

Given unary and binary causal constraints by the domain expert, our goal is to
exploit them when generating CF instances. To model binary constraints, we
define a vector $\bm{b}^{i}\in\mathbb{R}^{p}$ for each feature
$i\in\\{1,2,3\cdots,p\\}$ that represents causal relationships where $x_{i}$
is the upstream variable,

$\bm{b}^{i}_{j}=\begin{cases}-1&\text{if }\left(x_{i}\rightarrow
x_{j}\right)\mbox{ OR }\left(i==j\right)\\\ \;0&\text{otherwise}\end{cases},$
(6)

where $x_{i}\rightarrow x_{j}$ refers to the binary causal constraint where
any change (increase/decrease) in the value of the $i$th feature monotonically
changes the value of the $j$th feature.

Let $\mathbb{L}$ be the function to be minimized (using gradient descent) to
generate CF instances. At each iteration of the optimization,
$\nabla_{\bm{c}}\mathbb{L}=\left[\frac{\partial\mathbb{L}}{\partial\bm{c}_{1}},\frac{\partial\mathbb{L}}{\partial\bm{c}_{2}},\cdots,\frac{\partial\mathbb{L}}{\partial\bm{c}_{p}}\right]^{T}$
where the $i$th element of the gradient vector,
$\frac{\partial\mathbb{L}}{\partial\bm{c}_{i}}$, determines the direction
(sign of $\frac{\partial\mathbb{L}}{\partial\bm{c}_{i}}$) and the size
(magnitude of $\frac{\partial\mathbb{L}}{\partial\bm{c}_{i}}$) of the step to
be taken w.r.t. the $i$th feature to reach the minimum of $\mathbb{L}$. The
combination of $\nabla_{\bm{c}}\mathbb{L}$ and $\bm{b}$ tells whether or not
the gradient update violates existing binary causal constraints. This is
obtained by

$\bm{v}^{i}=\left(\bm{b}^{i}\circ\mbox{sgn}\left(\nabla_{\bm{c}}\mathbb{L}\right)\right)+\left(\mbox{sgn}\left(\frac{\partial\mathbb{L}}{\partial\bm{c}_{i}}\right).\mathbf{1}\right),\;\;\;\forall
i\in\\{1,\cdots,p\\}$ (7)

where $\circ$ is the element-wise multiplication, $\mbox{sgn}\left(.\right)$
is the sign function, and $\mathbf{1}\in\mathbb{R}^{p}$ is the vector of all
ones. Elements of $\bm{v}^{i}\in\mathbb{R}^{p}$ could have values in $\\{0,\pm
1,\pm 2\\}$ with the following interpretation:

  * •

when $v^{i}_{j}=0$, its either because there is no (monotonic) causal
relationship between the two features, or because the gradient signs satisfies
the causal relationship, or because $i=j$,

  * •

$v^{i}_{j}=\pm 2$ happens when there is a causal relationship between the two
features, but the gradient signs for the two features violates the causal
relationship,

  * •

finally, $v^{i}_{j}=\pm 1$ happens only when the gradient of the loss function
w.r.t one of the features is zero. If the gradient w.r.t. the downstream
feature, i.e., $x_{j}$, is zero, then the causal relationship is violated
since the upstream feature changes while the downstream feature remains
constant. Otherwise, the update is valid.

With the above interpretation and given the vectors
$\\{\bm{v}^{i}\\}_{i=1}^{p}$, one accepts or rejects gradient update w.r.t.
the $i$th feature as follows:

$\frac{\partial\mathbb{L}}{\partial\bm{c}_{i}}=\begin{cases}0&\text{if
}\left(v^{i}_{j}=\pm 2\right)\vee\left(v^{i}_{j}=\pm
1\wedge\frac{\partial\mathbb{L}}{\partial\bm{c}_{j}}=0\right)\\\
\frac{\partial\mathbb{L}}{\partial\bm{c}_{i}}&\text{otherwise}\end{cases},\;\;\;\forall
j\in\\{1,\cdots,p\\}.$ (8)

To apply unary constraints, we define two hypothetical features shown by
$U^{+}$ and $U^{-}$. These can only be used as upstream features;
$U^{+}\rightarrow x_{k}$ determines that $x_{k}$ can only increase, while
$U^{-}\rightarrow x_{k}$ states that the feature can only decrease. We define
vectors $\bm{b}^{+}\in\mathbb{R}^{p}$ and $\bm{b}^{-}\in\mathbb{R}^{p}$ to
model unary constraints where each row of these vectors determines whether or
not the feature of the corresponding row is downstream feature for $U^{+}$ or
$U^{-}$, respectively. With this definition, we have
$\bm{v}^{+}=\bm{b}^{+}\circ\mbox{sgn}\left(\nabla_{\bm{c}}\mathbb{L}\right)$
and
$\bm{v}^{-}=\bm{b}^{-}\circ\mbox{sgn}\left(\nabla_{\bm{c}}\mathbb{L}\right)$.
Both vector $\bm{v}^{+}$ and $\bm{v}^{-}$ have values in $\\{0,\pm 1\\}$. When
$\bm{v}^{+}_{j}=-1$, the unary constraint is violated by the update to the
$j$th feature, while $\bm{v}^{+}_{j}=0$ means that the update is zero and
$\bm{v}^{+}_{j}=+1$ represents that the gradient update satisfies the
constraint. For $\bm{v}^{-}$, if the value of an element equals $+1$, the
constraint is violated.

##  4 Empirical Evaluation

We conduct two user studies as explained in the following subsections to
verify how good are the explanations generated using our approach when
compared to DiCE [8].

###  4.1 Datasets and Model

The datasets used in the experiments are Adult and german credit datasets. We
adopted the pre-processed version of the adult dataset based on [3] with $8$
features, namely, age, work class, education level, marital status,
occupation, race, gender, and hours per week. For german credti, we consider
$9$ demographic and socio-economic features, including: duration in month,
credit history, credit amount, present employment since unemployed, sex, age,
job, and number of people liable to provide maintenance for. For both
datasets, categorical features are one-hot-encoded and continuous features are
scaled between $0$ and $1$. Datasets are divided into $80\%-20\%$ train and
test sets.

The black-box model is a single layer neural network trained for $20$ epochs
with learning rate $0.01$ using ADAM optimizer [5] to minimize cross entropy
loss. Accuracies of the trained neural network on held-out test set are $83\%$
and $73\%$ for Adult and German Credit, respectively. In all experiments, we
followed [8] and used Manhattan distance normalized by median absolute
deviation of features as distance metric for continuous features when
generating counterfactual instances using DiCE and our approach. Both
regularizers in Eq. 1, i.e., $\lambda_{1}$ and $\lambda_{2}$, are set to $1$
in the implementation. All experiments were run on a $1.9$ GHz CPU with $8$ GB
RAM.

###  4.2 Baselines

We focus on DiCE [8] as baseline, because our implementation is based on the
DiCE code and therefore comparison is easier. In the following subsection, for
the sake of brevity, we call our proposed method C-DiCE, i.e. DiCE that also
takes into account Causal constraints.

###  4.3 Feasibility Constraints and User Satisfaction

Most of the existing approaches for counterfactual explanation are validated
without user experiment. A key limitation of this approach is that predefined
metrics such as those adopted in [7, 8, 9], do not precisely capture the human
cognition when evaluating subjective criteria such as desirability of CF
instances. For this, we develop experiments with real users to study the
goodness of CF explanations considering following conditions,

  * •

C1: CF instances generated by DiCE, i.e., without feasibility constraints;

  * •

C2: CF instances generated by C-DiCE by taking into account only global
feasibility constraints;

  * •

C3: CF instances generated by C-DiCE considering both global and local
feasibility constraints.

Users are asked to rank CF instances generated with DiCE, C-DiCE:C2, and
C-DiCE:C3, for several real instances based on three criteria: (i) validity,
(ii) feasibility, and (iii) desirability. Validity determines whether or not
the generated CF instance flips the outcome of the original instance.
Feasibility states whether or not the values assigned to each feature or
combination of features are feasible considering real-world constraints.
Finally, desirability reflects users’ satisfacation on CF instances. Our
hypotheses were:

  * •

H1: CF instances generated by C-DiCE:C2 have, on average better ranks in terms
of the above-mentioned criteria, than those generated by DiCE;

  * •

H2: C-DiCE:C3 generates CF instances with better average ranks compared to
both DiCE and C-DiCE:C2 in terms of the introduced criteria.

We conduct two user studies: in the first one users were asked to rank CF
instances generated by DiCE and C-DiCE:C2 using Adult dataset, and in the
second study users were asked to rank CF instances generated by C-DiCE:C2 and
C-DiCE:C3 using German credit dataset.

Since we do not have the true causal constraints, we asked $4$ data scientists
who were familiar with the datasets to determine unary and binary causal
constraints. Extracted constraints highly agreed with agreement rate above
$90\%$ on average over all experts, however we only adopted constraints all
experts agreed on. Constraints defined for adult dataset are:
$education\rightarrow age$, $U^{+}\rightarrow age$, and $U^{+}\rightarrow
education$. For German credit dataset the feasibility constraints are:
$present\\_employment\\_since\\_unemployed\rightarrow age$, $job\rightarrow
age$, $U^{+}\rightarrow age$.

From the test set of each dataset, we randomly assign $5$ samples with
undesirable output to each participant. These samples are drawn randomly
without replacement to ensure that participants are assigned unique samples.
For each sample, we generate $10$ CF instances; in the first experiment, $5$
using DiCE and $5$ using C-DiCE:C2, and in the second experiment, $5$ using
C-DiCE:C2 and $5$ using C-DiCE:C3. These CF instances are shuffled and showed
to a participant together with their original instance, however, the
participant is unaware of the condition using which the CF instances are
generated. CF instances are visualized to the user in a spreadsheet.
Developing a more user-friendly tool for visualization and feedback collection
is out of the scope of this paper. The participant is asked to rank these CF
instances from $1$ to $10$, where $1$ is the best rank and $10$ is the worst.
Participants are not allowed to assign same ranks to several CF instances
unless they are identical. At the end of the user study, we conduct an
interview with the participants to discuss their ranks and ensure they
understood the task correctly.

---
(a)

---
(b)

---
(c)

---
(d)

Figure 1: Results of the user studies for comparing CF instances generated by
(a,b) DiCE and C-DiCE:C2 and (c,d) C-DiCE:C2 and C-DiCE:C3: (a,c) mean of the
ranks provided by each user, (b,d) mean and standard deviation of the ranks
provided by all users. In (a) and (b) results of User 3 are discarded.

We adopted a between subject study: $18$ participants for first experiment
($10$ male)333We hired $19$ participants for this study, but one of them (user
$3$) was discarded as s/he provided incorrect input. This was pointed out
during the interview., and $16$ participants for the second one ($9$ male).
Participants were from different backgrounds including computer science
($13$), engineering fields ($8$), mathematics ($5$), social sciences ($5$),
and management science ($3$). They were also from different education level
including Bachelor ($6$), Masters ($8$), PhD student in 2nd or 3rd year
($10$), and PhD graduated ($10$); were aged $20$-$45$. Each participant was
first introduced to the concept and the task using a toy example and then
asked to complete the test task. To generate CF instances using C-DiCE:C3, we
first ask users to provide their local feasibility constraints by assigning
feature perturbation values $\gamma_{j}$ to each feature. In our
implementation, we use $\gamma_{j}\in\left[1,5\right]$. Each experiment took
$\sim 1$ hour and a movie ticket was awarded.

Figure 1.(a,b) demonstrates results of the first user study with Adult
dataset. Figure 1.a compares the average ranks provided by each participant to
each method, DiCE and C-DiCE:C2. All users provide better ranks to CF
instances generated by C-DiCE:C2. Figure 1.b compares mean and standard
deviation of the ranks over all users. Note that average rank over all users
is a value in the range $\left[3,8\right]$.

To better compare the superiority of the CF instances generated by C-DiCE:C2
compared to DiCE, Table 1 demonstrates the ratio of the top-$k$ CF instances
generated by DiCE and C-DiCE:C2. According to the table, CF instances
generated by C-DiCE:C2 were the top-$1$ CF instance almost twice as those
generated by DiCE ($64.4\%$ compared to $35.6\%$). The difference is larger
for top-$2$ and top-$3$ CF instances.

Table 1: Ratio (%) of top-$k$ CF instances generated by each method to all CF instances for $k=1$, $2$, $3$ in Adult dataset. Values are averaged over $18$ users. | Top-1 | Top-2 | Top-3
---|---|---|---
C-DiCE:C2 | $\bm{64.4}$ | $\bm{70.6}$ | $\bm{74.2}$
DiCE | $35.6$ | $29.4$ | $25.8$

Figure 1.(c,d) shows the results of the second experiment. Figure 1.c
demonstrates that $10$ out of $16$ users gave better ranks to the CF instances
generated by C-DiCE:C3, where the differences between the average ranks
provided by User 27 and User 35 are very marginal. Figure 1.d shows that
average ranks over all users for the two approaches are very close to each
other.

Similar to the first experiment, Table 2 demonstrate the ratio of the top-$k$
ranked CF instances generated by C-DiCE:C3 compared to those generated by
C-DiCE:C2. The table demonstrates that users prefer CF instances generated by
C-DiCE:C3, however the differences in the top-$k$ values are not as large as
those shown in Table 1. In this experiment, there were many identical CF
instances generated by C-DiCE:C2 and C-DiCE:C3. Participants assigned same
ranks to identical CF instances; consequently summation of the values of
corresponding columns in Table 2 are not $1$.

Table 2: Ratio (%) of top-$k$ CF instances generated by C-DiCE under two different conditions to all CF instances for $k=1$, $2$, $3$ in German credit dataset. Values are averaged over $16$ users. | Top-1 | Top-2 | Top-3
---|---|---|---
C-DiCE:C2 | $55$ | $50.7$ | $53.5$
C-DiCE:C3 | $\bm{68.7}$ | $\bm{59.3}$ | $\bm{58}$

We use Bayesian T-test [6] to assess the significance of the differences
between average ranks provided by the participants to CF instances generated
by different approaches. Bayesian T-test constructs a distribution for the
mean and standard deviation for the group of ranks given by the users to the
CF instances generated by each approach. Then, it constructs a probability
distribution over the differences between the group-specific distributions
using MCMC estimation. Figure 2 demonstrate these distributions for the (a)
first experiment, where $\mu_{1}$ and $\mu_{2}$ refer to the average ranks of
the CF instances generated by C-DiCE:C2 and DiCE, respectively, and (b) second
experiment, where $\mu_{1}$ and $\mu_{2}$ refer to the average ranks given by
the participants to CF instances generated by C-DiCE:C3 and C-DiCE:C2,
respectively. Each distribution include the mean credible value as the best
guess of the actual difference and the $95\%$ Highest Density Interval (HDI)
as the range were the actual difference is with $95\%$ credibility. If the
$95\%$ HDI includes zero, the difference between the average ranks is not
significant. Figure 2.a, Figure 1, and Table 1 confirm that CF instances
generated by C-DiCE:C2 are significantly better than those generated by DiCE.
This difference is due to the adoption of global feasibility constraints. On
the other hand, in Figure 2.b, the $95\%$ HDI contains zero which states that
the difference between average ranks of CF instances generated by C-DiCE:C2
and C-DiCE:C3 is not statistically significant.

---
(a)

---
(b)

Figure 2: Results of the Bayesian t-test for the (a) first experiment where
$\mu_{1}$ and $\mu_{2}$ refer to the average ranks given to the CF instances
generated by C-DiCE:C2 and DiCE, respectively, and (b) second experiment where
$\mu_{1}$ and $\mu_{2}$ demonstrate average ranks of the CF instances
generated by C-DiCE:C3 and C-DiCE:C2, respectively. If the $95\%$ HDI contains
zero, then the difference between the average ranks of the CF instances
generated by the two approach is not statistically significant.

##  5 Conclusion and Discussion

Building upon prior work on CF explanations, we introduce a novel approach to
generate feasible and desirable CF explanation. We consider two levels of
feasibility, i.e., global and local. The former is defined in terms of causal
constraints among variables and is extracted from domain experts, while the
later captures end-user defined constraints. For global feasibility, the
current work accounts for unary and binary monotonic causal constraints as two
most common types of constraints and left the more complicated constraints for
future work. User studies demonstrate the effectiveness of the proposed
approach in increasing user satisfaction about CF explanation. Designing a
tool to help domain experts to provide global feasibility constraints
especially in datasets with large number of features is considered as a future
work.

## References

  * Afrabandpey et al. [2019] H. Afrabandpey, T. Peltola, and S. Kaski.  Human-in-the-loop active covariance learning for improving prediction in small data sets.  In _Proceedings of the 28th International Joint Conference on Artificial Intelligence_ , pages 1959–1966, 2019.
  * Burkart and Huber [2021] N. Burkart and M. F. Huber.  A survey on the explainability of supervised machine learning.  _Journal of Artificial Intelligence Research_ , 70:245–317, 2021.
  * Haojun [2016] Z. Haojun.  Predicting earning potential using the adult dataset.  _https://rpubs.com/H_Zhu/235617_ , 2016.
  * Karimi et al. [2020] A.-H. Karimi, G. Barthe, B. Balle, and I. Valera.  Model-agnostic counterfactual explanations for consequential decisions.  In _Proceedings of the 23rd International Conference on Artificial Intelligence and Statistics (AISTATS)_ , pages 895–905, 2020.
  * Kingma and Ba [2015] D. P. Kingma and J. L. Ba.  Adam: A method for stochastic gradient descent.  In _Proceedings of the 3rd International Conference on Learning Representations (ICLR)_ , 2015.
  * Kruschke [2013] J. K. Kruschke.  Bayesian estimation supersedes the t test.  _Journal of Experimental Psychology: General_ , 142(2):573–603, 2013.
  * Mahajan et al. [2019] D. Mahajan, C. Tan, and A. Sharma.  Preserving causal constraints in counterfactual explanations for machine learning classifiers.  _NeurIPS Workshop on Do the right thing: Machine learning and Causal Inference for improved decision making_ , 2019.
  * Mothilal et al. [2020] R. K. Mothilal, A. Sharma, and C. Tan.  Explaining machine learning classifiers through diverse counterfactual explanations.  In _Proceedings of the ACM Conference on Fairness, Accountability, and Transparency_ , pages 607–617, 2020.
  * Sharma et al. [2020] S. Sharma, J. Henderson, and J. Ghosh.  CERTIFAI: A common framework to provide explanations and analyse the fairness and robustness of black-box models.  In _Proceedings of the AAAI/ACM Conference on AI, Ethics, and Society_ , page 166–172, 2020.
  * Tjoa and Guan [2020] E. Tjoa and C. Guan.  A survey on explainable artificial intelligence (xai): Toward medical xai.  _IEEE transactions on neural networks and learning systems_ , 32(11):4793–4813, 2020.
  * Tolomei et al. [2017] G. Tolomei, F. Silvestri, A. Haines, and M. Lalmas.  Interpretable predictions of tree-based ensembles via actionable feature tweaking.  In _Proceedings of the 23rd ACM International Conference on Knowledge Discovery and Data Mining (SIGKDD)_ , pages 465–474, 2017.
  * Wachter et al. [2017] S. Wachter, B. Mittelstadt, and C. Russell.  Counterfactual explanations without opening the black box: Automated decisions and the gpdr.  _Harv. JL & Tech._, 31:841, 2017.


# Hi Imaging of a Blueberry Galaxy Suggests a Merger Origin

Saili Dutta1, Apurba Bera2,3, Omkar Bait4, Chaitra A. Narayan1, Biny
Sebastian5, Sravani Vaddi6
1 National Centre for Radio Astrophysics, Tata Institute of Fundamental
Research, Pune University Campus, Pune 411007, India
2 International Centre for Radio Astronomy Research, Curtin University,
Bentley, WA 6102, Australia
3 Inter-University Centre for Astronomy and Astrophysics, Pune University
Campus, Pune 411007, India
4 Observatoire de Genève, Université de Genève, 51 Ch. des Maillettes, 1290
Versoix, Switzerland
5 Department of Physics and Astronomy, University of Manitoba, Winnipeg, MB
R3T 2N2, Canada
6 Arecibo Observatory, NAIC, HC3 Box 53995, Arecibo, Puerto Rico, PR 00612,
USA

###### Abstract

Blueberry galaxies (BBs) are fainter, less massive, and lower redshift
counterparts of the Green pea galaxies. They are thought to be the nearest
analogues of the high redshift Lyman Alpha (Ly$\alpha$) emitters. We report
the interferometric imaging of Hi 21 cm emission from a Blueberry galaxy,
J1509+3731, at redshift, z = 0.03259, using the Giant Metrewave Radio
Telescope (GMRT). We find that this Blueberry galaxy has an Hi mass of
$M_{\text{H{\sc i}}}\approx 3\times 10^{8}\,M_{\odot}$ and an Hi-to-stellar
mass ratio $M_{\text{H{\sc i}}}/M_{*}\approx$ 2.4. Using SFR estimates from
the H$\beta$ emission line, we find that it has a short Hi depletion time
scale of $\approx 0.2$ Gyr, which indicates a significantly higher star-
formation efficiency compared to typical star-forming galaxies at the present
epoch. Interestingly, we find an offset of $\approx 2$ kpc between the peak of
the Hi 21 cm emission and the optical centre which suggests a merger event in
the past. Our study highlights the important role of mergers in triggering the
starburst in BBs and their role in the possible leakage of Lyman-$\alpha$ and
Lyman-continuum photons which is consistent with the previous studies on BB
galaxies.

###### keywords:

galaxies: dwarf — galaxies: high-redshift — galaxies: ISM — galaxies:
starburst — galaxies: star formation

##  1 Introduction

In the current understanding of the $\Lambda$-Cold Dark Matter ($\Lambda$CDM)
cosmology, the structure formation in the Universe follows a hierarchical
process. This bottom-up structure formation starts with small over-densities,
the dwarf galaxies. Later massive galaxies grow from assembling these smaller
galaxies through mergers and other accretion activities. A fraction of these
low-mass young star-forming galaxies are the Lyman Alpha (Ly$\alpha$) emitters
(LAEs); Lyman continuum (LyC) photons leaked by these LAEs are thought to be
the main contributors in reionizing the Universe (Bouwens et al., 2015;
Dressler et al., 2015; Drake et al., 2017; Dayal et al., 2020). Most known
LAEs are high redshift, chemically less evolved, compact galaxies with low
stellar mass, less dust, and a high ionization ratio. They are thought to be
the progenitors of present-day Milky Way like galaxies, and provide an
important probe to the earliest era of galaxy formation and the epoch of
reionization (Finkelstein et al., 2007; Gawiser et al., 2007; Nakajima et al.,
2012; Malhotra et al., 2012; Nakajima & Ouchi, 2014; Yang et al., 2017a;
Nakajima et al., 2020). Detailed study of the inter-stellar medium (ISM) of
LAEs is difficult because most LAEs are distant high redshift galaxies.
However, identification of low redshift LAE-analogues allows us to indirectly
probe the ISM properties of distant LAEs.

The Galaxy Zoo project (Lintott et al., 2008) identifies a class of green-
colored, small-sized objects from the Sloan Digital Sky Survey Data Release 7
(SDSS DR7; Abazajian, et al., 2009), known as the "Green Pea" galaxies (GPs).
These highly star-forming (SFR $\sim 10\;M_{\odot}yr^{-1}$) galaxies have high
specific star formation rates (sSFR $\sim 10$ Gyr-1) for their stellar masses
($\log[M_{*}/M_{\odot}]\sim 8-10$), low metallicities ($12+\log[O/H]\sim
8.05$) and high [O III]$\lambda$5007/[O II]$\lambda$3727 (O32) ratios
(Cardamone et al., 2009; Amorín, Pérez-Montero, & Vílchez, 2010; Amorín et
al., 2012; Jaskot & Oey, 2014; Izotov, Guseva, & Thuan, 2011). The bright
green color of these compact galaxies arises from very strong [O
III]$\lambda$5007 emission line with large equivalent widths ($\sim 1000$ Å)
which falls into the SDSS ‘r’ band. After the discovery of these GPs, several
studies have brought forward interesting similarities between these GPs and
the high redshift LAEs (Henry et al., 2015; Yang et al., 2016, 2017a; Izotov
et al., 2020). Jaskot & Oey (2014) suggest that GPs are optically thin
systems, which makes it possible for the LyC photons to escape. The results
from Henry et al. (2015) support this scenario indicating to low atomic
hydrogen (Hi) column densities in GPs. Outcomes of several studies on
Ly$\alpha$ profiles, Ly$\alpha$ and LyC escape fractions of GPs (Yang et al.,
2016, 2017a; Izotov et al., 2016a, b, 2018; Lofthouse, Houghton, & Kaviraj,
2017; McKinney et al., 2019) strengthened the hypothesis that GPs are local
analogues of high redshift LAEs. Importantly, it has been found that GPs can
leak a significant fraction of LyC photons, making them similar to early star-
forming galaxies believed to be responsible for reionization of the Universe
(e.g., Izotov et al., 2016a, b). Several high-redshift ($z>6$) galaxies,
recently detected by James Webb Space Telescope (JWST), show properties very
similar to GPs, thus increasing our confidence of them being good local
analogues to early star-forming galaxies (Schaerer et al., 2022; Rhoads et
al., 2023; Cameron et al., 2023; Matthee et al., 2022).

Yang et al. (2017b) (hereafter, Y17) present 40 starburst dwarf galaxies from
SDSS Data Release 12 (SDSS DR12; Alam et al., 2015) catalogue111
_https://www.sdss4.org/dr12/_ with redshifts $z\leq 0.05$ and similar
properties as GPs. The [O III]$\lambda$5007 emission lines from these galaxies
fall into the SDSS ‘g’ band (because of their lower redshifts), giving blue
colours to the composite images, and leading to the name “Blueberry” Galaxies
(BBs). Later it has become a convention to refer all these low redshift
analogues of LAEs as "Green Pea" galaxies. However, we follow the traditional
nomenclature and refer the lower redshift ($z<0.1$) galaxies with blue-
coloured composite images as "Blueberry" galaxies. These BBs represent
fainter, less massive, and lower redshift counterparts of the GPs, and are
thought to be the nearest analogues of the high redshift LAEs. A radio
continuum study of some of these BBs shows suppression in the radio-based SFR
compared to the optical SFR estimates (Sebastian & Bait, 2019).

Jiang et al. (2019) present a collection of 835 GPs from SDSS DR13 (Albareti
et al., 2017) within the redshift range $0.011<z<0.411$. Kanekar et al. (2021)
(hereafter, K21) select 40 galaxies from the sample of Jiang et al. (2019) at
$z<0.1$ and report the first detection of Hi 21 cm line emission in 19 GPs (we
refer to these as BBs in this article as per the traditional nomenclature)
along with estimates of their Hi masses using the Arecibo Telescope or the
Green Bank Telescope (GBT). They also find a lower Hi detection rate for
galaxies with higher O32 ratios. Purkayastha et al. (2022) present the first
interferometric Hi 21 cm imaging of a BB galaxy using the Giant Metrewave
Radio Telescope (GMRT), which suggests a galaxy-galaxy merger event to be the
reason for the starburst activity.

Here we present, Hi 21 cm imaging of a Blueberry galaxy, J1509+3731, at
redshift222We assume a flat $\Lambda$CDM cosmology with
($\Omega_{m},\Omega_{\Lambda},h$) = (0.315, 0.685, 0.674) (Planck
Collaboration et al., 2020), where $h$ is the dimensionless Hubble parameter
related to the Hubble Constant $H_{0}=100\,h\,$ km s-1 Mpc-1. z = 0.03259. We
describe the observations and data reduction process in Section 1. The
observed Hi properties of the Blueberry galaxy are described in Section 3. The
implications of the observed Hi properties, in the context of high redshift
LAEs, are discussed in Section 4. Finally, we summarize our results and
conclusions in Section 5.

##  2 Observations and data analysis

Figure 1: Radio continuum emission from BB10: The 1.4 GHz radio continuum
contours are overlaid on the DESI Legacy Survey g-band image. The contour
levels are at 4, 8, and 16 times the RMS noise of the radio image. The
synthesized beam of the radio image (2.5$\arcsec$ x 2.5$\arcsec$) is shown at
the bottom left corner. The radio emission is resolved with a source size of
$\approx 5\arcsec$, about the same as the optical diameter of the galaxy.

###  2.1 The target galaxy

We select J1509+3731 (hereafter BB10) with J2000 coordinates RA =
15h09m34.17s, DEC = +37d31$\arcmin$46.1$\arcsec$ from Y17 as the target galaxy
for Hi 21 cm observation with the GMRT. The redshift of BB10 is
$z=0.03259\,(\pm 0.00001)$, which corresponds to a distance of 144 Mpc. The
optical properties of BB10 are compiled in Table 1. It has a high specific
star formation rate ($\approx 13$ Gyr-1) and a very high O32 ratio $(\approx
15.7)$. In addition to these, its low metallicity and small optical size make
it an ideal candidate to study as a local analogue of high redshift LAEs.
Importantly, this galaxy shows strong Ly$\alpha$ emission and also evidence of
Ly$\alpha$ absorption in the Hubble Space Telescope (HST) Cosmic Origins
Spectrograph (COS) spectrum presented by Jaskot et al. (2019). The Ly$\alpha$
emission features show two peaks with a large Ly$\alpha$ peak velocity
separation of $\Delta v_{\mathrm{Ly\alpha}}$ $=400\pm 27$ km s-1and a weak
Ly$\alpha$ escape fraction of $f_{\mathrm{esc}}^{\mathrm{Ly\alpha}}$ $=0.05\pm
0.03$ (Jaskot et al., 2019). Previously, K21 targeted this galaxy using GBT
and provided an upper limit (3$\sigma$) on the Hi mass, $M_{\textsc{H{\sc
i}}}<4.4\times 10^{8}\,M_{\odot}$.

$M_{*}$ | SFR | O32 | $12+\log[O/H]$
---|---|---|---
$\times\;10^{8}\;M_{\odot}$ |  $M_{\odot}\;$ yr-1 |  |
1.26 | 1.61 | 15.67 | 7.87
Table 1: The properties of J1509+3731 (BB10) from Yang et al. (2017b). The
columns are i)stellar mass, ii) star formation rate, iii) O32 ratio, and iv)
metallicity.

Figure 2: [A, left] Hi 21 cm emission map of BB10 at a resolution of
$\mathbf{2.9\arcsec}$: The Hi contours, corresponding to the emission map
obtained by averaging six central velocity channels within the FWHM of the Hi
21 cm signal, are overlaid on the DESI Legacy Survey g-band image. The contour
levels are at 2, 3, and 4 times the RMS noise in the Hi 21 cm emission map.
The negative contours are shown in dashed lines. The synthesized beam is shown
at the bottom left corner. [B, right] Spatially integrated Hi 21 cm emission
spectrum of BB10 at a velocity resolution of 21.8 km s-1: The Hi 21 cm
spectrum has been integrated over the region within the outermost contour on
the left plot. The dotted line shows the spectral RMS noise.

###  2.2 Observations and Data Analysis

The target galaxy, BB10, was observed with the band-5 (1000-1450 MHz) system
of the upgraded GMRT (uGMRT), under the observing proposal DDTC095 (PI:
Chaitra A. Narayan), for a total on-source time of $\sim$9 hrs. The
observations were carried out with a bandwidth of 100 MHz (divided into 8192
channels) centered at 1380.00 MHz, using the GMRT Wideband Backend (GWB).
Radio sources 3C286 and 1602+334 were used as the flux and the phase
calibrator, respectively. The central square baselines, for the GWB system,
were affected by a power-offset333A detailed description of the issue can be
found at http://gmrt.ncra.tifr.res.in/gmrt_users/help/csq_baselines.html..
Hence the central square baselines were excluded from our analysis.

###  2.3 Radio continuum imaging

Radio interferometric imaging was carried out using the CASA (Common Astronomy
Software Applications) software package (McMullin et al., 2007). After gain
calibration and bandpass calibration using the flux calibrator, visibility
data were inspected for the presence of radio frequency interference (RFI).
Following the excision of RFI-affected data, visibilities were averaged to 781
kHz frequency channels for continuum imaging. Channels within $\pm 200$ km s-1
of the expected location of the Hi 21 cm line were not used for continuum
imaging. The continuum image was made using the task tclean, and was used to
self-calibrate the visibility data. The final continuum image has an RMS noise
of 12.4 $\mu$Jy beam-1 and an angular resolution of 2.5$\arcsec$ x
2.5$\arcsec$.

Radio continuum emission from BB10 was detected in our continuum image, as
shown in Figure 1 where we plot the 1.4 GHz radio continuum contours on the
optical g-band image 444https://www.legacysurvey.org/ from Dark Energy
Spectroscopic Instrument (DESI) Legacy imaging survey (Dey et al., 2019). We
used the CASA task IMFIT to estimate the properties of the emission source by
fitting a 2-dimensional elliptical Gaussian to the spatial brightness
distribution. A detailed description of the task, and the algorithms used, can
be found in the CASA documentation
555https://casadocs.readthedocs.io/en/latest/api/tt/casatasks.analysis.imfit.html.
The uncertainties on the flux density and the sky position of the peak are
estimated using the prescription of Condon et al. (1997).

The radio continuum emission is marginally resolved666The peak flux density
and the integrated flux density of the source differ by $\approx 2.5\sigma$.
It has a deconvolved size of 1.6″$\times$ 1.0″, which is smaller than the
synthesized beam. We verified these estimates by performing source extraction
using a different software package (PyBDSF; Mohan & Rafferty, 2015), which
yielded consistent results. with an integrated 1.4 GHz flux density of 340
$\pm$ 29 $\mu$Jy. Our estimate of the (integrated) flux density is consistent
with its previously published measurement ($345\pm 16\;\mu$Jy; Sebastian &
Bait, 2019). The peak flux density of the radio continuum emission (at an
angular resolution of 2.5″) is 261 $\mu$Jy beam-1, detected at a statistical
significance of (S/N)peak = 21 (i.e. peak / RMS noise).

The uncertainties on the sky position of the peak of the emission source is
given by (Condon et al., 1997, 1998)

$\rm\Delta RA\approx\Delta DEC\approx\frac{\Theta_{\rm FWHM}}{2(S/N)_{peak}}$
(1)

where $\Theta_{\rm FWHM}$ is the FWHM of the synthesized beam and
$\rm(S/N)_{peak}$ is the ratio of the peak flux density (not the spatially
integrated flux density) to RMS noise. For example, the sky position of a
source with a peak S/N = 5 can be measured with an accuracy of approximately a
tenth of the FWHM of the synthesized beam in each of RA and DEC (Condon et
al., 1998).

The peak of the radio continuum emission associated with BB10 is offset by
$0.7\arcsec\pm 0.1\arcsec$ (0.49 $\pm$ 0.06 kpc) from the peak of optical
emission (see Figure 1). The uncertainties on the position of the radio peak
is consistent with the relation in Equation 1. To investigate any systematic
spatial offset in our radio image, we compared the positions of bright compact
sources in our continuum image (within the FWHM of the GMRT primary beam) with
their reference positions in the NRAO VLA Sky Survey
(NVSS777https://www.cv.nrao.edu/nvss/; Condon et al., 1998) catalogue. The
individual offsets along RA and DEC, i.e. $\delta{\rm RA=RA-RA_{NVSS}}$ and
$\delta{\rm DEC=DEC-DEC_{NVSS}}$, are randomly distributed with both positive
and negative values, around their median values of $\delta{\rm
RA_{median}}=0.1\arcsec$ and $\delta{\rm DEC_{median}}=0.04\arcsec$. Note that
$\delta{\rm RA_{median}}$ and $\delta{\rm DEC_{median}}$, the median values of
the offsets along RA and DEC, represent the systematic offset between our
radio image and NVSS. This systematic offset is significantly smaller than the
offset between the radio and the optical peaks, and hence we do not apply any
correction for the systematic offset.

###  2.4 Spectroscopic imaging

The gain solutions from self-calibration were applied to visibilities at the
native frequency resolution (before channel averaging). After another round of
RFI excision, the calibrated visibilities were used to make Hi 21 cm data cube
with a velocity resolution of 21.8 km s-1. The Hi data cube was made using
visibilities with a UV upper limit of 60 k$\lambda$ (and Briggs weighting
scheme, with robust = 1.0), which yielded a spatial resolution of
$2.9\arcsec\times 2.9\arcsec$. The RMS noise at this angular resolution, per
21.8 km s-1 velocity channel, is $\approx 170\;\mu$Jy beam-1. A spectral
baseline was subtracted from each spatial pixel in the cube, where the
spectral baseline was obtained by fitting a second-order polynomial to the
spectra excluding channels within $\pm 200$ km s-1 of the expected location of
the Hi 21 cm line. Additional data cubes were made at lower spatial
resolutions, using smaller UV upper limits. For example, the cube made using
visibilities with a UV upper limit of 15 k$\lambda$ (and Briggs weighting
scheme, with robust = 1.0) has a spatial resolution of $7.4\arcsec\times
6.5\arcsec$ and an RMS noise of $\approx 280\;\mu$Jy beam-1 per 21.8 km s-1
velocity channel. Both the data cubes are CLEANed down to 1.5 times the rms
noise in each channel.

As mentioned earlier, the central square baselines (< 5 k$\lambda$) of GMRT
were not used for making the spectral cubes. The largest angular scale probed
by these cubes is thus $\approx 40\arcsec$. Hence, the Hi 21 cm emission can
be reliably mapped up to angular scales of $\lesssim 20\arcsec$ (spatial
scales of $\lesssim 15$ kpc at the redshift of BB10).

Figure 3: The Hi contours at a resolution of $7.4\arcsec\times 6.5\arcsec$
overlaid on the DESI Legacy Survey g-band image: The contour levels are at 2,
and 3 times the RMS noise in the Hi 21 cm map at this resolution. The negative
contours are shown in dashed lines. The synthesized beam is shown at the
bottom left corner.

##  3 Hi properties of BB10

Hi 21 cm emission from BB10 is detected in our data cube (at $2.9\arcsec\times
2.9\arcsec$ angular resolution). The Hi 21 cm spectrum was obtained by
spatially integrating over the source. The Hi 21 cm emission (moment-0) map
was obtained by averaging six velocity channels within the FWHM of the Hi 21
cm line (also above $1\sigma$). Figure 2[A] shows the Hi 21 cm emission map of
BB10, overlaid on the DESI Legacy Survey g-band image, and Figure 2[B] shows
its Hi 21 cm spectrum at a velocity resolution of 21.8 km s-1. The integrated
Hi 21 cm flux density is $S_{21}=60\pm 11$ mJy km s-1, which yields an Hi mass
of $M_{\text{H{\sc i}}}=(3.12\pm 0.57)\times 10^{8}\,M_{\odot}$. At an angular
resolution of $2.9\arcsec\times 2.9\arcsec$, the Hi 21 cm emission map has a
2$\sigma$ column density sensitivity of $N_{\text{H{\sc i}}}=1.9\times
10^{21}$ cm-2.

The Hi 21 cm spectra and emission maps obtained from data cubes made at
coarser spatial resolution yield Hi masses consistent with the above estimate.
The Hi 21 cm emission map at a spatial resolution of $7.4\arcsec\times
6.5\arcsec$ is shown in Figure 3, which has a 2$\sigma$ column density
sensitivity of $N_{\text{H{\sc i}}}=0.7\times 10^{21}$ cm-2. Table 2 compiles
the details of integrated Hi 21 cm flux density, Hi mass, and 2$\sigma$ Hi
column density for different beam sizes. As the central square baselines are
not included, it could not be possible to consider a bigger synthesized beam
and make data cubes with coarser resolution. Therefore, our estimate for the
Hi mass is formally a lower limit of the actual Hi mass of the galaxy.

Due to the lack of short baselines, our observations are not sensitive to Hi
emission on spatial scales $\gtrsim$ 15 kpc. However, we note that the non-
detection of Hi 21cm emission from BB10 with the GBT by K21 and the associated
upper limit on its Hi mass suggest that the galaxy does not have a significant
amount of Hi at larger spatial scales. The GBT upper limit — $M_{\text{H{\sc
i}}}<4.4\times 10^{8}M_{\odot}$ — and the (remarkably) tight Hi mass–size
relation (Wang et al., 2016) suggests a Hi disk size (diameter) of
$d_{HI}\lesssim 15$ kpc. It thus appears unlikely that the actual Hi mass of
BB10 is significantly larger than our estimate.

The Hi depletion timescale, $t_{\text{dep}}=M_{\text{H{\sc i}}}/SFR$,
indicates how long a galaxy can continue forming stars without refilling its
gas reservoir. BB10 has a short Hi depletion timescale of $\approx 0.2$ Gyr.

Beam | $S_{\text{int}}$ | $M_{\text{H{\sc i}}}$ | $N_{\text{H{\sc i}}}$
---|---|---|---
($\arcsec\times\arcsec$) | (mJy km s-1) | ($\times\;10^{8}\;M_{\odot}$) | ($\times 10^{21}$ cm-2)
2.9 $\times$ 2.9 | 60 $\pm$ 11 | 3.12 $\pm$ 0.57 | 1.9
7.4 $\times$ 6.5 | 68 $\pm$ 14 | 3.53 $\pm$ 0.73 | 0.7
Table 2: The values of Hi mass ($M_{\text{H{\sc i}}}$) and 2$\sigma$ Hi column
density sensitivity (integrated over the Hi 21 cm emission line with velocity
width $\approx$ 130.8 km s-1) at different angular resolutions.

An offset of $2.88\arcsec\pm 0.70\arcsec$, i.e. $1.94\pm 0.48$ kpc, is seen
between the peak of the optical emission and the peak of Hi distribution888The
peak of the Hi distribution was localized by fitting a 2D Gaussian to the
emission source in the moment-0 map. The Hi contours are elongated along the
optical extension(see Figure 2[A]). The offset is seen in the Hi 21 cm maps at
coarser spatial resolutions as well. This offset is significantly larger than
the offset between the peak of optical emission and that of radio continuum
emission (see Figure 1).

Figure 4: $\mathbf{M_{\text{H{\sc i}}}/M_{*}}$ [top panels] and Hi depletion
time scale [bottom panels] of galaxies from K21 and BB10 as functions of their
stellar mass [left] and O32 [right]: The grey filled circles show Hi 21 cm
detections of Blueberry galaxies and the grey arrows show 3$\sigma$ upper
limits for non-detections from Kanekar et al. (2021). The blue star represents
BB10.

##  4 Discussion

We compare the Hi properties of BB10 with the BBs from K21. We note that BB10
was also included in the sample of K21. Out of a sample of 40 BBs, K21 provide
estimates of Hi masses for 19 Hi detections and upper limits for the remaining
21 including BB10. Our Hi mass estimate for BB10 is consistent with the upper
limit provided by K21 ($4.4\times 10^{8}\,M_{\odot}$). Figure 4 shows the Hi-
to-stellar mass ratio, $M_{\text{H{\sc i}}}/M_{*}$, and the Hi depletion
timescale, $t_{\text{dep}}=M_{\text{H{\sc i}}}/SFR$, against the stellar mass
and O32 ratio for the BBs from the sample of K21 along with BB10. The grey
filled circles represent the Hi detections and the grey downward arrows show
the upper limits for the non-detections from K21. In all four panels, BB10 is
shown as a blue star. As can be seen in the top left and bottom left panels of
Figure 4, BB10 has a somewhat similar Hi-to-stellar mass ratio and Hi
depletion timescale as the BBs from K21 with similar stellar masses. However,
it is on the lower side of the Hi-to-stellar mass ratio and Hi depletion
timescale distribution. The top right and bottom right panels of Figure 4
show, as also noted by K21, that galaxies with higher O32 values have
preferentially lower Hi-to-stellar mass ratio, and thus also typically lower
depletion timescales. Our estimates of the Hi-to-stellar mass ratio and Hi
depletion timescale of BB10 are consistent with this scenario. We note that
BB10 has a higher O32 ratio than the other BBs with Hi 21 cm detections in the
K21 sample. However, more measurements of the Hi properties of galaxies in the
low $M_{*}$ \- high O32 space are needed to make a robust conclusion.

Purkayastha et al. (2022) find Hi 21 cm emission from a broken ring-like
structure around a BB galaxy and suggest a merger event, with a companion
galaxy separated by $\approx 4.7$ kpc, to be responsible for the starburst
activity in that BB galaxy. Recently, Le Reste et al. (2023) present Hi
imaging of a well-studied local LyC emitter (LCE), Haro 11. Both these studies
highlight that mergers can cause a significant displacement of the neutral gas
which can facilitate the escape of Lyman-$\alpha$ and possibly also Lyman-
continuum photons.

We do not find any clear evidence of a close companion galaxy (or a merging
galaxy) from the DESI Legacy Survey (Dey et al., 2019) optical image of BB10.
However, as mentioned earlier, we do notice a clear spatial offset between the
peak of the Hi 21 cm emission and the star-forming region (the brightest part
of the optical emission). In addition, the Hi distribution extends over the
star-forming region as well. We note that optical emission in the DESI image
also appears asymmetric as can be seen in Figure 1. While the brightest part
of the optical emission is offset from the densest part of the Hi spatial
distribution, an extended ‘tail’ of the optical emission coincides with the
peak of the Hi distribution (see Figure 2[A]). Both of these might be caused
by a dwarf-dwarf merger. And in that case, BB10 is perhaps in a late-stage
merger where the two galaxies have nearly coalesced with each other.
Interestingly, the COS-NUV image of BB10, which has a higher spatial
resolution than the DESI image, shows relatively extended emission with two or
more UV emitting knots within the bright optical emission (see Figure 3 in
Jaskot et al., 2019). It is well known that mergers can trigger a starburst
(e.g., Luo, Yang, & Zhang, 2014; Knapen & Cisternas, 2015; Stierwalt et al.,
2015; Moreno et al., 2021). Thus a dwarf-dwarf merger scenario can altogether
explain the origin of the offset in Hi 21 cm emission, multiple UV knots, and
the origin of the starburst in BB10.

Finally, compact starbursts can have strong outflows due to stellar/supernova-
driven feedback (e.g., Weiner et al., 2009; Heckman et al., 2011; Chisholm et
al., 2017). Such outflows can push significant amounts of gas from the star-
forming knot of BB10 leading to a paucity of gas in that region.

Using a sample of $z\sim 0.3$ high O32 star-forming galaxies, Izotov et al.
(2018) propose an empirical relation between $\Delta v_{\mathrm{Ly\alpha}}$
and the escape fraction of LyC photons ($f_{\mathrm{esc}}^{\mathrm{LyC}}$).
Further, Flury et al. (2022) using a multi-parameter analysis finds that
$\Delta v_{\mathrm{Ly\alpha}}$ is the best in-direct indicator of
$f_{\mathrm{esc}}^{\mathrm{LyC}}$. BB10 has properties which lie within the
range of galaxies studied in Izotov et al. (2018), in terms of the O32 ratio,
stellar mass, metallicity, SFR, etc. We can thus get an indirect estimate of
$f_{\mathrm{esc}}^{\mathrm{LyC}}$ from the observed $\Delta
v_{\mathrm{Ly\alpha}}$ $=400\pm 27$ km s-1from Jaskot et al. (2019) and the
empirical relation from Izotov et al. (2018). We find an indirect
$f_{\mathrm{esc}}^{\mathrm{LyC}}$ $\approx$ 0.034 for BB10. Thus BB10 is a
weak Ly$\alpha$ and possibly a weak LyC leaker. Indeed mergers can lead to
favourable conditions such as triggering an intense starburst (and thus
increasing the intrinsic LyC production) and simultaneously displacing a
significant amount of Hi gas through tidal interactions to open low column
density channels which can facilitate the leakage of Ly$\alpha$ and possibly
LyC photons (Witten et al., 2024). However, Hi 21 cm imaging of a sample of
BBs would be needed to help us obtain a better understanding of the role of
galaxy mergers, in particular the complex role of orientation and geometrical
effects, in the escape of LyC photons. Targeted simulations of merging
galaxies with realistic feedback will also be helpful in understanding the
role of mergers in LyC escape from these young star-forming galaxies (e.g.,
Renaud, Bournaud, & Duc, 2015).

##  5 Conclusion & Summary

In this article, we present the results obtained from Hi 21 cm imaging of a
Blueberry galaxy at a redshift of $z=0.03259$, using uGMRT band-5 receivers.
This Blueberry galaxy has an Hi mass of $M_{\text{H{\sc i}}}=(3.12\pm
0.57)\times 10^{8}\,M_{\odot}$, and an Hi-to-stellar mass ratio
$M_{\text{H{\sc i}}}/M_{*}\approx$ 2.4. Using SFR estimates from the H$\beta$
emission line (Yang et al., 2017b), we find that it has a short Hi depletion
time scale of $\approx 0.2$ Gyr. Comparing the Hi 21 cm emission map with the
DESI Legacy Survey optical image of BB10, we find an offset of $(1.94\pm
0.48)$ kpc between the centres of Hi 21 cm emission and star-forming region.
Such an offset could be a sign of a merger event which can also trigger a
starburst. Our study along with Purkayastha et al. (2022) and Le Reste et al.
(2023) highlight the role of dwarf galaxy mergers in the leakage of ionizing
photons, and thus their role in cosmic reionization.

## ACKNOWLEDGMENTS

We are grateful to the anonymous reviewer and the editor for feedback and
suggestions on the manuscript which have improved the article significantly.
We thank the staff of the GMRT who have made these observations possible. The
GMRT is run by the National Center for Radio Astrophysics of the Tata
Institute of Fundamental Research. O. B. is supported by the AstroSignals
Sinergia Project funded by the Swiss National Science Foundation. S.D and
C.A.N. also acknowledge the Department of Atomic Energy for funding support,
under project 12-R&D-TFR-5.02-0700.

The Legacy Surveys consist of three individual and complementary projects: the
Dark Energy Camera Legacy Survey (DECaLS; Proposal ID #2014B-0404; PIs: David
Schlegel and Arjun Dey), the Beijing-Arizona Sky Survey (BASS; NOAO Prop. ID
#2015A-0801; PIs: Zhou Xu and Xiaohui Fan), and the Mayall z-band Legacy
Survey (MzLS; Prop. ID #2016A-0453; PI: Arjun Dey). DECaLS, BASS and MzLS
together include data obtained, respectively, at the Blanco telescope, Cerro
Tololo Inter-American Observatory, NSF’s NOIRLab; the Bok telescope, Steward
Observatory, University of Arizona; and the Mayall telescope, Kitt Peak
National Observatory, NOIRLab. Pipeline processing and analyses of the data
were supported by NOIRLab and the Lawrence Berkeley National Laboratory
(LBNL). The Legacy Surveys project is honored to be permitted to conduct
astronomical research on Iolkam Du’ag (Kitt Peak), a mountain with particular
significance to the Tohono O’odham Nation.

NOIRLab is operated by the Association of Universities for Research in
Astronomy (AURA) under a cooperative agreement with the National Science
Foundation. LBNL is managed by the Regents of the University of California
under contract to the U.S. Department of Energy.

This project used data obtained with the Dark Energy Camera (DECam), which was
constructed by the Dark Energy Survey (DES) collaboration. Funding for the DES
Projects has been provided by the U.S. Department of Energy, the U.S. National
Science Foundation, the Ministry of Science and Education of Spain, the
Science and Technology Facilities Council of the United Kingdom, the Higher
Education Funding Council for England, the National Center for Supercomputing
Applications at the University of Illinois at Urbana-Champaign, the Kavli
Institute of Cosmological Physics at the University of Chicago, Center for
Cosmology and Astro-Particle Physics at the Ohio State University, the
Mitchell Institute for Fundamental Physics and Astronomy at Texas A&M
University, Financiadora de Estudos e Projetos, Fundacao Carlos Chagas Filho
de Amparo, Financiadora de Estudos e Projetos, Fundacao Carlos Chagas Filho de
Amparo a Pesquisa do Estado do Rio de Janeiro, Conselho Nacional de
Desenvolvimento Cientifico e Tecnologico and the Ministerio da Ciencia,
Tecnologia e Inovacao, the Deutsche Forschungsgemeinschaft and the
Collaborating Institutions in the Dark Energy Survey. The Collaborating
Institutions are Argonne National Laboratory, the University of California at
Santa Cruz, the University of Cambridge, Centro de Investigaciones
Energeticas, Medioambientales y Tecnologicas-Madrid, the University of
Chicago, University College London, the DES-Brazil Consortium, the University
of Edinburgh, the Eidgenossische Technische Hochschule (ETH) Zurich, Fermi
National Accelerator Laboratory, the University of Illinois at Urbana-
Champaign, the Institut de Ciencies de l’Espai (IEEC/CSIC), the Institut de
Fisica d’Altes Energies, Lawrence Berkeley National Laboratory, the Ludwig
Maximilians Universitat Munchen and the associated Excellence Cluster
Universe, the University of Michigan, NSF’s NOIRLab, the University of
Nottingham, the Ohio State University, the University of Pennsylvania, the
University of Portsmouth, SLAC National Accelerator Laboratory, Stanford
University, the University of Sussex, and Texas A&M University.

BASS is a key project of the Telescope Access Program (TAP), which has been
funded by the National Astronomical Observatories of China, the Chinese
Academy of Sciences (the Strategic Priority Research Program “The Emergence of
Cosmological Structures” Grant # XDB09000000), and the Special Fund for
Astronomy from the Ministry of Finance. The BASS is also supported by the
External Cooperation Program of Chinese Academy of Sciences (Grant #
114A11KYSB20160057), and Chinese National Natural Science Foundation (Grant #
12120101003, # 11433005).

The Legacy Survey team makes use of data products from the Near-Earth Object
Wide-field Infrared Survey Explorer (NEOWISE), which is a project of the Jet
Propulsion Laboratory/California Institute of Technology. NEOWISE is funded by
the National Aeronautics and Space Administration.

The Legacy Surveys imaging of the DESI footprint is supported by the Director,
Office of Science, Office of High Energy Physics of the U.S. Department of
Energy under Contract No. DE-AC02-05CH1123, by the National Energy Research
Scientific Computing Center, a DOE Office of Science User Facility under the
same contract; and by the U.S. National Science Foundation, Division of
Astronomical Sciences under Contract No. AST-0950945 to NOAO.

## Data Availability

Raw data used in this work are available in the GMRT Online Archive
(_http://naps.ncra.tifr.res.in/goa/_) under proposal ID: DDTC095. Reduced data
are available with the corresponding author and will be shared on reasonable
request.

## References

  * Abazajian, et al. (2009) Abazajian, K. N., et al. 2009,ApJS,182, 543
  * Alam et al. (2015) Alam S., Albareti F. D., Allende Prieto C., Anders F., Anderson S. F., Anderton T., Andrews B. H., et al., 2015, ApJS, 219, 12.
  * Albareti et al. (2017) Albareti F. D., Allende Prieto C., Almeida A., Anders F., Anderson S., Andrews B. H., Aragón-Salamanca A., et al., 2017, ApJS, 233, 25.
  * Amorín, Pérez-Montero, & Vílchez (2010) Amorín R. O., Pérez-Montero E., Vílchez J. M., 2010, ApJL, 715, L128.
  * Amorín et al. (2012) Amorín R., Pérez-Montero E., Vílchez J. M., Papaderos P., 2012, ApJ, 749, 185.
  * Bouwens et al. (2015) Bouwens R. J., Illingworth G. D., Oesch P. A., Trenti M., Labbé I., Bradley L., Carollo M., et al., 2015, ApJ, 803, 34.
  * Cameron et al. (2023) Cameron A. J., Saxena A., Bunker A. J., D’Eugenio F., Carniani S., Maiolino R., Curtis-Lake E., et al., 2023, arXiv, arXiv:2302.04298.
  * Cardamone et al. (2009) Cardamone C., Schawinski K., Sarzi M., Bamford S. P., Bennert N., Urry C. M., Lintott C., et al., 2009, MNRAS, 399, 1191.
  * Chisholm et al. (2017) Chisholm J., Orlitová I., Schaerer D., Verhamme A., Worseck G., Izotov Y. I., Thuan T. X., et al., 2017, A&A, 605, A67.
  * Condon et al. (1997) Condon J. J., 1997, PASP, 109, 166.
  * Condon et al. (1998) Condon J. J., Cotton W. D., Greisen E. W., Yin Q. F., Perley R. A., Taylor G. B., Broderick J. J., 1998, AJ, 115, 1693.
  * Dayal et al. (2020) Dayal P., Volonteri M., Choudhury T. R., Schneider R., Trebitsch M., Gnedin N. Y., Atek H., et al., 2020, MNRAS, 495, 3065.
  * Dey et al. (2019) Dey A., Schlegel D. J., Lang D., Blum R., Burleigh K., Fan X., Findlay J. R., et al., 2019, AJ, 157, 168.
  * Drake et al. (2017) Drake A. B., Garel T., Wisotzki L., Leclercq F., Hashimoto T., Richard J., Bacon R., et al., 2017, A&A, 608, A6.
  * Dressler et al. (2015) Dressler A., Henry A., Martin C. L., Sawicki M., McCarthy P., Villaneuva E., 2015, ApJ, 806, 19.
  * Finkelstein et al. (2007) Finkelstein S. L., Rhoads J. E., Malhotra S., Pirzkal N., Wang J., 2007, ApJ, 660, 1023.
  * Flury et al. (2022) Flury S. R., Jaskot A. E., Ferguson H. C., Worseck G., Makan K., Chisholm J., Saldana-Lopez A., et al., 2022, ApJ, 930, 126.
  * Gawiser et al. (2007) Gawiser E., Francke H., Lai K., Schawinski K., Gronwall C., Ciardullo R., Quadri R., et al., 2007, ApJ, 671, 278.
  * Heckman et al. (2011) Heckman T. M., Borthakur S., Overzier R., Kauffmann G., Basu-Zych A., Leitherer C., Sembach K., et al., 2011, ApJ, 730, 5.
  * Heintz et al. (2023a) Heintz K. E., Giménez-Arteaga C., Fujimoto S., Brammer G., Espada D., Gillman S., González-López J., et al., 2023, ApJL, 944, L30.
  * Heintz et al. (2023b) Heintz K. E., Watson D., Brammer G., Vejlgaard S., Hutter A., Strait V. B., Matthee J., et al., 2023, arXiv, arXiv:2306.00647.
  * Henry et al. (2015) Henry A., Scarlata C., Martin C. L., Erb D., 2015, ApJ, 809, 19.
  * Izotov, Guseva, & Thuan (2011) Izotov Y. I., Guseva N. G., Thuan T. X., 2011, ApJ, 728, 161.
  * Izotov et al. (2016a) Izotov Y. I., Orlitová I., Schaerer D., Thuan T. X., Verhamme A., Guseva N. G., Worseck G., 2016, Natur, 529, 178.
  * Izotov et al. (2016b) Izotov Y. I., Schaerer D., Thuan T. X., Worseck G., Guseva N. G., Orlitová I., Verhamme A., 2016, MNRAS, 461, 3683.
  * Izotov et al. (2018) Izotov Y. I., Schaerer D., Worseck G., Guseva N. G., Thuan T. X., Verhamme A., Orlitová I., et al., 2018, MNRAS, 474, 4514.
  * Izotov et al. (2020) Izotov Y. I., Schaerer D., Worseck G., Verhamme A., Guseva N. G., Thuan T. X., Orlitová I., et al., 2020, MNRAS, 491, 468.
  * Jaskot & Oey (2014) Jaskot A. E., Oey M. S., 2014, ApJL, 791, L19.
  * Jaskot et al. (2019) Jaskot A. E., Dowd T., Oey M. S., Scarlata C., McKinney J., 2019, ApJ, 885, 96.
  * Jiang et al. (2019) Jiang T., Malhotra S., Rhoads J. E., Yang H., 2019, ApJ, 872, 145.
  * Kanekar et al. (2021) Kanekar N., Ghosh T., Rhoads J., Malhotra S., Harish S., Chengalur J. N., Jones K. M., 2021, ApJL, 913, L15.
  * Knapen & Cisternas (2015) Knapen J. H., Cisternas M., 2015, ApJL, 807, L16.
  * Le Reste et al. (2023) Le Reste A., Cannon J. M., Hayes M. J., Inoue J. L., Kepley A. A., Melinder J., Menacho V., et al., 2023, MNRAS.tmp.
  * Lintott et al. (2008) Lintott C. J., Schawinski K., Slosar A., Land K., Bamford S., Thomas D., Raddick M. J., et al., 2008, MNRAS, 389, 1179.
  * Lofthouse, Houghton, & Kaviraj (2017) Lofthouse E. K., Houghton R. C. W., Kaviraj S., 2017, MNRAS, 471, 2311.
  * Luo, Yang, & Zhang (2014) Luo W., Yang X., Zhang Y., 2014, ApJL, 789, L16.
  * Malhotra et al. (2012) Malhotra S., Rhoads J. E., Finkelstein S. L., Hathi N., Nilsson K., McLinden E., Pirzkal N., 2012, ApJL, 750, L36.
  * Matthee et al. (2022) Matthee J., Mackenzie R., Simcoe R. A., Kashino D., Lilly S. J., Bordoloi R., Eilers A.-C., 2022, arXiv, arXiv:2211.08255.
  * McKinney et al. (2019) McKinney J. H., Jaskot A. E., Oey M. S., Yun M. S., Dowd T., Lowenthal J. D., 2019, ApJ, 874, 52.
  * McMullin et al. (2007) McMullin J. P., Waters B., Schiebel D., Young W., Golap K., 2007, ASPC, 376, 127
  * Mohan & Rafferty (2015) Mohan N., Rafferty D., 2015, ascl.soft. ascl:1502.007
  * Moreno et al. (2021) Moreno J., Torrey P., Ellison S. L., Patton D. R., Bottrell C., Bluck A. F. L., Hani M. H., et al., 2021, MNRAS, 503, 3113.
  * Nakajima et al. (2012) Nakajima K., Ouchi M., Shimasaku K., Ono Y., Lee J. C., Foucaud S., Ly C., et al., 2012, ApJ, 745, 12.
  * Nakajima & Ouchi (2014) Nakajima K., Ouchi M., 2014, MNRAS, 442, 900.
  * Nakajima et al. (2020) Nakajima K., Ellis R. S., Robertson B. E., Tang M., Stark D. P., 2020, ApJ, 889, 161.
  * Planck Collaboration et al. (2020) Planck Collaboration, Aghanim N., Akrami Y., Ashdown M., Aumont J., Baccigalupi C., Ballardini M., et al., 2020, A&A, 641, A6.
  * Purkayastha et al. (2022) Purkayastha S., Kanekar N., Chengalur J. N., Malhotra S., Rhoads J., Ghosh T., 2022, ApJL, 933, L11.
  * Renaud, Bournaud, & Duc (2015) Renaud F., Bournaud F., Duc P.-A., 2015, MNRAS, 446, 2038.
  * Rhoads et al. (2023) Rhoads J. E., Wold I. G. B., Harish S., Kim K. J., Pharo J., Malhotra S., Gabrielpillai A., et al., 2023, ApJL, 942, L14.
  * Sebastian & Bait (2019) Sebastian B., Bait O., 2019, ApJL, 882, L19.
  * Schaerer et al. (2022) Schaerer D., Marques-Chaves R., Barrufet L., Oesch P., Izotov Y. I., Naidu R., Guseva N. G., et al., 2022, A&A, 665, L4.
  * Stierwalt et al. (2015) Stierwalt S., Besla G., Patton D., Johnson K., Kallivayalil N., Putman M., Privon G., et al., 2015, ApJ, 805, 2.
  * Wang et al. (2016) Wang J., Koribalski B. S., Serra P., van der Hulst T., Roychowdhury S., Kamphuis P., Chengalur J. N., 2016, MNRAS, 460, 2143.
  * Weiner et al. (2009) Weiner B. J., Coil A. L., Prochaska J. X., Newman J. A., Cooper M. C., Bundy K., Conselice C. J., et al., 2009, ApJ, 692, 187.
  * Witten et al. (2024) Witten C., Laporte N., Martin-Alvarez S., Sijacki D., Yuan Y., Haehnelt M. G., Baker W. M., et al., 2024, NatAs, 8, 384.
  * Yang et al. (2016) Yang H., Malhotra S., Gronke M., Rhoads J. E., Dijkstra M., Jaskot A., Zheng Z., et al., 2016, ApJ, 820, 130.
  * Yang et al. (2017a) Yang H., Malhotra S., Gronke M., Rhoads J. E., Leitherer C., Wofford A., Jiang T., et al., 2017, ApJ, 844, 171.
  * Yang et al. (2017b) Yang H., Malhotra S., Rhoads J. E., Wang J., 2017, ApJ, 847, 38.