© ESO 2020

1. Introduction

What is the structure of orbits in the interstellar medium of disk galaxies? What is the shape of the gaseous velocity dispersion ellipsoid in the mid-plane of galactic disks?

A quick look at observations of nearby galaxies is helpful to guess that orbits should not be perfectly circular, owing to the presence of large-scale perturbations like bars, spiral arms, warps, or disk lopsidedness. Rotational motions are accompanied by radial and streaming motions, and orbits become asymmetric (Kalnajs 1973; Visser 1980). In a galactic disk, the structure of orbits can be studied by means of the shape of the velocity dispersion ellipsoid σ_R, σ_ϕ, and σ_z, which are respectively the radial, tangential, and vertical components of random motions in cylindrical coordinates. The ellipsoid is characterized by three axis ratios, one of which is the azimuthal-to-radial ratio that traces the degree of anisotropy in the mid-plane, as characterized by the azimuthal anisotropy parameter, β_ϕ = 1 − (σ_ϕ/σ_R)². Orbits that are biased tangentially have β_ϕ <  0, with β_ϕ → −∞ for circular orbits, whereas those biased radially have 0 <  β_ϕ <  1 (Binney & Tremaine 2008). This parameter is appropriate to explain the orbital structure of a kinematic tracer that is collisionless (stars). We would expect to observe such values within gaseous disks of galaxies, if interstellar gas behaved partly like a collisionless medium, as shown in numerical models (Bottema 2003; Agertz et al. 2009).

Interstellar gas is a medium in which collisions and shocks coming from all directions are thought to completely “isotropize” the motions. Assuming a velocity ellipsoid that is not tilted, the line-of-sight dispersion can be written

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{los}}=(({\sigma }_{\mathrm{R}}^2{sin}^2\varphi +{\sigma }_{\varphi }^2{cos}^2\varphi ){sin}^{2}i+{({\sigma }_{z}cosi)}^2+{\sigma }_{\mathrm{T}}^2+{\sigma }_{\mathrm{ins}}^2{)}^{1/2},\end{array}$$(1)

where i is the disk inclination, ϕ the azimuthal angle in the disk plane, σ_T the thermal (isotropic) component, and σ_ins the instrumental broadening. Isotropy greatly simplifies the modeling of the observed dispersion as it implies σ_R = σ_ϕ = σ_z = σ_iso, reducing Eq. (1) to

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{los}}={({\sigma }_{\mathrm{iso}}^2+{\sigma }_{\mathrm{T}}^2+{\sigma }_{\mathrm{ins}}^2)}^{1/2}.\end{array}$$(2)

Important values for galactic dynamics like the gas radial pressure support depending on σ_R (Dalcanton & Stilp 2010; Oh et al. 2015) or the vertical equilibrium depending on σ_z (Combes & Becquaert 1997; Koyama & Ostriker 2009; Bershady et al. 2010; Martinsson et al. 2013) can then be easily estimated directly from observations, without needing to deproject σ_los.

The analysis of σ_los is less straightforward in the hypothesis that part of the gas behaves as a collisionless medium. Equation (1) has to be modeled or inverted in order to fit or deproject (respectively) the data, and constrain the shape of the velocity ellipsoid. Moreover, Eq. (1) predicts that σ_los is a function of the azimuthal angle owing to the projection of the radial and tangential components into σ_los. In an axisymmetric approach of this kind, an anisotropic dispersion ellipsoid has a significant implication on σ_los: it has to be asymmetric, greater near the minor (major) axis for σ_ϕ <  σ_R (σ_ϕ >  σ_R, respectively). To our knowledge, such strong prediction has never been tested. Many observed velocity dispersion fields exhibit asymmetric σ_los dependent on the azimuth (Fig. 1, and Walter et al. 2008), and some of them are reminiscent of this particular anisotropic signature.

Fig. 1. Observational data of Messier 33. Left panel: HI column density map (VLA, 100 pc resolution, logarithmic stretch) with Spitzer/IRAC 3.6 μm stellar distribution overlaid (gray contours, showing the 0.1, 0.3, 0.5, 0.7, 0.9, 1.4, and 2.5 MJy sr−1 levels). Right panel: observed HI velocity dispersion (VLA, 100 pc resolution, not corrected for the instrumental dispersion), with HI column densities contours (1, 2, 3 × 1021 cm−2). The dashed line represents the location of the major axis of the inner unwarped disk, with position angle of 22.5°. Concentric ellipses show the projected locations of R = 2, 4, 6 kpc.

The two competing collisional and collisionless models for the interstellar medium thus predict different shapes for the velocity dispersion ellipsoid, respectively isotropic and anisotropic, that are important to study from an observational viewpoint. At the same time, the traditional assumption of axisymmetric motions under the isotropy argument seems not suitable to asymmetric observations as it is independent of azimuth, by construction. It is therefore important to evaluate the degree of asymmetry in velocity dispersion fields as well.

We want to address these problems in this work. Our long-term objective is to determine observationally whether collisional or collisionless models are appropriate to explain the structure of gaseous velocity dispersion fields. We want to study the asymmetries in velocity dispersion maps for the atomic and/or molecular gas in nearby galaxies, and estimate the azimuthal velocity anisotropy parameter that is needed to reproduce asymmetric observations. We start this study with the Local Group spiral Messier 33 (M 33). Its proximity (D = 840 kpc) has made it possible to obtain very high-quality data in the HI and CO lines, tracing the atomic and molecular gas components (Sect. 2). Velocity dispersion maps are used to quantify asymmetries assuming an isotropic velocity model (Sect. 3), and estimate β_ϕ and the corresponding structure of gas orbits for the atomic and molecular gas in M 33 assuming an anistropic velocity model (Sect. 4). A short discussion on the two competing models and comparisons with numerical models and the velocity anisotropy expected in the framework of the epicycle theory of collisionless orbits is also presented (Sect. 5).

2. Atomic and molecular gas data of M 33

