#  Simpson’s Paradox and Lagging Progress in Completion Trends of
Underrepresented Students in Computer Science

John Mason Taylor<EMAIL_ADDRESS>Furman
UniversityGreenvilleSCUSA29613 ,  Rebecca Drucker<EMAIL_ADDRESS>Furman UniversityGreenvilleSCUSA29613 ,  Chris Alvin<EMAIL_ADDRESS>Furman UniversityGreenvilleSCUSA29613 and  Syed Fahad Sultan
<EMAIL_ADDRESS>Furman UniversityGreenvilleSCUSA29613

(2023)

###### Abstract.

It is imperative for the Computer Science (CS) community to ensure active
participation and success of students from diverse backgrounds. This work
compares CS to other areas of study with respect to success of students from
three underrepresented groups: Women, Black and Hispanic or Latino. Using a
data-driven approach, we show that trends of success over the years for
underrepresented groups in CS are lagging behind other disciplines. Completion
of CS programs by Black students in particular shows an alarming regression in
the years 2011 through 2019. This national level decline is most concentrated
in the Southeast of the United States and seems to be driven mostly by a small
number of institutes that produce a large number of graduates. We strongly
believe that more data-driven studies in this area are necessary to make
progress towards a more equitable and inclusive CS community. Without an
understanding of underlying dynamics, policy makers and practitioners will be
unable to make informed decisions about how and where to allocate resources to
address the problem.

IPEDS Data, Underrepresented Groups, Analytics, Computing education,
Completions

††copyright: acmcopyright††journalyear: 2023††doi: XXXXXXX.XXXXXXX††price:
15.00††isbn: 978-1-4503-XXXX-X/18/06††ccs: Social and professional
topics††ccs: Social and professional topics Race and ethnicity††ccs: Social
and professional topics Gender††ccs: Applied computing Education

##  1\. Introduction

The Computer Science (CS) student body of today is going to build and shape
our futures for tomorrow. To ensure that the benefits of future technology
extend to all, it is important that the field of CS is inclusive, equitable
and reflects the diverse perspectives and talents of all individuals.

Diversity, Equity, and Inclusion (DEI) in CS is a topic of growing interest in
the United States (Denner and Campe, 2023; Sax et al., 2017; Ibe et al.,
2018). However, in our review of the literature, we found a dearth of DEI
studies in CS that take a data-driven approach. Most of the existing work is
qualitative and prescriptive in nature (Kulkarni et al., 2018; Larsen and
Stubbs, 2005). To address this gap in Computing Education literature, in this
work, we use 9 years of Integrated Post-secondary Education Data System
(IPEDS) Completions (Statistics, 2012) data to better understand key trends in
success or struggles of underrepresented students in Computer Science.

Our results show that the completion trends in CS are lagging behind those in
other areas of study (i.e., non-CS) for students from three key
underrepresented demographic groups: Women, Black and Hispanic or Latino. For
Black students, we observe an alarmingly sharp decline in CS completions
during the same period when no such decline is observed in other areas of
study. We show that this is despite the fact that most institutions are
individually exhibiting a positive trend in CS completions for Black students.
This is an example of Simpson’s Paradox: a phenomenon of different or opposite
patterns observed at two different scales. We explain that a few institutions
with a high number of total completions disproportionately influence national
level trends in CS completions for Black students. Negative trends in CS for
Black students are observed to be concentrated primarily in the Southern
United States with a large cluster in the Southeast region.

For Women, we observe that the completions in CS are far behind completions in
other areas of study but the trend over the years is consistently positive.
For Hispanic or Latino students, we report that even though completion trends
are positive, they are still lagging behind trends observed in other areas of
study. Altogether, our results strongly suggest a call to action for more
data-driven research of underrepresented groups in CS and the broader topic of
DEI in CS. Only with a better understanding of the underlying dynamics can we
hope to make progress toward a more equitable and inclusive CS community.

##  2\. Methodology

Data. The Integrated Postsecondary Education Data System (IPEDS) is a set of
interrelated surveys conducted annually by the National Center for Education
Statistics (NCES), a part of the Institute for Education Sciences within the
United States Department of Education. The _Completions survey_ within this
system is designed to collect information on the number and types of degrees
awarded by all accredited U.S. postsecondary institutions as well as
information on the characteristics of degree recipients. The data is publicly
available to download from the IPEDS website (Statistics, 2012).

In this work we analyze $9$ years of data (2011-2019) from the Completions
survey to understand the success of underrepresented students in CS. These
nine years were selected to minimize any confounding effects that might be
attributed to the financial crisis starting in 2008 and the COVID-19 pandemic
beginning in 2020.

The data we use in this paper are from the survey table titled _Awards/degrees
conferred by program, award level, race/ethnicity, and gender_. We
specifically focus on undergraduate level degrees (Award Level code $5$). We
define ‘Computer Science’ as degrees awarded under Classification of
Instructional Programs (CIP) code $11$ titled _Computer and Information
Sciences and Support Services_ in the data.

Analyses. Our analyses focus on the following three underrepresented groups:
Women (variable=CTOTALW), Black (variable= CBKAAT) and Hispanic or Latino
(variable=CHISPT). Values for these three variables are self-reported by
students and communicated to NCES by institutions and are in the IPEDS
Completions survey tables at the institutional level. With the gender
category, there was no way to report anything other than ‘Men’ and ‘Women’ in
IPEDS prior to $2022$. A question was added in the 2022-23 survey to capture
the total number of non-binary students and the total number of students for
whom gender is unknown (FAQ, 2023). We acknowledge the limitations of this
binary definition of gender for past data and welcome expansion of gender
categories going forward.

Counts were computed using a sum of the relevant variable over all
institutions for a given year. Percentages were computed by dividing relevant
sums by CTOTALT: total of students awarded degrees with CIPCODE in range
$[11,12)$.

We use the Pearson correlation coefficient $r$ between _count of completions_
with years $[2011,2019]$ to measure trends over time; similarly for
_percentage of completions_ over years. $r$ measures how two continuous
variables co-vary and indicate a linear relationship as a number between $-1$
(perfectly negatively correlated) to $0$ (not correlated) to $1$ (perfectly
positively correlated).

We use Kernel Density Estimation (KDE) to identify regions of high density in
bivariate visualizations. Given a random variable $v$, KDE estimates the
probability density function (PDF) non-parametrically. Conceptually, KDE is a
smoothed version of a histogram for $v$.

Source code. Complete analyses and source code is available at
https://github.com/REDACTED.

(a) Absolute count of completions. These data have been log-scaled and all
vertical axes use the same scale.

(b) Percentage of total completions

Figure 1. National completion trends in CS and other areas of study from 2011
to 2019 for Women, Black, and Hispanic or Latino students.

Figure 2. Comparing CS completion trends with completion trends in other areas
of study using $r_{CS}-r_{other}$. Each dashed line is a mean for the
corresponding distribution ($\mu_{r_{CS}-r_{other}}$) while each dotted line
corresponds to National statistics of $r_{CS}-r_{other}$ as reported in Figure
1(b). Frequencies of all institutions with $r_{CS}-r_{other}>0$ are colored
orange; $r_{CS}-r_{other}<0$ are blue.

##  3\. Results

In this section, we present our findings from analysis of 9 years of IPEDS
completions data in Computer Science. In subsection 3.1, we argue for the use
of relative counts (i.e. percentages) over absolute counts by demonstrating
how completion trends using percentages are more informative and meaningful.
In subsection 3.2, we demonstrate the lagging progress of Women, Black and
Hispanic or Latino students in CS when compared to other areas of study (i.e.,
non-CS).

We continue in subsection 3.3 by comparing trends observed at the national
level to trends observed at the level of individual institutions. Our results
demonstrate that for Black and Hispanic or Latino students, completion trends
disappear or reverse at these two scales thus indicating the occurrence of
Simpson’s paradox. We then attempt to explain underlying reasons behind the
paradox.

Finally, we present the geographic distribution of completion trends in CS
across the United States, reporting spatial clusters of declining trends for
students from underrepresented groups.

###  3.1. Absolute vs. Relative Counts

Table 1. National-level and Institutional-level completion trends for CS and other areas of study.  | Demographic | $r_{CS}$ | $r_{other}$ | $r_{CS}-r_{other}$
---|---|---|---|---
National | National | National | Institutional Average
1 | Women | 0.91 | 0.85 | 0.06 | 0.07
2 | Black | -0.93 | -0.06 | -0.87 | -0.06
3 | Hispanic or Latino | 0.97 | 1 | -0.03 | -0.45

In comparing CS with other areas of study, we are interested not only in the
gap between completions for students from underrepresented groups, but also in
the trends over time.

Figure 1 depicts national level trends in CS completions as they compare to
other areas of study (Non-CS). Figure 1(a) presents completion trends computed
as absolute counts while Figure 1(b) completion trends with percentages.
Absolute counts are computed as the sum of all completions for a given year
and demographic group. It can be observed that when computed as absolute
counts, completions across all underrepresented groups exhibit strong positive
trends for all demographic groups in CS as well as other areas of study.

Used as a measure of participation or success, absolute counts do not lead to
fair comparisons. In Figure 1(a), CS completions are compared to total
completions for all other areas of study. The wide gap between the CS line
(orange) and the other line (blue) is by definition; hence, the vertical axes
of all panels in Figure 1(a) are log-scaled. Figure 1(b) presents completion
trends as percentages, computed as the sum of completions by students from a
demographic group in a given year divided by the sum of all completions that
year. It is now the case in Figure 1(b) that the gap between the trend lines
for the two groups is meaningful. It is also the case that the trend lines
individually exhibit more interesting patterns in stark contrast to absolute
counts in Figure 1(a).

