© The Authors 2025

1. Introduction

There has been mounting observational evidence in the past few years that star-forming disc galaxies harbour a complex, filamentary network of dense gas and dust (Thilker et al. 2023, and references therein). These are elongated, kiloparsec-scale, dense features known to host star-forming regions, H II regions and young star clusters (Schinnerer et al. 2017; Williams et al. 2022). Not only are they found in external galaxies, but there has also been evidence of such large-scale filamentary structures in the Milky Way (Zucker et al. 2015; Pantaleoni González et al. 2021; Kuhn et al. 2021; Veena et al. 2021). External galaxies, however, give us the advantage of the ‘birds eye view’, which allows us to study filament structure and morphology. Because of their richness in dust, they are visible as attenuation features in Hubble Space Telescope (HST) images (Elmegreen 1980; La Vigne et al. 2006), or as emission features in the near/mid-infrared, for example in the 8 μm infrared band of Spitzer (Elmegreen et al. 2006, 2014, 2018; Elmegreen & Elmegreen 2019) and in recent JWST/Mid-Infrared Instrument (MIRI) observations (Williams et al. 2022; Thilker et al. 2023; Meidt et al. 2023), which have reached unprecedented resolutions of ∼10 pc. This enables us for the first time to compare the detailed morphology of dense gas in simulations with observations of nearby galaxies and to test the theories for the origins of filamentary structures.

The question of their origin is far from settled, and a number of different mechanisms have been proposed. Filaments have been most commonly found in grand-design spiral galaxies, and within spirals most frequently in Sb-Sc galaxies (La Vigne et al. 2006). This increased detection frequency in grand-design spirals has led some authors to posit that filament formation is driven by spiral arms, and to carry out linear stability analyses and two-dimensional (2D) simulations showing that gas traversing arms created by an external spiral potential is prone to a hydrodynamical instability that causes them to develop over-densities, which are then sheared into filamentary structures in the inter-arm regions (Wada & Koda 2004; Dobbs & Bonnell 2006; Lee & Shu 2012; Kim et al. 2014; Lee 2014; Kim et al. 2015; Sormani et al. 2017; Mandowara et al. 2022). Other authors have attributed filaments to Parker instability driven by large-scale magnetic fields (Körtgen et al. 2019) or supernova feedback from star-forming regions in the arms (Kim et al. 2020). Yet other authors suggest that gas self-gravity alone could potentially create filamentary structures without the need to invoke external potentials, supernova feedback or magnetic fields (Griv & Gedalin 2012; Meidt 2022).

At present, however, we lack a systematic, 3D, global, numerical study that aims to disentangle the contribution of different physical processes to dense structure formation in galaxies. Prior works that explicitly focus on gravitational instability in a galactic disc are either 2D (Shetty & Ostriker 2006; Griv & Wang 2014) or are low-resolution models with limited parameter space that were unable to resolve the spatial scales necessary for filament formation (Kim & Ostriker 2006; Wang et al. 2010). More recent and higher resolution simulations have different aims. Some are focused on studying Milky-Way analogues that have similar mid-plane pressures, distribution of the total mass into the halo, stellar bulge and the gas components as our own Galaxy (Kim et al. 2016; Jeffreson et al. 2020). Others include complex physical processes such as an external spiral potential and subgrid models for star formation and feedback that make it hard to separate effects of different physical processes (for example; Tress et al. 2020; Robinson & Wadsley 2024; Zhao et al. 2024; Khullar et al. 2024).

The aim for this study is to use numerical simulations to determine the conditions of filament formation in an idealised setting where we have control over the different potential drivers of filament formation and thereby can carry out clean numerical experiments to gain insight into the underlying physics. To this end, we carry out high-resolution 3D simulations of isothermal, self-gravitating isolated disc galaxies that are initialised in equilibrium. The main question we address is whether self-gravity by itself, in the absence of the other mechanisms described above (such as magnetic fields, stellar spiral potentials, supernovae), is sufficient to create the dense filamentary network observed in galaxies. We demonstrate that it is, and conduct a parameter study to quantify the morphological features of the filamentary structures that form in our simulations. Finally, we compare the physical spacing of the simulated filaments with recent observations.

The paper is organized as follows: Section 2 describes our simulation setup, initial conditions and the library of simulations used in this work. In Section 3, we discuss the basic morphology of our galaxies and the filaments that form in them. We further quantify the dependence of the morphology and filament formation timescales on the initial galaxy-wide parameters. In Section 4, we compare our results with observations and existing theoretical and numerical work. Finally, we summarise our conclusions in Section 5.

2. Methods

2.1. Physical setup

Our 3D disc galaxy simulations consist of self-gravitating isothermal gas in an analytical axisymmetric, time-steady dark matter plus stellar potential with a flat rotation curve. We initialise our disc galaxies in equilibrium using the “potential method” described by Wang et al. (2010). Here, we briefly describe the setup.

2.1.1. Basic setup

We seek to initialise our simulated galaxies in equilibrium, such that the forces due to thermal pressure, self-gravity, dark matter and rotation balance in both the radial and vertical directions. The surface density of the galaxies follows a modified exponential profile given by

$$\begin{array}{c}\hfill \mathrm{\Sigma }(R)={\mathrm{\Sigma }}_{°}exp[-\frac{R}{R_{\mathrm{d}}}-\beta exp(-\alpha \frac{R}{R_{\mathrm{d}}})],\end{array}$$(1)

where R is the galactic (cylindrical) radius, Σ_° is a scaling constant that will vary from run to run, R_d = 3 kpc is the disc scale length, and α = 2 and β = 1/2 are constants that determine the shape of the surface density profile near the centre of the galaxy. We use this function instead of the pure exponential in order to flatten the surface density profile close to the centre and guarantee that the scale height there is resolved well enough (≥4 cells) to avoid artificial numerical oscillations due to inadequate resolution of the vertical pressure gradient that balances self-gravity. Since we are not interested in the central part of the galaxy, this does not affect our results.

The rotation profile of the disc is given by

$$\begin{array}{c}\hfill v_{\mathrm{rot}}=v_{\mathrm{c}}\frac{R}{\sqrt{R^2+R_{\mathrm{c}}^2}},\end{array}$$(2)

where v_c is the saturated circular velocity of the gas and R_c = 2 kpc. The gas follows an isothermal equation of state with

$$\begin{array}{c}\hfill P=\rho {c_{\mathrm{s}}}^2,\end{array}$$(3)

where P is the gas pressure, ρ is the density, and c_s is the sound speed.

We use the constraint of radial equilibrium to determine the stellar plus dark matter potential ϕ_dm. For the exact expression and details regarding the procedure, the reader is referred to Appendix A. We then use this to place the disc in vertical equilibrium, where we balance the force due to the gas self-gravity and the gravity due to the stellar and dark matter components with the vertical pressure gradient of the gaseous disc at each galactocentric radius,

$$\begin{array}{c}\hfill \rho \frac{\partial }{\partial z}({\varphi }_{\mathrm{g}}+{\varphi }_{\mathrm{dm}})+{c_{\mathrm{s}}}^2\frac{\partial \rho }{\partial z}=0,\end{array}$$(4)

where ϕ_g is the gravitational potential of the self-gravitating gas (see Appendix A). This is subject to the constraint that

$$\begin{array}{c}\hfill 2{\displaystyle {\int }_0^{\infty }}\rho (R,z)\mathrm{d}z=\mathrm{\Sigma }(R),\end{array}$$(5)

where Σ(R) is the surface density profile (Equation (1)). Since ϕ_g and ϕ_dm do not depend on ρ in our simple approximation, we can find a vertical profile ρ(R, z) that satisfies these two constraints numerically. We do this by starting with an initial guess of the midplane density ρ(R, 0) and integrating Equation (4) to obtain the profile ρ(R, z), which is then rescaled by a constant factor in order to satisfy Equation (5). We point out that this is an approximate procedure rather than exact, since the gas gravitational potential ϕ_g should depend on ρ. However, we have found that the approximation is sufficiently accurate that our initial discs quickly relax to an exact equilibrium.

Finally, we immerse our galactic discs in a uniform hot circum-galatic medium (CGM) with a constant T_CGM = 10⁷ K and density. We transition from the cold galactic disc to the hot CGM when our discs reach a fixed transition density of ρ = 10⁻²⁸ cm⁻³, which is ∼5 orders of magnitude less than the central density of our models. The density of the CGM is such that there is no pressure gradient between the two media.

2.1.2. Turbulent velocity field

The rotation curve (Equation (2)) gives the mean initial velocity of the gas, but we also impose a perturbation on top of this to seed the instabilities that we strive to study. Our procedure for doing so mirrors that described in Arora et al. (2024): we impose a turbulent initial velocity field with a sonic Mach number ℳ = 0.5 and a Kolmogorov scaling of k^−5/3 on scales [50,  200] pc. These are close to the driving scales of turbulence of ∼140 ± 80 pc estimated in the Milky Way (Chepurnov et al. 2010). The turbulent velocity field contains a natural mixture of solenoidal and compressible modes, and we generate it with the methods described in Federrath et al. (2010), using the publicly available TurbGen code (Federrath et al. 2022).

