Turbulence statistics of H i clouds entrained in the Milky Way’s nuclear wind

ABSTRACT

The interstellar medium (ISM) is ubiquitously turbulent across many physically distinct environments within the Galaxy. Turbulence is key in controlling the structure and dynamics of the ISM, regulating star formation, and transporting metals within the Galaxy. We present the first observational measurements of turbulence in neutral hydrogen entrained in the hot nuclear wind of the Milky Way. Using recent MeerKAT observations of two extra-planar H i clouds above (gal. lat.|$\, \sim 7.0^{\circ }$|⁠) and below (gal. lat.|$\, \sim -3.9^{\circ }$|⁠) the Galactic disc, we analyse centroid velocity and column density maps to estimate the velocity dispersion (σ_v,3D), the turbulent sonic Mach number (⁠|$\mathcal {M}$|⁠), the volume density dispersion (⁠|$\sigma _{\rho /\rho _0}$|⁠), and the turbulence driving parameter (b). We also present a new prescription for estimating the spatial temperature variations of H i in the presence of related molecular gas. We measure these turbulence quantities on the global scale of each cloud, but also spatially map their variation across the plane-of-sky extent of each cloud by using a roving kernel method. We find that the two clouds share very similar characteristics of their internal turbulence, despite their varying latitudes. Both clouds are in the sub-to-trans-sonic Mach regime, and have primarily compressively driven (b ∼ 1) turbulence. Given that there is no known active star formation present in these clouds, this may be indicative of the way the cloud–wind interaction injects energy into the entrained atomic material on parsec scales.

turbulence, ISM: clouds, ISM: kinematics and dynamics, Galaxy: centre

1 INTRODUCTION

Turbulence is one of the main driving forces that shapes the evolution of the interstellar medium (ISM) in galaxies. To understand the mechanisms governing the life cycle of the ISM, we need to characterize the turbulence we observe in a diverse range of galactic environments through comparison with theoretical models and simulations (see review by Burkhart 2021). Statistical methods such as the spatial power spectrum (SPS) (e.g. Stanimirovic et al. 1999; Stanimirović & Lazarian 2001; Kowal & Lazarian 2007; Heyer et al. 2009; Chepurnov et al. 2015; Nestingen-Palm et al. 2017; Pingel et al. 2018; Szotkowski et al. 2019) and the probability distribution function (PDF; e.g. Burkhart et al. 2010; Patra, Chengalur & Begum 2013; Bertram et al. 2015; Maier et al. 2017; Nestingen-Palm et al. 2017) are important because they allow us to compare the distinctly two-dimensional (2D) information we can measure on the plane of the sky with the three-dimensional (3D) information we can access in increasingly sophisticated (magneto)hydrodynamic (MHD) simulations of the ISM (e.g. Kim & Ostriker 2017, 2018; Kim et al. 2023).

The density PDF of simulated isothermal, supersonic gas can be described by a log-normal function, meaning that the logarithm of the density follows a normal Gaussian distribution (Vazquez-Semadeni 1994; Padoan, Jones & Nordlund 1997; Passot & Vázquez-Semadeni 1998; Federrath, Klessen & Schmidt 2008; Hopkins 2013; Squire & Hopkins 2017; Beattie et al. 2021). The fluctuations of the density field are captured in the width of this PDF, a key parameter in starformation theories (Krumholz & McKee 2005; Hennebelle & Chabrier 2011; Padoan & Nordlund 2011; Federrath & Klessen 2012; Burkhart & Mocz 2019). Furthermore, studies by Price, Federrath & Brunt (2011), Konstandin et al. (2012), Molina et al. (2012a), Federrath & Banerjee (2015), Nolan, Federrath & Sutherland (2015), and Kainulainen & Federrath (2017) have shown that the width of the density PDF is proportional to the turbulent (sonic) Mach number. Formally, this is described by

$$\begin{eqnarray} \sigma _{\rho /\rho _0}= b\mathcal {M}, \end{eqnarray}$$

(1)

where |$\sigma _{\rho /\rho _0}$| is the standard deviation of the PDF of the 3D density field (ρ), scaled by the mean density (ρ₀), |$\mathcal {M}$| is the turbulent sonic Mach number, and b is a constant of proportionality known as the turbulence driving parameter (Federrath et al. 2008, 2010). This constant of proportionality is another key ingredient in star formation theories, and describes the ratio of compressive and solenoidal modes in the acceleration field of the gas, which is what drives the turbulence. We refer the reader to the large body of work testing the theoretically derived b parameter in numerical studies of driven turbulence (e.g. Federrath et al. 2008, 2010; Price et al. 2011; Molina et al. 2012a; Federrath & Banerjee 2015; Nolan et al. 2015; Beattie et al. 2021), as well as to a more in-depth discussion of the functionality of b in an observational context in Gerrard et al. (2023, hereafter G23). G23 described a new method for extracting |$\sigma _{\rho /\rho _0}$|⁠, |$\mathcal {M}$|⁠, and b (which we refer to as the main analysis quantities within this text) from observations, specifically from H i observations of the Small Magellanic Cloud (SMC).

Recent high-resolution observations of the large-scale wind activity in the centre of the Milky Way has revealed significant amounts of molecular and atomic gas entrained in the wind (Di Teodoro et al. 2018, 2020; Noon et al. 2023, hereafter N23). In this study, we will use the G23 techniques to probe the internal turbulence of two such extra-planar H i clouds. Extending our analysis to extra-planar clouds, positioned at significant angular distances from the Galactic plane, introduces an additional layer of complexity to our understanding of interstellar dynamics. These clouds present unique laboratories for studying turbulence in the ISM without the influence of ongoing internal star formation activity as is present in most cold clouds in the Galactic disc.