When measured using absolute counts, CS exhibits positive trends for all three
underrepresented groups as shown in Figure 1(a). However, when using relative
counts, the CS trends weaken, and in the case of Black students, a reversal
occurs. We conclude that even though CS is seeing an increase in the absolute
count of completions overall, the _share of underrepresented students in these
completions is not increasing at the same rate_ , if at all. In fact, as shown
in Table 1, the share of Black students in CS completions experienced a sharp
decline $r_{CS\wedge Black}=-0.93$ during the same time when their share in
other areas of study remained relatively stable $r_{other\wedge Black}=-0.06$.

Figure 3. Scatterplot of institutional completion values for CS (vertical
axis) compared to other areas of study (horizontal axis). National trends can
be further observed using the overlaying, nested regions generated using
Kernel Density Estimation (KDE). Figure 4. A guide to assist in interpreting
the plots in Figure 3: Cartesian plane restricted to $[-1,1]\times[-1,1]$ with
$r_{CS}$ as the vertical axis and $r_{other}$ horizontal axis.

###  3.2. Lagging Progress

National Trends. Figure 1(b) presents strong evidence of the large gap between
CS and other areas of study when it comes to participation and success of
students from underrepresented groups. For most underrepresented groups, these
gaps are only becoming wider with time. This should be alarming to the CS
community.

In 2011, women accounted for only 19% of the graduating class in CS compared
to other areas of study (58%; Figure 1b). While representation of women in the
graduating class remained relatively stable at 59% in non-CS fields in 2019,
women accounted for only 21% of the graduating class in CS. For black students
in CS, the trends are even more bleak. While CS has much work to do with
female student representation, the trend in CS has greater positive strength
than other areas of study ($r_{CS\wedge Women}=0.91>0.85=r_{other\wedge
Women}$ summarized in Table 1). For black students, completions in CS
_regressed_ sharply from 10% to 8% of the graduating class. During the same
period, the share of Black completions in other areas of study remained
relatively stable at around 9% of the graduating class. For reference, 14.2%
of the total US population is Black or African American alone, as per the US
Census Bureau (Bureau, 2020).

For Hispanic or Latino students, trends in both CS and other areas of study
are strongly positive, indicating good progress. However, the positive trend
is stronger in other areas of study as compared to in CS
$r_{other}=1.0>0.97=r_{CS}$ (Table 1). This indicates that CS is lagging
behind and not able to keep up with other disciplines in its progress for
Hispanic or Latino student representation and success. If these trends (albeit
positive) persist, the gap between CS and other areas of study would only
widen.

Institutional Trends. At the institutional level, we compare completion trends
in CS to other areas of study directly for each institution by taking the
difference of the two completion trend measures: i.e., $r_{CS}-r_{other}$.
Figure 2 depicts the distributions of $r_{CS}-r_{other}$ for each institution
and demographic group. The blue area of each distribution represents
institutions where CS completion trends are lagging behind completion trends
in other areas of study: $r_{CS}<r_{other}$. The orange part of each
distribution represents institutions where CS trends are better than other
areas of study: $r_{CS}>r_{other}$. Table 1 summarizes the distributions in
Figure 2.

For Women, we observe in Figure 2 that the frequency distribution is a bell-
shaped curve with the peak just to the right of zero. The mean difference
(dashed line) for Women is $0.07$ indicating that for an average institution,
the completion trends for Women is marginally better in CS than other areas of
study. For Black students, the frequency distribution is similarly bell-
shaped, but with a peak to the left of zero. However, the mean difference
(dashed line) for Black students is $-0.06$ indicating that for an average
institution, CS is lagging behind other disciplines with respect to completion
trends of Black students.

For Hispanic or Latino students, the frequency distribution is also bell-
shaped but the peak is far to the left of zero with mean difference of
$-0.45$. This suggests that for an average institution in the United States,
CS is lagging far behind other areas of study in terms of completions by
Hispanic or Latino students. We can also observe this result in the rightmost
panel of Figure 2 in which most of the distribution is blue representing the
fact that $r_{CS}<r_{other}$ for an overwhelming majority of institutions.

###  3.3. Simpson’s Paradox

In subsections 3.1 and 3.2, we compared completion trends in CS to other areas
of study for Women, Black and Hispanic or Latino students. Our analysis
considered data from a nationwide scale (dotted line in Figure 2) as well as
at an institutional scale (dashed line in Figure 2). In some cases, we observe
a congruence at both scales and for others, an alarming incongruence.

For Women, the trends at national and institutional scales are congruent: at
the national level $r_{CS}-r_{other}=0.06$ whereas, for an average
institution, $\mu_{r_{CS}-r_{other}}=0.07$, a marginal difference. However,
for Black and Hispanic or Latino students, the trends observed at national
level are incongruous with trends observed at institutional level. We
informally observe this incongruity by noting the wide gaps between dotted and
dashed lines in the center and right panels of Figure 2. Numerically, these
gaps are also reflected in divergent national and institutional values in
columns 3 and 4 of Table 1. These data are a clear example of _Simpson’s
Paradox_ : when a trend appears in several groups of data but disappears or
reverses when the groups are combined.

It is of particular interest that the effect of this paradox on CS completion
trends is different for Black students than for Hispanic or Latino students.
For Black students in CS, the national trend is strongly negative
($r_{CS}-r_{other}=-0.87$) but at the institutional level, the average trend
is relatively neutral ($\mu_{~{}r_{CS}-r_{other}}=-0.06$). In contrast, for
Hispanic or Latino students, the completion trend observed at the national
level is relatively flat ($r_{CS}-r_{other}=-0.03$) whereas the completion
trend observed at the institutional level is strongly negative
($\mu_{~{}r_{CS}-r_{other}}=-0.45$).

Unpacking Simpson’s Paradox. To better understand these observed discrepancies
between national and institutional trends, we attempt to graphically blend
national- and institutional-level data into a single plot for each
underrepresented group. Figure 3 depicts institutional distributions as a
scatterplot with respect to completion trends in CS, $r_{CS}$, and completion
trends in other areas of study, $r_{other}$. We simultaneously observe
national trends in each of plots in Figure 3 by considering densities
identified with the curved regions computed using Kernel Density Estimation
(KDE). To aid the reader, Figure 4 is a legend for how to interpret each plot
in Figure 3. Our plots are restricted to $[-1,1]\times[-1,1]$ of the Cartesian
plane with $r_{CS}$ as the vertical axis and $r_{other}$ horizontal axis
crossing at $(0,0)$. Quadrants are labeled with dotted boxes counterclockwise
$1$ through $4$ consistent with mathematical nomenclature. The colored arrows
indicate an upward (orange) or downward (blue) completion trend for a
particular institution. For example, quadrant $2$ indicates $r_{CS}<0$ and
$r_{other}>0$: the trend of CS completions at an institution is increasing
while completions for other areas of study is decreasing.

For all three demographic groups in Figure 3, the majority of the institutions
fall in the quadrant $1$ (top right) where $r_{CS}>0$ and $r_{other}>0$. The
KDE distribution for Women is nearly uniformly distributed with peaks in the
top-right of quadrant $1$. This aligns with the mean difference of $0.07$ for
Women students observed in Figure 2. The region of the plot corresponding to
$r_{CS}<-0.5$ is relatively sparse indicating that with respect to CS, very
few institutions exhibit strong negative trends. This is consistent with
Figure 1 where completions in other areas of study have stabilized at the
saturation point of around $60\%$ but CS still has a lot of ground to cover.

For Black students (center panel), Figure 3 shows a near bimodal distribution
with most institutions in either quadrant $1$ or quadrant $3$, with quadrants
$2$ and $4$ being relatively sparse. This indicates that institutions with a
positive trend in CS completions also exhibit a positive trend in other areas
of study and vice versa. Reconciling the institutional-level scatter-KDE plot
with the national level trends we observed in Figure 1(b) suggests that
_institutions with $r_{CS}<0$ for Black students must have a
disproportionately large influence to cause a sharply declining trend at the
national level_. We consider this further in subsection 3.4.

For Hispanic or Latino students, the KDE plot in Figure 3 (right panel) is far
more concentrated than the distributions for Women and Black students. This
grouping of institutions in quadrants $1$ and $4$ means the variance along the
$r_{CS}$ vertical axis is far greater than the variance along the $r_{other}$
horizontal axis. Thus, for Hispanic or Latino students, an overwhelming
majority of institutions have a strong positive trend for completions in other
areas of study while the trends of CS completions varies relatively more from
institution to institution. Reconciling the institutional-level scatter-KDE
plot with the national level trends we observed in Figure 2 suggests that
_institutions with $r_{CS}<0$ for Hispanic or Latino students must also have a
disproportionately large influence to cause a positive trend at the national
level but a neutral trend at the institutional level_.

Figure 5. Institutions divided into 4 $r$-categories based on completion
trends: values of $r_{CS}$ and $r_{other}$. The top panel reflects the number
of institutions in each $r$-category for each underrepresented, demographic
group. The bottom panel depicts the average number of CS completions for
institutions in each $r$-category broken out for each underrepresented,
demographic group.

Resolving Simpson’s Paradox. We attempt to resolve the question of why trends
observed at the national level are incongruous with those observed at the
institutional level for Black and Hispanic or Latino students. Figure 5
divides institutions in our dataset into the 4 $r$-categories corresponding to
the quadrants defined in Figure 4. The top panel of Figure 5 shows the
frequency distribution of institutions and the bottom panel is a frequency
distribution with respect to size of institution (where size is based on total
CS completions).