The HI data of M 33 come from three studies. Our main H I reference is the 25″ resolution (100 pc) Very Large Array (VLA) data from Gratier et al. (2010). We also use other recent VLA observations at 18″ resolution (70 pc) by Koch et al. (2018). Both data sets are very low noise, which is essential to calculate reliable second-moment maps (velocity width). Even though from the same telescope (the VLA), the data are from completely independent observations. We also use the low-resolution Dominion Radio Astrophysical Observatory (DRAO) data set by Kam et al. (2017, 2′, 490 pc) as a further check. Both are useful to estimate any link between resolution and the anisotropies presented here. The pixel size of the VLA velocity dispersion maps used in this study is 8″, and that of the DRAO map is 21.8″.

Velocity widths via calculation of the second moment of a spectrum are extremely subject to noise, particularly the noise in channels far from the line center. The velocity dispersion (second moment of a spectrum) is σ_los, as given by ${\sigma }_{\mathrm{los}}^2=\frac{{\int }_{{\mathit{v}}_{\min }}^{{\mathit{v}}_{\max }}T{(\mathit{v}-{\mathit{v}}_{\mathrm{cen}})}^2\mathrm{d}\mathit{v}}{{\int }_{{\mathit{v}}_{\min }}^{{\mathit{v}}_{\max }}T\mathrm{d}\mathit{v}}$, where v_cen is the central (intensity-weighted mean) velocity. The (v − v_cen)² term makes it critical to exclude noise beyond the line profile. Thus, we choose the low-noise HI data at 100 pc resolution presented in Gratier et al. (2010) to define the windows over which to calculate the moments. Over 90% of the surface of M 33 reaches a line temperature of 10 K in HI, corresponding to approximately four times the noise level in the data. We do not calculate the velocity dispersion where the line temperature does not reach 10 K because the uncertainty becomes too great, both on v_cen and thus on the dispersion. The 10 K cut applies only for the 100 pc resolution data.

The integration window is defined by taking the maximum and then descending to either side (higher and lower velocities) until a channel goes below zero. In less than 1% of the disk, a double-peaked profile can be observed, but at 100 pc resolution there are no pixels for which the line temperature goes below zero between the two peaks. A mask is created for each position, containing a value for v_min and v_max, and a flag indicating whether the line temperature reaches 10 K. This is the mask used in all the second-moment calculations of the HI gas. We note that the velocity dispersions for the 70 pc and 490 pc resolution data are ∼1 km s⁻¹ larger than those measured with the 100 pc resolution VLA data. As single-dish data have been merged to the 70 pc VLA and 490 pc DRAO interferometric data, but not to the 100 pc VLA data, combined data are more sensitive to larger scale structures that slightly widen the HI profile wings. This small difference has no impact on the results. Figure 1 shows examples of column density and velocity dispersion fields for the atomic gas (100 pc resolution VLA data).

The HI disk of M 33 is known to locally exhibit emission with anomalous velocities (Kam et al. 2017; Koch et al. 2018). This emission, which has typical characteristics of an extraplanar HI layer lagging the rotation of the disk (e.g., in the galaxy NGC 2403; Fraternali et al. 2001, 2002) with forbidden velocity gas and high-velocity components (see Sect. 3.2. of Kam et al. 2017), could participate in broadening the HI profiles. Quantifying the effect of such gas on the HI velocity dispersion models has not been attempted, however.

The CO observations are at 12″ resolution (50 pc) and were made at the 30 m telescope of the Institut de Radioastronomie Millimétrique (IRAM), as described in detail in Druard et al. (2014). The velocity dispersions for the molecular gas are not strictly derived the same way as for HI. The molecular gas of M 33 is organized into clouds. Although there is emission which appears diffuse, we do not know whether it is truly diffuse molecular gas or small clouds producing a weak but extended CO signal. Either way, the signal-to-noise ratio for the “diffuse” component is too low to try to calculate a velocity dispersion map from the IRAM datacube. Therefore, we use the positions and velocity dispersions of the 566 clouds presented in Corbelli et al. (2017) and Braine et al. (2018) rather than a continuous map of the velocity dispersion. The decomposition into clouds presented in these articles was made using CPROPS (Rosolowsky & Leroy 2006).

The instrumental dispersions of the HI and CO observations is σ_ins ≲ 1 km s⁻¹. In the kinematic modeling, the isotropic thermal component σ_T was chosen to be constant with radius. For the atomic gas, fits were performed assuming two extreme cases: HI seen as a cold or warm neutral medium (CNM and WNM, respectively with T ∼ 100 K and T ∼ 5000 K, Draine 2011). The CNM has a thermal dispersion of ∼0.9 km s⁻¹, which is comparable to the instrumental dispersion and significantly smaller than the observed HIσ_los. For simplicity, we thus set σ_T = 0 km s⁻¹ for the CNM case. The WNM has a thermal dispersion of 6.4 km s⁻¹, thus larger than the instrumental dispersion, and comparable to σ_los in many disk regions. For simplicity we thus set σ_T = 6 km s⁻¹ for the WNM case. Given that the majority of fits at σ_T = 6 km s⁻¹ failed (see below) most of results shown hereafter are those obtained for the CNM case, unless specified. The fits to the molecular gas dispersion were made assuming only a cold medium (σ_T = 0 km s⁻¹).

This work focuses only on the radial range R ≤ 7.5 kpc that is not affected by the warping of the gaseous disk (Corbelli & Schneider 1997; Corbelli et al. 2014; Kam et al. 2017). The projection of the dispersion model of Eq. (1) along the line of sight thus assumed a fixed inclination of 56°, and a fixed orientation of the major axis shown as a dashed line in Fig. 1, with a position angle of 22.5° (ϕ = 0 is chosen aligned with the semi-major axis of the approaching half to the northeast, and increases in the counterclockwise direction). The fits were performed at radii starting from R = 20″, 25″, and 120″ for the 70 pc, 100 pc, and 490 pc resolution data, respectively, with a radial bin width of 20″, 25″ and 120″. The fits were performed at radii starting from R = 60″ for the molecular gas, with a radial bin width of 120″. Fits were performed within the radial range R ≲ 1500″ for the molecular gas, using 553 from the initial 566 clouds. We also note that the sampling is different than for the VLA data, although the CO data were obtained at higher angular resolution. The smaller sampling for the discrete molecular gas clouds is necessary to have a number of degrees of freedom large enough to yield successful least-squares fits. For both components we thus excluded rings that had fewer than 10 points.

