H i Vertical Structure of Nearby Edge-on Galaxies from CHANG-ES

The vertical structure of H i gas disks is an important tracer of the galactic potential and dynamical effects in spiral galaxies. Under the assumption of hydrostatic equilibrium, the H i gas vertical structure exhibits a force balance between the self-gravity of the disk and the effective pressure of H i (Boulares & Cox 1990; Piontek & Ostriker 2007; Koyama & Ostriker 2009; Ostriker et al. 2010; Krumholz et al. 2018). The effective pressure, largely supported by the turbulence (Mac Low 1999; Tamburro et al. 2009), further reflects an energy balance between the radiative dissipation and the energy input that comes from stellar feedback and galactic-scale gas inflows driven by non-axisymmetric torques (Krumholz & Burkert 2010; Forbes et al. 2012, 2014). Real gas disks can be additionally influenced by gas accretion (Kereš et al. 2005; Dekel & Birnboim 2006; Oosterloo et al. 2007; Dekel et al. 2009; Zheng et al. 2017), and tidal interactions (Toomre & Toomre 1972; Barnes & Hernquist 1992; Di Matteo et al. 2007), as well as energetic feedback from active nuclei (Sijacki et al. 2007; Fabian 2012; Cicone et al. 2014).

Observationally, the vertical structure of an H i disk is quantified by the scale height to first order (Randriamampandry et al. 2021), which reflects properties consistent with the hydrostatic equilibrium model. The H i scale heights are nearly constant in the inner region of the galaxy and increase appreciably as a function of radius (the so-called flaring phenomenon) in the outer region between 5 and 35 kpc (Dickey & Lockman 1990; Kalberla & Kerp 2009). Narayan & Jog (2002) explained this radial trend of the Galactic H i scale heights and emphasized the importance of gravitational coupling between gas, stars and dark matter. They concluded that in addition, to the contributions from stars, gravity from the atomic and molecular gas helps set the scale heights in the inner region, while the atomic gas is more important than the molecular gas at intermediate radius; at large Galactic radius, the H i scale heights flare because the gravitational force decreases quickly and the dark matter halo dominates the gravity. The sharpness of the H i flares at large galactic radius is explained by the truncation of the stellar disks (Van Der Kruit 1988). H i flares are found to be common in spiral galaxies (Brinks & Burton 1984; Bigiel & Blitz 2012), though they are not as sharp as in the Milky Way, as truncated stellar disks are not observed in all galaxies (Bland-Hawthorn et al.2005).

The vertical hydrostatic equilibrium model of the H i disk serves as a useful tool to derive properties involved in this equilibrium. Typically, for face-on galaxies, observational studies easily obtain the gas velocity dispersion but rely on the equilibrium model when deriving the gas disk thickness, while for edge-on galaxies, the situation is the other way round. Olling (1995) and Narayan et al. (2005) developed a model using the flare of H i to constrain the shape, mass and size of the dark matter halo. Krumholz et al. (2018) identified the isothermal nature of gas based on the vertical force balance between the gravitational drag force of gas, stars and dark matter and the gas thermal, turbulent and magnetic pressure. Recent developments include Bacchini et al. (2019) deriving the H i volumetric density to study the volumetric star formation law. They found that the volumetric star formation law is much tighter than the surface star formation law, particularly, the previously known break in the slope of the surface star formation law (Bigiel et al. 2008; Leroy et al. 2008) is more likely due to disk flaring rather than a decrease of the star-forming efficiency at low surface densities (Bacchini et al.2019).

The vertical structure of the H i disk is more complex than described by this simple hydrostatic equilibrium model. A thick H i layer, also called extraplanar gas or H i halo in the literature, with a typical 1–2 kpc scale height (in contrast to 100–200 pc for the thin disk) is directly observed in edge-on galaxies, with a lag in rotation with respect to the thin H i disk (Oosterloo et al. 2007; Kamphuis et al. 2013). The Galactic intermediate-velocity clouds can be viewed as a form of extraplanar gas (Oort 1970; Wakker & van Woerden 1997). The extraplanar gas is also modeled and found to be prevalent in external galaxies that are not edge-on (Fraternali et al. 2001). Marasco et al. (2019) did a systematic study of extraplanar gas in 15 nearby late-type galaxies and concluded that both the mass and kinematics of extraplanar gas are in good agreement with the galactic fountain model which is powered by stellar feedback (Shapiro & Field 1976; Bregman 1980). In the fountain model, gas is pushed away from the disk by the supernova feedback and falls back with extra gas from the circumgalactic medium (CGM) after metal enrichment.

