# Turbulent Gas-rich Disks at High Redshift: Origin of Thick Stellar Disks Through 3D “Baryon Sloshing”

Joss Bland-Hawthorn, Thor Tepper-Garcia, Oscar Agertz, Christoph

Federrath, Misha Haywood, Paola Di Matteo, Timothy R Bedding, Takafumi

Tsukui, Emily Wisnioski, Melissa Ness, et al.

<figure><img src="image_1.png" /></figure>

## ► To cite this version:

Joss Bland-Hawthorn, Thor Tepper-Garcia, Oscar Agertz, Christoph Federrath, Misha Haywood, et al. Turbulent Gas-rich Disks at High Redshift: Origin of Thick Stellar Disks Through 3D “Baryon Sloshing”. The Astrophysical Journal, 2025, 994 (1), pp.22. 10.3847/1538-4357/ae0931 . hal-05495239

HAL Id: hal-05495239

https://hal.science/hal-05495239v1

Submitted on 5 Feb 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
<figure><img src="image_3.png" /></figure>

# Turbulent Gas-rich Disks at High Redshift: Origin of Thick Stellar Disks Through 3D “Baryon Sloshing”

Joss Bland-Hawthorn¹,², Thor Tepper-Garcia¹,², Oscar Agertz³, Christoph Federrath²,⁴, Misha Haywood⁵, Paola di Matteo⁵, Timothy R. Bedding¹, Takafumi Tsukui⁴, Emily Wisnioski⁴, Melissa Ness⁴, and Ken Freeman⁴

¹Sydney Institute for Astronomy, School of Physics, A28, The University of Sydney, NSW 2006, Australia; jonathan.bland-hawthorn@sydney.edu.au

²Centre of Excellence for All-Sky Astrophysics in Three Dimensions (ASTRO 3D), Australia

³Lund Observatory, Division of Astrophysics, Department of Physics, Lund University, Box 43, SE-221 00 Lund, Sweden

⁴Research School of Astronomy and Astrophysics, Australian National University, Canberra, ACT 2611, Australia

⁵GEPI, Observatoire de Paris, CNRS, Université Paris Diderot, 5 Place Jules Janssen, 92190 Meudon, France

Received 2025 January 31; revised 2025 September 9; accepted 2025 September 16; published 2025 November 12

## Abstract

In response to recent observations from JWST and Atacama Large Millimeter Array, we explore a new class of dynamically self-consistent models that mimics a plausible progenitor of the Milky Way over a wide range of disk gas fractions, $f_{gas}$. The high gas surface densities encourage vigorous star formation, which in turn couples with the gas to drive turbulence. We show that this coupling through momentum recoil drives a random walk of the baryonic potential minimum with respect to total gravitational potential, $\Phi_{tot}(R, \phi, z)$. The amplitude of the bulk motion depends on the feedback strength, which in turn is directly associated with $f_{gas}$. At its most extreme, when gas is the sole contributor to the disk potential ($f_{gas} = 100\%$), the amplitude of the walk can reach up to $R \approx 5$ kpc within $\Phi_{tot}$. The disk dominates over dark matter ($f_{disk} \gtrsim 50\%$) within $R_s = 2.2R_{disk}$, where $R_{disk}$ is the exponential disk scale length. For a lower $f_{disk}$ and/or $f_{gas}$, the 3D sloshing amplitude and velocity are reduced. The combination of strong feedback and sloshing leads to the newly formed stars being dynamically heated and settling to a more spatially extended disk population. The 3D heating process is roughly isotropic but its effects are more noticeable in $|z|$ due to the initial dynamical coldness of the star-forming disk. Such a disk has enhanced [α/Fe] stellar abundances and a vertical (but no radial) gradient in stellar age and metallicity, both consistent with the Milky Way’s thick stellar disk.

Unified Astronomy Thesaurus concepts: Galaxy disks (589); Milky Way dynamics (1051); Galaxy chemical evolution (580); Galaxy dynamics (591); Galaxy evolution (594)

## 1. Introduction

### 1.1. The Need for New Simulations

Most disk galaxies observed in the near field today went through an early phase ($z > 1$) when the disk mass was dominated by gas (L. J. Tacconi et al. 2020). In order for a gas disk to support itself against collapse, high star formation rates must have been commonplace, consistent with the high gas surface densities (R. Genzel et al. 2010). The energetic feedback presumably coupled to the gas with moderate efficiency (~10%) and provided a natural source of (turbulent) pressure to sustain the gas-rich medium (O. Agertz et al. 2021; J. Bland-Hawthorn et al. 2024). This was an intensely energetic era that can be traced to $z \sim 10$ by our most powerful telescopes, including the Atacama Large Millimeter Array (ALMA) and the James Webb Space Telescope (JWST). Today, Milky Way analogs have a low gas fraction ($f_{gas} \lesssim 10\%$), and so provide limited information on the gas processes that were active at those early times. Here, we seek to gain some physical insight about the nature of turbulent gas-rich disks at high redshift through controlled simulations made possible with the Nexus framework (T. Tepper-García et al. 2024).

Our new work is in response to mounting evidence for gas-dominant, star-forming disks in the redshift range $z \sim 1-8$ (S. C. Chapman et al. 2004; N. M. Förster Schreiber et al. 2006; R. Genzel et al. 2006; K. L. Shapiro et al. 2008; A. M. Swinbank et al. 2012; J. A. Hodge et al. 2019; T. Tsukui & S. Iguchi 2021; F. Rizzo et al. 2022; F. Roman-Oliveira et al. 2024; L. E. Rowland et al. 2024; T. Tsukui et al. 2024). There are clear cases where the gas entirely dominates the baryon disk fraction (e.g., F. Roman-Oliveira et al. 2024; L. E. Rowland et al. 2024; A. Amyrosiadis et al. 2025; J. A. Hodge et al. 2025). To date, there has been limited work on realistic high-resolution simulations of isolated gas-rich galaxies sustained by underlying star formation and turbulence (e.g., F. van Donkelaar et al. 2022). Here we provide a sophisticated simulation framework (Nexus; see T. Tepper-García et al. 2024) to explore gas-dominant disk systems within the context of multicomponent galaxy models (stellar disk and bulge, live dark matter halo, hot gaseous corona, multiphase ISM with star formation and metal production, etc.). Thus, it is now possible to explore disk evolution over the full range of $f_{gas} = 0\%–100\%$ in controlled experiments, as illustrated in Figure 1. (The parameters $f_{disk}$ and $f_{gas}$ are defined in Section 2.)

Our work has highlighted new processes in early gas-dominant disks that were missed in earlier (lower resolution) simulations. The Nexus simulations reveal large-scale instabilities that arise in turbulent gas disks, triggered by pressure waves from stochastic star formation (see Figure 2). These instabilities have important implications: as we now

<figure><img src="image_4.png" /></figure>

Original content from this work may be used under the terms of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
<figure><img src="image_5.png" /><figcaption>**Figure 1.** Edge-on gas surface density maps taken from our early disk simulations at $t_0 = 1$ Gyr. We adopt $f_{\text{disk}} = 0.5$ for the Milky Way progenitor at $z \sim 3$. The top panels are for five different initial gas fractions (as indicated) with no halo accretion; the bottom five are the matching simulations that include halo accretion from cooling hot gas. From left to right, for both rows, the galaxy's center falls at the right-hand side (RHS) edge in the first four panels, but is centered in the last panel, where the spatial scale is indicated. The extent of the disk-halo interaction is a reflection of the disk-wide star formation shown in Figure 2.</figcaption></figure>

show, they drive disk asymmetries and baryon sloshing modes within the smooth gravitational potential. Here, we focus on the latter and leave disk asymmetries to a companion paper (J. Bland-Hawthorn et al. 2025, in preparation). We show how baryon sloshing provides us with an alternative heating mechanism leading to vertically thickened stellar disks, which are observed in most disk galaxies (P. Yoachim & J. J. Dalcanton 2006; T. Tsukui et al. 2025).

## 1.2. A New Era in Disk Research

The formation of galactic disks has been one of the fastest developing areas of galactic astronomy in recent years. A host of observational results is revolutionizing this field of research, including the discovery of dynamically cold galactic disks at very high redshift ($1 < z < 8$) with ALMA and JWST (e.g., L. Ferreira et al. 2022, 2023; J. S. Kartaltepe et al. 2023; A. Pope et al. 2023, E. Nelson et al. 2024), the existence of already-mature disks containing bars (L. Costantin et al. 2023; B. Ratcliffe et al. 2025) and spiral structures (V. Kuhn et al. 2024) as early as 2 Gyr after the Big Bang. These discoveries concern objects that are expected to be massive, settling to disks earlier than Milky Way-like galaxies (disk downsizing, see T. Tsukui et al. 2025), in some cases evolving into ellipticals in the present-day Universe. Exactly when Milky Way-like galaxies began to build their disks is uncertain (H. J. Mo et al. 1998; I. B. Santistevan et al. 2020). A. Pope et al. (2023) present a galaxy at redshift 4.3 dominated by a stellar disk ($V_{\text{circ}}/\sigma_{\text{LOS}} \approx 5$) with baryon mass $M_* \approx 2 \times 10^9 M_\odot$, a value in line with the estimated stellar mass of Milky

Way (MW) precursors at such epochs (P. G. van Dokkum et al. 2013).

In our own Galaxy, the oldest and most metal-poor stars with disk kinematics are part of the thick disk and its metal-poor extension. The thick disk is traditionally defined in terms of its warm kinematics and its enhanced [α/Fe] abundances compared to the solar value; this population is largely confined to $R < 10$ kpc (M. R. Hayden et al. 2015). Thick disk stars are thus expected to hold the key to understanding and constraining the physical mechanisms that govern how and when these dissipative structures formed. Unfortunately, the epoch and duration of formation of this population cannot be determined precisely from stellar ages, even with recent advances in asteroseismology (S. Sharma et al. 2019). Some studies suggest thick disk formation lasted for 2–3 Gyr (M. Haywood et al. 2013; O. Snaith et al. 2015), while others suggest a significantly faster formation, on the order of 1 Gyr (A. Miglio et al. 2021). In all cases, the bulk of the stellar mass of this population was in place 10 Gyr ago (at $z \approx 2$), with $M_*$ in the range of 20% (J. Bland-Hawthorn & O. Gerhard 2016) to 50% (O. N. Snaith et al. 2014) of the total present-day stellar mass of the disk.

The existence of such stellar disks at high redshifts requires that massive and early inflows of gas occurred in the inner kiloparsec regions of galaxies. Our understanding of the mechanism by which gas fed the inner regions of early galaxies has changed significantly in the last 20 yr, with the identification of a cold mode of gas accretion in which gas is brought along a few dark matter filaments (Y. Birnboim & A. Dekel 2003; D. Kereš et al. 2005, 2009; A. Dekel & Y. Birnboim 2006). Although cold-mode accretion provides a
<figure><img src="image_6.png" /><figcaption>**Figure 2.** Periodograms drawn from our heavy disk ($f_{\text{disk}} = 0.5$) simulations that capture how the star formation rate (SFR) surface density (in units of $M_{\odot}\text{yr}^{-1}\text{kpc}^{-2}$) evolves as a function of radius and cosmic time: (left) without halo accretion; (right) with halo accretion. The five rows correspond to different gas fractions, in order from the top, $f_{\text{gas}} = 20\%, 40\%, 60\%, 80\%, 100\%$. All models commence with a short-lived, disk-wide “starburst” phase that moves outwards for the first 300 Myr. These stochastic starbursts are defined by SFR surface densities that exceed the average rate by more than an order of magnitude (i.e., blue regions). The higher $f_{\text{gas}}$ models sustain a higher SFR surface density over a larger radial range for both models, as do the halo accretion models in general when compared to the matching nonaccreting model.</figcaption></figure>

mechanism by which the early disk building could be fed, simulating disk formation with cold ($V_{\text{circ}}/\sigma_{\text{LOS}} \gtrsim 3$) stellar kinematics—such as the thick disk of the MW at epochs corresponding to $z > 2$—continues to be very challenging (O. Agertz et al. 2009; A. Pillepich et al. 2019). While more recent attempts appear promising (e.g., T. Tamfal et al. 2022; M. Kohandel et al. 2024), the formation of bimodal structures with thick and thin disks is not a natural outcome of these simulations.

Historically, the first thick disk formation scenarios focused on the response of a preexisting, thin stellar component to various heating mechanisms: the accretion of stars during a merger (A. C. Robin et al. 1996; M. G. Abadi et al. 2003), the heating of the disk by an interaction or merger event (P. J. Quinn et al. 1993; I. R. Walker et al. 1996; Á. Villalobos & A. Helmi 2008; Y. Qu et al. 2011), and the diffusion of stellar orbits by massive star clusters (P. Kroupa 2002; P. Assmann et al. 2011) or by radial migration (S. R. Loebman et al. 2011). Some studies also investigated the importance of a dissipative gaseous component in limiting the heating of a preexisting stellar disk (see, for example, B. P. Moster et al. 2010, 2011).