The telescope beams cover a larger area near the minor axis than on the major axis so the beam-smearing could induce a slightly higher velocity dispersion along the minor axis, and thus an apparent (but false) anisotropy (Chemin 2018). The estimation of the beam-smearing effect was performed using a high-resolution model of gas intensity and velocity of M 33 smoothed to the native resolution of observations (100 and 490 pc), following prescriptions given in Epinat et al. (2010). In particular, we used the rotation curve model given in Koch et al. (2018), and assumed a constant surface density (Kam et al. 2017) and constant velocity dispersion (σ_los = 8 km s⁻¹) within R ≤ 7.5 kpc. We found that the contamination of the VLA dispersion map by the smearing effect is completely negligible along the principal axes (≪1 km s⁻¹). The effect is stronger at a resolution of 490 pc, but the agreement between the resolutions shows that the beam-smearing is not an issue. For simplicity we assumed that artificial changes in dispersion with azimuth caused by other mechanisms, like possible finite disk thickness (Bacchini et al. 2019), are negligible.

Finally, we add that the models presented in Sects. 3 and 4 can be applied to linear or squared velocity dispersions. In the squared velocity case, Eqs. (1) and (2) resume to simple additions which could make the analysis easier at first glance. We thus performed the modeling of both σ_los and ${\sigma }_{\text{los}}^2$ and found no significant difference between results from the two approaches. Hereafter, the presented results are those obtained for the modeling of σ_los.

3. Evidence of asymmetric line-of-sight H I velocity dispersions in M 33

Asymmetries are obvious in the HI velocity dispersion map of M 33 (Figs. 1 and 2). The variations in σ_los with the azimuthal angle occur at all angular scales. There is for example a clear hint of a large angular scale variation for R >  4 kpc, with σ_los lower towards the major axis of the disk and larger towards the minor axis. In that region the difference between individual pixel values σ_los and the azimuthally averaged dispersion ⟨σ_los⟩_ϕ at the same location is easily > 25% of σ_los, and can sometimes reach up to 60%. A mean axisymmetric and isotropic dispersion is thus clearly ruled out to explain the observed variations occurring on larger angular scales (see also Sect. 4). The anisotropic and isotropic models presented below attempt to capture such large angular scale variations, not those occurring on angular scales smaller than a few degrees (like the pixel-to-pixel variations).

Fig. 2. Azimuth-velocity dispersion diagrams of atomic gas at selected radii in M 33. Each point corresponds to an individual measurement, i.e., the dispersion of each pixel within the map. Solid lines represent the results of a Fourier Transform model of the HI dispersion map (Sect. 3, assuming a cold neutral medium case). The dispersions are from the VLA 100 pc resolution data.

To investigate the properties of the asymmetric random motions of the atomic gas, we expanded the isotropic random motions σ_iso of Eq. (2) with Fourier coefficients, following σ_iso, asym = ∑_kσ_kcosk(ϕ − ϕ_k), where k is an integer. The isotropic axisymmetric component is σ₀ and σ_k and ϕ_k are the amplitude and phase of the asymmetric modes in the velocity dispersion map. Both least-squares (LSQ) fits and discrete fast Fourier transforms (FFTs) were performed at the radial rings defined in Sect. 2 to cross-validate the analysis. The minor differences between fits and FFTs are that fits provide formal errors for the phases and amplitudes, and are more time consuming than FFTs. To avoid possible degeneracies, divergence, and increased computational time during minimizations, the maximum order that was fit is k = 4 (harmonics of k = 2). This maximum order is much smaller than the last order yielded by FFTs, which is the number of points in the azimuthal dimension at a given radius. The instrumental and thermal contributions were subtracted from the original maps before derivations of the kinematic harmonics.

Figure 3 shows the results of the FFT derivation for the 100 pc resolution VLA observation in the CNM case, with the amplitudes (top panel) and phases (middle panel) of the first four Fourier modes. The comparison of the bisymmetric amplitude with the first-, third-, and fourth-order terms is shown in the bottom panel of Fig. 3. Results of the LSQ fits are not shown because they are extremely similar to the FFTs.

Fig. 3. Results of the discrete FFT of the HI velocity dispersion map of M 33 (100 pc resolution, in the CNM case). The amplitudes (upper panel) and phases (middle panel) of the Fourier modes are shown as black open circles (k = 0), green open circles (k = 1), filled circles (k = 2), blue upward filled triangles (k = 3), and red open diamonds (k = 4). The reference ϕ = 0 is chosen along the semi-major axis of the NE approaching half of M 33. Bottom panel: σ2/σk, using the same colors and symbols as above for the k = 1, 3, 4 modes. The range of ratios is chosen up to 10 for clarity, but the total range is 78 for σ2/σ1 and σ2/σ4, and 12 for σ2/σ3.

The amplitude of the axisymmetric term is unsurprisingly similar to the azimuthally averaged velocity dispersion. The other modes are weak, with median amplitudes of σ₁ = 0.6 km s⁻¹, σ₂ = 1.2 km s⁻¹, σ₃ = 0.6 km s⁻¹, and σ₄ = 0.7 km s⁻¹. The bisymmetry is the most interesting mode. It shows the largest variation, which can reach 25% of the axisymmetric value at around R = 1500″ (6 kpc). It is stronger than σ₁, σ₃, and σ₄ in ∼80%, 80%, and 65% of the cases, respectively. In the inner region (R <  1000″, or 4 kpc), these ratios fall to 65% for the k = 1 and k = 3 amplitudes, and to 40% for k = 4. We can thus define the locations where the bisymmetry is not stronger than other asymmetries as the locations where σ₂ is weaker than at least two other modes, that is 188″ ≤ R ≤ 237″, 513″ ≤ R ≤ 687″, and 737″ ≤ R ≤ 837″. These regions are highlighted in subsequent figures. In the outer region (R >  1000″) it significantly dominates the other modes, by factors of 5, 4, and 7 on average (respectively to the k = 1, 3, 4 modes), reaching up to 78 times the k = 1 and k = 4 terms, and 12 times the k = 3 term.