More complexities come from the fact that the edge-on views of H i disks are not necessarily flat but warped outside the edge of optical disks (Burke 1957; Kerr 1957; Sancisi 1976; Newton & Emerson 1977; Bosma 1978; Briggs 1990). Warps are found to be ubiquitous in disk galaxies (Sancisi 1976; García-Ruiz et al. 2002). The fact that they usually onset at the edge of optical disks suggests that the inner flat disk and the outer warped disk have different formation histories and probably involve different epochs (Van der Kruit 2007). The mainstream explanation for the formation of warps seems to be the accretion of material with an angular momentum vector misaligned with that of the main disk (Jiang & Binney 1999; Shen & Sellwood 2006; Roškar et al. 2010). Other explanations include the torques from the misaligned inner disk and the associated inner oblate halo (Dekel & Shlosman 1983; Toomre 1983; Sparke & Casertano 1988), and the tidal force from companion galaxies (Weinberg 1995; Weinberg & Blitz 2006).

An observational census would usefully gain more insights into the physics that shape the vertical extent of H i disks. We use the highly inclined galaxies from the Continuum Halos in Nearby Galaxies—an EVLA Survey (CHANG-ES) (Irwin et al. 2012; Wiegert et al. 2015; Irwin et al. 2019). The H i information of 19 edge-on galaxies in CHANG-ES is published in Zheng et al. (2022). Two galaxies are for the first time presented in H i interferometric images and twelve of the galaxies have better H i spatial resolutions and/or sensitivities of intensity maps than literature. The H i data are not well resolved kinematically, but we manage to derive the radially averaged scale height for the 15 most edge-on galaxies (inclination > 80°, Section 3) with a similar method as another CHANG-ES study Krause et al. (2018) used to measure the scale height of continuum halos. In the literature, the scale heights have been derived with sophisticated kinematical modelings based on data of much higher spectral resolutions, particularly when the galaxies are less inclined (Yim et al. 2014, 2020). Most of those modeling methods need to assume a quasi-equilibrium between the gravity and the gas pressure, and thus are most accurate for unperturbed galaxies. The procedure to derive the scale heights in this paper is mostly photometric using images, but we carefully correct the measurements for several types of observational artifacts, also called systematic biases or systematic broadening effects, including point spread function (PSF) smearing, planar projection and edge-on projection (Sections 3.2 and 3.3). Thus the procedure and measurements presented in this paper provide an alternative and convenient way to derive the scale heights without significant assumptions for the dynamic states of gas. It can easily be applied to statistical samples in the future, and may potentially provide an indicator for the dynamic states of H i gas when compared to results or expectations from kinematically derived scale heights. We investigate possible dependence of the H i scale height on other galactic properties (Section 4), and also more closely discuss the uncertainties and effects of the external environment (Section 5). We assume a ΛCDM cosmology (H₀ = 73 km s⁻¹ Mpc⁻¹) and Kroupa initial mass function (Kroupa 2001).

2.1. Sample and HI Data

2.2. SFR and Mass Densities

We take the star formation rate (SFR) from Vargas et al. (2019), who calculated SFR based on fluxes from narrow-band Hα images and fluxes from the Wide-field Infrared Survey Explorer 22 μm images. Because NGC 5084 lacks Hα observation, we exclude this galaxy from SFR related analysis. The SFR values are listed in Table 2, in which the SFR of NGC 5084 is only estimated by 22 μm data. The star formation surface density is estimated as ${{\rm{\Sigma }}}_{\mathrm{SFR}}=\mathrm{SFR}/(\pi * {R}_{25}^{2})$.

Table 2. The H i Scale Height and Other Physical Parameters

Galaxies	${\bar{h}}_{\mathrm{phot}}$	$\bar{h}$	$\bar{h}$	${\bar{h}}_{{\rm{Z}}21}^{a}$	${{\rm{\Sigma }}}_{\mathrm{Mtot},{\rm{r}}25}^{b}$	${{\rm{\Sigma }}}_{\mathrm{Mbaryon},{\rm{r}}25}^{c}$	SFR${}_{{{\rm{H}}}_{\alpha }+22\mu {\rm{m}}}^{d}$
	('')	('')	(kpc)	(kpc)	(107 M⊙ kpc−2)	(107 M⊙ kpc−2)	(M⊙ yr−1)