2.2. Parameter study

We characterise our disc galaxies with two dimensionless parameters,

$$\begin{array}{cc}\hfill ⟨Q⟩& =\frac{1}{R_{\max }}{\displaystyle {\int }_0^{R_{\max }}}Q(R)\mathrm{d}R=\frac{1}{R_{\max }}{\displaystyle {\int }_0^{R_{\max }}}\frac{\kappa c_{\mathrm{s}}}{\pi G\mathrm{\Sigma }}\mathrm{d}R,\hfill \end{array}$$(6)

$$\begin{array}{cc}\hfill {\mathcal{M}}_{\mathrm{c}}& =\frac{v_c}{c_{\mathrm{s}}},\hfill \end{array}$$(7)

where κ is the epicyclic frequency given by κ² = 4Ω² + 2RΩ(dΩ/dR), Ω = v_rot/r is the angular frequency, and we take R_max = 8 kpc. The first parameter above is the mean value of the Toomre-Q parameter Q(R) (Toomre 1964) over the inner 8 kpc of the galaxy; it quantifies the gravitational stability of our galaxies. The second parameter describes the ratio of rotation speed to sound speed in the region R ≫ R_c where the rotation curve reaches its full speed; it relates to the rotational and thermodynamic properties of the galaxy and controls its thickness¹. We show in Appendix B that these two dimensionless parameters, along with the auxiliary dimensionless ratios R_max/R_d, R_c/R_d, α, β, and q that we do not vary, fully determine both the initial state and the subsequent evolution of the system. In particular, we note that disc thickness is not an independent parameter, but is instead fully determined by ⟨Q⟩ and ℳ_c – primarily the latter.

The parameter space of our simulations can be split into two main groups. The first group varies Toomre-Q with ⟨Q⟩≃{1, 2, 3}. (The precise values are not exactly 1, 2, and 3, and are given in Table 1.) We show the radial profiles of Q(R) for these three cases in Figure 1. The solid lines show Q at the specified galactocentric radius and the dot-dashed lines indicate the average ⟨Q⟩ in the region 0 ≤ R/kpc ≤ 8. The second group has constant ⟨Q⟩ = 1, but varies ℳ_c with ℳ_c ≃ {14, 21, 29, 36}; in terms of dimensional parameters, we do this by varying v_c while keeping c_s = 7 km s⁻¹ fixed. The range of v_c is chosen to reflect the range of rotation speeds observed in galaxies, with v_c ∈ {100,150,200,250} km s⁻¹ (Levy et al. 2018; Lang et al. 2020). Our full parameter space is summarised in Table 1. For the sake of brevity, we refer to the different runs using the rounded-off values of the dimensionless parameters in the rest of the work.

Fig. 1. Initial radial Toomre-Q (see Equation (6)) profiles of our galaxies. The solid lines show Q as a function of galactocentric radius and the dot-dashed lines show ⟨Q⟩, the average of the solid curves over the disc region 0 ≤ R/kpc ≤ 8.

2.3. Numerical methods

We use the FLASH hydrodynamical code for performing the simulations (Fryxell et al. 2000; Dubey et al. 2008). The disc is initialised at the centre of a cuboidal box with side length L_x = L_y = 20 kpc in the plane of the disc and L_z = 2.5 kpc in the direction perpendicular to it. We use outflow boundary conditions as in earlier works (Körtgen et al. 2019; Arora et al. 2024), but the region we use for our analysis is limited to 2 ≤ R (kpc)≤4 and H ≤ 250 pc for our thickest galaxy, far enough from the edges of our box that the boundaries have negligible effects.

We use adaptive mesh refinement (AMR) for our simulations. Our base grid is 256 × 256 × 32 cells, giving a cell size ≈80 pc. However, we achieve a maximum effective resolution of 4096 × 4096 × 512, cells corresponding to a minimum cell size of ≈4.8 pc, by adding 4 additional levels of refinement. We apply this extra refinement to the same geometric region in all the simulations: a cylinder whose centre coincides with the centre of the galaxy and with radius R_cyl = 6 kpc and half-height |H_cyl| = 100 pc centred on the galactic plane. This ensures that all the simulations have identical spatial resolutions. For our thinnest disc, at this resolution we have ≥4 cells per galactic scale height at galactic centre, and ≥8 cells per scale height for our region of interest. Outside the cylindrical region where we enforce the highest refinement level, we refine and coarsen based on the condition that the Jeans length be resolved by at least 32 and not more than 64 grid cells (Federrath et al. 2011).

3. Results

3.1. Disc morphology

We begin our analysis by showing the time evolution of one example run, Q1_ℳ29 (that is, ⟨Q⟩≃1 and ℳ_c = 29), in Figure 2. Each panel depicts a time snapshot of the projected density of the disc in the z-plane. Time is expressed in units of t_rot = 2πR/v_rot, calculated at R = 3 kpc, increasing from left to right. The colourbar is normalised to the initial central surface density of the disc, that is, Σ_c0 = Σ(R = 0, t = 0). Following the evolution of the galaxy, we can see that it forms dense structures in less than one rotation period. These structures are ∼kiloparsec long, uniformly spaced in azimuth, and pitched back as a result of differential rotation. For reference, we will refer to these structures interchangeably as feathers or filaments for the rest of the work. They are most prominent at radii R ∼ 2–4 kpc, which Figure 1 shows is the region of maximum Toomre instability, where Q dips well below unity. These features develop even in the absence of magnetic fields, spiral arms or feedback from stars, establishing that self-gravity alone is sufficient to produce them.

Fig. 2. Projected density of the Q1_ℳ29 model in the z-plane. The projected density is normalised by the initial central projected density of the galaxy, shown at 0, 0.4, and 0.8 rotational times (trot), from left to right. We see the gradual emergence of filamentary structures around R = 2–4 kpc.

Interestingly, we do not see the development of feathers in the region R ≤ 2 kpc (at least up to the time we have run the simulations) even though it also has values of Q ≤ 1. One possible explanation for why not is the difference in the amount of shear in the rotation curve at these smaller radii. We can quantify this through the shear parameter q = −(dlnΩ/dlnR), which is equal to zero for solid body rotation and one for a flat rotation curve. The value of q only rises to a value of 0.5 at R = R_c = 2 kpc, and thus the lack of feathers at smaller radii is potentially a result of the low shear in this region.

3.1.1. Dependence on ⟨Q⟩

We next examine the dependence of feather formation on the initial ⟨Q⟩. Similar to Figure 2, we show the z projections of disc density in Figure 3 for the ⟨Q⟩ = 1,  2,  3 (from left to right) runs at fixed ℳ_c = 29 at identical times t = 0.80 t_rot. We can immediately see the stark difference between the run with ⟨Q⟩ = 1 and the ones that have a higher ⟨Q⟩. For ⟨Q⟩∈{2, 3}, we see diffuse and faint structures in the disc. This is in contrast to the dense network of filaments that emerge in the ⟨Q⟩ = 1 simulation. We do see development of dense features in the ⟨Q⟩ = 2 run at times later than that shown in the Figure, but the rate of growth of this run is an order of magnitude slower than the ⟨Q⟩ = 1 run, and thus it does not reach significant density contrasts until much later than the other runs. We discuss this case further in Appendix C, and for the remainder of the main body of this paper we exclusively focus on the ⟨Q⟩ = 1 runs that form dense structures within a single rotation.

Fig. 3. Same as Figure 2, but for runs with ⟨Q⟩ = 1,  2,  3 (from left to right) at t = 0.8 trot. We see that the filamentary structures only emerge for the ⟨Q⟩ = 1 run.

3.1.2. Dependence on ℳ_c

We now examine the runs with ⟨Q⟩ = 1 but varying ℳ_c. Figure 4 is similar to Figure 3 in that it shows the z projections of the disc density for different runs, but in this case with a common value of ⟨Q⟩ = 1 but varying ℳ_c ∈ {14, 21, 29, 36}. All runs shown are at time t = 0.8 t_rot except for the ℳ36 run, which we show at t = 0.54 t_rot because it evolves faster than the other cases and collapses beyond our ability to resolve by 0.8 t_rot. We see that all the runs form filaments, but that there is a systematic and striking change of feather morphology between the panels. With increasing ℳ_c, the feathers are thinner and more closely-spaced. The ℳ_c = 29 and 36 cases also exhibit formation of dense clumps, in contrast to the ℳ_c = 14 and 21 runs where we see uniform larger-scale radially extended filaments.

Fig. 4. Same as Figure 3, but for runs with ℳc = 14,  21,  29,  36 (from left to right) and fixed ⟨Q⟩ = 1. All the discs are unstable to filament formation, but filament morphology differs starkly between the different runs, with higher ℳc producing more and thinner filaments. The filaments also tend to be shorter at higher ℳc.