Theoretical studies encounter difficulties when attempting to simulate a fully entrained cold gas cloud in a hot wind before instabilities destroy the cloud (Schneider & Robertson 2017; Zhang et al. 2017; Girichidis et al. 2021), whilst other simulations find that it might be possible given specific criteria such as a large initial cloud mass or molecular gas forming via self-shielding in an entrained atomic cloud (Armillotta, Fraternali & Marinacci 2016; Armillotta et al. 2017; Gronke & Oh 2018; Farber & Gronke 2022; Chen & Oh 2023). Observational analysis of three such clouds has revealed that they were possibly molecular-dominated in the disc, but have since been accelerated in the nuclear wind, triggering the photo-dissociation of the molecular gas (N23). Di Teodoro et al. (2018) and Lockman, Di Teodoro & McClure-Griffiths (2020) use a kinematic, bi-conical wind model to determine the distances to each cloud from the Galactic centre (GC), revealing that the intermediate-latitude cloud also possesses broader CO line profiles, a disordered velocity field, and a wider full width at half-maximum (FWHM) range than the closest cloud. Further, the most distant cloud possesses no detectable molecular gas. The differences in cloud properties as a function of distance could indicate a cloud’s evolution in the Galactic wind, suggesting they slowly lose molecular content as they transit the wind. Alternatively, it could indicate varying localized conditions in the wind or stem from the relationship between cloud properties and survival in a wind. Probing the internal turbulence of these clouds could further reveal their history. Our aim is to see whether the cloud’s turbulence characteristics vary between the clouds or with distance from the GC.

The rest of this work is structured as follows: in Section 2, we outline our data preparation and the modifications we have made to the G23 methods, in Section 3, we discuss the results of our analysis and present some potential interpretations. We summarize our results in Section 4.

2 METHODS

In this section, we provide an overview of the methods used to extract and map the turbulent properties of the H i gas in two extra-planar clouds. Specifically, we recover the main analysis quantities |$\sigma _{\rho /\rho _0}$|⁠, |$\mathcal {M}$|⁠, and b, as well as the intermediate analysis quantities σ_{v, 3D}, |$\mathcal {R}^{1/2}$|⁠, and c_s (described in full in the following sections). Details of the methodology are outlined in G23. For this work, we have made some modifications to the method so that the pipeline can analyse data cubes that have a significant amount of empty pixels (see Section 2.2.1), as well as adding a low-pass filtering method (further discussed in Section 2.2.2).

2.1 Data

Here, we will outline the characteristics of the data relevant to the analysis presented in this paper. For a more detailed summary on the observation set-up and data reduction, see Di Teodoro et al. (2020) and N23. We analyse two of the three extra-planar clouds presented in N23, C1 (closest to the GC), and C3 (furthest from the GC). The intermediate-latitude cloud, C2, has a particularly complex velocity structure, as is likely a conglomerate of several distinct structures coincident in position–position–velocity (PPV) space. Our methods are not applicable to this kind of object, and as such, we have excluded it from the analysis.

We observed H i 1.42 GHz and ¹²CO(2→1) 230.538 GHz emission lines using the MeerKAT radio interferometer in late 2020 (SARAO project DDT-20201006-NM-01; principal investigator NMG) and the Atacama Pathfinder EXperiment (APEX) telescope (ESO project 0104.B-0106A; principal investigator EMDT), respectively. The H i observations have a spectral resolution of 5.5 km s⁻¹ and a spatial resolution of 24.9 arcsec × 21.7 arcsec with position angle −2.8°. The CO observations have a spectral resolution of 0.2 km s⁻¹ and a spatial resolution of 28 arcsec. The corresponding root-mean-square (rms) noise in the H i data cubes (σ_chn) are 304 and 365 mK for C1 and C3, respectively, per 5.5 km s⁻¹ channel. The rms noise of the CO data cubes is 65 mK per 0.25 km s⁻¹ channel for both clouds.

2.2 Data preparation

Before inputting the column density and centroid velocity maps of each cloud to the G23 pipeline, we first perform a few data-cleaning processes. The first step is to make a per-channel signal-to-noise ratio (SNR) cut at 3σ_chn. This removes much of the low-SNR emission from around the central bulk of the cloud, as well as spurious signals. We then create the moment-0 map (M₀), given by

$$\begin{eqnarray} \mathrm{M_0}(x,y) = \int T_b(x,y,v)\, \mathrm{ d}v, \end{eqnarray}$$

(2)

where T_b(x, y, v) is the brightness temperature (K) in the velocity channel v and in the pixel (x, y), and dv is the channel spacing (km s⁻¹). We calculate the column density in cm⁻² of each cloud by multiplying the M₀ by a factor of 1.823 × 10¹⁸ cm⁻² (K km s⁻¹)⁻¹, assuming the gas is optically thin (McClure-Griffiths, Stanimirović & Rybarczyk 2023).

The moment-1 (M₁) is the intensity-weighted centroid velocity, and is given by

$$\begin{eqnarray} \mathrm{M_1}(x,y) = \frac{\int v\, T_b(x,y,v)\, \mathrm{ d}v}{\int T_b(x,y,v)\, \mathrm{ d}v} =\frac{\int v\, T_b(x,y,v)\, \mathrm{ \mathrm{ d}}v}{\mathrm{M_0}}, \end{eqnarray}$$

(3)

where v is the velocity of each channel.

Post integration, we construct a map of the noise in the integrated intensity map to perform another SNR cut. We follow the method outlined in Lelli, Verheijen & Fraternali (2014, appendix A1, equation A5) for a uniform tapered cube that has been continuum-subtracted, defining

$$\begin{eqnarray} \sigma _{\mathrm{M0}} = \sqrt{1+0.25 N_{\mathrm{chn}} (1/N_1 + 1/N_2)}\times \sqrt{N_{\mathrm{chn}}}\times \sigma _{\mathrm{chn}}\mathrm{ d}v, \end{eqnarray}$$

(4)

where σ_M0 is the noise in the integrated map, σ_chn is the intrinsic noise per channel, N_chn is the number of channels integrated to make the moment-0 map, and N₁ and N₂ are the number of channels used during continuum subtraction from each end of the spectral axis of the cube (N₁ = 86, 133 and N₂ = 96, 49 for C1 and C3, respectively). Once we have constructed a map of the noise, we create an SNR map, and use it to threshold the moment-0 map. We choose to include everything above 5σ_M0 in the maps for the following analysis.