NGC 2683	22.64 ± 4.60	9.69 ± 2.18	0.29 ± 0.07	0.79 ± 0.27	29.84	6.70	0.25 ± 0.03
NGC 3003	26.77 ± 11.53	13.00 ± 5.94	1.60 ± 0.73	4.33 ± 2.08	3.64	1.69	1.56 ± 0.16
NGC 3044	11.43 ± 0.94	4.00 ± 0.63	0.39 ± 0.06	1.0 ± 0.28	12.18	3.11	1.75 ± 0.16
NGC 3079	20.26 ± 6.59	6.82 ± 2.52	0.68 ± 0.25	2.35 ± 1.15	7.01	4.58	5.08 ± 0.45
NGC 3556	25.82 ± 5.29	4.91 ± 1.12	0.34 ± 0.08	1.89 ± 0.45	7.51	5.15	3.57 ± 0.3
NGC 3877	15.90 ± 1.69	7.56 ± 1.00	0.65 ± 0.09	1.74 ± 0.2	10.15	5.77	1.35 ± 0.12
NGC 4096	23.93 ± 4.43	8.28 ± 1.71	0.41 ± 0.09	1.48 ± 0.28	13.42	3.82	0.71 ± 0.08
NGC 4157	17.51 ± 1.70	5.02 ± 0.64	0.38 ± 0.05	1.23 ± 0.28	12.62	4.90	1.76 ± 0.18
NGC 4217	9.85 ± 0.67	3.66 ± 0.56	0.37 ± 0.06	0.93 ± 0.17	14.61	5.15	1.89 ± 0.18
NGC 4302	9.89 ± 1.16	6.03 ± 1.76	0.57 ± 0.17	0.96 ± 0.21	14.36	4.16	0.92 ± 0.08
NGC 4565	16.76 ± 2.51	9.56 ± 1.76	0.55 ± 0.10	1.27 ± 0.22	11.15	3.03	0.96 ± 0.09
NGC 4631	18.66 ± 4.01	17.85 ± 4.45	0.64 ± 0.16	0.78 ± 0.27	6.99	2.47	2.62 ± 0.22
NGC 5084	21.45 ± 1.89	18.47 ± 2.43	2.10 ± 0.28	2.95 ± 0.35	14.30	1.98	0.1 ± 0.03*
NGC 5775	15.17 ± 1.86	3.14 ± 0.75	0.44 ± 0.11	2.01 ± 0.56	14.57	9.31	7.56 ± 0.65
UGC 10288	11.35 ± 0.29	7.19 ± 0.51	1.19 ± 0.09	1.91 ± 0.09	6.33	1.87	0.66 ± 0.07

Notes.

^a The Gaussian averaged scale height corrected using the method in Z21. ^b Total mass surface density ${{\rm{\Sigma }}}_{\mathrm{Mtot},{\rm{r}}25}={M}_{\mathrm{tot}}/(\pi * {R}_{25}^{2})$, where M_tot (Irwin et al. 2012) is the total mass within the blue isophotal radius R₂₅ in kpc, calculated using the velocity amplitude based on line width. ^c Baryon mass surface density ${{\rm{\Sigma }}}_{\mathrm{Mbaryon},{\rm{r}}25}={M}_{\mathrm{baryon}}/(\pi * {R}_{25}^{2})$, where M_baryon = 1.4 × M_HI + M_* is the baryon mass within the blue isophotal radius R₂₅ in kpc. ^d Star formation rates were estimated by combination of Hα and 22 μm data, taken from Paper XVII (Vargas et al. 2019), except NGC 5084. Lacking Hα observation, the star formation rate of NGC 5084 is only estimated by 22 μm data, taken from Paper IV (Wiegert et al. 2015).

Download table as:  ASCIITypeset image

The total (or dynamic) mass surface density within the optical radius Σ_Mtot,r25 is taken from Paper IX (Krause et al. 2018), who calculated it based on the inclination corrected line widths of the H i spectrum.

We also derive the baryonic mass surface density within R₂₅, Σ_Mbaryon,r25. We take the stellar mass M_* from Z21, who derived it based on the i-band luminosity and g − i color correlated mass-to-light ratio (Bell et al. 2003). We further measure the H i mass within R₂₅, M_HI,r25. The baryonic mass is calculated as the sum of the stellar mass and gas mass within R₂₅, M_baryon,r25 = M_*,r25 + 1.4M_HI,r25, where 1.4 is the standard correction factor to account for helium and metals. Then the baryonic mass surface density is calculated as ${{\rm{\Sigma }}}_{\mathrm{baryon}}\,={M}_{\mathrm{baryon},{\rm{r}}25}/(\pi * {R}_{25}^{2})$. These surface densities are listed in Table 2.

We use a photometric approach instead of kinematic approach to determine scale heights due to two characteristics of our sample. First the accuracy of three-dimensional (3D) kinematic modeling routines is challenged by highly inclined observations. Second, the low-velocity resolution of the CHANG-ES H i cubes does not permit accurate modeling.

We first derive the raw photometric H i scale height (h_phot) in a similar way as in Z21 (see Section 3.1), but skip their relatively simple correction for beam smearing and planar projection (caused by the not perfectly edge-on inclinations of galaxies). Instead, we investigate in detail and design a more robust procedure to correct for systematic biases including these two effects in Z21 and additional edge-on projection effect (mainly caused by flaring).

3.1. Deriving the Raw Photometric Scale Heights

The raw photometric scale heights h_phot are derived in two steps: deriving the vertical profile of surface densities in strips perpendicular to the disk midplane, and fitting the vertical distribution profiles with a Gaussian function to obtain the width. ⁸

