The WISDOM of power spectra: how the galactic gravitational potential impacts a galaxy’s central gas reservoir in simulations and observations

Observations indicate that the central gas discs are smoother in early-type galaxies than their late-type counterparts, while recent simulations predict that the dynamical suppression of star formation in spheroid-dominated galaxies is preceded by the suppression of fragmentation of their interstellar media. The mass surface density power spectrum is a powerful tool to constrain the degree of structure within a gas reservoir. Specifically here, we focus on the power spectrum slope and aim to constrain whether the shear induced by a dominant spheroidal potential can induce sufficient turbulence to suppress fragmentation, resulting in the smooth central gas discs observed. We compute surface density power spectra for the nuclear gas reservoirs of fourteen simulated isolated galaxies and twelve galaxies observed as part of the mm-Wave Interferometric Survey of Dark Object Masses (WISDOM) project. Both simulated and observed galaxies range from disc-dominated galaxies to spheroids, with central stellar mass surface densities, a measure of bulge dominance, varying by more than an order of magnitude. For the simulations, the power spectra steepen with increasing central stellar mass surface density, thereby clearly linking the suppression of fragmentation to the shear-driven turbulence induced by the spheroid. The WISDOM observations show a different (but potentially consistent) picture: while there is no correlation between the power spectrum slopes and the central stellar mass surface densities, the slopes scatter around a value of 2.6. This is similar to the behaviour of the slopes of the simulated galaxies with high central stellar mass surface densities, and could indicate that high shear eventually drives incompressible turbulence.


INTRODUCTION
How galaxies become quiescent and cease their star formation activity despite the presence of molecular gas reservoirs, remains an open question.Empirically, the total star formation rates (SFRs) of late-type, star-forming galaxies are proportional to the mass surface densities of their gas reservoirs (e.g.Kennicutt 1998;Bigiel et al. 2008;de los Reyes & Kennicutt 2019).By contrast, the ≈ 25 per cent of early-type galaxies (ETGs) which host molecular gas (e.g.Young et al. 2011;Davis et al. 2019) exhibit lower SFRs despite hosting gas reservoirs with similar mass surface densities (Davis et al. 2014).A similar behaviour is seen in the central, bulge-dominated region of the Milky Way (Longmore et al. 2013;Kruĳssen et al. 2014).Additionally, there is some evidence of intermediate-redshift ETGs that are quiescent despite hosting molecular gas reservoirs (Suess et al. 2017;Williams et al. 2021).
In recent years, the 'dynamical suppression' of star formation ★ E-mail: jindra.gensior@uzh.ch(also known as 'morphological quenching' 1 ; Martig et al. 2009) has emerged as a viable pathway to quiescence in bulge-dominated galaxies (e.g.Martig et al. 2009Martig et al. , 2013;;Gensior et al. 2020;Gensior & Kruĳssen 2021;Kretschmer & Teyssier 2020).Dynamical suppression describes the process whereby shear, induced by the spheroidal component of the gravitational potential, drives turbulence in the gas, thereby stabilising it against collapse and suppressing star formation.Although dynamical suppression requires low gas fractions to be effective (Martig et al. 2013;Gensior & Kruĳssen 2021), it offers a compelling explanation for the suppressed SFRs observed in ETGs with gas reservoirs.Furthermore, dynamical suppression is predicted to maintain quiescence in a spheroid-dominated galaxy when its cold interstellar medium (ISM) is replenished through gasrich minor mergers and stellar mass loss (e.g.Davis et al. 2019).
A key step for the dynamical suppression of star formation is the (dynamical) suppression of fragmentation.Gensior et al. (2020) simulated ten isolated galaxies with an identical stellar mass, but morphologically ranging from pure disc to spheroid.The galaxies with dominant bulges formed smooth, dense circumnuclear gas discs that increased in size with increasing central stellar mass surface density ( * , a proxy for the dominance of the spheroidal component).By comparison, the low- * galaxies in the Gensior et al. (2020) sample host ISM that are porous, sub-structured and clumpy throughout.Davis et al. (2022) used non-parametric morphological indicators to quantify the cold ISM structure of the galaxies from Gensior et al. (2020), the Physics at High Angular resolution in Nearby GalaxieS (PHANGS) survey (e.g Leroy et al. 2021b,c) and the mm-Wave Interferometric Survey of Dark Object Masses (WISDOM) 2 sample.They found strong anti-correlations between the central stellar surface density and the Gini, Smoothness and Asymmetry coefficients in both simulations and observations, demonstrating that galaxies with a deep central gravitational potential (high  * ) have a more smooth, symmetric and regular cold ISM (low Gini, Asymmetry, Smoothness).In this paper we wish to further quantify the effects the gravitational potential and galactic dynamics have on the gas reservoir of a galaxy and its turbulent state.To do so, we turn to an ISM statistic that has been used numerous times to analyse turbulence in galaxies: the spatial power spectrum of the ISM.
The power spectrum encodes information on the nature of the turbulence (e.g.Elmegreen & Scalo 2004).For example, the slopes of the density and velocity power spectra can be used to distinguish whether the turbulence is incompressible (Kolmogorov 1941), or compressible (e.g.Burgers 1948;Fleck 1996), and whether the turbulent forcing is solenoidal or compressive (e.g.Federrath 2013;Nandakumar & Dutta 2020, for an application to a spiral galaxy).Power spectra have been extensively studied within regions of the Milky Way (e.g.Crovisier & Dickey 1983;Green 1993;Stutzki et al. 1998;Lazarian & Pogosyan 2000;Miville-Deschênes et al. 2003;Swift & Welch 2008;Miville-Deschênes et al. 2010;Martin et al. 2015;Pingel et al. 2018).Over the last two decades this has been extended to external galaxies such as the Local Group (e.g.Stanimirovic et al. 1999 for the Small Magellanic Cloud and Elmegreen et al. 2001,Block et al. 2010and Colman et al. 2022 for the Large Magellanic Cloud, see also Koch et al. 2020 for both Magellanic Clouds, M31 and M33) and nearby spirals (e.g.Dutta et al. 2013) and dwarfs (e.g.Dutta et al. 2009b;Zhang et al. 2012).
Crucially for our purpose, the density power spectrum is very sensitive to the morphology of the tracer (Koch et al. 2020).Specifically, the absence of fragmentation shifts power from the smaller to the larger spatial scales (Renaud et al. 2013;Grisdale et al. 2017).Therefore, a suppression of fragmentation in the data will manifest itself as a steep power spectrum.However, there is some debate surrounding the shape (and slope) of extragalactic power spectra.Some studies find that their data are better represented by a broken powerlaw (e.g.Dutta et al. 2008Dutta et al. , 2009b;;Block et al. 2010;Combes et al. 2012), where the small-scale power-law has a steeper slope than the large-scale one.This break in the power spectrum is hypothesised to indicate the transition from three-dimensional turbulence on small scales to two-dimensional turbulence on large scales, the location of the break coinciding with the scale height of the galaxy's gas disc (e.g.Dutta et al. 2008).Alternatively, based on numerical simulations, Körtgen et al. (2021) proposed that the break scale in cold gas traces the average size of molecular clouds, which is not necessarily 2 https://wisdom-project.org/ equal to the disc scale height.Yet, a break in the power spectrum is far from ubiquitous.(e.g.Stanimirovic et al. 1999;Dutta et al. 2009a;Zhang et al. 2012;Dutta et al. 2013;Koch et al. 2020).Recently, Koch et al. (2020) found evidence that the break could be an artefact from not accounting for the shape of the point spread function of the data and can otherwise result from a bright, dominant source (see also Willett et al. 2005).Similarly, different methods of calculating the power spectrum (in particular, whether the uv-plane visibility or the zeroth-moment, mass surface density map is used in the computation3 ) and different tracers (e.g.Combes et al. 2012;Koch et al. 2020 and references therein) can yield very different power spectrum slopes.The above could imply that the density power spectrum is a poorly constrained probe of ISM turbulence, at least when considering whole galaxies.However, given its sensitivity to the morphology of the gas, and given the ability to systematically apply the same procedure to the data, the density power spectrum should be a good tool to quantify the impact of the gravitational potential on the turbulence within a galaxy.If the smooth, circumnuclear gas reservoirs of bulge-dominated galaxies are the result of shear-driven turbulence, we expect a continuous steepening of the power spectrum slope going from disc-dominated galaxies with a porous and sub-structured ISM to bulge-dominated ones.
Given the good agreement between non-parametric morphological indicators (and their trends with galaxy type) between the Gensior et al. (2020) simulated galaxies and observed galaxies in the local universe (see Davis et al. 2022), we wish to compare the power spectra findings of the simulations to observations.To do so, we turn to the WISDOM project (e.g.Onishi et al. 2017;North et al. 2019;Smith et al. 2019;Liu et al. 2021;Davis et al. 2022), which akin to the simulations includes galaxies with diverse stellar morphologies.The primary aim of WISDOM is to accurately measure the mass of each galaxy's central super-massive black hole (SMBH) using the molecular gas dynamics.As a result, exquisite, high-resolution CO maps of the circumnuclear gas reservoirs are available for a variety of galaxies, ranging from late to early types.This makes the WISDOM sample the ideal comparison sample for our purposes.
This paper is structured as follows.In Section 2 we give an overview of the Gensior et al. (2020) simulations and the WISDOM observations used, as well as of our method to measure the density power spectrum.The power spectra of simulations and observations are presented in Section 3 and their trends are discussed in Section 4. We summarise our findings and conclude in Section 5.