The phases of the asymmetries contain a wealth of information as well. They show extended angular ranges where they remain closer or aligned with one of the principal axes of the M 33 disk, modulo half or one period of each mode. In other words, ϕ₁ is consistent with the position of the minor and major axes (2.4 <  R <  3.4 kpc, 4 <  R <  5 kpc, R >  6 kpc) with some variation elsewhere, ϕ₂ is mostly aligned with the minor axis, except for R = 0.7 and 1.5−2.5 kpc, ϕ₃ remains close to 2π/3 in the inner disk half and to π/3 beyond R = 5.2 kpc (the period and half-period of the k = 3 mode, respectively), while ϕ₄ shows three steps, one close to the major axis for R <  3 kpc (ϕ₄ ∼ 0), one with the half-period of the k = 4 mode for 3 <  R <  5.5 kpc (ϕ₄ ∼ π/4) and another close to π/3 beyond 5.5 kpc. In the azimuth-dispersion diagrams drawn at selected radial rings of 0.2 kpc in width (Fig. 2), the solid lines show the results of the FFT models, highlighting the effects of the bisymmetry (all panels), the k = 1 mode (R = 3 kpc, middle panel), and the k = 4 mode (R = 1.5 kpc, left panel).

Finally, comparable results have been obtained with the 70 pc and 490 pc data. We also performed the same analysis in the WNM case for the 100 pc resolution data. The phases are comparable to the CNM case, but on average σ₀ drops by 25%, while the median σ₁, σ₂, σ₃, and σ₄ rises respectively by a factor of 1.2, 1.1, 1.4, and 1.6, though it remains small. The k = 2 mode thus still dominates the asymmetric random motions in the WNM case.

4. Velocity anisotropy of gas in M 33

We now consider that gas can behave like a collisionless medium. The axisymmetric and anisotropic velocity model of Eq. (1) can be recast in

$$\begin{array}{cc}\hfill {\sigma }_{\mathrm{los}}=& (0.5({\sigma }_{\varphi }^2-{\sigma }_{\mathrm{R}}^2){sin}^{2}icos2\varphi +0.5({\sigma }_{\varphi }^2+{\sigma }_{\mathrm{R}}^2){sin}^{2}i\hfill \\ \hfill +& {({\sigma }_{z}cosi)}^2+{\sigma }_{\mathrm{T}}^2+{\sigma }_{\mathrm{ins}}^2{)}^{1/2}.\hfill \end{array}$$(3)

The expression with cos2ϕ therefore implies that such an anisotropic velocity dispersion model looks like a bisymmetric mode that has a null phase, i.e., aligned with the major or minor axis (depending on the sign of ${\sigma }_{\varphi }^2-{\sigma }_{\text{R}}^2$). The significant bisymmetric perturbation mostly aligned with the minor axis found in the HI velocity dispersion map of M 33 could thus be interpreted as a signature of velocity anisotropy, except maybe for the clearly identified regions of weaker k = 2 mode (see Sect. 3). The anisotropic velocity model makes it possible to constrain the two planar components and study the structure of orbits of gas in M 33 through the azimuthal anisotropy parameter β_ϕ. We also apply this model to observations of the molecular gas (Sect. 4.2).

Non-linear LSQ fits of Eq. (1) to σ_los were performed using the program developed in Chemin (2018) for stellar disks. It fits radial, tangential, and vertical components of the random motion ellipsoid to a vector of observed dispersions at a given radius. A more developed model where the anisotropic dispersion components are asymmetric should also be considered to allow direct comparisons with Sect. 3. Yet such modeling is beyond the scope of the article. We refer to Chemin (2018) for more details of the minimization process and the validation of the methodology.

As σ_R, σ_ϕ, and σ_z are expected to be roughly comparable, degeneracies can occur if the model is fit blindly to the observations. To make the analysis possible, we chose to hold the vertical component constant and let the radial and tangential dispersions vary freely. The value of σ_z can then be varied to measure its effect on the values of σ_R and σ_ϕ. Fits with σ_z assumed constant as a function of R were performed, as well as others where σ_z varied with radius. In this latter case, it was fixed at the azimuthally averaged dispersion ⟨σ_los⟩_ϕ (dashed lines in Figs. 4, 5, and 9). In the case when σ_z was assumed constant, iterations in the range 4 ≤ σ_z ≤ 10 km s⁻¹ were made for the atomic gas, sampled every 0.2 km s⁻¹. This broad range of values contains the line-of-sight dispersion of ∼8 km s⁻¹ observed in nearly face-on HI disks (Shostak & van der Kruit 1984; van der Kruit & Shostak 1984), and which is a good proxy for σ_z. For the molecular gas, we explored 2 ≤ σ_z ≤ 5 km s⁻¹ because the velocity dispersions are lower. The resulting anisotropy parameters are not observed to be strongly dependent on the choice of σ_z (see also Chemin 2018). We define the quoted error on σ_R and σ_ϕ at a given value of σ_z and at a given radius as the standard deviation of the posterior distribution of the parameters, and that for β_ϕ as the propagation of the σ_R and σ_ϕ uncertainties.

Fig. 4. Results of the anisotropic and axisymmetric velocity model for the atomic gas in Messier 33 (100 pc resolution data). Results are those of the one-parameter model and show the profile of σR (σϕ, respectively) as filled green circles (open diamonds), obtained assuming σϕ = σz = ⟨σlos⟩ϕ (σR = σz = ⟨σlos⟩ϕ). A dashed blue line is for ⟨σlos⟩ϕ (corrected from instrumental dispersion). Results obtained assuming a null thermal component. The shaded areas highlight the regions where the k = 2 mode was found to be weaker than other dispersion asymmetries (see Sect. 3).

Fig. 5. Same as Fig. 4, but for the two-parameter model with σR and σϕ as free parameters.

4.1. Results for the atomic neutral hydrogen

4.1.1. One-parameter model

