# On the Origin of High Velocity Clouds in the Galaxy

Scott Lucchini Center for Astrophysics $|$ Harvard & Smithsonian, 60 Garden
Street, Cambridge, MA 02138, USA Scott Lucchini<EMAIL_ADDRESS>Jiwon Jesse Han Center for Astrophysics $|$ Harvard & Smithsonian, 60 Garden
Street, Cambridge, MA 02138, USA  Lars Hernquist Center for Astrophysics $|$
Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA  Charlie
Conroy Center for Astrophysics $|$ Harvard & Smithsonian, 60 Garden Street,
Cambridge, MA 02138, USA

###### Abstract

The origin of our Galaxy’s high-velocity clouds (HVCs) remains a mystery after
many decades of effort. In this paper, we use the TNG50 simulation of the
IllustrisTNG project to identify cool, dense clouds that match observations of
Galactic HVCs. We track these clouds back in time to determine their origin.
We find that only 17% of HVCs (by mass) can be tracked directly to the disk,
and 21% to satellite accretion. The majority of HVCs (62%) arise from warm and
hot circumgalactic gas that cools through thermal instability. They then
obtain their anomalous velocities through interactions with the turbulent
circumgalactic medium. At TNG50 resolution, we do not see evidence for HVCs
forming out of very low metallicity intergalactic material. Instead, low
metallicity HVCs are most likely associated with satellite accretion. These
results suggest that Galactic HVCs are highly heterogeneous in their origin,
and can provide insight into the physical processes that shape the
circumgalactic medium such as disk outflows, satellite accretion, and thermal
instabilities.

High-velocity clouds (735); Milky Way Galaxy (1054); Astronomical simulations
(1857); Circumgalactic medium (1879)

##  1 Introduction

The circumgalactic medium (CGM) plays a critical role in galaxy evolution – it
is the interface between a galaxy and its larger cosmological environment. The
CGM also contains a substantial reservoir of gas that can continue to feed the
galaxy and maintain star formation over Gyr timescales (Spitzer, 1956; Binney,
1977; Chevalier & Oegerle, 1979; Tumlinson et al., 2017). This star formation
and gas accretion subsequently induces strong feedback from supernovae and
active galactic nuclei returning gas and energy to the CGM and beyond.
However, this process of accretion and outflow – the “baryon cycle” – remains
poorly understood in detail. One of the most significant puzzles is how we can
connect our theoretical models and simulations back to observations of the CGM
(Peeples et al., 2019a; Faucher-Giguère & Oh, 2023).

Some of the earliest observations of circumgalactic gas come from our own
Galaxy in the form of “high-velocity clouds” (HVCs; van Woerden et al. 1957;
Muller et al. 1963; Wakker & van Woerden 1997; Putman et al. 2012). As
observations have improved we have uncovered more information about these
elusive clouds, but our understanding of their properties is still limited.
Modern H I 21-cm emission surveys give us an excellent view of their warm,
neutral gas structure and morphology (Kalberla & Haud, 2006; Moss et al.,
2013; Westmeier, 2018), and absorption spectroscopy can give us information
about their ionization state, composition, and distance (e.g., Smoker et al.,
2011; Fox et al., 2016; Lehner et al., 2022). However, these sightlines are
not as numerous as necessary to obtain strong constraints on the global
properties of HVCs.

A fundamental mystery surrounding the HVC population of the Milky Way (MW) is
its origins. Many models have been proposed, including satellite accretion
(Olano, 2008; Richter et al., 2018), galactic fountain recycling (Habe &
Ikeuchi, 1980; Bregman, 1980; Fox et al., 2023), accretion from the
intergalactic medium (Kereš & Hernquist, 2009), or the precipitation of small
cloudlets (Field, 1965; Maller & Bullock, 2004; Voit, 2021; Tripp, 2022).
While any single model is unlikely to explain all of the HVCs in the Galaxy,
it is important to understand the relative importance of each mechanism at
play in the Galactic CGM. Cosmological hydrodynamical simulations pose an
opportunity to tackle precisely this problem in a realistic galactic
environment.

The CGM is notoriously difficult to simulate due to its low density, extremely
large volume, and high resolution requirements (Nelson et al., 2016; van de
Voort et al., 2019; Hummels et al., 2019; Faucher-Giguère & Oh, 2023). It is
important to capture both the cosmological context (e.g., merger histories,
intergalactic medium enrichment) and the small-scale processes (e.g.,
turbulence, instabilities) to describe this complex environment (Lochhaas et
al., 2023). For this reason, the IllustrisTNG cosmological box simulations
provide an excellent avenue for exploration of the CGM. The TNG50 simulation
(Nelson et al., 2019a; Pillepich et al., 2019) has high baryonic resolution
and the publicly available subboxes provide high time cadence outputs perfect
for the analysis of HVC origins.

Previous investigations with cosmological simulations have explored the
origins of the CGM as a whole (Hafen et al., 2019), the cold component
(Decataldo et al., 2024), or the high-velocity cold component (Ramesh et al.,
2023). In this work, we offer a complementary view by identifying HVC analogs
in the TNG50 simulation based on criteria specifically motivated by Galactic H
I observations. By tracking these HVCs into the past, we identify their
origins, and investigate their relative contributions to the overall HVC
population.

This paper is organized as follows. In Section 2 we describe the simulations
and our analysis techniques for identifying and tracking the clouds. Section 3
compiles our distribution of origins and explores the specific histories of
several clouds. We discuss our results and conclude in Sections 4 and 5.

##  2 Methods

###  2.1 The TNG Simulations

IllustrisTNG is a suite of high resolution cosmological box simulations using
the magnetohydrodynamic code arepo (Weinberger et al., 2017; Pillepich et al.,
2018). The TNG suite consists of three volumes – TNG50, TNG100, and TNG300
corresponding to simulation cubes of 50, 100, and 300 physical Mpc on a side
(Nelson et al., 2019b). In this work, we use the TNG50-1 simulation which has
mass resolutions for dark matter and gas of $4.5\times 10^{5}$ and $8.5\times
10^{4}$ $\mathrm{M}_{\odot}$, respectively. TNG also produced several
“subboxes” in which a smaller region of the simulation is output at much
higher time cadence ($\sim$6 Myr). subbox0 contains a region of space with
moderate density which contains $\sim$ 6 MW-like galaxies. We select one such
MW-like galaxy corresponding to subhalo ID 537941 at $z=0$ in the full TNG50-1
simulation. This is the only galaxy in the TNG50 MW/M31 sample that is
included in a subbox from early times (Pillepich et al., 2023).