Simulations set-up
We analyse the ISM morphologies of the isolated galaxies presented in Gensior et al. (2020), that were simulated with a dynamicsdependent star formation efficiency (SFE), as well as the ISM morphology of four additional simulations.We briefly describe the initial conditions and simulation physics below, but refer the reader to Gensior et al. (2020) for a more complete description.
The Gensior et al. (2020) suite comprises ten simulations, where each isolated galaxy consists of a gaseous and a stellar component  2020), while the latter four have been added here.Column 1 lists the name of each simulation, following the naming convention of Gensior et al. (2020): the prefix refers to the presence of a bulge, the subsequent 'MX' and 'RY' to the bulge mass fraction in percent and the bulge scale radius in kpc, respectively.Column 2 lists the initial bulge mass, column 3 lists the bulge-to-disc stellar mass ratio, column 4 lists the bulge scale radius and column 5 lists the central stellar mass surface density  * .within a Hernquist (1990) dark matter halo.In all cases excepting one (run noB), the stellar contribution has been split into an exponential disc and a Hernquist (1990) bulge component.The total initial stellar mass is fixed at  * ≈ 4.71 × 10 10 M ⊙ , but the initial mass fraction of the bulge component,  b , is varied from 30 to 90 per cent, thereby varying the bulge-to-disc and bulge-to-total mass ratio between the simulations.Simultaneously the initial bulge scale radius,  b , is varied from 1 to 3 kpc across the simulations to create a range of different gravitational potentials.All galaxies are initialised with a gas-to-stellar mass ratio  gas / * = 0.05, in good agreement with the cold gas fraction of massive galaxies in the local Universe (e.g.Saintonge et al. 2017), and within the upper range of gas fractions of early-type galaxies (Young et al. 2011).
The central stellar mass surface density is defined as where  e is the effective radius.For the simulations, we use the stellar half-mass radius, i.e. the radius in the - plane of the galaxy which contains 50 per cent of the total stellar mass, to compute  * .Similar to the gas fraction and bulge-to-disc ratio, the central stellar mass surface density changes negligibly during runtime.Hence, we refer to all simulated galaxies with their initial conditions.Due to the sparsity of simulated galaxies with log( * /M ⊙ kpc −2 ) ≳ 8.75 in Gensior et al. (2020) compared to the WISDOM galaxies in our sample, we extend the Gensior et al. (2020) suite by performing four additional simulations.Initial conditions for the additional simulations were created following Springel et al. (2005) as outlined in the previous paragraph, again only changing the distribution of the stellar contribution to the potential.To add more high- * galaxies, we include one intermediate mass, compact galaxy ( b = 0.75 * ,  b = 1 kpc), an additional  b = 0.9 * galaxy with  b = 1.5 kpc, and two pure spheroids ( b =  * ) with scale radii of 1 and 2 kpc respectively.We summarise the properties of our initial conditions in Table 1.
The simulations are performed with the moving-mesh code Arepo (Springel 2010;Weinberger et al. 2020), which treats stars and dark matter as Lagrangian particles, while the hydrodynamics are solved on an unstructured mesh in form of a Voronoi tessellation.Gravity is solved with a tree-based scheme and we use an adaptive gravitational softening approach (Price & Monaghan 2007) for optimal gravitational resolution.With a mass resolution of ∼ 104 M ⊙ for gas and stars, and a mass resolution of ∼ 10 5 M ⊙ for the dark matter, we fix the minimum softening length at 12 pc for gas and stars and at 26 pc for dark matter, respectively.This ensures that the gas is well resolved spatially at our star formation density threshold n thresh = 1 cm −3 and higher.
Star formation is modelled using the Katz (1992) parametrisation for the volumetric SFR density  SFR =  ff  −1 ff , where  ff is the SFE per free-fall time,  the mass volume density and  ff = √︁ 3/32  the free-fall time of the gas.We introduce a dependence on the gas dynamics, using the virial parameter,  vir , dependent SFE of Padoan et al. (2012Padoan et al. ( , 2017)): This dynamics-dependent SFE reproduces the observed suppressed SFRs of ETGs (Gensior et al. 2020;Gensior & Kruĳssen 2021;Kretschmer & Teyssier 2020), but converges to the Kennicutt (1998) relation for late-type or gas-rich galaxies (e.g.Semenov et al. 2016;Gensior & Kruĳssen 2021).
The virial parameter is calculated from a cloud-like overdensity around each star-forming gas cell using the algorithm presented in Gensior et al. (2020).Apart from a minimum density threshold of n thresh = 1 cm −3 , the gas must also be cooler than 10 3 K to be eligible for star formation.We model the thermal state of the gas using the Grackle chemistry and cooling library 4 (Smith et al. 2017) with the six-species chemistry network enabled.Tabulated atomic metal line cooling at solar metallicity and heating from the interstellar radiation field using the Haardt & Madau (2012) UV-background are included as well.
We consider feedback from massive stars in the form of type II supernovae (SNe).These are modelled using the mechanical feedback formulation (Kimm & Cen 2014;Hopkins et al. 2014): numerically the energy, mass and metals injected by each SN are distributed to the 32 nearest neighbours using a smooth particle hydrodynamics kernel.We assume one SN per 100 M ⊙ of stellar mass formed (Chabrier 2003;Leitherer et al. 2014), where each SN injects 10 51 erg s −1 of energy and ejects 10.5 M ⊙ of mass and 2 M ⊙ of metals into the surrounding medium.The SNe are detonated with a delay time of 4 Myr, towards the upper end of the range derived from recent observations of feedback disruption (Kruĳssen et al. 2019;Chevance et al. 2020a,b).

Obtaining gas surface density maps
A total gas mass surface density (Σ gas ) map of each simulation snapshot is created using the ray-tracing functionality within AREPO.The primary focus of this work is to analyse the circumnuclear gas discs of galaxies, therefore we create face-on maps that focus on the central 2.25 kpc of each simulated galaxy.This is comparable in extent to the larger WISDOM maps, e.g. that of NGC 0383 (North et al. 2019, see also Davis et al. 2022).Each gas surface density map is computed by integrating the densities of all the gas cells that fall within the line of sight of each pixel, from 1 kpc below to 1 kpc above the galactic mid-plane.The pixel size of our original maps is 1.2 pc per pixel, comparable to the size of the densest Voronoi cells.
The size is chosen to ensure that all structures are sufficiently spatially resolved.To make this more comparable to the observations, we then smooth each map with a Gaussian kernel of full width at half maximum (FWHM) of 43 pc, equal to the geometric average of the NGC 0383 synthesised beam.Following this smoothing, each map is regridded to a coarser resolution of 10.8 pc per pixel so that the pixel grid nyquist samples the convolution kernel, while conserving surface density.
Although we use NGC 0383 as a point of reference, we tested the spatial resolution dependence of our results using different kernel sizes.We found that this size does not qualitatively affect our conclusions, because fits to the power spectra are restricted to spatial scales larger than the kernel size and we focus on the best-fitting power-law slopes during our main analysis (see Section 2.3 for more details on how we compute the power spectra, and Appendix A2 for details of the resolution test).While the surface density maps of the simulations include contributions from all gas cells within the field of view, which will predominantly be cold and dense, the WISDOM observations only map molecular gas.Hence, we tested how our results are affected by only including gas above different mass surface density thresholds.Using a threshold of 10 M ⊙ pc −2 , above which the gas is expected to be predominantly molecular (e.g.Krumholz et al. 2009), does not alter our conclusions.Additionally, we created mock radio interferometric observations of the simulations using the KINematic Molecular Simulation (KinMS) 5 tool (Davis et al. 2013b(Davis et al. , 2020a)).Similarly to the resolution test, we did not find any qualitative difference between the results of the Arepo projections and the KinMS maps.Therefore, we use the gas mass surface density projections as described above and refer the interested reader to Appendices A1 and B1 for a more quantitative discussion of these tests.
To ensure that we average over fluctuations due to cloud-scale rapid cycling of gas (e.g.Kruĳssen et al. 2019;Chevance et al. 2020b), we generate 8 surface density maps per simulated galaxy, starting from the 300 Myr snapshot, when the gas disc has settled into equilibrium.Subsequently a map is generated for each of seven additional snapshots separated by 100 Myr, approximately a galactic dynamical time, until the end of the simulation at 1 Gyr. Figure 1 shows the gas mass surface density map of the central region of each simulated galaxy, 600 Myr after the beginning of the simulation.One can clearly see the aforementioned differences of ISM structure: the maps of the disc-dominated galaxies (top-left panels) are dominated by dense clumps and filamentary structures, whereas the surface density maps of the bulge-dominated galaxies (bottom-right panels) appear very smooth, apart from the very bright central mass concentrations.