A first goal is to assess whether isotropy is present within the model of Eq. (1). In this context, we tried to maximize the likelihood of finding isotropy by fixing σ_R = σ_z = ⟨σ_los⟩_ϕ, with σ_ϕ as the only free parameter, and alternatively σ_ϕ = σ_z = ⟨σ_los⟩_ϕ with σ_R as free parameter. We would indeed expect σ_ϕ ≃ ⟨σ_los⟩_ϕ and σ_R ≃ ⟨σ_los⟩_ϕ, respectively, in those particular cases. The result is displayed in Fig. 4 (for the CNM case). The locations where σ_R or σ_ϕ is closer to ⟨σ_los⟩_ϕ are R ∼ 1 kpc, and 2 ≤ R ≲ 3.5 kpc. These radii match closely those where the bisymmetry was found dominated by other asymmetries (Sect. 3). The radial (azimuthal) component is observed to be larger (smaller, respectively) than ⟨σ_los⟩_ϕ in most disk regions. Beyond 4 kpc, the difference with ⟨σ_los⟩_ϕ is important, up to 2 and 3 km s⁻¹ for σ_R and σ_ϕ (respectively), which is significant compared to the formal errors of the fits (≪1 km s⁻¹). Comparisons between the two isotropic fits show that the model at free σ_R is significantly worse than the one at free σ_ϕ, particularly in the outer disk. Indeed the assumption σ_ϕ = ⟨σ_los⟩_ϕ automatically maximizes the tangential component, and in regions where the velocity dispersion map exhibits larger values near the minor axis, the model has no other choice than finding an even larger radial component than ⟨σ_los⟩_ϕ, hence yielding a best-fit σ_los model that overestimates the bulk of the observed σ_los. We also performed fits with σ_T = 6 km s⁻¹ (not shown), but they failed in more than 75% of the rings. Results for a one-parameter case that most closely approaches an axisymmetric and isotropic model thus confirm the results found in Sect. 3: axisymmetry with isotropy cannot apply to the HI random velocity field of M 33.

4.1.2. Two-parameter model

Figure 5 shows the results found for σ_R and σ_ϕ for the VLA 100 pc resolution map, assuming a null thermal component. The results corresponding to the warm neutral medium case are briefly discussed below. Results for σ_z = ⟨σ_los⟩_ϕ only are shown for clarity. Beyond R = 2 kpc, the azimuthal dispersion decreases by ∼6 km s⁻¹, while the radial dispersion tends to increase slightly, although bumps are observed. In the inner kiloparsec, the overall variation of σ_R reaches ∼5 km s⁻¹. For the most realistic σ_z values (σ_z <  10 km s⁻¹) the results yield 0.5 ≤ σ_z/σ_R ≤ 0.9.

The resulting anisotropy parameter in the disk mid-plane is shown in Fig. 6 for the three resolutions of the HI data, again for σ_z = ⟨σ_los⟩_ϕ (solid lines). The result for another value, σ_z = 4 km s⁻¹, for the 70 pc resolution data is also given (dotted line). It shows that the impact on β_ϕ of the choice of σ_z is negligible as differences in anisotropy parameter ≲0.15 are observed. A similar finding with stellar velocity anisotropy was presented in Chemin (2018).

The anisotropy parameters of the VLA 70 and 100 pc resolution H I data match perfectly. In the inner R = 3.5 kpc, β_ϕ is highly variable, sometimes corresponding to isotropic-to-radial motions (≲0.25), sometimes more radially biased (∼0.4), or showing dips down to −0.4 (R ∼ 0.7 and 2.2 kpc). Beyond R = 3.5 kpc, β_ϕ increases steadily reaching ∼0.8, showing that orbits become more radial at large radii.

The radially oriented orbits at large radius are also observed in the 490 pc resolution DRAO data (bottom panel of Fig. 6). The variations occurring on small angular scales are lost, however, because of the lower resolution. The dip of β_ϕ at R ∼ 0.7 kpc is not observed, while the second dip seems to be detected at ∼2.5 kpc, although at a lower (absolute) amplitude than at higher resolution, making the velocity ellipsoid more consistent with isotropy. Given the coincidence of the dips within or close to the regions where σ₂ was measured lower, we interpret such negative values as a failure of the anisotropic velocity model at these locations.

Fig. 6. Profiles of azimuthal velocity anisotropy βϕ of the atomic gas in M 33 at different angular resolutions. Solid lines show the anisotropy profiles obtained assuming σz = ⟨σlos⟩ϕ and a null thermal component. Short-dashed lines show the ±1 rms errors; a dotted line is an illustration of result choosing another vertical dispersion, σz = 4 km s−1 (70 pc resolution only); and a blue long-dashed line is the velocity anisotropy profile expected from the epicycle theory (Sect. 5). Shaded areas are as in Fig. 4.

The effects of isotropic and radial orbits on σ_los are shown in azimuth-dispersion diagrams of Fig. 7 (solid lines), as extracted from the 70 pc and 490 pc resolution HI observations. The widths of the radial rings were set to 0.1−0.2 kpc for the 70 pc resolution and 0.3−0.5 kpc for the 490 pc resolution. The left column illustrates locations where the anisotropic models found isotropy, the middle column is for gas orbits that are found more radial, and the right column illustrates one of the locations of a hypothetical tangentially biased orbit in the VLA data. While the model dispersion clearly varies with ϕ on a large angular scale, it does not account for the variations seen on smaller scales.

Fig. 7. Same as Fig. 2, but using HI data from VLA (70 pc resolution), DRAO (490 pc resolution), and CO data from IRAM 30 m antenna (50 pc resolution) to show results of the anisotropic velocity dispersion model (Sect. 4). The columns illustrate various cases of velocity anisotropy parameter found by the best-fit model: βϕ ∼ 0 (for isotropy, left ), 0 <  βϕ (for radial bias, middle), βϕ <  0 (for supposedly tangential bias, right). For each resolution, the selected radii are among the best positions illustrating each case of anisotropy parameter, following Fig. 6. Solid lines represent the best-fit dispersion model at the considered radii, assuming σz = ⟨σlos⟩ϕ.

As for the case of HI seen as a warm neutral medium (σ_T = 6 km s⁻¹), β_ϕ is shown in Fig. 8, again assuming σ_z = ⟨σ_los⟩_ϕ. The anisotropy could not be derived for the outer disk because there is little room left for both planar components. This shows that at least some of the HI is cool in the outer region within the anisotropy assumption. The remaining velocity anisotropy is stronger than in the case σ_T = 0 km s⁻¹, corresponding to even more radially biased gas orbits.

Fig. 8. Profile of azimuthal anisotropy βϕ of the velocity dispersion of the atomic gas in M 33 obtained using the 100 pc resolution VLA data by assuming σz = ⟨σlos⟩ϕ and HI as a warm neutral medium (thermal component of 6 km s−1). Shaded areas are as in Fig. 4.

4.2. Results for the molecular gas

