Large fluctuations and cusp-core transformation in dwarf galaxies

Main

Astronomical observations, interpreted in terms of general relativity and Newtonian dynamics, imply that our universe is dominated by a dark sector. In the standard $\Lambda$CDM scenario, this consists of dark energy, represented by a cosmological constant $\Lambda$, and cold dark matter (CDM) [1]. Phenomenologically, the dark energy is particularly simple; its density is homogeneous in space and time. The description of the dynamics of the dark matter is more involved, as it can collapse into nonlinear haloes with average overdensities of the order of a hundred relative to the cosmological background density. Nevertheless, the structure of such haloes is now well established through numerical simulations; in particular, it has long been known that their central density diverges as one zooms in towards the centre down to the resolution limit [2]. However, soon after the results of the first simulations displaying central density cusps became available, it was pointed out that the dynamics of dark matter dominated galaxies, as inferred from observations, may imply less steep central density profiles [3]. Indeed, although there is a diversity in galaxy rotation curves, perhaps partly born of systematic effects related to non-circular motions, at least some dwarfs appear to display non-singular central density profiles, and generally contain less dark matter in the central region than is expected from dissipationless cosmological simulations [4, 5].

In other words, the dark matter structures invoked in part to explain the outer dynamics of galaxies may not fit those of the inner parts. Solutions to this longstanding ’cusp/core’ problem have invoked variations on the dark matter model, in attempts to dynamically ‘heat’ the centres of haloes and decrease their concentrations. These include pre-heated warm dark matter [6], thermally conducting self-interacting dark matter [7], and zero-point quantum fluctuations in fuzzy dark matter [8]. But cold dark matter can itself be heated by gaining random energy from massive gaseous or stellar clumps via dynamical friction [9], or through sudden changes in the potential triggered by rapid feedback (e.g. starburst) driven gas causing outflows [10]. In the latter case, numerous outflows and inflows, from and towards the central region, are most effective in transforming the density cusp into a core [11]. The resulting irreversible dynamics resemble those stemming from massive clumps imparting random energy to the dark matter as they move through it [12, 13, 14, 15, 16]. Observationally, however, it may be possible to distinguish between the two mechanisms, through star formation histories and the associated energetics of the feedback that drives large scale gas motions. Thus, while both mechanisms may act in tandem [17], with dynamical friction from massive clumps being efficient in large galaxies, simultaneously forming halo cores and stellar bulges [18], feedback-driven fluctuations may be the dominant mechanism in dwarfs [19, 20]. But modelling the physical processes that may lead to the cusp/core transformation in this context is far more involved than simulating the dark matter, and the conclusions are consequently much less robust or secure.

In standard cosmology, galaxies form as gas settles into the centres of CDM haloes, becoming dynamically important and forming stars [21]. Beyond an early period of weak force mediated thermal equilibrium, which sets the scale of the smallest objects, the formation of CDM structures may be described in terms of Newtonian gravitational dynamics on scales significantly smaller than the cosmological horizon, relevant to galaxies. The baryons, on the other hand, can also interact electromagnetically and through the strong nuclear force. They are thus subject to complex hydrodynamic forces and thermal and feedback processes in a highly inhomogeneous multiphase medium. These may be modelled through large scale numerical simulations, but sub-grid models are required to describe star formation and stellar/AGN feedback, and their mode of energy injection and transfer [22].

Such simulations generally reproduce the observed low baryon to dark matter content in dwarf galaxies, resulting from a shallow potential well. This property renders dwarf galaxies dark matter dominated even in the central region, and consequently makes for ideal laboratories for testing predicted dark matter density profiles against observations [23, 24, 25, 26, 27, 28, 29]. However, whether cusps are transformed into cores in a given hydrodynamic simulation depends on the details of the implemented subgrid modelling. It has been long suspected that the key ingredient involved is the star formation threshold—with a high threshold leading to conditions favourable for core formation and a low threshold being unfavourable [13, 14, 30, 31, 32]. But it has recently been realised that additional details of the modelling of the subgrid physics of feedback and energy transfer may also be important [33].

The key dynamical condition for core formation in CDM-based simulations with feedback appears to be the presence of sustained fluctuations in the gravitational potential of the gas, which couple to the dynamics of the dark matter. Ref. [15] proposed that such a situation is amenable to analysis in terms of first dynamical principles if the fluctuations form a stochastic process described by a power spectrum. In this case, the dynamical ‘heating’ of the dark matter particles can be given a precise meaning, through the random energy they gain from stochastic fluctuations. The effect is predicted to depend primarily on two parameters: the gas density and the amplitude of the spectrum of density fluctuations. It thereby bypasses the complex details of the baryonic physics of galaxy formation. The assumptions and consequences of this theoretical model were tested against the results of a numerical simulation of a model dwarf galaxy in ref. [16].

Although gravitational potential fluctuations in the gas cannot be observed, the model maps these to the density field and determines the dynamical response of halo particles. It thus provides a pathway for directly testing whether the dynamics driven by these may indeed lead to the cusp-core transformation. In this context, observations of local galaxies should provide a lower limit of the most important parameters, as the gas density should on average be smaller than at earlier times (due to outflows, and star formation), and likewise with the amplitude of the fluctuations (as feedback is presumably stronger at earlier times). Measuring these quantities in local galaxies should thus provide an upper limit of the relaxation time for core formation. And with this comes a direct test of the proposal of cusp/core transformation through potential fluctuations, independently of the complexities involved in numerical and subgrid implementations inherent to simulations of the process. This is what we set out to do here.

Power spectra

The bulk of the gas mass in dwarf irregular galaxies is in atomic form [34]. In this medium, density fluctuations arise from the presence of gas clouds, feedback-driven voids, shells and turbulent structures. We choose a sample of such galaxies from the LITTLE THINGS HI galaxy survey of local ($< 10$ Mpc) dwarfs, which has sufficient resolution (110 pc at mean distance 3.7 Mpc), to capture the aforementioned structures [35]. The survey includes gas-rich galaxies with absolute visual magnitudes between -9 and -18. They are drawn from earlier databases where morphological inclusion criteria were also used [36, 37]. The total sample includes 37 dwarf irregulars and 4 blue compact dwarfs.

We produced power spectra of relative surface density fluctuations $\delta\Sigma/\Sigma= \Sigma(\mathbf{R})/\Sigma(R) - 1$ (as detailed in § 8), where $\Sigma(\mathbf{R})$ is the surface density of the gas at position $\mathbf{R}$ on the image of the galaxy and $\Sigma(R)$ refers to averages over rings at radius $\lvert\mathbf{R}\rvert$. This is a measure of the level of density perturbations that is independent of the absolute value of the average density. In addition to being theoretically relevant (as we will see below), it is a dimensionless measure facilitating comparison between different systems.