Observational Data
The molecular gas mass surface density (obtained from the zerothmoment) maps of twelve galaxies from the WISDOM project constitute the observational component of our comparison, making it the largest sample of CO extragalactic power spectra studied to date.This sample comprises five late-type and six early-type galaxies as well as one dwarf spheroidal galaxy, with a range of central stellar mass surface densities.We calculate  * following Equation 1, using stellar masses derived from dynamical models where available, or from the  = 0 Multiwavelength Galaxy Synthesis survey (z0MGS; Leroy et al. 2019) otherwise.The -band effective radius is used as  e in Equation 1 and taken from the Two Micron All Sky Survey 5 https://kinms.space/(2MASS; Skrutskie et al. 2006) whenever possible.Table 2 lists the selected galaxies, their associated properties and observational references.We assume uncertainties of 0.1 dex for gas and stellar mass, 1.5 ′′ for  e and 0.2 dex for the SFR.
All selected WISDOM galaxies were observed with the Atacama Large Millimetre/submillimetre Array (ALMA) in both a compact and an extended configuration, such that the data are sensitive to emission on spatial scales up to ≈11 arcsec in all cases.For ten galaxies the  = 2→1 transition of 12 CO was observed.For NGC 4429 and NGC 4826 the 12 CO(3-2) line was observed instead.The raw ALMA data for each galaxy were calibrated using the standard ALMA pipeline, as provided by ALMA regional centre staff.The Common Astronomy Software Applications (Casa;McMullin et al. 2007) package was then used to create the final data cube, from which the zerothmoment map was generated using a masked moment technique (e.g.Dame 2011).A more detailed description of the data reduction and calibration can be found in the relevant, associated WISDOM project papers (e.g.Davis et al. 2017Davis et al. , 2018Davis et al. , 2020b;;North et al. 2019North et al. , 2021;;Smith et al. 2019Smith et al. , 2021;;Liu et al. 2021Liu et al. , 2022;;Lu et al. 2022, andparticularly Davis et al. 2022).We use this zeroth-moment map to estimate the total molecular gas mass ( H2 ) listed in Table 2. To convert from brightness temperatures to surface densities we assume a CO-to-H 2 conversion factor  CO = 4.36 M ⊙ (K km s −1 ) −1 pc −2 (corresponding to  CO = 2 × 10 20 cm −2 (K km s −1 ) −1 , Dickman et al. 1986), a CO(2-1)/CO(1-0) intensity ratio (in temperature units) of 0.7 and a CO(3-2)/CO(1-0) intensity ratio of 0.3 (Leroy et al. 2021a).This is the same procedure as used in the companion paper on ISM morphologies Davis et al. 2022.Table 3 details the properties such as synthesised beam and pixel sizes for all galaxies in the sample.
We show the moment zero maps of all selected galaxies in Figure 2. Given that the primary aim of the WISDOM project is to dynamically measure SMBH masses in nearby galaxies, the resultant set of observed circumnuclear gas discs is quite heterogeneous.The molecular gas surface density sensitivity of the observations ranges from approximately 4 to 98 M ⊙ pc −2 , while the spatial resolution varies between approximately 1 to 80 pc.The extent and morphology of the circumnuclear gas reservoir also differ significantly across the twelve galaxies.As quantified with non-parametric morphological indicatiors by Davis et al. (2022), the circumnuclear gas discs shown in Figure 2 largely follow the prediction of Gensior et al. (2020), that a bulge-dominated galaxy with high central stellar mass surface density will have a smooth central gas reservoir.To quantify this further, we hereafter turn to computing the mass surface density power spectra for all simulation snapshots and selected WISDOM galaxies.

Computing power spectra
The power spectra are computed using the python package TurbuStat6 (Koch et al. 2020).For each map, the two-dimensional power spectrum is obtained by multiplying the Fourier transform of the mass surface density map by its complex conjugate.In this work, we focus on the one-dimensional (1D) power spectrum, which is generated by azimuthally averaging the two-dimensional power spectrum of each map.To simplify the comparison between simulations and observations, we fit a single power-law to each power spectrum: where  is the amplitude of the power-law,  the wavenumber, and  the power-law (i.e.power spectrum) slope.As discussed in Section 1, a broken power-law occasionally provides a better fit to the data (e.g.Dutta et al. 2009a;Grisdale et al. 2017), and we will discuss individual deviations from a single power-law briefly in the Section 3.2 and in more detail in Appendix D.
Individual bright emission regions can dominate a power spectrum on small scales, to the extent that they can induce a bump and break in the power spectrum (e.g.Willett et al. 2005;Koch et al. 2020).In our simulated galaxy sample, we observe this behaviour in all of our bulge-dominated galaxies, due to the extremely dense region at the centre.There, gas accumulates continuously because star formation is suppressed and the simulations do not include AGN feedback, which could remove (some of) the innermost gas (Gensior et al. 2020).Therefore, we mask the centre of each simulated galaxy using an inverted cosine bell window function and refer the reader to Appendix B2 for more details.To keep the analysis as consistent as possible, we similarly mask the centre of each WISDOM zerothmoment map, where the extent of the central region masked depends on the size of the map.To avoid contamination of the power spectra by Gibbs ringing (stripes in the two-dimensional power spectrum; artefacts caused by a sharp cut-off of the emission), we additionally use a Tukey filter with  = 0.1 to taper the edges of simulated and observed gas maps.The  value of the Tukey filter tapering was determined based on a visual inspection of the two-dimensional power spectrum, where  = 0.1 is the most conservative value that   93  9, 19, 24, 25  NGC 4501 3.3 (Spiral) 15.3 [4]  8.9 [9]  11.0 [14] 58 [15]  8.94  0.43 [14]  58.7  135 9  NGC 4826 2.2 (Spiral) 7.36 [4]  7.9 [9]  10.2 [14] 69 [15]  8.62  -0.71 [14]  59.5  100 9  NGC 5806 3.2 (Spiral) 21.4 [4]  9.0 [9]  10.6 [14] 30 [15]  8 70.0 115 9, 27 Notes: For each galaxy, column 1 lists its name and column 2 its morphological type according to the HyperLEDA database (Makarov et al. 2014).For the following properties, a reference to the source is indicated in square brackets and listed at the end of this description.The galaxy distance is listed in column 3, the molecular gas mass, estimated from the moment zero map as described in the text, in column 4. The total stellar mass and the -band effective radius are listed in columns 5 and 6, respectively.Column 7 lists the central stellar mass surface density, calculated from the stellar mass and effective radius as described in the text.fully removes the effects of Gibbs ringing.Changing the value of the Tukey filter  to a slightly larger value of 0.15 does not significantly affect the best-fitting power-law slopes.When varying ,  changes by ≲0.05, i.e. within the uncertainties on the power spectrum fit.
The fit to each power spectrum is limited by the beam size at small spatial scales and by the spatial scale within which all flux has been recovered at large scales.Specifically, here, we limit the fit to spatial scales larger than 3 FWHM Gauss / √︁ 8 log 2, where FWHM Gauss is the FWHM of the Gaussian smoothing kernel used on the simulations, or the geometric average of the beam for the WISDOM observations which has not been circularised to maintain resolution.This means we avoid contamination of the power spectrum by beam effects 7 and noise from pixelation (Koch et al. 2020).At large scales the extent of the fit is determined by the extent of the map (simulations and observations where the largest scale is smaller than 11 ′′ ), or by the spatial scale corresponding to 11 ′′ .The exact fit range is stated in Table 4 for the simulations and in Table 5 for the observations.

Simulations
We first show example power spectra, with the power-law fits of Equation 3 overlaid, for the suite of simulations at 600 Myr in Figure 3.In each panel (of fixed initial bulge mass, increasing from the left to the right panel), the power spectra of different initial bulge scale radii are shown as solid lines, colour-coded by their central stellar mass surface densities.On each power spectrum the power-law fit is overlaid as a black dashed line.The shaded area around each solid line indicates the standard deviation of the azimuthal power in the two-dimensional power spectrum at this spatial scale.The vertical dashed line indicates the smallest spatial scale considered for the fit.As detailed in Section 2.3, smaller scales are too affected by the beam to be included in the fit.
The panels in Figure 3 represent the different bulge mass fractions of 30 (left), 60 (centre left), 90 (centre right) and 100 (right) per cent.All simulated galaxies are well fitted by a single power-law within 7 We have verified that including a beam response term in the power spectrum model and fitting within the same limits does not significantly affect the slope of the best-fitting power-law, compared to the fit to Equation 3. 2.52 ± 0.05 0.14 83 -1500 Notes: For each simulation, column 1 lists the name, column 2 lists the time-averaged best-fitting slope and its uncertainty, column 3 lists the timeaveraged uncertainties on the slope listed in column 2, and column 4 lists the spatial scales across which a power-law is fitted to the 1D power spectrum.
the uncertainties of the azimuthal averages.In each panel, the power spectrum slopes increase from the more extended ( b = 3 kpc) toward the more compact ( b = 1 kpc) bulge.The slopes also increase from the left panel to right panel, demonstrating that a relatively more massive bulge (as the bulge-to-disc ratio also changes with increasing bulge mass) with the same scale radius results in a steeper power spectrum.This steepening of the power-law is expected if fragmentation is progressively suppressed (Grisdale et al. 2017).The trend shown in Figure 3 quantitatively confirms the qualitative impression from Figure 1 and the results of the non-parametric morphological indicator study (Davis et al. 2022) that galaxies with more compact and massive bulges, and thus with higher central stellar mass surface densities  * , host smoother central gas reservoirs in which fragmentation is suppressed, compared to the disc-dominated galaxies that have clumpier and more sub-structured ISM.
The power spectrum of the bulge-less galaxy is a mild outlier from the aforementioned trend, in that it has a slope that is steeper than that of two low- * , bulge-containing galaxy, even when accounting for the uncertainty of the fit.However, it is also the galaxy with the highest SFR of the simulated sample (Gensior et al. 2020), and therefore has the most SN feedback.This feedback contributes to the gas turbulence (e.g.Krumholz et al. 2018), leading to a steeper slope (see also e.g.Walker et al. 2014).
Table 4 summarises the time-averaged results of the power spectrum fits to all simulated galaxies.The trend observed in Figure 3, that galaxies with higher  * have steeper power spectra, also holds true when averaged across 700 Myr.In addition to the time-averaged best-fitting slope and its uncertainty, we also list the time-averaged uncertainty on the slope fit for each simulation.Table 4 also confirms that the average power spectrum slope of the noB run agrees with that of the runs B_M30_R3, B_M30_R2 and B_M60_R3 within their variations over time and/or the uncertainties on the fits.
To be more quantitative in linking the galactic gravitational potential to the power spectrum, we plot the power spectrum slopes as a function of the central stellar mass surface densities of all simulated galaxies (calculated as described in Section 2.1) in Fig- ure 4.There is a clear trend of increasing  with increasing  * for log( * /M ⊙ kpc −2 ) ≤ 8.75, which flattens at log( * /M ⊙ kpc −2 ) > 8.75.However, it is unclear whether this flattening is physical in origin, or whether it is due to the limited spatial resolution of the simu- lations, or is caused by other parameter choices.While the power-law fit is restricted to spatial scales larger than the minimum gravitational softening, the sub-structure on smaller scales might be affected by the gravitational softening choice.Small scale fragmentation might be even more suppressed in the highest  * galaxies which would lead to a steeper power spectrum, and thus a continuous trend of increasing  with  * , however higher resolution simulations would be required to test this.Interestingly, the flattening of the - * trend occurs at a slope of ∼2.6, which is close to that expected of two-dimensional Kolmogorov (1941) turbulence (8/3).If physical, this flattening could therefore indicate that high shear eventually drives incompressible turbulence.We use Spearman's rank correlation coefficient to assess the strength and statistical significance of this correlation.To take into account the uncertainties on the time-averaged slopes, we perform 1000 Monte-Carlo simulations assuming that the power spectrum slope of each data point is well described by a Gaussian distribution.This allows us to quote both Spearman's rank correlation coefficient and p-value with uncertainties, as given by the average and the 16thto-84th percentile of the Monte-Carlo simulation.A Spearman rank correlation coefficient of  = 0.81 +0.05 −0.05 with  = 0.001 +0.001 −0.001 confirms that the correlation between  and central stellar mass surface density is strong and statistically significant.This directly links the galactic gravitational potential to the suppression of fragmentation via the power spectrum slope.It also confirms that shear from a high- * potential does not only drive turbulence that suppresses star formation, but it also tears apart clouds and suppresses fragmentation.This is a strong prediction from the simulations, to which we can compare the WISDOM observations.

