BSM physics using photon-photon fusion
processes in UPC in Pb+Pb collisions
with the ATLAS detector

Klaudia Maj111On behalf of the ATLAS Collaboration

AGH University of Kraków

al. Mickiewicza 30, 30-059 Kraków, Poland

e-mail<EMAIL_ADDRESS>
> Relativistic heavy-ion beams at the LHC are accompanied by a large flux of
> equivalent photons, leading to multiple photon-induced processes. This
> proceeding presents searches for physics beyond the Standard Model enabled
> by photon-photon processes in both di-tau and diphoton final states. The
> tau-pair production measurements can constrain the tau lepton’s anomalous
> magnetic dipole moment (g-2), and a recent ATLAS measurement using muonic
> decays of tau leptons in association with electrons and tracks provides one
> of the most stringent limits available to date. Similarly, light-by-light
> scattering proceeds via loop diagrams, which can contain particles not yet
> directly observed. Thus, high statistics measurements of light-by-light
> scattering provide a precise and unique opportunity to investigate
> extensions of the Standard Model, such as the presence of axion-like
> particles.

> PRESENTED AT
>
>
>
>
> DIS2023: XXX International Workshop on Deep-Inelastic Scattering and Related
> Subjects,
> Michigan State University, USA, 27-31 March 2023
>

##  1 Introduction

In recent years, photon-induced processes in heavy-ion (HI) collisions have
emerged as a promising path for studying beyond the Standard Model (BSM)
physics. The ATLAS experiment [1] at the Large Hadron Collider (LHC) dedicates
part of its annual operational time to the HI physics including ultra-
peripheral collisions (UPC). UPC are a unique category of HI collisions, which
occur when the distance separating the interacting nuclei exceeds the sum of
their radii. The large electromagnetic fields generated by relativistic ions
can be considered as fluxes of photons - as described in the Equivalent Photon
Approximation formalism [2, 3]. Photon-photon interactions occur in both
proton-proton, $pp$, and HI collisions. However, in the latter, the cross-
sections for a specific process experience a significant increase due to the
$Z^{2}$ scaling of the photon fluxes ($Z$ being the atomic number).
Furthermore, HI collisions exhibit minimal hadronic pile-up, allowing for the
identification of exclusive events and triggering on low-$p_{\mathrm{T}}$
particles. This exceptional characteristics make UPC an excellent tool for
studying rare processes and searching for BSM phenomena. In this report two
results of the ATLAS experiment [1] are discussed: the observation of
$\gamma\gamma\rightarrow\tau^{+}\tau^{-}$ process with the measurement of
anomalous magnetic moment of the $\tau-$lepton [4] and the measurement of
light-by-light (LbyL, $\gamma\gamma\rightarrow\gamma\gamma$) scattering with
the search for axion-like particles (ALP) [5]. Both measurements utilize data
from the UPC Pb+Pb collisions and are potentially sensitive to BSM effects.

##  2 Exclusive $\tau^{+}\tau^{-}$ production and constraints on $a_{\tau}$

ATLAS provides the exclusive observation of
$\gamma\gamma\rightarrow\tau^{+}\tau^{-}$ process [4] using data from 2018
Pb+Pb collisions at $\sqrt{s_{\mathrm{NN}}}$= 5.02 TeV with an integrated
luminosity of 1.44 nb-1. The measurement of the exclusive production of
$\tau$-leptons is used to set new constraints on the anomalous magnetic moment
of the $\tau$-lepton, $a_{\tau}$. The theoretical SM prediction is:
$a_{\tau}^{\mathrm{SM}}=0.00117721(5)$ [6], which is remarkably smaller than
the currently available experimental bounds. Various BSM theories, such as
lepton compositeness, supersymmetry, and TeV-scale leptoquarks, etc.,
predicted modifications to the Standard Model (SM) value of $a_{\tau}$. The
most stringent limits on $a_{\tau}$ are currently provided by the DELPHI
experiment: $-0.052<a_{\tau}<0.013$ (95% CL) [7].

The identification techniques commonly employed in ATLAS cannot be used for
signal $\tau$-leptons due to its very low transverse momentum
($p_{\mathrm{T}}$) values. Instead, it is required that events considered in
the analysis contain one muon from $\tau$-lepton decay, and electron or
charged-particle track(s) from the other $\tau$-lepton decay. Three distinct
signal regions (SR) are defined: muon and electron ($\mu e$-SR), muon and one
track ($\mu$1T-SR), and muon and three tracks ($\mu$3T-SR). Candidate events
are selected with a single muon trigger requiring muon $p_{\mathrm{T}}$ above
4 GeV. To ensure the exclusivity of the selected events a veto on forward
neutron activity in the Zero Degree Calorimeter is imposed. Muons selected for
the analysis are required to have $p_{\mathrm{T}}>$ 4 GeV and $|\eta|<$ 2.4,
selected electrons have $p_{\mathrm{T}}>$ 4 GeV and $|\eta|<$ 2.47 and
selected tracks should have $p_{\mathrm{T}}>$ 100 MeV and $|\eta|<$ 2.5.
Events containing additional low-$p_{\mathrm{T}}$ tracks are rejected. Since
different background processes contribute to each signal category, further
requirements are introduced in the $\mu$1T-SR (muon and track system
$p_{\mathrm{T}}>$ 1 GeV) and the $\mu$3T-SR (mass of the three-track system
below 1.7 GeV). The main sources of background contributions arise from the
exclusive dimuon production with the final-state radiation (FSR) and
diffractive photonuclear interactions. The $\gamma\gamma\rightarrow\mu\mu$
background is constrained with a dimuon control region, 2$\mu$-CR. It requires
exactly two opposite-charge muons with invariant mass above 11 GeV to suppress
quarkonia backgrounds and no additional tracks separated from the muons.

Figure 1:  Muon $p_{\mathrm{T}}$ distributions in the (left) $\mu$1T-SR and
(right) 2$\mu$-CR [4]. Post-fit distributions are shown with the signal
contribution corresponding to the best-fit $a_{\tau}$ value
($a_{\tau}$=-0.042). For comparison, signal contributions with alternative
$a_{\tau}$ values are shown.

After applying the event selection, a total of 656 data events were observed
in three signal regions in which the analysis was performed. The fitted muon
$p_{\mathrm{T}}$ distributions for the $\mu$1T-SR and $2\mu$-CR are shown in
Figure 1. A very good data-to-prediction agreement is seen for the best-fit
value of the $a_{\tau}$. The $\gamma\gamma\rightarrow\tau^{+}\tau^{-}$ process
was observed with a significance exceeding 5 standard deviations, and a signal
strength of $\mu_{\tau\tau}$ = 1.03${}_{-0.05}^{+0.06}$ assuming the SM value
of $a_{\tau}$. To measure $a_{\tau}$, a fit to the muon $p_{\mathrm{T}}$
distribution is performed in the three SRs with $a_{\tau}$ being the parameter
of interest. Also a control region with events from the
$\gamma\gamma\rightarrow\mu^{+}\mu^{-}$ process is used in the fit to
constrain initial-photon fluxes. Figure 2 presents a comparison of the ATLAS
measurement of the anomalous magnetic moment of $\tau$-lepton in comparison to
previous results obtained at the LEP experiments. The precision of this
measurement is similar to the most precise single-experiment measurement by
the DELPHI Collaboration [7].

Figure 2: Measurements of $a_{\tau}$ from fits to individual signal regions
(including the $2\mu$-CR), and from the combined fit [4]. These are compared
with existing measurements from the OPAL, L3 and DELPHI experiments at LEP.

##  3 Light-by-light scattering and search for axion-like particles

The $\gamma\gamma\rightarrow\gamma\gamma$ scattering is a rare phenomenon
allowed by the quantum electrodynamics (QED) at the lowest order via a quantum
loop of virtual charged fermions or $W^{\pm}$ bosons. LbyL production can be
altered by various BSM phenomena: new particles entering the loop, Born-Infeld
extensions of the QED, space-time non-commutativity in the QED, extra spatial
dimensions, etc. Furthermore, the diphoton mass spectrum obtained from the
LbyL process can be explored to search for potential neutral axion-like
particles, ALP. ALP may contribute to the distribution as a narrow diphoton
resonance [8].

LbyL scattering was also measured by ATLAS in UPC Pb+Pb collisions at
$\sqrt{s_{\mathrm{NN}}}=$ 5.02 TeV using a combined 2015+2018 data sample with
an integrated luminosity of 2.2 nb-1 [5]. The signature of interest is the
exclusive production of two photons, each with transverse energy
$E_{\mathrm{T}}^{\gamma}>$ 2.5 GeV, pseudorapidity $|\eta^{\gamma}|<$ 2.4 and
diphoton invariant mass $m_{\gamma\gamma}>$ 5 GeV with transverse momentum
$p_{\mathrm{T}}^{\gamma\gamma}<$ 1 GeV. Any extra activity in the detector is
vetoed, in particular no reconstructed tracks originating from the nominal
interaction point with $p_{\mathrm{T}}>$ 100 MeV are accepted. The final state
photons are aligned in the azimuthal angle $\phi$. Back-to-back topology is
studied using diphoton acoplanarity defined as
$A_{\phi}=1-\frac{|\Delta\phi|}{\pi}$ . Event candidates are expected to have
$A_{\phi}<$ 0.01. The main background contribution originates from exclusive
production of the electron–positron pairs ($\gamma\gamma\rightarrow
e^{+}e^{-}$). In the measurement, the $\gamma\gamma\rightarrow e^{+}e^{-}$
background is suppressed with the requirement of no tracks and pixel-tracks
reconstructed in the Inner Detector. A remaining dielectron contribution is
evaluated using a data-driven method. The second significant background source
is gluon-induced central exclusive production (CEP) of photon pairs. The CEP
background is evaluated using a dedicated control region in data ($A_{\phi}>$
0.01) and then extrapolated to the LbyL signal region.

ATLAS established the observation of LbyL process with 97 events observed in
data, while the signal and background expectations are 45 events and 27 $\pm$
5 events, respectively. The integrated cross-section measured in the fiducial
phase space, defined by the requirements reflecting the event selection, is
$\sigma_{fid}=120\pm 17(stat.)\pm 13(syst.)\pm 4(lumi)$ nb. The presented
value can be compared with two theoretical predictions considered to be $78\pm
8$ nb from the SuperChic v3.0 MC generator [9] and $80\pm 8$ nb from [10]. In
addition to the integrated fiducial cross-section, ATLAS measured
$\gamma\gamma\rightarrow\gamma\gamma$ differential cross-sections involving
four kinematic variables of the final-state photons. In general, a good
agreement between the measurement and SM predictions is found.

ALP may be produced in the photon–photon fusion, $\gamma\gamma\rightarrow
a\rightarrow\gamma\gamma$, followed by the decay to the diphoton pair, where
$a$ denotes the ALP field. Thus, a diphoton invariant mass distribution,
$m_{\gamma\gamma}$, presented in Figure 3, can be interpreted for ALP
searches. The ALP production would result in a resonance peak with diphoton
mass equal to the mass of $a$. The diphoton mass distribution was examined for
a mass range between 6 and 100 GeV. No significant excess of events over
expected background was found in the analysis. The 95% CL limit was derived
for ALP production cross-section and ALP coupling to photons
$\frac{1}{\Lambda_{a}}$ as a function of ALP mass. A summary of exclusion
limits from different experiments together with the new ATLAS constraints is
shown in Figure 4. The new ATLAS analysis places the strongest limits on the
ALP production in the intermediate mass region to date.

Figure 3:  The diphoton invariant mass for
$\gamma\gamma\rightarrow\gamma\gamma$ event candidates [5]. Data are compared
with the sum of signal and background expectations. Figure 4:  Compilation of
exclusion limits at 95% CL in the ALP-photon coupling (1/$\Lambda_{a}$) versus
ALP mass ($m_{a}$) plane obtained by different experiments. The new limit
obtained by ATLAS is marked with label: ATLAS
$\gamma\gamma\rightarrow\gamma\gamma$ [5]

##  4 Summary

