Open Access
Issue
A&A
Volume 695, March 2025
Article Number A155
Number of page(s) 15
Section Extragalactic astronomy
DOI https://doi.org/10.1051/0004-6361/202453501
Published online 14 March 2025

© The Authors 2025

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.

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

Σ ( R ) = Σ ° exp [ R R d β exp ( α R R d ) ] , $$ \begin{aligned} \Sigma (R) = \Sigma _{\circ } \exp {\left[ -\frac{R}{R_{\rm d}} -\beta \exp \left(-\alpha \frac{R}{R_{\rm d}}\right) \right]}, \end{aligned} $$(1)

where R is the galactic (cylindrical) radius, Σ° is a scaling constant that will vary from run to run, Rd = 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

v rot = v c R R 2 + R c 2 , $$ \begin{aligned} { v}_{\rm rot} = {{ v}_{\rm c}}\frac{R}{\sqrt{R^{2} + R_{\rm c}^{2}}}, \end{aligned} $$(2)

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

P = ρ c s 2 , $$ \begin{aligned} P = \rho {c_{\rm s}} ^{2}, \end{aligned} $$(3)

where P is the gas pressure, ρ is the density, and cs 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,

ρ z ( ϕ g + ϕ dm ) + c s 2 ρ z = 0 , $$ \begin{aligned} \rho \frac{\partial }{\partial z}\left(\phi _{\rm g} + \phi _{\rm dm}\right) + {c_{\rm s}}^2 \frac{\partial \rho }{\partial z} = 0, \end{aligned} $$(4)

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