Observations
Figure 5 shows the power spectra of the WISDOM galaxies considered in this paper.All spectra are plotted on the same spatial scale, both to make comparison easier and to highlight the different spatial scales under consideration (that are function of the galaxy distance, map size and synthesised beam).In each panel, the navy dashed line shows the best-fitting power-law to the power spectrum, an inset shows the centre-masked zeroth-moment map from which the power spectrum was computed and the black vertical dashed line indicates the smallest spatial scale considered for the fit.We summarise the results of the 1D power spectrum fits of the WISDOM galaxies in Table 5.
Only half of the galaxy power spectra (FRL 0049, NGC 0404, NGC 1387, NGC 1574, NGC 4429 and NGC 4826) are well described by a single power-law.The other spectra show very pronounced bumps and/or breaks.In NGC 0383, NGC 0524, NGC 4501, NGC 5806, and NGC 6958 this is likely caused by ring structures.NGC 6753 possesses a very bright central region within a very low density region, itself surrounded by a low-density gas disc with weak spiral arms.These likely give rise to the complicated shape of the power spectrum, with multiple break points.
To assess the validity of our single power-law fits, we recompute .Gas mass surface density power spectra of the WISDOM galaxies considered in this paper (solid teal lines, with uncertainties as shading), plotted with unique spatial scale for ease of comparison.In each panel, the best-fitting power-law is shown as navy dashed line, the vertical black dashed line indicates the smallest spatial scale considered for the fit and the legend lists the best-fitting power-law slope with its 1 uncertainty.The inset shows the centre-masked moment zero map of the galaxy, from which the power spectrum was computed.42 -700 Notes: For each galaxy, column 1 lists the name, column 2 lists the best-fitting slope and the 1 uncertainty on the fit.Column 3 lists the spatial scales across which a powerlaw was fitted to the 1D power spectrum.the fit to the power spectra of the observations, allowing for a break, i.e. a dual power-law fit using the turbustat fitting routine and a segmented linear fit.Comparing the Bayesian information criterion of the fits with and without a break yields only three galaxies, NGC 4501, NGC 4826 and NGC 5806, whose power spectra are better fit by broken power-laws.A break in a surface mass density power spectrum is thought to indicate the scale height of the gas disc, as turbulence transitions from three-dimensional turbulence on small scales to two-dimensional turbulence on large scales (e.g.Dutta et al. 2009b).The break scale of the NGC 4826 power-law is at 27 pc, which would imply an extremely thin molecular gas disc.Yim et al. (2014) measured molecular gas scale heights of 31 ± 8 and 33 ± 8 pc for two other galaxies, NGC 5907 and NGC 4565, respectively, but CO scale heights of 100-200 pc are more common (Wilson et al. 2019).However, 27 pc approximately matches the size of the dense clouds seen in the zeroth-moment map, which may cause the break (Körtgen et al. 2021).The breaks in NGC 4501 and NGC 5806 occur at 239 pc and 191 pc, respectively, and could more plausibly be linked to the scale height of their molecular gas discs.Alternatively, the fits might be affected by the prominent ring structures.However, even the best-fitting single power-laws tend to be within the uncertainties associated with the 1D gas mass surface density power spectrum of these galaxies.Therefore, we use the best-fitting slope of the single power-law for all observed galaxies for consistency and we refer the reader to Appendix D for an in-depth discussion of the multicomponent power-law fits.In Figure 6 we overplot the WISDOM galaxies' best-fitting power spectrum slopes as a function of their central stellar mass surface densities on the simulation data.Contrary to the clear trend present in the simulations, the slopes of the observed galaxies have more scatter with respect to  * .This is a clear departure from the prediction of the simulations.There is neither a trend of power spectrum slope with central stellar mass surface density for the whole dataset of WISDOM galaxies ( = −0.14+0.24  −0.24 ,  = 0.55 +0.31 −0.34 ), nor for the sub sets of spiral galaxies ( = 0.02 +0.48  −0.52 ,  = 0.57 +0.30 −0.39 ) or ETGs ( = −0.09+0.33  −0.31 ,  = 0.57 +0.29 −0.31 ).The lack of a correlation with  * could in part stem from the unsuitability of some power spectra to being fit by a single power-law.Alternatively, the observations might follow the trend of the simulations (at log( * /M ⊙ kpc −2 ) > 8.75), but the measurement uncertainties on  and  * as well as potential scatter from physical effects not included in the simulations could obscure it.We will discuss the discrepancies between power spectra of the simulations and WISDOM galaxies in more depth in Section 4.2.However, before doing so, we wish to investigate whether the power spectrum slopes of the WISDOM galaxies are correlated with other quantities.

Correlations with the slopes of the WISDOM power spectra?
In an attempt to identify the driver of turbulence in the WISDOM galaxies, we have checked for correlations between the power spectrum slope and a large number of quantities.Some are observational (beam size, extent of the fit, lower spatial scale limit of the fit, sensitivity of observations) while others are galaxy properties (gas mass, mean central molecular gas mass surface density, stellar mass, gas fraction central gas-to-stellar mass surface density ratio, SFR, specific SFR, star formation efficiency, Gini, Asymmetry, Smoothness and the stellar-to-gas Toomre (1964) Q ratio).Table C1 lists all Spearman rank correlation coefficients and corresponding p-values.
If we define  ≤ 0.1 as being statistically significant, we find a single statistically-significant correlation (with mean central molecular gas mass surface density) for the observations.This statistically-significant correlation, as well as some other quantities of interest, are discussed in more detail below.We focus on the observations here, because the simulations only differ in the underlying stellar potential by construction.They are initialised with the same gas fraction and the same radial gas extent, all of which only evolve very slightly, thus resulting in a very narrow range of mean central molecular gas mass surface densities (see also Figure 4 in Davis et al. 2022).Differences in the simulated galaxies' (s)SFR are driven by changes in  * (see Figure 15 in Gensior et al. 2020), therefore any trend of the physical quantities investigated for the observations, if present in the simulations, will be a result of the correlation with the central stellar mass surface density discussed in Section 3.1.

Sensitivity of the observations
Koch et al. ( 2020) highlighted that observational effects can influence a power spectrum, even on spatial scales larger than the nominal spatial resolution (geometric beam size; see also Körtgen et al. 2021).With sensitivities (Σ gas,thresh ) that vary by more than an order of magnitude across WISDOM galaxies, from 4 (NGC 4501) to 97 M ⊙ pc −2 (NGC 6958), it is worth examining whether the sensitivities of the WISDOM galaxies affect the derived power spectra.The top left panel of Figure 7 reveals a putative trend between the molecular gas mass surface density sensitivity of the dataset and the best-fitting power-law slope of the resultant power spectrum.The power spectrum appears to steepen with worsening (i.e. higher surface density) sensitivity.However, the  − Σ gas,thresh trend has a Spearman rank correlation coefficient  = 0.52 +0.16 −0.17 and p-value  = 0.14 +0.14 −0.13 , so it is not statistically significant.To nevertheless test if any of our results are affected by this possible bias, we have performed the observational analysis again with all zeroth-moment maps clipped to a molecular gas mass surface density of 70 M ⊙ pc −2 , to mimic a uniform sensitivity.It turns out even this crude approximation only has a negligible effect on the correlation between best-fitting powerlaw slope and central stellar mass surface density shown in Figure 6.Despite this, we advise to be mindful of potential effects when comparing the power spectra of a heterogeneous set of observations.