These steps are similar to the "BoxModels" task in the new NOD3 program package (Müller et al. 2017). Krause et al. (2018) applied the "BoxModels" task to measure the radio halo scale height of the CHANG-ES sample. However, the H i disks are thin and asymmetric compared with the radio halos. Z21 thus designed additional steps specifically for measuring H i scale heights. They derived the PA of the H i plane instead of using the optical PA. They also measured the vertical distributions on both sides separately.

As in Z21, we only use the radially averaged h_phot within the optical radius R₂₅ (${\bar{h}}_{\mathrm{phot}}$) in the analysis. This is because the nearly edge-on view prevents us from recovering the flaring shape of the radial profile of scale heights. ${\bar{h}}_{\mathrm{phot}}$ should thus be only viewed as an indicator of the H i disk thickness. The error of each ${\bar{h}}_{\mathrm{phot}}$ combines the uncertainty from the model fitting and the scatter of h_phot in the profile.

Figure 1 shows the radial profile of h_phot and the value of ${\bar{h}}_{\mathrm{phot}}$ (black dashed line) in each galaxy. We present four profiles in each galaxy. The fitting results of the "up" and "low" represent two sides with respect to the midplane of a disk, while the "left" and "right" are two sides with respect to the minor axis of a disk. The "up" and "low" profiles here correspond to (but are not the same, as the deriving methods are different) the red and blue fitting curves in Figure A2 of Z21. We show these profiles mainly to demonstrate how ${\bar{h}}_{\mathrm{phot}}$ values are derived, but emphasize that only ${\bar{h}}_{\mathrm{phot}}$ values are considered reliable in the following analysis.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

3.2. New Procedure to Correct for Artificial Broadening

3.3. Investigation of Scale Height Broadening Based on Mock Cubes

We use mock cubes to quantify the extent of overestimating the scale heights due to effects of PSF smearing, planar projection and edge-on projection. The errors are calculated as 1σ.

We generate mock H i cubes with the GALMOD task of BBarolo (Di Teodoro & Fraternali 2015), using the H i scale height profiles derived in Bacchini et al. (2019) for eight THINGS galaxies, along with the velocity dispersion profiles, surface density profiles and rotation curves in that paper. We have excluded two galaxies from the original sample of Bacchini et al. (2019), DDO 154 and NGC 2976, because of their low stellar masses $\mathrm{log}{M}_{* }=7.1$ and 9.1, while the lowest stellar mass of the CHANG-ES edge-on sample is $\mathrm{log}{M}_{* }=9.96$. We call these eight galaxies the input sample.

The inclination i and maximum radius of ring ${r}_{\max }$ (larger than R₂₅ in all cases) are fixed for each model disk, but vary between model disks. Each model disk is built with rings of different radius r. For each r, the H i ring has scale height h_r , rotation velocity v_r , velocity dispersion σ_r and surface density Σ_r , with values determined by interpolating the related radial profiles of a galaxy from the input sample. The last data point in a radial profile is repeated if ${r}_{\max }$ exceeds that. The cube and resulting moment-0 map are first generated at the original THINGS resolution (i.e., the highest resolution) with GALMOD. Later, the spatial resolution of the output moment-0 maps can be changed through convolving with Gaussian kernels. We add Gaussian noise with sigma equivalent to the median rms of the CHANG-ES HI moment-0 maps to the model moment-0 maps. Similar to how we treat the CHANG-ES data, when deriving h_phot, we only utilize the data points with signal to noise ratio larger than 3 along each vertical strip to fit Gaussian models. We use the scale height averaged within the radius R₂₅, $\bar{h}$, as the parameter to be tested. We compare the photometrically measured $\bar{h}$ (${\bar{h}}_{\mathrm{phot}}$) to the real $\bar{h}$ of models (${\bar{h}}_{\mathrm{true}}$).

In the following, we first investigate separately each of the effects causing systematic biases, including PSF smearing, planar projection and edge-on projection. Then we take these effects together and correct ${\bar{h}}_{\mathrm{phot}}$ against ${\bar{h}}_{\mathrm{true}}$.

3.3.1. Effect of PSF Smearing

For the test of the PSF smearing effect, we fix the inclination of the mock galaxies to 90° and vary the beam's major axis (b_maj ) of the data cube from the highest resolution of b_maj = 15 (pixel size of the THINGS data) to the lowest resolution corresponding to the minimum value of the uncorrected ${\bar{h}}_{\mathrm{phot}}/{b}_{\mathrm{maj}}$ in the CHANG-ES edge-on sample, with a step of Δb_maj = 15. We generate 117 mock H i cubes. We quantify, with decreasing resolution, how ${\bar{h}}_{\mathrm{phot}}$ is increasingly overestimated with respect to the measurement at the best resolution ${\bar{h}}_{0}$. The result is shown in Figure 2. The median trend of the overestimation has a relatively small scatter of 1.85% and is fitted with the following equation

$$\begin{eqnarray}&&{\bar{h}}_{\mathrm{phot}}^{2}-{\bar{h}}_{0}^{2}=0.234\ {b}_{\mathrm{maj}}^{2},\end{eqnarray} \tag{ 1 }$$