In the top panel of Figure 5 we observe that for every underrepresented
demographic group, most institutions fall within quadrant $1$: $r_{CS}>0$ and
$r_{other}>0$. It is therefore the case that most institutions exhibit a
positive trend in both CS completions and completions in other areas of study.
However, the bottom panel of Figure 5 shows that for Black students, the
institutions where $r_{CS}<0$ (i.e., quadrants $3$ and $4$) also have the
largest graduating classes. This explains why the national-level trends are so
different from the institutional-level trends for Black students: _the few
institutions that have relatively large CS graduating classes are the ones
that exhibit a negative trend in CS completions_. This contrasts with most
institutions that have relatively smaller graduating classes, but exhibit a
positive trend in CS completions. The results are similar for Hispanic or
Latino students in Figure 5 where $r_{CS}<0$ category institutions have a
large number of total completions (graduates). However, because
$r_{CS}>0\wedge r_{other}>0$ category institutions have the largest size for
Hispanic or Latino students, the institution level trends are neutral or flat
as opposed to negative for Black students.

Figure 6. Spatial distribution over the United States of completion trends in
CS and other areas of study from 2011 to 2019 for Women, Black, and Hispanic
or Latino students. Each hexagon implies a geographic region and does _not_
correspond to a state.

###  3.4. Geographic Trends

Thus far, our discussion has focused on completion trends in CS and other
areas of study that are strictly temporal. We augment our analyses with the
added dimension of space to identify any spatiotemporal trends in completions
for underrepresented students. Figure 6 shows hexagon plots over the United
States for each demographic group color-coded by completion trends in CS and
other areas of study. White hexagons correspond to areas of the United States
with either no institution or institutions with no CS graduates belonging to
the given demographic group. The color of each hexagon is determined by the
mode of $r$-category of institutions within the region. The four categories
along with their corresponding color are given in the color bars to the right
and are consistent with the $r$-categories colors from Figure 5.

We observe that completion trends are spatially localized. For Women, negative
trends in CS and other areas of study are observed away from the East and West
coasts with a negative cluster near the Central West region. For Black
students, the negative completion trends in CS and other areas of study seem
to strongly concentrate in the South. In particular, we observe a large
cluster of negative completion trends in the Southeast. These negative trends
are alarming since the majority of the Black population in the United States
lives in these regions, as per the US Census Bureau (Bureau, 2020). A strong
but relatively smaller negative cluster can also be observed near Southern
California and adjacent areas. For Hispanic or Latino, the US census (Bureau,
2020) shows that the population is concentrated in the Southwest but in
contrast to Black students we do not observe any regions of concentrated
negative trends.

##  4\. Conclusions and Future Work

We draw the following conclusions from the analyses presented in this work:

  1. (1)

For Women in Computer Science, the completion rates are far behind when
compared to other areas of study. However, they are trending in the positive
direction.

  2. (2)

For Black students in Computer Science, we observe a sharp decline in
completion trends. This is particularly alarming because a similar decline was
not observed for black students in other areas of study. The negative trend in
completions for Black students at the national level occurs despite the fact
that most institutions individually exhibit positive trends. However, a few
institutions with large graduating classes in CS exhibit strong negative
trends and are driving most of the national level results. Furthermore, the
negative trends for Black students are concentrated generally in the South and
particularly in the Southeast of the United States.

  3. (3)

For trends with Hispanic or Latino students, completions in Computer Science
and other areas of study are both strongly positive. However, the positive
trend in CS is still lagging behind other areas of study. If this trend
continues, the gap between Computer Science and other areas of study will only
increase.

Taken together, the results presented in this paper serve as a call to action
for the CS Education community to include a more data-driven approach to the
problem of Diversity, Equity and Inclusion in CS. There are many potential
avenues for future work that expand on the results presented. For example,
while a completion implies that a student is able to overcome the challenges
of a computing program, it is still just one measure of success at the end of
a long and complex educational pipeline. Similarly, larger institutions, that
seem to have an out-sized influence on national trends in our results, are
more likely to be public with low tuition fees and high student-to-faculty
ratio.

In future work, we intend to include more nuanced categories of race and
ethnicity and American Indian and Alaska Natives, Native Hawaiian or Other
Pacific Islanders and Mixed-race students. While gender and sexual orientation
are important elements of diversity, these are unable to be studied in the
same way due to limitations in the way IPEDS data has been collected. The
IPEDS data is nonetheless a rich source of information that can be used to
better understand the low participation and success of underrepresented
students in CS. Our work only scratches the surface.

## References

  * (1)
  * Bureau (2020) US Census Bureau. 2020\.  Summary File 1.  (2020).
  * Denner and Campe (2023) Jill Denner and Shannon Campe. 2023.  Equity and Inclusion in Computer Science Education: Research on Challenges and Opportunities.  _Computer Science Education: Perspectives on Teaching and Learning in School_ (2023), 85\.
  * FAQ (2023) IPEDS FAQ. 2023\.  https://surveys.nces.ed.gov/ipeds/public/survey-materials/faq?faqid=11
  * Ibe et al. (2018) Nwannediya Ada Ibe, Rebecca Howsmon, Lauren Penney, Nathaniel Granor, Leigh Ann DeLyser, and Kevin Wang. 2018.  Reflections of a diversity, equity, and inclusion working group based on data from a national CS education program. In _Proceedings of the 49th ACM Technical Symposium on Computer Science Education_. 711–716.
  * Kulkarni et al. (2018) Anagha Kulkarni, Ilmi Yoon, Pleuni S Pennings, Kazunori Okada, and Carmen Domingo. 2018.  Promoting diversity in computing. In _Proceedings of the 23rd Annual ACM Conference on Innovation and Technology in Computer Science Education_. 236–241.
  * Larsen and Stubbs (2005) Elizabeth A Larsen and Margaret L Stubbs. 2005.  Increasing diversity in computer science: Acknowledging, yet moving beyond, gender.  _Journal of women and minorities in science and engineering_ 11, 2 (2005).
  * Sax et al. (2017) Linda J Sax, Kathleen J Lehman, Jerry A Jacobs, M Allison Kanny, Gloria Lim, Laura Monje-Paulson, and Hilary B Zimmerman. 2017.  Anatomy of an enduring gender gap: The evolution of women’s participation in computer science.  _The Journal of Higher Education_ 88, 2 (2017), 258–293.
  * Statistics (2012) NC Statistics. 2012\.  Integrated Postsecondary Education Data System.  (2012).  https://nces.ed.gov/ipeds/.


11institutetext: Departamento de Ciencias Físicas, Facultad de Ciencias
Exactas, Universidad Andres Bello, Fernández Concha 700, Las Condes, Santiago,
Chile 22institutetext: Vatican Observatory, Vatican City State, V-00120, Italy
33institutetext: Instituto de Astronomía, Universidad Católica del Norte, Av.
Angamos 0610, Antofagasta, Chile 44institutetext: Centro de Astronomía
(CITEVA), Universidad de Antofagasta, Av. Angamos 601, Antofagasta, Chile
55institutetext: Millennium Institute of Astrophysics , Nuncio Monseñor Sotero
Sanz 100, Of. 104, Providencia, Santiago, Chile 66institutetext: Max-Planck
Institute for Astronomy, Koenigstuhl 17, 69177 Heidelberg, Germany
77institutetext: North Optics Instrumentos Científicos, La Serena, Chile

# A new low-luminosity globular cluster discovered in the Milky Way with the
VVVX survey

E. R. Garro 11 D. Minniti 1122 M. Gómez 22 J. G. Fernández-Trincado 33 J.
Alonso-García 4455 M. Hempel 1166 R. Zelada Bacigalupo 77

(Received: February 17, 2022; Accepted: May 2, 2022)

###### Abstract

Context. The VISTA Variables in the Via Láctea Extended Survey (VVVX) allows
us to probe previously unexplored regions of the inner Milky Way (MW),
especially those that are affected by stellar crowding and strong extinction.

Aims. Our long-term goal is to identify new star clusters and investigate them
to reveal their true nature. In particular, we are looking for new candidate
globular clusters (GCs) located in the Galactic bulge and disk, with the aim
of completing the census of the MW GC system.

Methods. We searched and characterised new GCs, using a combination of the
near-infrared (IR) VVVX survey and Two Micron All Sky survey (2MASS) datasets,
and the optical Gaia Early Data Release 3 (EDR3) photometry and its precise
proper motions (PMs).

Results. We report the discovery of a new Galactic GC, named Garro 02,
situated at RA = 18:05:51.1, Dec = -17:42:02 and $l=12^{\circ}.042$,
$b=+1^{\circ}.656$. Performing a PM-decontamination procedure, we built a
final catalogue with all cluster members, on which we performed a photometric
analysis. We calculated a reddening of $E(J-K_{s})=1.07\pm 0.06$ mag and
extinction of $A_{Ks}=0.79\pm 0.04$ mag in the near-IR; while $E(BP-
RP)=2.40\pm 0.01$ mag and $A_{G}=4.80\pm 0.02$ mag in optical passbands. Its
heliocentric distance is $D=5.6\pm 0.8$ kpc, which places Garro 02 at a
Galactocentric distance of $R_{G}=2.9$ kpc and $Z=0.006$ kpc above the
Galactic plane. We also estimated the metallicity and age by comparison with
known GCs and by fitting PARSEC isochrones, finding [Fe/H] $=-1.30\pm 0.2$ dex
and age $=12\pm 2$ Gyr. We derived the mean cluster PM of
$(\mu_{\alpha}^{\ast},\mu_{\delta})=(-6.07\pm 0.62,-6.15\pm 0.75)$ mas yr-1.
We calculated the cluster luminosity in the near-IR of $M_{Ks}=-7.52\pm 1.23$
mag, which is equivalent to $M_{V}=-5.44$ mag. The core and tidal radii from
the radial density profile are $r_{c}=1.25\pm 0.27$ arcmin (2.07 pc) and
$r_{t}=7.13\pm 3.83$ arcmin (11.82 pc), respectively.

Conclusions. We confirm Garro 02 as a new genuine Galactic GC, located in the
MW bulge. It is a low-luminosity, metal-poor, and old GC, and it is a lucky
survivor of the strong dynamical processes that occurred during the MW’s
entire life.