This galaxy has a total mass of $7.64\times 10^{11}$ $\mathrm{M}_{\odot}$ at
$z=0.5$, increasing to $1.02\times 10^{12}$ $\mathrm{M}_{\odot}$ by $z=0$. Its
stellar mass is $\sim 3\times 10^{10}$ $\mathrm{M}_{\odot}$ and its star
formation rate is slightly higher than that of our own Galaxy at $\sim 5$
$\mathrm{M}_{\odot}$/yr.

To visualize the physical gas cell resolution, Figure 1 shows the mean
effective cell radius (the radius of a sphere with equivalent volume) as a
function of height above the plane for our $z=0.25$ snapshot. The bottom axis
shows height in physical kpc while the top axis shows the ratio of height over
virial radius of the halo ($r_{200}=186$ kpc). We include all gas within 200
comoving kpc of the galaxy in our cold cloud analysis (the “cutoff radius”;
237 physical kpc at $z=0.25$). While this resolution is not sufficient to
resolve the smallest clouds, it does capture thermal instabilities and cloud
growth on larger scales (see Section 4).

Figure 1: Mean gas cell spatial resolution as a function of $Z$ height above
the galactic disk. Distances and cell sizes are given in physical kpc. Figure
2: The locations of the identified HVC analogs overlaid on a projection of the
galaxy at our final time ($z\sim 0.25$). The solid circle in the right panel
shows the virial radius of the galaxy ($r_{200}=186$ kpc) and the dashed
circle shows the limit of our cloud identification region (“cutoff radius”,
$r_{\mathrm{cut}}=200\;\mathrm{ckpc}=237\;\mathrm{pkpc}$).

###  2.2 Cloud Identification

To identify contiguous groupings of cold gas cells, we regenerate the Voronoi
mesh for the cells in a given snapshot. The arepo output files contain gas
particles that represent the positions of each Voronoi cell. Thus the Voronoi
tessellation can be reconstructed from the snapshot output of the code. We
then check each cell within 200 comoving kpc to see if it meets our “cold
cloud” criteria. For this work, we are including any cells in which
$T<10^{4.5}$ K. Using a recursive neighbor-finding algorithm, we can obtain a
list of cell IDs which are contiguous and meet our temperature criteria.

This technique results in a list of “clouds” for each snapshot. Each cloud
consists of one or more gas cells with $T<10^{4.5}$ K. We can then calculate
mean properties for each of these clouds by averaging the properties of its
constituent cells (e.g. position, velocity, mass, density, temperature, etc.).
The detected cloud with the largest number of members will be the galactic
disk and we exclude it from our analysis. This method is based on techniques
published in Nelson et al. (2020) and Ramesh et al. (2023). A Cartesian
projection of the galaxy with the identified clouds is shown in Figure 2

###  2.3 Cloud Tracking