After the first structures were identified in high-redshift disks (e.g., B. G. Elmegreen & D. M. Elmegreen 2005; R. Genzel et al. 2008; K. L. Shapiro et al. 2008), the idea that massive gas clumps could provide a means of thickening the stellar component by diffusion of orbits was proposed by

F. Bournaud et al. (2009). These authors carried out sticky particle simulations rather than more realistic hydrodynamic simulations. As simulations of primordial disks improved (A. Dekel et al. 2009), the proposal that thick disks could arise from stars forming in a thick turbulent ISM could finally be tested. Other work suggested that stars forming in a turbulent, kinematically hot ISM generate a thick stellar disk, which can be further heated by mergers or interactions with satellites (G. S. Stinson et al. 2013; X. Ma et al. 2017; J. C. Bird et al. 2021; F. van Donkelaar et al. 2022). Here we challenge this claim by showing that the effect is not seen in higher-resolution *Nexus* simulations. Another scenario comes from X. Meng & O. Y. Gnedin (2021) who, using cosmological simulations, find that it is the accumulation of the successive generations of stars forming in thin disks of different orientations that gradually produces a thick disk.

Previous scenarios have generally focused on providing a way to thicken a disk, or have assumed the stars are born thick. However, the MW thick disk has important properties that were recently discovered, particularly from large spectroscopic surveys. First, the scale height of the thick disk is not constant, but increases from the most metal-poor and least α-enhanced stars to the most metal-poor and most α-enhanced (the so-called mono-abundance populations, see J. Bovy et al. 2012b; J. T. Mackereth et al. 2017). The α-rich disk is radially compact, with a scale length of about 2 kpc, which seems to be independent of the mono-abundant population (see J. Bovy
et al. 2012b, their Figure 5). It is unclear whether the thick disk has a synchronized, monotonic, chemical evolution, with little radial dependence (M. Haywood et al. 2013, 2015) or, conversely, an inside-out formation leading to a (steep) radial metallicity gradient (e.g., B. Ratcliffe et al. 2025). The former case is expected to lead to a tight age-metallicity relation, as observed in several studies (M. Haywood et al. 2013; M. Xiang & H.-W. Rix 2022; L. Casamiquela et al. 2024), while the latter implies uncorrelated age-chemistry distributions (A. Miglio et al. 2021; F. Anders et al. 2023). All these features constitute further constraints for any model of early disk formation in the Milky Way. The aim of our paper is to revisit some of these ideas and to explore a new mechanism for disk thickening.

The structure of the paper is as follows. In Section 2, we show how controlled simulations link to cosmological simulations, summarize the Nexus models and provide links to the simulation movies. In Section 3, we discuss the physical processes seen in the simulations, with a focus on the processes that lead to baryon sloshing. In Section 4, we ask whether the rapid growth of instabilities seen in the simulations are supported by theory. In Section 5, we discuss the implications of our new work and suggest future studies based on kinematic and morphological signatures. In Section 6, we place our new work in context and look forward to future improvements in the Nexus framework.

## 2. Nexus Simulations in Brief

### 2.1. Controlled Simulations in the Context of Cosmological N-body Simulations

In a seminal paper, H. J. Mo et al. (1998) introduced a simple theoretical framework describing how galaxy disks form from cooling hot gas within virialized dark matter halos, with their properties largely governed by angular momentum conservation (see also the updated treatment in J. Stern et al. 2020). This work laid the foundation for semianalytic models (SAMs) of galaxy formation. An unrecognized assumption in these models is that the circumgalactic medium (CGM) becomes structurally and dynamically stable (i.e., virialized) out to the virial radius at some early epoch. However, this assumption cannot hold in a strict sense, as the models also require that hot gas continues to cool within some inner radius.

Modern cosmological simulations diverge significantly in their predictions of CGM properties for Milky Way-like galaxies, mostly because of different implementations of star formation driven galactic winds (q.v. J. Stern et al. 2021; R. J. Wright et al. 2024). Controlled simulations are superior to cosmological simulations in tracing the evolution of galactic winds; the former include more realistic processes at much higher spatial resolution (q.v. B. Sike et al. 2025). Observationally, absorption-line studies of H I and heavy elements along quasar sightlines suggest that the CGM began to settle into a more stable state by around $z \gtrsim 2$ (H.-W. Chen 2012), consistent with CGM models showing self-similar behavior since $z \sim 4$ (A. Pallottini et al. 2014). This agreement supports our working assumption for comparisons with controlled simulations, such as those in the VINTERGATAN suite (O. Agertz et al. 2021). Specifically, we adopt characteristic CGM properties from VINTERGATAN progenitors at $z \approx 3$, a redshift by which the CGM appears largely stabilized (V. Rufo-Pastor et al. 2025, in preparation). In future work,

we will compare the evolutionary trajectories produced by these complementary approaches to modeling galaxy evolution.

### 2.2. Models and Movies

In earlier papers, we presented our new high-resolution (parsec scale) simulation framework (Nexus), which is ideally suited to studying how turbulent disks evolve dynamically (T. Tepper-García et al. 2024). These ecosystems manifest intermittent and long-lived behavior, including disk-halo interaction, bars and spiral arms, and bulge formation (J. Bland-Hawthorn et al. 2024). We were able to generate “bar-like” phenomena for the first time at all gas fractions, even in fully gas-dominated turbulent disks, including a remarkable “radial shear flow” with a characteristic signature that sets up before the central bulge is formed. As for the stellar disks, these bar-like signatures only arise in disks that dominate the local gravitational potential (J. Bland-Hawthorn et al. 2023).

The choice of parameters for our progenitor analogs is described in J. Bland-Hawthorn et al. (2024). This present work is focused on objects with halo masses of order $10^{11} M_{\odot}$, a dynamical mass that is characteristic of the recent JWST discoveries at $z \sim 1-6$ (L. Costantin et al. 2023; Y. Guo et al. 2023; Z. A. Le Conte et al. 2024). This is the expected mass of a Milky Way progenitor at $z \approx 3$ (Figure 1 in J. Bland-Hawthorn & O. Gerhard 2016). In high-redshift galaxies, the observed high gas fraction and enhanced gas surface densities (L. J. Tacconi et al. 2020) lead to elevated star formation rates, and the output energy and momentum couple with the gas to drive the turbulence (O. Agertz et al. 2009; E. Jiménez et al. 2023). The adopted star/galaxy formation physics is presented in O. Agertz et al. (2013) and O. Agertz et al. (2021). The star formation activity feeds back sufficient energy and momentum to maintain the turbulent support. (We note, in passing, that turbulent gas-rich systems have also been explored where the energy injection arises from self-gravity rather than star formation and feedback; see J. Fensch et al. 2023.)

Our first step is to consider an isolated galactic ecosystem in dynamical equilibrium, both with and without smooth accretion from the ambient hot corona. All computations were carried out with the RAMSES N-body/hydrodynamics code (R. Teyssier 2002) at benchmark ($N_{\text{io}} \sim 10^7$ elements) and high ($N_{\text{hi}} \sim 10^8$ elements) resolution, including star formation and metal production (see Table 1); the number of effective gas “cells” is a factor of 10 higher in both cases. In defining a Milky Way progenitor, we adopted a model with three key components: a live dark matter halo, a massive stellar/gaseous disk, and a hot coronal gas filling the live dark matter halo, which serves to supply the disk with a smooth flow of accreting gas after it cools.

For our Milky Way progenitor, we adopted halo parameters of $R_{\text{vir}} \approx 40$ kpc and $\log M_{\text{vir}}/M_{\odot} \approx 11$. Two other key parameters were (i) the disk mass fraction, $f_{\text{disk}}$, which determines whether the disk baryons dominate the underlying dark matter halo,

$$f_{\text{disk}} = \left( \frac{V_c, \text{disk}(R_e)}{V_c, \text{tot}(R_e)} \right)_{R_e=2.2R_{\text{disk}}}^2 . \quad (1)$$

Here, $V_c(R)$ is the circular velocity at a radius $R$, $R_{\text{disk}}$ is the exponential disk scale length, and $R_e = 2.2R_{\text{disk}}$ is the traditional scale length adopted in studies of disks as it
<table><caption>**Table 1**
Overview of Galaxy Models of the Initial Parameters</caption><thead><tr><th>M<sub>s</sub><br>(10<sup>8</sup> M<sub>☉</sub>)</th><th>M<sub>gas</sub><br>(10<sup>8</sup> M<sub>☉</sub>)</th><th>R<sub>disk</sub><br>(ckpc)</th><th>f<sub>disk</sub></th><th>f<sub>gas</sub></th><th>Label</th><th>High</th><th>Benchmark</th></tr></thead><tbody><tr><td>112</td><td>0.0</td><td>1.8</td><td>0.5</td><td>0.0</td><td>fd50_fg00_nac</td><td>0</td><td>1</td></tr><tr><td>88.8</td><td>22.2</td><td>1.8</td><td>0.5</td><td>0.2</td><td>fd50_fg20_nac</td><td>2</td><td>1</td></tr><tr><td>66.6</td><td>44.4</td><td>1.8</td><td>0.5</td><td>0.4</td><td>fd50_fg40_nac</td><td>2</td><td>3</td></tr><tr><td>44.4</td><td>66.6</td><td>1.8</td><td>0.5</td><td>0.6</td><td>fd50_fg60_nac</td><td>2</td><td>1</td></tr><tr><td>22.2</td><td>88.8</td><td>1.8</td><td>0.5</td><td>0.8</td><td>fd50_fg80_nac</td><td>0</td><td>1</td></tr><tr><td>0.0</td><td>112</td><td>1.8</td><td>0.5</td><td>1.0</td><td>fd50_fg100_nac</td><td>0</td><td>1</td></tr><tr><td>88.8</td><td>22.2</td><td>1.8</td><td>0.5</td><td>((0.2))</td><td>fd50_fg20_ac</td><td>0</td><td>1</td></tr><tr><td>66.6</td><td>44.4</td><td>1.8</td><td>0.5</td><td>((0.4))</td><td>fd50_fg40_ac</td><td>0</td><td>1</td></tr><tr><td>44.4</td><td>66.6</td><td>1.8</td><td>0.5</td><td>((0.6))</td><td>fd50_fg60_ac</td><td>0</td><td>1</td></tr><tr><td>22.2</td><td>88.8</td><td>1.8</td><td>0.5</td><td>((0.8))</td><td>fd50_fg80_ac</td><td>0</td><td>1</td></tr><tr><td>0.0</td><td>112</td><td>1.8</td><td>0.5</td><td>((1.0))</td><td>fd50_fg100_ac</td><td>0</td><td>1</td></tr></tbody></table>

Note. The DM host halo properties ($M_{\text{halo}} = 10^{11} M_{\odot}$; $R_{\text{vir}} = 37$ ckpc; $r_s = 9.2$ ckpc) are identical across models, and they are consistent with a MW-progenitor analog at $z \approx 3$. Table columns are as follows: (1) preexisting disk stellar mass (vertical scale height = 250 pc for all models); (2) disk gas mass; (3) disk scale length; (4) disk-to-total mass ratio (Equation (1)); (5) gas to total disk mass fraction at $t = 0$ (Equation (2))—the double brackets indicate accreting gas from a hot corona with a mass $M \approx 5.5 \times 10^9 M_{\odot}$ at $t = 0$; note that the hot halo mass is not accounted for in $f_{\text{gas}}$. (6) Model designation. (7)–(8) Number of high and benchmark resolution models with different random seeds.

corresponds approximately to the “turnover radius” in the rotation curve.<fnref n="6" /> (ii) The gas mass fraction within the same radial scale:

$$f_{\text{gas}} = \left( \frac{M_{\text{disk,gas}}}{M_{\text{disk}}} \right)_{R_s = 2.2 R_{\text{disk}}}, \qquad (2)$$

where $M_{\text{disk}}$ is the total disk mass and $M_{\text{disk,gas}}$ is the gas contribution.

The values for all parameters are given in Table 1. The simulations are grouped into different categories—benchmark resolution, high resolution, and repeated cases with different random seeds. There were six types of simulations distinguished by the different initial gas fractions ($f_{\text{gas}} = 0\%, 20\%, 40\%, 60\%, 80\%, 100\%$); these are run at benchmark resolution for a total of 2 Gyr because we are focused on the high-redshift Universe ($z \gtrsim 3$).