###### Key Words.:

Galaxy: bulge – Galaxy: stellar content – Stars Clusters: globular – Infrared:
stars – Surveys

##  1 Introduction

Recent studies (e.g. Baumgardt et al. 2019) have demonstrated that star
clusters with lower masses and lower luminosities than those of known Galactic
globular clusters (GCs) may be ubiquitous. We expect that these objects have
been dissolved due to past strong dynamical processes, such as tidal shocks,
evaporations, and galaxy mergers (e.g. Kruijssen et al. 2011), but also they
may be hidden in the Galactic bulge and disk (Minniti et al. 2017; Garro et
al. 2020). Therefore, their identification seems to be very complicated.
However, in the last decade, many new GC candidates have been discovered
across the Milky Way (MW), especially in the bulge regions where the
interstellar dust strongly hides the baryonic matter over the sky (e.g.
Gonzalez et al. 2012; Barbuy et al. 2018) and the stellar crowding is
substantially high. Both the VISTA variables in the Vía Láctea (VVV; Minniti
et al. 2010; Saito et al. 2012) and its eXtension (VVVX; Minniti 2018) have
substantially increased the number of stellar cluster candidates in the MW.
These studies use either visual inspection or photometric analysis to search
for henceforth unknown stellar clusters and they have been extremely
successful in that respect (e.g. Bica et al. 2019; Minniti et al. 2021b).
Although the number of star cluster candidates has increased substantially,
more work has to be done to understand their nature, for example to classify
them as GCs (as well as open clusters or star associations). The Gaia mission
(e.g. Cantat-Gaudin et al. 2019; Castro-Ginard et al. 2020) is extremely
helpful in that endeavour and has in fact confirmed many of the recently
discovered star clusters as GCs, some examples include the following: Garro 01
(Garro et al. 2020), FSR 1758 (Barbá et al. 2019; Villanova et al. 2019;
Romero-Colmenares et al. 2021), Patchick 99 (Garro et al. 2021c), Patchick 125
(Fernández-Trincado et al. 2022; Garro et al. 2021a, named also Gran 3 by Gran
et al. 2022), and VVV-GC160 (Minniti et al. 2021a; Garro et al. 2021d);
another nine GCs confirmed by Garro et al. (2021d), FSR 19 and FSR 25 (Obasi
et al. 2021), and Minni48 (Minniti et al. 2021c); and another five GCs
confirmed in Garro et al. (2021a). These studies as well as other ones are
fundamental since they can contribute, from a broad point of view, for one to
understand the formation and evolution of the MW GC system.

In this paper, we focus on the discovery and confirmation of a new GC
candidate in the MW (see Section 2), which we named Garro 02. In Section 3,
the datasets used to perform our analysis are presented. The decontamination
procedure and the selection of star cluster members are explained in Section
4. In Section 5, we derive its main physical parameters, such as reddening,
extinction, mean cluster proper motions (PMs), distance modulus, heliocentric
and Galactocentric distances, age, metallicity, total luminosity, and finally
structural parameters. A final summary and conclusions are given in Section 6.

##  2 Discovery of a new GC candidate

We report the discovery and confirmation of Garro 02, a new GC located in the
Galactic bulge at equatorial coordinates RA = 18:05:51.1 and Dec = -17:42:02
and Galactic coordinates l = 12.042∘ and b = +1.656∘. Following the same
procedure described in Minniti et al. (2017), we found this cluster using the
VVVX database as a clear over-density of red giant branch stars (RGBs), above
the background, in the tile $e0969$, as shown in Fig. 1. We also observed this
cluster in the optical wavelengths, using the Gaia EDR3 dataset. Also in this
case, we can see a clear excess of stars above the background, as the Gaia
EDR3 source density map (Fig. 1) is clearly shown in a nearly circular region.

However, as demonstrated by previous works, such as Gran et al. (2019), Palma
et al. (2019), and Minniti et al. (2021d), not all over-densities are real
clusters, but they may be simply a grouping of stars or statistical
fluctuations of the projected stellar density in the plane of the sky.
Therefore, we performed four independent tests to confirm Garro 02 as a
genuine GC. First, we visually analysed this cluster at different wavelengths,
using the Panoramic Survey Telescope and Rapid Response System (Pan-STARRS,
Chambers et al. 2016), the Two Micron All Sky Survey (2MASS, Skrutskie et al.
2006), and the Wide-field Infrared Survey Explorer (WISE, Wright et al. 2010)
images of the sky regions centred on Garro 02 as shown in Fig. 2. Both the
near-infrared (IR) 2MASS and mid-IR WISE surveys are sensitive to the
clusters’ bright RGBs; especially in the 2MASS image, a clear excess of RGB
stars is visible at the Garro 02 coordinates, suggesting the presence of a
cluster.

As a second test, we compared the spatial density distribution (Fig. 3) of our
cluster, including stars within $r<4^{\prime}$, with those of a sample of
field stars selected at $5^{\prime}\leq r\leq 8^{\prime}$ from the cluster
centre. As explained in our previous works (Garro et al. 2021a, d), we applied
the Gaussian kernel density estimate (KDE; e.g. Rosenblatt 1956; Parzen 1962)
in order to better distinguish over-densities from lower densities than by
using a visual inspection. As expected, we find that the spatial distributions
of the cluster and the field are very different since a central over-density
is clear for the star sample centred on Garro 02 coordinates, which decreases
towards larger radii. Instead, no central over-density is visible for the
field sample, but a lumpy and non-homogeneous distribution is seen for it.
However, from the contour figures, the maximum number density of the cluster
region is about six times the minimum level, which is not very different from
the referenced field. Therefore, we explored the possibility of having
statistical fluctuations by performing a third test. We constructed the radial
density profile (RDP) for both cluster and stellar field samples. As expected,
we find a flat radial density distribution for a stellar field, whereas a King
(1962) profile for the cluster sample. We postpone a detailed explanation of
this test to Sections 4 and 5.4.

A fourth test regards both the comparison of vector proper motion (VPM)
diagrams and also colour-magnitude diagrams (CMDs) for the cluster and star
field samples. The former gives us information about the dynamics, whereas the
latter allows the photometry of these two samples to be studied. Tests allow
the main differences between the cluster and field samples to be shown. We
postpone a detailed explanation of this test to Section 4.

Finally, these tests support the existence of an over-density of RGB stars,
which could be associated with a new GC candidate. Even so, a spectroscopic
follow-up will confirm the real nature of Garro 02. We performed a photometric
analysis, as explained in Sections 4 and 5.

Figure 1: Density maps to search for the cluster as an over-density. Left
panel: VVVX density map for the entire tile $e0969$ in equatorial coordinates.
The broad stripes seen in this image of the whole tile is the residual of the
combination of the different pawprints. Middle and right panels: Density maps
of the zoomed region around Garro 02 using near-IR VVVX data and optical Gaia
EDR3 data, respectively. We added a black circle and arrow in each panel in
order to highlight the over-density, which indicates the presence of the
cluster.  Figure 2: Pan-STARRS (a), 2MASS (b), VVVX (c), and WISE (d) images
centred on Garro 02. These fields of view are about $21’\times 14’$, oriented
along Galactic coordinates, with the longitude increasing to the left and the
latitude increasing to the top. We show the Pan-STARRS and WISE images for
comparison just to illustrate that the cluster is present at these
wavelengths. The over-density is certainly better seen in the near-IR images
(which is what we study in this paper), due to extinction and incompleteness
of RGB stars in the other bands.

Figure 3: Spatial distribution for Garro 02 selected at $r\lesssim 4^{\prime}$
(on the top). We also compare the cluster spatial distribution with its
relative field (on the bottom) selected at $5^{\prime}\lesssim r\lesssim
8^{\prime}$ from the cluster centre. We used the KDE technique in order to
better distinguish the over-densities (yellow and green colours) from lower
densities (blue colours). All samples were cleaned by nearby stars. These
panels show the following from left to right: the 2D histograms and the
calculated Gaussian KDE on the which we overplotted iso-density contours.

##  3 Datasets

We studied Garro 02 using cross-matched near-IR and optical datasets, adopting
a $0.5^{\prime\prime}$ matching radius. In particular, we used the near-IR
dataset from the VVVX survey, acquired with the VISTA InfraRed CAMera (VIRCAM)
at the 4.1m wide-field Visible and Infrared Survey Telescope for Astronomy
(VISTA; Emerson & Sutherland 2010). The point spread function (PSF) photometry
was extracted as described by Alonso-García et al. (2018). The astrometry was
calibrated to the Gaia Data Release 2 (DR2; Gaia Collaboration et al. 2018b)
astrometric reference system, whereas the photometry was calibrated into the
VISTA magnitude system (González-Fernández et al. 2018) against the 2MASS
using a globally optimised model of frame-by-frame zero points plus an
illumination correction. The advantage of using this public survey is due to
the fact that Garro 02 is located in an obscured region, which may be affected
by differential reddening and high stellar crowding. Therefore, VVVX can help
to see through the interstellar dust since it becomes more transparent at
these wavelengths. Moreover, VVVX is capable of resolving individual stars
better than other surveys, such as 2MASS (Hajdu et al. 2020). Additionally,
the VVVX passbands ($JHK_{s}$) allow one to study evolved stars in the CMDs,
such as red clump (RC) stars, which can be used to derive reliable distance
estimates.

We also included data from 2MASS, a near-IR survey centred on the $JHK$ bands.
We utilized these data to build more complete CMDs since the brightest stars
($K_{s}<11$ mag) in the VVV images are saturated (Saito et al. 2012). We
treated the 2MASS and VVVX data separately. Therefore, we transformed the
2MASS photometry into the VISTA magnitude scale (González-Fernández et al.
2018), since the magnitude scale is different in the two photometric systems,
noting very small differences of $\Delta JK_{s}<0.03$ mag.