which could be approximated as a convolution with an effective Gaussian kernel. When the beam is less than 25% of ${\bar{h}}_{\mathrm{phot}}$, the ${\bar{h}}_{\mathrm{phot}}$ will have an uncertainty due to PSF smearing of 0.19%.

3.3.2. Effect of Planar Projection

In the following, we investigate the planar projection effect, which however is likely to interfere with the edge-on projection effect. We use the highest resolution of 15, and vary the inclination of the galaxy disk from 80° to 90° with a step of 1°. We also vary the disk sizes in the mocks, r_max. The range of r_max is set such that ${r}_{\max }/{R}_{25}$ values are between 1 and maximum R_HI/R₂₅ of the CHANG-ES H i sample. We build 21 mock cubes for each inclination, thus in total 231 mock cubes.

We quantify, with decreasing inclination angle, how ${\bar{h}}_{\mathrm{phot}}$ is increasingly overestimated with respect to the measurement when the galaxy is perfectly edge-on ${\bar{h}}_{90}$. Figure 3(a) shows the result of the 231 mocks (points), with the median distribution (red line). At a given inclination, the planar projection effect has large scatter varying between different input galaxies. As we explain and justify below, the edge-on projection is likely responsible for this large scatter.

Zoom In Zoom Out Reset image size

3.3.3. Effect of Edge-On Projection

The edge-on projection is mainly caused by disk flaring. It interferes with the planar projection because each line of sight intercepts different (and more) disk rings in a flared disk from in a flat disk.

To parameterize the shape of underlying H i flaring, for each galaxy in the input sample, we fit the scale height radial profile with a 3rd order polynomial equation h(r) = a_h r³ +b_h r² + c_h r + d_h . The four coefficients of the polynomial equation, together with r_max, are parameters to quantify the effect of edge-on projection. We define ${\rm{y}}^{\prime} =({\bar{h}}_{\mathrm{phot}}\,-{\bar{h}}_{90})/{\bar{h}}_{\mathrm{phot}}$, so that the scatter of data points from the median curve (the red line in Figure 3(a)) is quantified as ${\rm{\Delta }}y^{\prime} =y^{\prime} -\mathrm{median}(y^{\prime} )$. We fit a linear relation of the five parameters describing the edge-on projection to predict ${\rm{\Delta }}y^{\prime}$ (dashed line in Figure 3(b)). The strong linear correlation with a small scatter of 5.32% implies that the effect of edge-on projection indeed accounts for a large fraction of the scatter around the median relation in Figure 3(a).

The analysis in Figures 3(a) and (b) demonstrates that the edge-on projection effect is indeed involved in the planar projection effect. The underlying flaring shape and r_max used to fit ${\rm{\Delta }}y^{\prime}$ in 3(b) are inputs of mocks which cannot be directly measured in observations. Motivated by the quasi-static equilibrium model of gas (based on which Bacchini et al. 2019 derived the scale heights of the input sample), we consider M_*, M_HI, R_HI, R₂₅ and SFR as candidate parameters that may mimic the combined effect of the underlying shape of H i flaring and r_max. We test different combinations of these candidate parameters and find that the combination of M_*, R_HI and R₂₅ is enough to derive a similarly tight relation with ${\rm{\Delta }}y^{\prime}$ as when using r_max and the four coefficients describing the underlying shape of H i flaring.

In addition to interfering with planar projection, the edge-on projection further causes a difference between the perfectly edge-on ${\bar{h}}_{90}$ and model input ${\bar{h}}_{\mathrm{true}}$. This is because in an edge-on view, the flaring outer disk can dominate the surface brightness at high z, where z is the vertical distance from the midplane. Thus, in the following, we consider the planar projection and edge-on projection effects together, and correct ${\bar{h}}_{\mathrm{phot}}$ to ${\bar{h}}_{\mathrm{true}}$.

3.3.4. Correcting for Both Planar and Edge-on Projection Effects Together

We correct two projection effects from ${\bar{h}}_{\mathrm{phot}}$ to derive ${\bar{h}}_{\mathrm{true}}$. We do it in a similar way as in the last section, with two major modifications: ${\bar{h}}_{90}$ is replaced by ${\bar{h}}_{\mathrm{true}}$, and the parameters describing the flaring profile shapes are replaced by M_*, R_HI and M_HI.

The median relation between $({\bar{h}}_{\mathrm{phot}}-{\bar{h}}_{\mathrm{true}})/{\bar{h}}_{\mathrm{phot}}$ and $\cos i$ is derived to quantify the effect of planar projection (Figure 3(c)). The median trend has a scatter of 10.04%. Interestingly, the scatter looks much smaller than in Figure 3(a), possibly because different systematic biases cancel out. A 3rd order polynomial equation is fitted to the median curve

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{{\bar{h}}_{\mathrm{phot}}-{\bar{h}}_{\mathrm{true}}}{{\bar{h}}_{\mathrm{phot}}}={c}_{3}{(\cos \ i)}^{3}\\ \quad +\,{c}_{2}{(\cos \ i)}^{2}+{c}_{1}(\cos \ i)+{c}_{0},\end{array}\end{eqnarray} \tag{ 2 }$$