In all models, the gas fraction declines as more of the mass is locked up in stars (see Figure 3 in J. Bland-Hawthorn et al. 2024). For models with $f_{\text{gas}} < 100\%$, we included a preexisting (“old”) stellar disk such that $f_{\text{old}} = 1 - f_{\text{gas}}$. For models with $f_{\text{gas}} \approx 100\%$, we include a tiny preexisting stellar population ($f_{\text{old}} \approx 1\%$) as a useful tracer of the galaxy’s gas-driven evolution. The total mass of stars produced in 2 Gyr from low to high gas fraction spans $2.6-6.2 \times 10^9 M_{\odot}$ for the accreting halo models, and $0.6-4.0 \times 10^9 M_{\odot}$ for the nonaccreting models. The baryon mass was preserved in both the accreting and nonaccreting halo models. In the former case, the total baryon mass of the disk increases with time. Note that the initial scale length of the disk (both gas and stars) was roughly maintained at $R_{\text{disk}} \approx 1.8$ kpc across models. The above parameters are broadly consistent with arguably the best Milky Way progenitor analog to date, the object CEERS-2112 at a photometric redshift of $z \approx 3.0$ (L. Costantin et al. 2023).

In Table 1, the filename convention is `fdXX_fgYY_nac` for “no halo accretion” models, and `fdXX_fgYY_ac` for “halo accretion” models. Here, `XX` is the disk mass fraction percentage and `YY` is the gas fraction percentage. If the simulation has multiple versions generated by different

random seeds, the file is referred to as `fdXX_fgYY_nac`. `ZZ`, where `ZZ` is the version number. At our website,<fnref n="7" /> we provide a series of movies (animations) that show the evolution of each of our model galaxies. There are two types of animations: (1) “on-the-fly” animations; and (2) “post-processed” animations. The “on-the-fly” animations were created at simulation runtime; these have a high time resolution ($\delta t \approx 1$ Myr), and they show the evolution of the newly formed stars, of the gas density, and of the gas temperature on a face-on projection. The “postprocessed” animations have been created from the simulation outputs; they have a lower time resolution ($\delta t \approx 10$ Myr), and they show the evolution of the surface density of all disk components (gas, preexisting stars, and—if available—newly formed stars, gas) along three orthogonal projections.

## 2.3. Reference Frames

An idealized, multicomponent galaxy system is made up of several distinct contributions to the gravitational potential. Since the initial galaxy is set up as an equilibrium figure using AGAMA (E. Vasiliev 2019), this inertial reference frame conserves momentum in all components, for both gas, stars, and dark matter. In N-body simulations, the different mass components can become separated, depending on the nature and dynamics of the dark matter halo. In one of their examples, R. Joshi & L. M. Widrow (2024) plot the separate trajectories for the center of mass of the halo, bulge, and disk. Once the baryonic component starts to slosh, since each subsystem has its own internal structure, the restoring force and the drag imposed by dynamical friction can be very different.

As the gas fraction increases, the system experiences stronger gravitational potential fluctuations, primarily driven by large-scale turbulent energy sustained by ongoing star formation. In particular, momentum recoil from sporadic galaxy-wide starbursts adds significantly to these background fluctuations. As we demonstrate, the cumulative effect is a pronounced sloshing of the baryonic matter within the overall gravitational potential well. We find it useful and revealing to

<footnote n="6">⁶ This radius contains about two-thirds of the disk’s total mass.</footnote>

<footnote n="7">⁷ http://www.physics.usyd.edu.au/turbo_discs/</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 3.** The evolution of the azimuthally averaged, radial profile of the disk potential $\Phi_{\text{bary}}(R, t)$ with cosmic time, normalized to the starting potential $\Phi_{\text{bary}}(R, 0)$ encoded in color for (top) no halo accretion and (bottom) halo accretion. Blue tracks are early in cosmic time and red tracks are later times, as indicated. For $f_{\text{gas}} < 0.5$, the disk potential is fairly constant but, above this limit, the loss of gas mass in circulating winds leads to a substantial weakening of the disk potential. Note that the vertical axis is a ratio of gravitational potentials, such that a curve moving downwards reflects a weaker potential.</figcaption></figure>

define two distinct centers that are tied to two different reference frames, with the first being defined by the minimum of the baryonic gravitational potential:

$$ \nabla\Phi_{\text{bary}}(R, \phi, z) \approx 0, \quad \nabla^2\Phi_{\text{bary}}(R, \phi, z) > 0, \quad (3) $$

with the Hessian matrix of partial second derivatives for $\Phi_{\text{bary}}$ being positive definite. This baryonic center ($x_{\text{bary}}$, $y_{\text{bary}}$, $z_{\text{bary}}$) is useful because it happens to coincide with the peak of the light (surface density) distribution throughout. This is most often how observers perceive the center of a galaxy.

The second reference frame defines the minimum of the total potential:

$$ \nabla\Phi_{\text{tot}}(R', \phi', z') \approx 0, \quad \nabla^2\Phi_{\text{tot}}(R', \phi', z') > 0, \quad (4) $$

with all partial second derivatives for $\Phi_{\text{tot}}$ being positive definite. This center of the total potential ($x_{\text{tot}}$, $y_{\text{tot}}$, $z_{\text{tot}}$) is not traced by light, typically, in sloshing disks, only inferred from the disk kinematics. In practice, the minima defined in Equations (3) and (4) are identified after convolving the gravitational field at each timestep with a window function defined by 50 nearest-neighbor particles. Given that the central baryons dominate over dark matter, we anticipate that both minima have similar sloshing amplitudes with respect to an inertial frame.

### 3. The Origins of Sloshing: Simulations

#### 3.1. Unbinding the Disk

In Figure 3, we show how the radial profile of the simulated disk’s gravitational potential $\Phi_{\text{bary}}(R, t)$ evolves with time $t$ as a function of the disk gas fraction $f_{\text{gas}}$. For the “no halo accretion” models, there is only a 15%–20% change for $f_{\text{gas}} < 0.5$. In the halo accretion models, there is substantial accretion to the center at low $f_{\text{gas}}$, but the overall radial profile is relatively unchanged for $f_{\text{gas}} < 0.5$. At higher $f_{\text{gas}}$, the gravitational potential decreases rapidly over the first 500 Myr for both accretion cases. Up to half of the gravitational potential is lost due to gas launched into the halo by strong disk-wide winds. (The same process has been discussed in the context of dwarf galaxies by H. Zhao 2002.) Over the next

<figure><img src="image_8.png" /><figcaption>**Figure 4.** The total gravitational potential for a *Nexus* galaxy simulation with $f_{\text{disk}} = 0.5$ and $f_{\text{gas}} = 0$. The equipotential surfaces are shown for the inner 10 kpc (comoving) diameter.</figcaption></figure>

1.5 Gyr, there is a slow reaccretion to the inner disk, but most of it is recirculated to the outer disk through recycling in the galactic halo. The slow spin of the halo ensures the gas acquires angular momentum. We will explore the physics of this process in a later paper.

A great deal of energy is needed to “unbind” the disk in this way, but the feedback requirements through active star formation are consistent with high-redshift ALMA and JWST observations. We can obtain a rough estimate from the following formula for the binding energy:

$$ E_{\text{bind}} = -\frac{3}{5}\left(1 + \frac{q^2}{2}\right)\frac{GM_{\text{disk}}^2}{R_{\text{disk}}}, \quad (5) $$

assuming the disk dominates the local gravitational potential. The factor $q$ is the degree of flattening of the oblate spheroid along the polar axis, or about $q \approx 0.7$ from Figure 4. For convenience, we modify our estimate in terms of the disk’s local surface density $\Sigma(R) = \Sigma_0 \exp(-R/R_{\text{disk}})$ such that

$$ E_{\text{bind}} \approx -3\pi^2 G \Sigma_0^2 R_{\text{disk}}^3. \quad (6) $$
In practice, this is an upper limit because the unbound material does not need to reach infinity.

Alternatively, the timescale $\Delta\tau$ needed to partially unbind (loosen) the disk is essentially the time it takes to lose half of its mass or, equivalently, its surface density (e.g., J. G. Hills 1980). Thus

$$ \Sigma(R, t) = \Sigma(R) - \dot{\Sigma} \Delta\tau, \quad (7) $$

where $\dot{\Sigma} = \epsilon\dot{\mu}_*$, for which $\dot{\mu}_*$ is the SFR surface density, and where $\epsilon$ (=0.1) is the coupling efficiency of the feedback energy and radiation with the gas (O. V. Egorov et al. 2023, see their Section 5.2). Here, $\Delta\tau$ is the duration of the starburst,

$$ \Delta\tau = \frac{\Sigma(R)}{2\epsilon\dot{\mu}_*(R)}. \quad (8) $$

The SFR surface density $\dot{\mu}_*$ is the quantity shown in Figure 2 averaged over annular bins at each radius. Therefore, after rearranging Equation (8) and inserting the simulated surface density profile $\Sigma(R)$,

$$ \dot{\mu}_*(R) = 27 \left( \frac{\Sigma_0 e^{-R/R_{\text{disk}}}}{5.45 \times 10^8 M_{\odot} \text{ kpc}^{-2}} \right) \left( \frac{0.1}{\epsilon} \right) \times \left( \frac{100 \text{ Myr}}{\Delta\tau} \right) M_{\odot} \text{ kpc}^{-2} \text{ yr}^{-1}. \quad (9) $$

This is about $\dot{\mu}_* = 0.3 \pm 0.2 M_{\odot} \text{ kpc}^{-2} \text{yr}^{-1}$ at $R \approx 8$ kpc where the disk is more reactive and easily unbound, consistent with Figure 2, where sustained bursts over 100 Myr are a regular occurrence, hence the choice of normalization for $\Delta\tau$ in the above equation. The simulated bursts can be much higher but these are occurring in compact regions and on shorter timescales than $\Delta\tau$. The inferred time-averaged rates are sufficiently impulsive to (partially) unbind the disk, particularly at large galactocentric radius. Here, we have taken the $R_{\text{disk}} (=1.7 \pm 0.3 \text{ kpc})$ and $\Sigma_0$ normalizations from our simulated models (see below). At these early times, the most powerful starbursts create of order $3 \times 10^9$ stars in a $\Delta\tau$ interval across the inner kiloparsec of the disk.

We see that the required SFR surface density is achieved for all halo accretion models in the first 500 Myr. But at low $f_{\text{gas}}$, the disk potential is dominated by the existing stellar populations and these remain unmoved by feedback processes. It is only at high values of $f_{\text{gas}}$ that the disk becomes marginally unbound.<fnref n="8" /> The consequences of the unbinding are discussed in the next section.

### 3.2. Stochastic Fluctuations and Increasing Entropy

In Figure 3, the unbinding of the disk at high $f_{\text{gas}}$ appears to show a smooth evolution but these profiles are radially averaged. Locally, the unbinding is highly impulsive. Because the disk mass is lost through radiation pressure and entrainment in winds, the stellar disk becomes inflated vertically and radially. Stars are now oscillating vertically and radially in a weakened potential while retaining the same orbital energy. As

we see later, this is not a smooth outward progression over hundreds of megayears. The stochastic nature of the starburst-driven fluctuations is crucial to the inflated disk remaining in that state.

In fact, we note that for high $f_{\text{gas}}$ models, the disk reaccretes at later times. One might imagine that the process is reversible if the gravitational potential evolved smoothly, allowing the orbital properties of stars to evolve adiabatically to their former state. But this can never happen. The inflated stellar orbits are progressively more radial during expansion (e.g., S. Yu et al. 2021). The randomized radial phases contribute to the overall increase in entropy and this is an irreversible process.

There is an interesting precedent for the disk inflation process in the context of dwarf galaxies (H. Zhao 2002). A more extensive and related literature focuses on unbinding the NFW cusp predicted in dark-matter-only simulations, but never observed. For example, A. Pontzen & F. Governato (2012) show how fluctuations in the central gravitational potential $\Phi_0$ irreversibly transfer energy into collisionless particles, thus removing any dark matter “cusp” (J. F. Navarro et al. 1996). They consider blowout-recollapse cycles driven by a central starburst (see also G. Ogiya & M. Mori 2011; C. Nipoti & J. Binney 2015). By analogy, consider a star that is oscillating vertically with orbital frequency $\nu$ under simple harmonic motion

$$ z(t) = a(t)\cos[\nu(R)t + \psi(t)], \quad (10) $$

where the orbital amplitude $a$ and phase $\psi$ can undergo slow or rapid changes. The orbital frequency can be determined from

$$ \nu^2(R) \approx 4\pi G\Sigma(R)/z_{\text{disk}}, \quad (11) $$

where $z_{\text{disk}}$ is the disk’s vertical scale height. (The formal definition is given below.) The key point is that a prompt change (i.e., faster than the orbital period) of $\nu$ at a fixed radius due to a sudden blowout and recollapse is now associated with a different amplitude $a_1$ and phase $\psi_1$ along the star’s orbit compared to the original amplitude $a_0$ and phase $\psi_0$ (see S. Yu et al. 2021). According to A. Pontzen & F. Governato (2012), the change in energy has an explicit dependence on the orbital phase, such that

$$ \Delta E = \left[ \left( \frac{\nu_1}{\nu_0} \right)^2 - 1 \right] \sin^2 \psi_0, \quad (12) $$