We also detected this cluster in the Gaia Early Data Release 3 (EDR3; Gaia
Collaboration et al. 2021) and in the Pan-STARRS datasets (Figs. 1 and 2).
However, we notice that the sampling of the area around Garro 02 is more
homogeneous for the Gaia dataset than for Pan-STARRS. Indeed, we find that a
greater number of stars are resolved in the Gaia survey; they also have more
precise and reliable parameters (astrometry, PMs, and magnitudes). Therefore,
we preferred using only the Gaia EDR3 dataset for our analysis in order to
benefit from its precise astrometry and exquisite PMs. We excluded all nearby
foreground stars with parallax $Plx>0.5$ mas, which is equivalent to $D<3$ kpc
(Bailer-Jones et al. 2018). We also used Gaia EDR3 PMs in order to perform the
PM-decontamination procedure (Section 4). Additionally, we made a cut using
those stars with renormalised unit weight error ruwe¡1.4 (Fabricius et al.
2021, Gaia Collaboration et al. 2021, and references therein), guaranteeing
the high quality of the astrometric solutions. No photometric colour or
magnitude cuts were applied to the Gaia data.

##  4 Decontamination procedure and cluster membership

When studying the inner regions of the MW, various aspects have to be taken
into consideration. First, the differential reddening and extinction, which
not only affects the photometric distance estimates, but also the derived
cluster age and metallicity, especially when isochrones are fitted. Also,
foreground and background field contamination can affect final results because
it can both alter the total luminosity estimate and increase the errors in the
cluster size derivation. However, we have used different procedures to
overcome those problems and reduce uncertainties.

First, we performed a decontamination procedure which allowed us to build
clean CMDs and construct final catalogues with probable star cluster members.
The decontamination of this dataset follows the one described by Garro et al.
(2020, 2021a, 2021c, 2021d). As mentioned in Section 3, we excluded all nearby
stars, through a cut in parallax. Thereafter, we included all stars within a
given radius from the cluster centre. We evaluated the Garro 02 dimension,
using both a visual inspection of the density diagram as a function of the sky
position and applying the Gaussian KDE, as shown in Fig. 3. Therefore, we used
the density map and the KDE iso-density contours in order to identify the more
suitable cluster size, selecting the higher density area around the Garro 02
centre we adopted a $r<5.4^{\prime}$. Finally, we preferred using the PM-
decontamination procedure for two reasons: (i) PMs ensure the stars in the
cluster are divided better from the stars in the field since a cluster is
expected to share coeval motions; (ii) as demonstrated by Garro et al. (2021c)
and Minniti et al. (2021c), this method substantially reduces the percentage
of contaminants better than the statistical decontamination procedure,
excluding stars with disk-like PMs. For this purpose, we constructed the VPM
diagrams as a 2D histogram shown in Fig. 4, in order to better identify the
cluster PM peak. Indeed, as we can see in the bottom-left panel of Fig.4, a
highly peaked distribution is concentrated at $\mu_{\alpha}^{\ast}\approx-0.5$
mas yr-1 and $\mu_{\delta}\approx-2.0$ mas yr-1; while a detached and lower
peaked distribution is located at $\mu_{\alpha}^{\ast}=-6.07\pm 0.62$ mas yr-1
and $\mu_{\delta}=-6.15\pm 0.75$ mas yr-1. The first is related to bulge+disk
field stars, while the second one indicates the presence of the cluster. We
derived the mean cluster PM (and its error) using the $\sigma$-clipping
technique. There is one concern that the PM dispersion is a bit large for a
normal GC, but this may be due to the combination of different effects: the
astrometric error, the residual contamination by field stars, or the possible
ongoing GC disruption. Also, we used this value in order to select stars
within 1 mas yr-1 from the mean cluster PM (see, top panels of Figure 4).
Additionally, we calculated the PM membership probability (Fig. 4), using a
radius $r<1$ mas yr-1 and $r<2$ mas yr-1, in order to exclude as many
contaminants as possible. We can note from Figure 4 that when we adopted $r<2$
mas yr-1, we surely included a higher number of bulge+disk stars, hence we
preferred adopting $r<1$ mas yr-1. This allowed us to select only stars that
are likely members of the cluster and thus construct the final catalogue111The
final catalogue for Garro 02 GC is available in upon request. Please, write to
<EMAIL_ADDRESS>including stars with $>50\%$ probability,
subsequently minimising the likelihood of embedding contaminants. Even so, we
expect that a very small fraction of contaminants may be present in the final
catalogues due to the fact that these regions are characterised by high
stellar crowding and by disk stars with similar PMs, as demonstrated by Garro
et al. (2021c). In fact, as we can see from Figure 5, the two Gaussians
(cluster and field samples) in the VPM diagram overlap. We have reduced the
percentage of contaminants including only stars with a membership probability
$>50\%$; however, we estimate that $<6\%$ of contaminants persist in our
catalogue after the cleaning.

Figure 4: Top. VPM diagram for all star sample (left panel), field sample
(middle panel), and Garro 02 PM members (right panel). In the right panel, the
blue and red circles indicate a radius within 2 and 1 mas yr-1, respectively.
As it is possible to see, we selected the most probable members ($>50\%$)
within 1 mas yr-1 from the mean cluster PMs. Bottom. Left panel: VPM diagram
for Garro 02 cluster as a 2D histogram, displaying both the PMRA (in blue) and
PMDEC (in red) histograms separately. The yellow and light green colours
depict over-densities, whereas the dark green and blue colours represent lower
densities. We indicate with dotted white circles the 1$\sigma$, 2$\sigma$, and
3$\sigma$ centred on the mean cluster PMs. Middle and right panels: VPM
diagrams, where the brightest cloud represents the highest PM probability,
within a radius $r<1$ mas yr-1 and $r<2$ mas yr-1.

Once the final catalogue is obtained, we can build the CMDs for our cluster.
However, before performing the photometric analysis, we carried out a third
and a fourth test, mentioned in Section 2. In this case, we compare the VPM
diagrams and the CMDs for both the cluster (using the decontaminated
catalogue) and the field samples. As shown in Fig. 5, we built the normalised
PM histograms and we fitedt the symmetric Gaussian distributions, which show
two very separated peaks ($\mu$) and very different dispersions ($\sigma$). We
find for Garro 02 cluster: $\mu^{Garro02}_{pmra}=-6.07$ mas yr-1 with
$\sigma^{Garro02}_{pmra}=0.62$ mas yr-1, and $\mu^{Garro02}_{pmdec}=-6.15$ mas
yr-1 with $\sigma^{Garro02}_{pmdec}=0.75$ mas yr-1. Whereas, for the field
sample: $\mu^{field}_{pmra}=-1.57$ mas yr-1 with $\sigma^{field}_{pmra}=2.56$
mas yr-1, and $\mu^{field}_{pmdec}=-3.60$ mas yr-1 with
$\sigma^{field}_{pmdec}=3.07$ mas yr-1.
However, it is notable that a symmetric Gaussian distribution is a rough fit
for the field sample, but this is expected since the field stars show some
structured motions, but their broad and skewed PM distributions are very
different from the cluster that is more concentrated and nearly Gaussian,
which is as expected for a GC.
Finally, we used the CMDs as excellent tools to understand the cluster nature.
In fact, comparing the cluster and field near-IR VVVX CMDs (Fig. 5), we can
see that the blue horizontal branch (HB) is visible in the cluster CMD,
suggesting a low metallicity, and also a narrow red giant branch (RGB)
sequence is evident. On the other hand, no clear sequences are displayed by
the field CMD, but a very elongated blue sequence is shown, indicating an
excess in the disk stars’ main sequence, as well as wide red sequence,
belonging to the RGB stars from the bulge. Unfortunately, fainter stars
belonging to the cluster MS turn off are missing in our compilation since they
just are below the magnitude detection limit.

The last test we carried out was to probe the presence or lack thereof of
statistical fluctuations in the Garro 02 regions is comparing the RDP of Garro
02 with the field sample. We have dedicated a section (see Sect. 5.4), where
we explain how we constructed the RDP for our target. Here, we just show the
RDP for the field sample. We preferred using a sub-sample of stars, selecting
only those within 1.5 mas yr-1 of the mean stellar field PMs, which were
obtained previously, in order to make a sensible comparison between the two
RDPs. Figure 6 shows the RDP for the PM-selected field sample. As expected for
these areas, we find a flat distribution, with a constant value of 5-8 stars
arcmin-2, while as explained in Section 5.4, a King profile represents a good
approximation for the Garro 02 RDP. Finally, these features confirm the
cluster nature of Garro 02.

Figure 5: Top panels: Normalised $\mu_{\alpha}^{\ast}$ (left) and
$\mu_{\delta}$ (right) histograms for the field sample (in red) and Garro 02
cluster (in blue). We fitted the symmetric Gaussian distributions to highlight
the differences in the peaks and in the dispersions of the two distributions.
Bottom panels: Near-IR VVVX CMDs for the field sample (on the left) and for
Garro 02 (on the right). The green-coloured areas are representative of over-
densities, while the blue ones are lower densities.  Figure 6: Flat radial
density profile for the PM-selected field stars. The black points correspond
to the field density profile obtained with no subtraction of the background
contribution.

##  5 Photometric characterisation of Garro 02

Once decontaminated catalogues with probable members are built, we can measure
the main physical parameters that characterise the cluster, such as the
reddening and extinction, the distance modulus and heliocentric distance, the
metallicity and age, the total luminosity, and also the main structural
parameters, such as the tidal radius ($r_{t}$) and core radius ($r_{c}$).

###  5.1 Reddening, extinction, and distance determinations

Although the advantage of using multi-wavelength photometry, as in our case,
is that it can help to reduce the sources of error substantially and reach
robust results, it is often tricky to find all the parameters which fit both
optical and near-IR CMDs simultaneously. For that reason, we derived the
reddenings and extinctions using different methods.