The resulting power spectra for the whole sample of galaxies are shown in Figure 1, where we also plot the time average of the corresponding spectra in the simulation analysed in ref. [16] (referred to as ISIM), which employs a feedback recipe with high star formation threshold (300 atoms cm$^{-3}$ and and efficiency $\epsilon_{\rm eff} = 0.1$ per local free fall time, allowing for bursty episodic feedback. The gas was simulated with maximum resolution of about 4 pc, to capture momentum driving from individual supernovae events. The recipe led to results that compared well with various properties of observed dwarf irregular galaxies, including rotation curves, stellar surface density, metallicity, gas and stellar velocity dispersion and star formation history [24, 40], despite the isolated setting of the simulation. On the other hand, this property renders the results ideal for testing the theoretical model, in order to capture and isolate the central physical and dynamical features of the core formation process in a controlled setting.

We have also examined the power spectra of gas embedded in simulated cosmological haloes, as illustrated in Figure 1. The FIRE-2 cosmological simulation [39] employs a feedback recipe with a high star formation threshold, broadly similar to ISIM, which allows for core formation [41]. For haloes that formed cores in FIRE, the resulting power spectra were likewise consistent with the data within the uncertainties. In contrast, results for the gas embedded in haloes from the Auriga simulation [38], were quite different. There, a low star formation threshold is adopted, and feedback is regulated by a long star formation timescale and a stiff equation of state. This results in less episodic star formation and, crucially, in the absence of large density fluctuations. Consequently, the spectra display substantial quantitative and qualitative differences, reflected in their normalization and tilt, respectively, as can be readily seen from the example spectrum included in Figure 1. We present further details regarding the analysis of these cosmological simulations in § 7.2.

In ref [16] it was shown that the gas density fluctuations in ISIM, as reflected in the power spectrum, can explain the cusp/core inferred transformation in the context of our model. Figure 1 shows that the normalizations and slopes of the spectra of the observed galaxies are distributed around the line corresponding to the simulation. This appears remarkable, given that the average in the simulation is taken with a starting time of 4 Gyr, when a significant core starts forming, to a Hubble time, while the power spectra inferred from the observed sample correspond to local galaxies. (Note that these are the measured spectra in physical units, with no scaling introduced.) This may suggest that the density and potential fluctuations may indeed be approximately represented by quasi-stationary processes, as is assumed in the model and approximately realized in the simulation.

As the density profiles of LITTLE THINGS galaxies are generally consistent with having halo cores rather than cusps (§ 10), the matching of their gas spectra with those of gas embedded in core-forming haloes (in ISIM and FIRE)—in contrast to the case of Auriga, where no cores were found to form [32]—naturally leads to the following question: have the fluctuations characterized by the inferred spectra led to a cusp/core transformation in the observed galaxies, in accordance with expectations from our model? We test for this in the coming section.

The relaxation time of core formation

The physical picture underlying our model is as follows. As the gas settles into the centre of a dark matter halo, a critical density threshold for star formation is reached. If that threshold is sufficiently large, star formation leads to a prolonged episodic feedback process, accompanied by repeated large scale outflows and inflows confined within the central region of the dark matter halo (as opposed to escaping the potential well). Within that region, the large scale flows carry smaller scale density fluctuations, which may be represented by a three dimensional (3D) power spectrum $P_{3D}(k) \equiv P(k)$. In analogy with compressible turbulence and observations of the interstellar medium, reference [15] assumed that $P(k)$ is a power law between maximal and minimal wavenumbers $k_x$ and $k_m$ respectively. As we saw, a power law within a certain interval is also an adequate description of the two dimensional (2D) power spectra of the HI gas in LITTLE THINGS galaxies.

A summary of the main features of the dynamical model developed in [15] is given in § 5. In its context, gaseous density perturbations are readily transformed into potential fluctuations through the Poisson equation. The resulting fluctuating force imparts random velocities to the dark matter particles moving through the gas, into higher energy orbits away from the centre. This reduces the central density, transforming the cusp into a core. The process takes place over a relaxation time, over which the random velocities gained by the particles become of the order of their orbital speed. This timescale is given by

Here, $n > 0$ is the index of the power law spectrum, $P(k) \propto k^{-n}$; $P(k_m)$ is the value of the 3D power spectrum at the minimal wavenumber $k_m$; $\rho_0$ is the characteristic 3D gas density, evaluated at the projected half-mass radius (§ 6); $v_p$ is the typical velocity of a dark matter particle affected by the stochastic fluctuation field; and $v_r$ is the velocity of the large scale gaseous flows carrying the density fluctuations. This equation is a generalization of Chandrasekhar’s two body relaxation formula when $n \ne0$. The latter corresponds to a white noise power spectrum, $P = 1/N$ (with $N$ the particle number density) and Coulomb logarithm $\Lambda= k_x/k_m$ [42]. For $n \ne0$, equation (1) replaces the usual formula, with $n/P(k_m)$ standing for $1/(P\ln\Lambda)$. It is customary in applying Chandrasekhar’s formula to realistic (inhomogeneous) systems, to assume a local density, possibly averaged within a certain radius (e.g., the half mass radius), while the maximum impact parameter is assumed to correspond to the entire size of the system. In the present case, while the characteristic density is defined in terms of the half mass radius, the spectra are evaluated at the 90% radius, in order to capture the power at the largest scales associated with $k_m$ and ensure convergence of the spectra (§ 8). Note that the dependence of $t_{\mathrm{relax}}$ on the largest fluctuations—when the spectrum displays a well-defined power-law region between $k_m$ and $k_x$—renders its estimation robust against small-scale noise and limited resolution.

As the index $n$ and normalization $P(k_m)$ refer to 3D spectra, to connect with observations we calibrate the 2D spectra into 3D ones using simulations and a projection model (as described in § 7). We find that, for the parameters of the galaxies at hand, the index $n_{3D} \equiv n \simeq1.1\,n_{2D}$. The maximal fluctuation scales (corresponding to $k_m$) are set through the maximal scales of HI cloud structures, estimated for each galaxy (§ 8). This has the advantage of being independent from the transformation of the power spectrum from 2D to 3D and insensitive to large scale deformations in the gas patterns (e.g., spiral structures).

For a subset of the galaxies in LITTLE THINGS, high resolution mass models are available [43], including fits to the halo mass distributions using isothermal profiles with quasi-constant density inside core radius $r_c$. We use these to estimate the characteristic velocities entering into equation (1), by fixing $v_p$ to the circular velocity at $r_c$. In this way $t_{\mathrm{relax}}(r_c)$ estimates the timescale over which a core of the observed size $r_c$ can be formed through the effect of potential fluctuations. The assumption here is that the core radius reflects the scale inside which an initial inner density cusp is significantly altered. For, although the existence of a precisely zero inner logarithmic slope for the density profile is generally difficult to confirm observationally (§ 10), for most galaxies in the sample the cored profile was found to be strongly favoured by the data, while a steep inner cusp is entirely inconsistent with it [43]. Note that because the circular speed increases with radius (as does the velocity dispersion in an initial cusp), the relaxation time satisfies $t_{\mathrm{relax}}(r \leq r_c) \leq t_{\mathrm{relax}}(r_c)$. Consequently, if gas fluctuations significantly affect particles at $r_c$, they will also affect all particles within $r_c$. The criterion is tested in Figure 5.

Following refs. [15] and [16], we estimate the velocity of large-scale gas flows—assumed to carry the density fluctuations—as $v_r = (d/2)/t_D(d/2)$, where $t_D(d/2) = 1/\sqrt{G\bar{\rho}_{\mathrm{DM}}(d/2)}$ (with $\bar{\rho}_{\mathrm{DM}}(d/2)$ being the average dark matter density inside $d/2$). This form of $v_r$ implies a characteristic diameter $d$ within which feedback-driven fountain flows are contained, maintaining a nearly constant total gas mass inside $d/2$ even as the inner mass fluctuates strongly (cf. Figure 2 of ref. [16]). Because feedback-driven fluctuations diminish rapidly beyond $d/2$ by definition, previous studies assumed this boundary to be the scale at which the power spectra are evaluated [15, 16]. We adopt the same assumption here; for our galaxy sample, we equate $d/2$ to the convergence scale of the power spectra, $R_{90}$. Evaluating $v_r$ at $R_{90}$ is conservative in the sense that $R_{90}$ typically lies in the flat region of the rotation curve, making $v_r$ insensitive to exact radial positioning. Conversely, evaluating it at smaller radii yields shorter relaxation times $t_{\mathrm{relax}}$ (as noted above in the case of $v_p$).

Images showing the density contrasts $\Sigma(\mathbf{R})/\bar{\Sigma}(R)$ in maps of the galaxy subsample for which dynamical modelling is available, are shown in Figure 2, along with the relaxation times. As can be seen, the latter are of the order of a Hubble time or less for most galaxies. This is explicitly illustrated in Figure 3, where the relaxation times are seen to cluster around a value somewhat smaller than a Hubble time. The conspicuous exception is DDO 47, which has an anomalously large $t_{\mathrm{relax}}$. But as discussed in § 10, it may also have a steep inner cusp.

We do not, of course, know the dynamical history of these galaxies; what we measure is the relaxation time inferred from current parameters. Strictly speaking, it is related to future dynamical evolution. However, again, unless the gas fluctuations can be characterized by a perfectly stationary stochastic process, all crucial parameters in equation (1) may collude towards smaller relaxation times, and thus more efficient cusp/core transformation in the past: the fluctuation levels (as measured by $P(k_{\mathrm{m}})$) should have been larger; the gas density $\rho_0$ should also have been larger; furthermore, the velocities should also have been smaller if the galaxy masses increased through past mergers. In this sense, the relaxation times calculated here are to be taken as upper limits. Alternatively, one can interpret the clustering of the results around a Hubble time as follows: suppose the process of stochastic fluctuations is indeed strictly stationary as the model assumes, then a core of the order of that observed is expected to form over a Hubble time, and consequently have an associated relaxation time of that order (as is the case with the simulation, where the relaxation time at the core radius at the end of the simulation was calculated to be of order a Hubble time [16]).

However, we find a clear correlation between the relaxation time and the virial mass of the system, as shown in Figure 4. Such scaling is expected, as the relaxation time increases as the cube of the velocity. This dependence arises from a $v_p^2$ term (due to larger particle kinetic energies in deeper potential wells) and a $v_r$ term (reflecting the fact that gravitational encounters between gas fluctuations and dark matter particles are weaker at larger relative velocities). Assuming a Tully-Fisher or Faber-Jackson type relation ($v \sim M^{1/4}$), the relaxation time scales with mass as $v^3 \sim M^{3/4}$.

For the points displayed with error bars in Figure 4, the increase in $t_{\mathrm{relax}}$ with mass is further reinforced by the gas density scaling as $\rho_0 \sim1/r_0$, where $\rho_0$ is given by equation (11) and displayed in Figure 7. This arises from the well-known mass-radius relation for HI gas and leads to a lower density in more massive systems. However, the HI density scaling may not universally hold; the total gas density may be nearly constant, or even increase with galaxy mass, when accounting for the increasing fractions of molecular gas in more massive galaxies. This may arise from an ascending molecular-to-atomic gas ratio [44, 45], potentially enhanced by undetected CO-dark molecular gas [46, 47]. Thus, for more massive galaxies in particular, our fiducial estimates of the relaxation times (displayed in Figure 11 and Figure 3) should be taken as upper limits if a significant molecular component is present. For comparison, we also show the relaxation times evaluated while assuming a constant $\rho_0 = 5 \times10^6\,M_\odot/\mathrm{kpc}^3$ in Figure 4, which are indeed generally shorter at higher masses.

On the other hand, the potential limitation of forming cores at higher mass may be compounded when accounting for the formation of dominant central stellar systems within haloes evolving in a cosmological setting [19, 48, 49, 50]. A related issue concerns whether cusps may reform during mergers, as cusps can represent equilibria of the merging process [51, 52]; massive haloes that grow through cusp-reforming processes may not maintain early-formed cores [41]. They are also less likely to form new ones given their longer relaxation times, unless they were subject to much higher density and potential fluctuations than those inferred here.

In summary, relaxation times of the order of a Hubble time or less indicate that past relaxation was sufficient to form significant cores in the present sample of galaxies. However, extrapolating the correlation in Figure 4 suggests that core formation through feedback-driven potential fluctuations may become less effective in more massive systems, particularly when assuming no significant contribution from molecular gas fluctuations. As a conservative conclusion, our results demonstrate that stochastic potential fluctuations can be an entirely viable mechanism for explaining halo cores within the standard CDM scenario, at least in systems with up to $M_{\rm vir}=10^{11}M_\odot$ that have not undergone recent cusp reformation. Further analysis of imminent observations, with larger, deeper, and higher-resolution samples, may track the genealogy of core formation and possible cusp reformation, potentially elucidating the observed diversity in the central mass distribution of galaxies.

Methods

Sketch of the theoretical model leading to the relaxation time (equation 1)

In this section, we outline the relationship between density and potential fluctuations, as well as between spatial and temporal ones, within our framework. We then show how the fluctuations drive a random increase in the velocity dispersion of particles forming a halo cusp. In the process, we delineate the simplifying assumptions of the model. We finally discuss the physical picture envisioned and verify the association of the relaxation time with core formation.

Following [15], we first consider a time-snapshot of a gaseous medium with average density $\rho_0$. The gas is driven by an energy source, causing large scale flows and an associated cascade of density and potential fluctuations within a volume $V$. To characterize this, we expand the density contrast $\delta=\frac{\rho(\mathbf{r})}{\rho_0}-1$ and potential $\Phi$ over a volume $V$, so that

and

The Fourier amplitudes are related through the Poisson equation, $\nabla^2\Phi=4\pi G\rho_0\delta$, by $\Phi_{\mathbf{k}}=-4\pi G\rho_0\delta_{\mathbf{k}}k^{-2}$, where $k=|\mathbf{k}|$. The power spectrum of density fluctuations is $P(k)=V\langle|\delta_{\mathbf{k}}|^2\rangle$, and that of the corresponding force fluctuations is $P_F(k)=(4\pi G\rho_0)^2k^{-2}P(k)$. The averages here are strictly speaking ensemble averages, but they may also be considered as average over all snapshots (in which case the process is said to be ergodic). The fluctuations $\delta(\mathbf{r})$ are statistically isotropic if $P(\mathbf{k})=P(k)$, implying that the variance in different spatial directions is the same: $\langle k_x^2\rangle=\langle k_y^2\rangle=\langle k_z^2\rangle=\frac{1}{3}\langle k^2\rangle$.

If the fluctuations are described by time independent statistics (a stationary stochastic process), the power spectra converge to the same values $P(k)=V|\delta_{\mathbf{k}}|^2$, for all snapshots, as $V$ increases. If the statistics of the fluctuations are independent of the point at which they are measured (a spatially stationary, or homogeneous, process), the spatial correlation functions are Fourier transforms of the power spectra. In particular, the force correlation function is

To transform this into a time correlation function, we employ the ’sweeping’ approximations long established in turbulence theory [53, 54, 55]: if each mode with wavenumber $k$ is transported (advected) with velocity $v_k$, such that $r_k(t) = v_k t$, one can rewrite the above equation as

where $p(v_k)$ is the probability distribution function of $v_k$. The derivation of this equation and application to the case of Maxwellian $p$ was studied in [42]. Here we follow [15] and [16] in assuming that the spectrum of fluctuations is swept by large scale fluid flows with a single characteristic speed $v_r$; so that $p(v_k) = \delta(v - v_r)$, and equation (5) reduces to (4) with $r(t) = v_r t$. In this way the spatial statistics of the fluctuating field are transformed into the time domain by being transported, frozen in, by large scale flows with characteristic speed $v_r$. The ergodic and sweeping assumptions were tested against the results of the numerical simulation primarily used here in [16] (Section 3.3).

Using the Newtonian equations of motion, one may deduce that, after time $T$, the gas fluctuations induce velocity changes on a dark matter particle with variance

This is proportional to the kinetic energy injection from the fluctuating gas to the halo particles. Note that as $P_F$ scales as the density squared, so will the force correlation function from equation (4), and consequently so will the kinetic energy injection rate.

For a power law spectrum $P_F \sim k^{-(n+2)}$ (with $n$ the index of density spectrum), one may evaluate this energy, assuming maximal and minimal cutoffs, $k_x \gg k_m$. When the diffusion limit ($v_r t \gg k_m^{-1}$) is taken, equation (1) follows by setting $\left\langle(\Delta v_p)^2\right\rangle= v_p^2$ as is usual in calculations of the relaxation time.

As the correlation function in (6) should be understood as taken along the path of a dark matter particle, $v_r$ is the typical speed of the large scale gas flow with respect to the particle. We follow [15] in associating this with the larger characteristic speed of repeated outflows and inflows confined within spatial scale $d \simeq V^{1/3}$.

The physical picture envisioned is that of persistent fountain flows convecting a cascade of multiscale fluctuations within $V$, rather than outflows permanently exiting $V$. In this context, $v_r$ refers to an average speed corresponding to the flows being in energetic equilibrium with the gravitational potential within $d$.

As the volume $V$ inside which this takes place is finite, an ideal statistical stationarity process as defined above cannot strictly apply, even in principle (due to random variations between snapshots). But a weaker condition of quasi-stationarity may still apply if the power spectra display a well defined mean and variance over time. This is reflected in Figure 1 by the limited variance in the simulated spectra. Physically, ’bursty’ feedback can prompt a quasi-stationary process if it leads to similar (in terms of energy input and effect on the gas) consecutive episodes of fountain flows within $V$, driving density and mass fluctuations that display similar statistics over time. In an isolated simulation this leads to the stability of the power spectra over time (as in Figure 1 and also Figs. 4 and 5 of ref. [16]). In cosmological simulations, the situation is complicated by mass changes arising from accretion and mergers. In that case, the process may be considered piecewise quasi-stationary over limited time (as discussed in § 7.2 below).

In ref. [16], a detailed model for the evolution of the central dark matter distribution under the influence of the gas potential fluctuation was presented. It fit the numerical results with no free parameters unconstrained by the simulation. As those parameters are not available from observations, we here associate the formation time of a core with radius $r_c$ with the relaxation time at $r_c$. This is estimated by setting the dark matter velocity to the circular velocity at $r_c$ in equation (1). As mentioned in the main text, the procedure may be justified by the fact $t_{\rm relax}$ is generally an increasing function of radius. We tested this criterion against the simulation for the same fiducial parameters used in ref. [16]. The result is shown in Figure 5.

The surface densities of observed galaxies

To produce the HI surface density maps of the LITTLE THINGS galaxies, we used the HI flux images created using the procedures described in Section 3 of ref. [56].¹ The HI data cubes, in the image-domain, are thus created with the ROBUST = 0.5 weighting scheme in the AIPS task IMAGER. The cubes are then cleaned with a threshold flux level of two times the RMS noise, as determined in the line-free channels. The data cubes are furthermore corrected for primary beam attenuation, and then de-projected using the photometry parameters (galactic centre, position angle, and axial ratio) determined from the optical V-band morphology [36, 56 table 6]. We compared the results with the original maps, with no de-projection. As the inclinations of the LITTLE THINGS galaxies are generally less than 30 degrees the differences in final surface density profiles were not generally large, although in some cases the central densities are larger by up to a factor of two in the original maps due to larger gas columns resulting from inclination effects. The use of the de-projected maps matches the requirements for evaluating the 3D density (§ 6) and projecting power spectra (§ 7.3), which assume cylindrical symmetry along a projection axis normal to the image.

The HI mass is obtained from the integrated HI flux ($S$) in the robust weighted maps using

where $D$ is the distance to the galaxy and $\Delta V$ is the velocity channel (ch) width. The surface density is then given by:

where $A_{\mathrm{pixel}}$ is the pixel area in $\mathrm{kpc}^2$. The HI surface density profiles $\Sigma(R)$ are then measured in annuli of width equal to the HI beam size.

The collected results for all galaxies in the sample are shown in Figure 6, along with the surface density obtained from the simulation ISIM, which is used to calibrate the fluctuation power spectra as described below. In presenting the results we scale our curves to the half mass and density, as this exploits the near homology reflected in the HI mass-radius relation [57]. This facilitates the comparison of the profiles of systems with different gas masses.

The characteristic density $\rho_0$

Equation (1) involves the 3D average gas density $\rho_0$. If the gas is homogeneous on large scales, $\rho_0$ tends to a definite value when measured over increasing scales $r_0$. In a realistic gas distribution, on the other hand, the density generally decreases with radius, falling off steeper beyond a certain scale. In this case, the relevant effective radius $r_0$ is one inside which $\rho_0(r_0)$ is large enough to ensure significant energy transfer from the fluctuations in the gas to the particles making the dark matter cusp. This radius was estimated for the simulation studied in ref. [16] by measuring the actual energy transfer from the gas to the dark matter within progressively larger radii (cf. Figure 7 of that work). In this case, it was found that it saturated beyond radius $r_{\mathrm{sat}} = 1.3$ kpc, which may in this case be identified with $r_0$. The scale $r_{\mathrm{sat}}$ cannot be deduced this way from data, but it may nevertheless be connected to the properties of the observed density profile.

In ref. [16] it was noted that $r_{\mathrm{sat}}$ approximately matched the radius at which the density fell off steeper than $\rho\sim1/r^2$. Since, as discussed in § 5, the rate of energy transfer from the gas fluctuations to the dark matter scales as the square of the density, this corresponds to a transfer rate rapidly declining at least as $1/r^4$. Furthermore, an $\sim1/r^2$ profile corresponds to surface density $\Sigma\sim1/R$, when projected, we may expect the projected radius where the surface density starts falling off as $1/R$ to also approximate $r_{\mathrm{sat}}$. We checked this by examining the 2D and 3D density slopes in our fiducial simulation. For each snapshot, we calculated the logarithmic slope of the 3D gas density $\alpha= \frac{d\ln\rho}{d\ln r}$ and evaluated the radius at which $\alpha\le b$. To reduce noise, we repeat the calculation for a range of values within $-2.5 \le b \le-2$. We then average the resulting cutoff radii over different $b$ and over simulation snapshots, to obtain the cutoff radius $r_{\alpha=-2}$. We do the same for the surface density, except that in this case $-1.5 \le b \le-1$, to obtain $R_{\alpha=-1}$. The results are shown in Figure 6 by the vertical lines. As can be seen, the two measures are consistent, within the error bars, and close to $r_{\mathrm{sat}}$ (indeed $R_{\alpha=-1}$ turns out to better approximate $r_{\mathrm{sat}}$).

Figure 6, also shows that all three aforementioned measures match the projected half mass radius (at unity). They are likewise close to $R_{\alpha=-1}$ for the galaxy gas density profiles (denoted by $R_{\alpha=-1}^{\mathrm{data}}$). The latter is shown by the vertical dotted line, with dispersion due to scatter between galaxies delineated by the surrounding vertical shaded region. Because the half-mass radius $R_{1/2}$ falls within this scatter, we will use it to define $r_0$ and $\rho_0(r_0)$. Choosing $R_{1/2}$ over $R_{\alpha=-1}$ is advantageous because it is readily defined and evaluated, avoiding the profile density differentiation that introduces significant uncertainties.

In estimating the 3D density at $R_{1/2}$, it is simplest to assume cylindrical symmetry, as the observed surface density is defined in terms of projections on cylindrical shells with radii $R$ on the image planes, along the normal axis $z$ (an analogous argument, assuming the density is stratified along spheroids leads to similar results). The density profile of the gas is thus assumed to be of the form $\rho(R,z) = \rho(R)\rho(z)$, which is constant along cylindrical shells with axis ratio $z/R = q \le1$. For the flattening of dwarf irregular galaxies a characteristic value is $q \simeq0.4$ [58, 59, 60].

The enclosed mass is given in terms of the average density inside the shells by

The average density at the half projected mass radius $R_{1/2}$ is

Here $M_{1/2}$ is the measured projected mass at $R_{1/2}$, and the factor $\eta=M(R_{1/2},qR_{1/2})/M(R_{1/2},qR_M)$ accounts for the fact that $M_{1/2}$ is defined in terms of the whole projected gas column density up to a total (maximum) height $z=qR_M$ normal to the image planes, while the average 3D density is defined for mass inside $z=qR_{1/2}$. For homogeneous density distributions (with sudden cut-off at $(R_M,qR_M)$), $\eta=R_{1/2}/R_M=1/\sqrt{2}$, and is generally closer to unity when the density gradually decreases outwards from the centre.² We then define $\rho_0(r_0)=\bar{\rho}_{1/2}(\eta,q)$ and an effective spherical radius

such that

Note that $r_0 \simeq R_{1/2}$; e.g., $r_0 = (0.8)^{1/3}R_{1/2}$ for $q = 0.4$ and $\eta= 0.75$. The densities $\rho_0(r_0)$, thus calculated, are shown in Figure 7. In this plot, as in all calculations using equation (11), the mean value of $q$ and the uncertainty on it is taken, for each galaxy, from [60], and a flat prior on $\eta$, such that $0.71 \leq\eta\leq0.95$. In calculating the mass, the mass of HI is multiplied by a factor 1.4 to account for helium and other elements.

Simulated power spectra and calibration

Isolated simulation (ISIM HEZFRC)

We are interested in the power spectra of fluctuations in the density contrasts, as their normalization and tilt enter into equation (1). We define these in 2D and 3D through

$\delta_{2D} = \Sigma(R)/\Sigma(R)-1$ and $\delta_{3D} = \rho(r)/\rho(r)-1$. Here $\Sigma(R)$ and $\rho(r)$ refer, respectively, to averages over circular rings and spherical shells at cylindrical radii $R = |\mathbf{R}|$, and spherical radii $r = |\mathbf{r}|$. Over a surface $S = d^2$ and volume $V = d^3$ the corresponding Fourier modes are

In practice, the power spectra are calculated through discretization of these equations and applying a fast Fourier transform (FFT) at a resolution defined by a pixel size $l \ll d$. In the convention used here, the power spectra are computed, at each snapshot, from $P(K) = S|\delta_{\mathbf{K}}|^2$ and $P(k) = V|\delta_{\mathbf{k}}|^2$. The isotropic spectra are obtained by averaging over shells or rings around radial wavenumbers $k$ and $K$ respectively.

On the spectra inferred from the simulation ISIM, we performed various convergence and robustness tests. In Figure 8, we show the effect of pixel size $l$ over which the discrete Fourier transform is evaluated. Except for an expected resolution effect at large wave numbers, the results are stable. This robustness of the spectra at small wavenumbers will be confirmed in further tests discussed below in connection to the observations. As our principal results depend on the normalization of the spectra at small wavenumbers—$P(k_m)$ in (1)—this property renders them also robust.

We also tested the effect of inclination on the 2D and 3D power spectra, by rotating the system relative to axis of angular momentum, and found the results to be insensitive to it. This reinforces the assumption of isotropic fluctuations (cf. § 5).

In the above tests, we used the fiducial box size $d = 5$ kpc (so that $d/2 \simeq2r_{\rm sat}$), following ref. [16]. Variation of the box size $d$ does not introduce very large differences as long as $d/2$ is large enough to incorporate the largest fluctuation scales, but not so large as to contain regions free of fluctuations—that is, beyond the region where episodic feedback-driven gas outflows and inflows are confined. From Fig. 2 of [16], this region satisfies $d/2 \gtrsim2$ kpc.

Thus, while in an idealized situation where fluctuations exist over an arbitrary volume $V$, the power spectra should converge as $V$ increases, when the fluctuations are confined to a finite $V$—as in a realistic simulation, and as assumed in the physical picture envisioned in our model (§ 5)—the spectra will weaken as regions essentially free of fluctuations are incorporated. A compromise is therefore necessary between selecting box sizes that capture large fluctuation scales and avoiding the inclusion of significant fluctuation-free regions.

With this in mind, in Figure 9 we show results for a range of $d$ that generously encompasses all relevant scales, from $d = 2.5$ kpc (which is smaller than the effective density scale) to $d = 6$ kpc (which encompasses a volume larger than that within which fluctuations are expected to be significant). The spectra in this latter case indeed fall below the other lines, which are more consistent within the ranges allowed by the box sizes and the uncertainties shown by the shaded regions. Nevertheless, as we will now see, for the purpose of transformation from 2D to 3D spectra, the differences are less pronounced; the ratios of the 3D to 2D spectra are on average less sensitive to the larger box size than the absolute values, and so the 6 kpc box will also be used for that purpose.

As $P(k_m)$ and $n$ in equation (1) refer to the 3D power spectrum, while only 2D spectra are evaluated from observations, we use snapshots from ISIM to calibrate the 3D spectra to the 2D ones. In the following two subsections we extend the results to systems with different effective gas radial scale ($\simeq r_0$), by examining higher mass haloes from the FIRE-2 simulation and through a simple projection model.

To extract spectral parameters, we define a power-law scaling domain: starting at the midpoint between wavenumbers of 1 and 50 kpc$^{-1}$, we calculate the logarithmic derivative, moving away from the midpoint until the slope deviates by 0.5 from the central value, defining this as the power law scaling region of the spectrum (the relatively large deviation allowed takes into account fluctuations in the logarithmic derivatives due to noise). Within this scaling window, we fit power laws of the form:

where $A_{2\mathrm{D}}$ and $A_{3\mathrm{D}}$ are dimensionless parameters. The results for box sizes between 3 and 6 kpc are shown in Figure 10. We find the 2D and 3D slopes to be nearly identical, with slightly steeper 3D slopes. The amplitudes $A_{3\mathrm{D}}$ are significantly larger than $A_{2\mathrm{D}}$, but their ratio is consistent with being approximately constant. Taken together, these properties imply a quasilinear relation between $P_{3\mathrm{D}}$ and $P_{2\mathrm{D}}$, which may be interpreted in terms of the projection model described below (§ 7.3).

The mean values of the indices are generally related by $n_{3\mathrm{D}} \equiv n \simeq1.1n_{2\mathrm{D}}$. For $d = 5$ kpc, the fiducial value used in [16] and Figure 1, the amplitudes are related through $A = A_{3\mathrm{D}}/A_{2\mathrm{D}} = 6.27$. For the purpose of the transformation of the observed 2D spectra into 3D ones, we use the ratios of the amplitudes and slopes of all simulation snapshots and box sizes used, and estimate the mean values and the associated errors (by evaluating the dispersion). This gives

Comparison with cosmological simulations

A detailed examination of the mechanism studied here in the context of cosmological simulations is better suited for a separate study. Nevertheless, to investigate the presence, or lack thereof, of the essential features required by our model, and to verify the calibration beyond the idealized simulation considered above, we examined the gas properties in haloes identified in two different cosmological simulations for which data are publicly available, and which employed substantially different star formation and feedback recipes.

This included several small haloes from the public release of the Auriga gravo-MHD cosmological zoom-in simulation [38]. The subgrid physics in this suite follows the prescription of Springel & Hernquist [61]. The gas is considered as a two-phase medium, with star formation occurring in the cooler component. The threshold for stochastic star formation is small (at number densities $0.13\,\mathrm{cm}^{-3}$), and the process is regulated by a long timescale (2.2 Gyr), a stiff equation of state, and the removal of gas through large scale winds in proportion to the star formation rate. This makes for more quiescent star formation, with bursty episodes accompanied by large density fluctuations being less likely.

The 2D power spectrum of gas embedded in an example halo from the Auriga Level 3 suite (resolution $\sim30\,\mathrm{pc}$) was already presented in Figure 1. Relative to the observations and our fiducial simulation ISIM, the slope and normalization of the spectrum represent quantitatively and qualitatively different gas properties; in terms of the level of the fluctuations and their distribution among scales, respectively. The central densities of haloes from this simulation were previously analysed, and density cusps, rather than cores, were invariably found [32, 62]. As we will see below, this is also the prediction of our model. Indeed, as we will find, the 3D gas spectra in Auriga lack extended power-law scaling regions. This reflects the absence of a scale-free cascade of density fluctuations that may be carried by larger scale flows, which is a central feature of our model (§ 5), In light of this, and the differences with the 2D spectra of the galaxies, we will not use the Auriga results in verifying the projection model that calibrates 2D to 3D spectra (§ 7.3).

We also examined the properties of the gas in haloes in publicly released data from the FIRE-2 zoom-in simulations [39]. These simulations have higher (pc scale) resolution and different recipe for the gas physics. In particular, the star formation threshold is far higher than in Auriga, and unlike the former is unregulated; its onset is set by the local Jeans scale and it occurs on the local cloud collapse time. This may lead to bursty episodic feedback, blowing holes and driving large scale flows, which can carry large density fluctuation cascades of the type envisioned in our model. As already seen in Figure 1, the spectra of the 2D fluctuations are also closer to those of the galaxies. We will thus use the results from the FIRE suite for verifying the spectral calibration process in § 7.3.

The issue of core formation in the FIRE simulations was first studied in [41]. Cores were found to be present in some haloes, though cusps were also found to reform through mergers, an effect we do not consider here. With the goal of verifying the projection and calibration model in mind, we consider instead the gas properties inside haloes where cores were identified in the FIRE-2 simulation [50], and do this within time intervals that are large enough to evaluate a statistical sample of spectra but which are still small compared to the timescale of halo mass evolution, which necessarily occurs through accretion and mergers in a cosmological context.³

This ensures quasi-stationarity during that period; by which we mean that the power spectra form a homogeneous sample, such that the average remain stable, during the period considered, with modest scatter (§ 5). After some trials, we choose time intervals of $\sim0.5$ Gyr, during periods when no major mergers occur. Such a time interval is significantly larger than the timescale of the largest scale fluctuations ($t_f \sim(k_m v_r)^{-1}$). The separation of timescales allows for multiple realizations of the fluctuation field while the global parameters of the system (like total mass) do not change significantly. Within the considered timescale, a halo particle thus experiences fluctuations with the same quasi-stationary statistics. (This has analogy for example in the application of two body relaxation time on evolving haloes considered in [64]). The total effect of the gas fluctuations, as modelled here, may be considered as a cumulative effect of such quasi-stationary time intervals.

In Figure 11, we compare the gas density contrasts in three haloes from the FIRE simulation with those in ISIM and halo 1 of the Auriga simulation. The lack of significant fluctuations in the latter relative to the FIRE and ISIM cases is clear. This can also be probed through measuring the temporal mass variation (Figure 12). Note that in the case of the Auriga halo we show the evolution over nearly a Hubble time, to illustrate that no significant evolution occurs other than due to mass accretion.

Figure 13 compares the power spectra associated with the density contrast of Figure 11. For further verification of the qualitative difference between the gas spectra in Auriga compared to FIRE and ISIM, we also include an additional halo 0 from Auriga, which has smaller mass but also exhibited qualitatively the same density contrast as halo 1. All spectra from the FIRE suite are much closer to those for ISIM (and to the galaxies) than in the Auriga simulation. The 3D spectra especially show the qualitative difference: the region where a power law may be identified is small, and its slope is extremely steep ($P_{3D}(k) \sim k^{-6}$). Indeed, the large scale power, which is an order magnitude smaller than in FIRE, is largely due to the spiral patterns in a disc-like system rather than the kind of scale-free energy cascade expected from feedback-driven turbulence. A naive application of our model leads to exceedingly large relaxation times. For halo 1, we may take $n = 6$, $v_c = 100$ km/s, $r_0 = 6.4$ kpc, $\rho_0(r_0) = 2 \times10^6$ $M_\odot$ kpc$^{-3}$. This gives $t_{\mathrm{relax}} \simeq1.5 \times10^{13}$ years at $r_0$, or about 1100 Hubble times. For the smaller halo 0—with $P(k_m) = 0.1$ kpc$^{-1}$, $n = 6$, $v_c = 35$ km/s, $r_0 = 5$ kpc, $\rho_0(r_0) = 5 \times10^5$ $M_\odot$ kpc$^{-3}$—one obtains $t_{\mathrm{relax}} \simeq2 \times10^{13}$ years at $r_0$. But the model does not even strictly apply, due to the absence of a wavenumber regime $k_m \ll k \ll k_x$ through which a scale-free cascade of fluctuations can be defined.

The relation between the 2D and 3D spectra for the FIRE haloes are shown in Figure 14. These are analogous to the ones shown in Figure 10 for ISIM, and also complementary, since the haloes in the FIRE simulation have larger effective gas radius. They will thus be used to verify our projection model described below.

Interpretation in terms of projection model

We consider the following projection [65]

where $k_z$ denotes the component of $\mathbf{k}$ that is vertical to the plane of projection (containing the radius vector $\mathbf{R}$), and $\widetilde{W}$ refers to a filter corresponding to a slab of thickness $2L$ (corresponding to the effective galaxy thickness) along that direction. We use a Gaussian filter, which averages the fluctuations such that $\delta_{2D} = \frac{1}{\sqrt{\pi}L}\int\delta_{3D}e^{-z^2/L^2}\,\mathrm{d}z$. The corresponding filter in Fourier space is $\widetilde{W}^{2} = e^{-k_z^2L^2}$. Note that $\delta_{2D}$, thus defined, corresponds to our definition of $\delta_{2D} = \delta\Sigma(\mathbf{R})/\Sigma(\mathbf{R})$ when the average density $\rho(r)$ is constant over scales $\sim L$ (cf. first paragraph of § 7.1 for definition). Thus, equation (16) is an estimate of the effect of projection that is particularly effective when most of the power comes from radii where the density is slowly varying.

Fig. 15 shows the result of using equation (16) for various values of $L$. For illustration, the 3D power law has index $-2$ but the following trends are generic. Two regimes may be isolated. In the limit $L \ll1/k$, the filter permits modes with $k_z \gg k$ to pass, rendering equation (1) equivalent to an unfiltered integral over wavenumbers. In this case, a spectrum $P_{3D} \propto k^{-n}$ is simply projected to $P_{2D} \sim K^{-(n-1)}$ (this is the counterpart, in wavenumber space, to the surface density scaling as $1/R$ when the 3D density goes as $1/r^2$, mentioned in § 6). In contrast, when the line-of-sight projection scale is large relative to the scale of interest ($L \gg1/k$), the filter restricts $k_z$ to small values ($k_z \ll k$), yielding $P_{3D} \sim L P_{2D}$ (this approximates the large $L$ limit of equation (16), which is $P_{3D} = 2\sqrt{\pi}L P_{2D}$). In this quasi-linear limit, $n_{3D}/n_{2D}$ is near unity, matching the results of the ISIM and FIRE simulations discussed above.

We now verify that it is the quasi-linear regime that is also relevant for converting the observed 2D power spectra of our sample of galaxies to 3D spectra. For this, we first need an estimate of the relevant $L$, which one may expect to be of the order of the effective scale of the density distribution, beyond which the gas density declines steeply. By inspecting Fig. 15, we note that a ratio of the 3D to 2D power spectra of order $A \approx6$—as found from the calibration using the fiducial simulation ISIM (eq (15))—corresponds to $L \gtrsim1$ kpc. This matches the energy transfer saturation scale $r_{\mathrm{sat}}$ and the closely related measures of the radial scale setting the effective density, as discussed in § 6 (Figure 6). It suggests that the prediction of the simple projection model matches that of the simulation when $L$ is indeed of the order of that effective scale. As we will see below, the results of the analysis of the gas in the FIRE-2 haloes are also consistent with this.

For galaxies, the measure we adopted for the effective gas radial scale was $r_0$, which was derived from the half mass projected radius. We will now assume that the scaling $L = r_0$ holds for the galaxies and examine whether that renders them in the quasilinear regime of conversion between 2D and 3D spectra. As noted above, significant departures from linearity in the conversion from 3D to 2D are clearest at smaller $k$. To verify their absence, we thus calculate the ratio $P_{3D}(k_m)/P_{2D}(K = k_m)$ as a function of $L = r_0$ for each galaxy, using equation (16) and the values of $k_m$ found in the subsection below. The results are shown in Figure 16. There we also plot two lines where equation (16) is used with large values of $k_m$, defining the quasi-linear regime, as reference. As smaller galaxies have smaller $L = r_0$ but larger $k_m$, the ratio between the 3D and 2D spectra invariably lies in the nearly linear segment of the relation ($r_0 \gg1/k_m$). We thus infer that the galaxies lie in the same quasilinear projection regime as the gas in ISIM and also the FIRE haloes, which are seen to be consistent with the projection model within the error bars.

Given the general consistency of the projection model with the simulation results (with $L = r_0$), and that the galaxies lie in the same quasilinear regime $P_{3D} \sim LP_{2D}$ as the simulations, we replace equation (16) by extrapolating the results of the calibration from the simulations as

where $P_{2\mathrm{D}}$ and $P_{3\mathrm{D}}$ are assumed to be given as in equations (14). This is readily inverted and solved for $P_{3\mathrm{D}}$, given $P_{2\mathrm{D}}$. The mean values and errors on $P_{2\mathrm{D}}$ (i.e on the normalization $A_{2\mathrm{D}}$ and index $n_{2\mathrm{D}}$) are obtained as described in § 8 and § 9 (item 2). The factor $\beta= r_{\mathrm{sat}}/r_0$ scales the estimated effective density scale of each galaxy relative to the effective density scale in ISIM. The means and errors of $r_0$ for each galaxy are obtained as described in § 6 and § 9. The weak departure from linearity is incorporated through the index $n_{3\mathrm{D}} - n_{2\mathrm{D}} = n_{2\mathrm{D}}(n_r - 1)$. As we have not comprehensively investigated the results of the FIRE simulations, compared to ISIM, we only use the latter for estimating $\beta$, as well $A$ and $n_r$ (as given by equations 15).

Observed spectra and minimal wavenumber

To produce power spectra from the galaxy sample data, we use the image maps (cf. § 5.1). The use of the image maps, rather than the interferometric data directly [67], is preferred because we are interested in quantifying relative fluctuations of the surface density $\delta\Sigma/\Sigma= \Sigma(R)/\bar{\Sigma}(R)-1$, and therefore positional information in real space is required (in order to obtain a density field $\Sigma(R)$ and scale it with $\bar{\Sigma}(R)$ before taking the spectrum). Discussion and comparison of the two methods for some of the galaxies studied here can be found in ref [68].

As with the simulation ISIM, we checked that varying the pixel size over which the power spectra are evaluated does not significantly affect the power spectra. Figure 17 compares the beam scale power spectra to the cleaned image pixel scale spectra. The fluctuations at the largest scales, which are crucial to our estimates of the relaxation time, again remain relatively stable and robust. On the other hand, small scale noise is smoothed and attenuated at large $K$ by the Multiscale CLEAN algorithm used to produce the images (§ 9).

We also varied the box size and studied the convergence of the spectra. As with the simulation, when fixing the box size a compromise is made. The box size $d$ has to be large enough to capture the largest fluctuations at minimal wavenumber $k_{\mathrm{m}}$. For the observations, the outer limit is effectively set by requiring $d$ not to be so large as to introduce spurious discreteness effects and noise near the observational minimal surface density limit of $10^{19}$ atoms per $\mathrm{cm}^{2}$. For our sample, this generally corresponds to radii $\lesssim2R_{1/2}$, which for ISIM (where $R_{1/2} \simeq1.5\,\mathrm{kpc}$), translates to $d \lesssim6\,\mathrm{kpc}$. For the galaxies we are therefore invariably within the box sizes considered in Figure 9, even when the FFT is taken over the whole observed gas distribution.

Figure 18 shows the convergence of the spectra for various radii comprising different percentages of the enclosed projected gas mass. We find general convergence when $d/2 \gtrsim R_{90}$, where $R_{90}$ is the cylindrical radius containing 90% of the observed gas mass. To test this, and obtain an error estimate, we estimated the mean and variance of the distribution of the power spectra fitting parameters for half box sizes between radii containing 85% and 95% of the mass, as described in § 9. The slopes and normalization of the power law parts of the galaxy spectra are determined starting at the midpoint between the two vertical lines shown in Figure 18 (i.e. between the minimal wavenumber and that corresponding to the beam scale), and applying the same procedure described above in relation to the simulations (equation 14).

In order to complete the quantification of the ingredients in equation (1), we still need to estimate the wavenumber associated with the largest fluctuation scale $k_m$. We do not have direct access to the 3D power spectra of the observed galaxies, and the turnover from the power law at large scales in the 2D spectrum can, in principle, result from projection. Moreover, power at small wavenumbers can arise from large scale asymmetries and spiral structures, rather than the scale-free cascade of fluctuations relevant to our model, as was apparent in the case of the Auriga simulation gas. It may therefore be more accurate to attempt to infer $k_m$ in equation (1) in an independent manner. To this end, we use the results of ref. [66], which determined maximal cloud sizes for galaxies in our sample. This has the added benefit that the largest sizes more likely correspond to actual hydrodynamic fluctuations rather than from structures that may arise from large scale deformations, which may not be as relevant to our model. In this way, we associate the radius $R_{\mathrm{cloud}}$ with $1/k_m$ (noting that the cloud radius corresponds to 3D fluctuations characterized by wavenumbers $k$, rather than projected fluctuations expressed in terms of $K$). We also note that, for the purpose of illustration, in Figure 18 we transform the inferred $k_m$ into $K_m$ assuming isotropic fluctuations ($\langle K^2\rangle= \langle\frac{2}{3}k^2\rangle$), such that the vertical lines correspond to $K_m = \sqrt{\frac{2}{3}}k_m$. But this transformation does not enter into the calculation of the relaxation time using equation (1), as this depends directly on the 3D spectrum (and thus involves an inverse transformation from $K$ to $k$). The relation between $P_{3D}(k_m)$ and $P_{2D}(K = k_m)$ over the power law part of the spectra is thus obtained using equation (17).

Error estimation and propagation

Equation (1) contains two terms related to the power spectra, namely $P(k_m)$ and $n$; another two related to the mass model ($v_p$ and $v_r$); in addition to the effective gas density $\rho_0$. We calculated error estimates on all these quantities, some of which were mentioned above when relevant. All errors were propagated through a Monte Carlo procedure. We adopt the Latin-hypercube MC sampling and perform non-order-specific error propagation via the Python package mcerp ⁴. The MC sampling is conducted over a number of points of order $10^6$.

We note that we do not place error estimates on the HI images themselves, as these are obtained from the ’dirty’ beams using Multiscale CLEAN [56], which is an inherently nonlinear procedure with no readily quantifiable errors. A reassuring result is that the final residuals are small and consist of nearly perfectly Gaussian white noise. The residual are furthermore subtracted in the noise-free images we use here, but adding them gives an idea of the errors involved. We estimate these from the robust data-cube signal-free channels (as described in CASA-imstat task [69]). A typical example is shown in Figure 19 where it is clear that the effect is minimal.

The error estimates that were included in the calculations are as follows.

1. The errors on the mass model are obtained from ref. [43], where Table 2 gives the core radius and characteristic density of the isothermal halo fits to the data. Monte Carlo distributions are produced for these quantities and the circular velocity at the core radius $v_c$ and $v_r = d\sqrt{G\bar{\rho}_{\rm DM}(d/2)}/2$ are estimated by propagating the errors as described above.

2. Errors on the observed (2D) power spectra include those on the normalization $A_{2D}$ and the index $n_{2D}$ of the power law fits. We also include in our estimate the error entering from our assumed convergence box size, which has sides $\simeq2 \times R_{90}$ (where $R_{90}$ is the radius containing 90 percent of the projected mass). To this end, we take a distribution of power spectra including radii between $R_{85}$ and $R_{95}$, fit power laws to this sample of spectra and calculate the mean values—which are statistically convergent and usually centred close to the value at $R_{90}$—as well as the variances of the normalization and index of the spectrum.

3. We use the isolated dwarf galaxy simulation (ISIM) to estimate the error in the conversion of the 2D spectra into 3D. Similar to the above procedure, we stack all power spectra from all snapshots and including all box sizes used (i.e., 3, 4, 5, and 6 kpc), and estimate the errors in the conversion fits to the normalization and index conversion factors $A_{3D}/A_{2D}$ and $n_r = n_{3D}/n_{2D}$ (Figure 10). This leads to equations (15).

4. We then invert equation (16) to obtain the 3D power spectra, including the errors in $n_{2D}$, calibration errors in $A$ and $n_r$, and also errors in $r_0$ (described below). This is done through equation (17), which is verified by the data and simulation points (including from FIRE-2) in Figure 16.

We finally introduce mean values and errors in $k_m$ from the maximum cloud size values, as given in [66]. These are propagated to obtain estimates of $P(k_m)$.

5. The effective radius $r_0$ and density $\rho_0(r_0)$, as given by equation (11), incorporate errors in the flattening $q$, which for each individual galaxy is obtained from ref. [60]. It moreover includes the uncertainty in the density distribution, parametrized by

$\eta$, for which we use a flat prior comprising reasonably expected values, $0.71 \le\eta\le0.95$, as discussed in § 6.

6. For Figure 4, the virial masses are estimated as in [43], by assuming an initial NFW [70] profile and associating its $V_{200}$ with the maximal speed of the isothermal halo fit, then using the empirical relation $M_{200} \simeq3.29 \times10^{5}V_{200}^{3}$ [71]. The errors on $M_{200}$ are propagated from those on $V_{200}$ (as given in Table 2 of [43]).

Table 1 displays the mean values with uncertainties of quantities entering into equation (1).

Galaxy	$r_{0}$	$\log(\rho_{0})$	$n_{\mathrm{Obs}}$	$n_{\mathrm{Obs+Calib}}$	$P^{\mathrm{Obs}}_{3D}\ (km)$	$P^{\mathrm{Obs+Calib}}_{3D}\ (km)$	$v_{p}$	$v_{r}$	$\log(t_{\mathrm{relax}})$
CVnIdwA	$0.73 \pm0.15$	$7.25^{+0.25}_{-0.17}$	$3.03 \pm0.28$	$3.03 \pm0.37$	$0.04^{+0.03}_{-0.02}$	$0.04^{+0.05}_{-0.02}$	$19.77^{+5.68}_{-5.35}$	$8.06^{+0.95}_{-1.08}$	$0.09^{+0.64}_{-0.76}$
DDO43	$2.14 \pm0.45$	$6.41^{+0.25}_{-0.17}$	$2.49 \pm0.23$	$2.49 \pm0.30$	$14.10^{+9.06}_{-5.62}$	$14.07^{+14.18}_{-7.11}$	$18.54^{+4.57}_{-4.26}$	$16.96^{+3.73}_{-3.61}$	$-0.59^{+0.58}_{-0.71}$
DDO47	$5.02 \pm0.89$	$5.67^{+0.22}_{-0.15}$	$2.68 \pm0.12$	$2.68 \pm0.24$	$17.79^{+5.42}_{-4.31}$	$17.73^{+13.20}_{-7.61}$	$44.95^{+6.54}_{-6.30}$	$33.14^{+3.40}_{-3.49}$	$1.90^{+0.45}_{-0.55}$
DDO50	$3.12 \pm1.03$	$6.46^{+0.33}_{-0.22}$	$1.96 \pm0.07$	$1.96 \pm0.16$	$2.94^{+0.73}_{-0.59}$	$2.94^{+2.09}_{-1.22}$	$8.34^{+8.69}_{-8.33}$	$8.66^{+8.88}_{-8.64}$	$-1.39^{+1.27}_{-0.00}$
DDO52	$3.99 \pm0.16$	$5.85^{+0.05}_{-0.05}$	$1.98 \pm0.42$	$1.98 \pm0.44$	$9.94^{+5.99}_{-3.77}$	$9.93^{+9.54}_{-4.87}$	$32.32^{+2.10}_{-2.05}$	$28.18^{+1.59}_{-1.58}$	$1.35^{+0.35}_{-0.35}$
DDO53	$1.10 \pm0.23$	$6.76^{+0.25}_{-0.17}$	$2.07 \pm0.35$	$2.07 \pm0.38$	$0.71^{+0.53}_{-0.31}$	$0.71^{+0.80}_{-0.38}$	$38.88^{+33.44}_{-26.19}$	$18.64^{+4.31}_{-8.43}$	$0.51^{+0.91}_{-1.47}$
DDO70	$1.32 \pm0.28$	$6.42^{+0.25}_{-0.17}$	$2.37 \pm0.18$	$2.37 \pm0.25$	$0.10^{+0.08}_{-0.04}$	$0.10^{+0.12}_{-0.05}$	$19.03^{+5.29}_{-4.89}$	$17.12^{+4.16}_{-4.04}$	$1.55^{+0.64}_{-0.77}$
DDO87	$3.96 \pm0.84$	$5.77^{+0.25}_{-0.17}$	$2.42 \pm0.14$	$2.42 \pm0.23$	$18.17^{+7.45}_{-5.46}$	$18.10^{+14.66}_{-8.15}$	$31.93^{+1.69}_{-1.66}$	$25.16^{+1.05}_{-1.05}$	$1.25^{+0.45}_{-0.59}$
DDO126	$1.84 \pm0.07$	$6.58^{+0.05}_{-0.05}$	$2.78 \pm0.46$	$2.78 \pm0.51$	$2.11^{+1.66}_{-0.95}$	$2.11^{+2.45}_{-1.15}$	$21.47^{+1.92}_{-1.87}$	$16.68^{+1.16}_{-1.17}$	$0.13^{+0.42}_{-0.42}$
DDO133	$1.72 \pm0.13$	$6.53^{+0.10}_{-0.09}$	$3.08 \pm0.56$	$3.08 \pm0.60$	$4.54^{+2.94}_{-1.79}$	$4.54^{+4.56}_{-2.28}$	$24.75^{+2.23}_{-2.16}$	$20.88^{+1.59}_{-1.58}$	$0.15^{+0.40}_{-0.41}$
DDO154	$3.14 \pm0.13$	$6.19^{+0.05}_{-0.06}$	$1.84 \pm0.27$	$1.84 \pm0.31$	$2.41^{+1.99}_{-1.09}$	$2.41^{+2.82}_{-1.30}$	$24.12^{+1.05}_{-1.04}$	$22.87^{+0.94}_{-0.93}$	$0.91^{+0.38}_{-0.39}$
DDO168	$1.80 \pm0.38$	$6.85^{+0.25}_{-0.17}$	$2.15 \pm0.09$	$2.15 \pm0.19$	$0.55^{+0.20}_{-0.15}$	$0.55^{+0.44}_{-0.25}$	$61.18^{+19.13}_{-18.37}$	$40.13^{+6.74}_{-8.11}$	$1.27^{+0.58}_{-0.74}$
DDO210	$0.28 \pm0.06$	$7.15^{+0.25}_{-0.17}$	$2.71 \pm0.58$	$2.71 \pm0.62$	$0.11^{+0.18}_{-0.07}$	$0.11^{+0.23}_{-0.08}$	$7.42^{+2.07}_{-1.95}$	$6.02^{+1.26}_{-1.31}$	$-1.22^{+0.78}_{-0.89}$
DDO216	$0.42 \pm0.07$	$7.08^{+0.20}_{-0.14}$	$2.81 \pm0.87$	$2.81 \pm0.89$	$0.21^{+0.80}_{-0.17}$	$0.21^{+0.92}_{-0.18}$	$5.74^{+1.63}_{-1.52}$	$5.55^{+1.47}_{-1.41}$	$-1.60^{+1.00}_{-1.06}$
F564-V3	$1.73 \pm0.08$	$6.12^{+0.25}_{-0.06}$	$1.80 \pm0.63$	$1.80 \pm0.65$	$2.69^{+4.37}_{-1.72}$	$2.69^{+5.09}_{-1.82}$	$15.36^{+7.44}_{-6.77}$	$13.96^{+6.44}_{-6.15}$	$0.31^{+1.05}_{-0.83}$
IC10	$1.29 \pm0.27$	$6.61^{+0.17}_{-0.15}$	$1.50 \pm0.19$	$1.50 \pm0.22$	$0.03^{+0.01}_{-0.01}$	$0.03^{+0.02}_{-0.02}$	$12.28^{+5.99}_{-5.79}$	$12.06^{+5.70}_{-5.15}$	$0.85^{+0.83}_{-1.16}$
IC1613	$1.47 \pm0.04$	$6.25^{+0.04}_{-0.05}$	$2.35 \pm0.22$	$2.35 \pm0.28$	$0.14^{+0.08}_{-0.05}$	$0.14^{+0.15}_{-0.07}$	$3.03^{+0.68}_{-0.65}$	$2.97^{+0.62}_{-0.60}$	$-0.56^{+0.46}_{-0.47}$
NGC1569	$1.82 \pm0.05$	$6.84^{+0.04}_{-0.04}$	$1.65 \pm0.15$	$1.65 \pm0.20$	$0.33^{+0.16}_{-0.11}$	$0.33^{+0.31}_{-0.16}$	$36.49^{+11.52}_{-11.09}$	$26.43^{+5.94}_{-1.03}$	$0.81^{+0.53}_{-0.38}$
NGC2366	$2.84 \pm0.26$	$6.70^{+0.13}_{-0.10}$	$1.69 \pm0.21$	$1.69 \pm0.25$	$0.73^{+0.26}_{-0.24}$	$0.73^{+0.35}_{-0.16}$	$27.90^{+1.22}_{-1.21}$	$25.64^{+1.03}_{-1.03}$	$0.54^{+0.38}_{-0.40}$
UGC8508	$0.62 \pm0.02$	$7.13^{+0.04}_{-0.04}$	$2.47 \pm0.55$	$2.47 \pm0.58$	$0.05^{+0.11}_{-0.04}$	$0.05^{+0.13}_{-0.04}$	$45.66^{+5.06}_{-5.01}$	$16.38^{+0.54}_{-0.57}$	$1.23^{+0.64}_{-0.61}$
WLM	$1.10 \pm0.05$	$6.94^{+0.06}_{-0.06}$	$1.56 \pm0.04$	$1.56 \pm0.13$	$0.06^{+0.02}_{-0.02}$	$0.06^{+0.05}_{-0.03}$	$19.53^{+0.37}_{-0.37}$	$16.67^{+0.29}_{-0.37}$	$0.60^{+0.28}_{-0.28}$
Haro29	$1.55 \pm0.33$	$6.40^{+0.25}_{-0.17}$	$1.69 \pm0.37$	$1.69 \pm0.39$	$0.29^{+0.33}_{-0.16}$	$0.29^{+0.43}_{-0.18}$	$15.97^{+4.07}_{-3.86}$	$16.04^{+3.96}_{-3.81}$	$0.77^{+0.71}_{-0.82}$

The inner halo density slopes

In the model used here the relaxation time goes to zero as one approaches the centre of an initial NFW [70] halo, where both circular speeds and velocity dispersions (and hence $v_{\mathrm{p}}$) go to zero. Whenever our mechanism is effective, a small core should therefore immediately form, expanding outwards with time as in Figure 5 (although its size is eventually limited by the energy transfer radius $r_{\mathrm{sat}}$, and also the NFW scale length $r_{\mathrm{s}}$; cf. [16] Appendix C).

Thus the expectation is to find nearly constant density cores, unless the relaxation time is exceedingly large (or recent core reformation occurred), so as to make a core undetectable at any resolvable radius.

In ref. [43], both the cored isothermal mass models and NFW halo profiles were tried in order to fit the data. It was found that, except for a few cases, the latter were entirely unsuitable (e.g., resulting in negative concentration parameters $c$). The inner density profiles obtained do vary significantly because of the observational errors, but the median profiles were indeed consistent with an inner nearly constant density core (Figure 5 in [43]). An attempt was also made to estimate the inner logarithmic slopes of the inner density profiles. This is necessarily sensitive to observational errors, as measuring the logarithmic slope effectively requires evaluating a second derivative of the rotation curve (and the estimation of the inner radius and its associated error involves a third derivative).

We plot the result of their analysis against the relaxation time in Figure 20. Except for DDO 47, the points cluster in a column, with larger error bars and central values of −0.28. Many are consistent with a flat core. Some of the more negative slopes can in principle result from resolution effects, as the transition from $\alpha= -2$ to $\alpha= 0$ occurs over more than an order of magnitude in radius [43, 72, 73]. Except for DD 47, the relaxation times are of the order of Hubble time or less.

The case of DDO 47 is interesting because the rotation curve from which the slope is inferred actually exhibits a nearly flat inner part, corresponding to a steep density cusp. This is followed, at larger radii, by an intermediate region where the inferred circular velocity rises in accordance with a nearly constant density distribution, before flattening again (cf. Figure 19 in [43]). In principle, the associated apparent flattening of the density distribution at intermediate radii could be due to non-circular motion, which may be in turn induced by a triaxial halo, a warp or a spiral structure [74]. The complex behaviour of the rotation curve of DDO 47 is not apparent in earlier work, where such complications were found to be minor and the curve well fitted by a model with a large core (of 7 kpc [75]). This core (which is larger still compared to the mean value $r_c = 3.72$ kpc used here) cannot be formed through the mechanism discussed in the present work, unless the fluctuation levels or the gas density were far larger in the past. If, on the other hand, the strong inner cusp with $\alpha\simeq-1.25$ is the more faithful description of the inner profile, this would be in accordance with the model presented here. For, using column 2 and 6 of Table 2 of [43], one derives a resolution of 0.26 kpc for this galaxy. On the other hand, to bring down the relaxation time to within a Hubble time we need to decrease $v_c^2$ in equation (1) by about 14 times. From Figure 19 of ref. [43], we find that this brings the core inside the resolution limit, and therefore cannot be detected.

Data availability

The observational HI data cubes and moment maps analyzed in this study are publicly available from the NRAO LITTLE THINGS repository (https://science.nrao.edu/science/surveys/littlethings/data/data). The cosmological simulation datasets are accessible through the public data releases of the FIRE-2 project (https://fire.northwestern.edu/data/) and the Auriga Project (https://wwwmpa.mpa-garching.mpg.de/auriga/data.html). All other intermediate data supporting the plots and conclusions of this study are available from the corresponding author upon reasonable request.

Code availability

A Python package to calculate the power spectra of observed and simulated galaxies is publicly available on GitHub at GalaxyPower.

Correspondence

Correspondence and requests for materials should be addressed to A. A. El-Zant (email: [email protected]).

Author contributions

MH undertook all analysis of simulation and observational data and wrote Section 2.1 of the Methods. AAEZ. conceived the project and led the interpretation of the results, in collaboration with M.H. and FC, and wrote the bulk of the paper. F.C. guided the analysis and interpretation of the observational data and contributed to the text and its finalization.

Funding

This work was supported by the Science and Technology and Innovation Funding Authority (STDF), Egypt, under grant number 48289.

Competing interests

The authors declare no competing interests.