Lastly, we use a watershed algorithm (Fiorio & Gustedt 1996; Wu, Otoo & Shoshani 2005) to find the largest region connected pixel-by-pixel in the moment-0 map and mask out any remaining pixels outside of that region, so as to exclude any small ‘islands’ of emission. This mask is also applied to the moment-1 maps. The resultant maps are presented in Fig. 1.

Figure 1.

Moment maps for each of the clouds analysed in this work. The left-hand panels show the logarithmic normalized column density, while the right-hand panels show the centroid velocity, with C1 on the top and C3 on the bottom. The orange circle in the bottom left-hand corner of each figure shows the FWHM of the kernel, while the small white circle inside the kernel represents the beam of the H i observations. The black contours in the top panels show H₂ column density at 5.0 × 10²¹ cm⁻². The scale bar shows 10 pc, assuming the physical distance to the clouds is 8.3 kpc. In each panel, the arrow points toward the Galactic Centre (GC), indicating the opposite direction of the wind’s flow.

2.2.1 Roving kernel

As in G23, we use a roving kernel to map the turbulent quantities across the on-sky extent of the clouds. This is a Gaussian kernel, truncated at 3σ. The method in G23 only computes the analysis quantities in kernels which are entirely filled with valid data, and assigns NaN to any kernel which has empty pixels inside of it. In this new iteration of the method, kernel windows are permitted to contain empty pixels, but must have a minimum number of valid pixels contained within the kernel FWHM to be processed. In essence, this allows us to compute statistics in kernels that contain a few NaN pixels, such as at the edges of the clouds, while still maintaining confidence in those statistics. As described in Sharda et al. (2018), the minimum number of resolution elements required to recover a reasonably converged standard deviation is of order 20, so we set a minimum of 20 beams worth of viable pixels inside the kernel’s FWHM for that kernel to be processed. We chose the actual size of the kernel such that the FWHM is the diameter of a circle that has an area 3 times larger than the minimum number of valid pixels. This is discussed further in Appendix  A. This change to the G23 method allows us to analyse objects with a comparable number of independent resolution elements to the minimum size of the kernel FWHM, which was not necessary in the original work because the kernel was small compared to the entire SMC data cube and the number of beams contained within it.

2.2.2 Low-pass filtering

In addition to the kernelled maps that allow us to spatially resolve turbulent properties of the clouds, we also compute global statistics for each cloud, meaning that we calculate the main and intermediate analysis quantities using the entire cloud to measure the variance in density and velocity. This allows us to make a comparison of the properties on two different spatial scales: the minimum size that the variance can be resolved in, and the largest scale, i.e. the entire cloud. In order to isolate the purely turbulent fluctuations in both density and velocity on the whole-cloud scale, we perform a Gaussian convolution of the column density and centroid velocity maps, to create a low-pass filter (LPF) map that we then subtract from the original maps. By using this convolution method to construct the large-scale variations across each cloud, we can account for their unique geometry and clumpy H i emission structure. The size of the convolution kernel is dictated by the size of the cloud: we take the total area of the cloud and use the radius of a circle with that same area as the FWHM of the Gaussian. We use these LPF-corrected maps to compute the global analysis quantities of each cloud, as well as the input to the roving kernel method, although for continuity of the method we still perform a linear gradient subtraction on the kernel scale. Whether or not we pass the LPF-corrected maps or the raw moment-0 and moment-1 maps does not effect the results of the kernelling method, however. In Appendix  B the results of the LPF-method can be seen for both density (Figs. B1 and B2) and velocity (Figs. B3 and B4), for each cloud.

2.3 Turbulence analysis

Our approach to extracting the density and velocity statistics for each cloud centres around isolating the turbulent fluctuations in the column density and centroid velocity maps, and using these plane-of-sky quantities to reconstruct the associated 3D quantities required for estimating the volume density dispersion, the 3D velocity dispersion, the Mach number, and the turbulence driving parameter for each cloud as per equation (1). We compute these quantities in a roving Gaussian kernel, where we subtract a kernel-scale gradient from the LPF-corrected moment-0 and moment-1 maps, so as to isolate the turbulent fluctuations before measuring the plane-of-sky dispersion. An in-depth description of the pipeline can be found in G23, but the main points are summarized in the following and changes with respect to G23 highlighted where applicable.

2.3.1 Density statistics

To convert the plane-of-sky column density dispersion |$\sigma _{N_{\mathrm{H I}}/N_0}$| to the volume density dispersion, |$\sigma _{\rho /\rho _0}$|⁠, we follow the method outlined in Brunt, Federrath & Price (2010). The relation between the 2D density power spectrum P_2D(k) and the 3D density power spectrum P_3D(k) is given by

$$\begin{eqnarray} P_\mathrm{3D}(k) = 2k P_\mathrm{2D}(k), \end{eqnarray}$$

(5)

where k is the wave number (Federrath & Klessen 2013). We exploit this relation to recover P_3D(k) by first computing P_2D(k) of the gradient-subtracted column density, which immediately gives us P_3D(k) of the quantity ρ/ρ₀−1 as per equation (5). The ratio of the sums over k-space of these two quantities gives the density variance ratio (Brunt et al. 2010),

$$\begin{eqnarray} \mathcal {R}^{1/2} = \frac{\sigma _{N_{\mathrm{HI}}/N_0}}{\sigma _{\rho /\rho _0}} = \left(\frac{\sum _k P_\mathrm{2D}(k)}{\sum _k P_\mathrm{3D}(k)}\right)^{1/2}, \end{eqnarray}$$

(6)

referred to as the ‘Brunt actor’, from which we can recover the volume density dispersion itself, |$\sigma _{\rho /\rho _0}$|⁠.

2.3.2 Velocity statistics

The plane-of-sky dispersion of the velocity centroid can be used to estimate the line-of-sight (LOS) velocity variance, from which we can reconstruct the 3D turbulent velocity dispersion, defined as

$$\begin{eqnarray} \sigma _{v,\mathrm{3D}}= \left(\sigma _{v_x}^2 + \sigma _{v_y}^2 + \sigma _{v_z}^2\right)^{1/2}, \end{eqnarray}$$

(7)