in which c₃ = −97.91, c₂ = 7.15, c₁ = 5.36 and c₀ = 0.18.

We calculate the offset of $({\bar{h}}_{\mathrm{phot}}-{\bar{h}}_{\mathrm{true}})/{\bar{h}}_{\mathrm{phot}}$ from the prediction of the best-fit relation, and denote it as Δy. Finally, we fit a linear relation of Δy as a function of M_*, R_HI and R₂₅.

$$\begin{eqnarray}\,\begin{array}{rcl}{\rm{\Delta }}y & = & 0.28\,\mathrm{log}{M}_{* }+0.68\,\mathrm{log}{R}_{{\rm{H}}{\rm{I}}}\\ & & -\,0.99\,\mathrm{log}{R}_{25}-2.69\end{array}.\end{eqnarray} \tag{ 3 }$$

We use Equations (1)–(3) derived above to correct for the three associated observational effects together from photometric scale heights. To test the goodness of such a methodology, we produce mock cubes with the median beam size (15'') and median inclination angle (85°) of the CHANG-ES edge-on sample. We obtain the photometric scale heights, and correct for the observational effects as described above. Figure 4 shows that the corrected scale heights agree very well with the true scale heights. The fit demonstrates that the uncertainty in the photometrically measured scale height is typically 9.36% of the true scale height, when the mock galaxies have median properties of CHANG-ES galaxies.

Zoom In Zoom Out Reset image size

We do not consider the projection effect of warps, for with the limited sensitivity and spectral resolution, as well as the nearly edge-on nature of the sample, it is difficult to reliably identify warps misaligned with respect to the inclination of the main disk. For example, most warps start near or beyond R₂₅ (Briggs 1990) but due to our limited sensitivity we are unable to detect much H i beyond R₂₅ (median R_HI/R₂₅ = 1.14 ± 0.33). Fortunately, this may reduce a putative warp's contribution to our photometric determination of ${\bar{h}}_{\mathrm{phot}}$. Additionally, the potential contamination from warps is difficult to quantify even when kinematical modeling is involved (Yim et al. 2014; Bacchini et al. 2019). For these reasons we leave the effect of warps as a caveat to be investigated in the future.

In summary, these mock tests justify that photometrically derived scale heights can indicate the true scale heights of galaxies after properly accounting for the observational effects of PSF smearing, planar projection and edge-on projection.

3.4. Deriving the Corrected Scale Heights

We present how the H i scale heights averaged within the optical radius, $\bar{h}$, correlate with other galactic properties in the left panels of Figures 5 and 6. As labeled in the figures, the Pearson correlation coefficient (R) is calculated based on the whole sample excluding NGC 5084.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

NGC 5084 is excluded from the correlation coefficient calculation (and fitting for linear relations as well) because of the globally peculiar structure of the whole disk. Most of its H i is distributed on an outer ring-like structure. This H i "ring" is following a faint optical outer ring and is tilted with respect to the main optical disk by about 5° (Gottesman & Hawarden 1986; Zeilinger et al. 1990). The origin of the misaligned H i ring is possibly external gas accretion from satellites, as NGC 5084 is one of the most massive galaxies in its group environment. As the H i ring is beyond the optical disk, its scale height is unlikely to be strongly related with the stellar properties. We keep it in the sample, mainly to demonstrate how far galaxies can deviate from major scaling relations when they are globally unsettled.

For the other galaxies, we use different colors to separate the tidally interacting systems (red) from the relatively unperturbed galaxies (blue).

The right panels of Figures 5 and 6 depict the corresponding previous results presented in Z21 as comparisons. The averaged Gaussian H i scale heights from the method of Z21 are labeled as ${\bar{h}}_{\ {\rm{Z}}21}$ (see Table 2). In Z21, those figures with exponential scale height were presented in Figures 7 and 8 but not discussed in detail, for the systematic broadening effects had not been investigated and quantified. In this paper, we can look into the trends in figures in more detail.

Zoom In Zoom Out Reset image size

Figure 5(a) features a correlation between $\bar{h}$ and R_HI. The black dashed line shows the best-fit linear relation $\mathrm{log}\ \bar{h}=0.49\ \mathrm{log}{R}_{{\rm{H}}{\rm\small{I}}}-0.89$. The most significant outliers include two tidally interacting systems (NGC 3003 and NGC 5775) and the peculiar galaxy N5084. From Figure 5 in Z21, NGC 3003 presents a signature of recent gas-rich minor merger in the southwest corner. From literature studies, NGC 5775 seems to be in the early stage of a major merger and H i masses are likely transferred from its neighbor galaxy NGC 5774. This result shows that H i disks grow thicker when their diameters are larger, which has been found for disks in other wavelengths, for example the optical (de Grijs & Peletier 1997; Zasov et al. 2002) and the radio continuum (Krause et al. 2018). In Z21, we obtained a similar relation between ${\bar{h}}_{\ {\rm{Z}}21}$ and R_HI but with a different slope, for ${\bar{h}}_{\ {\rm{Z}}21}$ values there were not sufficiently corrected for artificial broadening.