The report highlighted the significance of UPC in exploring rare SM processes
and searching for BSM phenomena. The $\gamma\gamma\rightarrow\tau^{+}\tau^{-}$
process has been observed in Pb+Pb UPC by the ATLAS experiment, surpassing a
5$\sigma$ significance. The signal strength is consistent with the Standard
Model expectations. The new constraints on the $a_{\tau}$ have been set, and
are competitive to the best limits obtained during the LEP era. With the
upcoming Run-3 data, an improvement in precision is anticipated. Additionally,
the ATLAS experiment has established the presence of
$\gamma\gamma\rightarrow\gamma\gamma$ scattering, with the results consistent
with the Standard Model prediction. The measured invariant mass of the
diphoton system was used to set new exclusion limits on axion-like particles.
This measurement provides the strongest constraints on the ALP production in
the mass region of 6–100 GeV to date.

## Acknowledgments

This work was partly supported by program ”Excellence initiative – research
university” for the AGH University of Kraków, by the National Science Centre
of Poland under grant number UMO-2021/40/C/ST2/00187 and by PL-GRID
infrastructure.
Copyright [2023] CERN for the benefit of the ATLAS Collaboration. CC-BY-4.0
license.

## References

  * [1] ATLAS Collaboration, JINST 3 (2008) S08003.
  * [2] E. Fermi, Nuovo Cim. 2 (1925) 143–158.
  * [3] E. J. Williams, Phys. Rev. 45 (1934) 729–730.
  * [4] ATLAS Collaboration, arXiv:2204.13478 [hep-ex], accepted by PRL.
  * [5] ATLAS Collaboration, JHEP 03 (2021) 243.
  * [6] S. Eidelman and M. Passera, Mod. Phys. Lett. A 22 (2007) 159–179.
  * [7] DELPHI Collaboration, J. Abdallah et al., Eur. Phys. J. C 35 (2004) 159–170.
  * [8] S. Knapen, T. Lin, H. K. Lou, and T. Melia, Phys. Rev. Lett. 118 (2017) 171801.
  * [9] L. A. Harland-Lang, V. A. Khoze, and M. G. Ryskin, Eur. Phys. J. C 79 no. 1, (2019) 39.
  * [10] M. Kłusek-Gawenda, P. Lebiedowicz, and A. Szczurek, Phys. Rev. C 93 (2016) 044907.


The Age–Metallicity–Specific Orbital Energy Relation for the
Milky Way's Globular Cluster System Confirms the Importance of Accretion
for Its Formation

Turner Woody

0000-0002-0721-6715]Turner Woody
Department of Physics and Astronomy

Johns Hopkins University

3400 N Charles St

Baltimore, MD 21218, USA

0000-0001-5761-6779]Kevin C. Schlaufman
Department of Physics and Astronomy

Johns Hopkins University

3400 N Charles St

Baltimore, MD 21218, USA

Globular clusters can form inside their host galaxies at high redshift
when gas densities were higher and gas-rich mergers were common.  They can
also form inside lower-mass galaxies that have since been accreted and
tidally disrupted, leaving their globular cluster complement bound to
higher-mass halos.  We argue that the age–metallicity–specific orbital
energy relation in a galaxy's globular cluster system can be used to
identify its origin.  Gas-rich mergers should produce tightly bound
systems in which metal-rich clusters are younger than metal-poor clusters.
Globular clusters formed in massive disks and then scattered into a halo
should have no relationship between age and specific orbital energy.
Accreted globular clusters should produce weakly bound systems in
which age and metallicity are correlated with each other but inversely
correlated with specific orbital energy.  We use precise relative ages,
self-consistent metallicities, and space-based proper motion-informed
orbits to show that the Milky Way's metal-poor globular cluster system
lies in a plane in age–metallicity–specific orbital energy space.
We find that relatively young or metal-poor globular clusters are weakly
bound to the Milky Way, while relatively old or metal-rich globular
clusters are tightly bound to the Galaxy.  While metal-rich globular
clusters may be formed either in situ or ex situ, our results suggest that
metal-poor clusters formed outside of the Milky Way in now-disrupted dwarf
galaxies.  We predict that this relationship between age, metallicity,
and specific orbital energy in a $L^{*}$ galaxy's globular cluster system
is a natural outcome of galaxy formation in a $\Lambda$CDM universe.

§ INTRODUCTION

Globular cluster formation is one the enduring unsolved problems
in astrophysics.  Globular clusters can form either in situ
inside their parent galaxy or ex situ in lower-mass galaxies that
have been accreted and tidally disrupted by a more massive halo
<cit.>.  In either case, the ancient
ages inferred for most globular clusters indicate that they formed long
ago when the universe was very different from its current $z = 0$ state.
Indeed, the interval spanning the 16th and 84th percentiles of the Milky
Way globular cluster system's age distribution is $11~\text{Gyr} \lesssim
\tau \lesssim 13~\text{Gyr}$ or $2.4 \lesssim z \lesssim 7.3$ in terms
of redshift <cit.>.
Both in situ and ex situ formation channels are theoretically
expected to contribute to the Milky Way's globular cluster system

In the in situ formation scenario, globular clusters
form inside their parent galaxy in gas-rich galaxy mergers
and/or in gas-rich disks
Galaxies in general and galaxies destined to be similar to the Milky Way
in particular are observed to be smaller at high redshift when globular
clusters form in this scenario <cit.>.
Consequently, globular clusters formed in this way should be tightly
bound to their now more-massive and larger parent galaxies at $z = 0$
<cit.>.  Another robust prediction of the gas-rich
merger scenario is that metal-rich clusters should be systematically
younger than metal-poor clusters [Li & Gnedin, 2014, Li & Gnedin, 2019].

In the ex situ scenario, globular clusters form
outside their $z = 0$ host halo in parent galaxies
that have since been accreted and tidally disrupted
While in most cases it is difficult to conclusively associate a candidate
accreted globular cluster with its parent galaxy, the globular clusters
Arp 2, NGC 6715 (M 54), Pal 12, Terzan 7, and Terzan 8 have all been
securely associated with the Sagittarius dwarf spheroidal (dSph) galaxy
<cit.>.  Given the collisionless nature of
their accretion, globular clusters formed ex situ should be on average
less tightly bound to their host galaxy at $z = 0$ than clusters formed
in situ in gas-rich conditions.

The observed properties of accreted globular clusters can be used to
set limits on the properties of their parent galaxies.  Likewise, the
ensemble properties of a galaxy's accreted globular cluster system
can be used to explore its formation.  The maximum metallicity
realized in a dwarf galaxy's globular cluster system should not
exceed the dwarf galaxy's typical metallicity.  This is the case in
both the Fornax and Sagittarius dSph galaxies.  The Fornax dSph has
$\langle[\text{Fe/H}]\rangle = -1.04$ [Kirby et al., 2013], while its
five globular clusters span the range $-2.5 \lesssim [\text{Fe/H}]
\lesssim -1.4$ <cit.>.
While the Sagittarius dSph has $\langle[\text{Fe/H}]\rangle
\approx -0.5$ <cit.>,
the five globular clusters securely associated with it span
the metallicity range $-2.3 \lesssim [\text{Fe/H}] \lesssim -0.6$
A lower limit on the stellar mass of a galaxy necessary for the formation
of a globular cluster of a given metallicity can therefore be derived
from a stellar mass–metallicity relationship.  Kirby et al., 2013 showed
that there exists a linear relationship between log stellar mass and
metallicity for dwarf galaxies over five orders of magnitude such that
the typical metallicity of a dwarf galaxy rises linearly with log stellar
mass from $[\text{Fe/H}] \approx -2.5$ at $M_{\ast} \sim 10^{3}~M_{\odot}$
to $[\text{Fe/H}] \approx -1.0$ at $M_{\ast} \sim 10^{8}~M_{\odot}$.  So
while both metal-rich and metal-poor globular clusters can form in massive
dSph galaxies, low-mass dSph galaxies will only be able to form metal-poor
clusters.  More sophisticated arguments based on cosmological dark
matter-only or hydrodynamical simulations have reached similar conclusions

The mass of an accreted globular cluster's parent galaxy will also
affect its orbit.  A dwarf galaxy orbiting inside the dark matter halo
of a more massive galaxy will lose angular momentum due to dynamical
friction and its orbit will decay on a timescale that is proportional
to its mass <cit.>.  Massive dSph galaxies
will therefore be dragged into the inner regions of their host halo
and tidally disrupted much more quickly than low-mass dSph galaxies.
In particular, a massive dSph galaxy resembling Fornax with a total
mass $M_{\mathrm{tot}} \sim 10^{8}~M_{\odot}$ capable of forming a
globular cluster with $[\mathrm{Fe/H}] \approx -1.0$ will be tidally
disrupted 100 times faster than a low-mass dSph galaxy with a total
mass $M_{\mathrm{tot}} \sim 10^{6}~M_{\odot}$ that could not form a
globular cluster more metal-rich than $[\mathrm{Fe/H}] \approx -2.0$
<cit.>.  The metal-rich globular
clusters deposited deep inside the more massive halo when their parent
massive dSph galaxy is tidally disrupted will be much more tightly bound
to the massive halo than any metal-poor globular clusters left behind
when their parent low-mass dSph galaxy is tidally disrupted much later.

The ages of accreted globular clusters should be related to the masses
and therefore the metallicities of their parent galaxies too.  In similar
environments, galaxies destined to be massive form in higher $\sigma$
peaks in the primordial matter density distribution than galaxies destined
to be low mass <cit.>.  Since the dynamical time
$t_{\mathrm{dyn}}$ scales like $t_{\mathrm{dyn}} \propto \rho^{-1/2}$,
in similar environments more massive galaxies form stars earlier and
therefore the ages of the oldest stellar populations in galaxies should
be correlated with their masses.  The implication is that in similar
environments the age of the oldest globular cluster formed in a dwarf
galaxy should be correlated with its mass.  This is the case for the
$M_{\mathrm{tot}} \approx 9 \times 10^{7}~M_{\odot}$ Fornax dSph galaxy,
the $M_{\mathrm{tot}} \gtrsim 4 \times 10^{8}~M_{\odot}$ Sagittarius
dSph galaxy, and the $M_{\mathrm{tot}} \sim 10^{11}~M_{\odot}$ Large
Magellanic Cloud <cit.>.
The Fornax globular clusters 1, 3, and 5 are comparable
in age to NGC 4590 (M 68) while cluster 4 is somewhat younger
<cit.>.  On the other hand, Terzan 8
is the oldest globular cluster with a confirmed Sagittarius association
and is older than NGC 4590 <cit.>.
Likewise, the oldest of the Large Magellanic Cloud's globular
clusters that are not affected by crowding, Hodge 11, is older than
Terzan 8 but younger than the Milky Way's oldest globular clusters
[Wagner-Kaiser et al., 2017].  Analyses based on scaling relations
calibrated with cosmological hydrodynamical simulations reach similar
conclusions <cit.>.

The relationship between age and metallicity in the Galaxy's globular
cluster system has long been used to constrain Milky Way formation.
As we described above, analyses of the relationship between age,
metallicity, and orbital properties in the Milky Way's globular cluster
system can be even more informative.  If globular clusters formed in
situ are produced in gas-rich mergers, then metal-rich clusters should
be younger than metal-poor clusters and both metal-rich and metal-poor
clusters should be relatively tightly bound to their parent galaxy.
If globular clusters are formed in situ in gas-rich disks and subsequently
scattered into the halo, then there should be no relationship between
age and specific orbital energy.  In the ex situ formation and accretion
scenario, older, metal-rich globular clusters should be tightly bound to
their $z = 0$ host galaxies while younger, metal-poor globular clusters
should be weakly bound to their $z = 0$ host.  These straightforward
arguments have been confirmed by cosmological hydrodynamical simulations

In this paper, we calculate the orbits of the Milky Way's globular
clusters with precise space-based astrometry.  We use those orbits
along with precise relative ages and metallicities from the literature
to quantify the relationship between age, metallicity, and specific
orbital energy in the Milky Way's globular cluster system.  We describe
the construction of our globular cluster sample in Section <ref>.
We detail in Section <ref> our orbital integrations and
statistical analyses of the age–metallicity–specific orbital energy
relation in our complete sample of globular clusters as well as its
metal-rich/metal-poor subsamples.  We review our results and their
implications for the formation of the Milky Way's globular cluster
system specifically and globular cluster formation in general in Section
<ref>.  We conclude by summarizing our findings in Section