Once the clouds in each snapshot have been identified, we track them in time
to generate a merger tree for each cloud (inspired by Jeffreson et al. 2021).
For every cloud we propagate each of its member cells forward in time based on
their current positions and velocities. We then compute the convex hull of the
future cell positions and compare against all of the clouds in the next
snapshot. If more than 10% of the gas cells in any of the detected clouds in
this next snapshot are contained within the convex hull of the propagated
cloud, we link these two clouds. The merger tree itself is generated using a
graph matrix (specifically, the python implementation
NetworkX111https://networkx.org). Each cloud is a node in the graph and clouds
are linked with graph “edges.”

We then prune each connected component of the graph to determine a single
evolutionary track for each cloud. By removing dead end branches and taking
the largest clouds upon splitting, we are left with one cloud per snapshot
that we can follow in the analysis below.

This is distinct from previous works in that we are tracking the observable
cold cloud itself in time as opposed to using tracer particles to follow the
motion of individual gas parcels. Due to gas mixing driven by radiative
cooling (e.g., Gronke & Oh, 2018), the actual gas parcels within the cloud
will change even if the global properties of the cloud are relatively
constant. In this paper, we want to track these physical clouds to identify
how their initial seeds formed. However each cloud is formed at a certain time
and before that time, the above tracking method doesn’t work.

The first snapshot in the evolutionary track of each cloud is its initial
snapshot. Prior to this time, the cloud either has too few members or its
temperature has risen enough that it is no longer identified using the
algorithm described above (Section 2.2). In order to track these clouds back
in time beyond this point, we resort to two imperfect methods. The most
natural choice would be to use the Monte-Carlo tracer particles that reside in
the cloud at this initial snapshot (Genel et al., 2013). However, due to the
overdiffusive nature of the tracers in arepo, we have also tracked these cloud
precursors using gas cell IDs. Obviously using the IDs does not consider the
gas flux between cells during their evolution, however out in the CGM, we
expect these fluxes to be relatively small and this will allow us to explore
the nondiffusive limit when compared with the overdiffusive tracers. As we
show below, we obtain consistent results when using both methods.

Figure 3: The distribution of high velocity gas. On the left, we show two
perspectives in cartesian coordinates with the approximate extent of the disk
outlined in white. On the right, we show the high velocity gas projected onto
the sky assuming a solar-like position within the galaxy. The several large
complexes near the plane are physically connected to the disk and are not
included in our HVC analysis. The top panel shows column density and the
bottom shows GSR velcity with the identified HVCs outlined in white.

###  2.4 HVC Analogs

In order to determine those clouds that are most similar to HVCs observed
around the MW, we use a mock observation technique following Westmeier (2018).
We first determine the expected line of sight velocity due to disk rotation
using the following formula.

$v_{\mathrm{disk}}(\ell,b,R)=\left(v_{\mathrm{rot}}(R)\frac{R_{0}}{R}-v_{\mathrm{rot}}(R_{0})\right)\sin(\ell)\cos(b)$
(1)

Here, $\ell$ and $b$ are Galactic longitude and latitude, respectively, $R$ is
cylindrical radius, $v_{\mathrm{rot}}$ is the circular rotation curve for the
galaxy, and $R_{0}$ is the Galactocentric radial distance to the “sun” (taken
to be $R_{0}=8$ kpc). The angular position of the sun was randomly determined,
however we found that this did not play a significant role in the distribution
of HVC origins (see Appendix A). We calculate the rotation curve using the
formula

$v_{\mathrm{rot}}(R)=\sqrt{\frac{GM_{\mathrm{enc}}(R)}{R^{2}}}$ (2)

where $G$ is Newton’s gravitational constant and $M_{\mathrm{enc}}(R)$ is the
total mass enclosed within radius $R$.

Finally, we transform the positions and velocities of all the cold clouds
identified above (Section 2.2) from Galactocentric coordinates into spherical
Galactic coordinates. We then compare the line of sight radial velocity with
the calculated $v_{\mathrm{disk}}$ at that position. Our HVC-analogs are those
cold clouds for which the difference is larger than the deviation velocity of
70 km s-1 (as used in Westmeier 2018).

This same methodology is used on each individual gas cell in the simulation in
order to reproduce an all-sky map of the high velocity gas. As in the
observations (Westmeier, 2018), we calculate the range of possible disk
velocities for each line of sight by stepping out in $R$ (assuming the disk
extends to cylindrical $R=20$ kpc, and $|Z|=5$ kpc). We then exclude all gas
particles along that line of sight within 70 km s-1 of this range. Figure 3
shows the 3D Cartesian representation of this high velocity gas edge-on and
face-on in the left two panels. The white outline represents the approximate
location of the galactic disk. The top right panel shows how this would appear
on the sky in column density from an observer at the solar location within
this galaxy. The bottom right panel shows this perspective now colored by line
of sight velocity (in the Galactic standard of rest, GSR, frame). The
lightness represents column density. This panel also outlines the projected
positions of the detected HVC analogs shown in Figure 2. We note that there
are several large complexes near the equator that are not outlined. This is
because they are connected to the galactic disk and are thus excluded from our
HVC identification algorithm; however, this gas still appears at anomalous
velocities due to the warp of the disk, thus appearing in this figure.

##  3 Results

Our fiducial snapshot in which we identify our HVC analogs is snapshot 3230
($z=0.245$). We find 70 cold clouds with high deviation velocities comprised
of more than 10 Voronoi cells below $10^{4.5}$ K. We trace these 70 clouds
back in time using the method described in Section 2.3 and track their
properties in order to determine their origin.

###  3.1 Comparison with the MW

Before tracking our HVCs back to study their origins, we compared the
simulated population against current observations of the MW. As stated above,
the global properties of the TNG galaxy presented here agree well with those
of the MW. The properties of the high-velocity gas also provide good
agreement. At its limiting column density ($2.3\times 10^{18}$ cm-2), the
HI4PI survey shows a covering fraction of 11% for the gas with deviation
velocities $>70$ km s-1(Westmeier, 2018) while our galaxy has a covering
fraction of 7.5% using the same criteria (Figure 3).

Furthermore, we compare directly to the discrete HVC population identified by
Moss et al. (2013) in the Galactic All-Sky Survey (GASS; McClure-Griffiths et
al. 2009). Restricting to the parameter space accessible by both the survey
and our simulation (peak $N_{\mathrm{HI}}>10^{18}$ cm-2 and cloud area $>2$
deg2) we find 41 simulated HVCs compared with 30 HVCs in the observational
sample (covering the southern sky). And while the median of the peak column
density distributions agree well at $10^{19.4}$ and $10^{19.6}$ cm-2 for the
simulated and observed HVCs, respectively, the areas of the simulated HVCs are
generally larger with a distribution median at 35 deg2 compared with 6 deg2
for the GASS sample. We believe this is a resolution effect and plan to
explore the distribution of HVC projected areas in future works. These results
align with previous works that have generally found good agreement between
high velocity gas in the TNG simulations and observations (Nelson et al.,
2020; Ramesh et al., 2023).

Figure 4: Properties of the HVCs used in classification of their origins. From
left to right, we sequentially categorize each HVC based on four properties
plotted with respect to initial 3D radial position of the clouds. First, we
use the distance to the closest gas-rich subhalo at $z=0.5$ ($\Delta
r_{\mathrm{sat,0}}$). Any point within the orange shaded region is thus
classified as a satellite origin HVC. Next, we define disk origin HVCs as any
cloud with high maximum circularity ($(L_{z}/|L|)_{\mathrm{max}}$) or positive
initial radial velocity ($v_{\mathrm{rad,0}}$), again with the shaded regions
denoting our selection. Finally, the remaining HVCs are classified as forming
through thermal instability or out of the enriched CGM based on a cut in
maximum temperature reached over the lifetime of the cloud
($T_{\mathrm{max}}$). The bottom right plot shows the maximum value reached
for the ratio of the cooling time to the free-fall time
($(t_{\mathrm{cool}}/t_{\mathrm{ff}})_{\mathrm{max}}$) to confirm that the
thermally unstable clouds have values less than 10. The points are also
colored and shaped based on their determined origin throughout this process –
stars for satellites, ovals for disk, dots for thermal instability, and plus
symbols for enriched CGM.

###  3.2 Identifying Origins

In order to identify the origins of these clouds, we locate each HVC analog in
a six dimensional space of its properties. These individual components are
calculated by taking the mean value across all the cloud cells at each
snapshot. Several are the initial values corresponding to the property at
$z=0.5$ (snapshot 2900), while others are the maximum value attained
throughout their history.

  * •

$r_{0}$ – The initial 3D radius of the cloud cells.

  * •

$(L_{z}/|L|)_{\mathrm{max}}$ – The maximum circularity of the cloud in its
history. The $Z$ direction is aligned with the angular momentum of the
disk.222Throughout this paper, we will be considering specific angular
momentum, namely $\mathbf{L}=\mathbf{r}\times\mathbf{v}$.

  * •

$T_{\mathrm{max}}$ – The maximum temperature reached by the cloud cells in
their history.

  * •

$\Delta r_{\mathrm{sat},0}$ – The minimum separation between the cloud
position and all gas-rich satellites (identified using the SUBFIND algorithm)
at the initial snapshot.

  * •

$v_{\mathrm{rad},0}$ – The initial radial velocity of the cloud.

  * •

$(t_{\mathrm{cool}}/t_{\mathrm{ff}})_{\mathrm{max}}$ – The maximum value of
the ratio of the cooling time to the free fall time over the cloud’s history.

Satellite origin HVCs must have a small separation from an identified subhalo.
Disk origin HVCs exhibit high circularities and positive initial radial
velocities. Thermal instability origins must have reached high maximum
temperatures and maintained cooling times short relative to their freefall
time. While in ideal scenarios this ratio should be less than 1, simulations
have shown that cooling can occur for values less than 10 in realistic
environments (McCourt et al., 2012; Sharma et al., 2012). The remaining clouds
are all relatively metal-enriched at 0.1 $Z_{\odot}$ and found at all initial
radii and radial velocities. We therefore classify these clouds as the
“enriched CGM.” Our selection criteria are listed below and shown in Figure 4
as shaded regions colored by the different origins.

Satellite:

$(r_{0}>r_{\mathrm{vir}})\quad\&\quad(\Delta
r_{\mathrm{sat},0}<50\;\mathrm{kpc})$

Disk:

$(r_{0}<50\;\mathrm{kpc})\quad\&\\\
\Big{(}((L_{z}/|L|)_{\mathrm{max}}>0.7)\quad|\quad(v_{\mathrm{rad},0}>0\;\mathrm{km/s})\Big{)}$

Thermal instability:

$(T_{\mathrm{max}}>T_{\mathrm{vir}}/2)\quad\&\quad\Big{(}(t_{\mathrm{cool}}/t_{\mathrm{ff}})_{\mathrm{max}}<10\Big{)}$

Enriched CGM:

$(T_{\mathrm{max}}<T_{\mathrm{vir}}/2)\quad\&\quad(\Delta
r_{\mathrm{sat},0}>50\;\mathrm{kpc})$

In addition to using $\Delta r_{\mathrm{sat},0}$ to identify satellite origin
HVCs, we isolated the tracer particles associated with each subhalo at $z=0.5$
and propogated them forwards to our final snapshot. We see significant overlap
between satellite particles and the largest HVC (located in the southern
hemisphere) and nearby clouds. The stellar component of this satellite has not
yet merged with the galaxy, lying at a galactocentric radius of 56 kpc. It is
spatially separated from the gas by a distance of 59 kpc due to ram pressure
and drag forces. The stellar component is still colocated with an identified
SUBFIND subhalo.

Interestingly, we also see two other satellites in different stages of
disruption. For the first, both the stars and the gas have been integrated
into the galactic disk. For the second, the gas content of the satellite now
resides within the central galaxy’s disk, however the stellar component (and
DM halo) is still undergoing tidal disruption and tidal tails can clearly be
seen in the distribution of stars.

###  3.3 Distribution of Origins

Tracking gas cell IDs, we find 12 HVCs originating from the galactic disk (12%
by mass), 15 stripped from satellite galaxies (36% by mass), 16 forming
through thermal instability (18% by mass), and 27 that are part of the
enriched CGM (34% by mass). These categories are discussed further below, and
this distribution is shown in the top panels of Figure 5.

Figure 5 also shows the distribution of origins using properties of the tracer
particles instead of the gas cell IDs in the bottom panels. Using tracer
particles we find many more clouds classified as formed via thermal
instability because the maximum temperatures reached by these particles are on
average higher due to their increased diffusivity. Since the hot phase of the
CGM is the volume-filling component, the tracers are more likely to diffuse
into this phase while the gas cells themselves are more likely to remain in
the dense, mass-dominant phase. This change has only a very minor effect on
the disk and satellite origin HVCs. The specific cuts and results for the
tracer analysis are discussed in Appendix B.

Figure 5: The distribution of HVC origins. On the left we count up the number
of individual clouds each coming from the disk, thermal instability, the IGM,
and satellite accretion. While on the right we sum up the total mass of all
the clouds in each category. In the top plots, the origins were determined
based on properties of the tracked particle IDs. While in the bottom plots, we
used the properties of the tracer particle locations.

Figure 6: Properties of example clouds originating from the disk and a
satellite. For each example cloud, we show four panels on the left depicting
the history of the HVC analog and its precursor particles or tracers. The top
half of each panel shows the history by tracking particle IDs, while the
bottom half shows the history of the tracer particles. At the dashed vertical
line, the cloud is identified as a contiguous group of cold Voronoi cells and
is tracked using the method described in Section 2.3. To the left of the
dashed vertical line, the blue distribution shows the history of each cell or
tracer particle while the black line denotes the mean of all the cells. From
the top left, clockwise, the panels depict temperature, the ratio of
$t_{\mathrm{cool}}$ to $t_{\mathrm{ff}}$, the circularity ($L_{z}/|L|$), and
the 3D radial velocity. The four panels on the right show the appearance of
the cloud at the final snapshot ($z\sim 0.25$). The top two panels show a
slice through the domain centered on the HVC analog and show 3D mass density
on the left, and temperature on the right. An approximate demarkation of the
galactic disk is shown as a rectangle (if applicable). The bottom two panels
show projected gas density in the background with an overlay of the HVC
trajectory shown in blue (the dotted line denotes when it is being tracked by
particle IDs, while the solid line is once it is identified as a cold cloud)
and the cloud’s final position as a red plus symbol.

Figure 7: Continuation of Figure 6 showing the properties of example clouds
originating via thermal instability and the enriched CGM.

###  3.4 Cloud Histories

We are also able to track the individual HVCs and their precursors throughout
cosmic time. Figures 6 and 7 show four example clouds, one for each origin
pathway. The panels on the left show the evolution of the cloud’s temperature,
cooling time ($t_{\mathrm{cool}}/t_{\mathrm{ff}}$), angular momentum
($L_{z}/|L|$), and radial velocity. The black lines on the right sides of
these plots (in the region labelled “Cold Cloud”) track the mean properties of
the identified HVC as a collection of Voronoi cells with $T<10^{4.5}$ K.

However, as we move backward in time, at some point the cloud heats up and we
can no longer track it this way. Here we show two different methods of
tracking the HVC’s precursor gas in the two sub-panels in each plot. On the
top we show the evolution of the specific gas cells that comprise the HVC at
first identification (using the cell IDs). The black line is the mean value
and the blue lines are the individual gas cells. As discussed in Section 2.3,
this method does not fully capture the gas flows between cells, however it
provides a limit in the minimum diffusion case.

Figure 8: The on-sky and radial distributions of the HVCs by origin. The left
panel shows all the high-velocity gas (as in Figure 3) with each of the HVCs
outlined in a color representing their origin. The right panel shows the
distribution of HVCs as a function of Galactocentric radius (with each curve
offset in y for clarity).

In the bottom sub-panel we show the evolution of the gas tracer particles that
were in the HVC at identification. Again, the black line is the mean, but the
blue background depicts a histogram showing the distribution of tracer
properties (since there were so many tracers that individual lines became
unreadable). As discussed in Genel et al. (2013), the tracers diffuse more
quickly than the gas and thus this panel gives us a limit in the maximum
diffusion case. By comparing these two methods, we can be sure that our
analysis is consistent and robust.

The four panels on the right depict trajectory of the HVC as well as its final
position and appearance. The top panels show the density and temperature of
the HVC (located in the center of each plot) with the galactic disk shown
outlined in a white rectangle. The bottom panels show two zoomed out views of
the galaxy with the HVC’s trajectory shown in blue and its final position
shown as a red plus. The blue dashed line is the mean trajectory of the
precursor traced by gas cell IDs, while the solid blue line is the mean
trajectory of the tracked cold cloud.

For the disk origin example (Figure 6, top), we can see that it remains cold
throughout its history and has a circularity of $\approx$1 for the first
$\sim$half of its evolution. For the satellite origin HVC (Figure 6, bottom),
it also remains relatively cool however it has a consistent angular momentum
and radial velocity as it approaches the galaxy from large distances. For the
example cloud of thermal instability origin (Figure 7, top), we see a steady
decrease in temperature until the cloud forms with
$t_{\mathrm{cool}}/t_{\mathrm{ff}}$ remaining below 1 throughout. Finally, the
enriched CGM origin (Figure 7, bottom) shows a relatively constant, cool
temperature with slowly varying radial velocity and smooth trajectory. This
category has several different populations of clouds with a variety of initial
radii and radial velocities. However, they all have cool max temperatures and
metallicities $\sim 0.1$ Z⊙.

###  3.5 Observable Properties

Figures 8 and 9 show the distributions and properties of the HVC sample by
origin. Figure 8 shows the onsky (left) and radial (right) distributions for
each of the different formation pathways. Most of the satellite origin clouds
are in the southern hemisphere close to the location of the infalling
satellite. And the inner-most clouds are generally from satellites or the
disk.

Figure 9 shows the distributions of metallicities (left) and galactocentric
radial velocities (right) for our HVC sample. Each panel shows the
distribution for a different origin (classified as shown in Figure 4). In the
left panels, the vertical dashed lines denote the mean metallicities of the
source gas (disk, satellite ISM, and circumgalactic gas). Additionally, the
grey background shows the initial distribution of the cells comprising the
HVCs at $z=0.5$. We see that most of the HVCs have metallicities of $z\sim
0.1$ Z⊙, and are infalling towards the galaxy. We discuss the implications of
these results for observations below.

##  4 Discussion

While we have focused on the origin of HVC analogs, there have been several
other works exploring the origin of the CGM as a whole. Hafen et al. (2019)
use the FIRE simulations (Hopkins et al., 2014, 2018) to study the CGM of MW-
mass halos. For $z\sim 0.25$, they find 16% of the CGM comes from satellites,
20% from recycling of the galaxy’s own winds, and 53% originates in the IGM.
Decataldo et al. (2024) also identify the origins of the cold CGM using the
Eris2K simulations. For $z\sim 0.5$, they find 1.5% entering the CGM from the
IGM in the cold phase, 2.5% from satellites, 10% from the galactic disk, and
84% cooled from hot gas while bound to the galaxy.

These classifications vary slightly across the literature, but there are a few
key differences with our results here. Notably the HVC analog mass fraction
originating from satellites is 36% compared with 2.5% of the cold CGM in
general (Decataldo et al., 2024) and 16% of the CGM in total (Hafen et al.,
2019). This shows that stripped material from satellites is more likely to
appear at high velocities relative to the galactic disk rotation than
thermally unstable CGM gas or galactic outflows.

Figure 9: The final snapshot properties of the HVCs separated by origin. The
left four panels show the final snapshot metallicity distribution in color and
the initial distribution of the cloud progenitor cells in grey. The vertical
lines denote the mean metallicity of the disk, CGM, IGM, and satellite
galaxies for reference. On the right we show the distribution of 3D radial
velocities with positive velocities denoting outflow. As expected, disk origin
HVCs are the only outflowing clouds.

More consistent with the global CGM is the distribution of the remaining
material. Roughly 15% of the HVC analogs (by mass) originate in the disk,
while the remaining material cools in the CGM ($\sim$50%). This is comparable
to the previous results listed above.

However, it has also been shown that the gaseous histories of galaxies depend
greatly on numerical effects (e.g. Nelson et al., 2013; Hobbs et al., 2013).
Early smoothed-particle hydrodynamics (SPH) codes were unable to accurately
model various instabilities which led to the increased stability of cold
gaseous filaments in the CGM. Using more modern simulation techniques, we have
seen that the fragmentation of circumgalactic gas is ubiquitous with
resolution and numerical methods playing an important role in the realism of
the results (van de Voort et al., 2019; Peeples et al., 2019b; Suresh et al.,
2019; Ramesh & Nelson, 2024).

On smaller scales, wind tunnel simulations have shown that cold clouds moving
through a hot gaseous medium can survive and grow (depending on the initial
properties; Kwak et al. 2011; Gronke & Oh 2018; Gronke et al. 2022). The
criteria for cloud survival is often given in terms of critical size; large
clouds will grow, and small clouds will shatter. Pulling values from our
simulation, we see that our clouds have overdensities of $4-10$ and Mach
numbers of $0.2-0.9$. Using equation (2) from Gronke et al. (2022), we find
critical radii on the order of a few parsecs. Due to the resolution of the
simulation, all our clouds have sizes $>1$ kpc, so we would expect them to
grow. They also specify a resolution requirement to resolve the cloud by at
least $\sim$ 16 cells (Tan et al., 2021), and our clouds contain $13-137$
cells at identification which should be sufficient in most cases.

Moreover, we can compare the simulation resolution to an estimate of the
fragmentation length to confirm that we are able to resolve the thermal
instability. Maller & Bullock (2004) derive the size of clouds formed in a
two-phase pressure confinement model of thermal instability driven
fragmentation. Using their equation (26), we can calculate the expected sizes
of clouds that will form in the CGM of our galaxy. By using the mass of cold
connected cells and the temperature of the surrounding ambient gas, we find
cloud sizes of $1-5$ kpc which is in agreement with their actual sizes
($r_{\mathrm{HVC}}/\ell_{\mathrm{frag}}\sim 0.5-1.75$). This reassures us that
at the resolution of TNG50-1, we will be able to see clouds condense due to
thermal instability in the CGM.

Finally, we return to the observability of these different origin pathways. As
shown in Figure 9, we see that most of the present-day HVCs have metallicities
of $\sim 0.1$ Z⊙, presumably due to mixing. However, we do see a few lower
metallicity clouds originating from satellite accretion. While we do not know
the origins of many of the MW HVCs, this does seem to be in agreement based on
the metallicity observations that we do have. The Smith Cloud and Complex C
have metallicities of $0.3-0.5$ Z⊙ (Fox et al., 2016, 2023) while the
Magellanic Stream and Leading Arm (material tidally stripped from the Large
and Small Magellanic Clouds) have lower metallicities of $\sim 0.1$ Z⊙ (Fox et
al., 2018).

In the right panels of Figure 9, we see that nearly all the simulated HVCs are
infalling (negative $v_{\mathrm{rad}}$). The only formation pathway that leads
to significant outflowing gas is disk origin, as expected. In our
observational sample of HVCs, we also see the majority of the gas infalling
towards us, however transforming from Heliocentric coordinates to
Galactocentric coordinates without proper motions for these gas clouds can be
imprecise (Wakker & van Woerden, 1997). Most of the positive velocity clouds
are very small, which will allow them to be more easily pushed and pulled by
turbulence in the inner CGM.

###  4.1 Further Considerations

Gas mixing plays a key role in these processes and as stated above, has a
strong dependence on numerical methodologies (Springel, 2010; Vogelsberger et
al., 2012; Sijacki et al., 2012; Hayward et al., 2014). It has been shown that
SPH codes are unable to produce certain instabilities which results in a lack
of mixing and cooling while adaptive-mesh refinement (AMR) and Eulerian codes
suffer from overcooling and preferred directions. arepo’s moving-mesh scheme
based on Voronoi tessellation mitigates these issues and allows for
appropriate gas mixing and cooling. While a full comparison between
cosmological scales and those of the turbulent mixing layers governing gas
mixing should be performed (Fielding et al., 2020), arepo’s numerical model is
theoretically able to capture instabilities and hydrodynamic mixing correctly.
However, resolution of the simulations plays an important role in this issue
as well.

When considering simulations, especially studying the CGM, resolution and
computational feasability are continuously at odds. IllustrisTNG was designed
to be a set of cosmological box simulations, so while they are cutting edge,
the baryonic resolution is limited to $4.5\times 10^{5}$ $\mathrm{M}_{\odot}$.
However, as stated above, the spatial resolution is sufficient to resolve the
fragmentation length for thermal instabilities (Maller & Bullock, 2004). This
is consistent with previous works that have studied the thermally unstable CGM
in TNG50 (Nelson et al., 2020). A full convergence study is required and we
will cover this in future work utilizing cosmological zoom-in simulations with
higher resolution.

##  5 Conclusions

In this paper, we have identified 70 HVC analogs in a TNG50-1 galaxy. Using
the high time cadence of subbox0, we are able to track these HVCs in time to
determine their origins.

  * •

We find a realistic distribution of high velocity gas which is well captured
by our algorithm.

  * •

By identifying cold clouds at each snapshot and propagating them forward in
time, we can track the transient HVC as gas cycles in and out of it.

  * •

For the first time, we are able to identify HVC analogs formed via thermal
instability. Due to the intense turbulence in the CGM, these clouds can be
pushed to high deviation velocities consistent with the observed HVCs in the
MW. 23% of our HVC analogs are of thermal instability origin.

  * •

We also find 17% originating from the disk, 21% coming from satellites, and
39% comprised of relatively high metallicity circumgalactic gas at a variety
of radii (the enriched CGM).

  * •

We obtain similar results by using both gas cell IDs and Monte Carlo tracer
particles to track the HVC analogs back before they were cold.

  * •

Based on this analysis, low metallicity HVCs are most likely to be associated
with a satellite infall, and positive radial velocity clouds are most likely
of disk origin.

As we increase the resolution in the CGM, more and more structure continues to
appear. By applying this technique to new, cosmological zoom and CGM-
refinement simulations, we will be able to get a better understanding of the
low-mass end of the HVC population. Moreover, we will be able to explore the
future of these clouds as they interact with the galactic disk and facilitate
star formation and galaxy growth.

The authors thank Sarah Jeffreson, Jake Bennett, and Cassi Lochhaas for useful
discussions in the development of this paper. The computations in this paper
were run on the FASRC cluster supported by the FAS Division of Science
Research Computing Group at Harvard University.

Figure 10: Similar to Figure 5 showing the number of identified HVCs by origin
for different solar positions. The left-most panel matches that in Figure 5
and the subsequent three panels show the distributions after rotating the
solar position by 90, 180, and 270 degrees.

##  Appendix A Azimuthal Solar Position

Our galaxy was initially rotated from its position within the TNG box such
that the disk’s angular momentum was aligned with the $z$-axis. Specifically,
it was rotated by the angle $\cos^{-1}(\hat{L}\cdot\hat{z})$ around the vector
$\hat{L}\times\hat{z}$. The fiducial azimuthal position of the “Sun” is
therefore $(x,y,z)=(-8,0,0)$ kpc in this rotated coordinate system.

We also performed the full analysis after rotating the galaxy by 90, 180, and
270 degrees to determine the impact, if any, of the azimuthal solar position
on our results. Using the full set of 211 identified cold clouds at $z=0.25$,
we calculated the deviation velocities for each of these different viewing
angles. We found the number of clouds identified as “high velocity” ranged
from 68 to 78 with 41 consistent across all four angles. Furthermore, Figure
10 shows the distribution of origins for each of these solar positions. While
there is some slight variation, the overall relative distributions are
consistent.

##  Appendix B Tracer Particle Results

Figure 5 shows the distribution of origins determined for both particle ID
tracking and tracer particle tracking. However, due to the inflated
temperatures seen by the tracer particles, a few modifications to the
selection criteria were required. The cutoff temperature between thermal
instability and the enriched CGM origins was increased from
$T_{\mathrm{vir}}/2$ to $T_{\mathrm{vir}}$, and the minimum separation to a
satellite, $\Delta r_{\mathrm{sat},0}$, was also increased from 50 to 60 kpc.
The resulting property distributions are shown in Figure 11.

Figure 11: As Figure 9 except with the origins determined using the tracer
particle values. This leads to an increase in the number of disk origin clouds
and thermal instability clouds, while decreasing the number of IGM clouds.
However, the overall trends are unchanged. Figure 12: The metallicity
distribution of gas within and beyond 300 kpc as a function of radial
velocity. On the left we show $r<300$ kpc gas, and on the right $r>300$ kpc
gas. The shaded curves on the far left and right are the total metallicity
histograms for the $r<300$ kpc and $r>300$ kpc gas, respectively. A vertical
line denotes zero velocity with inflows as negative and outflows as positive.

##  Appendix C Metallicity of the IGM

Interestingly, we did not find any HVCs originating from very low metallicity
IGM gas. All of the clouds that form at large radii have relatively high
metallicities of $\sim 0.1$ Z⊙. This seems to be due to very strong outflows
coming from the galaxy’s supermassive black hole which is propelling high
metallicity material out to hundreds of kpc away from the galaxy (Suresh et
al., 2015). These outflows prevent IGM gas from falling directly into the
galaxy’s CGM, and mixes strongly with the low and intermediate metallicity
material. Additionally, small-scale simulations show that high velocity gas
experiences very strong mixing and can raise or lower a cloud’s metallicity
dramatically within a few hundred Myrs (Gritton et al., 2014).

Figure 12 shows the distribution of gas within (left) and beyond (right) 300
kpc in metallicity vs radial velocity. Inflows correspond to negative
$v_{\mathrm{rad}}$ while outflows are positive. On the left, we can see
inflowing low metallicity material and outflowing high metallicity material.
On the right, we see a similar trend even at these very large radii. There are
metal enriched outflows at the level of $>0.1$ Z⊙ in addition to low
metallicity infalling material. The integrated histograms for these two radial
cuts are shown on the left and right (for the $r<300$ kpc gas and the $r>300$
kpc gas, respectively). Interestingly, the very low metallicity gas ($Z<0.01$
Z⊙) does not survive within 300 kpc. This is consistent with our findings that
the HVC analogs that we detect originating at large radii have $Z\sim 0.1$ Z⊙
indicating that they were formed out of this metal rich outflow rather than
the pristine intergalactic medium.

## References

  * Binney (1977) Binney, J. 1977, ApJ, 215, 483, doi: 10.1086/155378
  * Bregman (1980) Bregman, J. N. 1980, ApJ, 236, 577, doi: 10.1086/157776
  * Chevalier & Oegerle (1979) Chevalier, R. A., & Oegerle, W. R. 1979, ApJ, 227, 398, doi: 10.1086/156744
  * Decataldo et al. (2024) Decataldo, D., Shen, S., Mayer, L., Baumschlager, B., & Madau, P. 2024, A&A, 685, A8, doi: 10.1051/0004-6361/202346972
  * Faucher-Giguère & Oh (2023) Faucher-Giguère, C.-A., & Oh, S. P. 2023, ARA&A, 61, 131, doi: 10.1146/annurev-astro-052920-125203
  * Field (1965) Field, G. B. 1965, ApJ, 142, 531, doi: 10.1086/148317
  * Fielding et al. (2020) Fielding, D. B., Ostriker, E. C., Bryan, G. L., & Jermyn, A. S. 2020, ApJ, 894, L24, doi: 10.3847/2041-8213/ab8d2c
  * Fox et al. (2023) Fox, A. J., Cashman, F. H., Kriss, G. A., et al. 2023, ApJ, 946, L48, doi: 10.3847/2041-8213/acc640
  * Fox et al. (2016) Fox, A. J., Lehner, N., Lockman, F. J., et al. 2016, ApJ, 816, L11, doi: 10.3847/2041-8205/816/1/L11
  * Fox et al. (2018) Fox, A. J., Barger, K. A., Wakker, B. P., et al. 2018, ApJ, 854, 142, doi: 10.3847/1538-4357/aaa9bb
  * Genel et al. (2013) Genel, S., Vogelsberger, M., Nelson, D., et al. 2013, MNRAS, 435, 1426, doi: 10.1093/mnras/stt1383
  * Gritton et al. (2014) Gritton, J. A., Shelton, R. L., & Kwak, K. 2014, ApJ, 795, 99, doi: 10.1088/0004-637X/795/1/99
  * Gronke & Oh (2018) Gronke, M., & Oh, S. P. 2018, MNRAS, 480, L111, doi: 10.1093/mnrasl/sly131
  * Gronke et al. (2022) Gronke, M., Oh, S. P., Ji, S., & Norman, C. 2022, MNRAS, 511, 859, doi: 10.1093/mnras/stab3351
  * Habe & Ikeuchi (1980) Habe, A., & Ikeuchi, S. 1980, Progress of Theoretical Physics, 64, 1995, doi: 10.1143/PTP.64.1995
  * Hafen et al. (2019) Hafen, Z., Faucher-Giguère, C.-A., Anglés-Alcázar, D., et al. 2019, MNRAS, 488, 1248, doi: 10.1093/mnras/stz1773
  * Hayward et al. (2014) Hayward, C. C., Torrey, P., Springel, V., Hernquist, L., & Vogelsberger, M. 2014, MNRAS, 442, 1992, doi: 10.1093/mnras/stu957
  * Hobbs et al. (2013) Hobbs, A., Read, J., Power, C., & Cole, D. 2013, MNRAS, 434, 1849, doi: 10.1093/mnras/stt977
  * Hopkins et al. (2014) Hopkins, P. F., Kereš, D., Oñorbe, J., et al. 2014, MNRAS, 445, 581, doi: 10.1093/mnras/stu1738
  * Hopkins et al. (2018) Hopkins, P. F., Wetzel, A., Kereš, D., et al. 2018, MNRAS, 480, 800, doi: 10.1093/mnras/sty1690
  * Hummels et al. (2019) Hummels, C. B., Smith, B. D., Hopkins, P. F., et al. 2019, ApJ, 882, 156, doi: 10.3847/1538-4357/ab378f
  * Jeffreson et al. (2021) Jeffreson, S. M. R., Keller, B. W., Winter, A. J., et al. 2021, MNRAS, 505, 1678, doi: 10.1093/mnras/stab1293
  * Kalberla & Haud (2006) Kalberla, P. M. W., & Haud, U. 2006, A&A, 455, 481, doi: 10.1051/0004-6361:20054750
  * Kereš & Hernquist (2009) Kereš, D., & Hernquist, L. 2009, ApJ, 700, L1, doi: 10.1088/0004-637X/700/1/L1
  * Kwak et al. (2011) Kwak, K., Henley, D. B., & Shelton, R. L. 2011, ApJ, 739, 30, doi: 10.1088/0004-637X/739/1/30
  * Lehner et al. (2022) Lehner, N., Howk, J. C., Marasco, A., & Fraternali, F. 2022, MNRAS, 513, 3228, doi: 10.1093/mnras/stac987
  * Lochhaas et al. (2023) Lochhaas, C., Tumlinson, J., Peeples, M. S., et al. 2023, ApJ, 948, 43, doi: 10.3847/1538-4357/acbb06
  * Maller & Bullock (2004) Maller, A. H., & Bullock, J. S. 2004, MNRAS, 355, 694, doi: 10.1111/j.1365-2966.2004.08349.x
  * McClure-Griffiths et al. (2009) McClure-Griffiths, N. M., Pisano, D. J., Calabretta, M. R., et al. 2009, ApJS, 181, 398, doi: 10.1088/0067-0049/181/2/398
  * McCourt et al. (2012) McCourt, M., Sharma, P., Quataert, E., & Parrish, I. J. 2012, MNRAS, 419, 3319, doi: 10.1111/j.1365-2966.2011.19972.x
  * Moss et al. (2013) Moss, V. A., McClure-Griffiths, N. M., Murphy, T., et al. 2013, ApJS, 209, 12, doi: 10.1088/0067-0049/209/1/12
  * Muller et al. (1963) Muller, C. A., Oort, J. H., & Raimond, E. 1963, Academie des Sciences Paris Comptes Rendus, 257, 1661
  * Nelson et al. (2016) Nelson, D., Genel, S., Pillepich, A., et al. 2016, MNRAS, 460, 2881, doi: 10.1093/mnras/stw1191
  * Nelson et al. (2013) Nelson, D., Vogelsberger, M., Genel, S., et al. 2013, MNRAS, 429, 3353, doi: 10.1093/mnras/sts595
  * Nelson et al. (2019a) Nelson, D., Pillepich, A., Springel, V., et al. 2019a, MNRAS, 490, 3234, doi: 10.1093/mnras/stz2306
  * Nelson et al. (2019b) Nelson, D., Springel, V., Pillepich, A., et al. 2019b, Computational Astrophysics and Cosmology, 6, 2, doi: 10.1186/s40668-019-0028-x
  * Nelson et al. (2020) Nelson, D., Sharma, P., Pillepich, A., et al. 2020, MNRAS, 498, 2391, doi: 10.1093/mnras/staa2419
  * Olano (2008) Olano, C. A. 2008, A&A, 485, 457, doi: 10.1051/0004-6361:20077556
  * Peeples et al. (2019a) Peeples, M., Behroozi, P., Bordoloi, R., et al. 2019a, BAAS, 51, 368, doi: 10.48550/arXiv.1903.05644
  * Peeples et al. (2019b) Peeples, M. S., Corlies, L., Tumlinson, J., et al. 2019b, ApJ, 873, 129, doi: 10.3847/1538-4357/ab0654
  * Pillepich et al. (2018) Pillepich, A., Springel, V., Nelson, D., et al. 2018, MNRAS, 473, 4077, doi: 10.1093/mnras/stx2656
  * Pillepich et al. (2019) Pillepich, A., Nelson, D., Springel, V., et al. 2019, MNRAS, 490, 3196, doi: 10.1093/mnras/stz2338
  * Pillepich et al. (2023) Pillepich, A., Sotillo-Ramos, D., Ramesh, R., et al. 2023, arXiv e-prints, arXiv:2303.16217, doi: 10.48550/arXiv.2303.16217
  * Putman et al. (2012) Putman, M. E., Peek, J. E. G., & Joung, M. R. 2012, ARA&A, 50, 491, doi: 10.1146/annurev-astro-081811-125612
  * Ramesh & Nelson (2024) Ramesh, R., & Nelson, D. 2024, MNRAS, 528, 3320, doi: 10.1093/mnras/stae237
  * Ramesh et al. (2023) Ramesh, R., Nelson, D., & Pillepich, A. 2023, MNRAS, 522, 1535, doi: 10.1093/mnras/stad951
  * Richter et al. (2018) Richter, P., Fox, A. J., Wakker, B. P., et al. 2018, ApJ, 865, 145, doi: 10.3847/1538-4357/aadd0f
  * Sharma et al. (2012) Sharma, P., McCourt, M., Quataert, E., & Parrish, I. J. 2012, MNRAS, 420, 3174, doi: 10.1111/j.1365-2966.2011.20246.x
  * Sijacki et al. (2012) Sijacki, D., Vogelsberger, M., Kereš, D., Springel, V., & Hernquist, L. 2012, MNRAS, 424, 2999, doi: 10.1111/j.1365-2966.2012.21466.x
  * Smoker et al. (2011) Smoker, J. V., Fox, A. J., & Keenan, F. P. 2011, MNRAS, 415, 1105, doi: 10.1111/j.1365-2966.2011.18647.x
  * Spitzer (1956) Spitzer, Lyman, J. 1956, ApJ, 124, 20, doi: 10.1086/146200
  * Springel (2010) Springel, V. 2010, MNRAS, 401, 791, doi: 10.1111/j.1365-2966.2009.15715.x
  * Suresh et al. (2015) Suresh, J., Bird, S., Vogelsberger, M., et al. 2015, MNRAS, 448, 895, doi: 10.1093/mnras/stu2762
  * Suresh et al. (2019) Suresh, J., Nelson, D., Genel, S., Rubin, K. H. R., & Hernquist, L. 2019, MNRAS, 483, 4040, doi: 10.1093/mnras/sty3402
  * Tan et al. (2021) Tan, B., Oh, S. P., & Gronke, M. 2021, MNRAS, 502, 3179, doi: 10.1093/mnras/stab053
  * Tripp (2022) Tripp, T. M. 2022, MNRAS, 511, 1714, doi: 10.1093/mnras/stac044
  * Tumlinson et al. (2017) Tumlinson, J., Peeples, M. S., & Werk, J. K. 2017, ARA&A, 55, 389, doi: 10.1146/annurev-astro-091916-055240
  * van de Voort et al. (2019) van de Voort, F., Springel, V., Mandelker, N., van den Bosch, F. C., & Pakmor, R. 2019, MNRAS, 482, L85, doi: 10.1093/mnrasl/sly190
  * van Woerden et al. (1957) van Woerden, H., Rougoor, G. W., & Oort, J. H. 1957, Academie des Sciences Paris Comptes Rendus, 244, 1691
  * Vogelsberger et al. (2012) Vogelsberger, M., Sijacki, D., Kereš, D., Springel, V., & Hernquist, L. 2012, MNRAS, 425, 3024, doi: 10.1111/j.1365-2966.2012.21590.x
  * Voit (2021) Voit, G. M. 2021, ApJ, 908, L16, doi: 10.3847/2041-8213/abe11f
  * Wakker & van Woerden (1997) Wakker, B. P., & van Woerden, H. 1997, ARA&A, 35, 217, doi: 10.1146/annurev.astro.35.1.217
  * Weinberger et al. (2017) Weinberger, R., Springel, V., Hernquist, L., et al. 2017, MNRAS, 465, 3291, doi: 10.1093/mnras/stw2944
  * Westmeier (2018) Westmeier, T. 2018, MNRAS, 474, 289, doi: 10.1093/mnras/stx2757