In Figure 5(b), $\bar{h}$ is significantly anticorrelated with total mass surface density Σ_Mtot,r25. The globally peculiar galaxy NGC 5084 is again an outlier. Interestingly, most of the tidally perturbed galaxies in our sample do not behave as outliers in the relation. It implies that normal tidal interaction effects are not dominant in determining the thickness of those H i disks. To confirm this speculation, we also test and find no clear correlation between $\bar{h}$ and the local number density of galaxies ρ (taken from Paper I; Irwin et al. 2012). This result holds even when we select galaxies with low mass densities. The data points from Z21 show similar strong anticorrelation between ${\bar{h}}_{\ {\rm{Z}}21}$ and Σ_Mtot,r25, indicating that the trend is strong enough to show itself above the systematic biases.

In Figure 5(c), there is a moderate anticorrelation between $\bar{h}$ and baryonic mass surface density Σ_Mbaryon,r25. The correlation is weaker than that of Σ_Mtot,r25 (implications are discussed in Section 5). Different from the anticorrelation presented here, the trend of Z21 shows no correlation with R value equal to −0.19 with a scatter of 0.28. This inconsistent result implies that the analysis of the systematic biases in this work is useful and a correction is necessary.

In Figure 6(a), for our whole sample, there is no clear correlation between $\bar{h}$ and star formation rate surface density Σ_SFR. We also study the possible relation of $\bar{h}$ with the ratio between Σ_SFR and Σ_Mtot,r25 in Figure 6(b) and no correlation is observed. This is consistent with Z21 in which no correlation was found between ${\bar{h}}_{\ {\rm{Z}}21}$ and Σ_SFR, or between ${\bar{h}}_{\ {\rm{Z}}21}$ and Σ_SFR/Σ_Mtot,r25, when the broadening effects are not as properly corrected for as in this paper.

One of the major results in this paper is that we find the averaged H i scale height $\bar{h}$ to be significantly anticorrelated with the total mass surface density Σ_Mtot,r25. Such a relation seems largely expected from the vertical hydrostatic equilibrium model of H i, as gravity restores the H i to the midplane of the disk. However, in theory, only the mass within the disk layers contributes effectively to the gravity that restores the H i (Forbes et al. 2012; Krumholz et al. 2018), which our mass density parameters do not strictly indicate. For Σ_Mtot,r25, the dark matter certainly extends beyond the disk layers, though there has been evidence that dark matter halos can be oblate in shape and thus concentrate more mass in the disk layers (Navarro et al. 1996, 1997). On the other hand, as typically observed in other massive galaxies (McGaugh et al. 2000; Catinella et al. 2012; Richards et al. 2018), the baryonic-to-total mass ratio is relatively high in our edge-on sample (∼0.35 ± 0.15), leaving space for uncertainties in the dark matter geometry. Σ_Mtot,r25 can be viewed as an upper limit of the gravitational mass restoring the H i disk in the vertical direction. The weaker anticorrelation between $\bar{h}$ and baryonic mass surface density Σ_Mbaryon,r25 implies that Σ_Mbaryon,r25 is not a perfect indicator of the restoring gravity when averaging over the radial range of the whole optical disk. The total mass is not commonly available for extra-galaxies, and thus as a result it is typically assumed that dark matter is not important and the baryonic mass is sufficient to explain most of the effects related to gravity within the optical radius of massive spiral galaxies. However, our results suggest that the effect of dark matter cannot be ignored at least for H i thickness. Typically, stars are dynamically much hotter than the neutral gas (Habets & Heintze 1981), particularly so in massive galaxies (Weidner & Vink 2010; Tan et al. 2014; Krumholz 2015). We may hence overestimate the gravitational contribution from stars, but on the other hand, we have neglected the molecular gas which may have compensated the overestimated stellar contribution. There are also other uncertainties in the relation between $\bar{h}$ and Σ_Mtot,r25 (Σ_Mbaryon,r25). For example, the Poisson equation in cylindrical coordinates indicates that the vertical gravitational acceleration not only depends on the local mass density within a given thickness, but also depends on the changing rate of the disk rotation curve (e.g., Equation (13) of Krumholz et al. 2018). However, our galaxies have high M_* and massive galaxies typically have steeply rising and then flattened rotation curves (de Blok et al. 2008), minimizing the effective gravitational contribution from this term. Despite all the caveats discussed above, $\bar{h}$ anticorrelates with Σ_Mbaryon,r25 and particularly strongly with Σ_Mtot,r25, indicating the CHANG-ES H i disks are strongly regulated by gravity.