First, we used the reddening maps of Schlafly & Finkbeiner (2011), recovering
an optical extinction of $A_{V}=8.57\pm 0.02$ mag. Therefore, adopting the
known extinction relation $A_{Ks}\leavevmode\nobreak\ =\leavevmode\nobreak\
0.11\leavevmode\nobreak\ \times\leavevmode\nobreak\ A_{V}$ (Schlafly &
Finkbeiner 2011), we find in the near-IR $A_{Ks}=0.94\pm 0.02$ mag. On the
other hand, we used two extinction coefficients for deriving the excess
colour: $R_{Ks}=0.428$ (Alonso-García et al. 2017) was obtained for the inner
bulge and $R_{Ks}=0.75$ (Cardelli et al. 1989) was mostly used for regions in
the outer bulge; they allowed us to calculate an excess colour of
$E(J-Ks)=2.20\pm 0.05$ mag and $E(J-Ks)=1.27\pm 0.03$ mag, respectively. In
the first case, we measured a distance modulus of $(m-M)=14.16\pm 0.05$ mag,
equivalent to $D=6.8\pm 0.2$ kpc, whereas in the second case the distance
modulus is $(m-M)=13.64\pm 0.03$ mag, which is equivalent to a heliocentric
distance of $D=5.4\pm 0.2$ kpc.

Figure 7: Comparison between the near-IR CMDs for Garro 02 (grey dots) and NGC
6642 (grey diamonds). We shifted Garro 02 by 0.8 in $(J-K_{s})$ colour to
avoid overlapping points.

Second, we built near-IR CMDs via comparison with the well-known Galactic GC
NGC 6642, which is located at a distance of $D=8.1\pm 0.2$ kpc with a distance
modulus of $(m-M)=14.53\pm 0.17$ mag (Balbinot et al. 2009; Baumgardt &
Vasiliev 2021). We compared these two clusters (see Fig. 7), because they show
similar RGB and HB morphologies. Clearly, we can see a complete CMD, showing
very narrow and well-defined sequences for NGC 6642. Whereas, for Garro 02, we
can see a wider RGB due to a higher reddening in this field, as well as a
scarcely populated HB due to the fact that Garro 02 is fainter than NGC 6642.
Indeed, we observed that their sequences are separated by
$\Delta(J-K_{s})=0.86$ mag, with Garro 02 being more reddened, and
$\Delta(K_{s})=0.09$, with Garro 02 being slightly fainter. Taking into
account that the NGC 6642 reddening is $E(J-K_{s})=0.21\pm 0.06$ mag, this
yields a reddening $E(J-K_{s})=1.07\pm 0.06$ mag. The differential comparison
with NGC 6642 implies a distance modulus of $(m-M)=13.79\pm 0.06$ mag,
equivalent to $D=5.7\pm 0.7$ kpc.

Independently, we derived reddening and extinction in the optical passbands.
Since clusters with blue HBs usually do not show a defined RC, as Garro 02, we
derived its colour excess, comparing the cluster HB and RGB colours with the
absolute colours by Gaia Collaboration et al. (2018a). We obtained a reddening
of $E(BP-RP)=2.40\pm 0.01$ mag and an extinction $A_{G}=4.80\pm 0.02$ mag,
from which we calculated a distance modulus of $(m-M)_{0}=13.32\pm 0.02$ mag,
which is equivalent to a heliocentric distance of $D=4.6\pm 0.2$ kpc. Knowing
that $A_{G}=10.116\times A_{Ks}$, we derived the near-IR extinction
$A_{Ks}=0.47\pm 0.01$ mag, which allowed us to obtain a distance modulus of
$(m-M)_{0}=14.07\pm 0.01$ mag, equivalent to $D=6.5\pm 0.2$ kpc.
Hence, we find that all these distance estimates are all in agreement within
errors ($<2\sigma$). We adopted a mean value of $D=5.6\pm 0.8$ kpc, as the
cluster distance, thus we placed Garro 02 at a Galactocentric distance of
$R_{G}=2.9$ kpc, assuming $R_{\odot}=8.2$ kpc (Gravity Collaboration et al.
2019) and at a distance above the Galactic plane of $Z=0.006$ kpc, using the
relation $Z=D\times\sin(b)$ and assuming $Z_{\odot}=0$ kpc. Finally, once
distance and mean cluster PMs were evaluated, we could derive the tangential
velocities, finding $V_{T}^{RA}=-163.87\pm 33.62$ km s-1 and
$V_{T}^{Dec}=-166.26\pm 16.15$ km s-1. Therefore, the position and motions
indicate that Garro 02 is a new bulge cluster.

###  5.2 Metallicity and age estimates

In the previous section (Sect. 5.1), we derived the reddening and the
extinction with different methods. Therefore, we used these values in order to
fit a family of isochrones in the near-IR and optical CMDs. We found that the
parameters that allow us to better reproduce the evolutionary sequences in the
CMDs are those reported in Table 1, which were derived using the comparison
with the GC NGC 6642. Once these parameters were fixed, we performed the
isochrone-fitting methods in order to derive the metallicity and age of our
target. From the CMD morphology of this cluster and the similarities with the
GC NGC 6642 CMD, we expect that Garro 02 is a metal-poor and old cluster.
We downloaded the PARSEC+COLIBRI isochrones222http://stev.oapd.inaf.it/cgi-
bin/cmd (Bressan et al. 2012; Marigo et al. 2013). We first generated a family
of isochrones, considering a metallicity range of $-1.0\leqslant$ [Fe/H]
$\leqslant-2.0$ dex, and we superimposed these tracks on the CMDs both in
near-IR and optical passbands (Fig. 8). After that, we followed the same
procedure with ages, generating isochrones with ages between 8 and 13 Gyr.
Hence, comparing the evolutionary sequences in the Garro 02 CMDs with the
isochrones generated for different ages and metallicities, we visually
selected the best fit, yielding a [Fe/H] = $-1.30\pm 0.2$ dex and Age $=12\pm
2$ Gyr. We provide an estimation of the age and metallicity errors by changing
the parameters until the isochrones are not adequate to reproduce all
evolutionary sequences both in Gaia and VVVX+2MASS CMDs. Recovering the
metallicity is simpler than the age since the more evolved sequences, such as
the HB and the RGB, are more affected by variations in metallicity, whereas
deriving a precise age estimate is a tricky task to perform, especially when
stars fainter than the main sequence turn-off (MSTO) are below the detection
limit (see Figure 9). Despite this, the vertical (as well as the horizontal)
method (e.g. Gratton & Ortolani 1987; Sarajedini & King 1989; Marín-Franch et
al. 2010) is largely used for the relative age determination. The estimate of
the magnitude difference between the TO and the HB level is a good relative
age indicator for a cluster. The advantage of using this method is that it is
independent of the reddening, distance modulus, and metallicity. We can see
from Figures 8 and 9 that the difference in magnitude between the HB and MSTO
is $\Delta K_{s}(HB-MSTO)>3$ mag and $\Delta G(HB-MSTO)>2.5$ mag, indicating
that our cluster is older than 10 Gyr. This method helped us to limit the age-
metallicity degeneracy, since we know that ages ¡ 10 Gyr are not suitable for
Garro 02. However, this represents a lower limit since our CMDs do not reach
the MSTO. From our analysis, we can conclude that Garro 02 is a metal-poor and
old bulge GC.

Figure 8: Optical Gaia EDR3 (right panel) and near-IR VVVX+2MASS (left panel)
CMDs for the GC Garro 02. We have included only PM-selected members,
considering also those stars with a membership probability greater than 50%.
The isochrone-fitting method yields an Age = 12 Gyr and [Fe/H] = -1.30 dex. We
have superimposed the PARSEC model represented as red dotted lines in both
panels.

Figure 9: Optical Gaia EDR3 (middle and right panels) and near-IR VVVX+2MASS
(middle and left panels) CMDs for the GC Garro 02. We have included only PM-
selected members, considering also those stars with a membership probability
greater than 50%. We used a family of isochrones, changing metallicity and age
in turn, in order to limit the age-metallicity degeneracy and derive their
errors. For our convenience, we specify here that the unit of the age is ‘Gyr’
and that of the metallicity is ‘dex’ in the legend.

###  5.3 Total luminosity derivation

As mentioned in Section 5.2, MS faint stars are missing in our compilation;
therefore, it is complicated to derive a precise estimate of the total
luminosity for this GC. However, as demonstrated in our previous works (Garro
et al. 2021a, b, c, d), we can derive a good estimation of this parameter,
comparing the luminosity of Garro 02 with those of known and well-
characterised GCs, as long as they have a similar metallicity333We consider
the metallicity as good parameter for comparing GCs because it has an effect,
as seen in theoretical isochrones for example. In general, metal-poor GCs
should be bluer, and therefore brighter, than metal-rich GCs at the same
mass.. We first derived the absolute magnitude, measuring the cluster total
flux. We obtain $M_{Ks}=-7.52\pm 1.23$ mag, which is equivalent to
$M_{V}=-5.24$ mag, if we assume the typical GC colour of $(V-K_{s})=2.5$ (e.g.
Barmby et al. 2000; Conroy & Gunn 2010). Then, we applied the procedure
described by Garro et al. (2021b) in order to quantify the fraction of
luminosity missed. In summary, we first calculated the absolute magnitude in
$K_{s}$-band for the Galactic GCs, using the VVV/VVVX catalogues, thus in the
same magnitude ranges as Garro 02 ($10.2\lesssim K_{s}\lesssim 17.2$ mag).
After that, the obtained $M_{Ks}$ was converted into the absolute magnitude in
V-band $M_{V}$, using again the typical GC colour. Finally, we scaled the
$M_{V}$ values to the absolute magnitude listed in the 2010 version of the
Harris (1996) catalogue, estimating the fraction of luminosity that comes from
the faintest stars for the known GCs. In particular, we compared Garro 02 with
NGC 6642 ([Fe/H] $=-1.26$, Barbuy et al. 2006) and NGC 6626 ([Fe/H] $=-1.32$,
Kerber et al. 2018). Therefore, we calculated the average of $\Delta(M_{V})$
fractions for these two clusters and finally we added this average to the
Garro 02 luminosity. We estimate that this fraction of missed luminosity is
$\sim 0.4$ mag (on average), and therefore fainter stars do not strongly
change the total luminosity. In this way, we have obtained an empirical
correction to our cluster’s luminosity under the assumption of similarity
among Galactic GCs. We obtain a total luminosity in V-band of $M_{V}=-5.44$
mag, thus finding a low luminosity GC, that falls in the fainter tail of the
MW GC luminosity function (Rejkuba 2012). It is notable that the total
luminosity should also be uncertain because the luminosity function of a
cluster may be altered within a strong tidal field. However, as demonstrated
by Kharchenko et al. (2016), we expect that the first 10–12 brighter members
accumulate more than half of the integrated luminosity of a cluster;
therefore, our luminosity estimate is a lower limit of the real total
luminosity of the GC since both optical and near-IR catalogues are incomplete.