§ DATA

In an analysis of the age–metallicity–specific orbital energy
relation in the Milky Way's globular cluster system, the limiting
factors are precise cluster ages and proper motions.  Most of
our proper motion measurements come from Gaia Data Release 2 <cit.>.
We supplement the globular cluster proper motions published in
Gaia Collaboration et al., 2018 with the Gaia DR2-based proper motions for NGC 6584 and
NGC 6723 from Baumgardt et al., 2019 that have proper motion precisions
$\mu/\sigma_{\mu} > 85$.[We use this threshold because all
of the globular cluster proper motions from Gaia Collaboration] et al., 2018 have
$\mu/\sigma_{\mu} > 85$ in the Baumgardt et al., 2019 catalog. We also
add proper motions measurements for NGC 1261, NGC 4147, NGC 6101, Terzan
7, Arp 2, Terzan 8, NGC 6934, and Pal 12 from the High-resolution Space
Telescope Proper Motion Collaboration <cit.>.  We
list our input globular cluster proper motions in Table <ref>.

Globular cluster age determinations are inherently difficult to put
on an absolute age scale and can be non-trivially affected by even
small differences in analysis techniques.  Since the analyses we
will describe in Section <ref> only make use of relative
age differences, we choose to use the homogeneous relative globular
cluster ages presented in Marín-Franch et al., 2009 to avoid the difficulties
associated with absolute ages.  Marín-Franch et al., 2009 calculated
relative clusters ages $\tau = \text{age}/\text{12.8~Gyr}$ for two
different metallicity scales [Zinn & West, 1984, Carretta & Gratton, 1997]
as reported in Rutledge et al., 1997, Rutledge et al., 1997
and four different sets of theoretical stellar models
[Bertelli et al., 1994, Girardi et al., 2000, Pietrinferni et al., 2004, Dotter et al., 2007].
We list our adopted globular cluster metallicities and normalized ages
in Table <ref>.  While we focus on the relative ages
calculated assuming the more modern Carretta & Gratton, 1997 metallicity
scale and Dotter et al., 2007 models given in Table <ref>,
the analysis we describe in Section <ref> reaches the same
conclusion regardless of the assumed metallicity scale or model grid.