or simply, the sum in quadrature of the velocity variance in each spatial direction x, y, z. Of course we cannot measure the variance directly from observations, so instead we assume that the fluctuations in the x- and y-directions are the same as the LOS fluctuations. Whilst this assumption may not be reasonable for a cloud trapped in an accelerating flow originating from the GC, we minimize any potential variance by removing systematic motions through our gradient and large-scale subtraction methods. Given the lack of spectral resolution in these data, we cannot directly measure the LOS dispersion (the second moment), and so we make a further assumption that the variance in velocity centroid (first moment) in the plane-of-sky is proportional to the variation along the LOS. To recover σ_v,3D we follow the methods developed by Stewart & Federrath (2022), who show that the 3D turbulent velocity dispersion can be recovered from PPV space using the standard deviation of the gradient-corrected moment-1 map together with a correction factor¹ of |$C_\mathrm{(c-grad)}^\mathrm{\, any}=3.3\pm 0.5$|⁠.

The sonic Mach number (⁠|$\mathcal {M}$|⁠) of the turbulent component of the velocity field is given by

$$\begin{eqnarray} \mathcal {M}= \frac{\sigma _{v,\mathrm{3D}}}{c_\mathrm{s}}, \end{eqnarray}$$

(8)

and the sound speed, c_s is defined as

$$\begin{eqnarray} c_\mathrm{s}= \left(\frac{\gamma k_\mathrm{B} T}{\mu m_\mathrm{H}}\right)^{1/2}, \end{eqnarray}$$

(9)

where γ = 5/3 is the adiabatic index, k_B is the Boltzmann constant, T is the gas temperature, μ = 1.4 is the mean particle weight of H i (Kauffmann et al. 2008), and m_H is the mass of a hydrogen atom.

2.3.3 Temperature and sound speed

Calculating the Mach number requires us to know the sound speed in the gas, as per equation (8), which in turn requires an estimate of the temperature of the gas. In these particular data, we are constrained by the coarse spectral resolution of 5.5 km s⁻¹, which means we cannot directly measure LOS velocity fluctuations smaller than this, and we cannot perform a Gaussian decomposition on the PPV cubes to estimate the fraction of cold gas in these clouds, and therefore estimate an average temperature for the emitting H i. In G23, we assumed a constant sound speed for the whole SMC, based on the assumption that most of the emitting gas was warm neutral medium (WNM). Here, we must consider that the cold gas fraction is likely much more dominant, as there is a relatively large fraction (⁠|$\sim 55~{{\ \rm per\ cent}}$|⁠) of molecular material in C1 (but none in C3) (N23). In this case, we use the location of the molecular gas to estimate the temperature variance of the H i across the clouds, assuming that where the molecular gas dominates the combined (molecular and atomic) column density of the cloud, the spatially corresponding H i is all cold neutral medium (CNM). In regions where there is H i without H₂, we assume a mixture of the three H i phases (e.g. Wolfire et al. 2003), with a significant percentage of the gas existing as unstable neutral medium (UNM).

To begin with, we determine how much molecular column density is in each pixel compared to the total column density (atomic plus molecular), i.e. the molecular fraction

$$\begin{eqnarray} f_{\mathrm{mol}} = \frac{N_{\mathrm{H}_{2}}}{N_{\mathrm{H I}} + N_{\mathrm{H}_{2}}}. \end{eqnarray}$$

(10)

Once we have created a map of f_mol (Fig. 2), we need to assign temperatures to each pixel based on the ratio of molecular to atomic gas. Noting that we always observe H i gas, we need to estimate the H i temperature, and we do so via the following equation:

$$\begin{eqnarray} T = T_{\mathrm{CNM}}\times f_{\mathrm{mol}} + (1-f_{\mathrm{mol}})\times \overline{T}_{\mathrm{H I}}, \end{eqnarray}$$

(11)

where the average H i temperature (⁠|$\overline{T}_{\mathrm{HI}}$|⁠) is given by

$$\begin{eqnarray} \overline{T}_{\mathrm{HI}} = \mathrm{f_{CNM}}\times T_{\mathrm{CNM}} + \mathrm{f_{UNM}}\times T_{\mathrm{UNM}} +\mathrm{f_{WNM}}\times T_{\mathrm{WNM}}, \end{eqnarray}$$

(12)

with the CNM, UNM, and WMN mass fractions and temperatures of each phase listed in Table 1. We use the best observational estimates we have for the mass fractions of each phase, taken from the 21-SPONGE survey (Murray et al. 2018). These values are derived from the solar-neighbourhood ISM, and although it is possible that the actual mass fractions of each phase differ from these estimates in the extra-planar clouds, we are working under the assumption that the cloud material was launched from the Galactic plane, resembling the H i phase composition of the local ISM. Estimates for the temperature range of each phase are taken from McClure-Griffiths et al. (2023), where theoretical and observational estimates for T_CNM, T_UNM, and T_WNM are collated.

Figure 2.

A map of the fraction of molecular column density, f_mol for C1. The colourbar shows the variation of f_mol where there is molecular material, although the amount of H₂ only varies between ∼0.95 and 1., meaning that the black regions regions are pure H i (f_mol = 0). The red contours shows H₂ column density at 5.0 × 10²¹ cm⁻². The scale bar shows 10 pc, assuming the physical distance to the clouds is 8.3 kpc. The arrow in each panel points towards the GC, which is a proxy for the opposite direction in which the Galactic wind is flowing.

Table 1.

Mass fractions and temperatures of CNM, UNM, and WNM.

	Phase	CNM	UNM	WNM
(1)	Fraction fphase	0.28	0.20	0.52
(2)	Minimum Tphase	25 K	250 K	4000 K
(3)	Maximum Tphase	250 K	4000 K	8000 K

Note. Row (1) is the mass fraction for each phase taken from the 21-SPONGE survey (Murray et al. 2018). Rows (2) and (3) are the minimum and maximum temperature for each phase, taken from McClure-Griffiths et al. (2023).