where $\nu_0$ and $\nu_1$ are the vertical frequencies immediately before and after blowout or recollapse. This tends to extract energy from a star’s orbit at the blowout phase and return energy to the orbit at recollapse, but the overall effect is an energy gain due to the randomized orbital phases.<fnref n="9" /> As demonstrated by an idealized simulation (A. Pontzen & F. Governato 2012), the irreversibility thus arises not from dynamical differences but statistical differences between the blow-out and recollapse events. Any amount of reaccretion or recollapse does not return a disk to its (vertically) thinner state.

<footnote n="8">⁸ There is a strong similarity here with the dissolution of open star clusters (J. Bland-Hawthorn et al. 2010). At the cluster’s birth, the hot young stars drive the gas away and, if the surviving gas dominates the cluster’s gravitational potential, the cluster becomes unbound. This is the fate of most open clusters since star formation efficiencies are consistently low in most environments, as observed in our models.</footnote>

<footnote n="9">⁹ The analogy here is pushing a child on a swing. If the applied force is introduced slowly in phase with the swing, it can be made to swing higher or lower. If the applied force is out of phase, the swing motion is disrupted and becomes erratic. In the Pontzen mechanism, in general, the swing gains energy in a random walk; the swing’s energy can never be reduced to zero.</footnote>
<figure><img src="image_9.png" /><figcaption>**Figure 5.** The upper figure in each panel shows evolution of the dark matter density profile as a function of radius *R* for the halo accretion models. The lower figure in each panel presents the ratio of the post-starburst to original density profiles. At the highest gas fractions, the central NFW cusp is suppressed by a factor of 2, but not destroyed. In all plots, the color coding indicates the time evolution of the density profile, as indicated, with the reddest profile indicating the end of the simulation.</figcaption></figure>

Figure 5 investigates whether the density cusp residing in the NFW dark matter halo, already weakened by the dominant central disk, can be further softened or removed by the Pontzen mechanism. The evolution of the density profile (upper plots) is shown for the usual gas fractions. Many of the time curves overlap. The effect is really only detectable at the highest gas fractions, in the sense that the ratio of the density profile, after being divided by the starting profile, shows a clear depression at late times near $R=0$; even then, the suppression is at worst a factor of 2. So total destruction of the NFW cusp is not

achieved, unlike the case for dwarf galaxies with their much weaker central potentials. This finding is compatible with the absence of cusps at our targeted halo mass ($10^{11} M_{\odot}$) in modern cosmological simulations of galaxy formation (T. K. Chan et al. 2015; E. Tollet et al. 2016).

### 3.3. Sloshing Behavior

In the Nexus simulations, the powerful starburst events drive turbulence across the gas-rich disk through coupling of the mechanical and radiative luminosity arising from star
<figure><img src="image_10.png" /><figcaption>**Figure 6.** Nexus simulation for an $f_{disk} = 0.5$ stellar disk with no cold gas component ($f_{gas} = 0\%$). The top left and bottom left panels are the $x-y$ and $x-z$ slosh projections, respectively, at a fixed timestep. The baryon sloshing is defined in terms of the separation of the minimum of the total gravitational potential ($\Phi_{tot}(0, 0)$ indicated by the cross at the origin) and the baryon potential ($\Phi_{bary}(0, 0)$). (The boxed information concerns the disk asymmetry and concentration, a topic discussed in our companion paper, J. Bland-Hawthorn et al. 2025, in preparation). The physical scales of the three circles are indicated in the lower key). The top right and bottom right panels are the $x-y$ and $x-z$ slosh projections, respectively, but highly magnified compared to the LHS panels. A small amount of “baryon sloshing” ($\lesssim 50$ pc) is expected due to the discreteness of the disk, as demonstrated. The color coding indicates the baryon centroid shifts with the passage of time, as indicated.</figcaption></figure>

formation. These scattered events, both in time and location, drive global oscillations throughout the disk potential.<fnref n="10" /> The global pressure modes driven by strong impulsive fluctuations arising in our simulations are not “normal” modes that lead to long-lived standing waves. In galaxies, these oscillations would decay in $\lesssim 30$ Myr without the turbulence maintained by ongoing feedback. Interestingly, L. M. Widrow et al. (2012) has used the neologism “galactoseismology” to describe global oscillation modes of the Milky Way,<fnref n="11" /> but this is a specific reference to the response of the stellar disk. These so-called van Kampen modes (J. Binney 2025) do not exhibit collective oscillations with a single frequency, as normal modes do, but rather form a continuum of phase-space oscillations. These modes can occasionally reinforce each other (e.g., Gaia phase

spiral; T. Antoja et al. 2018) or, alternatively, phase mix away any coherent behavior.

Small-scale baryonic oscillations with respect to the dark matter halo have long been recognized (q.v. M. Kuhlen et al. 2013; R. Joshi & L. M. Widrow 2024). These can be induced by mergers through momentum recoil, or excessively clumpy, spinning or tumbling dark matter halos. Furthermore, as recognized by J. Y. Lau & J. Binney (2021), an N-body model initialized from the distribution function of an equilibrium model quickly develops macroscopic fluctuations in excess of what is expected from Poisson noise. We confirm this behavior on 50 pc scales in Figure 6 for a Nexus simulation that incorporates a stellar disk component free of gas ($f_{gas} = 0\%$).

The baryon sloshing is defined in terms of the separation of the minimum of the total gravitational potential ($\Phi_{0,tot}$ indicated by the cross) and the baryon potential ($\Phi_{0,bary}$) given by Equations (4) and (3) respectively. In Figure 7, we repeat the exercise carried out in Figure 6 to show the sloshing behavior for $f_{gas} = 20\%, 40\%, 60\%, 80\%, 100\%$. The “postprocessed” movies at our website http://www.physics.usyd.edu.au/turbo_discs/ clearly show the baryon sloshing in the inertial frame. See, for example, the $x-y$ projection of the newly formed stars (left middle panel) in fd50_fg100_ac; the emerging stellar bar soon drifts away from the origin and follows a complex track with multiple loops.

In Figure 7, notice that the vertical migration is half that of the radial migration in all cases, including the $f_{gas} = 0\%$ case

<footnote n="10">¹⁰ There is a useful analogy here with the Sun. The turbulent motions in its convective zone create random pressure perturbations, which act like “kicks” to the surrounding plasma (G. Houdek & M.-A. Dupret 2015) The randomness of the turbulence excites a broad range of discrete modes simultaneously. The acoustic waves generated by convection can reinforce each other through resonance, leading to standing waves that persist. These standing waves correspond to eigenmodes of oscillation with specific frequencies and spatial patterns. Each mode is also damped over time by dissipative processes, but because of the constant churning of convection, new energy is continually injected into the modes, maintaining their amplitudes over long periods. Detectable oscillations induced by these p-modes form the basis of helioseismology and, for stars more generally, asteroseismology (C. Aerts 2021). Similar processes occur in gas-rich galaxies, but these excite a continuum of (van Kampen) modes rather than standing waves excited by discrete modes.</footnote>

<footnote n="11">¹¹ For a brief history, see J. Bland-Hawthorn & T. Tepper-García (2021).</footnote>
<figure><img src="image_11.png" /><figcaption>**Figure 7.** Nexus simulations in the presence of halo accretion for an $f_{\text{disk}} = 0.5$ stellar disk for different values of $f_{\text{gas}}$, as indicated. Left panel: $x - y$ (left column) and $x - z$ (right column) projections of the sloshing action. Note that the physical scale of each plot is growing with increasing $f_{\text{gas}}$. Right panel: $x - y$ (left column) and $x - z$ (right column) projections of the sloshing action. Note that the physical scale of each plot is fixed for all $f_{\text{gas}}$. The baryon sloshing is defined in terms of the separation of the minimum of the total gravitational potential $\Phi_{\text{tot}}(0, 0, 0)$ (indicated by a cross) and the baryon gravitational potential, $\Phi_{\text{bary}}(0, 0, 0)$. The color coding indicates the passage of time, as indicated by the RHS bar, extending from $z \approx 3.3$ to $z \approx 1.5$. Internal to each plot, the horizontal and vertical scales are always matched.</figcaption></figure>

in Figure 6. In Figure 4, the massive disk flattens the total gravitational potential along the vertical axis. The vertical restoring force is stronger than the radial force, and thus provides more resistance to the vertical drift of the baryons. The total (effective) gravitational potential (Figure 4) can be approximated with

$$ \Phi_0(R, z) \approx \Phi_1(0, 0) + \frac{1}{2}\kappa^2 R^2 + \frac{1}{2}\nu^2 z^2, \quad (13) $$

where the natural harmonic frequencies of the smooth potential are evaluated at $z = 0$ in both instances. For simplicity, we have dropped the $\phi$ ordinate even though the simulated

potential is rarely axisymmetric. Thus

$$ \kappa^2(R) = \left( \frac{\partial^2 \Phi_0}{\partial R^2} \right)_{z=0}, \quad (14) $$

$$ \nu^2(R) = \left( \frac{\partial^2 \Phi_0}{\partial z^2} \right)_{z=0}. \quad (15) $$

For small perturbations, the restoring force is

$$ F_{\text{rest}} = M_{\text{halo}} \frac{dV_{\text{disk}}}{dt} \approx -M_{\text{halo}}(\kappa^2 R \mathbf{i} + \nu^2 z \mathbf{k}), \quad (16) $$

where $M_{\text{halo}}$ is the total halo mass (dark matter + hot coronal gas), and $V_{\text{disk}}$ is the disk's velocity with respect to the total
<figure><img src="image_12.png" /><figcaption>**Figure 8.** The 2D slosh projections in Figure 7 can be rendered in 3D to reveal the erratic movement. The tracks all refer to the $f_{gas} = 100\%$ halo accretion case: (left) motion of $\Phi_{tot}(0, 0, 0)$ with respect to the inertial reference frame; (middle) motion of $\Phi_{bary}(0, 0, 0)$ with respect to the inertial reference frame; (right) motion of $\Phi_{bary}(0, 0, 0)$ with respect to $\Phi_{tot}(0, 0, 0)$. The slosh axes are in units of kiloparsecs in a comoving frame. Evidently, for the baryons, there are frequent major reversals over the 2 Gyr simulation. The color coding indicates the passage of time, from blue to red over this timespan.</figcaption></figure>

gravitational potential. Thus, as expected, since $\nu > \kappa$, the sloshing perturbations for all simulations is more compressed along the vertical axis when compared with in-plane motions.

This fact has important implications for the physical consequence of sloshing. For example, the momentum recoil kicks that push the gas around are essentially isotropic to first order, with no preferred direction. While we do see occasional evidence for oscillatory behavior, the overall effect is more chaotic due to strong stochastic impulses acting in three dimensions. R. Joshi & L. M. Widrow (2024) incorporate an additional damping term to the restoring force, i.e.,

$$F_{\text{damp}} = -\beta M_{\text{halo}} V_{\text{disk}}, \quad (17)$$

to account for the dynamical friction from the unseen dark matter and hot baryonic halo. The factor $\beta \approx 4\pi G^2 \rho_{\text{halo}}(R)M_{\text{halo}}b/\sigma_{\text{halo}}(R)^3$ where the variable $b \approx \frac{1}{4}$ near $R \approx 0$ and goes to zero at infinity, and $\rho_{\text{halo}}$ and $\sigma_{\text{halo}}$ are the density and velocity dispersion of the halo in the vicinity of the disk. We arrive at

$$F_{\text{rest}} \approx -M_{\text{halo}}(\kappa^2 R \mathbf{i} + \nu^2 z \mathbf{k} + \beta V_{\text{disk}}). \quad (18)$$

This can explain why, in most cases, the disk centroid appears to settle to the $z=0$ plane of the total gravitational potential. Furthermore, in some of the simulations (e.g., $f_{gas} = 40\%$ in Figure 7), the baryonic center executes a near-circular wide orbit about the total potential, which may arise from the orbit-circularizing effect of dynamical friction (I.-G. Jiang & J. Binney 2000). We return to these issues in the next section.

In Figure 8 (left), the 2D slosh projections shown in Figure 7 are rendered in 3D to reveal the erratic (Brownian) nature of the movement. Here, we show the halo accretion case for $f_{gas} = 100\%$. Evidently, there are frequent major reversals over the 2 Gyr simulation. A gallery of radial and kinematic signatures for all $f_{gas}$ is presented in the Appendix.

# 4. The Origins of Sloshing

## 4.1. Theory versus Simulations

The dynamics of stellar disks is a remarkably complex process, with only restricted analytic cases that are well understood (J. Binney 2025; C. Hamilton et al. 2025). So to

