The structure and characteristic scales of the HI gas in galactic disks

The spatial distribution of the HI gas in galaxies holds important clues on the physical processes that shape the structure and dynamics of the interstellar medium (ISM). In this work, we quantify the structure of the HI gas in a sample of 33 nearby galaxies taken from the THINGS Survey using the delta-variance spectrum. The THINGS galaxies display a large diversity in their spectra, however, there are a number of recurrent features. In many galaxies, we observe a bump in the spectrum on scales of a few to several hundred pc. We find the characteristic scales associated with the bump to be correlated with galactic SFR for values of the SFR>0.5 M$_{sol}$ yr$^{-1}$ and also with the median size of the HI shells detected in those galaxies. On larger scales, we observe the existence of two self-similar regimes. The first one, on intermediate scales is shallow and the power law that describes this regime has an exponent in the range [0.1-1] with a mean value of 0.55 which is compatible with the density field being generated by supersonic turbulence in the cold phase of the HI gas. The second power law is steeper, with a range of exponents between [0.5-1.5] and a mean value of 1.5. These values are associated with subsonic turbulence which is characteristic of the warm phase of the HI gas. The spatial scale at which the transition between the two regimes occurs is found to be $\approx 0.5 R_{25}$ which is similar to the size of the molecular disk in the THINGS galaxies. Overall, our results suggest that on scales<$0.5 R_{25}$, the structure of the ISM is affected by the effects of supernova explosions. On larger scales (>0.5 $R_{25}$), stellar feedback has no significant impact, and the structure of the ISM is determined by large scale processes that govern the dynamics of the gas in the warm neutral medium such as the flaring of the HI disk and the effects of ram pressure stripping.