Equation (11) describes three cases: pixels where there is more H₂ than H i in the column are the CNM temperature, pixels which have a mixture of H₂ and H i have a temperature which is a proportional mixture of T_CNM and the |$\overline{T}_{\mathrm{HI}}$|⁠, and pixels where there is only H i have the |$\overline{T}_{\mathrm{HI}}$| temperature.

Lastly, we use the minimum and maximum temperatures in Table 1 to construct a minimum and maximum temperature map for C1 via equations (12) and (11), and then use equation (9) to convert them into maps of the sound speed. We do this because there are several orders of magnitude in the temperature variations, and because the sound speed goes as the square root of the temperature, we are then able to take the arithmetic mean of the minimum and maximum maps as our final map to be used in the analysis, as the sound speed spans a smaller range. Using the lower temperature limits, we find that |$\overline{T}_{\mathrm{HI}} \sim 2100$| K, and with the upper limits |$\overline{T}_{\mathrm{HI}} \sim 5000$| K.

The above prescription is applied to both clouds, but as there is no molecular material in C3, equation (11) collapses to a constant temperature of |$\overline{T}_{\mathrm{HI}}$| (and therefore sound speed) across the cloud.

3 RESULTS

In this section, we present the results of the turbulence analysis described in the previous section. We discuss how the mapped turbulence statistics compare to the physical structure and orientation of the clouds, as well as comparing the median values of the maps with the clouds’ global turbulence quantities.

3.1 Main analysis quantities

Our methods are designed to measure the volume density dispersion (⁠|$\sigma _{\rho /\rho _0}$|⁠, Section 2.3.1), the Mach number (⁠|$\mathcal {M}$|⁠, Section 2.3.2), and the turbulence driving parameter (b, equation 1) which we refer to as the main analysis quantities, maps of which are shown in Fig. 3. The first column shows that the variation of |$\sigma _{\rho /\rho _0}$| is more pronounced in C3, while C1 has a smoother distribution. C1 has overall higher values of |$\sigma _{\rho /\rho _0}$|⁠, with a slight correlation between regions of high |$\sigma _{\rho /\rho _0}$| and the position of the molecular gas clumps at the ‘head’ of the cloud, when compared to its H i-only ‘tails’. While we cannot directly map the density itself, higher |$\sigma _{\rho /\rho _0}$| means there is a higher chance of regions of large absolute density, possibly signifying a phase transition in the atomic material that is required for molecule formation and maintenance (Krumholz, McKee & Tumlinson 2009; Banda-Barragán et al. 2021; Girichidis et al. 2021).

Figure 3.

Maps of C1 (top) and C3 (bottom), of the spatially resolved analysis quantities. From left to right, we show the volume density dispersion (⁠|$\sigma _{\rho /\rho _0}$|⁠), the Mach number (⁠|$\mathcal {M}$|⁠), and the turbulence driving parameter (b). The orange circle in the bottom left-hand corner of each figure shows the FWHM of the kernel, while the small white circle inside the kernel. The black contours show the integrated H₂ intensity 5.0 × 10²¹ cm⁻² in C1. The scale bar shows 10 pc at a distance of 8.3 kpc, and the arrows point towards the GC. For comparison, the colourbars are constant across clouds.

The Mach number map of C1 shows that there is a trend towards higher |$\mathcal {M}$| on the outside edges of the cloud, while the Mach number in the inner region is lower. There is some suggestion of higher Mach numbers near the molecular gas, as is to be expected from our construction of the sound speed map, in keeping with colder gas being more supersonic than warm gas. In C3, we have a great deal of variation in Mach numbers, although there is no obvious trend in said variations in relation to the head or tail of the cloud. If we subscribe to the interpretation of C3 as having been entrained in the wind for much longer than C1 (N23), we would expect it to be far more fragmented and disorganized (Armillotta et al. 2016; Zhang et al. 2017; Banda-Barragán et al. 2019; Gronke & Oh 2020; Schneider et al. 2020). From the aforementioned simulation work, the initial entrainment phase is marked by a distinct head–tail cloud structure, later evolving into an extended structure aligned with the wind, but where turbulent properties are likely to be homogeneous across the cloud, since there is no net shear with the local environment. This is borne out in our maps, as there is no trend in variations of the three main analysis quantities with position or orientation within the cloud. Overall, C3 is more subsonic than C1, although with pockets of higher Mach. This is expected, as C3 has no detectable molecular gas, and as such, we expect the average gas temperature to be warmer in C3 than in C1, and hence have a lower average Mach number.

Finally, taking the ratio of columns 1 and 2 gives us maps of the turbulence driving parameter (column 3). In C1, we can see that the tail of the cloud is driven more solenoidally than the head. This could be due to shearing at the interface of the atomic cloud material and the ionized wind. Overall however, C1 is compressively driven, particularly in the regions containing molecular gas, as is anticipated given the compression of the head of the cloud by the wind. C3 displays an interesting amount of variation across the cloud between solenoidal, mixed, and compressive driving, but the driving does not seem to be correlated with the head or tail of the cloud, as with the other two main quantities. It is possible that the compressive borders between the clumpy structures in C3 are cloudlet–cloudlet collisions.

Fig. 4 shows the PDF of the Mach number, density dispersion, and driving parameter, as mapped in Fig. 3. We see that C1 is distributed above |$\sigma _{\rho /\rho _0}= 1.0$|⁠, whereas the median value in C3 is 1.01. In the second column, we see that all the Mach numbers in C1 are in the trans to supersonic regime, whereas C3 has a significant portion of the distribution below |$\mathcal {M}=1$|⁠. In the third column, we see that the distribution of the turbulence driving parameter is slightly narrower in C1 as compared to C3, but that in the latter distribution the median value is skewed lower by roughly 10 per cent. This reflects the larger amount of variation seen in the map of b for C3. C1 and C3 both have significant portions of their distributions above b = 1.0, which is the upper limit describing purely compressive driving. While this presents as ‘unphysical’, as discussed in G23, these spatial measurements of the driving parameter should be taken as an indication of which parts of the cloud are relatively more or less compressively driven. We discuss the uncertainty associated with our measurement of b quantitatively in Section 3.3.

Figure 4.