Star formation rate
Star formation is one of the primary drivers of feedback and thus turbulence (e.g.Mac Low & Klessen 2004;Hennebelle & Falgarone 2012).Here we examine whether there is a correlation between the power spectrum slopes and the SFRs of the selected WISDOM galaxies.With our sample comprising ETGs and spirals, we cover several decades in SFR, thus expect to have sufficient dynamic range to identify a correlation, should it exist.However, the bottom left panel of Figure 7 shows that there is no trend between the power spectrum slope and the SFR of a galaxy.This is confirmed by a Spearman rank correlation coefficient of  = 0.35 +0.17 −0.16 with  = 0.31 +0.25 −0.23 .

Gas fraction
Dynamical suppression has a strong dependence on the gas fractions of galaxies (Martig et al. 2013;Gensior & Kruĳssen 2021).Additionally, Gensior & Kruĳssen (2021) demonstrated that the morphology of the circumnuclear gas reservoir can be strongly affected by the gas-to-stellar mass ratio of the galaxy.At fixed (high) central stellar mass surface density, the smooth circumnuclear region of the gas disc decreases in extent until finally the entire ISM is porous and sub-structured for increasing initial cold gas-to-stellar mass ratio   from 1 to 20 per cent.For simplicity we have restricted ourselves in this study to the analysis of the Gensior et al. ( 2020) simulations with a constant initial cold gas-to-stellar mass ratio of 0.05, but the observed galaxies have a range of gas fractions.Therefore, we investigate the molecular gas-to-stellar mass ratios of the WISDOM galaxies in relation to their power spectrum slopes in the top right of Figure 7.However, as evidenced by the Spearman rank correlation coefficient  = 0.40 +0.17 −0.16 with a  = 0.25 +0.21 −0.20 , there is no statistically-significant trend.There is also no secondary trend with  * , in contrast to the prediction from the simulations.

Mean central molecular gas mass surface density
The bottom right panel of Figure 7 shows the power spectrum slopes of the WISDOM galaxies as a function of their central mean molecular gas mass surface density, ΣH 2 .Following Davis et al. (2022), ΣH 2 is the average molecular gas mass surface density within an ellipse of semi-major axis 1 kpc, that is centred on the centre of the galaxy and oriented according to the galaxy's position angle.There is a trend of steeper power spectra with increasing mean central molecular gas mass surface density with a Spearman rank correlation coefficient  = 0.61 +0.13  −0.13 and  = 0.06 +0.05 −0.05 .Therefore, this trend is statistically significant for our definition of  < 0.1.Davis et al. (2022) found similar correlations between the non-parametric morphological pa-rameters and ΣH 2 , in that galaxies with higher central gas surface density had a smoother, less sub-structured ISM.A correlation with ΣH 2 suggests that the self-gravity of the gas could play an important role in setting the structure of (and driving turbulence in) the ISM.However, naively one would expect the strong self-gravity to affect the ISM structure in the opposite way, through fragmentation, from which one would expect an anti-correlation between  and ΣH 2 .Therefore, it is unclear whether self-gravity is the dominant driver of turbulence in these galaxies.Davis et al. (2022) analysed the stellarto-gas Toomre (1964) Q ratio to determine whether the stability of the system is dominated by the stellar or the gaseous component.They found that in ≈75 per cent of galaxies the stability of the gas disc is dominated by the stellar component, suggesting that the ΣH 2 trend is a secondary correlation with the gravitational potential.However, we find no trend between the power spectrum slopes and the stellar-togas Q ratios of the galaxies, nor with any other more direct tracers of the potential, implying that while the mean central gas mass surface density might weakly depend on the potential, it seems to be the more relevant quantity for the power spectrum.Another possibility is that the  correlates with ΣH 2 , as a tracer for the turbulent energy per unit area of gas, which depends on the gas surface density (e.g.Krumholz et al. 2018).This potentially hints at gravitational momentum transport or accretion driven turbulence as the drivers of turbulence in the observed galaxies, explaining the correlation with ΣH 2 .However, the differences in the sensitivities of the observations (see Table 3 and Section 4.1.1)may affect the positive correlation between power spectrum slope and mean central molecular gas mass surface density.Therefore, a large(r), more homogeneous sample of observations is required to confirm this trend.The simulations do not show any trend with ΣH 2 due to the small dynamic range of mean central gas mass surface densities they have by construction.
All data points in Figure 7 are colour-coded by their central stellar mass surface density, to see if any secondary correlations with  * exist.However, contrary to expectations from the simulations and the results of Davis et al. (2022), no such secondary correlations exist, implying that the stellar potential has a smaller effect on the power spectra (and possibly driving turbulence) of the observed galaxies, than on the ISM structure as quantified through Gini, Asymmetry and Smoothness.Henshaw et al. (2020) found that analysis of the gas velocity structure, in addition to the gas density, is required to fully understand the density structure in different galactic environments.Therefore, the velocity spectra of the WISDOM galaxies could provide more insight into the power spectrum slope discrepancies of galaxies with similar discs (as per non-parametric morphology indicators, see Davis et al. 2022) and galactic properties -this will be pursued in future works.

Differences between observations and simulations
Figure 6 shows the power spectrum slopes of both simulated and real galaxies as a function of their central stellar mass surface densities.Shear induced by the deep gravitational potentials is the main driver of turbulence in the simulations, as indicated by the trend of steeper power spectra as a function of  * for these.The observations with central stellar mass surface densities in the range 8.3 ≤ log( * /M ⊙ kpc −2 ) ≤ 8.8 follow the simulated data.However, as pointed out in Section 3.2, instead of even steeper power spectra or a flattening of the trend for WISDOM galaxies with larger central stellar mass surface densities, the other eight galaxies are scattered in  −  * space.Nonetheless, they remain broadly consistent with the simulations.This could indicate that the flattening of  to ∼ 2.6 at log( * /M ⊙ kpc −2 ) ≥ 8.75 for the simulated galaxies is physical in nature, and in that case, given the close correspondence with the slope of two-dimensional Kolmogorov (1941) turbulence, that high shear might eventually drive incompressible turbulence.In the remainder of this subsection, we will further discuss the differences between the simulations and observations, that may contribute to a more nuanced scenario and explain the larger scatter of the observations.
It should be kept in mind that the simulations are idealised galaxies in isolated boxes, where the only difference between the individual simulated galaxies is the distribution of the stellar matter.Therefore, their only sources of turbulence are shear, (self-)gravity and SN feedback.The simulated galaxies have not evolved over a Hubble time in a cosmological setting, i.e. they have never accreted gas from the cosmic web, undergone mergers or interactions with other galaxies, and by construction do not contain black holes.In contrast, each of the WISDOM galaxies hosts a SMBH, which is active in some of them (e.g.FRL 0049, NGC 0383 and NGC 6753).To our knowledge, the impact of an active galactic nucleus on the power spectrum of a galaxy has not yet been investigated.The radius of the SMBH sphere of influence tends to be comparable to the lower spatial-scale limit of the power spectrum fit, so the SMBHs direct influence on the potential should have little impact on the fit.However, it remains difficult to estimate the effect of an active SMBH on the power spectrum, as its feedback could drive turbulence in the gas (e.g.Venturi et al. 2021;Tamhane et al. 2022).
Cluster membership and related effects on galaxies can also affect their power spectra.Investigating the power spectrum of the Virgo cluster galaxy NGC 4254, Dutta et al. (2010) found that the power spectrum slope is different in the inner and outer regions of the galaxy (i.e. for different velocity channels), as well as in the velocity-integrated power spectrum.Dutta et al. (2010) attributed this difference to galaxy harassment.As most of the WISDOM galaxies considered here are members of a galaxy group or cluster, this could contribute to the scatter of their power spectrum slopes and the poor quality of the fit for some.It is also important to capture the effects of galaxy interactions and cosmological evolution in simulations, as demonstrated by comparing the the cold gas mass surface density power spectrum slope between observation and simulation of the Small Magellanic Cloud (SMC): With a (large-scale) slope of 1.7, Grisdale et al. (2017) found a much shallower power-law for their simulated SMC than that of Stanimirovic et al. (1999), who measured a slope of 3.04 from Hi observations.Grisdale et al. (2017) argued that this discrepancy is the result of running an isolated galaxy simulation lacking tidal effects from galaxy interactions.Another potentially important source of turbulence is gas accretion both on the larger spatial scale for accretion onto the gas disc of the galaxy, but also on the smaller spatial scale accretion of gas onto molecular clouds (e.g.Klessen & Hennebelle 2010), which also impacts the power spectrum slope (Körtgen et al. 2021).
In short, while there is broad agreement between the simulations and the WISDOM observations, the discrepancies are likely the result of several factors.In particular, the heterogeneity of the observations (e.g.sensitivity, see Section 4.1) could dilute the effect of the stellar potentials on the power spectrum slopes, and there are many effects not captured by the simulations (e.g.AGN feedback).A larger, more homogeneous sample of observations should be able to address some of these observational issues.Similarly, high-resolution cosmological zoom-in simulations should produce more realistic comparisons.