The filaments we see in Figure 4 are comparable to the dense filamentary network observed in nearby spiral galaxies with JWST/MIRI observations (Thilker et al. 2023; Meidt et al. 2023) and the lattice or elongated extinction features reported in HST observations (La Vigne et al. 2006). In addition to the lattice-like morphology, where the filaments are present throughout the azimuth of the galaxy, the ℳ_c = 36 case has the filaments arranged in the pattern of spiral arms. This arrangement of filaments is similar to “spurs” (see La Vigne et al. 2006) that are filamentary features connected to spiral lanes. A notable feature of spurs is that they emanate from spiral arms whose pitch angles are lower compared to the pitch angle of spurs. This is also reproduced in the ℳ_c = 36 case. We defer a more quantitative comparison to observations to Section 4.2.

3.2. Feather formation timescales

After establishing that the galaxies with ⟨Q⟩ = 1 forms feathers, we now quantify the timescales of their formation. For this purpose we define the logarithmic surface density,

$$\begin{array}{c}\hfill \eta =ln\left(\frac{\mathrm{\Sigma }}{⟨\mathrm{\Sigma }⟩}\right),\end{array}$$(8)

where ⟨Σ⟩ is the average surface density of some region of interest. We take our region of interest for this purpose to be an annulus 1 kpc thick centred at R = 3 kpc, which is the region where feathers appear (as seen in Figure 4). We use the logarithmic surface density in order to give equal weight to the under- and over-dense regions that differ by a similar factor relative to the average. The variance in the value of η over the region of interest, which we denote σ_η, therefore describes the strength of density variations in that region. We show the time evolution of this quantity in Figure 5. The left panel is for runs with varying ⟨Q⟩ and constant ℳ_c = 29, while the right panel is for the group with ⟨Q⟩ = 1 but varying ℳ_c. The solid lines in both the panels depict the data from the simulations and the dashed lines in the second panel are empirical fits to the data that we describe below.

Fig. 5. Time evolution of the variance in the logarithmic surface density, ση, where η is defined in Equation 8. The panel on the left is for runs with varying ⟨Q⟩ and fixed ℳc and the panel on the right for runs with varying ℳc but fixed ⟨Q⟩, as indicated in the legends. The solid lines represent the simulation data and the dotted lines are piecewise fits with a constant and an exponential part. We see that the ⟨Q⟩> 1 runs are stable, while there is exponential growth in all the runs with ⟨Q⟩ = 1. For the unstable runs, we see that growth is faster for higher ℳc. We also observe a decrease in the time of the onset of the instability with increasing ℳc. We do not provide an analytical fit to the ℳc = 14 case because of its complex behaviour.

Starting with the left panel in Figure 5 which shows runs with varying ⟨Q⟩, we see that the variation in the clumping factor rises initially due to the turbulent velocity field causing fluctuations in the density field, and that this rise is similar in all runs. After this initial rise the value of σ_η stabilises, but at this point the runs begin to diverge. In the ⟨Q⟩ = 1 run the variations begin to rise again after a short stagnation period, marking the onset of instability. By contrast the runs with ⟨Q⟩ = 2,  3 show gradual decline followed by a phase of slow growth after around a rotation period – see Appendix C for further discussion of the long-time behaviour of the ⟨Q⟩ = 2 run. For the group of runs with ⟨Q⟩ = 1 but different ℳ_c (right panel), we see an evolution similar to the Q1_ℳ29 case for all but one of the runs. That is, an initial jump, followed by a period during which σ_η remains roughly constant, and then an epoch of exponential growth. The ℳ_c = 14 differs slightly from this pattern, showing a dip in the midst of its steady period before recovering and then embarking on exponential growth.

While the runs show similar qualitative behaviour, there are quantitative differences. First, we see that the disc develops stronger density fluctuations (higher σ_η) with increasing ℳ_c, despite the initial velocity perturbations being unchanged. Second, the exponential growth phase of the instability is systematically earlier for increasing ℳ_c. Third, the slopes of the rising curves are higher with increasing ℳ_c. We quantify these trends by carrying out a simple χ² fit for the time-evolution of σ_η at times t ≥ 0.15 t_rot, which is approximately the time of the initial relaxation phase, to an empirically-motivated piecewise function of the form