PDFs of the analysis quantities |$\sigma _{\rho /\rho _0}$|⁠, |$\mathcal {M}$|⁠, and b (left to right), for C1 (top) and C3 (bottom). The dashed line shows the median of each distribution, while the dotted lines show the 16th and 84th percentiles, respectively.

3.2 Intermediate analysis quantities

Fig. 5 shows the maps of the intermediate analysis quantities that are needed to construct the main analysis quantities, but that are none the less important in their own right. In the first column, we see that the Brunt factor (⁠|$\mathcal {R}^{1/2}$|⁠) is predominantly of order 0.3∼0.4 for both clouds, which is what we expect from previous findings (Brunt 2010; Ginsburg, Federrath & Darling 2013; Federrath et al. 2016; Menon et al. 2021; Sharda et al. 2022; Gerrard et al. 2023). Comparing the map of the Brunt factor to |$\sigma _{\rho /\rho _0}$| in Fig. 3, we see that the variations in |$\mathcal {R}^{1/2}$| do not dominate the variations in |$\sigma _{\rho /\rho _0}$|⁠, such that the volume density dispersion is primarily a consequence of the observed column density. We discuss the errors introduced by the Brunt method further in Section 3.3.

Figure 5.

Same as Fig. 3 but showing the Brunt factor (⁠|$\mathcal {R}^{1/2}$|⁠), the 3D velocity dispersion (σ_v,3D), and the sound speed c_s. Because there is no molecular gas in C3, the sound speed is a constant value everywhere in this map.

The variation in the 3D velocity field is quite smooth across C1, although it exhibits higher σ_v,3D towards the edges of the tail of the cloud, which is reflected in the Mach number map in Fig. 3. It appears that there is a slight trend towards lower σ_v,3D in the molecular regions. C3 has the lowest σ_v,3D values of the two clouds, with a similar amount of variation across the cloud as in its |$\sigma _{\rho /\rho _0}$| map, which is again reflected in the |$\mathcal {M}$| map. By construction, the sound speed map (column 3) of C1 has the lowest values closest to where the molecular clumps are. The map has been processed through the same kernel method as all the other quantities, which has smoothed out any steep transitions between high and low sound speeds. Of course the sound speed of C3 is a constant value as it has no molecular component, and as such our simple model for the temperature does not include any variation in this case.

Fig. 6 shows PDFs of the intermediate analysis quantities. C3 has a broader distribution of Brunt factors than C1, and an overall lower median value. C1 has a larger value of σ_v,3D than in C3, as well as a more asymmetrical distribution. Dividing σ_v,3D by the sound speed (third column) gives us the Mach number in Fig. 4. By construction, C3 has only one value for the sound speed, while C1 has a range of values but the bulk of the cloud is at the H i average sound speed (as outlined in Section 2.3.3).

Figure 6.

Same as Fig. 4 but showing the intermediate analysis |$\mathcal {R}^{1/2}$|⁠, σ_v,3D, and c_s.

3.3 Uncertainties

There are several sources of systematic error associated with our methods. The error associated with the Brunt method is quoted as |$\sim 10~{{\ \rm per\ cent}}$| in Brunt et al. (2010), when the Fourier image of the column density is axisymmetric. In Appendix  C, we discuss the level of anisotropy in each cloud and conclude that the upper limit on the error introduced via the Brunt method is 40 per cent in this instance.

Converting the centroid velocity variance to σ_v,3D via the Stuart method introduces an error of |$\sim 10~{{\ \rm per\ cent}}$| (Stewart & Federrath 2022), which is carried through to our estimates of the Mach number. Lastly, we must consider the error introduced by our temperature method. For this, we will take the upper and lower bound on the possible values that the temperature and therefore c_s can take depending on the temperature ranges of each H i phase, as described in Section 2.3.3. For C1, the mean of the sound speed map constructed using the maximum temperatures of each phase is 6.7 km s⁻¹, and the mean of the same map constructed with the minimum temperatures of each phase is 4.4 km s⁻¹. For C3 we have 7.0 and 4.6 km s⁻¹, respectively. Comparing these values to the sound speeds derived using the mean of the minimum and maximum temperature maps gives us a relative error on the global sound speed of both maps of |$\sim 20~{{\ \rm per\ cent}}$|⁠.

We collect the relative errors outlined above and propagate them through the global quantities for each cloud. This results in a relative error of |$\sim 22~{{\ \rm per\ cent}}$| for the Mach number, and |$\sim 46~{{\ \rm per\ cent}}$| for the driving parameter. Table 2 summarizes the bulk characteristics of each cloud, like H i mass, molecular mass, age, and distance from the galactic centre, all from N23. To add to this, we list the median values of the mapped quantities, with error bars representing the 16th and 84th, as well as the global values for each cloud which we compute directly on the LPF-corrected moment-0 and moment-1 maps. We find that when comparing the mapped and global quantities, the latter is lower than the former across all values in C1, but the opposite is true (with the exception of |$\mathcal {R}^{1/2}$|⁠) for C3. This is likely due to the size of the analysis kernel compared to the cloud in each case, as C3 is a much larger object and therefore the smallest resolvable kernel is small compared to its on-sky extent. This means that the scale probed by the mapped quantities is smaller by 1–2 orders of magnitude (depending on the cloud) than that probed by the global quantities. However, comparing the clouds and the scales on which the analysis quantities are calculated on, all the values agree within 1σ, except for σ_v,3D in C3. Theoretically, the velocity variance should scale with the length scale on which it is measured as a power law with an exponent of p ∼ 0.5 in the supersonic regime, and p ∼ 0.4 in the subsonic regime (e.g. Larson 1981; Ossenkopf & Mac Low 2002; Heyer & Brunt 2004; Federrath et al. 2021), so it is expected that the global quantities of the clouds are higher than median values measured on the kernel scale.

Table 2.