It is well established that turbulence, which is ubiquitously observed in all phases of the gas, is one of the primary regulators of the ISM structure and dynamics of local disk galaxies. It is therefore responsible for setting the self-similar behavior of A&A proofs: manuscript no. dv_h1_v3 many of the physical quantities that are used to describe it (e.g., Elmegreen & Scalo 2004). In the warm (T ≈ 10 4 K) neutral medium (WNM), turbulence is transonic or possibly subsonic, while in the cold (T ≈ 100K) neutral medium (CNM), it is supersonic. Turbulent motions in the WNM and CNM phases can be sustained by a variety of instabilities and energy and momentum injection mechanisms, both internal and external to the galaxy. The spatial scales associated with the fastest-growing modes of these instabilities and those associated with direct energy and momentum injection mechanisms can break the self-similarity of the gas. Some of these scales could be detected as characteristic scales in the ISM (Dib et al. 2009;Eden et al. 2020;Dib et al. 2020). Internal processes include stellar feedback from massive stars, that is, ionizing radiation, radiation pressure, stellar winds, and supernova explosions, which impart significant amounts of energy and momentum to the ISM on intermediate scales, that is, ≈ 50 − 1000 pc (e.g., Heiles 1979, Ehlerova & Palous 1996de Avillez & Breitschwerdt 2005;Dib et al. 2006;Hodge & Deshpande 2006;Shetty & Ostriker 2008;Dib et al. 2011Dib et al. ,2013Gent et al. 2013;Agertz et al. 2013;Hony et al. 2015;Suad et al. 2019;Chamandy & Shukurov 2020;Pokhrel et al. 2020, Bacchini et al. 2020. Large-scale gravitational instabilities due to the combined action of gas and stars (Jog & Solomon 1984;Elmegreen 2011;Shadmehri & Khajenabi 2012;Dib et al. 2017;Marchuk 2018;Marchuk & Sotnikova 2018) can also drive turbulence in galactic disks. Other internal mechanisms of the galaxy that can perturb the self-similar nature of the gas and shape its spatial structure include stellar spiral density waves (e.g., Lin & Shu 1966;Guibert 1974;Adler & Westpfahl 1996;Tosaki et al. 2007;Khoperskov & Bertin 2015;Wang et al. 2015), the Parker instability (e.g., Parker 1967;Franco et al. 2002;Hanasz & Lesch 2003;Rodrigues et al. 2016;Mouschovias et al. 2009;Heintz et al. 2020), and the impact of high-velocity clouds on the galactic disk (e.g., Santillán et al. 1999;Boomsma et al. 2008;Heitsch & Putman 2009;Park et al. 2016). External mechanisms can also impart energy and momentum to the gas on galactic scales. They include ram pressure stripping (e.g., Clemens et al. 2000;Marcolini et al. 2003;Vollmer et al. 2004;Freeland et al. 2010;Steyrleithner et al. 2020) and tidal stripping in interacting systems (e.g., Combes et al. 1988;Marziani et al. 2003;Mayer et al. 2006;Holwerda et al. 2013;Lipnicki et al. 2018;Fattahi et al. 2018). In galaxy mergers, galactic disks can experience strong compressions due to tides, and these compressions can significantly affect the structure and dynamical properties of the gas in the interacting galactic disks (e.g., Renaud et al. 2009).
In this paper, we quantify the structure of the H i gas for a number of nearby galaxies using the ∆-variance spectrum (Stutzki et al. 1998;Ossenkopf et al. 2008). The ∆-variance spectrum is another expression of the power spectrum, and it has been employed successfully to characterize the self-similar structure of the gas as well as to uncover the existence of characteristic scales (e.g., Elmegreen et al. 2001;Dib et al. 2020). In addition to quantifying the structure of the H i gas, we aim to relate features that are observed in the ∆-variance spectra to physical processes that may affect the spatial distribution of the gas. In §. 2 we summarize the sample of galaxies we used in this work, which are taken from the THINGS survey of nearby galaxies . The ∆-variance method is discussed in §. 3, and its application to the THINGS maps is presented in §. 4. In §. 5 we interpret our results using simple models and results from a cosmological zoom-in simulation of a star-forming disk galaxy. We also explore the correlations that exist between characteristic scales detected in the ∆-variance spectrum of the galaxies and their star formation rate (SFR). In §. 6 we discuss our results and compare them to previous work, in particular, to results obtained using the identification of H i shells in the THINGS sample of galaxies. In §. 7 we discuss our results and conclude.

Data: The HI Nearby Galaxy Survey
We used the moment-0 (integrated intensity) H i maps from The H i Nearby Galaxy Survey (THINGS; Walter et al. 2008) 1 . THINGS is a homogeneous survey in the 21 cm H i emission line for 34 nearby galaxies. The observations, performed with the NRAO Very Large Array (VLA), have an angular resolution of ≈ 6 ′′ . At the distances of these galaxies (D gal ≈ 2 − 15 Mpc), this corresponds to spatial resolutions of a few to several hundred parsecs. The galaxies were mapped with various configurations, and the integrated H i intensity maps have a total 1024 × 1024 or 2048 × 2048 pixels. Each pixel represents an angular size of 1 ′′ to 1.5 ′′ , depending on the galaxy. The sample of galaxies spans a range of morphological types, metallicities, total H i mass, and star formation rates extending from low-mass, metal-poor, only weakly star-forming dwarf galaxies to metal-rich massive spiral galaxies with high star formation rates. Galaxies in the THINGS survey have a wide range of inclinations, and it is imperative to correct for the effect of inclination in order to minimize the projection effects. We deprojected all galaxies using the inclinations measured by de Blok et al. (2008) using the H i data alone (i H i ). For the few galaxies (NGC 1569, NGC 3077, and NGC 4449) for which no such measurement is reported in de Blok et al. (2008), we used values of the inclination of the optical disk that are reported in the LEDA database (Paturel et al. 2003). For the position angles (PA) needed to deproject the maps, we used the values listed in Walter et al. (2008). The adopted inclinations are listed in Tab. 1. For NGC 3031, we removed the very central nuclear H i ring (the 50 × 50 inner pixels) from the original data. This ring dominates the signal. We also discarded the first 600 pixels in each direction as they are affected by edge effects.

Method: ∆ Variance spectrum
In order to quantify the structure of the H i gas, we used the ∆-variance spectrum method, originally introduced in Stutzki et al. (1998) and Zielinsky et al. (1999). In this work, we used an improved version of the method presented in Ossenkopf et al. (2008) 2 . For a two-dimensional field A(x, y), the ∆-variance on a scale L is defined as being the variance of the convolution of A with a filter function ⊙ L , such that For the filter function, Ossenkopf et al. (2008) recommend the use of a Mexican hat function, which is defined as where the two terms on the right side of Eq. 2 represent the core and the annulus of the Mexican-hat function, respectively, and v is the ratio of their diameters (we used a value of v = 1.5). The maps have a resolution of 1000 × 1000 pixels. The 2D Gaussian functions all have an aspect ratio ( f = σ 1 /σ 2 = 1) and a contrast between the peak of the Gaussian and the mean value in the map of δ c = 3. The standard deviations of the 2D Gaussians are σ 1 = σ 2 = 10 pixels. All maps are normalized to their mean value. The maps correspond to the case of a single 2D Gaussian (bottom left), a number of 2D Gaussians (top left), an inverted 2D Gaussian (bottom mid), a number of inverted 2D Gaussians (top middle), and a mix of 2D Gaussians and inverted 2D Gaussians (top right). The corresponding ∆-variance functions calculated for each case are displayed in the bottom right subpanel, and they are compared to the ∆-variance function of the underlying fBm image as well to the case of the same fBm smoothed with a Gaussian beam whose FWHM is 6 pixels.
For a faster and more efficient computation of Eq. 1, Ossenkopf et al. (2008) performed the calculation as a multiplication in Fourier space, and thus, the ∆-variance is given by where P is the power spectrum of A, and⊙ L is the Fourier transform of the filter function. If β is the exponent of the power spectrum, then a relation exists between the slope of the ∆variance and β (Stützki et al. 1998). This is given by The slope of the ∆-variance can be inferred from the range of spatial scales over which it displays a self-similar behavior. It can be tied to the value of β. Characteristic scales are scales at which there are breaks of the self-similarity and that show up in the ∆-variance plots as break points or inflection points. The error bars of the ∆-variance are computed from the counting error determined by the finite number of statistically independent measurements in a filtered map and the variance of the variances, that is, the fourth moment of the filtered map. The ∆-variance has been employed to analyze the structure of observed molecular clouds (e.g., Bensch et al. 2001;Sun et al. 2006;Rowles & Froebrich 2011;Elia et al. 2014;Dib et al. 2020) as well as simulated molecular clouds (e.g., Ossenkopf et al. 2001;Bertram et al. 2015). Elmegreen et al. (2001) in the only work that used the ∆-variance spectrum to characterize the structure of the H i gas. However, their study was limited to the Large Magellanic Cloud (LMC). Coherent (i.e., nonhierarchical) structures in a larger selfsimilar medium generate a bump in the ∆-variance spectrum. This is true regardless of whether the structure is an overdensity (i.e., a clump in a column density map) or a region of low column density (i.e., a hole or void). The reason is that the ∆variance measures the variance of an image over a given scale. This issue was discussed in detail in Dib et al. (2020). Here, we show a limited number of examples for completeness. Figure  1 displays five realizations of a fractal Brownian motion (fBm) image with an exponent of β = 2.4, on which we superimpose a single clump (bottom left), a number of identical clumps (top left), a single void (bottom middle), a number of identical voids (top middle), and a mixture of clumps and voids (top right). The size of each image is 1000 × 1000 pixels, and the clumps and voids are represented by 2D generalized Gaussian functions and inverted 2D Gaussians, respectively. In the examples displayed in Fig. 1, the standard deviation of the Gaussian functions in each direction is 10 pixels. The ∆-variance spectra for all cases are displayed in the bottom right subpanel of Fig. 1. Additionally, we show the ∆-variance spectrum of the pure fBm image. In the latter case, the ∆-variance spectrum is a power law with an exponent α = β − 2 = 0.4. As demonstrated in Dib et al. (2020), coherent structures in a self-similar medium increase the ∆-variance on all spatial scales. However, the point of maximum Article number, page 3 of 21 A&A proofs: manuscript no. dv_h1_v3 deviation from the spectrum of the underlying fBm (i.e, the peak of the bump) occurs on scales that are ≈ 4 √ σ 1 σ 2 as this is where most of the signal of the 2D Gaussian lies and where σ 1 and σ 2 are the standard deviations in both directions. It is important to mention here that the scale at which the maximum deviation (∆(σ 2 ∆ ) max ) occurs between the ∆-variance spectrum in the pres-ence of added structure and the spectrum of the underlying fBm does not necessarily correspond to the position of the peak, and it is generally smaller (see the schematic representation in Fig. 4).
The results discussed so far relate to the case of a pure fBm image and to cases in which discrete coherent structures are overlaid on the fBm image. We have also calculated a case of Fig. 3. Delta-variance spectra for the galaxies classified as dwarfs in the THINGS sample. The vertical dashed black line in each subpanel indicates the spatial resolution for each galaxy, and the vertical dash-dotted orange line corresponds to the optical radius of the galaxy. The values of R 25 are adopted from Walter et al. (2008). The spectra are normalized by their respective mean values. Fig. 4. Schematic figure representing the shape of the ∆-variance spectrum for the THINGS galaxies. In some galaxies, some of the features of the spectrum such as the bump at small scales or the presence of two distinct power laws are not observed. The dashed line represents the extrapolation of the first power law down to smaller scales. The quantity ∆ σ 2 ∆ represents the maximum deviation between the bump and the extrapolated power law. The value of L s f represents the physical scale at which this maximum deviation occurs. As illustrated, this scales does not correspond to the peak of the bump and is generally smaller. The quantity L tr is the physical scales at which a transition is observed between the first and second power law. a pure fBm (with β = 2.4) in which the image was smeared with a Gaussian beam whose full width at half maximum (FWHM) is D beam = 6 pixels. The ∆-variance spectrum corresponding to this case is also displayed in Fig. 1 (lower right subpanel). The effect of the reduced resolution is to cause a depression in the ∆-variance spectrum on scales 2D beam . The effects of beam smearing can extend to scales larger than D beam . However, on scales 1.5D beam , the effects of the beam do not exceed a few tens down to a few percent on larger scales. This effect has been presented and discussed in earlier studies using the ∆-variance (e.g., Bensch et al. 2001;Dib et al. 2020). The important aspect of this is that the effects of beam smearing do not generate any spurious bump or other features in the ∆-variance spectrum.

Results
We calculated the ∆-variance spectrum for all galaxies in the THINGS sample that are available in the online database (33 galaxies) after deprojecting them. The results are displayed in Fig. 2 and Fig. 3 for spiral and dwarf galaxies, respectively. The ∆-variance spectra for the spiral and dwarf galaxies exhibit a variety of features. However, some features are the same in many galaxies. In some galaxies, a bump in the spectrum is observed on scales of a few to several hundred parsec (≈ 100 − 850 pc). With the exception of NGC 2077, this feature is more commonly observed in spiral galaxies than in dwarfs. A second feature that can be observed in the ∆-variance spectra of the majority of THINGS galaxies in Fig. 2 and Fig. 3 is the existence of two selfsimilar regimes where the ∆-variance can be described by two power laws, σ 2 ∆ ∝ L α 1 on intermediate spatial scales (i.e., one to several kpc) and σ 2 ∆ ∝ L α 2 on larger scales (i.e., a few to several kpc). Figure 4 shows a schematic sketch of the common features that are observed in the ∆-variance spectra of the THINGS galaxies. We determined the boundaries of each self-similar regime by visual inspection and avoided any overlap with other features of the spectra (i.e., bumps, dips, and inflection points). We also used a different approach in which we computed the first-order derivative of the spectrum in order to evaluate in which range it is constant. We find that there is no particular advantage in following this approach because the value of the slope is never exactly a constant, and the range at which the spectrum starts to deviate from a power law is also not unambiguously determined. We find that the visual inspection with a careful selection of the ranges on which the spectrum is assumed to be a power law is as accurate as an automated selection.
We fit the self-similar regimes and determined for each galaxy the values of α 1 and α 2 . In a few galaxies (i.e., NGC 2976, NGC 3184, NGC 3351, and NGC 7331), only one selfsimilar regime can be identified, and given its steepness and the involved spatial ranges, we categorized it as being described by the second power-law function, whose exponent is α 2 . Following the logic described in §. 3, we determined the position at which the maximum deviation occurs, L s f , using the following procedure: We extrapolated the power law that comes after the bump down to spatial scales where the bump is observed, and we measured the difference between the observed spectrum and the extrapolation of the first power law on scales where the bump is located (see Fig. 4). The upper limit for the scales that were considered in this subtraction is the lower limit of where the first power law is assumed to be valid. The difference in the spectra was then fit with a Gaussian function, and the position of the peak of the Gaussian function was assumed to be the position where the maximum deviation occurs. In all cases, we found that the value of L s f is lower than the position of the peak of the bump. A note of caution is probably due. By performing the extrapolation of the first power law down to physical scales where the bump is observed and by subtracting it from the observed spectra on those scales, we assumed that the first power is the underlying slope of the spectra if there are no discrete coherent structures in the H i on these scales. This is not entirely guaranteed, and there is a possibility that the process that generated the bump in the spectrum (both position and width) can also affect the value of the first power-law slope in the spectra. On the other hand, it is important to stress that the values of L s f measured in this way are very close to the position of the bump in each spectrum. This implies that any misassumption on what the true slope of the spectrum might be at small spatial scales has a very minor effect on our results and conclusions. We also point out that the wing on the left-hand side of the bump is likely to be affected by the resolution of the observations, and as a consequence, both the amplitude and width of the bump are reduced. However, this  Notes. Columns represent the (1) name of the galaxy, (2) adopted inclination, (3) value of the exponent of the first power low α 1 , (4) 1σ uncertainty on α 1 , (5) value of the exponent of the second power law α 2 , (6) 1σ uncertainty on α 2 , (7) spatial range over which the first power law was fitted, (8) spatial range over which the second power law was fitted, (9) position of the characteristic scales, L s f , (10) 1σ uncertainty on L s f , and the (11) position of the transition between the first and second power laws L tr . The upper and lower groups of galaxies are the spirals and dwarfs, respectively.
should affect the position of the bump only marginally (see the application to the simulated galaxy in §. 5.1). Finally, the position of the transition point between the first and second power law (when present), L tr , is estimated by eye, and given the high uncertainties that affect the spectra at these large scales, there is little ground to expect that any automated procedure will yield more accurate estimates. In some galaxies, this transition appears as an inflection point (e.g., NGC 3521 and NGC 3621), whereas in others, a dip can be observed (e.g., NGC 7793 and Ho II). We defined the value of L tr as being either the position of the inflection point or the deepest position in the dip, when present. A conservative estimate of the uncertainty on L tr is about 10%. The values of L s f , α 1 , α 2 , and L tr for all of the THINGS galaxies are reported in Tab. 1, along with the uncertainties measured for L s f , α 1 , and α 2 and the spatial ranges over which every powerlaw fit was performed.
Most H i holes are unlikely to be circular as they are affected by local inhomogeneities in the local velocity and density field and by the effects of galactic shear (Dib et al. 2006, Bagetakos et al. 2011, Ohlin et al. 2019, Aouad et al. 2020, and thus, L s f is a measure of the effective radius. Figure 5 (bottom left panel) displays the distribution of L s f (dN/dL s f ) for the spiral and dwarf galaxies and for the combined sample. In dwarf galaxies, L s f 250 pc, while in spiral galaxies, L s f can be as large as ≈ 850 pc. The distributions of α 1 and α 2 ((dN/dα 1 ) and (dN/dα 2 ), respectively) are displayed in Fig. 5 (top subpanels). The mean values of α 1 are 0.79 ± 0.23, 0.43 ± 0.20, and 0.55 ± 0.27 for the dwarfs, spirals, and for the entire sample, respectively. The ∆-variance spectra are steeper on large spatial scales, and the mean values of α 2 are 1.56 ± 0.38, 1.42 ± 0.04, and 1.47 ± 0.46 for the dwarfs, spirals, and the entire sample, respectively. The distributions of the transition scale between the first and second power laws in the spectra, (dN/dL tr ), are displayed in the bottom right subpanel of Fig. 5. The distribution of L tr in dwarf galaxies peaks at ≈ 1 kpc and extends to ≈ 3 kpc. For spirals, the distribution of L tr is broader, with values that fall in the range of 2 to 12 kpc (Tab. 1).
Article number, page 7 of 21 A&A proofs: manuscript no. dv_h1_v3 Fig. 6. Angular size of L s f plotted against the projection-corrected angular beam size λ beam (top subpanel) and the ratio of these two quantities plotted against the distance of the galaxy (bottom subpanel).
As discussed above and in §. 3, the shape of the ∆-variance spectrum might be affected by resolution effects on sizes of about the beam size and smaller. In most galaxies, the peak of the bump and its right wing are well resolved, whereas the left wing is more affected by the resolution of the observations. We have shown that resolution issues cannot cause the occurrence of a bump similar to the one that is observed in many galaxies in Fig. 2 and Fig. 3. Here, we test further the effects of resolution by comparing the angular size of L s f (θ S F ) with the inclinationcorrected angular beam size, λ beam . Similarly to what was presented in Li et al. (2021), we calculated these two quantities as being and where 6 ′′ is the angular resolution of the beam, D gal is the distance of the galaxy, and i is the inclination angle. The values of D gal were taken from Walter et al. (2008) and the adopted values of i are those described in §. 2. The top subpanel in Fig. 6 indicates that there is no obvious correlation between θ s f and λ beam , and the bottom subpanel shows that the ratio θ s f /λ beam is not dependent on the distance to the galaxy. The Pearson correlation coefficient for the points in the (θ s f − λ beam ) scatter plot is P ≈ −0.0015. This clearly indicates that there is no linear correlation between these two quantities and that the determination of L s f is largely unaffected by the beam size. Furthermore, the ratio θ s f /λ beam is not constant and varies by a factor of ≈ 4 at given distance. This rules out that L s f could be an artifact of the data reduction that would be affecting the THINGS data on small spatial scales.

Interpretation
The bumps in the ∆-variance spectra that are observed in most spiral galaxies and in some of the dwarf galaxies at small spatial scales (i.e., a few to several hundred parsec; see Fig. 2 and Fig. 3) might be due to large H i complexes and to H i holes that are created either by feedback from massive stars or by other mechanisms, such as large-scale thermal and gravitational instabilities (Kim et al. 1999;Dib & Burkert 2004Silich et al. 2006;Weisz et al. 2009;Bagetakos et al. 2011, Cannon et al. 2012. While H i complexes may make a certain contribution to the bump, it is unlikely that they are the main source of the signal that is observed at these scales. This is simply because H i complexes are themselves self-similar in nature and are a consequence of the large-scale turbulence cascade. On the other hand, supernova remnants and bubbles are filled mostly with hot rarefied gas and are devoid of any significant H i gas emission (e.g., Walter et al. 2008, Bagetakos et al. 2011. Thus, H i holes are more similar to the coherent structures described in §. 3, and their sizes (or distribution of sizes) can have a direct imprint on the shape of the ∆-variance spectrum. The bump can be described by three quantities: its amplitude, the position of the point of maximum deviation from the underlying power law (L s f ), and its width. The scale at which the bump joins the first power law is related to the H i hole separation, as shown in Dib et al. (2020) and in Fig. 1. Both the amplitude of the bump and the value of L s f could be related to the SFR, which given a certain IMF, sets the frequency of type II supernova explosions in the disks. However, given that the first power-law slope might itself be perturbed by the existence of the bump, the value of L s f is likely to have a tighter correlation with the SFR than the bump amplitude.
In order to validate our results and gain more insight into the features we observe in the ∆-variance spectra of the THINGS galaxies, we measured the ∆-variance spectrum of a simulated galaxy. We tested the validity of the star formation activity related scenario by exploring the connection between the features observed in the ∆-variance spectra of the THINGS galaxies and indicators of their star formation activity. We also explore the possible origin of the broken power law that is observed in the ∆-variance spectra and interpret the transition point between the two self-similar regimes.

Insight from numerical simulations of whole galaxies
We used the VINTERGATAN cosmological zoom-in simulation of a Milky Way-like galaxy (Agertz et al. 2021, Renaud et al. 2021a. The simulation reaches a resolution of 20 pc in the densest medium and includes prescriptions for ultraviolet back-  ground radiation, atomic and molecular cooling lines in the form of tabulated data of the Sutherland & Dopita (1993) and Rosen and Bregman (1995) cooling curves, and a prescription for star formation. Stellar feedback from massive stars is accounted for in the form of stellar winds, radiation pressure, and type II and type Ia supernovae (see Agertz et al. 2021 for details). The global properties of the simulated galaxy agree with measurements of the Galactic mass, the surface density profiles of its baryonic components, the rotation curve, and the chemical bimodality of the stellar populations.
The analysis shown here was conducted at a look-back time of 3.5 Gyr, corresponding to a redshift of z ≈ 0.3 (we refer to the time of the canonical snapshot with time t = t 0 ), when the stellar mass of the galaxy was ≈ 6×10 10 M ⊙ and the mass of the atomic gas component was ≈ 5 × 10 9 M ⊙ . At this epoch, the effects of the last major merger had faded, and the galaxy was in its phase of secular evolution, with an SFR of ≈ 9 M ⊙ yr −1 . In order to measure the mass of H i in each simulation cell, we first solved the local Saha equation using the cell temperature and density. This allowed us to obtain the mass of the ionized gas. By subtracting it from the total gas mass, we measured the total mass of neutral gas. Using the prescription of Krumholz et al. (2009) 3 , we then computed the mass of molecular hydrogen at the local gas metallicity and removed it from the neutral gas mass, which gives the mass of atomic hydrogen. The surface density maps of H i were then computed in face-on projections of the simulation cells, were remapped at a uniform resolution of 50 pc, and covered a surface area of 50 kpc x 50 kpc. For each snapshot, we also generated a number of H i column density maps with various inclinations. The effect of inclination on the ∆-variance spectrum is discussed in Appendix A. Figure 7 displays the column density maps of the face-on simulated galaxy at t = t 0 (middle panel), at t = t 0 − 160 Myr (left panel) and at t = t 0 + 160 Myr (right panel). The spatial resolution in all three maps is 50 pc. The corresponding ∆-variance spectra for these three cases are displayed in Fig. 8. Figure 8 provides clear evidence that the galaxy is in its phase of secular evolution because the ∆-variance spectra of the three snapshots are nearly identical. The spectra display a prominent bump at ≈ 0.8−1 kpc and are similar to those that are observed in some of the THINGS galaxies, such as NGC 2976, NGC 3351, and NGC 7331, whose spectra also display a prominent bump and an absence of a first power law on intermediate spatial scales. In Fig. 8 we compare the ∆-variance spectra of the models to that of the galaxy NGC 7731 (L s f = 305 pc, α 2 = 1.42 ± 0.024). The spectra are very similar, but we had to multiply the spectrum of NGC 7331 by a factor of 3. This factor is simply due to the difference in the mean surface density of the H i maps between the simulation galaxy and NGC 7331. A fit to the self-similar regime of the models in the scale range [10−20] kpc yields a value of the slope α 2 = 1.34 ± 0.08, which is very similar to the value derived for NGC 7331. Applying the same procedure as in the observations, we derive a value of L s f ≈ 890 pc from the ∆-variance spectra of the simulation. Because most THINGS galaxies have a spatial resolution higher than 50 pc, we generated H i surface den-3 The implementation of the Krumholz et al. (2009) method to compute f H 2 in the version of the RAMSES code used in this work is described in detail in Agertz & Kravtsov (2015) (Eqs. 2, 3, and 6). We used σ d,−21 /R −16.5 = 1, where σ d,−21 is the dust cross-section per hydrogen nucleus to radiation at 1000 Å normalized to 10 −21 cm −2 , and R −16.5 is the rate for H 2 formation on dust grains, normalized to the Milky Way value of 10 −16.5 cm 3 s −1 (Wolfire et al. 2008). Both quantities are directly proportional to the dust abundance and thus to the gas-phase metallicity, which is tracked in each cell of the simulation. As in Agertz & Kravtsov (2015), we adopted a value of 3 for the parameter φ CNM and calculate the optical depth τ c that appears in Eq. 3 of Agertz & Kravtsov (2015) as τ c = ρ cell s cell σ d , where s cell and ρ cell are the size of the cell and its density, respectively. A&A proofs: manuscript no. dv_h1_v3  sity maps in which the 50 pc resolution map is convolved with a beam whose FWHM is 150 pc and 300 pc. The three maps with different resolution for the fiducial timestep (i.e., t = t 0 ) are displayed in Fig. 9. The corresponding ∆-variance spectra for these cases are displayed in Fig. 10. The loss of resolution does not affect the value of the power law that characterizes the shape of the spectrum at large spatial scales. However, the reduced resolution affects the shape of the bump and its position. Increased smoothing reduces the variance in the map and thus the amplitude of the bump decreases, and as discussed earlier, the left-hand wing of the bump also becomes increasingly affected and the bump becomes narrower. The position of the bump, and consequently, the value of L s f are also shifted to higher values (L s f ≈ 965 pc at the resolution of 300 pc). However, this effect, which is present in the ∆-variance spectra of the observed galaxies shown in Fig. 2 and Fig. 3, is not dramatic. The increase in the position of the bump and in the value of L s f is only ≈ 10%.

Relation between the ∆-variance spectra and galactic star formation
If the bump that is observed in the ∆-variance spectra is connected with stellar feedback, in particular, with supernova ex-plosions that can carve holes in the H i gas, then we might expect a correlation between the value of L s f and the global galactic SFR. This correlation is expected on the basis of the known empirical relation between the global SFR of a galaxy and the maximum mass of clusters that form within it (Weidner et al. 2004;Gónzalez-Lópezlira et al. 2012Schulz et al. 2015 ). More massive clusters will statistically harbor more massive stars, and the bubble that can form when massive stars explode as supernova would therefore be larger, leading to a correlation between the characteristic size of H i holes in galaxies and the global SFR. Figure 11 (top left panel) displays the value of L s f plotted as a function of the global galactic SFR. The values of the SFRs for the ensemble of the THINGS galaxies are taken from Leroy et al. (2008) and Walter et al. (2008). Figure 11 shows that for low values of the SFR (i.e., SFR 0.5 M⊙ yr −1 ), L s f is quasi-constant or weakly dependent on the SFR. This is expected if H i holes in these galaxies are caused by one (or a few) supernovae remnants. For higher SFRs (SFR 0.5 M ⊙ yr −1 ), we observe a correlation between the SFR and L s f . We also include the point of the simulated galaxy smoothed at a resolution that is close to the observational resolution of the observations (i.e., 300 pc). An empirically motivated second-order polynomial fit to the L s f -SFR data yields the following result: log(L s f ) = 2.48 + 0.23 log(SFR) + 0.05 (log(SFR)) 2 , and the fit is overplotted on the data in Fig. 11   The same procedure was adopted by Leroy et al. (2008) to derive the H 2 masses. The H 2 surface densities are those obtained by the HERACLES survey (Leroy et al. 2009), the BIMA SONG survey (Helfer et al. 2003) by Walter et al. (2001) for NGC 3077, and by Bolatto et al. (2008) for NGC 4449. Leroy et al. (2008) adopted a constant CO-to-H 2 conversion factor of X CO = 2×10 20 cm −2 (K km −1 s −1 ) −1 , which is about the mean value found in the Milky Way. The 1σ uncertainty on this quantity is ≈ 30% (Bolatto et al. 2013). Because the galactic SFR and stellar mass M * are correlated (e.g., Lara-López et al. 2013), L s f has the same dependence on M * as on the SFR (see details in Appendix B). Hence, L s f is independent of the sSFR, but the scatter between the two quantities is large due to the combined uncertainties on the SFR and M * . The correlation between L s f and SFE g is unclear as well because this last quantity is only loosely related to the star-forming gas. However, a weak anticorrelation is observed between L s f and SFE m . A power-law fit to the L s f -SFE m data points yields the following result: The anticorrelation between L s f and S FE m is largely due to the low-metallicity galaxies, which have a higher SFE m . Whether low-metallicity galaxies, such as the dwarf galaxies in the THINGS sample, have a higher SFE m is still a matter of debate. The difficulty in measuring SFE m in subsolar metallicity galaxies is the determination of the appropriate ratio between the column density of H 2 and the CO intensity (N(H 2 )/I CO ) (X CO factor). The X CO factor is expected to increase as the metallicity decreases due notably to increased photodissociation of CO (Bolatto et al. 2013 and references therein). Thus, a CO flux at low metallicity corresponds to a higher H 2 mass as compared to the same flux from a higher metallicity environment. A careful comparison of the CO intensity with the column density of H 2 as estimated via dust emission in the nearby galaxies M33 (Gardan et al. 2007;Braine et al. 2010;Gratier et al. 2017) and NGC 6822  led to a higher X CO than in the Milky Way, but nonetheless to a higher SFE than in large spirals (e.g., Murgia et al 2002). An obvious reason for this is that the conversion of H i into H 2 occurs at higher density when fewer dust grains are present (Hollenbach et al 1971;Braine et al 2001), as is the case in low-metallicity environments. This reduces the freefall time in the molecular component 4 . A second, more subtle, mechanism is that the weaker stellar winds in low-metallicity environments expel gas less efficiently from protocluster-forming clouds, such that a higher stellar mass can be formed for a given molecular gas mass, which leads to a higher SFE m Dib 2011, Dib et al 2013. These effects are cumulative. 4 The conversion rate is proportional to Z n 2 H i where Z is the metallicity and n H i is the number density of neutral hydrogen. This means that as Z decreases, the rate of conversion from H i into H 2 occurs at a higher density. This is valid under the assumption that the dust number density n dust ∝ n H , which is correct down to Z ≈ 0.1Z ⊙ ). The free-fall time is t f f ∝ n −0.5 H , which implies that t f f ∝ Z 0.25 .
Article number, page 11 of 21 Fig. 12. Relation between the galactic SFR and the slope of the first self-similar regime in the ∆-variance spectrum, α 1 (top subpanel), and the slope of the second self-similar regime (bottom subpanel).

Two self-similar regimes and the transition point
With the exception of a few galaxies (NGC 2976, NGC 3184, NGC 3351, and NGC 7331), the ∆-variance spectrum of most galaxies in the THINGS sample displays two distinct power-law regimes. On scales 0.5R 25 , the ∆-variance can be described by a power law with an exponent that varies between 0.1 and 1.16 and whose mean value is α 1 ≈ 0.5 . On larger scales ( 0.5R 25 ), it is described by a second power law, whose exponent can be as large as 2.3 and has a mean value of ≈ 1.5 (Tab. 1). A similar result was obtained for the LMC by Elmegreen et al. (2001). In terms of the exponent of the power spectrum, this would correspond to exponents of β 1 ≈ 2.5 and β 2 ≈ 3.5. A detailed explanation of the specific values of α 1 and α 2 , and consequently of L tr , for each individual galaxy in the THINGS sample is beyond the scope of this work. This would require comparisons with numerical simulations of a cosmological volume that resembles the ensemble of nearby galaxies. Intuitively, the shape of the ∆-variance spectra that are observed for the THINGS samples of galaxies might be thought to be the result of an exponential disk. While we show in Appendix C that an exponential disk could indeed generate a spectrum with a broken power law, we discard this hypothesis on the basis that no exponential disks are observed in H i. This fact was noted earlier by other authors (e.g., Casasola et al. 2017). Instead, H i disks are observed to be nearly flat (i.e., nearly constant column density), with radial variations by a factor of ≈ 2 and in some cases a depression toward the inner regions of the galaxy, where most of the hydrogen gas becomes molecular Leroy et al. 2008).
The mean values of α 1 and α 2 are so different that the structure of the ISM in the range of spatial scales they represent must originate from different physical processes or correspond to different phases of the gas with different compression levels. The values of α 1 we find in this work, with a mean value of ≈ 0.5, are very similar to those found for molecular clouds using either molecular transitions or cold dust emission (e.g., Stutzki et al. 1998;Bensch et al. 2001;Miville-Deschênes et al. 2010;Dib et al. 2020). They are also similar to those found in H i seen in absorption which is tracing the cold (≈ 100 K) component of the H i gas (Deshpande et al. 2000). These values are also consistent with those found in numerical simulations of supersonic magnetohydrodynamic turbulence where the gas is compressed into smaller pockets (e.g., Kowal et al. 2007;Dib et al. 2008) 5 . In contrast, the range of values of α 2 is consistent with the values observed for the H i in emission toward diffuse regions that are dominated by the warm component of the H i gas both in the Galaxy (e.g., Miville-Deschênes et al. 2003;Chepurnov et al. 2010) and in the LMC (Elmegreen et al. 2001). The mean value of α 2 ≈ 1.5, which corresponds to an exponent of the power spectrum of β ≈ 3.5, is consistent with the picture in which turbulence in the warm neutral medium (WNM) is subsonic to transonic (e.g., Burkhart et al. 2013).
As discussed above, a bump on scales of a few to several hundreds of parsecs perturbs the underlying self-similar regime on these scales. Although we have avoided any overlap with the bump when we performed the fit for the first power law, we explored whether any correlation exists between the SFR and the exponent of the first power law, α 1 . Figure 12 (top panel) shows a weak anticorrelation between the SFR and α 1 . Weakly starforming dwarf galaxies have a systematically steeper spectrum on spatial scales that are covered by the first self-similar regime. A power-law fit to the SFR-α 1 data points yields The SFR-α 1 anticorrelation suggests that the star formation activity in galaxies shapes the structure of the gas distribution on scales larger than those associated with the sizes of individual supernova remnants or larger superbubbles, up to the transition scale L tr . Lower values of α 1 imply that more substructure is present in galaxies with a higher star formation rate. In contrast, the exponent of the power law that describe the second selfsimilar regime (α 2 ) is independent of the SFR (Fig. 12, middle panel), indicating that the dynamics of the gas on large scales, and consequently its structure, are shaped by processes that act on scales larger than those associated with feedback from supernova explosions.
The transition between the two regimes is observed in Fig. 2  and Fig. 3 as a dip or an inflection point. As stated in §. 4, we adopted as the value of L tr the position of the inflection point or the deepest position of the dip, when present. The values of the derived L tr are reported in Tab. 1, and the distributions of L tr for dwarf and spiral galaxies are displayed in Fig. 5 (bottom  right panel). The top panel of Fig. 13 displays the dependence of L tr on the galactic SFR. The correlation between the L tr and the galactic SFR (and with M * , see Appendix B) is clear. Figure  13 (bottom panel) displays the distribution of the ratio of L tr to the galactic optical radius, R 25 . While there is some scatter, most values of this ratio lie around (L tr /R 25 ) ≈ 0.4 − 0.5. Interestingly, the value of ≈ (0.4 − 0.5)R 25 is very similar to the size of the molecular disk in the THINGS galaxies . Our results do not support the idea that L tr is connected to the scale height of the H i disk, h H i . Recent derivations of h H i for a number of the THINGS galaxies clearly indicate that the H i disks are flared and the values of h H i vary radially from about 100 pc in the inner region of the disk to ≈ 1 kpc in its outer region (Bacchini et al. 2019;Patra 2020a,b). These values are lower than any of the values of L tr derived in this work.
Because the exponential disk scenario can be discarded, a different physical mechanism must cause the broken power law and the transition point that are observed in the H i ∆-variance spectra of most THINGS galaxies. The evidence gathered from the distributions of α 1 , α 2 , and L tr and the anticorrelation between α 1 and the galactic SFR appears to point out to the following scenario: The ∆-variance spectrum in the first self-similar regime is dominated by emission from the CNM component of the H i gas. The range of values found for α 1 around a mean value of ≈ 0.5 is compatible with the occurrence of compressible supersonic turbulence that governs the dynamics of cold gas, that is, the cold component of the H i gas (e.g., Bensch et al. 1998;Bertram et al. 2015;Dib et al. 2020). The anticorrelation between α 1 and the SFR is compatible with the idea that the dynamics of the gas is increasingly dominated by compressive motions for an increasing SFR, leading to a shallower spectrum (i.e., a lower value of α 1 ). The reason is that compressive turbulence for the same Mach number can compress gas to higher overdensities than solenoidal modes (e.g., Federrath et al. 2008). On large spatial scales (i.e., 0.5R 25 ), the signal is dominated by the contribution from the external regions of the galaxy where the WNM phase of the H i is larger. At large galactocentric radius, where the H i is more flared, the gas is more easily affected by ram pressure stripping and the heating by the extragalactic background UV field. The combination of these processes keeps the gas warm and diffuse, and thus its dynamics is governed by subsonic to mildly supersonic turbulence with little connection to the galactic SFR.
The signature of the H i gas in the CNM phase extends only to scales L tr ≈ (0.4 − 0.5) R 25 because most of the cold H i resides in the inner region of the galaxy, on scales smaller than and up to the size of the molecular disk, and there is little or no cold H i in the outer regions of the galaxy (e.g., Braun 1997;Zhang et al. 2012). In contrast, the warm component of the H i gas is likely to be present everywhere in the disk, but in smaller proportion (in terms of total local mass) in the inner regions and dominant in the outer regions. The warm H i dominates the emission on larger scales and has less substructure on smaller scales. It can still contribute to the signal on small scales, but most of the variance on those scales is dominated by the cold component. In order to illustrate this idea, we show in Fig. 14 (left subpanel) a toy model in which the inner parts of the disk are described by a fBm whose β = 2.4, characteristic of cold gas, and the outer regions are described by as second fBm with a steeper spectrum (β = 3.4), which is characteristic of warm gas. Both fBms are normalized by their mean values. The size of the map is 1000 × 1000 pixels and the inner fBm is contained within a region whose radius is 300 pixels. The corresponding ∆-variance spectrum of this model displays a broken power law with a transition point located at a scale of ≈ 50 − 100 pixels. The break point is smaller than the imposed transition radius of 300 pixels because the low-intensity regions in both fBm have the same values, and thus the outer fBm has a non-negligible contribution to the signal on small scales because it covers a large surface area. Fig. 14 (middle and right subpanels) shows that depressing the value of the outer fBm by a certain factor (here 2 and 4, respectively) results in a spectrum in which the transition between the two regimes in the ∆-variance spectra is sharper and the transition point nearer to the radius of the inner fBm. Despite being overly simplistic in comparison to a real galactic disk in which the stable cold and warm gas phases can locally coexist with gas in the unstable regime (e.g., Dib et al. 2006), this toy model shows that a dominant cold H i component in the inner region of a galactic disk and a warm H i component that dominates the emission in the outer regions of the disk can explain the broken power law that is observed in the ∆-variance spectra of the H i 21 cm emission line in the THINGS galaxies.

Discussion and connection to previous work
Several other studies have explored the structure of the diffuse ISM on the scale of entire galaxies either using the H i 21 cm line emission or other tracers of the diffuse gas, such as dust mid-to A&A proofs: manuscript no. dv_h1_v3 The maps have a resolution of 1000 × 1000 pixels. The inner fBm, residing within a circle with a radius of 300 pixels, has β = 2.4, and the outer fBm has β = 3.4. The top left subpanel corresponds to the fiducial case, and the middle and right top subpanels correspond to cases in which the inner and outer fBms have been divided by a factor of 2 and 4, respectively. The lower subpanels display the corresponding ∆-variance spectra. The values of α = 0.4 and 1.4 are not fits to the spectra, but are shown as a reference to guide the eye.
far-infrared emission (e.g., Koch et al. 2020). These studies can be sorted into two main categories. On the one hand, there are studies that used isotropic methods in order to characterize the structure of the H i gas distribution in galactic disks, such as the calculation of the power spectra of the H i intensity map or of the line-of-sight velocity fluctuations (e.g., Begum et al. 2006;Szotkowski et al. 2019;Nandakumar & Dutta 2020), the auto-correlation function of the H i intensity map  or the ∆-variance spectrum (Elmegreen et al. 2001). Elmegreen et al. (2001) computed the power spectrum and the ∆-variance spectrum of the H i intensity for the LMC. They found two distinct power-law regimes with a transition that occurs at ≈ 250 − 300 pc. They were unable to explore the dependence of the shape of the ∆-variance spectra on the galactic star formation activity as their study was restricted to the case of a single galaxy.  and  measured the power spectrum of the H i intensity for a number of the THINGS galaxies. While their approach and data sample overlaps with ours, they favored fitting the entire spectrum of each galaxy with a single power-law function. Their approach could be entirely valid over specific spatial ranges in each galaxy, however, there are instances where a single power-law fit cannot be justified (e.g., see the case of NGC 3184 in Fig. 2 of . Furthermore,  did not find any correlation between the exponent of the power spectra and some of the galactic properties they have considered, such as the inclination, the H i and dynamical masses of the galaxy, and the surface density of the SFR. At first glance, our results might seem to contradict those of Combes et al. (2012) for M33 and Szotkowski et al. (2019) in the Large and Small Magellanic Clouds. Those authors found a steep spectrum on "small" scales and a shallower spectrum at "larger" scales. However, because of the relative proximity of these galaxies, the small scales in these works refer to scales that are not resolved in the THINGS sample. The shallow slopes they find at "larger" scales are similar to those we measure in our study over the same range of spatial scales.
Another approach for studying the structure of the ISM relies on the identification of discrete structures in galactic disks and on quantifying their statistical properties. This approach has been employed to detect H i holes in the Milky Way (Ehlerová & Palouš 2005. Ehlerová & Palouš (2013) measured the size distribution of shells in the Leiden, Argentina, Bonn H i survey and found that it can be fit with a power-law function (dN/dR shell ) ∝ R −ξ≈−2.6 shell . Oey & Clarke (1997) noted that if the H i shells are the results of feedback from massive stars, then a relation exists between ξ and the exponent of the power-law function that describes the luminosity distribution of OB associations, φ (L) ∝ L −η , such that ξ = 2η − 1. Ehlerová & Palouš (2013) already noted that the value of ξ ≈ 2.6 they have derived implies that η ≈ 1.8, which is close to the value that is derived from observations (≈ 2, McKee & Williams 1997). Dib et al. (2009) found that the orientations of the main axis of molecular clouds in the outer Galaxy are correlated on spatial scales that are approximately the expected sizes of supernova remnants that are found in these regions of the Galaxy. The results of Dib et al. (2009) and Ehlerová & Palouš (2013) clearly suggest that feedback processes from massive stars play an important role in shaping the structure of the ISM in the Milky Way. In nearby galaxies, Bagetakos et al. (2011) searched for H i shells in a subsample of the THINGS galaxies using both a simple morphological selection and a selection based on localized expansion velocity. They found that the size distribution of the H i shells in each of the THINGS galaxies peaks at a few to several hundred parsecs (see Fig. 3 in their paper), which is broadly similar to the position and width of the bump observed in the ∆-variance spectrum displayed in Fig. 2 and Fig. 3. Using numerical simulations, Yadav et al. (2017) showed that multiple supernovae remnants can merge to form large bubbles with sizes that range between ≈ 100 pc and 700 pc, which is very similar to the range of L s f values we find in this work. Figure 15 (top subpanel) displays the median size of the H i shells (R sh,med ) from Bagetakos et al. (2011) plotted as a function of the galactic SFR. The similarity between the SFR-R sh,med and the SFR-L s f relations is striking. Like in the case of the L s f , R sh,med shows no significant dependence on the SFR for value of the SFR 0.5 M ⊙ yr −1 and a positive correlation at higher values. The lower subpanel of Fig. 15 displays R sh,med plotted as a function of L s f . A clear one-to-one correlation exists between these two quantities, and this confirms the supernovae feedback origin of the characteristic scale L s f that is measured from the ∆-variance spectra. With the exception of one outlier galaxy, the values of R sh,med seem to be systematically higher than those of L s f . A plausible explanation is that because the H i holes in Bagetakos et al. (2011) were entirely identified by eye, the selection favored the identification of the largest holes, and some of the smaller H i holes went undetected, especially those that could be elongated or deformed.

Conclusions
We analyzed the structure of the H i gas using the order zeromaps of 33 galaxies taken from the THINGS survey . In order to characterize the H i gas structure, we calculated the ∆-variance spectrum (Stützki et al. 1998;Ossenkopf et al. 2008). Most spectra possess common features that include a bump at scales of a few to several hundred parsec, a first selfsimilar regime at intermediate spatial scales, and a transition to a second, steeper, self-similar regime at larger spatial scales. When extrapolating the first power law to smaller scales and subtracting it from the observed spectra, we were able to measure the position of the maximum deviation between the spectra and the underlying power law, L s f . We find that L s f , whose values range from one to several hundred parsecs, correlates with the galactic SFR for SFR values 0.5 M ⊙ yr −1 . Below this value, L s f is independent of the SFR. We also find a strong correlation between the value of L s f for each galaxy and the median size of H i shells measured by Bagetakos et al. (2011). Both findings clearly suggest that L s f is a measure of the characteristic size of the H i shells in the THINGS galaxy. The first similar regime is observed to extend from beyond the bump up to a spatial scale of ≈ 0.5R 25 , and it can be described by a power law whose exponent ranges from ≈ 0.1 to 1 with a mean value of ≈ 0.55. This exponent is compatible with the occurrence of compressible supersonic turbulence, which governs the dynamics of the cold component of the H i gas. On larger spatial scales (i.e., 0.5R 25 ), the structure of the H i gas can be characterized by a second power law whose exponent is found to vary between ≈ 0.5 and 2.5 with a mean value of 1.47. We find that the values of α 2 do not correlate with the galactic SFR. This regime corresponds to the dynamics of the gas being governed by subsonic to transonic turbulence. It therefore is a signature of emission from the warm component of the H i gas 6 . The transition point between the two selfsimilar regimes, L tr , is found to correspond to a spatial scale of ≈ 0.4 − 0.5R 25 . Interestingly, in most THINGS galaxies, this scale is about the size of the molecular disk, and this is probably an indication of where most of the cold H i gas resides.
Earlier work on the scale of molecular clouds (scales of 5 to 50 pc) using the ∆-variance technique has allowed us to uncover characteristic scales of ≈ 1 pc in massive star-forming regions such as Cygnus X (Dib et al. 2020). These scales are thought to be associated with the sizes of hubs where stellar clusters form. With the current data from the THINGS survey, it is not possible to probe the connection between what we observe in the H i, and particularly, its cold component with the structure of molecular clouds. Future observations with the Square Kilometer Array will allow us to start probing scales that are about 4-20 pc for galaxies located at distances of ≈ [1 − 7] Mpc (Tolstoy et al. 2010). Combined with high-resolution observations of molecular clouds both in the Galaxy and in nearby galaxies using the Atacama Large Millimeter Array (ALMA), we will be able to A&A proofs: manuscript no. dv_h1_v3 probe the link between the structure observed in the cold H i gas and that seen in the submillimeter and in molecular line transitions. This will allow us to explore the effects of feedback in the pre-supernova phase on the structure of the ISM in greater detail.