###  5.4 Radial density profile

We computed the RDP for Garro 02 with the purpose of determining its true
physical size. The old centre was initially found by visual inspection of the
available near-IR images and density maps, and it should therefore be more
uncertain than the new centre determination. Therefore, we checked that the
centre coordinates of the system were accurately determined for a proper
determination of the RDP. Various papers, presented in the literature, explain
different methods for deriving the cluster centre, some of them are based on
the stellar counts (e.g. Lanzoni et al. 2019) or the position of the maximum
of the surface brightness (Trager et al. 1995). Since our compilation includes
only evolved stars, we used the maximum spatial density determined through a
2D Gaussian KDE density estimator in order to assess the location among our
points (or stars) where the maximum density is found. In this way we can avoid
any possible bias induced by the presence of a few bright stars, which would
alter the position of the surface brightness maximum.
We first divided the distribution of points in a grid both in RA and Dec,
between 0∘ (coincident with the cluster centre) and 0.08∘, with a step of
0.005∘. For estimating the point density using a kernel, we set two functions.
The first one performs two tasks at once: it takes two points defined by RA1,
DEC1, RA2, and DEC2; calculates their Euclidean distance; and uses that
distance to evaluate a Gauss kernel’s function value. The second function
transforms the Gauss kernel to approximately yield $1$ when the distance is
$0$ (or very small), and to have 2$\sigma$ fixed at the previously calculated
radius. After that, the code identifies the maximum value in each grid and
gives us the maximum of the density distribution of our cluster. In this way,
we find that the new centre is $RA=$ 18:05:54.0 (equivalent to
$271.475^{\circ}$) and $Dec=$ -17:42:00 (equivalent to $-17.700^{\circ}$) in
the equatorial system, and $l=12.048^{\circ}$ and $b=1.646^{\circ}$ in the
galactic system, which slightly differ from the old centre of $\Delta
RA=-0.0121^{\circ}$ and $\Delta Dec=-0.0002^{\circ}$.

We used this new centre in order to construct the RDP for our target. We
divided our sample into 14 circular annuli, increasing radii of $0.3^{\prime}$
out to an outer radius of $4.8^{\prime}$. Although our statistics are low
because they are only related to evolved stars, this allowed us to include a
good number of stars in each bin ($N\geqslant 7$ stars). We calculated the
density per bin as the number of stars (N) divided by the area (A). The
errorbars were computed as $e=\sqrt{N}/A$. We also derived the true RDP
correcting our profile subtracting the background level, fixed at 0.5
stars/arcmin2. Thereafter, the RDP for Garro 02 was plotted as function of the
mean distance of the circular annulus to the cluster centre (see Fig. 10).
One of the most used models to fit the cluster density profile is the King
(1962) model. Even so, the true shape should likely be more complex; since we
base our study on evolved stars which may be more concentrated, we find that
the King (1962) model is a good representation of the projected number
density. Also, it allows us to provide the cluster main structural parameters.
The best-fit King model yields a core radius of $r_{c}=1.25\pm 0.27$ arcmin
and a tidal radius of $r_{t}=7.13\pm 3.83$ arcmin, which are equivalent to
$r_{c}=2.04\pm 0.44$ pc and $r_{t}=11.61\pm 6.24$ pc at the cluster distance
of $D=5.6$ kpc. These values are consistent with the typical MW GC radii
listed in the 2010 Harris (1996) compilation. However, owing to the low
density and also no well-resolved profile due to the low statistics especially
at smaller radii, we are not able to suggest if Garro 02 is a pre-, post- or
core-collapse GC. However, calculating the concentration parameter
$c=\log(r_{t}/r_{c})=1.74$ and comparing this value with GCs listed in the
2010 Harris (1996) compilation, we exclude that Garro 02 is a core-collapse
cluster.

Figure 10: Left panel: Determination of the Garro 02 centre using the KDE
technique. As specified by the legend, the blue points are the PM selected
cluster members, while we highlight the old centre with a magenta cross, and
the new centre with a red plus sign. The yellow colour depicts higher
densities, which decreasing to lower densities becomes orange, blue, and
black. Right panel: The black points correspond to the cluster density profile
obtained after subtraction of the background contribution. The red line shows
the best-fit King model profile, with the corresponding values of the core
radius of $r_{c}=1.3^{\prime}$ and a tidal radius of $r_{t}=7.1^{\prime}$.

##  6 Summary and conclusions

We report the discovery of a new GC candidate located in the MW bulge, named
Garro 02. We found this cluster as a clear over-density of red giant stars in
the VVVX density map (Fig. 1), and also the cluster PMs peak differ from the
bulk of the field PMs (Fig. 4). We confirmed the cluster nature performing
four independent tests, making a comparison between the spatial distributions,
radial density distributions, VPM diagrams, and CMDs between the cluster and
the field samples. Thereafter, we carried out a multi-band photometric
analysis, using a combination of three surveys: VVVX and 2MASS in the near-IR,
and Gaia EDR3 in the optical passbands. We measured the cluster’s physical
parameters, including reddening, extinction, distance, metallicity, age, total
luminosity, and also structural parameters deriving core and tidal radii. Our
main results for the parameters of Garro 02 are listed in Table 1.
Garro 02 is a genuine GC located in the Galactic bulge at a heliocentric
distance of $5.6$ kpc and a galactocentric distance of $2.9$ kpc. It is a
metal-poor ([Fe/H] = $-1.30$) and old (age $\sim 12$ Gyr) GC. Additionally, it
is a low-luminosity object with $M_{V}=-5.44$ mag, which is fainter than the
MW GC luminosity function peak ($M_{V}^{peak}=-7.4\pm 0.2$, Harris 1991;
Ashman & Zepf 1998, and in agreement with the value found by Garro et al.
2021b of $M_{V}^{peak}=-7.46\pm 1.04$ mag), suggesting that this cluster may
have survived a strong dynamical process.
Finally, deeper observations need to be done in order to reach the main
sequence to provide a more precise age estimate and complete the Garro 02-star
catalogue, including fainter stars, to confirm its structural parameters and
better calculate the total luminosity for this cluster. Spectroscopic analysis
needs to be performed in order to confirm the metallicity found by the present
work, and also to derive radial velocities and chemical abundances.

Table 1: Main physical parameters for Garro 02 determined using near-IR and optical datasets. Physical parameters | Garro 02
---|---
RA (J2000) | 18h05m54.0s
DEC (J2000) | -17d42m00s
Latitude | $12^{\circ}.048$
Longitude | $+1^{\circ}.646$
$A_{K_{s}}$ [mag] | $0.79\pm 0.04$
$A_{G}$ [mag] | $4.80\pm 0.02$
$E(J-K_{s})$ [mag] | $1.07\pm 0.06$
$E(BP-RP)$ [mag] | $2.40\pm 0.01$
$D_{mean}$ [kpc] | $5.6\pm 0.8$
$R_{G}$ [kpc] | $2.9$
$Z$ [kpc] | $0.006$
$M_{K_{s}}$ [mag] | $-7.52\pm 1.23$
$M_{V}$ [mag] | $-5.44$
[Fe/H] [dex] | $-1.30\pm 0.2$
Age [Gyr] | $12.0\pm 2$
$\mu_{\alpha}^{\ast}$ [mas/yr] | $-6.07\pm 0.62$
$\mu_{\delta}$ [mas/yr] | $-6.15\pm 0.75$
$V_{T}^{RA}$ [km/s] | $-163.87\pm 33.62$
$V_{T}^{Dec}$ [km/s] | $-166.26\pm 16.15$
$r_{c}$ [arcmin] |  $1.25\pm 0.27$ ($2.04$ pc)
$r_{t}$ [arcmin] |  $7.13\pm 3.83$ ($11.61$ pc)

###### Acknowledgements.