Simulations
Over the past decade, power spectra have been used to assess the necessity for and quality of feedback models in simulations.Most of these simulated galaxies have breaks in their power spectra at scales of several tens of pc to 1 kpc (e.g.Bournaud et al. 2010;Combes et al. 2012;Renaud et al. 2013;Walker et al. 2014;Grisdale et al. 2017; although the power spectrum break in Renaud et al. 2013 is at a much smaller scale of about 1 pc).Unlike these, the power spectra of our simulations are well described by single power-laws, comparable to the two runs with less energetic feedback of Walker et al. (2014).The slopes (at large scales in the case of broken powerlaw fits) range from 0.67 (Renaud et al. 2013) to 4.5 (Pilkington et al. 2011), although they are mostly in the range 1.3 -2.6 (Bournaud et al. 2010;Combes et al. 2012;Walker et al. 2014;Grisdale et al. 2017) for runs with stellar feedback.The latter range is in good agreement with that of our own simulated power spectra (1.8-2.7).A common trend in the aforementioned simulations is a steepening of the power spectra in simulations with (increasing) feedback.This may seem counter-intuitive when compared to our result of steeper power spectra with increasing suppression of fragmentation (and therefore suppression of star formation and subsequent feedback).However, there is a simple explanation for this: both strong feedback and shear from the galactic gravitational potential have the same effect.Sheardriven turbulence shifts the power in the gas from small to large scales by suppressing fragmentation and the associated formation of GMCs.
Similarly, strong feedback shifts the power from small-scale to large scale by destroying clouds more effectively (e.g.Pilkington et al. 2011;Walker et al. 2014).Some evidence of the effects of feedback are also present in our simulated galaxies.Indeed, Table 4 shows that the simulated bulge-less galaxy has a steeper power spectrum than a couple of other simulated low- * galaxies, a direct consequence of its higher SFR.

Observations
A variety of tracers have been used in studies of extragalactic power spectra.However, only M 31 (Koch et al. 2020) and M 33 (Combes et al. 2012;Koch et al. 2020) have been previously analysed in CO.With slopes of 1.59 for M31 and 1.5 (Combes et al. 2012) or 0.91 (Koch et al. 2020) for M 33, these power spectra are much shallower than those measured here for the WISDOM galaxies.If ISM turbulence is truly scale-free, then it should not matter that the length scales considered are vastly different.The M 31 and M 33 data probe scales ranging from hundreds of pc to several tens of kpc, whereas the circumnuclear gas discs of some of our galaxies only span a few hundred parsecs, and even in the most extreme case the largest length scale probed is 3.4 kpc.However, like Koch et al. (2020), we find that the morphology of the molecular gas disc is important to determine its power spectrum slope and shape.With often only a small disc detected, it is difficult to compare8 .Encouragingly Miville-Deschênes et al. ( 2010), and Pingel et al. (2018) found power spectrum slopes of 2.7 when using dust and CO emission to study the Polaris flare and Perseus molecular cloud regions within the Milky Way, in better agreement with the slopes of the WISDOM galaxies than those of Combes et al. (2012) and Koch et al. (2020).
The power spectrum slopes of our galaxies are more similar to those often measured from Hi , which tend to be in the range 1.5 -3.0 (e.g.Dutta et al. 2009a;Zhang et al. 2012;Dutta et al. 2013;Walker et al. 2014;Grisdale et al. 2017;Koch et al. 2020, and references therein).This raises an interesting question, namely whether the Hi power spectra of the WISDOM galaxies would be steeper or shallower than their CO counterparts.Turbulence theory predicts that the power spectrum of the molecular gas should have a steeper slope than that of its atomic counterpart (Romeo et al. 2010), although for M31 and M33 Koch et al. (2020) found that the Hi spectra are steeper.With the extremes of the Hi power spectrum slopes at 0.3 (Dutta et al. 2013) and 4.3 (Zhang et al. 2012), both steeper and shallower Hi power spectra for the WISDOM galaxies would not be without precedent9 .
We do not have evidence of ubiquitous breaks in the power spectra, that would indicate transitions from two-dimensional turbulence in the planes of the galaxies to three-dimensional turbulence at the scale heights of the galaxies (e.g.Dutta et al. 2009b).Instead, morphological features specific to each galaxy (such as rings) show up as bumps or wiggles in the power spectrum.Only three galaxies are better described by broken power-laws (than by single power-laws).For NGC 4826 the break scale is at roughly the size of the brightest clumps in the gas disc, making it more likely that the break in the power spectrum indicates a transition from turbulence in the largescale ISM to that in the largest fragments, in good agreement with the findings of Körtgen et al. (2021).The break scale in NGC 4501 and NGC 5806 may be at around around the disc scale height, but it might also be caused by the effect of the dominant ring on the power spectrum.In general, we find that the morphology of the gas is important to determine the power spectrum shape, in good agreement with Koch et al. (2020).
The lack of a correlation between the power spectrum slopes of the WISDOM galaxies and their SFRs matches the results in the literature.Zhang et al. (2012) found a tentative anti-correlation between SFR and power spectrum slope for dwarf galaxies with an absolute -band magnitude lower than -14.5, but none for the entire sample of dwarf galaxies.Similarly, Dutta et al. (2009a) found tentative evidence for a correlation of power spectrum slope with SFR for five dwarf galaxies, but only scatter for a larger sample of 18 late-type galaxies (Dutta et al. 2013).
The trend of steeper power spectrum slopes with increasing mean central molecular gas mass surface densities for the WISDOM galaxies seems somewhat more at odds with past results.However, while Koch et al. (2020) argued that the shallow power spectrum they measured for M 33 is caused by the central concentration of CO (within the inner 2 kpc), the behaviour of the WISDOM galaxies matches that of the Walker et al. (2014) simulated galaxies.In the latter, both of the galaxies with an exponential Hi mass surface density profile have a power spectrum steeper than that of the galaxy that evolved to have a flat Hi mass surface density profile (from the same initial conditions; Walker et al. 2014).Furthermore, Davis et al. (2022) found that the mean central molecular gas mass surface density is one of the best predictors (next to  * ) for the morphology of the central ISM.
Lastly, as previously highlighted by Koch et al. (2020), we find that the observational parameters are important.The power spectrum slopes have a trend with the sensitivities of our observations, which may also affect the correlation with the central molecular gas mass surface densities (through its impact on the measured molecular gas masses).To further quantify the impact of observational parameters on the measured power spectrum slope, we performed tests with the simulations and refer to Appendices A and B for more detail.Even though the synthesised beam sizes (zeroth-moment map spatial resolutions), inclinations and mass surface density thresholds do not qualitatively affect the trend with the central stellar mass surface densities of the simulated galaxies' power spectra, they do make small quantitative differences.Similarly, the discrepancy of the power spectrum slopes when smoothing the simulated maps with Gaussian beams or mock interferometric beams highlights how sensitive the power spectrum slope is to how the data were obtained and the data analysis methods.This makes it doubly difficult to compare to literature results, and draw sound conclusions, as some discrepancies might be purely methodology dependent (see also Zhang et al. 2012;Körtgen et al. 2021).

SUMMARY AND CONCLUSIONS
In this paper, we have analysed the cold gas mass surface density power spectra of the circumnuclear gas discs of 14 isolated, simulated galaxies from Gensior et al. (2020) and 12 galaxies from the WIS-DOM project (e.g.Onishi et al. 2017), making it the largest sample of CO extragalactic power spectra studied to date.The morphologies of the simulated galaxies range from pure disc to spheroid, with a fixed initial gas-to-stellar mass ratio of 5 per cent, and the circumnuclear gas discs become increasingly smooth with increasing central stellar mass surface density ( * ), i.e. as the bulge dominance increases.The WISDOM galaxies comprise a mix of late-and early-type galaxies, some of which have very sub-structured, and others very smooth gas reservoirs.We computed the azimuthally-averaged, 1D power spectra of the gas mass surface densities, which are sensitive to the turbulent forcing of the gas.Our main findings are summarised below.
(i) There is a strong correlation between the power spectrum slopes () and the central stellar mass surface densities of the simulated galaxies; the power spectra being steeper in galaxies with higher  * (and smoother central gas reservoirs).This confirms that the shear induced by the gravitational potential leads to dynamical suppression of fragmentation.(ii) The flattening of the power spectrum slopes at log( * /M ⊙ kpc −2 ) ≥ 8.75 (to a constant ∼2.6) could be physical or could be an artefact caused by the limited spatial resolution of the simulations and other parameter choices.If the flattening is real, the close agreement with the slope of two-dimensional Kolmogorov (1941) turbulence (8/3) could indicate that high shear eventually drives incompressible turbulence.(iii) Contrary to the simulations, the observations do not show any trend between power spectrum slope and central stellar mass surface density, whether they are being considered as a whole, or are divided into sub-samples of early-and late-type galaxies.However, the majority of the WISDOM galaxies have log( * /M ⊙ kpc −2 ) ≥ 8.75, and have power spectrum slopes that scatter about the roughly constant slope of the simulations at such high  * .This suggests that the flattening to a constant  could be real and thus physical in origin.The large scatter of the observational slopes likely results from a combination of physical effects not captured by the simulations, the heterogeneity of the observations and the different gas morphologies (e.g.rings and spiral arms) affecting the (quality of the) power spectrum fits.(iv) Of all the observation and galaxy properties we investigated, we only found a statistically significant correlation between the power spectrum slopes and the mean central molecular gas mass surface densities of the WISDOM galaxies.Although this correlation may depend on the sensitivities of the observations, this suggests that the dominant driver of turbulence in the ISM is related to the gas density, while the gravitational potential (of the stars) might play a less important role.
We have thus demonstrated using isolated galaxy simulations that the shear and turbulence induced by a spheroidal component that dominates the gravitational potential can suppress fragmentation in an embedded cold gas disc, resulting in a steep mass surface density power spectrum.The observations of 12 WISDOM galaxies however reveal a more nuanced reality, and the galaxies' power spectrum slopes depend most strongly on the mean central molecular gas mass surface densities.However, the heterogeneity of the observations and the occasionally poor description of a galaxy's power spectrum by a single power-law make it difficult to reach any strong conclusion.Velocity power spectra of the WISDOM galaxies may further elucidate the turbulence mechanisms and explain the different power spectrum slopes of galaxies with similar properties.Power spectra of a larger, more homogeneous sample of observed galaxies are required to firmly confirm whether the tentative correlation we discovered for the WISDOM sample galaxies holds true.In addition, high resolution cosmological (zoom) simulations, that resolve the cold ISM in sufficient detail, could shed more light on mass surface density power spectra from a theoretical perspective.
Table A1.Spearman rank correlation coefficients and p-values of the correlations between power spectrum slope and central stellar mass surface density of the simulated galaxies, when using different Gaussian FWHM for the spatial smoothing kernel applied to the gas mass surface density maps.jections (1.77 -1.89).This is the result of the synthesised KinMS non-normalised interferometric beam smoothing out some of the small-scale structures, thereby steepening the power spectrum slopes.However, we can reproduce the steeper slopes of the low- * simulated galaxies (within the uncertainties) using the original, non-convolved, Arepo ray-tracing maps and convolving them with the KinMS beam kernel.Since both methods yield qualitatively the same results, in this paper we use the Arepo mass surface density projections convolved with a Gaussian kernel, which stay truer to the underlying gas distributions.However, this implies that the power spectra of interferometric observations could be artificially steepened by suppressing power on small scales.