what extent can we trust N-body simulations to trace the growth of instabilities and nonlinear perturbations realistically? One of the most powerful analytic theories to address this issue is the (local) shearing-sheet model presented by W. H. Julian & A. Toomre (1966) for a razor-thin 2D disk. Here, a perturbation in the distribution function (DF) can be followed in both density space and velocity space. A local disturbance generates a leading wake that is “swing-amplified” into a trailing wake. These density perturbations grow as a consequence of gravity, shear, and rotation (A. Toomre 1964), but eventually dissolve unless sustained by a perturbing (e.g., external) force. Even when the density enhancement has disappeared, its influence is seen in velocity space long afterwards but eventually fades from view due to phase mixing.

In the shearing-sheet approximation, the swing-amplified modes are analyzed in a frozen, corotating frame where both the global and local rotation (due to shear and vorticity) are removed. The model treats stars that are on near-circular orbits, although corrections can be made for mildly noncircular orbits (J. Binney 2020). In such a frame, all curvature is removed and spiral features can be treated as time-dependent plane waves, typically via a local density (or surface density) —gravitational potential pair of coupled functions.

## 4.2. *Jolting the Disk—The Importance of Self-gravity*

So how does the shearing sheet respond to a powerful jolt? This line of enquiry is revealing because it shows that an impulse experienced by a stellar disk is amplified by its self-gravity for all $Q_*(=\kappa\sigma_R/3.36G\Sigma_0)$ values relevant to Nexus galaxies. Here, $\Sigma_0$ is the local undisturbed surface density, $\kappa$ is the epicyclic frequency and $\sigma_R$ is the radial velocity dispersion. $Q_*$ is a dimensionless quantity that measures the local axisymmetric stability in the stellar disk against collapse or fragmentation; nonaxisymmetric growing modes are still allowed. The amplifying effect is strongest at intermediate wavelengths (see below), becoming weaker at long wavelengths, while short wavelengths are suppressed below a critical wavelength $\lambda_{crit} = 2\pi/k_{crit}$, for which $k_{crit} = \kappa^2/2\pi G\Sigma_0$. The characteristic epicyclic radius is roughly $\sigma_R/\kappa$. These physical dimensions are on kiloparsec
<figure><img src="image_13.png" /><figcaption>**Figure 9.** Evolution of wave amplitudes $A(t) = \tilde{\Sigma}_1(t)$ in the presence of an impulse, obtained by solving the iterative integral in Equation (23), where the tilde indicates a surface density amplitude normalized to evolution without gravity. The horizontal axis shows time evolution in units of $\kappa/\pi$, consistent with shearing-sheet analysis (W. H. Julian & A. Toomre 1966), starting at $\kappa t_i = -\frac{3\pi}{2}$. (Left) The response of a Mestel disk to a sharp impulse for $Q = 1.2, 1.3, 1.4$, demonstrating the increasing importance of self-gravity to the disk's response (J. Binney 2025). The red curves illustrate the weak disk response when self-gravity is neglected. (Right) This repeats the left panel but for a longer-wavelength disturbance; only wavelengths shorter than $k = k_{crit}$ are stabilized by the initial $Q$ value. The range of $Q$ values shown describes the Nexus disks for most of their evolution.</figcaption></figure>

scales, e.g., $k_{crit} \approx (2 \pm 1 \text{ kpc})^{-1}$, such that, for example, $\lambda_{crit} \sim 10 \text{ kpc}$ near the Solar circle today.

In the shearing-sheet approximation (W. H. Julian & A. Toomre 1966), the coordinate system is centered on a spiraling wave's corotation radius. In Galactic cylindrical coordinates $(R, \phi, z)$, individual stars have velocities $V = (V_R, V_\phi, V_z)$ and oscillation frequencies $\Omega = (\Omega_R, \Omega_\phi, \Omega_z) = (\kappa, \omega, \nu)$. Following J. Binney (2020), we introduce a Cartesian wavevector $\mathbf{k}$ that corotates around a point defined by the radius vector $\mathbf{R}$ such that

$$|\mathbf{k}| = |k_x \mathbf{i} + k_y \mathbf{j}| = k(t) = k_y \sqrt{1 + (2A(t - t_0))^2}, \quad (19)$$

where $k_y$ is the tangential wavenumber, and $A = -\frac{1}{2} \frac{d\Omega_\phi}{d \ln r}$ is Oort's constant that quantifies the shear rate. The wavevector ensures that the perturbation's wave crests track from leading ($t < t_0$, $k_x < 0$) to radial ($t = t_0$, $k_x = 0$) to trailing waves ($t > t_0$, $k_x > 0$), where the amplitude peaks near $t = t_0$. A suitable potential—surface density pair is

$$\Sigma(\mathbf{x}, t) = \Sigma_0 + \Sigma_1 \cos(\mathbf{k} \cdot \mathbf{x}) \quad (20)$$

$$\Phi_1(\mathbf{x}, z, t) = -2\pi G \frac{\Sigma_1}{k(t)} \cos(\mathbf{k} \cdot \mathbf{x}) e^{-k(t)/|z|}, \quad (21)$$

where $\mathbf{x}$ describes the $x-y$ plane, and $\Sigma_0$ and $\Sigma_1$ are the undisturbed and wave-disturbed density fields.

In Figure 9, we show the response of a razor-thin 2D disk to a strong local impulse for different values of the $Q$ parameter, from $Q = 1$ (marginally stable, high self-gravity) to $Q = 1.4$ (most stable, low self-gravity). We initiated the impulse at a time $t = t_i$ before $t = t_0$; the perturbed distribution function is $f_1(t_i) = 0$ at the start. The impulse is such that the surface density from the perturbation is $\tilde{\Sigma}_e(t, t_i) = \delta(t - t_i)/\kappa$. Time

is needed for the perturbation to build, but if that time is set too far in the past, the disturbance can phase mix away before it has a chance to build. From the so-called JT equation, the renormalized surface density evolution is given by

$$\tilde{\Sigma}_1(t) = \int_{t_i}^{t} \mathcal{K}(t, t') (\tilde{\Sigma}_1(t') + \tilde{\Sigma}_e(t')) dt' \quad (22)$$

$$= \mathcal{K}(t, t_i) + \int_{t_i}^{t} \mathcal{K}(t, t') \tilde{\Sigma}_1(t') dt'. \quad (23)$$

The JT kernel is

$$\mathcal{K}(t, t_i) = \frac{4\mathcal{F}(t, t_i)}{k(t)\exp(0.57Q^2b^2)}, \quad (24)$$

for which $b$ depends on $k_y/k_{crit}$, and $\mathcal{F}(t, t_i)$ depends on geometric factors involving $A$, $\kappa$, and $\Omega$ (J. Binney 2020). The first RHS term in Equation (23) arises from the impulse and the second integral term is the temporal, accumulated response to the disk's self-gravity. As observed in the simulations, the strong impulses excite a continuum of long-wavelength waves that carry energy and momentum throughout the stellar disk, the gas disk and the dark matter halo.

In the Nexus simulations, away from the galaxy's center, $Q$ is in the range 1.2–1.5, leading to swing-amplification factors $A$ as high as 10–60 compared to disks without self-gravity (see M. Aumer et al. 2016). J. Binney (2020) showed that the biggest amplification occurs around $\lambda/\lambda_{crit} \approx 2 \pm 1$ (see also A. Toomre 1981). Thus, it is entirely reasonable that a succession of disk jolts leads to a strong reaction from the large-scale stellar disk (e.g., spiral arm excitation, internal disk heating) on tens of kiloparsec scales (q.v. C. Hamilton et al.
<figure><img src="image_14.png" /><figcaption>**Figure 10.** Two different perspectives of the same Brownian motion track in a moderately flattened, spherical gravitational potential. The slosh coordinates (in kiloparsecs) are equivalent to those presented in Figure 7. The color scheme once again shows the evolution over 2 Gyr (blue → red). The track is similar in character to the baryon sloshing tracks in Figures 7 and 8.</figcaption></figure>

2025). Once the momentum-injecting starbursts die down, the disk ceases to slosh, as observed in the simulations.

### 4.3. Brownian Motion

We separate the phenomenon and consequences of baryon sloshing in the total gravitational potential into two parts, which can both be explored with the Fokker-Planck equation using a probability density function $P(R - R_0, t)$, for which $R$ is the radius vector from an initial position $R_0$. First, the baryon potential minimum moves around stochastically with what is effectively damped Brownian motion. We demonstrate this property with a simple model. Second, the internal response of the baryons is subject to energy loss and energy gain in different measures due to the changing underlying potential. In a later paper, we treat the internal disk evolution (heating and cooling) separately as a diffusion problem in action space.

For the simple model, there are “coasting” phases where the simulated galaxy can move on a trajectory, but with associated noise due to the stochastic nature of the starburst activity. The overall damping is due to the simulated galaxy moving in a deep potential well and dynamical friction due the disk’s motion in a live dark matter halo. On occasion, the entire galaxy can bob up and down with damped quasiperiodic motion before reverting to an erratic trajectory. Here we stress the Brownian motion behavior because it explains why the disk heating is essentially an isotropic process (Section 5.3).

To simulate the track of a particle, Langevin's equation is used to describe the motion subject to (overdamped) Brownian motion, i.e., $\gamma\dot{R} \approx F(R) + \eta(t)$. Here, $F$ is the continuously applied force (e.g., galactic potential gradient), $\eta$ supplies a random force term and $\gamma$ is a damping constant. In our N-body simulations, the damping is mostly due to dynamical friction, i.e., energy loss to the live dark matter halo. To add an extra layer of realism, we computed the baryon motion in a moderately flattened, spherical gravitational potential ($F = -\nabla\Phi_0$) that mimics the potential in Figure 4. This

weakly biases each step $dR$ along the radius vector toward the center and compresses the motion in $z$. An example of such a track is shown in Figure 10 from two perspectives; most realizations capture the essence of Figures 7 and 8. In the absence of external forces, we can mimic isotropic diffusion in three dimensions using a 3D Gaussian random variate such that $\eta(t) = G[R(t) - R_0, \sigma]$, where $|\sigma| \propto \sqrt{D_2t}$. For an ensemble, we find (i) the first moment is $\langle\eta_i(t)\rangle = 0$ and (ii) the second moment is $\langle\eta_i(t)\eta_j(t')\rangle = 2D_2\delta_{ij}\delta(t-t')$, where $\delta_{ij}$ is the Kronecker delta and $\delta$ is the Dirac delta function. The latter ensures the covariance of the (Wiener) noise process is uncorrelated. $D_2$ is a diffusion coefficient that we chose arbitrarily, along with the timesteps, to match the amplitude of motion in kiloparsecs seen in the N-body simulations.

## 5. Discussion: The Consequences of Sloshing

Here we discuss the main consequences of baryon sloshing, at least as we understand them today. The *Nexus* simulations are only run for a few billion years at present, long enough to explore the formation and evolution of the α-rich disk. In a later paper, we will trace the evolution of baryon sloshing over the lifetime of the Milky Way. Here, we look at how the vertical and radial scale lengths evolve and discuss in the context of the sloshing action. We also look at the predicted metallicity trends throughout the α-rich disk.

### 5.1. Disk Thickening

The formation of early thick stellar disks has been a challenge for theory and simulations alike. One class of models invoke stars forming in a turbulent, kinematically hot ISM to generate a thick stellar disk, which can be further heated by mergers or interactions with satellites (G. S. Stinson et al. 2013; X. Ma et al. 2017; J. C. Bird et al. 2021; F. van Donkelaar et al. 2022). However, this effect is not seen in our higher-resolution *Nexus* simulations. In Section 5.4, we present evidence *against* disk thickening arising from in situ star formation in a turbulent gas disk. Although star formation
<figure><img src="image_15.png" /><figcaption>**Figure 11.** The rms z-scale height for preexisting (old) stars as a function of radius *R* (in kiloparsecs) for each value of *f<sub>gas</sub>*. The left panel is for no halo accretion cases; the right panel is for halo accretion models. In all plots, the color coding indicates the time evolution, as indicated. The scale height of the α-rich thick disk today at *R*₀ = 8.2 kpc is indicated by the cross (J. Bland-Hawthorn & O. Gerhard 2016). Note that the vertical scale has a fourfold stretch compared to the horizontal axis. Gas-rich models are more able to reproduce the α-rich disk within 1.5–2 Gyr.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>**Figure 12.** (Top) The radial scale length for preexisting stars as a function of time (in megayears) for each value of $f_{\text{gas}}$ indicated in the key. The left panel is for no halo accretion models; the right panel is for halo accretion models. (Right) The same quantity but for newly formed stars as a function of time (in megayears) for each value of $f_{\text{gas}}$ indicated in the key. In all panels, the gray band indicates the radial scale length of the thick disk today, i.e., $R_{\text{disk}}^{T} = 2.0 \pm 0.2$ kpc (J. Bland-Hawthorn & O. Gerhard 2016).</figcaption></figure>

and turbulence are still involved, our mechanism is very different.