We adopt the globular cluster distances and radial
velocities listed in Table <ref> from
the December 2010 revision of the Harris, 1996
compilation.[<https://physwww.mcmaster.ca/ harris/mwgc.ref>]
Since the December 2010 revision of the Harris, 1996 compilation
does not provide uncertainties for its distance estimates, we estimate
distance uncertainties by first transforming the distance inferences into
distance moduli.  We next assume a uniform distance modulus uncertainty
of 0.05 mag and sample each cluster's distance modulus from a normal
distribution centered on its distance modulus with standard deviation
0.05 mag.  We then transform the simulated distance moduli calculated
in this way back into distances and uncertainties.  We will use the
57 Milky Way globular clusters listed in Tables <ref> and
<ref> with homogeneous relative ages and metallicities plus
precise space-based proper motions to explore the relationship between
age, metallicity, and specific orbital energy in Section <ref>.

§ ANALYSIS

We use the right ascensions $\alpha$, declinations $\delta$, proper
motions $\mu_{\alpha}$ & $\mu_{\delta}$, distances $d$, and radial
velocities $v_{r}$ given for the 57 globular clusters listed in Table
<ref> to calculate their orbits and specific orbital energies.
We use a Monte Carlo approach to average over the uncertainties in
our input proper motions, distances, and radial velocities.  On each
iteration, we sample the necessary input data for an orbit integration
($\alpha$, $\delta$, $\mu_{\alpha}$, $\mu_{\delta}$, $d$, and $v_{r}$)
from the data and uncertainties in Table <ref>.  When possible,
we account for the covariance between simultaneously inferred proper
motion and parallax inferences.  For the globular clusters with proper
motions from Gaia Collaboration et al., 2018, we account for the covariance between
proper motions and parallaxes even though we use Harris, 1996-based
distances instead of the Gaia DR2 parallax-based distances.  For the
globular clusters with proper motions from Baumgardt et al., 2019,
we account for the covariance between the proper motion components.
Sohn et al., 2018 did not provide estimates of the covariance between
proper motion components, so we sample each proper motion component
independently from within its quoted uncertainty.  For distances and
radial velocities, we sample from normal distributions with the means
and standard deviations given in Table <ref>.  We repeat this
process 100 times to generate 100 sets of initial conditions for each
globular cluster.

For each Monte Carlo initial condition realization, we integrate a
cluster's orbit in model Milky Way potentials using the
python module[<https://github.com/jobovy/galpy>]
[Bovy, 2015].  We calculate orbits from the present five Gyr
into the past with a time step of 200,000 years assuming a solar
motion with respect to the local standard of rest ($U_{\odot}$,
$V_{\odot}$, $W_{\odot}$) = (11.1, 12.24, 7.25) km s$^{-1}$

We calculate orbits in three different model Milky Way potentials:
the default  from Bovy, 2015, a scaled
version of  with its halo mass increased by 50%,
and  from McMillan, 2017.  In the default
potential, the virial mass and radius are
$M_{\text{vir}} = 0.8 \times 10^{12}~M_{\odot}$ and $r_{\text{vir}}
= 245$ kpc respectively.  Its bulge is parameterized as a power-law
density profile with exponent $-1.8$ that is exponentially cut-off at
1.9 kpc.  Its disk is represented by a Miyamoto–Nagai potential with
a radial scale length of 3 kpc and a vertical scale height of 280 pc
[Miyamoto & Nagai, 1975].  Its halo is modeled as a Navarro–Frenk–White
(NFW) halo with a scale length of 16 kpc [Navarro et al., 1996].  The scaled
version of the  potential differs from the default
potential in that we have increased its halo
mass by 50%.  This change increases the virial mass to $M_{\text{vir}} =
1.2 \times 10^{12}~M_{\odot}$.  In both cases, we set the solar distance
to the Galactic center to $R_{0}=8.122$ kpc, the circular velocity at
the Sun to $V_{0}=238$ km s$^{-1}$, and the height of the Sun above the
plane to $z_{0}=25$ pc [Jurić et al., 2008, Bland-Hawthorn & Gerhard, 2016, Gravity Collaboration et al., 2018].

The  potential has $M_{\text{vir}} = 1.4 \times
10^{12}~M_{\odot}$.  Its bulge is parameterized as a spherical power
law density profile with exponent $-1.8$ that is finite at $r = 0$ and
has an exponential cut-off at 2.1 kpc.  Its disk is split into four
components: a thin stellar disk, a thick stellar disk, an H1
gas disk, and an H$_{2}$ gas disk.  Its two stellar disks are modeled
as exponential in both the cylindrical radius $R$ and height $z$, with
scale lengths $R_{\text{d,thin}}=2.6$ kpc, $R_{\text{d,thick}}=3.6$ kpc,
$z_{\text{d,thin}}=300$ pc, and $z_{\text{d,thick}}=900$ pc.  Its two
gas disks are modeled as exponential in $R$ with a central hole and a
modified exponential profile in the $z$ coordinate.  Its halo is modeled
as a NFW halo with a scale length of 19.6 kpc.  The
potential uses a solar distance and circular velocity of 8.21 kpc and
233.1 km s$^{-1}$.

The Milky Way's mass within 25 kpc $M_{25}$ is especially relevant
for this problem, and estimates of $M_{25}$ for the Milky Way
based on globular cluster kinematics suggest that $M_{25} \approx
0.26_{-0.06}^{+0.10} \times 10^{12}~M_{\odot}$ [Eadie & Jurić, 2019].
The $M_{25}$ values for the default , the
scaled , and the  potentials
are $M_{25} = 0.25 \times 10^{12}~M_{\odot}$, $M_{25} = 0.29 \times
10^{12}~M_{\odot}$, and $M_{25} = 0.28 \times 10^{12}~M_{\odot}$
respectively.  All three potentials are therefore consistent the
best available Milky Way constraints in the same volume.  While
both the default  and
potentials are generally consistent with literature estimates of
the Milky Way's virial mass, most recent mass estimates making use
of Gaia DR2 proper motions suggest virial masses in the range $1.1
\times 10^{12}~M_{\odot} \lesssim M_{\text{vir}} \lesssim 1.4 \times
10^{12}~M_{\odot}$ <cit.>.  While the virial mass
of the  potential is in this interval, the virial
mass of the default  potential is below this
range.  In addition, the  potential provides a more
comprehensive model of the Milky Way's disk than either the default or
scaled  potentials.  This difference could be
important for the relatively nearby clusters that dominate our input
sample.  For all of these reasons, the  potential
likely provides a better approximation to the true Milky Way potential
than the default and scaled  potentials.

Our orbit integrations produce 100 plausible orbits for each of
the 57 globular clusters listed in Table <ref> in the
default , scaled , and
potentials.  We calculate the 16th, 50th, and 84th
percentiles of the specific orbital energy (SOE) distribution for every
cluster and use those data to define each cluster's specific orbital
energy and uncertainty in each potential.  We list the specific
orbital energies derived in this way for each cluster in Table

We use the  python module [Seabold & Perktold, 2010] to fit
an ordinary least squares linear regression of the form $\mathrm{SOE}
= \beta_0 + \beta_{\tau} \tau + \beta_{M} \mathrm{[M/H]}$ to the
age–metallicity–specific orbital energy distribution produced by
each set of initial conditions in all three potentials.  We aggregate
the results and calculate the 16th, 50th, and 84 percentiles of each
distribution.  We report the results in Table <ref>.  In all
three potentials, we find evidence of an age–metallicity–specific
orbital energy relationship in which relatively old, metal-rich globular
clusters have more negative specific orbital energies and are therefore
more tightly bound to the Milky Way than relatively metal-poor,
young globular clusters.  This relationship is especially strong in
the  potential, as $\beta_{\tau}$ and $\beta_{M}$
in that potential are significant at approximately the 5-$\sigma$ level.

We separately study the age–metallicity–specific orbital energy
relationship for metal-rich and metal-poor clusters.  We split the
complete sample into metal-rich/metal-poor subsamples at $[\text{M/H}]
= -0.8$ because that cleanly separates the two peaks of the complete
sample's bimodal metallicity distribution.  We find similar results
for both the complete sample and the metal-poor subsample of globular
clusters that we define as the 45 globular clusters with $[\text{M/H}]
\leq -0.8$.  The $t$ values and therefore the statistical significances
of the coefficients defining the relation are generally larger in
absolute value for the metal-poor subsample.  For the 12 metal-rich
clusters with $[\text{M/H}] > -0.8$, the age–metallicity–specific
orbital energy relationship is much weaker.  Indeed, if we exclude the
two metal-rich globular clusters associated with the Sagittarius dSph,
then in both the default and scaled  potentials
the coefficients defining the age–metallicity–specific orbital relation
energy change signs such that older, more metal-rich clusters are less
tightly bound to the Milky Way.  In this case, there is no significant
age–metallicity–specific orbital energy relationship for the the

Median Regression Model Coefficients and $t$ Values

$t$ Value
$t$ Value
$t$ Value

Complete     -5.39 $\pm$ 2.32     -2.33     -4.12 $\pm$ 2.54     -1.62     -1.16 $\pm$ 0.63     -1.84     0.081

Metal-poor     -2.21 $\pm$ 3.07     -0.72     -8.64 $\pm$ 3.62     -2.39     -1.98 $\pm$ 0.92     -2.16     0.146

Metal-richa     -11.26 $\pm$ 5.19     -2.17     -1.38 $\pm$ 4.20     -0.33     -7.45 $\pm$ 8.09     -0.92     0.087

Metal-richb     -22.52 $\pm$ 10.02     -2.25     19.30 $\pm$ 10.80     1.79     9.34 $\pm$ 10.16     0.92     0.317

Complete     -10.17 $\pm$ 2.44     -4.16     -3.95 $\pm$ 2.68     -1.47     -1.40 $\pm$ 0.67     -2.10     0.089

Metal-poor     -5.89 $\pm$ 3.19     -1.85     -9.53 $\pm$ 3.76     -2.54     -2.15 $\pm$ 0.95     -2.26     0.160

Metal-richa     -18.20 $\pm$ 5.42     -3.36     -0.12 $\pm$ 4.38     -0.03     -9.53 $\pm$ 8.46     -1.13     0.131

Metal-richb     -30.52 $\pm$ 9.72     -3.14     22.91 $\pm$ 10.48     2.19     9.52 $\pm$ 9.85     0.97     0.407

Complete     1.39 $\pm$ 3.95     0.35     -24.52 $\pm$ 4.33     -5.66     -4.85 $\pm$ 1.08     -4.49     0.435

Metal-poor     2.80 $\pm$ 5.23     0.54     -28.18 $\pm$ 6.17     -4.57     -6.32 $\pm$ 1.56     -4.05     0.381

Metal-richa     -3.86 $\pm$ 9.24     -0.42     -24.76 $\pm$ 7.47     -3.32     -15.79 $\pm$ 14.41     -1.10     0.551

Metal-richb     -38.61 $\pm$ 19.16     -2.02     15.39 $\pm$ 20.67     0.745     -5.92 $\pm$ 19.43     -0.31     0.123
Median age–metallicity–specific orbital energy
model coefficients, coefficient uncertainties, $t$ values, and $R^2$
values calculated from 100 Monte Carlo iterations for the complete,
metal-poor, and metal-rich subsamples (both including and excluding
clusters associated with the Sagittarius dSph) assuming all three
aAll metal-rich globular clusters
bMetal-rich globular clusters excluding those associated
with the Sagittarius dSph

We use $F$ tests to compare the full age–metallicity–specific orbital
energy model we derived above with each of its nested submodels to verify
that the full age–metallicity–specific orbital energy model best
describes the data in Table <ref>.  For each subsample
in each potential, we tested three null hypotheses: $\beta_{\tau}=0$,
$\beta_M=0$, and $\beta_{\tau}=\beta_M=0$.  We report the results of
those model comparisons in Table <ref>.  For the
metal-poor subsample in all three potentials, the $p$ values for
each test are small enough to confidently reject the null hypotheses.
The implication is that the full linear model is a better description of
the metal-poor subsample than any of its nested submodels.  In the case
of the  potential, the full linear model is a better
description of the complete sample than any of its nested submodels.
On the other hand, for the metal-rich subsample there is no reason to
prefer the full age–metallicity–specific orbital energy model over any
of its nested submodels.  We also use the Bayesian information criterion
(BIC) to compare models.  Smaller BIC values indicate better models,
and the results in Table <ref> show that the full
age–metallicity–specific orbital energy model is preferred to any of its
nested submodels for both the complete sample and metal-poor subsample
in the  potential.  All models are approximately
similar in the default and scaled  potentials.

Median Model Comparison $F$ Statistics

$H_0:\beta_{\tau} = 0$
$H_0:\beta_M = 0$
$H_0:\beta_{\tau} = \beta_M = 0$

Complete     2.63     $1.1\times10^{-1}$     3.27     $7.2\times10^{-2}$     2.38     $1.0\times10^{-1}$

Metal-poor     5.71     $2.1\times10^{-2}$     4.66     $3.7\times10^{-2}$     3.58     $3.7\times10^{-2}$

Metal-richa     0.11     $7.5\times10^{-1}$     0.85     $3.8\times10^{-1}$     0.43     $6.6\times10^{-1}$

Metal-richb     3.19     $1.1\times10^{-1}$     0.85     $3.9\times10^{-1}$     1.62     $2.6\times10^{-1}$

Complete     2.17     $1.5\times10^{-1}$     4.40     $4.1\times10^{-2}$     2.65     $8.0\times10^{-2}$

Metal-poor     6.44     $1.5\times10^{-2}$     5.12     $2.9\times10^{-2}$     4.00     $2.6\times10^{-2}$

Metal-richa     0.00     $9.8\times10^{-1}$     1.27     $2.9\times10^{-1}$     0.68     $5.3\times10^{-1}$

Metal-richb     4.78     $6.5\times10^{-1}$     0.93     $3.7\times10^{-1}$     2.40     $1.6\times10^{-1}$

Complete     32.05     $5.9\times10^{-7}$     20.13     $3.8\times10^{-5}$     20.81     $2.0\times10^{-7}$

Metal-poor     20.89     $4.2\times10^{-5}$     16.38     $2.2\times10^{-4}$     12.90     $4.3\times10^{-5}$

Metal-richa     11.00     $9.0\times10^{-3}$     1.20     $3.0\times10^{-1}$     5.52     $3.0\times10^{-2}$

Metal-richb     0.55     $4.8\times10^{-1}$     0.09     $7.7\times10^{-1}$     0.49     $6.3\times10^{-1}$
$F$ statistics and associated $p$ values comparing
the full three parameter age–metallicity–specific orbital energy
relation derived for each potential to all possible lower-dimensional
submodels.  The large $F$ statistics and small $p$ values indicate
that the full age–metallicity–specific orbital energy relation is
superior to all lower-dimensional nested models for the metal-poor
aAll metal-rich globular clusters
bMetal-rich globular clusters excluding those associated
with the Sagittarius dSph

Median Model Comparison BIC Statistics

Full Model
$\beta_{\tau} = 0$
$\beta_M = 0$
$\beta_{\tau} = \beta_M = 0$

Complete     241     240     240     238

Metal-poor     188     190     189     188

Metal-richa     55     53     54     52

Metal-richb     43     44     41     42

Complete     247     245     248     244

Metal-poor     192     194     193     192

Metal-richa     57     54     56     53

Metal-richb     42     45     41     43

Complete     302     324     316     326

Metal-poor     236     251     247     250

Metal-richa     69     76     68     74

Metal-richb     55     54     53     52
BIC values comparing the full three parameter
age–metallicity–specific orbital energy relation derived for each
potential to all possible lower-dimensional submodels.  The small BIC
values indicate that the full age–metallicity–specific orbital energy
relation is superior to all lower-dimensional nested models for both the
complete sample and its metal-poor subsample in the
aAll metal-rich globular clusters
bMetal-rich globular clusters excluding those associated
with the Sagittarius dSph

We plot the distribution of the Milky Way's globular cluster systems
in the age–metallicity–specific orbital energy space for all three
potentials in Figure <ref>.  It is clear that the 57 globular
clusters listed in Tables <ref> and <ref>
lie in a plane in age–metallicity–specific orbital energy
space in the  potential, exactly as predicted
by our linear modeling.  We plot two-dimensional projections of
age–metallicity–specific orbital energy space for the default
, scaled , and
potentials in Figures <ref>, <ref>,
and <ref>.  There remains significant dispersion about the
linear model visualized in Figures <ref>, <ref>,
and <ref>, so the age–metallicity–specific orbital energy
relation inference we described in Section 1 does not fully explain the
highly nonlinear process of globular cluster formation.  Nevertheless, the
meaningfully non-zero $R^2$ values presented in Table <ref> show
that a simple three-parameter linear model explains at least 15%—and
possibly as much as 38%—of the variance in the metal-poor globular
cluster age–metallicity–specific orbital energy relation.

Age–metallicity–specific orbital energy distribution of
the Milky Way globular clusters listed in Tables <ref> and
<ref>.  We plot metal-poor clusters with $[\text{M/H}] \leq
-0.8$ as blue circles and metal-rich clusters with  $[\text{M/H}] > -0.8$
as green triangles.  Downward-pointing triangles are metal-rich clusters
associated with the Sagittarius dSph galaxy, while upward-pointing
triangles are metal-rich clusters unassociated with the Sagittarius
dSph galaxy.  Each plot is looking down the plane defined by the median
best-fit age–metallicity–specific orbital energy relation for the
metal-poor subsample as reported in Table <ref>.  Top: specific
orbital energies calculated assuming the default
potential from Bovy, 2015.  Middle: specific orbital
energies calculated assuming the scaled
potential.  Bottom: specific orbital energies calculated assuming the
potential from McMillan, 2017.  It is visually
clear in the  potential that the metal-poor subsample
lies in a plane in age–metallicity–specific orbital energy space as
indicated by the highly significant coefficients $\beta_{\tau}$ and
$\beta_{M}$ listed in Table <ref> and the small $F$ test $p$
values listed in Table <ref>.

Two-dimensional projections of the age–metallicity–specific
orbital energy distribution of the 57 Milky Way globular clusters listed
in Tables <ref> and <ref> assuming the default
potential [Bovy, 2015].  We plot metal-poor
clusters as blue circles and metal-rich clusters as green triangles.
Downward-pointing triangles are metal-rich clusters associated with the
Sagittarius dSph galaxy, while upward-pointing triangles are metal-rich
clusters unassociated with the Sagittarius dSph galaxy.  The sizes of
the points corresponds to the relative values of the suppressed axis.
In blue we indicate the two-dimensional projections of the best-fit
age–metallicity–specific orbital energy relation in the metal-poor
subsample for each of the 100 Monte Carlo trials.  Top: metallicity–age
relation.  Middle: metallicity–specific orbital energy relation.
Bottom: age–specific orbital energy relation.

Two-dimensional projections of the age–metallicity–specific
orbital energy distribution of the 57 Milky Way globular clusters
listed in Tables <ref> and <ref> assuming
the scaled  potential.  We plot metal-poor
clusters as blue circles and metal-rich clusters as green triangles.
Downward-pointing triangles are metal-rich clusters associated with the
Sagittarius dSph galaxy, while upward-pointing triangles are metal-rich
clusters unassociated with the Sagittarius dSph galaxy.  The sizes of
the points corresponds to the relative values of the suppressed axis.
In blue we indicate the two-dimensional projections of the best-fit
age–metallicity–specific orbital energy relation in the metal-poor
subsample for each of the 100 Monte Carlo trials.  Top: metallicity–age
relation.  Middle: metallicity–specific orbital energy relation.
Bottom: age–specific orbital energy relation.

Two-dimensional projections of the age–metallicity–specific
orbital energy distribution of the 57 Milky Way globular clusters
listed in Tables <ref> and <ref> assuming
potential [McMillan, 2017].  We plot metal-poor
clusters as blue circles and metal-rich clusters as green triangles.
Downward-pointing triangles are metal-rich clusters associated with the
Sagittarius dSph galaxy, while upward-pointing triangles are metal-rich
clusters unassociated with the Sagittarius dSph galaxy.  The sizes of
the points corresponds to the relative values of the suppressed axis.
In blue we indicate the two-dimensional projections of the best-fit
age–metallicity–specific orbital energy relation in the metal-poor
subsample for each of the 100 Monte Carlo trials.  Top: metallicity–age
relation.  Middle: metallicity–specific orbital energy relation.
Bottom: age–specific orbital energy relation.

§ DISCUSSION

We find a significant age–metallicity–specific orbital energy
relationship for the metal-poor subsample of 45 globular clusters
with $[\text{M/H}] \leq -0.8$ listed in Tables <ref>
and <ref>.  In particular, old, metal-rich globular
clusters tend to have more negative specific orbital energies
and are therefore more tightly bound to the Milky Way than young,
metal-poor globular clusters.  The metal-poor subsample prefers the
full age–metallicity–specific orbital energy relation over any of its
nested submodels.  These results are independent of assumed potential,
although the relationship is much stronger in the
potential from McMillan, 2017 than either the default or scaled
potentials from Bovy, 2015.  We find
no clear age–metallicity–specific orbital energy relation in the
metal-rich subsample of 10 globular clusters with $[\text{M/H}] > -0.8$
that are not associated with the Sagittarius dSph galaxy.  We also find
in this metal-rich subsample that an age–metallicity–specific orbital
energy relation is not clearly preferred over lower-dimensional models
relating age, metallicity, or specific orbital energy.

As we argued in Section <ref>, accreted metal-rich globular
clusters can only be formed in massive dSph galaxies that experience
stronger dynamical friction and form their first stars more quickly
than lower-mass dSph galaxies.  These two properties imply that in
a globular cluster system with many accreted globular clusters, age
and metallicity should be correlated with each other but inversely
correlated with specific orbital energy.  Our observation that old,
metal-rich globular clusters are more tightly bound to the Milky
Way than young, metal-poor globular clusters confirms the importance
of accretion for the formation of the Milky Way's globular cluster
system.  Our inference also supports previous proposals that a large
fraction of the Milky Way's metal-poor globular clusters were accreted

All but two of the metal rich clusters in our sample have tightly bound
orbits and old ages, broadly consistent with an in situ origin for
metal-rich globular clusters.  These two clusters—Terzan 7 and Pal
12—are both young and on energetic orbits.  They have been securely
associated with the Sagittarius dSph galaxy though, and therefore have
been accreted <cit.>.  Because we do not
observe a strong age–metallicity–specific orbital energy for metal-rich
globular clusters unassociated with the Sagittarius dSph, our inferences
do not support an accretion origin for most metal-rich globular clusters.
As there are only 10 metal-rich clusters unassociated with the Sagittarius
dSph in our sample, we are unable to choose between gas-rich mergers or
formation in gas-rich disks for the origin of the Milky Way's metal-rich
globular clusters.

In our analysis, we used relative ages from Marín-Franch et al., 2009 calculated
assuming the Carretta & Gratton, 1997 metallicity scale and the theoretical
model grid from Dotter et al., 2007.  To ensure that our results are
insensitive to the assumed metallicity scale or theoretical model grid
underlying the Marín-Franch et al., 2009 age inferences, we performed the same
analysis using all of the relative ages from Marín-Franch et al., 2009 calculated
assuming both the Zinn & West, 1984 and Carretta & Gratton, 1997 metallicity
scales and all of the Bertelli et al., 1994, Girardi et al., 2000,
Pietrinferni et al., 2004, and Dotter et al., 2007 theoretical model
grids.  In all cases, our analysis produces the same trends with
comparable statistical significance.  This should not be surprising,
as Marín-Franch et al., 2009 found that their relative ages inferences were
insensitive to their choice of metallicity scale or theoretical model
grid.  This consistency supports the robustness of our findings.

As a consequence of the simple physical explanation we outline to explain
the age–metallicity–specific orbital energy relation we observe in the
Milky Way's globular cluster system, we predict that the globular cluster
systems of other $L^{\ast}$ galaxies should evince similar relations.
That is, we predict that other $L^{\ast}$ galaxies in isolation or in
low-density group environments like the Local Group should have globular
cluster systems in which old, metal-rich globular clusters have more
negative specific orbital energies and are therefore more tightly bound to
their host halo than young, metal-poor globular clusters.  We argue that
this is a natural outcome of galaxy formation in a $\Lambda$CDM universe.

While we infer the importance of accretion for the origin of
the Milky Way's globular cluster system from an analysis of its
ensemble properties, our approach cannot associate individual
globular clusters with individual accretion events.  Detailed
simulations of Milky Way analogs have shown how a globular cluster
systems properties can be used to reveal that halo's accretion history
When applied to the Milky Way, that formalism can quantify the properties
of massive satellite galaxies that have since merged with the Milky Way
Moreover, efforts have been made to match the kinematics of clusters to
stellar streams associated with known accretion events.  For example,
Massari et al., 2019 employed this technique in combination with
in-group age–metallicity relations to assign candidate accreted globular
clusters to specific progenitors.  However, our analysis shows that it is
beneficial to use similar techniques on the entire cluster population,
thereby capturing the general imprint that accretion has had on the
Milky Way's globular cluster system's age–metallicity–specific orbital
energy relation.

Because proper motions are more difficult to infer for distant clusters,
our input sample lacks outer halo clusters.  This bias is unlikely to
affect our conclusion, as the only plausible explanation for weakly
bound outer halo clusters is accretion from dwarf galaxies.  Our input
sample also lacks highly extincted disk and bulge clusters.  This bias
is unlikely to affect our conclusion either, as disk/bulge clusters are
metal-rich and would not appear in our metal-poor subsample even if they
were present in our input sample.

The age–metallicity–specific orbital energy relation we observe in
our metal-poor subsample of globular clusters is much stronger in the
potential than in the
potential.  As discussed in Section <ref>, the default
potential represents a comparatively
smaller, lower-mass Milky Way than the  potential.
In addition, the more massive  potential appears
to be a better match to Gaia DR2-informed Milky Way virial mass
inferences.  Under the default  potential,
most of the weakly bound clusters are relatively less bound than they
are in the  potential, though this difference is
diminished in the scaled .  This results in
a relatively larger spread in orbital energy, which might explain the
weaker age–metallicity–specific orbital energy relation.  Nevertheless,
the age–metallicity–specific orbital energy relation has the same form
in all three potentials and we are confident its existence is robust to
the choice of potential.

§ CONCLUSION

We find that in the metal-poor subsample of 45 Milky Way globular
clusters with $[\text{M/H}] \leq -0.8$ in Tables <ref> and
<ref>, relatively young or metal-poor globular clusters
are weakly bound to the Milky Way, while relatively old or metal-rich
globular clusters are tightly bound to the Galaxy.  We argue that this
relationship is naturally explained by the accretion of globular clusters
from now-disrupted dwarf galaxies.  We propose that this observation is
a consequence of the combined effects of the stellar mass–metallicity
relationship for dwarf galaxies, the dependence of dynamical friction
on mass, and the relationship between a galaxy's dynamical time and its
central density.  Accreted metal-rich globular clusters can only form
in massive dSph galaxies that can attain high metallicities and will be
strongly affected by dynamical friction.  These more massive dSph galaxies
capable of forming metal-rich clusters will also produce their oldest
stellar populations before lower-mass dSph galaxies that form in lower
$\sigma$ peaks in the universe's primordial matter density distribution.
These lower-mass dSph galaxies will be much less affected by dynamical
friction and can only produce low-metallicity globular clusters.
Because these properties are independent of the specific accretion
history of the Milky Way, we assert that they are a natural outcome
of galaxy formation in a $\Lambda$CDM universe.  We predict that the
globular cluster systems of other $L^{\ast}$ galaxies in isolation
or in low-density groups like the Local Group will lie in a plane in
age–metallicity–specific orbital energy space.

We thank Brendan Griffen and David Nataf for insightful suggestions that
improved our analyses.  This work has made use of data from the European
Space Agency (ESA) mission Gaia (<https://www.cosmos.esa.int/gaia>),
processed by the Gaia Data Processing and Analysis Consortium (DPAC,
<https://www.cosmos.esa.int/web/gaia/dpac/consortium>).  Funding for
the DPAC has been provided by national institutions, in particular the
institutions participating in the Gaia Multilateral Agreement.  This
research has made use of NASA's Astrophysics Data System Bibliographic
Services.  This research has made use of the SIMBAD database, operated
at CDS, Strasbourg, France [Wenger et al., 2000].  This research has
made use of the VizieR catalogue access tool, CDS, Strasbourg, France.
The original description of the VizieR service was published in A&AS 143,
23 [Ochsenbein et al., 2000].

[Bovy, 2015],
[van der Walt et al., 2011],
[McKinney, 2010],
[Seabold & Perktold, 2010]

[Arenou et al., 2018]
Arenou, F., Luri, X., Babusiaux, C., et al. 2018, , 616, A17

[Ashman & Zepf, 1992]
Ashman, K. M., & Zepf, S. E. 1992, , 384, 50

[Baumgardt et al., 2019]
Baumgardt, H., Hilker, M., Sollima, A., & Bellini, A. 2019, ,
482, 5138

[Bertelli et al., 1994]
Bertelli, G., Bressan, A., Chiosi, C., Fagotto, F., & Nasi, E. 1994,
, 106, 275

[Binney & Tremaine, 2008]
Binney, J., & Tremaine, S. 2008, Galactic Dynamics: Second Edition

[Bjork & Chaboyer, 2006]
Bjork, S. R., & Chaboyer, B. 2006, , 641, 1102

[Bland-Hawthorn & Gerhard, 2016]
Bland-Hawthorn, J., & Gerhard, O. 2016, , 54, 529

[Bovy, 2015]
Bovy, J. 2015, , 216, 29

[Buonanno et al., 1999]
Buonanno, R., Corsi, C. E., Castellani, M., et al. 1999, , 118, 1671

[Buonanno et al., 1998]
Buonanno, R., Corsi, C. E., Zinn, R., et al. 1998, , 501, L33

[Carretta et al., 2014]
Carretta, E., Bragaglia, A., Gratton, R. G., et al. 2014, , 561,

[Carretta & Gratton, 1997]
Carretta, E., & Gratton, R. G. 1997, , 121, 95

[Carretta et al., 2010]
Carretta, E., Bragaglia, A., Gratton, R. G., et al. 2010, , 520,

[Chaboyer et al., 2001]
Chaboyer, B., Fenton, W. H., Nelan, J. E., Patnaude, D. J., & Simon,
F. E. 2001, , 562, 521

[Choksi et al., 2018]
Choksi, N., Gnedin, O. Y., & Li, H. 2018, , 480, 2343

[Chou et al., 2007]
Chou, M.-Y., Majewski, S. R., Cunha, K., et al. 2007, , 670, 346

[Cohen, 2004]
Cohen, J. G. 2004, , 127, 1545

[Crowley et al., 2016]
Crowley, C., Kohley, R., Hambly, N. C., et al. 2016, , 595, A6

[Dotter et al., 2007]
Dotter, A., Chaboyer, B., Jevremović, D., et al. 2007, , 134,

[Eadie & Jurić, 2019]
Eadie, G., & Jurić, M. 2019, The Astrophysical Journal, 875, 159

[Erkal et al., 2019]
Erkal, D., Belokurov, V., Laporte, C. F. P., et al. 2019, , 487,

[Fabricius et al., 2016]
Fabricius, C., Bastian, U., Portell, J., et al. 2016, , 595, A3

[Forbes, 2020]
Forbes, D. A. 2020, , 493, 847

[Forbes & Bridges, 2010]
Forbes, D. A., & Bridges, T. 2010, , 404, 1203

[Forbes et al., 2018]
Forbes, D. A., Bastian, N., Gieles, M., et al. 2018, Proceedings of the
Royal Society of London Series A, 474, 20170616

[Gaia Collaboration et al., 2016]
Gaia Collaboration, Prusti, T., de Bruijne, J. H. J., et al. 2016,
, 595, A1

[Gaia Collaboration et al., 2018]
Gaia Collaboration, Helmi, A., van Leeuwen, F., et al.
2018a, , 616, A12

[Gaia Collaboration et al., 2018]
Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al.
2018b, , 616, A1

[Girardi et al., 2000]
Girardi, L., Bressan, A., Bertelli, G., & Chiosi, C. 2000, , 141,

[Gravity Collaboration et al., 2018]
Gravity Collaboration, Abuter, R., Amorim, A., et al. 2018, , 615,

[Griffen et al., 2010]
Griffen, B. F., Drinkwater, M. J., Thomas, P. A., Helly, J. C., &
Pimbblet, K. A. 2010, , 405, 375

[Hambly et al., 2018]
Hambly, N. C., Cropper, M., Boudreault, S., et al. 2018, , 616, A15

[Harris, 1996]
Harris, W. E. 1996, , 112, 1487

[Hasselquist et al., 2017]
Hasselquist, S., Shetrone, M., Smith, V., et al. 2017, , 845, 162

[Hughes et al., 2019]
Hughes, M. E., Pfeffer, J., Martig, M., et al. 2019, , 482, 2795

[Jurić et al., 2008]
Jurić, M., Ivezić, Ž., Brooks, A., et al. 2008, ,
673, 864

[Keller et al., 2020]
Keller, B. W., Kruijssen, J. M. D., Pfeffer, J., et al. 2020, ,
495, 4248

[Kim et al., 2018]
Kim, J.-h., Ma, X., Grudić, M. Y., et al. 2018, , 474, 4232

[Kirby et al., 2013]
Kirby, E. N., Cohen, J. G., Guhathakurta, P., et al. 2013, , 779,

[Kravtsov & Gnedin, 2005]
Kravtsov, A. V., & Gnedin, O. Y. 2005, , 623, 650

[Kruijssen, 2015]
Kruijssen, J. M. D. 2015, , 454, 1658

[Kruijssen, 2019]
—. 2019, , 486, L20

[Kruijssen et al., 2019]
Kruijssen, J. M. D., Pfeffer, J. L., Crain, R. A., & Bastian, N.
2019a, , 486, 3134

[Kruijssen et al., 2019]
Kruijssen, J. M. D., Pfeffer, J. L., Reina-Campos, M., Crain, R. A., &
Bastian, N. 2019b, , 486, 3180

[Kruijssen et al., 2020]
Kruijssen, J. M. D., Pfeffer, J. L., Chevance, M., et al. 2020, ,
498, 2472

[Larsen et al., 2012]
Larsen, S. S., Brodie, J. P., & Strader, J. 2012, , 546, A53

[Law & Majewski, 2010]
Law, D. R., & Majewski, S. R. 2010, , 714, 229

[Leaman et al., 2013]
Leaman, R., VandenBerg, D. A., & Mendel, J. T. 2013, , 436, 122

[Letarte et al., 2006]
Letarte, B., Hill, V., Jablonka, P., et al. 2006, , 453, 547

[Li & Gnedin, 2014]
Li, H., & Gnedin, O. Y. 2014, , 796, 10

[Li & Gnedin, 2019]
—. 2019, , 486, 4030

[Lindegren et al., 2018]
Lindegren, L., Hernández, J., Bombrun, A., et al. 2018, , 616,

[Luri et al., 2018]
Luri, X., Brown, A. G. A., Sarro, L. M., et al. 2018, , 616, A9

[Mackey & Gilmore, 2004]
Mackey, A. D., & Gilmore, G. F. 2004, , 355, 504

[Marín-Franch et al., 2009]
Marín-Franch, A., Aparicio, A., Piotto, G., et al. 2009, ,
694, 1498

[Massari et al., 2019]
Massari, D., Koppelman, H. H., & Helmi, A. 2019, , 630, L4

[McKinney, 2010]
McKinney, W. 2010, in Proceedings of the 9th Python in Science Conference

[McMillan, 2017]
McMillan, P. J. 2017, , 465, 76

[Miyamoto & Nagai, 1975]
Miyamoto, M., & Nagai, R. 1975, , 27, 533

[Mo et al., 2010]
Mo, H., van den Bosch, F. C., & White, S. 2010, Galaxy Formation and

[Mottini et al., 2008]
Mottini, M., Wallerstein, G., & McWilliam, A. 2008, , 136, 614

[Muratov & Gnedin, 2010]
Muratov, A. L., & Gnedin, O. Y. 2010, , 718, 1266

[Myeong et al., 2018]
Myeong, G. C., Evans, N. W., Belokurov, V., Sanders, J. L., &
Koposov, S. E. 2018, , 863, L28

[Navarro et al., 1996]
Navarro, J. F., Frenk, C. S., & White, S. D. M. 1996, , 462, 563

[Ochsenbein et al., 2000]
Ochsenbein, F., Bauer, P., & Marcout, J. 2000, , 143, 23

[Pascale et al., 2018]
Pascale, R., Posti, L., Nipoti, C., & Binney, J. 2018, , 480,

[Patel et al., 2013]
Patel, S. G., van Dokkum, P. G., Franx, M., et al. 2013a,
, 766, 15

[Patel et al., 2013]
Patel, S. G., Fumagalli, M., Franx, M., et al. 2013b,
, 778, 115

[Pfeffer et al., 2018]
Pfeffer, J., Kruijssen, J. M. D., Crain, R. A., & Bastian, N. 2018,
, 475, 4309

[Pfeffer et al., 2020]
Pfeffer, J. L., Trujillo-Gomez, S., Kruijssen, J. M. D., et al. 2020,
, 499, 4863

[Pietrinferni et al., 2004]
Pietrinferni, A., Cassisi, S., Salaris, M., & Castelli, F. 2004, ,
612, 168

[Renaud et al., 2017]
Renaud, F., Agertz, O., & Gieles, M. 2017, , 465, 3622

[Rutledge et al., 1997]
Rutledge, G. A., Hesser, J. E., & Stetson, P. B. 1997a,
, 109, 907

[Rutledge et al., 1997]
Rutledge, G. A., Hesser, J. E., Stetson, P. B., et al.
1997b, , 109, 883

[Sbordone et al., 2007]
Sbordone, L., Bonifacio, P., Buonanno, R., et al. 2007, , 465, 815

[Schönrich et al., 2010]
Schönrich, R., Binney, J., & Dehnen, W. 2010, , 403, 1829

[Seabold & Perktold, 2010]
Seabold, S., & Perktold, J. 2010, in 9th Python in Science Conference

[Searle & Zinn, 1978]
Searle, L., & Zinn, R. 1978, , 225, 357

[Simon, 2019]
Simon, J. D. 2019, , 57, 375

[Sohn et al., 2018]
Sohn, S. T., Watkins, L. L., Fardal, M. A., et al. 2018, , 862, 52

[Trujillo-Gomez et al., 2021]
Trujillo-Gomez, S., Kruijssen, J. M. D., Reina-Campos, M., et al. 2021,
, 503, 31

[van der Walt et al., 2011]
van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, Computing in
Science and Engineering, 13, 22

[VandenBerg et al., 2013]
VandenBerg, D. A., Brogaard, K., Leaman, R., & Casagrand e, L. 2013,
, 775, 134

[Vasiliev & Belokurov, 2020]
Vasiliev, E., & Belokurov, V. 2020, , 497, 4162

[Wagner-Kaiser et al., 2017]
Wagner-Kaiser, R., Mackey, D., Sarajedini, A., et al. 2017, ,
471, 3347

[Wang et al., 2020]
Wang, W., Han, J., Cautun, M., Li, Z., & Ishigaki, M. N. 2020,
Science China Physics, Mechanics, and Astronomy, 63, 109801

[Wenger et al., 2000]
Wenger, M., Ochsenbein, F., Egret, D., et al. 2000, , 143, 9

[West et al., 2004]
West, M. J., Côté, P., Marzke, R. O., & Jordán, A. 2004,
, 427, 31

[Wright, 2006]
Wright, E. L. 2006, , 118, 1711

[Zinn & West, 1984]
Zinn, R., & West, M. J. 1984, , 55, 45

Input Astrometric, Distance, and Radial Velocity Data

[mas yr$^{-1}$]
[mas yr$^{-1}$]

[km s$^{-1}$]
NGC 104     47 Tuc     6.0194     -72.0821     5.2477^+0.0016_-0.0016      -2.5189^+0.0015_-0.0015     0.1959^+0.0002_-0.0002     -0.06     -0.01     -0.01     4.3^+0.1_-0.1     -18.7^+0.2_-0.2

NGC 288     $\ldots$     13.1879     -26.5858     4.2385^+0.0035_-0.0035      -5.6470^+0.0026_-0.0026     0.1401^+0.0021_-0.0021     0.25     0.15     -0.13     8.1^+0.2_-0.2     -46.6^+0.4_-0.4

NGC 362     $\ldots$     15.8099     -70.8489     6.6954^+0.0045_-0.0045      -2.5184^+0.0034_-0.0034     0.0788^+0.0012_-0.0012     -0.09     -0.04     -0.12     8.3^+0.2_-0.2     223.5^+0.5_-0.5

NGC 1261     $\ldots$     48.067543     -55.216224     1.6900^+0.0400_-0.0400      -2.1100^+0.0400_-0.0400     . $\ldots$     . $\ldots$     . $\ldots$     . $\ldots$     16.0^+0.4_-0.4     68.2^+4.6_-4.6

NGC 1851     $\ldots$     78.528     -40.0456     2.1308^+0.0037_-0.0037      -0.6220^+0.0040_-0.0040     0.0298^+0.0011_-0.0011     -0.09     0.06     -0.07     12.2^+0.3_-0.3     320.9^+1.0_-1.0

NGC 2298     $\ldots$     102.2464     -36.0046     3.2762^+0.0060_-0.0060      -2.1913^+0.0061_-0.0061     0.0791^+0.0019_-0.0019     0.07     0.08     -0.07     10.6^+0.2_-0.2     148.9^+1.2_-1.2

NGC 2808     $\ldots$     138.0071     -64.8645     1.0032^+0.0032_-0.0032      0.2785^+0.0032_-0.0032     0.0560^+0.0006_-0.0006     -0.08     0.05     -0.01     9.3^+0.2_-0.2     93.6^+2.4_-2.4

NGC 3201     $\ldots$     154.3987     -46.4125     8.3344^+0.0021_-0.0021      -1.9895^+0.0020_-0.0020     0.1724^+0.0006_-0.0006     0.12     0.04     -0.02     5.1^+0.1_-0.1     494.0^+0.2_-0.2

NGC 4147     $\ldots$     182.52626     18.542638     -1.7800^+0.0400_-0.0400      -2.1000^+0.0400_-0.0400     . $\ldots$     . $\ldots$     . $\ldots$     . $\ldots$     18.8^+0.4_-0.4     183.2^+0.7_-0.7

NGC 4590     M 68     189.8651     -26.7454     -2.7640^+0.0050_-0.0050      1.7916^+0.0039_-0.0039     0.0664^+0.0025_-0.0025     -0.29     -0.0     0.13     10.1^+0.2_-0.2     -95.2^+0.4_-0.4

NGC 4833     $\ldots$     194.8978     -70.8718     -8.3147^+0.0036_-0.0036      -0.9366^+0.0029_-0.0029     0.1163^+0.0010_-0.0010     0.06     0.05     0.11     5.9^+0.1_-0.1     200.2^+1.2_-1.2

NGC 5024     M 53     198.2262     18.1661     -0.1466^+0.0045_-0.0045      -1.3514^+0.0032_-0.0032     0.0143^+0.0018_-0.0018     -0.28     -0.12     0.08     18.4^+0.4_-0.4     -79.1^+4.1_-4.1

NGC 5053     $\ldots$     199.1124     17.7008     -0.3591^+0.0071_-0.0071      -1.2586^+0.0048_-0.0048     0.0064^+0.0040_-0.0040     -0.32     0.13     -0.17     16.2^+0.4_-0.4     44.0^+0.4_-0.4

NGC 5139     $\omega$ Cen     201.7876     -47.4515     -3.1925^+0.0022_-0.0022      -6.7445^+0.0019_-0.0019     0.1237^+0.0011_-0.0011     -0.03     -0.04     0.17     5.1^+0.1_-0.1     232.3^+0.5_-0.5

NGC 5272     M 3     205.5486     28.376     -0.1127^+0.0029_-0.0029      -2.6274^+0.0022_-0.0022     0.0265^+0.0010_-0.0010     -0.03     -0.01     -0.05     10.0^+0.2_-0.2     -148.6^+0.4_-0.4

NGC 5286     $\ldots$     206.6136     -51.3723     0.1836^+0.0076_-0.0076      -0.1477^+0.0068_-0.0068     0.0168^+0.0025_-0.0025     -0.01     -0.02     0.08     10.7^+0.2_-0.2     58.3^+3.2_-3.2

NGC 5466     $\ldots$     211.3614     28.5331     -5.4044^+0.0042_-0.0042      -0.7907^+0.0041_-0.0041     0.0210^+0.0021_-0.0021     0.07     0.04     0.15     16.6^+0.4_-0.4     107.7^+0.3_-0.3

NGC 5904     M 5     229.6394     2.0766     4.0613^+0.0032_-0.0032      -9.8610^+0.0029_-0.0029     0.1135^+0.0010_-0.0010     0.03     -0.07     0.09     7.3^+0.2_-0.2     52.1^+0.5_-0.5

NGC 5927     $\ldots$     232.0065     -50.6694     -5.0470^+0.0060_-0.0060      -3.2325^+0.0055_-0.0055     0.0996^+0.0021_-0.0021     -0.08     -0.0     -0.01     7.4^+0.2_-0.2     -115.7^+3.1_-3.1

NGC 5986     $\ldots$     236.5211     -37.7826     -4.2217^+0.0084_-0.0084      -4.5515^+0.0065_-0.0065     0.0718^+0.0031_-0.0031     -0.18     -0.1     -0.01     10.3^+0.2_-0.2     88.9^+3.7_-3.7

NGC 6093     M 80     244.2564     -22.9723     -2.9469^+0.0090_-0.0090      -5.5613^+0.0073_-0.0073     0.0558^+0.0030_-0.0030     0.01     -0.1     0.06     8.7^+0.2_-0.2     9.3^+3.1_-3.1

NGC 6121     M 4     245.8976     -26.5279     -12.4956^+0.0033_-0.0033      -18.9789^+0.0030_-0.0030     0.5001^+0.0007_-0.0007     0.08     -0.08     0.03     2.2^+0.1_-0.1     70.2^+0.3_-0.3

NGC 6101     $\ldots$     246.450485     -72.202194     1.7300^+0.0500_-0.0500      -0.4700^+0.0500_-0.0500     . $\ldots$     . $\ldots$     . $\ldots$     . $\ldots$     15.1^+0.3_-0.4     361.4^+1.7_-1.7

NGC 6144     $\ldots$     246.8061     -26.0301     -1.7646^+0.0085_-0.0085      -2.6371^+0.0063_-0.0063     0.0668^+0.0040_-0.0040     0.08     -0.15     0.15     10.1^+0.2_-0.2     188.9^+1.1_-1.1

NGC 6171     M 107     248.135     -13.057     -1.9359^+0.0064_-0.0064      -5.9487^+0.0048_-0.0048     0.1480^+0.0026_-0.0026     0.03     -0.13     0.17     6.3^+0.1_-0.1     -33.8^+0.3_-0.3

NGC 6205     M 13     250.4217     36.4596     -3.1762^+0.0027_-0.0027      -2.5876^+0.0030_-0.0030     0.0801^+0.0007_-0.0007     0.19     -0.04     0.04     7.0^+0.2_-0.2     -246.6^+0.9_-0.9

NGC 6218     M 12     251.8101     -1.951     -0.1577^+0.0040_-0.0040      -6.7683^+0.0027_-0.0027     0.1563^+0.0013_-0.0013     0.29     -0.06     0.11     4.7^+0.1_-0.1     -42.1^+0.6_-0.6

NGC 6254     M 10     254.2861     -4.0981     -4.7031^+0.0039_-0.0039      -6.5285^+0.0027_-0.0027     0.1511^+0.0014_-0.0014     0.21     -0.07     0.1     4.3^+0.1_-0.1     75.8^+1.0_-1.0

NGC 6304     $\ldots$     258.637     -29.4816     -3.9478^+0.0095_-0.0095      -1.1248^+0.0069_-0.0069     0.1077^+0.0034_-0.0034     0.16     -0.04     0.13     6.0^+0.1_-0.1     -107.3^+3.6_-3.6

NGC 6341     M 92     259.2821     43.1352     -4.9367^+0.0040_-0.0040      -0.5559^+0.0040_-0.0040     0.0564^+0.0008_-0.0008     0.11     -0.04     -0.0     8.1^+0.2_-0.2     -121.6^+1.5_-1.5

NGC 6352     $\ldots$     261.3739     -48.427     -2.1889^+0.0046_-0.0046      -4.4209^+0.0036_-0.0036     0.1543^+0.0018_-0.0018     0.22     -0.12     0.13     5.6^+0.1_-0.1     -120.9^+3.0_-3.0

NGC 6366     $\ldots$     261.9393     -5.0752     -0.3835^+0.0054_-0.0054      -5.1309^+0.0044_-0.0044     0.2292^+0.0022_-0.0022     0.27     -0.08     0.14     3.6^+0.1_-0.1     -122.3^+0.5_-0.5

NGC 6362     $\ldots$     262.9772     -67.0492     -5.5014^+0.0028_-0.0028      -4.7417^+0.0032_-0.0032     0.0974^+0.0011_-0.0011     0.06     -0.06     0.06     7.5^+0.2_-0.2     -13.1^+0.6_-0.6

NGC 6388     $\ldots$     264.0654     -44.7423     -1.3548^+0.0072_-0.0072      -2.7144^+0.0061_-0.0061     0.0482^+0.0034_-0.0034     0.11     -0.12     0.18     11.5^+0.3_-0.3     81.2^+1.2_-1.2

NGC 6397     $\ldots$     265.1697     -53.6773     3.2908^+0.0026_-0.0026      -17.5908^+0.0025_-0.0025     0.3781^+0.0007_-0.0007     0.1     -0.05     0.07     2.2^+0.1_-0.1     18.9^+0.1_-0.1

NGC 6441     $\ldots$     267.5540     -37.066     -2.5394^+0.0070_-0.0070      -5.3010^+0.0057_-0.0057     0.0403^+0.0037_-0.0037     0.17     -0.18     0.17     9.7^+0.2_-0.2     18.3^+1.7_-1.7

NGC 6496     $\ldots$     269.7677     -44.2660     -3.0290^+0.0057_-0.0057      -9.1971^+0.0050_-0.0050     0.0803^+0.0031_-0.0031     0.07     -0.15     0.17     11.6^+0.3_-0.3     -112.7^+5.7_-5.7

NGC 6535     $\ldots$     270.959     -0.2953     -4.2101^+0.0115_-0.0115      -2.9461^+0.0108_-0.0108     0.1294^+0.0047_-0.0047     0.1     0.1     -0.02     6.8^+0.2_-0.2     -215.1^+0.5_-0.5

NGC 6541     $\ldots$     271.9827     -43.7144     0.2762^+0.0054_-0.0054      -8.7659^+0.0048_-0.0048     0.1139^+0.0025_-0.0025     -0.03     -0.04     0.08     7.4^+0.2_-0.2     -156.2^+2.7_-2.7

NGC 6584     $\ldots$     274.656646     -52.215778     -0.0500^+0.0100_-0.0100      -7.2200^+0.0100_-0.0100     . $\ldots$     -0.15     . $\ldots$     . $\ldots$     13.0^+0.3_-0.3     222.9^+15.0_-15.0

NGC 6637     M 69     277.8342     -32.3565     -5.0669^+0.0104_-0.0104      -5.8017^+0.0094_-0.0094     0.0746^+0.0032_-0.0032     0.23     -0.01     -0.05     8.2^+0.2_-0.2     39.9^+2.8_-2.8

NGC 6656     M 22     279.1048     -23.9102     9.8019^+0.0036_-0.0036      -5.5643^+0.0034_-0.0034     0.2602^+0.0009_-0.0009     0.15     -0.05     0.02     3.2^+0.1_-0.1     -148.9^+0.4_-0.4

NGC 6681     M 70     280.802     -32.2892     1.3853^+0.0076_-0.0076      -4.7174^+0.0065_-0.0065     0.1096^+0.0038_-0.0038     0.26     -0.24     0.1     8.7^+0.2_-0.2     218.6^+1.2_-1.2

NGC 6723     $\ldots$     284.888123     -36.632248     1.0000^+0.0100_-0.0100      -2.4200^+0.0100_-0.0100     . $\ldots$     0.11     . $\ldots$     . $\ldots$     8.6^+0.2_-0.2     -94.5^+3.6_-3.6

NGC 6752     $\ldots$     287.7175     -59.9833     -3.1908^+0.0018_-0.0018      -4.0347^+0.0020_-0.0020     0.2310^+0.0011_-0.0011     0.19     -0.29     0.03     3.9^+0.1_-0.1     -24.5^+1.9_-1.9

NGC 6779     M 56     289.148     30.184     -2.0092^+0.0051_-0.0051      1.6553^+0.0056_-0.0056     0.0702^+0.0015_-0.0015     0.03     -0.05     -0.03     9.9^+0.2_-0.2     -135.7^+0.8_-0.8

Terzan 7     $\ldots$     289.432983     -34.657722     -3.0400^+0.0600_-0.0600      -1.7100^+0.0500_-0.0500     . $\ldots$     . $\ldots$     . $\ldots$     . $\ldots$     23.0^+0.5_-0.5     166.0^+4.0_-4.0

Arp 2     $\ldots$     292.183807     -30.355638     -2.4000^+0.0400_-0.0400      -1.5400^+0.0300_-0.0300     . $\ldots$     . $\ldots$     . $\ldots$     . $\ldots$     27.6^+0.6_-0.6     115.0^+10.0_-10.0

NGC 6809     M 55     295.0046     -30.9621     -3.4017^+0.0031_-0.0031      -9.2642^+0.0028_-0.0028     0.1707^+0.0011_-0.0011     0.18     -0.03     0.0     5.3^+0.1_-0.1     174.8^+0.4_-0.4

Terzan 8     $\ldots$     295.435028     -33.999474     -2.9100^+0.0800_-0.0800      -1.6300^+0.0600_-0.0600     . $\ldots$     . $\ldots$     . $\ldots$     . $\ldots$     25.4^+0.6_-0.6     130.0^+8.0_-8.0

NGC 6838     M 71     298.4427     18.779     -3.3842^+0.0027_-0.0027      -2.6528^+0.0028_-0.0028     0.2252^+0.0010_-0.0010     0.11     -0.11     -0.01     3.8^+0.1_-0.1     -22.9^+0.2_-0.2

NGC 6934     $\ldots$     308.547393     7.404472     -2.6700^+0.0400_-0.0400      -4.5200^+0.0500_-0.0500     . $\ldots$     . $\ldots$     . $\ldots$     . $\ldots$     15.2^+0.3_-0.4     -411.4^+1.6_-1.6

NGC 6981     M 72     313.3662     -12.5386     -1.2488^+0.0089_-0.0089      -3.3117^+0.0068_-0.0068     0.0225^+0.0063_-0.0063     0.26     -0.38     -0.13     16.8^+0.4_-0.4     -345.1^+3.7_-3.7

NGC 7078     M 15     322.4949     12.1661     -0.6238^+0.0041_-0.0041      -3.7960^+0.0039_-0.0039     0.0568^+0.0014_-0.0014     -0.04     -0.02     -0.15     10.2^+0.2_-0.2     -107.5^+0.3_-0.3

NGC 7089     M 2     323.3497     -0.8177     3.4911^+0.0077_-0.0077      -2.1501^+0.0071_-0.0071     0.0591^+0.0035_-0.0035     -0.04     -0.14     -0.14     11.4^+0.3_-0.3     -5.3^+2.0_-2.0

NGC 7099     M 30     325.0888     -23.1792     -0.7017^+0.0063_-0.0063      -7.2218^+0.0055_-0.0055     0.0746^+0.0040_-0.0040     0.3     -0.29     -0.27     7.9^+0.2_-0.2     -184.2^+1.0_-1.0

Pal 12     $\ldots$     326.661804     -21.252611     -3.0600^+0.0500_-0.0500      -3.3200^+0.0500_-0.0500     . $\ldots$     . $\ldots$     . $\ldots$     . $\ldots$     18.7^+0.4_-0.4     27.8^+1.5_-1.5
Input globular cluster astrometry, distance,
and radial velocity data.  While most cluster astrometric data
come from Gaia Collaboration et al., 2018, the data for NGC 6584 and NGC 6723
come from Baumgardt et al., 2019 while the data for NGC 1261, NGC
4147, NGC 6101, Terzan 7, Arp 2, Terzan 8, NGC 6934, and Pal 12
come from Sohn et al., 2018.  Distance and radial velocity data
come from the December 2010 revision of the Harris, 1996

Age, Metallicity, and Specific Orbital Energy Inferences

NGC 104     47 Tuc     1.02^+0.07_-0.07     -0.64^+0.05_-0.05     -7.64^+0.06_-0.05     -11.82^+0.05_-0.05     -16.62^+0.05_-0.05

NGC 288     $\ldots$     0.83^+0.03_-0.03     -0.92^+0.05_-0.05     -9.04^+0.08_-0.04     -13.06^+0.07_-0.05     -16.00^+0.11_-0.07

NGC 362     $\ldots$     0.81^+0.04_-0.04     -0.87^+0.05_-0.05     -8.66^+0.16_-0.18     -12.84^+0.17_-0.17     -16.20^+0.17_-0.20

NGC 1261     $\ldots$     0.80^+0.04_-0.04     -0.86^+0.05_-0.05     -8.20^+0.07_-0.06     -12.04^+0.06_-0.08     -13.43^+0.13_-0.09

NGC 1851     $\ldots$     0.78^+0.04_-0.04     -0.81^+0.05_-0.05     -7.57^+0.15_-0.10     -11.34^+0.14_-0.09     -13.47^+0.14_-0.09

NGC 2298     $\ldots$     0.99^+0.05_-0.05     -1.49^+0.05_-0.05     -8.04^+0.19_-0.10     -11.85^+0.17_-0.13     -14.15^+0.22_-0.15

NGC 2808     $\ldots$     0.85^+0.03_-0.03     -0.89^+0.05_-0.05     -8.12^+0.08_-0.07     -12.24^+0.13_-0.08     -15.33^+0.11_-0.10

NGC 3201     $\ldots$     0.80^+0.03_-0.03     -1.02^+0.05_-0.05     -3.13^+0.19_-0.25     -7.21^+0.19_-0.25     -11.01^+0.23_-0.26

NGC 4147     $\ldots$     0.89^+0.04_-0.04     -1.28^+0.05_-0.05     -7.34^+0.04_-0.04     -11.08^+0.06_-0.04     -12.33^+0.11_-0.09

NGC 4590     M 68     0.90^+0.04_-0.04     -1.78^+0.05_-0.05     -3.94^+0.56_-0.59     -8.17^+0.71_-0.56     -11.41^+0.51_-0.55

NGC 4833     $\ldots$     0.98^+0.05_-0.05     -1.49^+0.05_-0.05     -9.37^+0.10_-0.09     -13.58^+0.08_-0.08     -18.66^+0.12_-0.11

NGC 5024     M 53     0.99^+0.05_-0.05     -1.64^+0.05_-0.05     -7.44^+0.88_-0.53     -11.04^+0.97_-0.70     -12.26^+0.90_-0.53

NGC 5053     $\ldots$     0.96^+0.04_-0.04     -1.76^+0.05_-0.05     -7.75^+0.73_-0.60     -11.68^+0.92_-0.60     -12.97^+0.75_-0.60

NGC 5139      $\omega$ Cen     0.90^+0.05_-0.05     -1.13^+0.05_-0.05     -9.29^+0.11_-0.06     -13.51^+0.07_-0.08     -18.60^+0.12_-0.07

NGC 5272     M 3     0.89^+0.04_-0.04     -1.12^+0.05_-0.05     -7.55^+0.13_-0.12     -11.61^+0.13_-0.12     -14.06^+0.13_-0.12

NGC 5286     $\ldots$     0.98^+0.04_-0.04     -1.19^+0.05_-0.05     -8.06^+0.07_-0.08     -12.48^+0.11_-0.07     -15.88^+0.16_-0.13

NGC 5466     $\ldots$     1.06^+0.05_-0.05     -1.98^+0.05_-0.05     -3.35^+0.65_-0.52     -7.27^+0.44_-0.53     -9.00^+0.73_-0.46

NGC 5904     M 5     0.83^+0.03_-0.03     -0.90^+0.05_-0.05     -4.44^+0.59_-0.39     -8.76^+0.50_-0.41     -12.72^+0.68_-0.43

NGC 5927     $\ldots$     0.99^+0.07_-0.07     -0.50^+0.05_-0.05     -7.18^+0.09_-0.10     -11.62^+0.10_-0.09     -18.37^+0.10_-0.07

NGC 5986     $\ldots$     0.95^+0.04_-0.04     -1.13^+0.05_-0.05     -10.81^+0.19_-0.16     -15.47^+0.21_-0.12     -20.28^+0.29_-0.22

NGC 6093     M 80     0.98^+0.04_-0.04     -1.25^+0.05_-0.05     -10.84^+0.07_-0.05     -15.42^+0.08_-0.07     -21.78^+0.17_-0.12

NGC 6121     M 4     0.98^+0.05_-0.05     -0.83^+0.05_-0.05     -9.72^+0.02_-0.01     -13.93^+0.02_-0.01     -19.69^+0.05_-0.04

NGC 6101     $\ldots$     0.98^+0.04_-0.04     -1.54^+0.05_-0.05     -5.04^+0.13_-0.11     -9.60^+0.13_-0.11     -10.27^+0.22_-0.26

NGC 6144     $\ldots$     1.08^+0.05_-0.05     -1.52^+0.05_-0.05     -9.58^+0.18_-0.18     -14.31^+0.15_-0.13     -19.82^+0.18_-0.20

NGC 6171     M 107     1.09^+0.06_-0.06     -0.81^+0.05_-0.05     -10.00^+0.09_-0.09     -14.44^+0.08_-0.07     -21.27^+0.11_-0.08

NGC 6205     M 13     0.91^+0.04_-0.04     -1.11^+0.05_-0.05     -9.32^+0.10_-0.07     -13.51^+0.10_-0.05     -17.34^+0.11_-0.08

NGC 6218     M 12     0.99^+0.03_-0.03     -0.92^+0.05_-0.05     -9.01^+0.15_-0.18     -13.32^+0.18_-0.19     -19.47^+0.17_-0.20

NGC 6254     M 10     0.89^+0.04_-0.04     -1.03^+0.05_-0.05     -8.82^+0.10_-0.08     -13.14^+0.11_-0.07     -19.37^+0.12_-0.08

NGC 6304     $\ldots$     1.06^+0.08_-0.08     -0.52^+0.05_-0.05     -8.08^+0.23_-0.21     -12.52^+0.22_-0.22     -22.01^+0.27_-0.34

NGC 6341     M 92     1.03^+0.04_-0.04     -1.94^+0.05_-0.05     -9.00^+0.10_-0.11     -13.16^+0.13_-0.10     -16.63^+0.09_-0.11

NGC 6352     $\ldots$     0.99^+0.07_-0.07     -0.56^+0.05_-0.05     -7.69^+0.17_-0.15     -12.06^+0.13_-0.19     -20.12^+0.18_-0.16

NGC 6366     $\ldots$     1.04^+0.13_-0.13     -0.59^+0.05_-0.05     -8.53^+0.19_-0.10     -12.78^+0.15_-0.11     -19.31^+0.20_-0.10

NGC 6362     $\ldots$     1.06^+0.05_-0.05     -0.85^+0.05_-0.05     -8.97^+0.08_-0.07     -13.36^+0.07_-0.09     -19.02^+0.09_-0.11

NGC 6388     $\ldots$     0.94^+0.08_-0.08     -0.63^+0.05_-0.05     -11.34^+0.21_-0.16     -16.13^+0.27_-0.18     -20.99^+0.26_-0.29

NGC 6397     $\ldots$     0.99^+0.04_-0.04     -1.54^+0.05_-0.05     -8.24^+0.04_-0.06     -12.44^+0.05_-0.06     -18.30^+0.04_-0.06

NGC 6441     $\ldots$     0.88^+0.07_-0.07     -0.46^+0.05_-0.05     -11.26^+0.09_-0.08     -16.03^+0.11_-0.08     -24.53^+0.40_-0.33

NGC 6496     $\ldots$     0.97^+0.07_-0.07     -0.56^+0.05_-0.05     -5.86^+0.60_-0.60     -10.94^+0.69_-0.43     -15.40^+0.77_-0.79

NGC 6535     $\ldots$     0.82^+0.09_-0.09     -1.29^+0.05_-0.05     -9.38^+0.09_-0.09     -13.84^+0.12_-0.11     -20.81^+0.08_-0.08

NGC 6541     $\ldots$     1.01^+0.04_-0.04     -1.31^+0.05_-0.05     -7.52^+0.11_-0.13     -12.05^+0.15_-0.14     -20.74^+0.12_-0.13

NGC 6584     $\ldots$     0.88^+0.03_-0.03     -1.10^+0.05_-0.05     -4.81^+0.78_-0.75     -9.54^+0.64_-0.73     -11.81^+0.92_-0.91

NGC 6637     M 69     1.02^+0.07_-0.07     -0.64^+0.05_-0.05     -10.39^+0.17_-0.17     -14.94^+0.17_-0.16     -24.16^+0.21_-0.17

NGC 6656     M 22     0.99^+0.05_-0.05     -1.27^+0.05_-0.05     -5.39^+0.15_-0.16     -9.63^+0.19_-0.19     -16.48^+0.15_-0.17

NGC 6681     M 70     1.00^+0.04_-0.04     -1.13^+0.05_-0.05     -7.31^+0.42_-0.35     -11.96^+0.45_-0.34     -19.70^+0.45_-0.34

NGC 6723     $\ldots$     1.02^+0.05_-0.05     -0.82^+0.05_-0.05     -9.40^+0.43_-0.27     -14.01^+0.42_-0.29     -21.21^+0.49_-0.35

NGC 6752     $\ldots$     0.92^+0.04_-0.04     -1.02^+0.05_-0.05     -8.18^+0.13_-0.09     -12.41^+0.08_-0.10     -18.53^+0.13_-0.10

NGC 6779     M 56     1.07^+0.05_-0.05     -1.50^+0.05_-0.05     -8.00^+0.17_-0.14     -12.27^+0.17_-0.16     -15.73^+0.18_-0.17

Terzan 7     $\ldots$     0.57^+0.04_-0.04     -0.42^+0.05_-0.05     -9.39^+0.21_-0.19     -14.83^+0.22_-0.18     -9.46^+0.34_-0.35

Arp 2     $\ldots$     0.85^+0.06_-0.06     -1.23^+0.05_-0.05     -9.59^+0.43_-0.33     -14.87^+0.35_-0.33     -9.10^+0.44_-0.36

NGC 6809     M 55     0.96^+0.04_-0.04     -1.32^+0.05_-0.05     -8.19^+0.05_-0.03     -12.56^+0.05_-0.04     -18.88^+0.05_-0.04

Terzan 8     $\ldots$     0.95^+0.04_-0.04     -1.58^+0.05_-0.05     -8.12^+0.40_-0.37     -13.22^+0.34_-0.39     -8.35^+0.56_-0.35

NGC 6838     M 71     1.07^+0.08_-0.08     -0.59^+0.05_-0.05     -7.69^+0.04_-0.05     -11.88^+0.05_-0.05     -17.25^+0.04_-0.04

NGC 6934     $\ldots$     0.87^+0.04_-0.04     -1.10^+0.05_-0.05     -4.97^+0.20_-0.24     -9.37^+0.22_-0.23     -10.43^+0.33_-0.37

NGC 6981     M 72     0.85^+0.02_-0.02     -0.99^+0.05_-0.05     -7.58^+0.67_-0.45     -12.04^+0.71_-0.45     -12.66^+0.79_-0.51

NGC 7078     M 15     1.01^+0.04_-0.04     -1.80^+0.05_-0.05     -9.24^+0.09_-0.06     -13.43^+0.09_-0.06     -16.15^+0.13_-0.13

NGC 7089     M 2     0.92^+0.04_-0.04     -1.09^+0.05_-0.05     -7.14^+0.42_-0.39     -11.54^+0.63_-0.42     -13.97^+0.38_-0.39

NGC 7099     M 30     1.01^+0.05_-0.05     -1.70^+0.05_-0.05     -9.24^+0.21_-0.12     -13.56^+0.23_-0.13     -17.53^+0.29_-0.17

Pal 12     $\ldots$     0.69^+0.09_-0.09     -0.69^+0.05_-0.05     -5.12^+0.24_-0.26     -9.48^+0.27_-0.24     -9.18^+0.34_-0.35
Table <ref> is ordered by R.A.  The first
column is normalized age $\tau = \text{age}/12.8\text{~Gyr}$ from
Marín-Franch et al., 2009 calculated using the Dartmouth Stellar Evolution
Program [Chaboyer et al., 2001, Bjork & Chaboyer, 2006, Dotter et al., 2007] assuming the
metallicities presented in Rutledge et al., 1997, Rutledge et al., 1997 on the
Carretta & Gratton, 1997 scale.  The second column metallicity [M/H] is from
Rutledge et al., 1997, Rutledge et al., 1997 on the Carretta & Gratton, 1997 scale.
The third, fourth, and fifth columns are globular cluster specific orbital
energies calculated using the data in Table <ref> assuming the
default, , the scaled
[Bovy, 2015], and the  [McMillan, 2017]
potentials respectively.