A2 Gaussian FWHM
The beam of observations and the smoothing kernel used to generate the maps of the simulated galaxies can affect the shape of a power spectrum up to a spatial scale equal to a few times the FWHM Gauss (e.g.Grisdale et al. 2017;Koch et al. 2020).When fitting the power spectra, we therefore only consider spatial scales in excess of 3FWHM Gauss / √︁ 8 log 2, which should mitigate this effect.Nonetheless, we test the robustness of our results to changes in Gaussian smoothing kernel size here, using the simulated galaxies.The effects should be similar when varying the synthesised beam size of a real observation.We vary the FWHM Gauss we apply to the gas mass surface density projections from 20 to 70 pc, in 10 pc increments and subsequently re-grid each map such that FWHM Gauss is always sampled with approximately 4 pixels.When fitting the power spectra, we account for the differences of the Gaussian FWHM by adjusting the lower spatial limit of the fits accordingly.
Figure A3 shows the time-averaged power spectrum slopes of the simulated galaxies as a function of their central stellar mass surface densities, colour-coded by the Gaussian FWHM sizes.The error bars denote the uncertainty on the time-averaged best-fitting power-law slopes, while the grey-shaded regions indicate the average uncertainty of the individual power-law fits.There are some changes of the power spectrum slope with varying FWHM Gauss size, however these are mostly over slope ranges comparable to the uncertainty of the powerlaw fits and there is no clear pattern of how the FWHM Gauss size affects the power spectrum slope.
We list the Spearman rank correlation coefficient and the probability that the trend of each Gaussian kernel data set is statistically significant in Table A1.For all FWHM Gauss considered here, the correlation between power spectrum slope and central stellar mass surface density remains strong and statistically significant.In short, although there are some small quantitative changes of power spectrum slope as a function of the Gaussian smoothing kernel size, these differences are only of the order of the uncertainties of the fits.Qualitatively, at a given  * , there is good agreement between all the power spectra.

A3 Inclination
All selected WISDOM galaxies have inclinations () ranging from 20 to 70 • .The power spectra are computed from (molecular) gas mass surface density maps, and a large part of the analysis focusses on the central stellar mass surface densities of the galaxies, i.e. using projected quantities.Therefore, the inclination of the galaxies could affect both the 1D power spectra (and thus ), as well as  * .Here, we systematically test the impact of inclinations ranging from 0 to 90 • , in 10 • increments, on our simulated galaxy results, by repeating the analysis after inclining the simulated galaxies.We use KinMS to produce the inclined total gas mass surface density projections adopting the same set-up described in Appendix A1.We also recompute the stellar half-mass radius and  * for each inclined stellar particle distribution.The top panel of Figure A4 shows the time-averaged power spectrum slopes of the simulated galaxies as a function of their central stellar mass surface densities, colour-coded by the inclination of the simulated galaxy disc from which the gas mass surface density map was produced.The error bars denote the uncertainty on the time-averaged best-fitting power-law slopes, while the grey-shaded regions indicate the average uncertainty of the individual power-law fits.We only show the points for every third galaxy at full opac- ity for better readability of the plot.For all galaxies with B/T < 1  * increases with , by up to 0.5 dex for the disc-dominated cases and by ≲0.15 dex for the more massive and compact bulges.The best-fitting power-law slope also changes with inclination.However, apart from the extreme,  = 90 • case, the slopes are generally consistent with each other within the uncertainties on the time-averaged slopes, and they are always within the uncertainties of the individual power-law fits.The correlation between  and  * is preserved for  ≤ 80 for  = 0 • ), due to the increase of  * with .By contrast, when the simulated galaxies are viewed edge-on, all power spectra are very steep.For low- * galaxies, some steepening already occurs at  ≥ 70 • .Thus, caution is advised when considering the power spectra of galaxies whose gas reservoirs have inclinations above 70 • .This is in good agreement with Grisdale et al. (2017), who compared the power spectra of their simulated galaxies at  = 0 • and  = 40 • , and found that the main difference is the total power in the power spectra, while the power-law slopes change by very little.Since the WISDOM galaxies to which we compare the simulated galaxies span a range of inclinations, we perform a second test where we use bootstrapping to determine median  and  * , as well as their 16th-to-84th percentile variations, randomly drawing an inclination between 0 − 70 • each time.This is shown as the dark red crosses in the bottom panel of Figure A4.The resultant trend of increasing  with  * (at log( * /M ⊙ kpc −2 ) ≤ 8.75) is more pronounced for the bootstrapped data, and has a Spearman rank correlation coefficient  = 0.83 +0.05 −0.05 with a p-value  = 0.0005 +0.009 −0.00002 .This further confirms that the power-law slope and the trend the simulations show with  * is not strongly affected for  ≲ 70 • .

APPENDIX B: SIMULATION PARAMETER TESTS B1 Surface density threshold
We use the total gas mass surface density projections of the simulated galaxies' gas discs to compute the power spectrum in the main analysis, as we do not use a chemical model that models molecular gas or tracks CO specifically.In contrast, we use molecular gas mass surface density maps to compute power spectra for the observed galaxies.Here, we therefore quantify how the best-fitting power-law slopes of the power spectra are affected by using total gas mass surface density maps clipped to surface density thresholds (Σ thresh ) of up to 100 M ⊙ pc −2 , to test how robust the results of the main analysis remain when using a density threshold that mimics molecular gas.We use the Arepo-generated total gas mass surface density maps, masking all surface densities smaller than our threshold prior to computing the power spectra.
Figure B1 shows the time-averaged power spectrum slopes of the simulated galaxies as a function of their central stellar mass surface densities, colour-coded by the Σ thresh of the gas mass surface density maps from which the power spectra were computed.The error bars denote the uncertainty on the time-averaged best-fitting power-law slopes, while the grey-shaded regions indicate the average uncertainty of the individual power-law fits.For Σ thresh < 10 M ⊙ pc −2 , the power-law slopes remain unaffected, the data points systematically fall on top of each other.For Σ thresh = 10 M ⊙ pc −2 , the power-law slopes are in good agreement with the ones of the unclipped maps, in some cases within the uncertainties of the time-averaged slopes and always within the average uncertainties of the individual fits.Notable differences are however present for Σ thresh ≥ 30 M ⊙ pc −2 , i.e. for the darker blue and purple diamonds in Figure B1.At those high Σ thresh , the power spectra become much shallower as the threshold is increased for all simulated galaxies.This is the result of ignoring low-density regions, or conversely focussing on increasingly small (very) high-density peaks.The latter are preferentially found in the low- * galaxies, which have clumpier ISM with starker density contrasts.Therefore, the differences of the power-law slopes are larger for the high- * galaxies.Having said that, these high surface density thresholds are somewhat unphysical for our simulations, as the gas that remains looks like a noise map in each case.Furthermore, most observations of CO in galaxies now routinly reach surface densities smaller than 10 M ⊙ pc −2 .As before, we list the Spearman rank correlation coefficients and the p-values of the correlations for different thresholds in Table B1.In summary, our results are robust even when introducing a reasonable gas mass surface density threshold that mimics (high-density) molecular gas, such as Σ thresh ≈ 10 M ⊙ pc −2 .