Results of the two-parameter model for the molecular gas, derived from the CO data, are shown in Fig. 9, obtained assuming σ_z = ⟨σ_los⟩_ϕ. Examples of azimuth-dispersion diagrams with the σ_los models are shown in Fig. 7 (ring of 0.35 kpc in width). Similarly to the atomic gas, the variation of the molecular gas velocity anisotropy as a function of σ_z is negligible within the spanned range of σ_z. The velocity anisotropy parameter is highly scattered (⟨β_ϕ⟩∼0, on average, throughout the disk).

Fig. 9. Results for the molecular gas in Messier 33. Top: profiles of σR and σϕ of CO gas, obtained assuming σz = ⟨σlos⟩ϕ (dashed blue line, corrected for instrumental dispersion) and a null thermal component. Filled circles are for σR and open diamonds are for σϕ. Bottom: profiles of azimuthal velocity anisotropy βϕ of CO gas. The solid line shows the anisotropy obtained assuming σz = ⟨σlos⟩ϕ and a null thermal component, and dashed lines show the ±1 rms errors. For comparison, the anisotropy of the atomic gas as derived from all pixels in the 100 pc resolution velocity dispersion map is shown as a dotted line, and that derived using only the pixels at the positions of the molecular clouds as open triangles. Filled circles are values derived by masking a few deviant dispersions of CO clouds (σlos >  5 km s−1). Shaded areas are as in Fig. 4.

We inspected the origin of the significant β_ϕ troughs at R = 1.7, 4.6, and 5.1 kpc and found they are caused by the presence of a few deviant observed dispersions, namely σ_los ∼ 6.2 km s⁻¹ along ϕ = 0 (R = 1.7 kpc), σ_los ∼ 7.2 km s⁻¹ along ϕ = 0 (R = 4.6 kpc), or σ_los ∼ 6.4 km s⁻¹ along ϕ = π (R = 5.1 kpc), while at these radii the other points are almost exclusively below 4.5 km s⁻¹. The masking of such outlying points made the anisotropy parameter closer to 0 (filled circles in Fig. 9). A close inspection of σ_los at R = 5.1 kpc shows no clear sine pattern that could yield β_ϕ <  0. Therefore, as for the HI, there is no evidence of a strong tangential bias of velocity anisotropy in the molecular gas in M 33.

The CO velocity anisotropy parameter is roughly in agreement with that of the atomic gas (dotted line) in the inner region, with the exception of R = 1.7 kpc. Interestingly, at larger radius where the HI velocity dispersion has become strongly anisotropic, the orbital structure of CO and HI in the mid-plane differ fundamentally. To verify whether these differences may be artifacts, we measured the velocity anisotropy of the atomic gas in a comparable way to that of the CO gas. This was done by using only the HI velocity dispersions at the locations of the discrete molecular clouds, instead of the whole HI dispersion map. The result shown as open triangles in Fig. 9 indicates that this HI anisotropy profile perfectly agrees with the profile inferred using the whole velocity dispersion field, including the outer regions with stronger HI radial bias. This suggests that the sparser distribution of points in the disk for the molecular gas than for the atomic gas is not the cause of the molecular-atomic difference.

In summary, we find no compelling evidence for a velocity ellipsoid of the molecular clouds being aligned systematically towards (or perpendicular to) the direction of the galactic centre of M 33. This result may indicate that the dynamics of clouds is locally dominated by the cloud gravitational potential. It also highlights the need for velocity dispersion maps of molecular gas in galaxies rather than cloud-based measurements to make the comparison with the HI gas more appropriate.

5. Discussion

This work is the first to our knowledge that examines the effects of the collisionless medium hypothesis for gas on the structure of velocity dispersions, and the implied azimuthal-to-radial axis ratio of the velocity ellipsoid in galactic disks. Therefore, no fully appropriate comparisons with other observational studies of the gas component are available.

Comparisons can be made with the stellar collisionless kinematic tracer, however. Observations of radially biased stellar random motions is not rare among nearby spiral galaxies. Taking the example of our Galaxy, the velocity anisotropy of stars in the disk of the Milky Way derived by Chemin (2018) using stellar dispersions from the Second Gaia Data Release published in Gaia Collaboration (2018) shows quite similar values to the atomic gas of M 33 inside R = 4 kpc, as well as the increase towards large radii. A possible mechanism for the origin of stellar anisotropic orbits in the Galaxy may be radial migration induced by the dynamics of spiral arms (Roškar et al. 2012; Grand et al. 2014). The similarity of the velocity anisotropy of HI gas in M 33 to that of the stars in the Milky Way may provide clues for the interpretation of the results presented here. There are clear spiral-like structures at large radius in M 33 (Fig. 1), with location and shape correlated with the bisymmetry in the velocity dispersion. Although non-axisymmetric perturbations are not necessary to have an anisotropic velocity dispersion ellipsoid, we can speculate that the spiral arms in M 33 could enhance the velocity anisotropy parameter in the outer disk regions if part of the HI gas in M 33 behaves like a collisionless medium. We also note that disks of similar stellar mass to that of M 33 show −0.1 ≲ β_ϕ ≲ 0.2 (Chemin 2018), thus stellar orbits in these galaxies are more isotropic.

The analysis also shows that it is not possible to strongly constrain σ_z given the small scatter of β_ϕ as a function of σ_z. It is nevertheless worth mentioning that under realistic hypotheses on σ_z, the range of σ_z/σ_R found for gas in M 33 is consistent with the value of ∼0.65 found for late-type stellar disks by Pinna et al. (2018).