The authors are grateful to the anonymous referee for providing helpful
comments and suggestions, which have improved the content of this paper. We
gratefully acknowledge the use of data from the ESO Public Survey program IDs
179.B-2002 and 198.B-2004 taken with the VISTA telescope and data products
from the Cambridge Astronomical Survey Unit. ERG acknowledges support from
ANID PhD scholarship No. 21210330. D.M. gratefully acknowledges support by the
ANID BASAL projects ACE210002 and FB210003, and Fondecyt Project No. 1220724.
This work has made use of data from the European Space Agency (ESA) mission
Gaia (http:/www.cosmos.esa.int/gaia), processed by the Gaia Data Processing
and Analysis Consortium (DPAC,
http://www.cosmos.esa.int/web/gaia/dpac/consortium). J.G.F-T gratefully
acknowledges the grant support provided by Proyecto Fondecyt Iniciación No.
11220340, and also from ANID Concurso de Fomento a la Vinculación
Internacional para Instituciones de Investigación Regionales (Modalidad corta
duración) Proyecto No. FOVI210020, and from the Joint Committee ESO-Government
of Chile 2021 (ORP 023/2021). J.A.-G. acknowledges support from Fondecyt
Regular 1201490 and from ANID – Millennium Science Initiative Program –
ICN12_009 awarded to the Millennium Institute of Astrophysics MAS.

## References

  * Alonso-García et al. (2017) Alonso-García, J., Minniti, D., Catelan, M., et al. 2017, ApJ, 849, L13
  * Alonso-García et al. (2018) Alonso-García, J., Saito, R. K., Hempel, M., et al. 2018, A&A, 619, A4
  * Ashman & Zepf (1998) Ashman, K. M. & Zepf, S. E. 1998, Cambridge Univ. Press, Cambridge (Cambridge Astrophysics Series; 30)
  * Bailer-Jones et al. (2018) Bailer-Jones, C. A. L., Rybizki, J., Fouesneau, M., Mantelet, G., & Andrae, R. 2018, AJ, 156, 58
  * Balbinot et al. (2009) Balbinot, E., Santiago, B. X., Bica, E., & Bonatto, C. 2009, Monthly Notices of the Royal Astronomical Society, 396, 1596
  * Barbá et al. (2019) Barbá, R. H., Minniti, D., Geisler, D., et al. 2019, ApJ, 870, L24
  * Barbuy et al. (2006) Barbuy, B., Bica, E., Ortolani, S., & Bonatto, C. 2006, A&A, 449, 1019
  * Barbuy et al. (2018) Barbuy, B., Chiappini, C., & Gerhard, O. 2018, ARA&A, 56, 223
  * Barmby et al. (2000) Barmby, P., Huchra, J. P., Brodie, J. P., et al. 2000, The Astronomical Journal, 119, 727
  * Baumgardt et al. (2019) Baumgardt, H., Hilker, M., Sollima, A., & Bellini, A. 2019, MNRAS, 482, 5138
  * Baumgardt & Vasiliev (2021) Baumgardt, H. & Vasiliev, E. 2021, MNRAS, 505, 5957
  * Bica et al. (2019) Bica, E., Pavani, D. B., Bonatto, C. J., & Lima, E. F. 2019, AJ, 157, 12
  * Bressan et al. (2012) Bressan, A., Marigo, P., Girardi, L., et al. 2012, MNRAS, 427, 127
  * Cantat-Gaudin et al. (2019) Cantat-Gaudin, T., Krone-Martins, A., Sedaghat, N., et al. 2019, A&A, 624, A126
  * Cardelli et al. (1989) Cardelli, J. A., Clayton, G. C., & Mathis, J. S. 1989, ApJ, 345, 245
  * Castro-Ginard et al. (2020) Castro-Ginard, A., Jordi, C., Luri, X., et al. 2020, A&A, 635, A45
  * Chambers et al. (2016) Chambers, K. C., Magnier, E. A., Metcalfe, N., et al. 2016, arXiv e-prints, arXiv:1612.05560
  * Conroy & Gunn (2010) Conroy, C. & Gunn, J. E. 2010, The Astrophysical Journal, 712, 833
  * Emerson & Sutherland (2010) Emerson, J. & Sutherland, W. 2010, The Messenger, 139, 2
  * Fabricius et al. (2021) Fabricius, C., Luri, X., Arenou, F., et al. 2021, A&A, 649, A5
  * Fernández-Trincado et al. (2022) Fernández-Trincado, J. G., Minniti, D., Garro, E. R., & Villanova, S. 2022, A&A, 657, A84
  * Gaia Collaboration et al. (2018a) Gaia Collaboration, Babusiaux, C., van Leeuwen, F., et al. 2018a, A&A, 616, A10
  * Gaia Collaboration et al. (2018b) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018b, A&A, 616, A1
  * Gaia Collaboration et al. (2021) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2021, A&A, 649, A1
  * Garro et al. (2021a) Garro, E. R., Minniti, D., Alessi, B., et al. 2021a, arXiv e-prints, arXiv:2112.13591
  * Garro et al. (2021b) Garro, E. R., Minniti, D., Gómez, M., & Alonso-García, J. 2021b, A&A, 654, A23
  * Garro et al. (2020) Garro, E. R., Minniti, D., Gómez, M., et al. 2020, A&A, 642, L19
  * Garro et al. (2021c) Garro, E. R., Minniti, D., Gómez, M., et al. 2021c, A&A, 649, A86
  * Garro et al. (2021d) Garro, E. R., Minniti, D., Gómez, M., et al. 2021d, arXiv e-prints, arXiv:2111.08317
  * Gonzalez et al. (2012) Gonzalez, O. A., Rejkuba, M., Zoccali, M., et al. 2012, A&A, 543, A13
  * González-Fernández et al. (2018) González-Fernández, C., Hodgkin, S. T., Irwin, M. J., et al. 2018, MNRAS, 474, 5459
  * Gran et al. (2019) Gran, F., Zoccali, M., Contreras Ramos, R., et al. 2019, A&A, 628, A45
  * Gran et al. (2022) Gran, F., Zoccali, M., Saviane, I., et al. 2022, MNRAS, 509, 4962
  * Gratton & Ortolani (1987) Gratton, R. G. & Ortolani, S. 1987, A&AS, 67, 373
  * Gravity Collaboration et al. (2019) Gravity Collaboration, Abuter, R., Amorim, A., et al. 2019, A&A, 625, L10
  * Hajdu et al. (2020) Hajdu, G., Dékány, I., Catelan, M., & Grebel, E. K. 2020, Experimental Astronomy, 49, 217
  * Harris (1991) Harris, W. E. 1991, ARA&A, 29, 543
  * Harris (1996) Harris, W. E. 1996, The Astronomical Journal, 112, 1487
  * Kerber et al. (2018) Kerber, L. O., Nardiello, D., Ortolani, S., et al. 2018, The Astrophysical Journal, 853, 15
  * Kharchenko et al. (2016) Kharchenko, N. V., Piskunov, A. E., Schilbach, E., Röser, S., & Scholz, R. D. 2016, A&A, 585, A101
  * King (1962) King, I. 1962, AJ, 67, 471
  * Kruijssen et al. (2011) Kruijssen, J. M. D., Pelupessy, F. I., Lamers, H. J. G. L. M., Portegies Zwart, S. F., & Icke, V. 2011, MNRAS, 414, 1339
  * Lanzoni et al. (2019) Lanzoni, B., Ferraro, F. R., Dalessandro, E., et al. 2019, ApJ, 887, 176
  * Marigo et al. (2013) Marigo, P., Bressan, A., Nanni, A., Girardi, L., & Pumo, M. L. 2013, MNRAS, 434, 488
  * Marín-Franch et al. (2010) Marín-Franch, A., Cassisi, S., Aparicio, A., & Pietrinferni, A. 2010, ApJ, 714, 1072
  * Minniti (2018) Minniti, D. 2018, in Astrophysics and Space Science Proceedings, Vol. 51, The Vatican Observatory, Castel Gandolfo: 80th Anniversary Celebration, ed. G. Gionti & J.-B. Kikwaya Eluo, 63
  * Minniti et al. (2021a) Minniti, D., Fernández-Trincado, J. G., Gómez, M., et al. 2021a, A&A, 650, L11
  * Minniti et al. (2017) Minniti, D., Geisler, D., Alonso-García, J., et al. 2017, The Astrophysical Journal, 849, L24
  * Minniti et al. (2021b) Minniti, D., Gómez, M., Alonso-García, J., Saito, R. K., & Garro, E. R. 2021b, A&A, 650, L12
  * Minniti et al. (2010) Minniti, D., Lucas, P. W., Emerson, J. P., et al. 2010, New A, 15, 433
  * Minniti et al. (2021c) Minniti, D., Palma, T., Camargo, D., et al. 2021c, A&A, 652, A129
  * Minniti et al. (2021d) Minniti, D., Palma, T., & Claria, J. J. 2021d, BAAA, 107M
  * Obasi et al. (2021) Obasi, C., Gómez, M., Minniti, D., & Alonso-García, J. 2021, A&A, 654, A39
  * Palma et al. (2019) Palma, T., Minniti, D., Alonso-García, J., et al. 2019, Monthly Notices of the Royal Astronomical Society, 487, 3140
  * Parzen (1962) Parzen, E. 1962, The Annals of Mathematical Statistics, 33, 1065
  * Rejkuba (2012) Rejkuba, M. 2012, Ap&SS, 341, 195
  * Romero-Colmenares et al. (2021) Romero-Colmenares, M., Fernández-Trincado, J. G., Geisler, D., et al. 2021, A&A, 652, A158
  * Rosenblatt (1956) Rosenblatt, M. 1956, The Annals of Mathematical Statistics, 27, 832
  * Saito et al. (2012) Saito, R. K., Hempel, M., Minniti, D., et al. 2012, A&A, 537, A107
  * Sarajedini & King (1989) Sarajedini, A. & King, C. R. 1989, AJ, 98, 1624
  * Schlafly & Finkbeiner (2011) Schlafly, E. F. & Finkbeiner, D. P. 2011, ApJ, 737, 103
  * Skrutskie et al. (2006) Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, 131, 1163
  * Trager et al. (1995) Trager, S. C., King, I. R., & Djorgovski, S. 1995, AJ, 109, 218
  * Villanova et al. (2019) Villanova, S., Monaco, L., Geisler, D., et al. 2019, ApJ, 882, 174
  * Wright et al. (2010) Wright, E. L., Eisenhardt, P. R. M., Mainzer, A. K., et al. 2010, AJ, 140, 1868