Here we only make a preliminary discussion and hypothesis on the results with Σ_SFR, due to our limited sample size. Stronger conclusions should be made through a larger data set in the future. The supernova feedback and streaming motions of gas are both important channels to produce turbulence (Elmegreen et al. 2003; Schnorr Müller et al. 2011; Schmidt et al. 2016; Ramón-Fox & Bonnell 2017) and thus raise the H i disk thickness, and both channels are expected to be associated with star formation. Yet, we do not find a significant relation between $\bar{h}$ and SFR surface densities. A possible reason is that the intrinsic trend is weak and buried under uncertainties, as SFR can be fueled in other ways than gas streaming motions, supernova feedback comes from stars formed earlier than reflected by the current SFR and a large portion of energy from the supernova feedback may be distributed into gas of other phases. The magnetic pressure may also support the disk vertical structures (Krumholz et al. 2018); although magnetic properties have been derived for several galaxies in the CHANG-ES program (Stein et al. 2019, 2020; Krause et al. 2020), directly quantifying the magnetic pressure remains to be done.

Our results of H i scale heights appear to correspond with the study of Krause et al. (2018) (K18 hereafter) about the scale heights in the radio continuum. K18 found the radio scale heights correlate with the halo radius (see Figures 12, 13 and 14 in K18), anticorrelate with mass surface densities (see Figures 15 and 16 in K18), and do not correlate with SFR surface densities (see Figure 19 in K18) both in C- and L-band. These consistent results reflect the similar way of organizing their vertical structures for the different constituents of the interstellar medium (ISM). In addition, we provide strong support for arguments in K18 that the underlying gravitational potential plays a more important role than star formation in determining disk scale heights. Star formation may be important in determining scale heights and structures in individual locations, but globally, the potential appears to dominate.

The really strong correlation between $\bar{h}$ and R_HI raises the possibility that other dependencies of $\bar{h}$ shown in this study are artificially caused by a mutual dependence on R_HI. Following the study in Z21, we thus test by showing the correlation between the normalized H i scale heights $\bar{h}/{R}_{{\rm{H}}{\rm\small{I}}}$ and other galactic properties in Figure 7. The $\bar{h}/{R}_{{\rm{H}}{\rm\small{I}}}$ shows no correlation with Σ_Mtot,r25, Σ_Mbaryon,r25, Σ_SFR or Σ_SFR/Σ_Mtot,r25. The same results were presented in Z21 and the right column of Figure 8, where the systematic biases are not corrected. These results imply that it is indeed possible that the relation between $\bar{h}$ and Σ_Mtot,r25 is caused by the more intrinsic correlation between $\bar{h}$ and R_HI and are not independent or new. On the other hand, this dependence may have been blurred by other dependencies (e.g., $\bar{h}$ versus Σ_SFR) which have contrary slopes. More data will be needed to test both hypotheses in the future.

In the end, we reiterate that the biggest caveat of our analysis is that, due to the lack of kinematic information, we only obtain a rough measure of H i scale heights averaged over a large radial range, and miss the typical flaring feature of galactic gas disks. The systematic uncertainty due to averaging needs to be quantified in the future, based on edge-on H i disks with properly modeled 3D distributions.

We measure H i scale height from 15 highly inclined (inclination > 80°) galaxies from the CHANG-ES sample. We build mock data cubes and images for the H i, to explore the feasibility of deriving disk thickness in the photometric way. We show that when the inclinations of disks are not perfectly 90°, the directly measured disk scale heights can be systematically overestimated as a result of planar projection, beam smearing and disk flaring. We quantify the trends of these overestimating effects, which are used as correcting equations for the directly measured disk scale heights. We derive the radially averaged Gaussian scale height within the optical radius for each galaxy and apply these corrections. We find a significant anticorrelation of the corrected scale height with the total mass surface density and the baryonic mass surface density. Compared with the results in Z21, these results imply that either the systematic broadening effects are corrected well in our method, or the underlying physical driver is strong enough to show itself above the biases.

The result is consistent with predictions from the quasi-equilibrium model (Poisson equation) of H i vertical distribution, where gravity is the primary force restoring the H i to the midplane of the disk (Krumholz et al. 2018).

We thank L. Bai, K. Zheng, X. Feng and Y. Fu, and especially thank J. Zhu, Y. Yang, N. Yu, Z. Liang and H. Xu. J. Wang acknowledges support from the National Natural Science Foundation of China (NSFC, Grant Nos. 12073002 and 11721303). Parts of this research were supported by the High-performance Computing Platform of Peking University. This research made use of Photutils, an Astropy package for detection and photometry of astronomical sources (Bradley et al. 2020), as well as NOD3 (Müller et al. 2017), Photutils (Bradley et al. 2020) and BBarolo (Di Teodoro & Fraternali 2015).