2 0 ρ ( R , z ) d z = Σ ( R ) , $$ \begin{aligned} 2 \int _0^\infty \rho (R,z) \, \mathrm{d}z = \Sigma (R), \end{aligned} $$(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 TCGM = 107 K and density. We transition from the cold galactic disc to the hot CGM when our discs reach a fixed transition density of ρ = 10−28 cm−3, 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,

Q = 1 R max 0 R max Q ( R ) d R = 1 R max 0 R max κ c s π G Σ d R , $$ \begin{aligned} \langle Q\rangle&= \frac{1}{R_{\rm max}}\int _0^{R_{\rm max}} Q(R) \, \mathrm{d}R = \frac{1}{R_{\rm max}}\int _0^{R_{\rm max}}\frac{\kappa {c_{\rm s}}}{\pi G \Sigma }\, \mathrm{d}R, \end{aligned} $$(6)

M c = v c c s , $$ \begin{aligned} {\mathcal{M} }_{\rm c}&= \frac{{ v}_c}{{c_{\rm s}}}, \end{aligned} $$(7)

where κ is the epicyclic frequency given by κ2 = 4Ω2 + 2RΩ(dΩ/dR), Ω = vrot/r is the angular frequency, and we take Rmax = 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 ≫ Rc where the rotation curve reaches its full speed; it relates to the rotational and thermodynamic properties of the galaxy and controls its thickness1. We show in Appendix B that these two dimensionless parameters, along with the auxiliary dimensionless ratios Rmax/Rd, Rc/Rd, α, β, 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 vc while keeping cs = 7 km s−1 fixed. The range of vc is chosen to reflect the range of rotation speeds observed in galaxies, with vc ∈ {100,150,200,250} km s−1 (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.

Table 1.

Parameter space of our simulations.

thumbnail 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 Lx = Ly = 20 kpc in the plane of the disc and Lz = 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 Rcyl = 6 kpc and half-height |Hcyl| = 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 trot = 2πR/vrot, 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.

thumbnail 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 = Rc = 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 trot. 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.

thumbnail 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 trot except for the ℳ36 run, which we show at t = 0.54 trot because it evolves faster than the other cases and collapses beyond our ability to resolve by 0.8 trot. 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.

thumbnail 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,

η = ln ( Σ Σ ) , $$ \begin{aligned} \eta = \ln \left(\frac{\Sigma }{\langle \Sigma \rangle } \right), \end{aligned} $$(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.

thumbnail 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 χ2 fit for the time-evolution of ση at times t ≥ 0.15 trot, which is approximately the time of the initial relaxation phase, to an empirically-motivated piecewise function of the form

σ η ( t ) = { A e ω t onset t t onset , A e ω t t > t onset , $$ \begin{aligned} \sigma _\eta (t) = {\left\{ \begin{array}{ll} A e^{\omega t_{\rm onset}}&t\le t_{\rm onset}, \\ A e^{\omega t}&t > t_{\rm onset}, \end{array}\right.} \end{aligned} $$(9)

where A sets the amplitude of ση during the steady phase and ω and tonset 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 tonset. We do not fit for the ℳc = 14 run since it shows a more complex behaviour compared to the other three runs.

Table 2.

Fit parameters for the evolution of the variance in the logarithmic surface density ση (Equation (8)).

For our best fits, we find t onset M c 1 $ 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 ω = 3.06 ± 0.08 t rot 1 = 0.49 ± 0.01 Ω $ \omega=3.06 \pm 0.08 \,{t_{\mathrm{rot}}}^{-1}=0.49 \pm 0.01\,\Omega $ for the ℳc = 36 case and the slowest growth rate ω = 2.41 ± 0.04 t rot 1 = 0.38 ± 0.01 Ω $ \omega=2.41 \pm 0.04 \,{t_{\mathrm{rot}}}^{-1}=0.38 \pm 0.01 \,\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 log10ση = −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.

thumbnail 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

thumbnail 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.

A ( m , p ) = u min u max 0 2 π η ( θ , u ) e i ( m θ + p u ) d θ d u , $$ \begin{aligned} A(m,p) = \int ^{u_{\rm max}} _{u_{\rm min}} \int ^{2\pi } _{0} \eta (\theta , u) e^{i(m\theta + pu)} \, \mathrm{d}\theta \, \mathrm{d}u, \end{aligned} $$(10)

where u = lnR, η is the logarithmic surface density (defined in Equation (8)), and (umin, umax) denote the radial bounds of our annuli. The term ei( + 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 A ¯ ( m , p ) = A ( m , p ) d m d p / d m d p $ \bar{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.

thumbnail 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−1(−m/p). Quantitatively, for each simulation we define mmax as the value of m at which the maximum of A(m, p) occurs. We then fit A(mmax, 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(mmax, p) versus p, and our best Gaussian fits to it, for the same four snapshots shown in Figure 7.

thumbnail 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−1(−mmax/μ) 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

Δ α = | d ( tan 1 ( m max / p ) ) d ( m max / p ) Δ ( m max / p ) | = | cos 2 α m max μ 2 σ | . $$ \begin{aligned} \Delta \alpha = \left|\frac{\mathrm{d}(\tan ^{-1} (-m_{\rm max}/p))}{\mathrm{d}(-m_{\rm max}/p)} \Delta (-m_{\rm max}/p)\right| = \left|\cos ^2 \alpha \frac{m_{\rm max}}{\mu ^{2}} \sigma \right|. \end{aligned} $$(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.

thumbnail 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 tonset, 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,

A ( m ) = A m m σ m 2 π exp [ ( log m μ m ) 2 2 σ m 2 ] , $$ \begin{aligned} A(m) = \frac{A_m}{m\sigma _{m} \sqrt{2\pi }} \exp \left[ \frac{(\log m - \mu _{m})^{2}}{2\sigma _{m}^{2}} \right], \end{aligned} $$(12)

where Am 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 mfil = 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.

thumbnail 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.

Table 3.

Fit parameters for the lognormal-like fit (see Equation (12)) on A(m) from our simulations at the times shown in Figure 11.

We show mfil 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 mfil, which is λfil = 2πRgal/mfil, where Rgal = 3 kpc is the mean radius of the annulus we use in our analysis. We find that the mfil varies approximately linearly with ℳc. A χ2 fit to this relationship of the form

m fil = a M c + b , $$ \begin{aligned} m_{\rm fil} = a{\mathcal{M} }_{\rm c} + b, \end{aligned} $$(13)

thumbnail 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 log10 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.

thumbnail 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 = (cs/4Δxmin)2(π/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Δxmin). Here, cs is the sound speed and Δxmin 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 = 2cs2/GΣ, λJeans = (πcs2/mid)1/2 and λSwing = λJeansXJ, respectively, where ρmid is the mid-plane density of the galaxy and XJ 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 XJ = 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

λ = 1 R out R in R in R out λ ( R , t = 0 ) d R , $$ \begin{aligned} \langle \lambda \rangle = \frac{1}{R_{\rm out}-R_{\rm in}} \int ^{R_{\rm out}} _{R_{\rm in}} \lambda (R, t=0)\, \mathrm{d}R, \end{aligned} $$(14)

where Rin = 2 kpc and Rout = 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

σ λ = { 1 R out R in R in R out [ λ ( R , t = 0 ) λ ] 2 d R } 1 / 2 . $$ \begin{aligned} \sigma _\lambda = \left\{ \frac{1}{R_{\rm out}-R_{\rm in}} \int ^{R_{\rm out}} _{R_{\rm in}} \left[\lambda (R, t=0)-\left\langle \lambda \right\rangle \right]^2\, \mathrm{d}R \right\} ^{1/2}. \end{aligned} $$(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.

thumbnail 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 mfil at all radii, which we extract by carrying out a simple χ2 fit for mfil to the radius-dependent values of λfil reported in Mandowara et al.’s Table 4; we then compute λfil = 2πR/mfil at R = 3 kpc. For ℳc in both sets of observations, we assume a fixed cs = 7 km s−1 for all galaxies, characteristic of the warm neutral medium (McClure-Griffiths et al. 2023), and take our circular velocities vc 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 vc, assuming no uncertainty in cs. 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 vc = 120 ± 20.6 km s−1, 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 (Mstellar), are tabulated in Table 4. The SFR and the Mstellar 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.

Table 4.

Physical properties of the galaxies we use for comparison with our simulations.

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 M c 1 $ {\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 M c 1 $ {\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


1

Our ℳc parameter is similar to the quantity vc/σ used in studies of galaxy kinematics to distinguish rotation-dominated from dispersion-dominated galaxies. Note that these are not fully identical since in observations σ contains both thermal and non-thermal contributions to the velocity dispersion. By contrast, for ℳc our dispersion is purely thermal.

2

Formally the expression that Wang et al. (2010) provide is for a purely exponential disc, without the central flattening we have introduced. However, since our disc differs in mass from a pure exponential disc by only ≈5%, we expect the difference in potential to be small.

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

  1. Arora, R., Federrath, C., Banerjee, R., & Körtgen, B. 2024, A&A, 687, A276 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  2. Binney, J., & Tremaine, S. 1987, Galactic dynamics (Princeton, N.J.: Princeton University Press) [Google Scholar]
  3. Chepurnov, A., Lazarian, A., Stanimirović, S., Heiles, C., & Peek, J. E. G. 2010, ApJ, 714, 1398 [NASA ADS] [CrossRef] [Google Scholar]
  4. Dobbs, C. L., & Bonnell, I. A. 2006, MNRAS, 367, 873 [NASA ADS] [CrossRef] [Google Scholar]
  5. Dubey, A., Fisher, R., Graziani, C., et al. 2008, in Numerical Modeling of Space Plasma Flows, eds. N. V. Pogorelov, E. Audit, & G. P. Zank, ASP Conf. Ser., 385, 145 [NASA ADS] [Google Scholar]
  6. Elmegreen, D. M. 1980, ApJ, 242, 528 [Google Scholar]
  7. Elmegreen, B. G., & Elmegreen, D. M. 2019, ApJS, 245, 14 [NASA ADS] [CrossRef] [Google Scholar]
  8. Elmegreen, D. M., Elmegreen, B. G., Kaufman, M., et al. 2006, ApJ, 642, 158 [NASA ADS] [CrossRef] [Google Scholar]
  9. Elmegreen, D. M., Elmegreen, B. G., Erroz-Ferrer, S., et al. 2014, ApJ, 780, 32 [Google Scholar]
  10. Elmegreen, B. G., Elmegreen, D. M., & Efremov, Y. N. 2018, ApJ, 863, 59 [NASA ADS] [CrossRef] [Google Scholar]
  11. Federrath, C., Roman-Duval, J., Klessen, R. S., Schmidt, W., & Mac Low, M. 2010, A&A, 512, A81 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  12. Federrath, C., Sur, S., Schleicher, D. R. G., Banerjee, R., & Klessen, R. S. 2011, ApJ, 731, 62 [NASA ADS] [CrossRef] [Google Scholar]
  13. Federrath, C., Roman-Duval, J., Klessen, R. S., Schmidt, W., & Mac Low, M. M. 2022, Astrophysics Source Code Library [record ascl:2204.001] [Google Scholar]
  14. Fryxell, B., Olson, K., Ricker, P., et al. 2000, ApJS, 131, 273 [Google Scholar]
  15. Fuchs, B. 2001, A&A, 368, 107 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  16. Griv, E., & Gedalin, M. 2012, MNRAS, 422, 600 [NASA ADS] [CrossRef] [Google Scholar]
  17. Griv, E., & Wang, H.-H. 2014, New Astron., 30, 8 [NASA ADS] [CrossRef] [Google Scholar]
  18. Jeans, J. H. 1902, R. Soc. Lond. Philos. Trans. Ser. A, 199, 1 [CrossRef] [Google Scholar]
  19. Jeffreson, S. M. R., Kruijssen, J. M. D., Keller, B. W., Chevance, M., & Glover, S. C. O. 2020, MNRAS, 498, 385 [Google Scholar]
  20. Khullar, S., Matzner, C. D., Murray, N., et al. 2024, ApJ, 973, 40 [NASA ADS] [CrossRef] [Google Scholar]
  21. Kim, W.-T., & Ostriker, E. C. 2002, ApJ, 570, 132 [Google Scholar]
  22. Kim, W.-T., & Ostriker, E. C. 2006, ApJ, 646, 213 [NASA ADS] [CrossRef] [Google Scholar]
  23. Kim, W.-T., Kim, Y., & Kim, J.-G. 2014, ApJ, 789, 68 [NASA ADS] [CrossRef] [Google Scholar]
  24. Kim, Y., Kim, W.-T., & Elmegreen, B. G. 2015, ApJ, 809, 33 [NASA ADS] [CrossRef] [Google Scholar]
  25. Kim, J.-H., Agertz, O., Teyssier, R., et al. 2016, ApJ, 833, 202 [NASA ADS] [CrossRef] [Google Scholar]
  26. Kim, W.-T., Kim, C.-G., & Ostriker, E. C. 2020, ApJ, 898, 35 [Google Scholar]
  27. Körtgen, B., Banerjee, R., Pudritz, R. E., & Schmidt, W. 2019, MNRAS, 489, 5004 [CrossRef] [Google Scholar]
  28. Kuhn, M. A., Benjamin, R. A., Zucker, C., et al. 2021, A&A, 651, L10 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  29. Lang, P., Meidt, S. E., Rosolowsky, E., et al. 2020, ApJ, 897, 122 [CrossRef] [Google Scholar]
  30. La Vigne, M. A., Vogel, S. N., & Ostriker, E. C. 2006, ApJ, 650, 818 [Google Scholar]
  31. Lee, W.-K. 2014, ApJ, 792, 122 [NASA ADS] [CrossRef] [Google Scholar]
  32. Lee, W.-K., & Shu, F. H. 2012, ApJ, 756, 45 [NASA ADS] [CrossRef] [Google Scholar]
  33. Leroy, A. K., Sandstrom, K. M., Lang, D., et al. 2019, ApJS, 244, 24 [Google Scholar]
  34. Leroy, A. K., Schinnerer, E., Hughes, A., et al. 2021, ApJS, 257, 43 [NASA ADS] [CrossRef] [Google Scholar]
  35. Levy, R. C., Bolatto, A. D., Teuben, P., et al. 2018, ApJ, 860, 92 [NASA ADS] [CrossRef] [Google Scholar]
  36. Mandowara, Y., Sormani, M. C., Sobacchi, E., & Klessen, R. S. 2022, MNRAS, 513, 5052 [NASA ADS] [CrossRef] [Google Scholar]
  37. McClure-Griffiths, N. M., Stanimirović, S., & Rybarczyk, D. R. 2023, ARA&A, 61, 19 [NASA ADS] [CrossRef] [Google Scholar]
  38. McGaugh, S. S., & Schombert, J. M. 2015, ApJ, 802, 18 [NASA ADS] [CrossRef] [Google Scholar]
  39. Meidt, S. E. 2022, ApJ, 937, 88 [NASA ADS] [CrossRef] [Google Scholar]
  40. Meidt, S. E., & van der Wel, A. 2024, ApJ, 966, 62 [NASA ADS] [CrossRef] [Google Scholar]
  41. Meidt, S. E., Leroy, A. K., Rosolowsky, E., et al. 2018, ApJ, 854, 100 [NASA ADS] [CrossRef] [Google Scholar]
  42. Meidt, S. E., Rosolowsky, E., Sun, J., et al. 2023, ApJ, 944, L18 [CrossRef] [Google Scholar]
  43. Pantaleoni González, M., Maíz Apellániz, J., Barbá, R. H., & Reed, B. C. 2021, MNRAS, 504, 2968 [CrossRef] [Google Scholar]
  44. Parker, E. N. 1966, ApJ, 145, 811 [NASA ADS] [CrossRef] [Google Scholar]
  45. Puerari, I., Elmegreen, B. G., & Block, D. L. 2014, AJ, 148, 133 [NASA ADS] [CrossRef] [Google Scholar]
  46. Robinson, H., & Wadsley, J. 2024, MNRAS, 534, 1420 [NASA ADS] [CrossRef] [Google Scholar]
  47. Romeo, A. B., & Falstad, N. 2013, MNRAS, 433, 1389 [NASA ADS] [CrossRef] [Google Scholar]
  48. Romeo, A. B., & Wiegert, J. 2011, MNRAS, 416, 1191 [NASA ADS] [CrossRef] [Google Scholar]
  49. Savchenko, S. S. 2012, Astrophys. Bull., 67, 310 [NASA ADS] [CrossRef] [Google Scholar]
  50. Schinnerer, E., Meidt, S. E., Colombo, D., et al. 2017, ApJ, 836, 62 [NASA ADS] [CrossRef] [Google Scholar]
  51. Shetty, R., & Ostriker, E. C. 2006, ApJ, 647, 997 [Google Scholar]
  52. Sormani, M. C., Sobacchi, E., Shore, S. N., Treß, R. G., & Klessen, R. S. 2017, MNRAS, 471, 2932 [NASA ADS] [CrossRef] [Google Scholar]
  53. Thilker, D. A., Lee, J. C., Deger, S., et al. 2023, ApJ, 944, L13 [NASA ADS] [CrossRef] [Google Scholar]
  54. Toomre, A. 1964, ApJ, 139, 1217 [Google Scholar]
  55. Tress, R. G., Sormani, M. C., Glover, S. C. O., et al. 2020, MNRAS, 499, 4455 [NASA ADS] [CrossRef] [Google Scholar]
  56. Truelove, J. K., Klein, R. I., McKee, C. F., et al. 1997, ApJ, 489, L179 [CrossRef] [Google Scholar]
  57. Veena, V. S., Schilke, P., Sánchez-Monge, Á., et al. 2021, ApJ, 921, L42 [NASA ADS] [CrossRef] [Google Scholar]
  58. Wada, K., & Koda, J. 2004, MNRAS, 349, 270 [Google Scholar]
  59. Wang, H.-H., Klessen, R. S., Dullemond, C. P., van den Bosch, F. C., & Fuchs, B. 2010, MNRAS, 407, 705 [CrossRef] [Google Scholar]
  60. Williams, T. G., Sun, J., Barnes, A. T., et al. 2022, ApJ, 941, L27 [NASA ADS] [CrossRef] [Google Scholar]
  61. Zhao, B., Pudritz, R. E., Pillsworth, R., Robinson, H., & Wadsley, J. 2024, ApJ, submitted [arXiv:2405.18474] [Google Scholar]
  62. Zucker, C., Battersby, C., & Goodman, A. 2015, ApJ, 815, 23 [NASA ADS] [CrossRef] [Google Scholar]

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

ρ R ( ϕ dm + ϕ g ) + c s 2 ρ R ρ Ω 2 R = 0 , $$ \begin{aligned} \rho \frac{\partial }{\partial R}\left(\phi _{\rm dm} + \phi _{\rm g}\right) + {c_{\rm s}}^2 \frac{\partial \rho }{\partial R} - \rho \Omega ^2 R = 0, \end{aligned} $$(A.1)

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

ϕ g = 4 π G Σ ° R d y 2 [ I 0 ( y ) K 0 ( y ) I 1 ( y ) K 1 ( y ) ] , $$ \begin{aligned} \phi _{\rm g} = 4\pi G \Sigma _\circ {R_{\rm d}} y^2 \left[I_0(y) K_0(y) - I_1(y) K_1(y)\right], \end{aligned} $$(A.2)

where y = R/2Rd and I0, K0, I1, and K1 are the modified Bessel functions of the first and second kind, of order zero and one.2 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/Rd). 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 ≪Rd; 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 Ω = vrot/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 = Ω2R. For this purpose we adopt the flattened centrifugal potential form suggested by Binney & Tremaine (1987),

ψ = v c 2 2 ln { 1 R c 2 [ R 2 + R c 2 + ( z q ) 2 ] } , $$ \begin{aligned} \psi = \frac{{v_{\rm c}} ^{2}}{2} \ln \left\{ \frac{1}{R_{\rm c} ^{2}} \left[ R^2 + R_{\rm c}^2 + \left( \frac{z}{q} \right)^2 \right] \right\} , \end{aligned} $$(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 = Ω2R. Following Binney & Tremaine’s estimates for realistic galactic potentials, we adopt q = 0.7. With these assumptions, our equation of force balance reduces to

ρ R ( ϕ dm + ϕ g c s 2 R R d ψ ) = 0 , $$ \begin{aligned} \rho \frac{\partial }{\partial R}\left(\phi _{\rm dm} + \phi _{\rm g} - {c_{\rm s}}^2 \frac{R}{R_d} - \psi \right) = 0, \end{aligned} $$(A.4)

which has the immediate solution

ϕ dm = ψ ϕ g + c s 2 R R d . $$ \begin{aligned} \phi _{\rm dm} = \psi - \phi _{\rm g} + {c_{\rm s}}^2 \frac{R}{{R_{\rm d}}}. \end{aligned} $$(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 – Rc/Rd, Rmax/Rd, q, α, and β – fully determine the evolution of the system. To do so, we non-dimensionalise the fluid variables via the transformations

r = x R d , $$ \begin{aligned} \boldsymbol{r}&= \boldsymbol{x} {R_{\rm d}}, \end{aligned} $$(B.1)

t = τ R d c s , $$ \begin{aligned} t&= \tau \frac{{R_{\rm d}}}{{c_{\rm s}}}, \end{aligned} $$(B.2)

v = u c s , $$ \begin{aligned} \boldsymbol{v}&= \boldsymbol{u}{c_{\rm s}}, \end{aligned} $$(B.3)

ρ = a Σ ° R d , $$ \begin{aligned} \rho&= a\frac{\Sigma _{\circ }}{{R_{\rm d}}}, \end{aligned}$$(B.4)

ϕ = ϕ c s 2 . $$ \begin{aligned} {\phi }&= \frac{\phi }{{c_{\rm s}}^{2}}. \end{aligned} $$(B.5)

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

a τ + 1 M c x ( a u ) = 0 , $$ \begin{aligned} \frac{\partial a}{\partial \tau } + \frac{1}{{\mathcal{M} }_{\rm c}} \boldsymbol{\nabla _{x}}(a\boldsymbol{u}) =&~0, \end{aligned} $$(B.6)

M c u τ + ( u · x ) · u = x a a x ( ϕ g + ϕ dm ) , $$ \begin{aligned} {\mathcal{M} }_{\rm c}\frac{\partial \boldsymbol{u} }{\partial \tau } + (\boldsymbol{u} \cdot \nabla _{x}) \cdot \boldsymbol{u} =&-\frac{\nabla _{x} a}{a} - \nabla _{x} \left({\phi }_{\rm g} + {\phi }_{\rm dm} \right), \end{aligned} $$(B.7)

x 2 ϕ = ξ M c Q a , $$ \begin{aligned} \nabla ^{2}_{x} {\phi } =&\xi \frac{{\mathcal{M} }_{\rm c}}{\langle Q\rangle } a , \end{aligned} $$(B.8)

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

ξ = 4 2 X max 0 X max 1 exp [ X β exp ( α X ) ] X 2 + 2 X c 2 X 2 + X c 2 d X 19.43 . $$ \begin{aligned} \xi&= \frac{4 \sqrt{2}}{X_{\rm max}} \int _0^{X_{\rm max}} \frac{1}{\exp \left[-X-\beta \exp \left(-\alpha X\right)\right]} \frac{\sqrt{X^2+2X_{\rm c}^2}}{X^2 + X_{\rm c}^2} \, dX \nonumber \\&\approx 19.43. \end{aligned} $$(B.9)

In the expression above, we have defined X = R/Rd as the dimensionless cylindrical radius, Xc and Xmax 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 ϕ $ \tilde{\phi} $ using the same non-dimensionalisation procedure. Substituting into equation A.5 we have

ϕ dm = ψ ϕ g + X , $$ \begin{aligned} {\phi }_{\rm dm}&= {\psi } - {\phi }_{\rm g} + X, \end{aligned} $$(B.10)

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

ϕ g = M c Q ξ X 2 [ I 0 ( X 2 ) K 0 ( X 2 ) I 1 ( X 2 ) K 1 ( X 2 ) ] , $$ \begin{aligned} {\phi }_{\rm g}&= \frac{{\mathcal{M} }_{\rm c}}{\langle Q\rangle \xi } X^2 \left[ 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)\right], \end{aligned} $$(B.12)

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

Appendix C: Toomre ⟨Q⟩ = 2 run

thumbnail 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 ∼trot timescales as the ⟨Q⟩ = 1 runs presented in the main body. We continue this run to time t = 2.5 trot, 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 trot 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 mfil ≃ 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.

thumbnail 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 trot, 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 ωtrot = 0.31, which is ≈8 times smaller than Q1_ℳ29 case.

Appendix D: Resolution study

thumbnail 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

Table 1.

Parameter space of our simulations.

Table 2.

Fit parameters for the evolution of the variance in the logarithmic surface density ση (Equation (8)).

Table 3.

Fit parameters for the lognormal-like fit (see Equation (12)) on A(m) from our simulations at the times shown in Figure 11.

Table 4.

Physical properties of the galaxies we use for comparison with our simulations.

All Figures

thumbnail 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
thumbnail 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
thumbnail 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
thumbnail 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
thumbnail 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
thumbnail 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
thumbnail 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
thumbnail 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
thumbnail 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
thumbnail 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
thumbnail 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
thumbnail 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
thumbnail 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
thumbnail 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
thumbnail 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
thumbnail 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
thumbnail 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

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.