$$\begin{array}{c}\hfill {\sigma }_{\eta }(t)=\{\begin{array}{cc}A{e}^{\omega t_{\mathrm{onset}}}\hfill & t\le t_{\mathrm{onset}},\hfill \\ A{e}^{\omega t}\hfill & t>t_{\mathrm{onset}},\hfill \end{array}\end{array}$$(9)

where A sets the amplitude of σ_η during the steady phase and ω and t_onset quantify the growth rate and the time of the onset of the instability, respectively. We report the best-fitting parameters and their uncertainties in Table 2 and plot the fits as the dashed lines in the right-hand panel of Figure 5; the figure demonstrates that these functional forms fit the data very well, particularly in the time after t_onset. We do not fit for the ℳ_c = 14 run since it shows a more complex behaviour compared to the other three runs.

For our best fits, we find $t_{\mathrm{onset}}\propto {\mathcal{M}}_{\mathrm{c}}^{-1}$ to within 5%. Compared to this dependence of the onset time, we find that the growth rate only weakly depends on ℳ_c, with the highest growth rate $\omega =3.06\pm 0.08{t_{\mathrm{rot}}}^{-1}=0.49\pm 0.01\mathrm{\Omega }$ for the ℳ_c = 36 case and the slowest growth rate $\omega =2.41\pm 0.04{t_{\mathrm{rot}}}^{-1}=0.38\pm 0.01\mathrm{\Omega }$ for the ℳ_c = 21 run. A decrease in the growth rate is expected as the galaxy gets thicker since the gravitational potential is smeared out in the thick disc case when compared to the thin disc case of equal Q value (Meidt 2022). To check the robustness of this result against resolution, we carry out a resolution study on the Q1_ℳ29 run, which we discuss in Appendix D. The analysis there shows that the run is converged with respect to both the onset time of the instability and the growth rate.

Since the growth rate of density contrast depends on ℳ_c, it is natural to ask whether the difference in morphology visible in Figure 4 is solely the result of the instability being at different evolutionary stages. To answer this question in Figure 6 we again show projections of all the ⟨Q⟩ = 1 runs, but now choosing times when the different runs are at a common value of log₁₀σ_η = −0.3. The figure shows that, even though we are now examining comparable amplitudes of the instability, substantial morphological differences remain between the runs with differing ℳ_c. We quantify these morphological differences in the next subsection.

Fig. 6. Same as Figure 4, but now showing the different runs at different times. The time is chosen such that the standard deviation of the logarithmic surface density log10ση = −0.3 in all runs. We see systematic differences in the filament morphology even when the filament forming regions have similar density contrasts. The filaments are large scale and uniform in the ℳc = 14 and 21 runs compared to the more complex and finely spaced in the ℳc = 29 and 36 cases.

3.3. Feather geometry

In order to better visualise the morphology of the filamentary structures formed in our galaxies, we project the galaxy in the (lnR, θ) plane, where R is the galactocentric radius in kpc and θ is the azimuth. This allows easy identification of spiral features, since spirals of the form R = e^θtanα appear as straight lines of slope tanα in the (lnR, θ) plane. It also makes the azimuthal periodicity of the features more apparent to the eye. We show this projection in Figure 6 for the same runs and times as in Figure 4. Since our galaxy rotates clockwise, the positively-sloped features visible in the plot correspond to trailing spiral features (tips pointing opposite to galactic rotation). We can immediately see that the feathers in the ℳ_c = 14 and 21 runs appear as straight lines of nearly constant slope with a periodicity in the azimuthal direction. The complex morphology of the ℳ_c = 29 and 36 runs appear as thinner spirals with a larger variation in their orientations and lengths compared to the ℳ_c = 14 and 21 cases. In these runs we also see distinctive “spur”-like arrangements of the feathers located periodically along spirals, which are visible for example in the distinctive alternating pattern from θ ≈ −1 to 0 in the ℳ_c = 36 run.

To quantify these morphological differences, we take “tilted” 2D Fourier transforms (Savchenko 2012; Puerari et al. 2014) of the windows shown in Figure 7. We define

Fig. 7. Normalised surface density Σ/Σc0 in the (lnR, θ) plane for the same runs and times as shown in Figure 6. The radial range in this figure is chosen to highlight the annulus from 2 − 4 kpc where feather formation is strongest. We can see the feathers appear as sloped lines that are periodic in azimuth, with length, thickness, orientation and spacing that vary with ℳc. The associated time evolution movies are available online.

$$\begin{array}{c}\hfill A(m,p)={\displaystyle {\int }_{u_{\min }}^{u_{\max }}}{\displaystyle {\int }_0^{2\pi }}\eta (\theta ,u)e^{i(m\theta +pu)}\mathrm{d}\theta \mathrm{d}u,\end{array}$$(10)

where u = lnR, η is the logarithmic surface density (defined in Equation (8)), and (u_min, u_max) denote the radial bounds of our annuli. The term e^{i(mθ + pu)} can be viewed as a filter function that picks up m-armed logarithmic spirals with a pitch angle of tanα = −m/p. We show the amplitude |A(m, p)| normalised to the mean value $\overline{A}(m,p)=\iint A(m,p)\mathrm{d}m\mathrm{d}p/\iint \mathrm{d}m\mathrm{d}p$ in Figure 8, plotted with m on the abscissa and p on the ordinate. Different panels show runs with varying ℳ_c, increasing from left to right. We see that there is maximum in all the runs, and from the position of this maximum we can read off the preferred azimuthal scale and orientation of the filaments. Shifts in the position of the maximum and the shape of the distribution around it reflect the morphological differences seen between the feathers of different runs. We see that as ℳ_c increases, power shifts to higher m and more positive p, and that the distribution in the (m, p) plane broadens. These shifts mirror the more tightly packed and morphologically-complex filamentary pattern we see in going from low to high ℳ_c in Figure 6.

Fig. 8. Normalised amplitude of the 2D Fourier transform of the galactic annuli shown in Figure 7. Panels show runs with varying ℳc (increasing from left to right). We can see that as ℳc increases the maximum shifts to higher m and the width of the distribution increases.

3.3.1. Pitch angle

The 2D Fourier transform contains information on both the arm number and the pitch angle of the perturbations. To extract the latter we follow the approach of Savchenko (2012) by finding the m-mode of maximum power and then examining the range of p value, recalling that the pitch angle α = tan⁻¹(−m/p). Quantitatively, for each simulation we define m_max as the value of m at which the maximum of A(m, p) occurs. We then fit A(m_max, p) with a Gaussian functional form in p, which provides a reasonably good description of the data. We demonstrate this in Figure 9, which shows A(m_max, p) versus p, and our best Gaussian fits to it, for the same four snapshots shown in Figure 7.

Fig. 9. Function A(mmax, p) for runs with ⟨Q⟩ = 1 and varying ℳc (as indicated in the panels), along with their Gaussian fits. Black points show the amplitude of the Fourier transform A(mmax, p) as a function of p (cf., Equation (10)), where mmax is the value of m for which A(m, p) has its maximum. Orange lines show Gaussian fits to the data, with the amplitude A, mean μ, and dispersion σ of each fit as indicated in the legend.

We carry out this procedure for every snapshot of runs Q1_ℳ14 to Q1_ℳ36, extracting the amplitude A, mean μ, and dispersion σ of the best-fitting Gaussian at each time. We then take α = tan⁻¹(−m_max/μ) as our best estimate for the pitch angle. For the uncertainty in α, we take the dispersion of our fit as an estimate for the uncertainty in p, that is Δp = σ, and then propagate it to the uncertainty in pitch angle via standard error propagation, which gives

$$\begin{array}{c}\hfill \mathrm{\Delta }\alpha =\left|\frac{\mathrm{d}({tan}^{-1}(-m_{\max }/p))}{\mathrm{d}(-m_{\max }/p)}\mathrm{\Delta }(-m_{\max }/p)\right|=\left|{cos}^2\alpha \frac{m_{\max }}{{\mu }^2}\sigma \right|.\end{array}$$(11)

We show the pitch angle as a function of time after the onset of the instability that we derive from this procedure in Figure 10. The upper panel is for the ℳ_c = 14 and 21 runs and the lower panel for the ℳ_c = 29 and 36 cases. We see that, independent of ℳ_c, the average pitch angle of all the runs decreases with time. This is also seen visually in the time evolution movies of the surface densities in the (lnR, θ) plane (available online). This is expected for material arms being sheared by the differential rotation of the galaxy, since more rapid rotation of the inner regions winds up the arms. In terms of Figure 7, this corresponds to a systematic decrease in the slopes of the spiral features with time. This decrease is much more uniform in the ℳ_c = 14 and 21 cases when compared to the other two. Moreover, the spread in pitch angle systematically decreases with time for the ℳ_c = 14 and 21 runs but remains nearly constant, and a factor of two larger, for the ℳ_c = 29 and 36 runs. This difference reflects the more complex and clumpy morphology visible in the higher ℳ_c runs.

Fig. 10. Time evolution of the pitch angle, α = tan−1(−m/p), of the filaments in the runs with ⟨Q⟩ = 1 and varying ℳc. The upper panel shows the ℳc = 14 and 21 runs, and the lower shows ℳc = 29 and 36. All panels show results after the onset of the instability, and all times are normalised by the rotation period of the galaxy at R = 3 kpc. Circular points represent the mean pitch angle and error bars show the standard deviation. The dashed line shows α = 90°, that marks the transition between the trailing spirals with a positive pitch angles to leading spiral with negative ones. We see that all the runs show a decreasing trend in the pitch angle as the feathers get sheared due to differential rotation of the galaxy.

3.3.2. Spacing

To quantify feather spacing we integrate the 2D transform over the p direction to yield A(m) = ∫A(m, p) dp. We find that A(m) does not vary appreciably after t_onset, and thus for simplicity we carry out our analysis only on the snapshots shown in Figure 7; results for other snapshots are qualitatively identical. As with our analysis of the pitch angle, we fit an empirically chosen functional form to A(m); since m is a positive-definite quantity, we use a lognormal-like rather than a normal form,

$$\begin{array}{c}\hfill A(m)=\frac{A_m}{m{\sigma }_m\sqrt{2\pi }}exp\left[\frac{{(logm-{\mu }_m)}^2}{2{\sigma }_m^2}\right],\end{array}$$(12)

where A_m is the amplitude, μ_m is the location of the peak in log m, and σ_m is the width in log m. In Figure 11 we show A(m) along with our best empirical fits, the parameters for which we report in Table 3. The figure demonstrates that our chosen functional form provides a reasonable description of the data. From the empirical fit we extract the median, given by m_fil = exp(μ_m), as a representation of the dominant azimuthal mode, and the distribution’s 16th and 84th percentiles given by exp(μ_m ∓ σ_m), as its lower and upper spread.

Fig. 11. Amplitude of the 2D Fourier transform A(m, p) of the logarithmic surface density (see Equation (10)) in the region 2–4 kpc, integrated over p. Different panels are for runs with ⟨Q⟩ = 1 and varying ℳc, evaluated at the same times as in Figure 6. The grey curves are simulation data and the orange ones are the empirical fits to the data (see text for details). We can see how the power shifts toward higher values of m with increasing ℳc.

We show m_fil as a function of ℳ_c for our four runs in Figure 12. The right vertical axis of this plot shows the physical spacing corresponding to this value of m_fil, which is λ_fil = 2πR_gal/m_fil, where R_gal = 3 kpc is the mean radius of the annulus we use in our analysis. We find that the m_fil varies approximately linearly with ℳ_c. A χ² fit to this relationship of the form

$$\begin{array}{c}\hfill m_{\mathrm{fil}}=a{\mathcal{M}}_{\mathrm{c}}+b,\end{array}$$(13)

Fig. 12. Azimuthal mode of the filaments (mfil) and various instabilities, shown with the ℳc ( ∼ vc/σ) of our galaxies. Black circular points with error bars show the mean and 16th to 84th percentile range of the mfil and corresponding physical spacing λfil for the snapshots shown in Figure 7. The dashed line is a linear fit to the median values of mfil (Equation (13)). Purple, red, and orange points with error bars represent the radial averages ⟨λ⟩ and standard deviations σλ (Equation (14) and Equation (15)) of the wavelength of the most unstable Toomre mode, the Jeans length and the Swing amplification length – see main text for details. These points are slightly offset from their overlapping ℳc values for visibility. The simulations show an increase in the median value of mfil and a broadening of the 16th to 84th percentile range with increasing ℳc. The median filament spacing agrees well with the Toomre length of the ℳc = 21, 29 and 36 runs, and with the swing length for ℳc = 14 run.

yields a = 2.13 ± 0.02 and b = −15.5 ± 0.6, and we show this fit via the black dashed line in Figure 12.

3.4. Feather density structure and evolution

Finally, we examine the three-dimensional density structure of the feathers. Our goal here is understanding how this depends on ℳ_c, and to see whether the feathers are likely to be sites for star formation. Figure 13, where we shows the mass-weighted probability density function (PDF) of log₁₀ of ρ/ρ_mid, where ρ is the density of the gas and ρ_mid is the radial average of the initial density of the region of interest. We do this for the cylindrical region bounded by 1 kpc thick radial annuli centred at 3 kpc. The left panel is for all the “feather forming” runs with ⟨Q⟩ = 1 and varying ℳ_c at times when they are similarly evolved; see snapshots in Figure 7. In the right panel we show the PDF at several different times for the Q1_ℳ36 run.

Fig. 13. Mass-weighted probability density function (PDF) of the log10 of the density normalised by the radially averaged mid-plane density, for the snapshots shown in Figure 7. The grey dashed line depicts the Jeans density (see text for details), at which we expect numerical effects to dominate. The left panel shows all the runs at times similar to Figure 7, and the right panel shows the time evolution of the Q1_ℳ36 case leading up to the time in the left panel. We see that all the runs in the left panel form high-density tails that represent feathers and that the densest end of the tail is reaching the limit imposed by the resolution of our simulations.

In this plot, the stratified medium of the background galaxy corresponds to the parts of the PDF with ρ/ρ_mid ≤ 1. Filaments and further dense regions inside them are represented in the higher-density tails. Looking at the left panel, we can see that all of our runs form the high-density tail, which is expected from feather-forming galaxies. However, we also see a systematic trend that the higher ℳ_c runs have longer tails extending towards higher densities. This densest part contains little mass and is largely de-coupled from the initial phase of the instability, and is instead driven by self-gravity.

To confirm the importance of self-gravity, in the right panel we add a vertical dashed line marking the Jeans density for the Q1_ℳ36 run, defined as ρ_Jeans = (c_s/4Δx_min)²(π/G), which is the density at which self-gravitating structures become unresolved on the grid scale (Truelove et al. 1997, formally, our expression is the density for which the Jeans length is equal to 4Δx_min). Here, c_s is the sound speed and Δx_min is the length of the grid cells. We see that the high-density tail of the PDF does not reach a steady-state, and instead rises until it reaches the Jeans density at our resolution. This likely indicates that the high-density tail of the PDF is made up of gas in a state of runaway collapse, where star formation is likely to occur in the future. A detailed study of how this star formation proceeds is beyond the scope of this work.

4. Discussion

In this section we discuss our results in the light of existing observational, numerical and theoretical works and point out limitations of our approach. In Section 4.1 we comment on the nature of instability that is responsible for feather formation, and in Section 4.2 we compare the filament spacing found in our simulations with observations of nearby galaxies. Lastly, in Section 4.3, we outline the limitations of our simulations and hint on future avenues.

4.1. Nature of instability

In order to better understand the physical mechanism responsible for feather formation in our simulations, we compare the mean azimuthal spacing of the filaments calculated in Section 3.3.2 to the length scales of various instabilities expected to be present in the disc: the wavelength of the most unstable Toomre mode (Toomre 1964), the Jeans length (Jeans 1902) in the galactic mid-plane, and the length scale of the Swing amplification mechanism (Fuchs 2001; Kim & Ostriker 2002; Meidt & van der Wel 2024). These are given by λ_Toomre = 2c_s²/GΣ, λ_Jeans = (πc_s²/Gρ_mid)^1/2 and λ_Swing = λ_JeansX_J, respectively, where ρ_mid is the mid-plane density of the galaxy and X_J measures the characteristic length of the swing amplification in terms of the Jeans length and has been found to be ≈3 for a wide range of perturbing wavelengths (Meidt & van der Wel 2024); we adopt X_J = 3 exactly for our comparison.

These length scales all depend on the local properties of the disc, and thus vary with position and time within our analysis region. In order to reduce this complexity to a single parameter for each run, we take averages over our region of interest using the initial configurations in the simulations. Specifically, for each length scale λ, we defined

$$\begin{array}{c}\hfill ⟨\lambda ⟩=\frac{1}{R_{\mathrm{out}}-R_{\mathrm{in}}}{\displaystyle {\int }_{R_{\mathrm{in}}}^{R_{\mathrm{out}}}}\lambda (R,t=0)\mathrm{d}R,\end{array}$$(14)

where R_in = 2 kpc and R_out = 4 kpc are the bounding radii of our annulus and λ(R, t = 0) represents the length scale of interest evaluated for the initial disc state at radius R. We also compute the standard deviations of λ over this radial range, defined as

$$\begin{array}{c}\hfill {\sigma }_{\lambda }={\left\{\frac{1}{R_{\mathrm{out}}-R_{\mathrm{in}}}{\displaystyle {\int }_{R_{\mathrm{in}}}^{R_{\mathrm{out}}}}{[\lambda (R,t=0)-\langle \lambda \rangle ]}^2\mathrm{d}R\right\}}^{1/2}.\end{array}$$(15)

We plot ⟨λ⟩ for each of our candidate length scales – Jeans, Toomre, and Swing – on top of values of λ_fil measured from our simulations in Figure 12; in this plot, the errorbars on the points are the standard deviations of the corresponding quantities. Comparing the theoretically-predicted length scales with the measured values of λ_fil, first we note that except for the λ_Jeans in the ℳ_c = 14 case, all lie within the spread of the measured filament spacings, but that λ_Jeans is consistently the farthest from the measured value of λ_fil. The other two length scales are closer to the measurements, with λ_Toomre lying very close to λ_fil for the ℳ_c = 29 and 36 runs, and λ_Swing matching it very closely for the ℳ_c = 14 case. The ℳ_c = 21 run lies at a transition for these two behaviours.

An appealing interpretation of this trend is that filament formation is driven by the Toomre instability in the two higher ℳ_c cases, while Swing amplification dominates at lower ℳ_c. Such a transition is expected theoretically, because the threshold value of Q required to trigger Toomre instability increases with disc thickness (Wang et al. 2010), and the low ℳ_c discs are thicker than the high ℳ_c ones. Thus for our series of runs at fixed Q, we might expect Toomre instability to dominate in high ℳ_c discs that are thin enough to trigger Toomre instability, and for the somewhat slower-growing Swing mechanism to dominate in thicker discs that are Toomre-stable. Our results also appear to be quantitatively consistent with this hypothesis. While the exact Q threshold for Toomre instability in a finite-thickness disc is somewhat sensitive to the exact vertical distribution of the gas, Wang et al. (2010) find typical thresholds of ≃0.7 − 0.8 for discs with ℳ_c in the range 16–24. Examining Figure 1, it is clear that our discs do not quite reach such small Q values in the feather forming region from 2–4 kpc, and thus it seems likely that our ℳ_c = 14 and 21 cases are Toomre stable, consistent with our finding that Swing amplification dominates in these.

Swing amplification also offers an explanation for the structure formation seen at very late times in the ⟨Q⟩ = 2 run, which has a Q value well above the threshold for Toomre instability at all radii (Appendix C).

In observations it can be even more challenging to single out the mechanism that is directly responsible for filament formation since we do not have the option to simply “switch-off” physical processes that affect gas flows. Still, there is one observational study that analyses the emergence of galactic filaments as a consequence of gravitational instability (Meidt et al. 2023). They find that their four filament-forming galaxies have Q ≥ 1, that is, above the more conservative Toomre threshold for razor-thin discs. While the galaxies could still form filaments via the Swing amplification, we point out that the inclusion of the major stellar component (Leroy et al. 2021) will bring down this Q value potentially making them gravitationally unstable (Romeo & Falstad 2013). Meidt et al. (2023) also find that the filament spacing in these galaxies is more comparable to the estimated Jeans length, rather than the Toomre length, which however, will also change with the inclusion of the stellar component. Given the high degree of complexity of the filaments observed in these galaxies, we believe that it is more likely that they are formed due to the Toomre-like galactic scale instability. We compare the filament spacings from observations with our simulations in the next section.

4.2. Comparison with observations

We now compare the filament spacing extracted from our simulations with those observed in nearby galaxies, adopting a fiducial radius of 3 kpc. We plot the filament spacing in our simulations together with observed filament spacings as a function of ℳ_c in Figure 14. In this plot we show the measured filament spacings from our simulations as black points with error bars, and the values shown match the right-axis values indicated in Figure 12; similarly, the dashed line shows the empirical curve given by Equation (13), and is identical to the dashed black line in Figure 12.

Fig. 14. Physical spacing of filaments as a function of the rotational Mach number, ℳc (∼vc/σ), in our simulations (black circles; cf., Section 3.3.2) compared with observations (yellow and blue diamonds). The yellow diamonds are JWST/MIRI observations (Meidt et al. 2023) and the blue diamonds are HST observations (La Vigne et al. 2006). We see that, except for NGC 1365 and NGC 4548, the filament spacing extracted from our simulations agrees within the uncertainty ranges of the ones reported in the observations.

We draw the observations (shown as diamonds in Figure 14) from two sources: the JWST/MIRI 11.3 μm observations reported by Meidt et al. (2023, their Table 1), which we plot in yellow, and the HST archival data of La Vigne et al (2006, shown in blue), and with the feather spacing of these sources taken from the tabulation of Mandowara et al. (2022, their Table 4). Thilker et al. (2023) show that there is good agreement between filament catalogues derived from optical and IR data, and thus for the analysis that follows we will not differentiate between the two data sources. We extract the filament spacings and values of ℳ_c that we plot as follows. For the JWST observations, our central points correspond to the reported 50th percentile of the filament spacing, and our lower and upper error bars to the 16th and 84th percentiles, respectively, for the region R ≤ 4 kpc; however, the choice of radius has little effect since Meidt et al. find no appreciable variation in filament spacing with galacto-centric radius. For the HST observations, Mandowara et al. report filament spacings at various galacto-centric radii, so to estimate the spacing at R = 3 kpc we assume that each galaxy is characterised by a single, global m_fil at all radii, which we extract by carrying out a simple χ² fit for m_fil to the radius-dependent values of λ_fil reported in Mandowara et al.’s Table 4; we then compute λ_fil = 2πR/m_fil at R = 3 kpc. For ℳ_c in both sets of observations, we assume a fixed c_s = 7 km s⁻¹ for all galaxies, characteristic of the warm neutral medium (McClure-Griffiths et al. 2023), and take our circular velocities v_c from the saturated circular velocities of all galaxies calculated from CO (2–1) emission maps taken from Lang et al. (2020, their Table 4); we derive error bars in ℳ_c from their quoted uncertainties in v_c, assuming no uncertainty in c_s. The only exceptions to this are the galaxies IC 5332 and NGC 5055, which are not included in this sample. For NGC 5055, we use the H I rotation curve velocity from McGaugh & Schombert (2015) and for IC 5332 we adopt v_c = 120 ± 20.6 km s⁻¹, derived from a semi-empirical model that links the rotation curve to the stellar mass (Meidt et al. 2018), where the spread is propagated from the uncertainty in the galaxy’s inclination angle (Leroy et al. 2021). These properties, along with the star formation rates (SFR) and total stellar masses (M_stellar), are tabulated in Table 4. The SFR and the M_stellar are taken from Leroy et al. (2021) for all but NGC 5055, which is taken from Leroy et al. (2019). Our sample consists of late-type main sequence star-forming galaxies.

From Figure 14, we can see that the filament spacing from our simulations agrees well with the observations despite the simplicity of our models. All but NGC 1365 fall within a factor of ∼2 of the median of our filament spacing. The trend of decreasing filament spacing with increasing ℳ_c predicted by the simulations is reproduced in the observational data, albeit with a fair amount of scatter since the data are sparse at both high and low ℳ_c. The comparison is currently limited by the unavailability of kinematic data for many galaxies whose filament spacings have already been tabulated.

4.3. Limitations and future work

Our simulations consist of self-gravitating, isothermal gas. As a result, we miss some key physical processes pervasive in the ISM that are dynamically relevant on galactic scales. These include magnetic fields, cosmic rays, and feedback processes such as supernovae and winds from massive stars. Moreover, we use static analytical potentials to approximate the gravitational interaction of stars and dark matter with the gas. Therefore, we miss the back-reaction of the former with the latter and vice-versa.

The inclusion of each of these aforementioned processes may affect feather formation and structure, acting as an additional source or damping mechanism for the perturbations or directly impacting the nature of the instability. For example, magnetic fields can have a stabilising effect due to the additional magnetic pressure. However, they are also known to destabilise the disc via Parker instability (Parker 1966; Körtgen et al. 2019; Arora et al. 2024) and the magneto-Jeans instability (Kim & Ostriker 2002). Further, since we expect star formation in feathers, energy and momentum injection due to feedback from stars may also affect their morphology. Finally, the presence of a live stellar component may destabilise the disc (Romeo & Wiegert 2011), which might lead to a shift in the Q thresholds for feather formation.

We point out, however, that the omission of these physical processes is intentional and motivated by our goal of conducting a clean numerical experiment to determine the role of gravity, rotation, and pressure alone in the formation of filaments, and to determine their morphological properties. Future studies will aim to incorporate more physical processes in order to study feather formation in more realistic galaxies. Future studies may also investigate the distribution of filament properties such as their masses, their velocity dispersion, and aspect ratios. The relative orientation of the filaments with respect to the local magnetic field is also of interest. These properties could help us gain insight into the variations of the initial conditions of star-forming regions residing inside feathers (for example; Thilker et al. 2023).

5. Summary

We simulate isothermal, self-gravitating, isolated disc galaxies that are initialised in equilibrium, with the aim of capturing structure formation via global gravitational instability. Our models are intentionally simplified, allowing us to isolate the dependence of structure formation on two dimensionless parameters: the radially averaged Toomre parameter ⟨Q⟩ and the rotational Mach number, ℳ_c. Our main conclusions are summarised below.

• Galaxies with ⟨Q⟩ = 1 in the region R ≥ 2 kpc form dense filaments and feathers within a single rotation, while galaxies with ⟨Q⟩ = 2, 3 fail to form them within similar evolution times. This demonstrates that filaments can form solely due to galactic-scale instability driven by self-gravity overcoming the stabilising effect of rotation and thermal pressure. Thus filament formation occurs without the need for an external stellar spiral potential, magnetic fields, or supernova feedback. This does not prove that these other mechanisms cannot alternatively generate filaments, or that they do not modify the properties of filaments formed via gravitational instability; it simply establishes that gravitational instability alone is sufficient.

• The dense filaments that form in our simulations are kiloparsec long and semi-regularly spaced in azimuth. Their morphology and their density contrasts vary systematically with ℳ_c: morphological complexity increases with ℳ_c and filament spacing varying as roughly ${\mathcal{M}}_{\mathrm{c}}^{-1}$. Filament pitch angle decreases with time in all runs due to shearing by differential rotation, but the dispersion in pitch angle is larger at higher ℳ_c as well. Filaments with a higher ℳ_c also show density contrasts with longer tails that extend towards higher densities.

• The time required for onset of the filament forming instability varies ${\mathcal{M}}_{\mathrm{c}}^{-1}$ as well; by contrast, the growth rate of variations in the density depends only weakly on ℳ_c, and is typically 30–50% of the galactic angular velocity. Highest for the run with higher ℳ_c.

• The filament spacings and their dependence on ℳ_c found in our simulations are in good agreement with the characteristic length scales of the Toomre instability at higher ℳ_c, and with the length scale of the Swing amplification mechanism for lower ℳ_c. They are also consistent within observational uncertainties with values measured in nearby galaxies.

Our results demonstrate that the filamentary ISM structures revealed by modern telescopes and facilities such as the JWST and ALMA in nearby galaxies with unprecedented resolution can at least potentially be understood as arising from very simple physical processes. Our numerical experiments show that self-gravity alone is sufficient to produce structures whose characteristics are in quantitative agreement with the observations. However, we are yet to see how this might change with the inclusion of other mechanisms such as magnetic fields, cosmic rays, or stellar feedback. In a forthcoming companion study we will quantify some of these effects.

Data availability

Movies associated to Fig. 7 are available at https://www.aanda.org

Acknowledgments

R. A. acknowledges funding by the Hamburg State Graduate Funding Program (HmbNFG) and heartfully thanks Sharon Meidt, Shivan Khullar, and Chris Matzner for insightful discussions and their enthusiasm. C. F. acknowledges funding provided by the Australian Research Council (Discovery Project grants DP230102280 and DP250101526), and the Australia-Germany Joint Research Cooperation Scheme (UA-DAAD). M. R. K. acknowledges support from the Australian Research Council through Laureate Fellowship FL220100020. We further acknowledge high-performance computing resources provided by the Leibniz Rechenzentrum and the Gauss Centre for Supercomputing (grants pr32lo, pr48pi and GCS Large-scale project 10391), the Australian National Computational Infrastructure (grants ek9 and jh2) and the Pawsey Supercomputing Centre (projects pawsey0807 and pawsey0810) in the framework of the National Computational Merit Allocation Scheme and the ANU Merit Allocation Scheme. The simulation software, FLASH, was in part developed by the Flash Centre for Computational Science at the University of Chicago and the Department of Physics and Astronomy at the University of Rochester.

References

Appendix A: Dark matter potential

Here, we determine the stellar plus dark matter potential ϕ_dm such that our disc is in radial equilibrium. In a reference frame that co-rotates with the disc angular frequency Ω at a given position (R, z), radial equilibrium demands that the sum of the gravitational, pressure, and centrifugal forces vanish. We can express this condition as

$$\begin{array}{c}\hfill \rho \frac{\partial }{\partial R}({\varphi }_{\mathrm{dm}}+{\varphi }_{\mathrm{g}})+{c_{\mathrm{s}}}^2\frac{\partial \rho }{\partial R}-\rho {\mathrm{\Omega }}^{2}R=0,\end{array}$$(A.1)

where ϕ_g is the potential due to disc self-gravity, which is given by Wang et al. (2010, their equation 29)

$$\begin{array}{c}\hfill {\varphi }_{\mathrm{g}}=4\pi G{\mathrm{\Sigma }}_{°}{R}_{\mathrm{d}}{y}^2[I_0(y)K_0(y)-I_1(y)K_1(y)],\end{array}$$(A.2)

where y = R/2R_d and I₀, K₀, I₁, and K₁ are the modified Bessel functions of the first and second kind, of order zero and one.² To proceed we must now make some assumptions about the vertical structure of the disc and potential. We start by assuming that the radial variation in density is approximately exponential, ρ(R)∝exp(−R/R_d). This assumption is exact only if the disc vertical profile is independent of radius, but is true approximately as long as the disc scale height does not vary on scales ≪R_d; since this term is highly sub-dominant in any event (roughly 100× smaller than the centrifugal term), we do not seek a more accurate approximation. Finally, at the midplane we require that Ω = v_rot/R, but we must make an assumption about the variation of Ω in the z direction, or, equivalently, assume a shape for the centrifugal potential ψ, where ∂ψ/∂R = Ω²R. For this purpose we adopt the flattened centrifugal potential form suggested by Binney & Tremaine (1987),

$$\begin{array}{c}\hfill \psi =\frac{{v_{\mathrm{c}}}^2}{2}ln\left\{\frac{1}{R_{\mathrm{c}}^2}[R^2+R_{\mathrm{c}}^2+{\left(\frac{z}{q}\right)}^2]\right\},\end{array}$$(A.3)

where q is a dimensionless parameter that determines the amount of flattening. One can readily verify that this potential satisfies our condition that ∂ψ/∂R = Ω²R. Following Binney & Tremaine’s estimates for realistic galactic potentials, we adopt q = 0.7. With these assumptions, our equation of force balance reduces to

$$\begin{array}{c}\hfill \rho \frac{\partial }{\partial R}({\varphi }_{\mathrm{dm}}+{\varphi }_{\mathrm{g}}-{c_{\mathrm{s}}}^2\frac{R}{R_d}-\psi )=0,\end{array}$$(A.4)

which has the immediate solution

$$\begin{array}{c}\hfill {\varphi }_{\mathrm{dm}}=\psi -{\varphi }_{\mathrm{g}}+{c_{\mathrm{s}}}^2\frac{R}{R_{\mathrm{d}}}.\end{array}$$(A.5)

We have therefore succeeded in identifying the stellar plus dark matter potential required to place the disc in equilibrium.

Appendix B: Two dimensionless parameters

In this appendix we prove the assertion made in the main text that the fluid equations that we solve depend only upon the two dimensionless parameters ⟨Q⟩ and ℳ_c. And that these together with the additional dimensionless parameters describing the initial conditions that we do not vary – R_c/R_d, R_max/R_d, q, α, and β – fully determine the evolution of the system. To do so, we non-dimensionalise the fluid variables via the transformations

$$\begin{array}{cc}\hfill \mathit{r}& =\mathit{x}{R}_{\mathrm{d}},\hfill \end{array}$$(B.1)

$$\begin{array}{cc}\hfill t& =\tau \frac{R_{\mathrm{d}}}{c_{\mathrm{s}}},\hfill \end{array}$$(B.2)

$$\begin{array}{cc}\hfill \mathit{v}& =\mathit{u}{c}_{\mathrm{s}},\hfill \end{array}$$(B.3)

$$\begin{array}{cc}\hfill \rho & =a\frac{{\mathrm{\Sigma }}_{°}}{R_{\mathrm{d}}},\hfill \end{array}$$(B.4)

$$\begin{array}{cc}\hfill \varphi & =\frac{\varphi }{{c_{\mathrm{s}}}^2}.\hfill \end{array}$$(B.5)

where x, τ, u, a, and $\stackrel{\sim }{\varphi }$ are dimensionless length, time, velocity, density and the gravitational potential, respectively, scaled by the appropriate combinations of the disc scale length R_d, the sound speed c_s and the surface density scaling of our galaxies Σ_°. The equations of mass and momentum conservation plus the Poisson equation under this transformation are

$$\begin{array}{cc}\hfill \frac{\partial a}{\partial \tau }+\frac{1}{{\mathcal{M}}_{\mathrm{c}}}{\mathbf{\nabla }}_{\mathit{x}}(a\mathit{u})=& 0,\hfill \end{array}$$(B.6)

$$\begin{array}{cc}\hfill {\mathcal{M}}_{\mathrm{c}}\frac{\partial \mathit{u}}{\partial \tau }+(\mathit{u}·{\mathrm{\nabla }}_x)·\mathit{u}=& -\frac{{\mathrm{\nabla }}_{x}a}{a}-{\mathrm{\nabla }}_x({\varphi }_{\mathrm{g}}+{\varphi }_{\mathrm{dm}}),\hfill \end{array}$$(B.7)

$$\begin{array}{cc}\hfill {\mathrm{\nabla }}_x^2\varphi =& \xi \frac{{\mathcal{M}}_{\mathrm{c}}}{⟨Q⟩}a,\hfill \end{array}$$(B.8)

where ∇_x indicates gradient with respect to the dimensionless position variable x and ξ is a dimensionless constant given by

$$\begin{array}{cc}\hfill \xi & =\frac{4\sqrt{2}}{X_{\max }}{\displaystyle {\int }_0^{X_{\max }}}\frac{1}{exp[-X-\beta exp(-\alpha X)]}\frac{\sqrt{X^2+2X_{\mathrm{c}}^2}}{X^2+X_{\mathrm{c}}^2}dX\hfill \\ \hfill & \approx 19.43.\hfill \end{array}$$(B.9)

In the expression above, we have defined X = R/R_d as the dimensionless cylindrical radius, X_c and X_max are defined analogously, and the numerical evaluation is for our fixed values of these quantities, α, and β. Note that we have used the isothermal equation of state to substitute for the pressure (equation 3). Finally, we can rewrite the various components of the gravitational potential that contributed to $\stackrel{\sim }{\varphi }$ using the same non-dimensionalisation procedure. Substituting into equation A.5 we have

$$\begin{array}{cc}\hfill {\varphi }_{\mathrm{dm}}& =\psi -{\varphi }_{\mathrm{g}}+X,\hfill \end{array}$$(B.10)

$$\begin{array}{cc}\hfill \psi & =\frac{{\mathcal{M}}_{\mathrm{c}}^2}{2}ln\left\{X_{\mathrm{c}}^{-2}[X^2+X_{\mathrm{c}}^2+{\left(\frac{Z}{q}\right)}^2]\right\},\hfill \end{array}$$(B.11)

$$\begin{array}{cc}\hfill {\varphi }_{\mathrm{g}}& =\frac{{\mathcal{M}}_{\mathrm{c}}}{⟨Q⟩\xi }{X}^2[I_0\left(\frac{X}{2}\right)K_0\left(\frac{X}{2}\right)-I_1\left(\frac{X}{2}\right)K_1\left(\frac{X}{2}\right)],\hfill \end{array}$$(B.12)

where Z = z/R_d is the dimensionless vertical height. This completes the non-dimensionalisation of the system, and demonstrates that the results depend only on ⟨Q⟩, ℳ_c, X_c, X_max, q, α, and β as claimed.

Appendix C: Toomre ⟨Q⟩ = 2 run

Fig. C.1. Same as Figure 2, but for the longer time evolution of the Q = 2 run. We see that there are faint structures at early times in the first two panels. The third panel shows a structure with somewhat more density contrast, visually confirming the slow growth of an instability.

In this Appendix we present the long-term evolution of the ⟨Q⟩ = 2 run (Q2_ℳ29), which does not produce significant structure on the same ∼t_rot timescales as the ⟨Q⟩ = 1 runs presented in the main body. We continue this run to time t = 2.5 t_rot, and show the projected surface density as a function of time in Figure C.1. We see that the galaxy does develop structures with noticeable density contrast by t = 2.5 t_rot at R = 3 kpc, and that as in the ⟨Q⟩ = 1 cases, these structures are regularly spaced in the azimuth. However, there are also significant differences. First, the structures in the ⟨Q⟩ = 2 run clearly have lower pitch angles and density contrasts. Second, their azimuthal spacing is also different - counting the number of dense structures along the azimuth we estimate m_fil ≃ 15 to be the dominant mode of the filaments. This is a factor of two less than the Q1_ℳ29 run, that has the same ℳ_c value as the ⟨Q⟩ = 2 run. Moreover, the number of filaments in the Q2_ℳ29 case is around the value expected from the swing amplification mechanism, which is = 13.

Fig. C.2. Same as Figure 5, but only for Q1_ℳ29 and Q2_ℳ29 runs, where the latter was evolved to a later time than shown in the main text. We see the oscillatory peaks in the Q = 2 run, but also a steady, yet slow growth over time. This growth is visibly smaller than in the Q1_ℳ29 run.

To compare the growth rate of filaments of both the cases, we plot the time evolution of the variation of the logarithmic density contrast (σ_η) of the Q1_ℳ29 run and the Q2_ℳ29 run together in Figure C.2. This is similar to the left panel of the Figure 5, except for the x-axis that extends to longer times. The dashed lines show the empirical fit on the two runs. The Q1_ℳ29 run has the piecewise-fit of the form equation 9, identical to the one shown in the right panel of the Figure 5. We fit only the exponential part of the fit for the Q2_ℳ29 run after t = 0.55 t_rot, since we are interested exclusively in the rising part of the curve. As seen in the main text, the Q2_ℳ29 shows an initial decline in σ_η. However, from t ≈ 0.55 onwards, it shows periodic oscillations with a rising mean. From our empirical fit to this part of the curve, we find ωt_rot = 0.31, which is ≈8 times smaller than Q1_ℳ29 case.

Appendix D: Resolution study

Fig. D.1. Same as Figure 5, but for the Q1_ℳ29 run at three different resolutions. The solid line depicts the resolution used in the main text (Δx = 4.8 pc) and the other two lines are runs with resolutions reduced by a factor of two each.

For a resolution study, we perform two runs with resolutions that are a factor of two (Δx = 9.6 pc) and four (Δx = 19.2 pc) lower than the one used in the main text. Similar to Figure 5, we show the time evolution of the standard deviation of the logarithmic surface density (defined in Equation (8)) of the annuli R = 2 − 4 kpc of the three runs in Figure D.1. We see similar qualitative evolution in all the three cases as seen for all the ⟨Q⟩ = 1 runs. The steady value of σ_η after the initial relaxation relaxation phase increases with increasing resolution. This is expected since we resolve the initial turbulent velocity fluctuations and the resulting density fluctuations with more cells. However, this change in σ_η between runs with consecutive increasing resolutions is less than a factor of two. The time of the onset of exponential growth, as well as the slopes of the lines following the onset are very similar, especially for the two runs with the highest resolutions.

All Tables

All Figures

	Fig. 1. Initial radial Toomre-Q (see Equation (6)) profiles of our galaxies. The solid lines show Q as a function of galactocentric radius and the dot-dashed lines show ⟨Q⟩, the average of the solid curves over the disc region 0 ≤ R/kpc ≤ 8.
In the text

	Fig. 2. Projected density of the Q1_ℳ29 model in the z-plane. The projected density is normalised by the initial central projected density of the galaxy, shown at 0, 0.4, and 0.8 rotational times (trot), from left to right. We see the gradual emergence of filamentary structures around R = 2–4 kpc.
In the text

	Fig. 3. Same as Figure 2, but for runs with ⟨Q⟩ = 1,  2,  3 (from left to right) at t = 0.8 trot. We see that the filamentary structures only emerge for the ⟨Q⟩ = 1 run.
In the text

	Fig. 4. Same as Figure 3, but for runs with ℳc = 14,  21,  29,  36 (from left to right) and fixed ⟨Q⟩ = 1. All the discs are unstable to filament formation, but filament morphology differs starkly between the different runs, with higher ℳc producing more and thinner filaments. The filaments also tend to be shorter at higher ℳc.
In the text

Fig. 5. Time evolution of the variance in the logarithmic surface density, ση, where η is defined in Equation 8. The panel on the left is for runs with varying ⟨Q⟩ and fixed ℳc and the panel on the right for runs with varying ℳc but fixed ⟨Q⟩, as indicated in the legends. The solid lines represent the simulation data and the dotted lines are piecewise fits with a constant and an exponential part. We see that the ⟨Q⟩> 1 runs are stable, while there is exponential growth in all the runs with ⟨Q⟩ = 1. For the unstable runs, we see that growth is faster for higher ℳc. We also observe a decrease in the time of the onset of the instability with increasing ℳc. We do not provide an analytical fit to the ℳc = 14 case because of its complex behaviour.

In the text

Fig. 6. Same as Figure 4, but now showing the different runs at different times. The time is chosen such that the standard deviation of the logarithmic surface density log10ση = −0.3 in all runs. We see systematic differences in the filament morphology even when the filament forming regions have similar density contrasts. The filaments are large scale and uniform in the ℳc = 14 and 21 runs compared to the more complex and finely spaced in the ℳc = 29 and 36 cases.

In the text

Fig. 7. Normalised surface density Σ/Σc0 in the (lnR, θ) plane for the same runs and times as shown in Figure 6. The radial range in this figure is chosen to highlight the annulus from 2 − 4 kpc where feather formation is strongest. We can see the feathers appear as sloped lines that are periodic in azimuth, with length, thickness, orientation and spacing that vary with ℳc. The associated time evolution movies are available online.

In the text

	Fig. 8. Normalised amplitude of the 2D Fourier transform of the galactic annuli shown in Figure 7. Panels show runs with varying ℳc (increasing from left to right). We can see that as ℳc increases the maximum shifts to higher m and the width of the distribution increases.
In the text

Fig. 9. Function A(mmax, p) for runs with ⟨Q⟩ = 1 and varying ℳc (as indicated in the panels), along with their Gaussian fits. Black points show the amplitude of the Fourier transform A(mmax, p) as a function of p (cf., Equation (10)), where mmax is the value of m for which A(m, p) has its maximum. Orange lines show Gaussian fits to the data, with the amplitude A, mean μ, and dispersion σ of each fit as indicated in the legend.

In the text

Fig. 10. Time evolution of the pitch angle, α = tan−1(−m/p), of the filaments in the runs with ⟨Q⟩ = 1 and varying ℳc. The upper panel shows the ℳc = 14 and 21 runs, and the lower shows ℳc = 29 and 36. All panels show results after the onset of the instability, and all times are normalised by the rotation period of the galaxy at R = 3 kpc. Circular points represent the mean pitch angle and error bars show the standard deviation. The dashed line shows α = 90°, that marks the transition between the trailing spirals with a positive pitch angles to leading spiral with negative ones. We see that all the runs show a decreasing trend in the pitch angle as the feathers get sheared due to differential rotation of the galaxy.

In the text

Fig. 11. Amplitude of the 2D Fourier transform A(m, p) of the logarithmic surface density (see Equation (10)) in the region 2–4 kpc, integrated over p. Different panels are for runs with ⟨Q⟩ = 1 and varying ℳc, evaluated at the same times as in Figure 6. The grey curves are simulation data and the orange ones are the empirical fits to the data (see text for details). We can see how the power shifts toward higher values of m with increasing ℳc.

In the text

Fig. 12. Azimuthal mode of the filaments (mfil) and various instabilities, shown with the ℳc ( ∼ vc/σ) of our galaxies. Black circular points with error bars show the mean and 16th to 84th percentile range of the mfil and corresponding physical spacing λfil for the snapshots shown in Figure 7. The dashed line is a linear fit to the median values of mfil (Equation (13)). Purple, red, and orange points with error bars represent the radial averages ⟨λ⟩ and standard deviations σλ (Equation (14) and Equation (15)) of the wavelength of the most unstable Toomre mode, the Jeans length and the Swing amplification length – see main text for details. These points are slightly offset from their overlapping ℳc values for visibility. The simulations show an increase in the median value of mfil and a broadening of the 16th to 84th percentile range with increasing ℳc. The median filament spacing agrees well with the Toomre length of the ℳc = 21, 29 and 36 runs, and with the swing length for ℳc = 14 run.

In the text

Fig. 13. Mass-weighted probability density function (PDF) of the log10 of the density normalised by the radially averaged mid-plane density, for the snapshots shown in Figure 7. The grey dashed line depicts the Jeans density (see text for details), at which we expect numerical effects to dominate. The left panel shows all the runs at times similar to Figure 7, and the right panel shows the time evolution of the Q1_ℳ36 case leading up to the time in the left panel. We see that all the runs in the left panel form high-density tails that represent feathers and that the densest end of the tail is reaching the limit imposed by the resolution of our simulations.

In the text

Fig. 14. Physical spacing of filaments as a function of the rotational Mach number, ℳc (∼vc/σ), in our simulations (black circles; cf., Section 3.3.2) compared with observations (yellow and blue diamonds). The yellow diamonds are JWST/MIRI observations (Meidt et al. 2023) and the blue diamonds are HST observations (La Vigne et al. 2006). We see that, except for NGC 1365 and NGC 4548, the filament spacing extracted from our simulations agrees within the uncertainty ranges of the ones reported in the observations.

In the text

	Fig. C.1. Same as Figure 2, but for the longer time evolution of the Q = 2 run. We see that there are faint structures at early times in the first two panels. The third panel shows a structure with somewhat more density contrast, visually confirming the slow growth of an instability.
In the text

	Fig. C.2. Same as Figure 5, but only for Q1_ℳ29 and Q2_ℳ29 runs, where the latter was evolved to a later time than shown in the main text. We see the oscillatory peaks in the Q = 2 run, but also a steady, yet slow growth over time. This growth is visibly smaller than in the Q1_ℳ29 run.
In the text

	Fig. D.1. Same as Figure 5, but for the Q1_ℳ29 run at three different resolutions. The solid line depicts the resolution used in the main text (Δx = 4.8 pc) and the other two lines are runs with resolutions reduced by a factor of two each.
In the text