B2 Centre masking
Lastly, we elaborate on the need for masking the central regions of the gas mass surface density maps of the simulated galaxies prior to our computation of the power spectrum.Figure B2 provides an example of just how extreme the effect of the very dense central regions in bulge-dominated galaxies can be, by contrasting the power spectra of the most disc-dominated (noB) and most bulge-dominated (B_M100_R1) simulated galaxies without central masking (colour) Figures D2 and D3 show the power spectrum slopes of either the best-fitting single power-laws or large scale slopes of the best-fitting broken power-laws of the WISDOM galaxies as a function of  * (Figure D2, cf. Figure 6), surface density sensitivity, molecular gasto-stellar mass ratio, star formation rate and average central molecular gas surface density (Figure D3, cf. Figure 7).The discrepancies between the trend between  and  * of the simulated and observed .Gas mass surface density power spectra of the three WISDOM galaxies NGC 4501, NGC 4826 and NGC 5806 (solid teal lines, with uncertainties as shading), plotted with unique spatial scale for ease of comparison.In each panel, the best-fitting single power-law is shown as navy dashed line, the large scale (small scale) best-fitting broken power-law is shown as a purple solid (dotted) line, the vertical black dashed (purple dotted) line indicates the smallest spatial scale (break scale) considered for the fit and the legend lists the best-fitting (large scale) power-law slopes with their 1 uncertainty.The inset shows the centre-masked moment zero map of the galaxy, from which the power spectrum was computed.
galaxies are larger when considering the large scale slopes of the best-fitting broken power-laws.The large scale slopes are (much) shallower than the slopes of the best-fitting single power-law for the simulated galaxies in the same central stellar mass surface density regime.In contrast to the single power-law lopes, the combination of best-fitting single and large scale slopes shows a statistically significant correlation with the surface density sensitivity of the observations.There are still no trends between the best-fitting slopes and star formation rate or molecular gas-to-stellar mass fraction.However, while there still exists a putative trend of increasing power-law slope with central stellar surface density, it is no longer statistically significant (Spearman rank correlation coefficient  = 0.45 +0.10  −0.10 with a p-value  = 0.16 +0.10  −0.10 ) when considering the large scale slopes of the best-fitting broken power-laws for NGC 4501, NGC 4826 and NGC 5806.This paper has been typeset from a T E X/L A T E X file prepared by the author.

Figure 1 .
Figure 1.Gaussian-smoothed and regridded mass surface density projection of the central gas reservoir of each simulation considered in this paper, 600 Myr after the start of the simulation.The simulations are ordered by increasing central stellar mass surface density, from the lowest (noB, top-left) to highest (B_M100_R1, bottom-right).Each map has a linear extent of 2.25 kpc and the purple solid circle in the bottom-left corner indicates the size of the Gaussian beam.The central regions of the bulge-dominated galaxies (bottom right) are much smoother and largely devoid of sub-structure, compared to the clumpy ISM found in the central regions of the disc-dominated galaxies (top left).

Figure 2 .
Figure2.Moment zero maps of all WISDOM galaxies considered in this paper.In the top-right corner of each map we show the mass surface density scaling (conversion from original brightness units as described in Section 2.2).The scale bar in the bottom-right corner of each map gives and indication of the size of the visible circumnuclear gas disc.The purple ellipse in the bottom-left corner of each map shows the size of the synthesised beam.A black dashed circle indicates the 11 ′′ large scale limit of the power spectrum fit.

Figure 3 .
Figure 3. Gas mass surface density power spectra (solid lines colour-coded by the central stellar mass surface density of the simulated galaxy, with uncertainties as shading), computed from the snapshots 600 Myr after the start of the simulations.Each panel compares the effect of changing bulge scale radius (1 -3 kpc) at constant bulge mass, with panels from left to right showing simulations with bulges containing 30, 60, 90 and 100 per cent of the initial stellar mass.The one bulge-less galaxy is shown with the leftmost panel and the single  b = 0.75 * bulge is shown in the middle left panel.In each panel, the best-fitting power-laws are shown as black dashed lines, the vertical black dashed line indicates the smallest spatial scale considered for the fit, and the legend lists the best-fitting power-law slopes with their 1 uncertainties.All galaxies have a similar amount of power (arbitrary unit), thus individual power spectra are vertically offset from each other for visibility.

Figure 4 .
Figure 4. Time-averaged slopes () of the power-laws best-fitting the power spectra of the simulated galaxies as a function of the galaxies' central stellar mass surface densities ( * ).Error bars indicate the uncertainties on the timeaveraged power-law slopes.As indicated by the Spearman rank correlation coefficient and p-value listed in the bottom-right corner, there is a statisticallysignificant correlation between the best-fitting slope of the power spectrum and the central stellar mass surface density of a galaxy.
Figure5.Gas mass surface density power spectra of the WISDOM galaxies considered in this paper (solid teal lines, with uncertainties as shading), plotted with unique spatial scale for ease of comparison.In each panel, the best-fitting power-law is shown as navy dashed line, the vertical black dashed line indicates the smallest spatial scale considered for the fit and the legend lists the best-fitting power-law slope with its 1 uncertainty.The inset shows the centre-masked moment zero map of the galaxy, from which the power spectrum was computed.

Figure 6 .
Figure6.Slopes () of the power-laws best-fitting the power spectra of the simulated (teal) and WISDOM (purple) galaxies as a function of their central stellar mass surface densities ( * ).The vertical error bars indicate the uncertainties on the time-averaged best-fitting power-law slopes (simulations) or the 1 uncertainties on the best-fitting power-law slopes (observations).

Figure 7 .
Figure 7. Slopes () of the power-laws best-fitting the power spectra of the WISDOM galaxies as a function of the galaxies' mass surface density sensitivities (Σ gas,thresh , top left), molecular gas fraction (  H 2 , top right), SFR (bottom left) and average molecular gas mass surface density within the central kpc ( ΣH 2 , bottom right).The data points are colour-coded by their central stellar mass surface densities ( * ).Vertical error bars indicate the 1 uncertainties on the best-fitting power-law slopes and horizontal error bars show the uncertainty on the respective quantity.The Spearman rank correlation coefficient and p-value for each correlation are shown in each panel.There only statistically significant correlation is that between power spectrum slope and the average central molecular gas mass surface density of the observations.

RFigure A1 .Figure A2 .
Figure A1.Total gas mass surface density maps generated using the KinMS tool, for the three galaxies with a bulge stellar mass of 0.3 * (see the top row of Figure1), 600 Myr after the start of the simulations.Each map has a linear extent of 2.25 kpc and the solid purple circle in the bottom-left corner indicates the size of the synthesised beam.There is good agreement between the Arepo-generated mass surface density projections and the KinMS maps, although the mock observations with KinMS tend to wash out small-scale structures.

Figure A3 .
Figure A3.Time-averaged slopes () of the power-laws best-fitting the power spectra of the simulated galaxies as a function of the galaxies' central stellar mass surface densities ( * ).The data points are colour-coded by the Gaussian FWHM of the smoothing kernels with which the maps were spatially smoothed prior to computing the power spectra.The error bars indicate the uncertainties on the time-averaged power-law slopes while the grey-shaded regions indicate the average uncertainties of the power-law fits.

Figure A4 .
Figure A4.slopes () of the power-laws best-fitting the power spectra of the simulated galaxies as a function of the galaxies' central stellar mass surface densities ( * ).Top panel: the data points are colour-coded by the inclinations () of the simulated galaxies, only showing every third data point with full opacity.The error bars indicate the the uncertainties on the time-averaged power-law slopes while the grey-shaded regions indicate the average uncertainties of the power-law fits.Bottom panel: The red crosses show the median  and  * for  ranging from 0 -70 • , error bars indicating the 16th-to-84th percentile variation determined via bootstrapping, overplotted on the  = 0 data points (teal diamonds).
Figure D1.Gas mass surface density power spectra of the three WISDOM galaxies NGC 4501, NGC 4826 and NGC 5806 (solid teal lines, with uncertainties as shading), plotted with unique spatial scale for ease of comparison.In each panel, the best-fitting single power-law is shown as navy dashed line, the large scale (small scale) best-fitting broken power-law is shown as a purple solid (dotted) line, the vertical black dashed (purple dotted) line indicates the smallest spatial scale (break scale) considered for the fit and the legend lists the best-fitting (large scale) power-law slopes with their 1 uncertainty.The inset shows the centre-masked moment zero map of the galaxy, from which the power spectrum was computed.

Figure D2 .
Figure D2.Slopes of the power-laws best-fitting the power spectra of the simulated (teal) and WISDOM (purple) galaxies as a function of their central stellar mass surface densities ( * ), including the large scale slopes of the broken power-laws for NGC 4501, NGC 4826 and NGC 5806.The empty symbols show the result of the best-fitting single power-law to the aforementioned galaxies.The vertical error bars indicate the uncertainties on the time-averaged best-fitting power-law slopes (simulations) or the 1 uncertainties on the best-fitting power-law slopes (observations).

Table 1 .
Initial conditions of the simulations.The top ten simulations were presented inGensior et al. (

Table 3 .
Zeroth-moment map-related quantities of the selected WISDOM galaxies.For each galaxy, columns list the galaxy name, synthesised beam size (FWHM), pixel size and the molecular gas mass surface density sensitivity (i.e.root-mean-square noise of the observations).

Table 4 .
Results of the power-law fits to the 1D power spectra of all simulations

Table 5 .
Results of the power-law fits to the 1D power spectra of all selected WISDOM galaxies • , it even becomes stronger for  ≤ 60 ( = 0.90+0.04