In Figure 11, the rate of disk heating increases dramatically with $f_{\text{gas}}$. This is shown for preexisting stars as these are a good tracer of the internal heating. (We present the matching plot for newly formed stars in Section 5.4.) The internal heating continues until the gas fraction drops and the star formation rate follows suit, at which point the internal heating becomes saturated. The high $f_{\text{gas}}$ models achieve the $\alpha$-rich thick disk values today. The $f_{\text{gas}} = 80\%$ and 100% models overshoot the measured values today, but this is likely to have taken place in the early Milky Way. Once the younger, massive thin disk forms, the preexisting $\alpha$-rich thick disk must have been compressed by the increased density of the $\alpha$-poor disk.

Another consequence of disk heating is the rise of the vertical stellar dispersion $\sigma_z$ in both preexisting and newly formed stars with cosmic time. The relevant figures are presented in our earlier paper (see Figures 6 and 20 in J. Bland-Hawthorn et al. 2024). In Section 5.3, we present a simple model for how this heating arises from the fluctuating gravitational potential. Evidence for quasi-isotropic heating is the roughly preserved ratios $\sigma_z/\sigma_R$ and $\sigma_\phi/\sigma_R$ across all radii

for the preexisting stars throughout the simulation. Initially, the situation is different for the newly formed stars because their radial distribution exhibits inside-outside growth (Figure 12), as observed in our movies. But the heating rate of the younger stars matches that of the preexisting stars once the formed stars extend across the preexisting disk (see Section 5.3).

In Figure 12, the evolution of the radial scale length is presented over the full simulation for the preexisting (top) and newly formed stars (bottom). To combat noisy data, we determine the disk scale length $R_{\text{disk}}$ from the radius enclosing half the disk's mass; we divide the latter by 1.67, which is appropriate for a 3D exponential disk. For the preexisting stars, we see the disk grow radially by no more than about 20% of its starting dimension. For newly formed stars, the growth is more dramatic. All halo accretion models achieve the radial scale length observed for the thick disk today (J. Bland-Hawthorn & O. Gerhard 2016; J. T. Mackereth et al. 2017), indicated by the gray band. For the nonaccreting models, a very high starting $f_{\text{gas}}$ value is required. In contrast, the vertical disk scale height grows by a factor of 1.5–8 depending on the model (Figure 11). This statement is not inconsistent with the
<figure><img src="image_17.png" /></figure>

**Figure 13.** (Left) The $f_{\text{gas}} = 100\%$ accreting halo case is repeated from Figure 11 and highlights the different heating phases due to rapid reversals along the sloshing track. In practice, we allow for 1% of the disk mass to be an existing stellar population so that we can trace the impact of the heating. The z-height is starting to saturate after 1.5 Gyr. (Right) This is a repeat of the LHS plot but for stars forming in the gas rather than preexisting stars. The heating rate is essentially the same. The bottom 10 tracks (100 Myr apart) are removed because of large intrinsic scatter; this is a consequence of small number statistics and early populations being born out of dynamical equilibrium. Note for both plots that the vertical scale has a fourfold stretch compared to the horizontal axis. The cross indicates the vertical scale height of the Milky Way's $\alpha$-rich disk at the solar circle today.

claim of isotropic heating: the young disk is initially very cold in the vertical direction ($\sigma_z \sim 5 \text{ km s}^{-1}$) compared to the radial direction ($\sigma_R \sim 60 \text{ km s}^{-1}$).

In Figure 13, we highlight a causal connection between the baryon sloshing motion and the progressive vertical heating of the disk. The vertical growth goes through different phases of slow and fast growth. With reference to the 3D sloshing track for this specific model (Figure 8), we should not be surprised that the overall heating is roughly isotropic when averaged over time. In what follows (Section 5.3), a dynamical argument is given in support of this statement.

## 5.2. Chemical Signatures

We provide preliminary results on the predicted stellar chemistry for the early thickened disk (see T. Tepper-García et al. 2024; O. Agertz et al. 2021, for details). Some of our “results” are not particularly compelling because, like all simulators, we set up the initial conditions for success. Specifically, the initial star formation rates were very high and this guaranteed the production of high-mass stars and supernovae under starburst conditions. These, in turn, produced $\alpha$-enhanced yields ensuring the first generations of stars are $\alpha$-enriched. This signature remains as some of the best supporting evidence that early disks were gas-rich (K. Freeman & J. Bland-Hawthorn 2002; J. Bland-Hawthorn & S. Sharma 2016). Simulators often present their “predicted” old $\alpha$-rich stars as a success of the model, but the first generations tend to emerge as a direct consequence of the initial setup conditions. Arguably, it is the subsequent generations that are more useful.

In particular, the left panel in Figure 14 shows a clear trend in [Fe/H] as time passes. This is significant because the thick disk has been described as a snap-frozen, mono-age stellar population (see Section 1 and K. Freeman & J. Bland-Hawthorn 2002), which is inconsistent with our model. Well defined age-metallicity trends have been discussed extensively in the

context of the $\alpha$-rich disk (M. Haywood et al. 2013; M. R. Hayden et al. 2014; J. T. Mackereth et al. 2019; S. Sharma et al. 2021; M. Xiang & H.-W. Rix 2022). Most recently, this is precisely what is seen for the $\alpha$-rich population in the GALAH/Gaia (C. L. Sahlholdt et al. 2022) and APOGEE/Gaia surveys (V. Cerqui et al. 2025), both in terms of timescale and metallicity range.

The right panel in Figure 14 show $R-z$ cross-sections of the thickened stellar disk at 2 Gyr ($z \approx 1.5$) after the disk formed at $z \approx 3.3$. The upper panel exhibits a clear vertical gradient in stellar age, with the oldest stars furthest from the plane. The lower panel shows a clear vertical stratification in stellar metallicity with the most metal-poor stars furthest from the plane. These general trends (albeit described by different functions of $R$ and $z$), including the flaring at large radius due to the declining disk potential, are well established for both the $\alpha$-poor and $\alpha$-rich disks (J. Bovy et al. 2012a; J. T. Mackereth et al. 2017).

In both panels, the radial trends are very weak, unlike what is observed in the mass-dominant $\alpha$-poor disk. Radial gradients have been searched for repeatedly in APOGEE data, for example, but never established (F. Anders et al. 2014; M. R. Hayden et al. 2014). Currently, we are unable to say whether the lack of a radial gradient reflects the thick disk's formation epoch, or its subsequent evolution. For example, was there an early radial gradient that was washed away by subsequent bar formation? The lack of radial structure makes it difficult to determine a likely birth radius for individual stars using chemical tagging arguments. This question requires more investigation.

We would argue that, on the evidence provided, the correct vertical age-metallicity trends and the nonexistent radial trends are successes of the baryon-sloshing model. More extensive analysis is required to establish whether the baryon-sloshing model has unique chemodynamic signatures that we can search for.
<figure><img src="image_18.png" /><figcaption>Figure 14. Predicted trends at 2 Gyr ($z \approx 1.5$) after disk onset for the halo accreting $f_{gas} = 100\%$ model, representative of the high $f_{gas}$ models in general. (Left) Age-metallicity trend in the simulated $\alpha$-rich disk in the simulated $R-z$ plane, i.e., within a vertical slice that is 100 pc thick along the line of sight. (Top right) Age distribution in the same plane showing a weak trend in radius but a stronger vertical gradient, such that the youngest stars lie near the plane. (Bottom right) Metallicity distribution in the same plane showing a weak trend in radius but a stronger vertical stratification, such that the most metal-rich stars lie near the plane. For the RHS figures, representative stellar density contours within the vertical slice are overlaid.</figcaption></figure>

<figure><img src="image_19.png" /></figure>

### 5.3. Internal Heating of the Stellar Disk

Contemporary discussions of Galactic heating tend to focus on either vertical or in-plane processes, or heat energy lost to the dark matter halo, although 3D disk heating has been discussed (S. Ida et al. 1993). This is because the discussion centers on an observed phenomenon and a specific process, e.g., radial abundance gradients and in-plane migration. But here we stress that baryon sloshing must lead to isotropic heating, at least to first order, such that all Galactic coordinates are affected simultaneously. Figures 11 and 12 are evidence enough that vertical and in-plane heating occur concurrently in the simulations. In a later paper, we return to a detailed discussion of 3D heating using a Fokker-Planck diffusion analysis in action space.

Instead, we use a simple model for how vertical heating of the stellar disk arises from sloshing in a time-dependent background potential. The connection between fluctuating gravitational fields and internal dissipation or heating is a classical problem in statistical physics (M. S. Green 1954; R. Kubo et al. 1957), and has also been examined in astrophysical dynamics (e.g., R. W. Nelson & S. Tremaine 1999; M. D. Weinberg 2001). To build intuition, consider a galactic disk oscillating vertically within an external potential. If the fluctuation $\delta\Phi$ has a zero mean around the equilibrium position ($z=0$), the equation of motion for a stellar mass $m$ is $m\dot{v}_z = F(z, t) = -\nabla\delta\Phi$. For simplicity, we set $m=1$ since the mass cancels out in the dimensionless form of the analysis. For a stochastic process, the diffusion coefficient is defined as

$$ \mathcal{D} = \frac{1}{2} \frac{d}{dt} \langle (v_z(t) - v_z(0))^2 \rangle \quad (25) $$

$$ = \frac{d}{dt} \int_{0}^{t} (t - \tau) C(\tau) d\tau, \quad (26) $$

where the force autocorrelation function is $C(\tau) = \langle F(0)F(\tau) \rangle$. As $t$ goes to infinity, the diffusion coefficient simplifies to

$$ \mathcal{D} = \int_{0}^{\infty} C(\tau) d\tau \approx \langle F^2 \rangle \tau_{\text{corr}}, \quad (27) $$

assuming $C$ decays over a correlation time $\tau_{corr}$ and has a variance $\langle F^2 \rangle$. The baryon sloshing decays over timescales of 400–700 Myr for the different models as the gas is converted into stars. We approximate $\langle F^2 \rangle = (\delta\Phi_z/Z)^2$, where $Z$ is the vertical distance over which the potential is varying. For $\Phi_z \approx \sigma_z^2$ and $\tau_{dyn} \approx Z/\sigma_z$, assuming the simulation time $T_o$ (2 Gyr) is much longer than $\tau_{dyn}$ ($\approx 600/2\pi$ Myr), we can approximate the change in the stellar velocity dispersion to be

$$ \frac{\delta\sigma_z^2}{\sigma_z^2} \approx \frac{2}{3} \left( \frac{\delta\Phi_z}{\Phi_z} \right)^2 \frac{\tau_{corr} T_o}{\tau_{dyn}^2}. \quad (28) $$

as a function of the change in the underlying gravitational potential.

The clearest demonstration of how stellar populations heat with time is presented in J. Bland-Hawthorn et al. (2024, Figure 6), drawn from the same data used in Figure 11. For the highest gas-fraction models, the preexisting stars increase their dispersions by factors of 1.5–2, requiring gravitational fluctuations of 10%–20%. The young stars, at least those born in a dynamically stable disk, have a colder intrinsic dispersion ($\sigma_z \approx 6–8$ km s$^{-1}$): these increase their dispersions by factors of 3–4 requiring fluctuations of 50%. This is an order of magnitude too high to be explained by the underlying dark matter potential or by stars being born out of equilibrium. Instead, we believe the extra heating arises from the large $\delta\Phi_z$ fluctuations due to the unbinding of the disk. We refer the reader to Figure 3 and the discussion in Section 3.1 where 30%–50% fluctuations are seen in the most extreme cases. How the unbinding works with the sloshing action is unclear
and deserves further study. This work has started and will be reported on in a future paper.

## 5.4. The Case Against Thick Disks ($|z| \approx 1$ kpc) from In Situ Star Formation

In Sections 5.1 and 5.3, we presented evidence for disk heating being a gradual accumulation of energy from baryon sloshing within the dark matter potential. Star formation feedback and turbulence are implicated in the thickening of disks, but *not* from stars being born in situ far from the cold molecular disk. This idea first arose from nonphysical, sticky particle simulations where there is no equation of state, and therefore no pressure balance, or even feedback from star formation processes to support the gas. In Figure 15, most stars in dynamical equilibrium are born close to the plane within the scale height of the cold gas. The remaining 20%–30% of stars are injected under nonequilibrium conditions. The LHS figures are similar in their temporal behavior for all $f_{\text{gas}}$, unlike the RHS figures. In the RHS panel, the disk heating rate for formed stars in time increases rapidly with $f_{\text{gas}}$, all starting from a low base set by the maximum dispersion allowed to turbulent-pressure supported molecular gas ($\le 20–30 \text{ km s}^{-1}$; S. M. R. Jeffreson et al. 2022). This is because the Brownian motion arising from baryon sloshing heats the stellar disk essentially isotropically, and the amplitude of this motion is a strong function of $f_{\text{gas}}$. Eventually, as the younger stellar component settles, its heating rate emerges from the noise and matches the preexisting population. This is evident by matching the different curves one for one in Figure 13, starting with the oldest tracks and moving backwards in time.