Comparisons can be made with numerical simulations as well. Hydrodynamical modeling of gas in simulated spiral galaxies shows velocity dispersion components that are anisotropic (Bottema 2003; Agertz et al. 2009). In these numerical models, the tangential component is smaller than the radial dispersion. This is in agreement with most of our measurements. The simulations of Agertz et al. (2009) are interesting for our study because they are supposed to simulate a disk with similar physical properties to M 33. They showed that the planar dispersion, σ_P, given by the root mean squared value of σ_R and σ_ϕ, is twice larger than σ_z. If we restrict the comparison to the radial range 4−7.5 kpc where M 33 shows a more significant anisotropy parameter in the framework of an axisymmetric and anisotropic velocity model, the simulations of Agertz et al. (2009) show σ_z within 3.5−7 km s⁻¹ once the spiral-like features are well defined in the simulated density map. Our models show that σ_P ∼ 2σ_z in M 33 for σ_z ∼ 6 km s⁻¹. Therefore, comparable planar and vertical dispersions are found in both the observations and simulations. More broadly, within the range of vertical dispersions that has been investigated here, we find σ_P from ∼10 km s⁻¹ (σ_z = 10 km s⁻¹) to ∼13 km s⁻¹ (σ_z = 4 km s⁻¹), hence σ_P/σ_z from ∼0.9−1 to ∼3−3.2, respectively. The vertical motion is thus the main driver of the ratio. For σ_z = ⟨σ_los⟩_ϕ (∼8 km s⁻¹), σ_P/σ_z ∼1.4−1.5, which is close to the value expected for an isotropic velocity ellipsoid ($\sqrt{2}$).

Agertz et al. (2009) then found that σ_ϕ/σ_R is roughly consistent with the expectation of the epicyclic approximation (EA) which considers collisionless orbits that only slightly deviate from circularity. The EA stipulates that σ_ϕ/σ_R is related to the slope of the circular velocity v_c (e.g., Binney & Tremaine 2008):

$$\begin{array}{c}\hfill {({\sigma }_{\varphi }/{\sigma }_{\mathrm{R}})}_{\mathrm{EA}}^2=\frac{1}{2}(1+{\displaystyle \frac{\mathrm{d}ln{v}_{\mathrm{c}}}{\mathrm{d}lnR}})·\end{array}$$(4)

The epicycle anisotropy of M 33, ${\sigma }_{\varphi }^2-{\sigma }_{\text{R}}^2{\beta }_{\text{EA}}=1-{\left({\sigma }_{\varphi }/{\sigma }_{\text{R}}\right)}_{\text{EA}}^2$ is shown as a dashed blue line in Fig. 6. It was derived using the model of the HI rotation curve of M 33 by Koch et al. (2018) as a proxy for the circular velocity. The assumption that the circular velocity can be approximated by the tangential velocity is reasonable for gas, except maybe in lower mass disks (Dalcanton & Stilp 2010). The value of β_EA shows very little variation as a function of radius. Interestingly, the inner R = 4 kpc of M 33 is the range of radii where the anisotropy is mild (Fig. 6) and differs by ≲0.25 from that derived from the epicyclic approximation (except for the two outlying dips at R ∼ 0.7 and 2.1 kpc). The velocity anisotropy of gas is systematically larger than expected from the epicyclic approximation beyond R ∼ 5 kpc, however. Equation (4) cannot allow the observed strongly radial HI orbits because the rotation curve of M 33 is barely rising at these radii (β_EA ∼ 0.4). We would need a decreasing rotation curve to get β_EA >  0.5 compatible with the observed velocity anisotropy in the plane. Our results thus do not fully agree with those found with the numerical simulations of Agertz et al. (2009) in the disk regions of stronger radial bias. On the other hand, if the ellipsoid of velocity is isotropic, then the epicyclic approximation is violated as well because a linearly rising velocity curve with radius is strictly required for an EA–isotropy agreement. This indicates a failure of the epicycle approximation of orbits for the gas component. Interestingly, this discrepancy is reminiscent of the result found in Chemin (2018) for stellar disks where the diversity of stellar orbits could not be reproduced by the theory.

Is the HI velocity dispersion ellipsoid of Messier 33 anisotropic or isotropic? Choosing a side for this question is not an easy task as the two models both explain the major asymmetry. The Fourier analysis has the flexibility to probe a large number of modes, hence isotropic and asymmetric models are unsurprisingly more accurate in modeling σ_los. An intriguing result of this work is the observation that the orientation of kinematic perturbations are often aligned with the principal disk axes (modulo half and full periods) as if a projection effect affected significantly the velocity dispersion map of M 33. Such coincidences may be fortuitous in M 33 given that asymmetries in the gas density often appear coincident with the dispersion asymmetries (Fig. 1). This also suggests that further analyses of phases of asymmetric modes inside velocity dispersion maps are promising to assess the nature of gaseous velocity ellipsoids. If the phase angles were found systematically near the principal axes, at the positions of supposedly stronger velocity anisotropy, then it would rule out isotropic velocity ellipsoids since that occurrence should occur only rarely for a random distribution of phase angles of perturbations in the isotropic scenario. Data from deep surveys like THINGS (Walter et al. 2008), HALOGAS (Heald et al. 2011), and LittleTHINGS (Hunter et al. 2012) for the neutral atomic gas, or PHANGS-ALMA (Sun et al. 2018) for the molecular gas will be helpful to study that problem.

6. Summary

Messier 33 has a non-axisymmetric distribution of observed random motions of HI gas. There is a prominent pattern that makes the HI velocity dispersion weaker near the major axis and stronger near the minor axis of the galaxy. The velocity dispersion of the R >  4 kpc disk can locally be larger by up to 60% than the azimuthally averaged value. Hypotheses allying axisymmetry and isotropy are ruled out to explain the variations of the velocity dispersion.

Among the models presented in this study a Fourier transform has shown that bisymmetric random motions having an amplitude of up to 2 km s⁻¹ (25% of the axisymmetric value) must be invoked to explain the discrepancy while maintaining isotropy. It dominates the harmonic asymmetries in the random motions, and first-, third-, and fourth-order motions were found mostly weaker. The phase angles of the asymmetries are often seen close to the principal axes of the HI disk, and particularly the bisymmetry, which is aligned with the minor axis. The asymmetries coincide well with the non-axisymmetric spiral-like distribution of HI gas in M 33.

Another model was to consider that the velocity dispersion ellipsoid is axisymmetric but anisotropic, acting as if part of gas behaved like a collisionless medium. That led us to constrain the radial and tangential components of the ellipsoid (σ_R and σ_ϕ) at fixed vertical dispersion σ_z, and from which the azimuthal velocity anisotropy parameter β_ϕ = 1 − (σ_ϕ/σ_R)² could be measured. In the framework of this axisymmetric model, β_ϕ is mostly positive and maximum at R ∼ 6 kpc, indicating orbits of the atomic gas that are strongly radial. The perturbed dynamics in the spiral-like structure could be responsible for the velocity anisotropy. It was also found that while anisotropic velocity dispersions could be measured when the HI gas is treated like a CNM, it is not the case in the outer regions of stronger velocity anisotropy when the gas is seen as a WNM. A high thermal component does not leaves as much room for the planar motions.