Summary of key quantities for C1 and C3. The top section of the table shows the kinematic wind model quantities (age and distance from the GC) from Di Teodoro et al. (2020) and Lockman, Di Teodoro & McClure-Griffiths (2020), assuming a distance to the GC of 8.2 kpc and the subsequent derived masses (H i and H₂ masses) from N23.The derived masses are calculated using an arbitrary boundary of each cloud, both spectrally and spatially as described in N23. The middle section shows the median analysis quantities, with error bars of the 16th and 84th percentiles of each distribution. The bottom section shows the analysis quantities calculated for the entire cloud, with error bars as outlined in Section 3.3.

	C1	C3	Reference
Kinematic model			Di Teodoro et al. (2020).;
quantities			Lockman, Di Teodoro & McClure-Griffiths (2020)
Distance (kpc)	0.8	2.1
Age (Myr)	3.5	8.3
Derived masses			N23
M|$_{\rm{\small HI}}\, (\mathrm{M}_{\odot })$|	385	3426
M|$_{\mathrm{H}_{2}}\, (\mathrm{M}_{\odot })$|	443	–
Mapped quantities			This work
|$\mathcal {R}^{1/2}$|	|$0.34^{+0.02}_{-0.02}$|	|$0.31^{+0.03}_{-0.05}$|
|$\sigma _{\rho /\rho _0}$|	|$1.4^{+0.2}_{-0.3}$|	|$1.0^{+0.2}_{-0.2}$|
σv, 3D (km s−1)	|$7.9^{+0.9}_{-1.5}$|	|$6.4^{+1.3}_{-1.5}$|
cs (km s−1)	|$5.8^{+0.6}_{-0.04}$|	|$5.8^{+0.0}_{-0.0}$|
|$\mathcal {M}$|	|$1.4^{+0.1}_{-0.2}$|	|$1.1^{+0.2}_{-0.3}$|
b	|$1.0^{+0.2}_{-0.3}$|	|$0.9^{+0.2}_{-0.3}$|
Global quantities			This work
|$\mathcal {R}^{1/2}$|	0.4 ± 0.1	0.2 ± 0.1
|$\sigma _{\rho /\rho _0}$|	1.5 ± 0.6	1.8 ± 0.7
σv,3D (km s−1)	9.6 ± 1.0	11.0 ± 1.1
cs (km s−1)	5.5 ± 1.1	5.8 ± 1.2
|$\mathcal {M}$|	1.7 ± 0.4	2.0 ± 0.4
b	0.9 ± 0.4	0.9 ± 0.4

3.4 Context

As outlined in the introduction, there have been several measurements of the turbulence driving parameter from observational data, mostly in star-forming molecular clouds. Fig. 7 summarizes these previous measurements, with the addition of C1 and C3, showing that these clouds sit between previous measurements of the diffuse H i and the dense molecular clouds on the |$\sigma _{\rho /\rho _0}-\mathcal {M}$| diagram. As we have shown, C3 has a lower Mach number and density dispersion, and so sits below C1 in the sub-to-transonic regime on this figure. More high-resolution observations of these extra-planar H i clouds will allow us to further populate this figure with these kind of objects, but from this investigation it seems that H i clouds of this type have lower |$\mathcal {M}$| and |$\sigma _{\rho /\rho _0}$| values than their molecular counterparts, although the ratio still results in high levels of compressive driving which is comparable to the molecular pillars in NGC 3372 (Menon et al. 2021).

Figure 7.

A summary of the available observational estimates for the density dispersion–Mach number relation in different environments. The y-axis shows the volume density dispersion (⁠|$\sigma _{\rho /\rho _0}$|⁠), and the x-axis shows the turbulent Mach number (⁠|$\mathcal {M}$|⁠), including a factor involving plasma β (ratio of thermal to magnetic pressure), as some of the literature values shown (Taurus and G0.253+0.016) have been calculated using the magnetized version of the |$\sigma _{\rho /\rho _0}-\mathcal {M}$| relation in equation (1) (see Molina et al. 2012b). The three diagonal lines show the theoretical limits for compressive (b = 1.0, dotted), mixed (b = 0.38, solid), and solenoidal (b = 0.33, dashed) driving of the turbulence (Federrath et al. 2010). The hexagons show the median values for C1 (red) and C3 (orange). The error bars on these points show the 16th to 84th percentile on each axis. For context, we include a variety of sources from the literature: Taurus (dark blue star; Brunt 2010), which includes magnetic field estimates and revised Mach number estimations from Kainulainen & Tan (2013), using ¹³CO line imaging observations; IC5146 (blue cross; Padoan et al. 1997), using ¹²CO and ¹³CO observations; GRSMC 43.30−0.33 (aqua plus; Ginsburg et al. 2013), observed in H₂CO absorption and ¹³CO emission; ‘The Brick’ (G0.253+0.016, teal square; Federrath et al. 2016), using HNCO observations; ‘The Pillars of Creation’ (NGC 3372 pillars, magenta diamonds; Menon et al. 2021), from ¹²CO, ¹³CO, and C¹⁸O; ‘The Papillon Nebula’ (LMC N159E, pink circle; Sharda et al. 2022), again in ¹²CO, ¹³CO, and C¹⁸O; the WNM in the MW (lilac triangle; Marchal & Miville-Deschênes 2021) (H i observations); and the WNM in the SMC (purple triangle) (G23).

4 SUMMARY

In this study, we employ the G23 methods for calculating the properties of the internal turbulence of two extra-planar H i clouds, and develop a novel approach for modelling the spatially varying temperature of a mixture of atomic and molecular hydrogen. Cloud C1 contains a significant amount (∼55 per cent by mass) of molecular material, is small and compact, at b ∼ −3.9° below the disc. Cloud C3 is much larger and more spatially dispersed, it contains no detectable molecular material, and is at a latitude of b ∼ 7.0° above the disc. Despite their physical differences and distance from the plane of the Milky Way, the internal turbulence statistics of these clouds are very similar. We find that both clouds are in the sub-to-trans-sonic regime with |$\mathcal {M}~\sim 1.4^{+0.1}_{-0.2}$| in C1 and |$\mathcal {M}~\sim 1.1^{+0.2}_{-0.3}$| in C3. We find that the turbulence in both clouds is predominately compressively-driven with |$b ~\sim 1.0^{+0.2}_{-0.3}$| in C1 and |$b ~\sim 0.9^{+0.2}_{-0.3}$| in C3. The volume density contrast in C1 is |$\sigma _{\rho /\rho _0}\sim 1.4^{+0.2}_{-0.3}$| and in C3 |$\sigma _{\rho /\rho _0}\sim 1.0^{+0.2}_{-0.2}$|⁠. Comparing the median of the spatially varying analysis quantities and the global quantities for each cloud, there is good 1σ agreement across all quantities and clouds. We find that there are no obvious trends in the spatial variation of the analysis quantities with orientation or position of C3, but C1 shows signs of the cloud-wind interaction at the head of the cloud leading to more compressively-driven turbulence in that region. This is consistent with the idea that the wind compressing the tip of the clouds drives compressive turbulence, but that this effect of the cloud-wind interaction may become less focused as the cloud disperses during its journey away from the Galactic disc.