In summary, although star formation and turbulence are involved, our mechanism is very different. The energy to thicken the stellar disk is originally from star formation and then converted to bulk kinetic motion of the baryonic gravitational potential with respect to an inertial frame (approximated by the dark matter frame), which is subsequently diffused into the thermal kinetic energy of the stars.

## 6. Looking Forward

The formation and evolution of disk galaxies is recognized as a difficult problem that is far from resolution. Once again, observations lead theory with the remarkable revelation (Section 1) of the early onset of disks from JWST and ALMA observations, in particular. However, more observational progress is needed: there is a distinct lack of galaxies with both high-quality JWST *and* ALMA data.

Here we do not consider the early onset problem. We focus on the simple fact that most disks must have gone through a gas-rich phase given that half of all Milky Way stars were born before $z \sim 1$ when the Universe was less than a third of its current age (K. Freeman & J. Bland-Hawthorn 2002; K. El-Badry et al. 2018). Very little is known about the evolution of gas-rich disks, and *Nexus* was developed with the specific goal of shedding light on this important epoch. What little we do know has come from cosmological simulations at much lower resolution. Here we demonstrate the benefits of increased resolution by revealing new physical processes in these early disks.

In our first *Nexus* paper (J. Bland-Hawthorn et al. 2024), we showed how bulges, bars, and spiral arms emerge even in turbulent, gas-rich environments, consistent with the

remarkable new JWST results (Y. Guo et al. 2023; Z. A. Le Conte et al. 2024). In the present work, we focus on dynamical processes that unbind the disk and cause it to slosh around. The sloshing leads to a thickened stellar disk with properties remarkably consistent with our Galactic $\alpha$-enriched thick disk, in particular, the prediction of vertical gradients in age and metallicity, with only a weak to no radial gradient. These are complex nonlinear processes that can only be addressed through numerical simulations presently. Analytic models that describe the evolution of fast and slow fluctuations are under way but at present they have been applied to only a few restricted cases (C. Hamilton et al. 2024; J. Binney 2025).

There is little doubt that gas-dominated disks were abundant at early times. This has been evident for the past two decades, starting with early molecular gas observations and ground-based integral field studies of ionized gas (S. C. Chapman et al. 2004; N. M. Förster Schreiber et al. 2006; R. Genzel et al. 2006; K. L. Shapiro et al. 2008; A. M. Swinbank et al. 2012). The most obvious manifestations of turbulent galaxies at high redshift are the observed high gas fractions and large gas clouds (e.g., J. Fensch & F. Bournaud 2021), elevated star formation rates and broadened gas kinematics seen in ALMA and JWST surveys (Section 1). The enhanced star-forming disks at high redshift exhibit larger velocity dispersions than for local disks, regardless of the emission diagnostic used (T. Ejdetjärn et al. 2022), and presumably these reflect enhanced levels of turbulent energy.

A strong prediction of our work is that early disks undergo “baryon sloshing” with a spatial and kinematic amplitude that gets larger with increasing $f_{\text{gas}}$. Can we find direct evidence of this process? What are the observable manifestations of the predicted sloshing action? There are already tantalizing clues. B. P. Venemans et al. (2020) has shown that quasars display $\sim 20 \text{ km s}^{-1}$ velocity offsets and $\sim 1 \text{ kpc}$ spatial offsets in projection with respect to their accompanying gas-rich disks. These amplitudes are consistent with our high $f_{\text{gas}}$ models. From JWST imaging, A. Le Bail et al. (2024) find that galaxy disk asymmetries dominate beyond $z \sim 1$ becoming less prevalent at lower redshift; in addition, this behavior is seen in the absence of companions and in all environments. These recent results may have other explanations but they are entirely consistent with our expectations.

In a companion paper (J. Bland-Hawthorn et al. 2025, in preparation), we explore both morphological and kinematic signatures of baryon sloshing. We highlight specific signatures in low-order ($m = 1, 2, 3$) Fourier components and show that these effects may already be evident in high-redshift ALMA observations (e.g., T. Tsukui & S. Iguchi 2021; T. Tsukui et al. 2024). We focus our study on both the morphology and kinematics of the starlight and the cool gas for high-redshift galaxies up to $z \sim 4$. There are now dozens of galaxies that have been mapped by ALMA in cool gas out to this redshift, as discussed in Section 1. The notion of stellar kinematics beyond $z \sim 1$ may seem fanciful, except that at least two sources now have (sufficiently) resolved stellar velocity fields in the range $z = 3–4$ (F. D’Eugenio et al. 2024; P. G. Pérez-González et al. 2025), with more to come. This was made possible by the remarkable sensitivity of the JWST NIRSpec integral field spectrograph. These are direct line-of-sight observations unaided by gravitational lensing. Reconstructed velocity fields from lensing are also now possible in the same redshift window. Consistent with our models, these stellar disks are
<figure><img src="image_20.png" /><figcaption>**Figure 15.** (Left) Where new stars are born in time as a function of vertical height $|z|$. The left-hand side (LHS) panels look very similar for all $f_{\text{gas}}$, quite unlike the RHS panels. (Right) The rms z-scale height for newborn stars as a function of radius $R$ (in kiloparsecs) for each value of $f_{\text{gas}}$ for the halo accretion cases. In all plots, the color coding indicates the time evolution, as indicated. The bottom 10 tracks (100 Myr apart) are removed because of large intrinsic scatter; this is a consequence of small number statistics and early populations being born out of dynamical equilibrium. The scale height of the $\alpha$-rich thick disk today at $R = 8$ kpc are indicated by the blue star (J. Bland-Hawthorn & O. Gerhard 2016). Note that the vertical scale has a fourfold stretch compared to the horizontal axis. The gas-rich models are able to reproduce the thickness of the $\alpha$-rich disk within 1.5–2 Gyr.</figcaption></figure>
relatively young (300–700 Myr) and thus detectable in Balmer absorption lines. For both disks, the star formation rates are of order $10 M_{\odot} yr^{-1}$ and considered “quiescent” for such an early epoch, but these rates were an order of magnitude higher at even earlier times.

In future work, already under way, we will look at how merger activity affects our results during the gas-rich phase. It will be important to separate and quantify the respective roles of baryon sloshing and merger activity. We will evolve the Nexus models over the lifetime of the simulated Milky Way–like galaxies. What are the likely chemodynamic signatures that survive? We conjecture that the different baryon components will have different heating signatures, induced by the sloshing, compared to other scenarios. While such orbits can arise from past mergers (e.g., C. Gallart et al. 2019; J. A. S. Amarante et al. 2020), there may be chemodynamic signatures that separate these two distinct processes. If most galaxies are ultimately found to have chemically bimodal disks, given that most disks go through a gas-rich phase in the early Universe, the sloshing mechanism is arguably the most natural process to explain their prevalence.

We are upgrading Nexus with more complex gas processes, including MHD and cosmic-ray heating, AGN heating, and molecule and dust formation primarily. This will allow our models to be compared to a wider class of objects at different epochs, with varying properties (e.g., total mass, star formation rate) observed across the full electromagnetic spectrum. These are important steps in understanding the full narrative that takes us from early disk systems to modern-day Milky Way counterparts.

## Acknowledgments

J.B.H. acknowledges a Université Paris Science et Lettres (PSL) professorial fellowship that partially supported the writing stage of this work. He is indebted to the graduate students at Paris and Meudon, for their engagement in the graduate course and for inspiring conversations on early disk formation.

J.B.H. is indebted to Merton College and the Beecroft building for providing a stimulating research environment. We are particularly indebted to James Binney for an advanced copy of his new book Stellar Dynamics (2025), and for its innovative analysis of dynamical fluctuations that was pivotal for this work. We also thank him for providing the code for the calculation in Section 4.2. Moreover, we received valuable feedback from seminars delivered at Marseille, Paris,

Shanghai, UCL, Cambridge, and Oxford. Important questions were raised by Lia Athanassoula, Chris Hamilton, Francoise Combes, Bruce Elmegreen, Walter Dehnen, Wyn Evans, and Vasily Belokurov, and we have sought to answer them here. We are grateful to Eugene Vasiliev and Romain Teyssier for continued assistance and adaptions with AGAMA and RAMSES, respectively. We are also indebted to a thorough and perceptive referee.

T.T.G. acknowledges financial support from the Australian Research Council (ARC) through Australian Laureate Fellowships awarded to JBH (FL140100278) and TRB (FL220100117), and for partial funding through the James Arthur Pollock memorial fund awarded to the School of Physics, University of Sydney. O.A. acknowledges support from the Knut and Alice Wallenberg Foundation, the Swedish Research Council (grant 2019-04659) and the Swedish National Space Agency (SNSA Dnr 2023-00164). C.F. acknowledges funding provided by the Australian Research Council (Discovery Project DP230102280 and DP250101526), and the Australia-Germany Joint Research Cooperation Scheme (UA-DAAD).

The computations and data storage were enabled by two facilities: (i) the National Computing Infrastructure (NCI) Adapter Scheme, provided by NCI Australia, an NCRIS capability supported by the Australian Government; and (ii) LUNARC, the Centre for Scientific and Technical Computing at Lund University (resource allocations LU 2023/2-39 and LU 2023/12-6). We further acknowledge high-performance computing resources provided by the Leibniz Rechenzentrum and the Gauss Centre for Supercomputing (grants pr32lo, pr48pi, pn76ga and GCS Large-scale project 10391), the Australian National Computational Infrastructure (grant ek9) and the Pawsey Supercomputing Centre (project pawsey0810) in the framework of the National Computational Merit Allocation Scheme and the ANU Merit Allocation Scheme.

# Appendix

## Radial and Kinematic Signatures of Baryon Sloshing

In Section 3.3, we presented initial evidence of the 3D baryon sloshing, with kinematic and radial information compared in Figure 7 (Right) for a few specific cases. The full gallery of radial and kinematic signatures is presented in Figure 16. As expected, the radial and kinematic excursion amplitudes increase with $f_{gas}$. Note also the decreasing frequency of oscillation as the baryon centroid wanders further and further from the center of total gravitational potential. This
<figure><img src="image_21.png" /></figure>

**Figure 16.** Baryon sloshing presented as the change in radius $R_{\text{slosh}} = \sqrt{x_{\text{slosh}}^2 + y_{\text{slosh}}^2 + z_{\text{slosh}}^2}$ with cosmic time (upper panels) in the reference frame of the total potential for $f_{\text{gas}} = 0\%, 20\%, 40\%, 60\%, 80\%, 100\%$. The lower panels are the derivative $\dot{R}_{\text{slosh}}$ of the upper panels. The red lines indicate the changing physical scales between panels, where the overall trend is increasing with $R$, as expected from Equation (15).

is a natural consequence of the radial trends in Equation (15) (see also Equation (11)).

## ORCID iDs

Joss Bland-Hawthorn https://orcid.org/0000-0001-7516-4016

Thor Tepper-Garcia https://orcid.org/0000-0002-1081-883X

Oscar Agertz https://orcid.org/0000-0002-4287-1088

Christoph Federrath https://orcid.org/0000-0002-0706-2306

Misha Haywood https://orcid.org/0000-0003-0434-0400

Paola di Matteo https://orcid.org/0000-0002-5213-4807

Timothy R. Bedding https://orcid.org/0000-0001-5222-4661

Takafumi Tsukui https://orcid.org/0000-0002-1499-6377

Emily Wisnioski https://orcid.org/0000-0003-1657-7878

Melissa Ness https://orcid.org/0000-0001-5082-6693

Ken Freeman https://orcid.org/0000-0001-6280-1207

## References

Abadi, M. G., Navarro, J. F., Steinmetz, M., & Eke, V. R. 2003, ApJ, 597, 21

Aerts, C. 2021, RvMP, 93, 015001

Agertz, O., Kravtsov, A. V., Leitner, S. N., & Gnedin, N. Y. 2013, ApJ, 770, 25

Agertz, O., Renaud, F., Feltzing, S., et al. 2021, MNRAS, 503, 5826

Agertz, O., Teyssier, R., & Moore, B. 2009, MNRAS, 397, L64

Amarante, J. A. S., Beraldo e Silva, L., Debattista, V. P., & Smith, M. C. 2020, ApJL, 891, L30

Amvrosiadis, A., Lange, S., Nightingale, J. W., et al. 2025, MNRAS, 537, 1163

Anders, F., Chiappini, C., Santiago, B. X., et al. 2014, A&A, 564, A115

Anders, F., Gispert, P., Ratcliffe, B., et al. 2023, A&A, 678, A158

Antoja, T., Helmi, A., Romero-Gómez, M., et al. 2018, Natur, 561, 360

Assmann, P., Fellhauer, M., Kroupa, P., Brüns, R. C., & Smith, R. 2011, MNRAS, 415, 1280

Aumer, M., Binney, J., & Schönrich, R. 2016, MNRAS, 462, 1697

Binney, J. 2020, MNRAS, 496, 767