As for the CO gas traced by a collection of a few hundreds of molecular clouds, β_ϕ is highly scattered and did not allow us to draw firm conclusions about the shape of CO cloud orbits. Although the comparison with HI gas remains limited because cloud-based dispersions were used in this analysis, unlike HI dispersions this result is not surprising if velocity dispersions of molecular clouds are driven by local cloud dynamics.

These results were found only marginally dependent on the assumptions made for σ_z (chosen as representative of values observed in face-on nearby disks). In future works, we will pursue the analysis of the properties of asymmetries in gas velocity dispersions by means of a larger sample of disk galaxies via further sensitive HI and CO measurements, and investigate to what extent anisotropic velocity ellipsoids can still explain the asymmetric gas random motions.

Acknowledgments

We are very grateful to an anonymous referee for insightful propositions which improved the analysis and the content of the article. This research was supported by the Comité Mixto ESO-Chile and the DGI at University of Antofagasta, and J. Braine by the MINEDUC-UA project code ANT 1755.

References

All Figures

Fig. 1. Observational data of Messier 33. Left panel: HI column density map (VLA, 100 pc resolution, logarithmic stretch) with Spitzer/IRAC 3.6 μm stellar distribution overlaid (gray contours, showing the 0.1, 0.3, 0.5, 0.7, 0.9, 1.4, and 2.5 MJy sr−1 levels). Right panel: observed HI velocity dispersion (VLA, 100 pc resolution, not corrected for the instrumental dispersion), with HI column densities contours (1, 2, 3 × 1021 cm−2). The dashed line represents the location of the major axis of the inner unwarped disk, with position angle of 22.5°. Concentric ellipses show the projected locations of R = 2, 4, 6 kpc.

In the text

	Fig. 2. Azimuth-velocity dispersion diagrams of atomic gas at selected radii in M 33. Each point corresponds to an individual measurement, i.e., the dispersion of each pixel within the map. Solid lines represent the results of a Fourier Transform model of the HI dispersion map (Sect. 3, assuming a cold neutral medium case). The dispersions are from the VLA 100 pc resolution data.
In the text

Fig. 3. Results of the discrete FFT of the HI velocity dispersion map of M 33 (100 pc resolution, in the CNM case). The amplitudes (upper panel) and phases (middle panel) of the Fourier modes are shown as black open circles (k = 0), green open circles (k = 1), filled circles (k = 2), blue upward filled triangles (k = 3), and red open diamonds (k = 4). The reference ϕ = 0 is chosen along the semi-major axis of the NE approaching half of M 33. Bottom panel: σ2/σk, using the same colors and symbols as above for the k = 1, 3, 4 modes. The range of ratios is chosen up to 10 for clarity, but the total range is 78 for σ2/σ1 and σ2/σ4, and 12 for σ2/σ3.

In the text

Fig. 4. Results of the anisotropic and axisymmetric velocity model for the atomic gas in Messier 33 (100 pc resolution data). Results are those of the one-parameter model and show the profile of σR (σϕ, respectively) as filled green circles (open diamonds), obtained assuming σϕ = σz = ⟨σlos⟩ϕ (σR = σz = ⟨σlos⟩ϕ). A dashed blue line is for ⟨σlos⟩ϕ (corrected from instrumental dispersion). Results obtained assuming a null thermal component. The shaded areas highlight the regions where the k = 2 mode was found to be weaker than other dispersion asymmetries (see Sect. 3).

In the text

	Fig. 5. Same as Fig. 4, but for the two-parameter model with σR and σϕ as free parameters.
In the text

Fig. 6. Profiles of azimuthal velocity anisotropy βϕ of the atomic gas in M 33 at different angular resolutions. Solid lines show the anisotropy profiles obtained assuming σz = ⟨σlos⟩ϕ and a null thermal component. Short-dashed lines show the ±1 rms errors; a dotted line is an illustration of result choosing another vertical dispersion, σz = 4 km s−1 (70 pc resolution only); and a blue long-dashed line is the velocity anisotropy profile expected from the epicycle theory (Sect. 5). Shaded areas are as in Fig. 4.

In the text

Fig. 7. Same as Fig. 2, but using HI data from VLA (70 pc resolution), DRAO (490 pc resolution), and CO data from IRAM 30 m antenna (50 pc resolution) to show results of the anisotropic velocity dispersion model (Sect. 4). The columns illustrate various cases of velocity anisotropy parameter found by the best-fit model: βϕ ∼ 0 (for isotropy, left ), 0 <  βϕ (for radial bias, middle), βϕ <  0 (for supposedly tangential bias, right). For each resolution, the selected radii are among the best positions illustrating each case of anisotropy parameter, following Fig. 6. Solid lines represent the best-fit dispersion model at the considered radii, assuming σz = ⟨σlos⟩ϕ.

In the text

	Fig. 8. Profile of azimuthal anisotropy βϕ of the velocity dispersion of the atomic gas in M 33 obtained using the 100 pc resolution VLA data by assuming σz = ⟨σlos⟩ϕ and HI as a warm neutral medium (thermal component of 6 km s−1). Shaded areas are as in Fig. 4.
In the text

Fig. 9. Results for the molecular gas in Messier 33. Top: profiles of σR and σϕ of CO gas, obtained assuming σz = ⟨σlos⟩ϕ (dashed blue line, corrected for instrumental dispersion) and a null thermal component. Filled circles are for σR and open diamonds are for σϕ. Bottom: profiles of azimuthal velocity anisotropy βϕ of CO gas. The solid line shows the anisotropy obtained assuming σz = ⟨σlos⟩ϕ and a null thermal component, and dashed lines show the ±1 rms errors. For comparison, the anisotropy of the atomic gas as derived from all pixels in the 100 pc resolution velocity dispersion map is shown as a dotted line, and that derived using only the pixels at the positions of the molecular clouds as open triangles. Filled circles are values derived by masking a few deviant dispersions of CO clouds (σlos >  5 km s−1). Shaded areas are as in Fig. 4.

In the text