ACKNOWLEDGEMENTS

The authors acknowledge Interstellar Institute’s programme ‘II6’ and the Paris-Saclay University’s Institut Pascal for hosting discussions that nourished the development of the ideas behind this work. The MeerKAT telescope is operated by the South African Radio Astronomy Observatory, which is a facility of the National Research Foundation, an agency of the Department of Science and Innovation. IAG would like to thank the Australian Government and the financial support provided by the Australian Postgraduate Award. CF acknowledges funding by the Australian Research Council (Discovery Projects grant DP230102280), and the Australia-Germany Joint Research Cooperation Scheme (UA-DAAD). CF further acknowledges high-performance computing resources provided by the Leibniz Rechenzentrum and the Gauss Centre for Supercomputing (grants pr32lo, pr48pi, and GCS Large-scale project 10391), the Australian National Computational Infrastructure (grant ek9) and the Pawsey Supercomputing Centre (project pawsey0810) in the framework of the National Computational Merit Allocation Scheme and the ANU Merit Allocation Scheme, through which the data analyses presented in this paper were performed. EDT was supported by the European Research Council (ERC) under grant agreement no. 10104075. This research was partially funded by the Australian Government through an Australian Research Council Australian Laureate Fellowship (project number FL210100039) to NMcG.

DATA AVAILABILITY

The data underlying this article along with a general implementation of the code used to process the data cubes is available via zenodo at https://doi.org/10.5281/zenodo.8060960.

Footnotes

1

This correction factor is the mean of the p₀ values in lines 4–6 of table E1 of Stewart & Federrath (2022). We choose the gradient-subtracted statistics and choose the mean of those values, which are independent of the LOS orientation with respect to the rotation axis of the cloud/kernel region.

APPENDIX A: KERNEL SIZE

The size of the kernel FWHM is set as the diameter of a circle with an area 3 times the minimum viable area required to recover the statistics in each kernel window. Limiting the size of the kernel as much as possible so that we recover as much spatial detail as possible must be balanced with making the kernel large enough so that the edges are not eaten away during the processing, and we recover as much of the area covered by the raw input data as possible. We explored this for C1, as shown in Fig. A1, comparing the area covered by the input and output maps, depending on how large we made the kernel in comparison to the minimum viable area. We see that by making the kernel 3 times as large as the minimum area we recover 95 per cent of the input data, and that larger areas don’t drastically increase the amount of area recovered. For each cloud, therefore, we make the FWHM of the kernel such that the area contained within it is 3 times the minimum viable area, which is to say that for any given window, at least a third of the area inside the kernel FWHM must be viable data. This is all ultimately set by the beam size, as we are ensuring that at least 20 beams worth of area is contained in this region. The major axis of the beam for each PPV cube is about ∼25 arcsec.

Figure A1.

The percentage area of the raw data recovered by the pipeline for C1. Using a multiplying factor of 3 recovers 95 per cent of the raw data, which is what we have chosen to use throughout the analysis to set the kernel FWHM.

APPENDIX B: LOW-PASS FILTERING METHOD

As discussed in Section 2.2.2. to extract 3D density statistics from the column density maps, we first subtract the low-pass filter from each map, such that the resulting column density map approximates a log-normal distribution. The results of this are shown in Fig. B1 and B2.

Figure B1.

C1 column density. The left-hand panel shows the log of the column density map prior to background subtraction above, and a PDF of the column density below. The middle panel shows the map of the low-pass background filter and its PDF, while the right-hand panel shows the background-subtracted map. On each PDF, a solid green line shows the fitted Gaussian, along with its σ value. The circle in the top right-hand corner of each panel represents the FWHM of the convolution kernel.

Figure B2.

Same as Fig. B1, but for C3.

We apply the LPF method to the centroid velocity maps in the same way as the column density, which can be seen in Figs B3 and B2.

Figure B3.

C1 centroid velocity. The panels represent the same steps as in Fig. B1, and the overlays are the same.

Figure B4.

Same as Fig. B3, but for C3.

APPENDIX C: BRUNT METHOD AND ANISOTROPY IN COLUMN DENSITY POWER SPECTRA

The Brunt method, described in Section 2.3.1, relied on the assumption of the column density fluctuations being isotropically distributed in k-space. To ascertain the amount of uncertainty introduced by using this method to obtain |$\sigma _{\rho /\rho _0}$|⁠, we fit ellipses to the Fourier image of C1 and C3 (Fig. C1) to investigate the symmetry in the 2D power spectra. In C1, the largest ellipticity is e = 1.47, while in C3 it is e = 1.21, which translates to a 47 per cent deviation from circular in C1 and a 21 per cent deviation in C3. This deviation can be used as a proxy for the amount of anisotropy in the images. Following Federrath et al. (2016), we conclude that the error associated with the conversion of the 2D-to-3D density contrast conversion is of order <40 per cent.

Figure C1.

2D power spectra of the normalized column density of the two clouds. Left: C1 with contours fitted at 10⁻⁵, 10⁻⁶, 10⁻⁷ corresponding to ellipticities of |$e=1.04,\, 1.17,\, 1.20$|⁠. Right: C3 with contours fitted at 10⁻⁴, 10⁻⁵, 10⁻⁶ corresponding ellipticities of |$e=1.24,\, 1.47,\, 1.00$|⁠.