Binney, J. 2025, Stellar Dynamics: First Edition (Princeton, NJ: Princeton Univ. Press)

Bird, J. C., Loebman, S. R., Weinberg, D. H., et al. 2021, MNRAS, 503, 1815

Birnboim, Y., & Dekel, A. 2003, MNRAS, 345, 349

Bland-Hawthorn, J., & Gerhard, O. 2016, ARA&A, 54, 529

Bland-Hawthorn, J., Krumholz, M. R., & Freeman, K. 2010, ApJ, 713, 166

Bland-Hawthorn, J., & Sharma, S. 2016, AN, 337, 894

Bland-Hawthorn, J., & Tepper-García, T. 2021, MNRAS, 504, 3168

Bland-Hawthorn, J., Tepper-Garcia, T., Agertz, O., & Federrath, C. 2024, ApJ, 968, 86

Bland-Hawthorn, J., Tepper-Garcia, T., Agertz, O., & Freeman, K. 2023, ApJ, 947, 80

Bournaud, F., Elmegreen, B. G., & Martig, M. 2009, ApJL, 707, L1

Bovy, J., Rix, H.-W., & Hogg, D. W. 2012a, ApJ, 751, 131

Bovy, J., Rix, H.-W., Liu, C., et al. 2012b, ApJ, 753, 148

Casamiquela, L., Reese, D. R., Lebreton, Y., et al. 2024, A&A, 692, A243

Cerqui, V., Haywood, M., Snaith, O., Di Matteo, P., & Casamiquela, L. 2025, A&A, 699, A277

Chan, T. K., Kereš, D., Oñorbe, J., et al. 2015, MNRAS, 454, 2981

Chapman, S. C., Smail, I., Windhorst, R., Muxlow, T., & Ivison, R. J. 2004, ApJ, 611, 732

Chen, H.-W. 2012, MNRAS, 427, 1238

Costantin, L., Pérez-González, P. G., Guo, Y., et al. 2023, Natur, 623, 499

Dekel, A., & Birnboim, Y. 2006, MNRAS, 368, 2

Dekel, A., Sari, R., & Ceverino, D. 2009, ApJ, 703, 785

D'Eugenio, F., Pérez-González, P. G., Maiolino, R., et al. 2024, NatAs, 8, 1443

Egorov, O. V., Kreckel, K., Glover, S. C. O., et al. 2023, A&A, 678, A153

Ejdetjärn, T., Agertz, O., Östlin, G., Renaud, F., & Romeo, A. B. 2022, MNRAS, 514, 480

El-Badry, K., Bland-Hawthorn, J., Wetzel, A., et al. 2018, MNRAS, 480, 652

Elmegreen, B. G., & Elmegreen, D. M. 2005, ApJ, 627, 632
Fensch, J., & Bournaud, F. 2021, MNRAS, 505, 3579

Fensch, J., Bournaud, F., Brucy, N., et al. 2023, A&A, 672, A193

Ferreira, L., Adams, N., Conselice, C. J., et al. 2022, ApJL, 938, L2

Ferreira, L., Conselice, C. J., Sazonova, E., et al. 2023, ApJ, 955, 94

Förster Schreiber, N. M., Genzel, R., Lehnert, M. D., et al. 2006, ApJ, 645, 1062

Freeman, K., & Bland-Hawthorn, J. 2002, ARA&A, 40, 487

Gallart, C., Bernard, E. J., Brook, C. B., et al. 2019, NatAs, 3, 932

Genzel, R., Burkert, A., Bouché, N., et al. 2008, ApJ, 687, 59

Genzel, R., Tacconi, L. J., Eisenhauer, F., et al. 2006, Natur, 442, 786

Genzel, R., Tacconi, L. J., Gracia-Carpio, J., et al. 2010, MNRAS, 407, 2091

Green, M. S. 1954, JChPh, 22, 398

Guo, Y., Jogee, S., Finkelstein, S. L., et al. 2023, ApJL, 945, L10

Hamilton, C., Modak, S., & Tremaine, S. 2024, arXiv:2408.03366

Hamilton, C., Modak, S., & Tremaine, S. 2025, arXiv:2507.16950

Hayden, M. R., Bovy, J., Holtzman, J. A., et al. 2015, ApJ, 808, 132

Hayden, M. R., Holtzman, J. A., Bovy, J., et al. 2014, AJ, 147, 116

Haywood, M., Di Matteo, P., Lehnert, M. D., Katz, D., & Gómez, A. 2013, A&A, 560, A109

Haywood, M., Di Matteo, P., Snaith, O., & Lehnert, M. D. 2015, A&A, 579, A5

Hills, J. G. 1980, ApJ, 235, 986

Hodge, J. A., da Cunha, E., Kendrew, S., et al. 2025, ApJ, 978, 165

Hodge, J. A., Smail, I., Walter, F., et al. 2019, ApJ, 876, 130

Houdek, G., & Dupret, M.-A. 2015, LRSP, 12, 8

Ida, S., Kokubo, E., & Makino, J. 1993, MNRAS, 263, 875

Jeffreson, S. M. R., Sun, J., & Wilson, C. D. 2022, MNRAS, 515, 1663

Jiang, I.-G., & Binney, J. 2000, MNRAS, 314, 468

Jiménez, E., Lagos, C. d. P., Ludlow, A. D., & Wisnioski, E. 2023, MNRAS, 524, 4346

Joshi, R., & Widrow, L. M. 2024, MNRAS, 527, 7781

Julian, W. H., & Toomre, A. 1966, ApJ, 146, 810

Kartaltepe, J. S., Rose, C., Vanderhoof, B. N., et al. 2023, ApJL, 946, L15

Kereš, D., Katz, N., Fardal, M., Davé, R., & Weinberg, D. H. 2009, MNRAS, 395, 160

Kereš, D., Katz, N., Weinberg, D. H., & Davé, R. 2005, MNRAS, 363, 2

Kohandel, M., Pallottini, A., Ferrara, A., et al. 2024, A&A, 685, A72

Kroupa, P. 2002, MNRAS, 330, 707

Kubo, R., Yokota, M., & Nakajima, S. 1957, JPSJ, 12, 1203

Kuhlen, M., Guedes, J., Pillepich, A., Madau, P., & Mayer, L. 2013, ApJ, 765, 10

Kuhn, V., Guo, Y., Martin, A., et al. 2024, ApJL, 968, L15

Lau, J. Y., & Binney, J. 2021, MNRAS, 507, 2241

Le Bail, A., Daddi, E., Elbaz, D., et al. 2024, A&A, 688, A53

Le Conte, Z. A., Gadotti, D. A., Ferreira, L., et al. 2024, MNRAS, 530, 1984

Loebman, S. R., Roškar, R., Debattista, V. P., et al. 2011, ApJ, 737, 8

Ma, X., Hopkins, P. F., Wetzel, A. R., et al. 2017, MNRAS, 467, 2430

Mackereth, J. T., Bovy, J., Schiavon, R. P., et al. 2017, MNRAS, 471, 3057

Mackereth, J. T., Schiavon, R. P., Pfeffer, J., et al. 2019, MNRAS, 482, 3426

Meng, X., & Gnedin, O. Y. 2021, MNRAS, 502, 1433

Miglio, A., Chiappini, C., Mackereth, J. T., et al. 2021, A&A, 645, A85

Mo, H. J., Mao, S., & White, S. D. M. 1998, MNRAS, 295, 319

Moster, B. P., Macciò, A. V., Somerville, R. S., Johansson, P. H., & Naab, T. 2010, MNRAS, 403, 1009

Moster, B. P., Macciò, A. V., Somerville, R. S., Naab, T., & Cox, T. J. 2011, MNRAS, 415, 3750

Navarro, J. F., Frenk, C. S., & White, S. D. M. 1996, ApJ, 462, 563

Nelson, E., Brammer, G., Giménez-Arteaga, C., et al. 2024, ApJL, 976, L27

Nelson, R. W., & Tremaine, S. 1999, MNRAS, 306

Nipoti, C., & Binney, J. 2015, MNRAS, 446, 1820

Ogiya, G., & Mori, M. 2011, ApJL, 736, L2

Pallottini, A., Gallerani, S., & Ferrara, A. 2014, MNRAS, 444, L105

Pérez-González, P. G., D'Eugenio, F., Rodríguez del Pino, B., et al. 2025, NatAs, 9, 1240

Pillepich, A., Nelson, D., Springel, V., et al. 2019, MNRAS, 490, 3196

Pontzen, A., & Governato, F. 2012, MNRAS, 421, 3464

Pope, A., McKinney, J., Kamieneski, P., et al. 2023, ApJL, 951, L46

Qu, Y., Di Matteo, P., Lehnert, M. D., & van Driel, W. 2011, A&A, 530, A10

Quinn, P. J., Hernquist, L., & Fullagar, D. P. 1993, ApJ, 403, 74

Ratcliffe, B., Khoperskov, S., Minchev, I., et al. 2025, A&A, 698, A267

Rizzo, F., Kohandel, M., Pallottini, A., et al. 2022, A&A, 667, A5

Robin, A. C., Haywood, M., Creze, M., Ojha, D. K., & Bienayme, O. 1996, A&A, 305, 125

Roman-Oliveira, F., Rizzo, F., & Fraternali, F. 2024, A&A, 687, A35

Rowland, L. E., Hodge, J., Bouwens, R., et al. 2024, MNRAS, 535, 2068

Sahlholdt, C. L., Feltzing, S., & Feuillet, D. K. 2022, MNRAS, 510, 4669

Santistevan, I. B., Wetzel, A., El-Badry, K., et al. 2020, MNRAS, 497, 747

Shapiro, K. L., Genzel, R., Förster Schreiber, N. M., et al. 2008, ApJ, 682, 231

Sharma, S., Hayden, M. R., & Bland-Hawthorn, J. 2021, MNRAS, 507, 5882

Sharma, S., Stello, D., Bland-Hawthorn, J., et al. 2019, MNRAS, 490, 5335

Sike, B., Thomas, T., Ruszkowski, M., Pfrommer, C., & Weber, M. 2025, ApJ, 987, 204

Snaith, O., Haywood, M., Di Matteo, P., et al. 2015, A&A, 578, A87

Snaith, O. N., Haywood, M., Di Matteo, P., et al. 2014, ApJL, 781, L31

Stern, J., Faucher-Giguère, C.-A., Fielding, D., et al. 2021, ApJ, 911, 88

Stern, J., Fielding, D., Faucher-Giguère, C.-A., & Quataert, E. 2020, MNRAS, 492, 6042

Stinson, G. S., Bovy, J., Rix, H. W., et al. 2013, MNRAS, 436, 625

Swinbank, A. M., Smail, I., Sobral, D., et al. 2012, ApJ, 760, 130

Tacconi, L. J., Genzel, R., & Sternberg, A. 2020, ARA&A, 58, 157

Tamfal, T., Mayer, L., Quinn, T. R., et al. 2022, ApJ, 928, 106

Tepper-García, T., Bland-Hawthorn, J., Vasiliev, E., et al. 2024, MNRAS, 535, 187

Teyssier, R. 2002, A&A, 385, 337

Tollet, E., Macciò, A. V., Dutton, A. A., et al. 2016, MNRAS, 456, 3542

Toomre, A. 1964, ApJ, 139, 1217

Toomre, A. 1981, in Structure and Evolution of Normal Galaxies, ed. S. M. Fall & D. Lynden-Bell (Cambridge: Cambridge Univ. Press), 111

Tsukui, T., & Iguchi, S. 2021, Sci, 372, 1201

Tsukui, T., Wisnioski, E., Bland-Hawthorn, J., & Freeman, K. 2025, MNRAS, 540, 3493

Tsukui, T., Wisnioski, E., Bland-Hawthorn, J., et al. 2024, MNRAS, 527, 8941

van Dokkum, P. G., Leja, J., Nelson, E. J., et al. 2013, ApJL, 771, L35

van Donkelaar, F., Agertz, O., & Renaud, F. 2022, MNRAS, 512, 3806

Vasiliev, E. 2019, MNRAS, 482, 1525

Venemans, B. P., Walter, F., Neeleman, M., et al. 2020, ApJ, 904, 130

Villalobos, Á., & Helmi, A. 2008, MNRAS, 391, 1806

Walker, I. R., Mihos, J. C., & Hernquist, L. 1996, ApJ, 460, 121

Weinberg, M. D. 2001, MNRAS, 328, 311

Widrow, L. M., Ryu, D., Schleicher, D. R. G., et al. 2012, SSRv, 166, 37

Wright, R. J., Somerville, R. S., Lagos, C. d. P., et al. 2024, MNRAS, 532, 3417

Xiang, M., & Rix, H.-W. 2022, Natur, 603, 599

Yoachim, P., & Dalcanton, J. J. 2006, AJ, 131, 226

Yu, S., Bullock, J. S., Klein, C., et al. 2021, MNRAS, 505, 889

Zhao, H. 2002, MNRAS, 336, 159