© The Authors 2023

1. Introduction

In discs of galaxies, large-scale perturbations such as bars, spiral arms, warps, or lopsidedness make the orbits non-circular. They generate streaming motions, which are asymmetric radial and tangential velocities in the stellar and gaseous media in galactic discs (e.g. Visser 1980). This is beautifully illustrated from an observational viewpoint via the kinematics of millions of individual stars in the Large Magellanic Cloud (LMC), for which the Gaia Collaboration used astrometry to deduce the velocity fields of both the radial and tangential components of a disc galaxy (Gaia Collaboration 2016, 2018a,b, 2021a), and from which significant asymmetries are observed in the velocity maps, as caused by the bar and spiral arms of the LMC (Gaia Collaboration 2021b).

The study of kinematic asymmetries has become relevant in galactic dynamics, as it impacts our knowledge of the structure of dark matter halos and the comparison with simulations made in a cosmological context. Interstellar gas is thought to follow circular orbits as it undergoes dissipative collisions. If asymmetries are important, then there is no guarantee that rotation curves perfectly trace the mass distribution of galaxies. Asymmetries could affect the inner slope of gaseous rotation curves and require a model including asymmetries to determine the scale parameters of dark matter halos (Hayashi & Navarro 2006; Oman et al. 2019). Strategies have been developed to quantify asymmetries in gaseous velocity fields in the last 25 years. One approach is to measure harmonic velocity components that overlap with the axisymmetric motions (Schoenmakers et al. 1997; Krajnović et al. 2006). Early works focussed on the large-scale kinematic lopsidedness in H I velocity fields, a perturbation of first order that causes the kinematic centre to drift with radius. In this case, rotation curves from the approaching and receding halves are not consistent with each other. This could be due to a misalignment between the axes of the disc and the host dark matter halo (Swaters et al. 1999). Trachternach et al. (2008) estimated asymmetries for a large sample of H I velocity fields to study the elongation of the gravitational potential. They found that the asymmetries do not alter the mass distribution inferred from the axisymmetric rotational motions. Spekkens & Sellwood (2007) measured bisymmetric perturbations in CO and Hα velocity fields to highlight the importance of bisymmetric flows caused by a stellar bar on the inner shape of the rotation curve. Chemin et al. (2016) modelled the harmonics observed in the Hα velocity field of a grand design spiral through direct derivation of gravitational potentials and assessed the impact of asymmetries on the structure of the host dark matter halo.

While studying streaming motions in ordered velocity fields has become routine, little is known about asymmetries in maps of apparently random motions of interstellar gas, as traced by the line-of-sight velocity dispersion (σ_los). The velocity dispersion can be used to assess the dynamical heating of discs, and to estimate the support due to pressure and how it compares with rotational motions (e.g. Combes & Becquaert 1997; Koyama & Ostriker 2009; Bershady et al. 2010; Oh et al. 2015). The main reason dispersion has not been studied extensively is that σ_los mixes turbulence, asymmetric drift, shear, superposition of multiple elliptical orbits, and unresolved motions on relatively large scales. The local velocity dispersion ellispoid of gas is also usually considered isotropic because the gaseous component dissipates energy through collisions between clouds in all directions. Within this hypothesis, no projection effect occurs in the random motions of gas, as traced by σ_los, unlike that of collisionless stars. Furthermore, the projection of the velocity ellipsoid can be degenerate if the ellipsoid is anisotropic. The finite velocity and spatial resolution further complicate the interpretation of maps of σ_los.

In a recent study of the Local Group spiral M 33, Chemin et al. (2020) found perturbations in the random motions through Fourier transforms of σ_los maps of 21-cm H I data at various angular resolutions. Bisymmetry dominated in the outer half of the disc with a phase aligned with the minor axis of M 33. The question of whether this is a typical characteristic of spiral galaxies in general remains open. One possibility could be that these observations are due to the influence of an anisotropic velocity ellipsoid, which harbours a dominant radial component, as suggested by these authors. The goal of this study is to determine whether what was observed in M 33 is a general feature of spiral discs.

Our goal is to systematically search for asymmetries in velocity dispersion maps, identify their origin, and evaluate the consequences for galactic dynamics. In this paper, we present the first census of asymmetries in the velocity dispersion of H I gas in local, massive discs. The samples used for this analysis are described in Sect. 2. The asymmetries are measured through fast Fourier transforms (FFTs) of velocity dispersion maps from interferometric data with a careful treatment of the effect of beam smearing (BS hereafter), as detailed in Sect. 3. The properties of the asymmetries are presented in Sect. 4, and the discussion of their origin is given in Sect. 5. Finally, we provide a synthesis and conclusions in Sect. 6.

2. Selection of a working sample

Several factors limit the modelling of spatially resolved velocity dispersions. The smearing induced by the finite telescope beam can produce systematic asymmetries in velocity dispersion maps (see Sect. 3). For a given instrumental configuration of spatial and spectral resolution and sampling, the BS effect increases with the distance and inclination of the galaxies.

For these reasons, our sample is from The high-resolution H I Nearby Galaxy Survey (THINGS, Walter et al. 2008). With galaxy distances between 2 and 15 Mpc, THINGS yields linear resolutions up to ∼900 pc. The velocity resolution is 2.6 or 5.2 km s⁻¹, and the pixel size 1.5″, except for the galaxy NGC 2403 (1″). We have used the three moment maps of the integrated H I emission (0th moment), line-of-sight velocity and velocity dispersion (1st and 2nd moments, respectively), which are made available by the THINGS collaboration¹. In particular, we used the data obtained with natural weighting. We however verified that the conclusions of this analysis are not changed by using higher resolution kinematics, from maps obtained with robust weighting. A discussion of the impact of the angular resolution can be found in Sect. 5.2.3.

To perform the Fourier analysis of velocity dispersion fields, robust constraints on the geometry of the discs are necessary, and particularly the variation of the inclination and the position angle of the discs as a function of galactocentric radius. The disc warping parameters for 19 of the 34 THINGS galaxies have been measured with tilted-ring models of the velocity fields by de Blok et al. (2008). Among several methods applied to obtain velocity fields from the H I data cubes, these authors adopted results from Hermite h3 polynomials fittings, due to their stability at low signal-to-noise ratio (S/N). They explained that two masks were successively applied to filter low-quality regions from the velocity fields. The first one consisted in rejecting H I profiles (i) for which the fitted maximum intensity was lower than 3σ_ch, with σ_ch being the average noise in the profile outside of the emission line, and (ii) for which the dispersion of the fitted function was lower than the channel separation. The second mask consisted in a sigma-clipping on the H I column density maps to suppress noise pixels and to exclude regions inside r_min and outside r_max, which are respectively the innermost and outermost radii of the tilted-ring models fitted of the final velocity fields. They finally measured the warp parameters and rotation curves, adopting a sampling of two points per synthesised beam size.

We masked the moment maps following de Blok et al. (2008) and adopted the parameters of their tilted-ring models. This makes our gaseous and kinematic distributions, and the galaxy cylindrical frames fully consistent with their analysis. A minor difference with their study is that we considered sampling the radial profiles with one data-point per beam only, to minimise the correlation between adjacent rings. Finally, among the 19 galaxies, we selected the 15 massive, regular galaxies, excluding dwarfs or strongly disturbed galaxies (NGC 2366, DDO 154, IC 2574 and NGC 4826).

Chemin et al. (2020) showed that the angular resolution did not affect the detection of the bisymmetry and/or anisotropy observed in H I velocity dispersion maps of M 33, as it was seen at 70, 100, and 490 pc resolution. To increase the size of our sample, we included sources from the Westerbork H I Survey of Spiral and Irregular Galaxies (WHISP; van der Hulst et al. 2001). The reduction pipeline of the WHISP observations provides data at resolutions of 14″ × 14″/sinδ, where δ is the galaxy declination, 30″ × 30″, and 60″ × 60″, for a spectral resolution of 2.06 to 16.5 km s⁻¹. At distances out to more than 100 Mpc, the implied linear resolution for many WHISP galaxies is beyond the kiloparsec scale. Considering the balance between spatial resolution and sensitivity, and, due to the necessity to work with an appropriate spatial resolution, we only selected the sources with a physical resolution better than 1.5 kpc (Adamczyk 2021). This corresponds to lower S/N measurements than for THINGS targets, and consequently to a smaller extent of the gas distribution and kinematics. We further required the neutral gas to cover seven times the size of the beam (as suggested in Bosma 1978, to derive rotation curves), in the same inclination range as THINGS, and with no obvious signs of tidal interactions. In the end, 17 WHISP galaxies were added to the THINGS sample. We produced moment maps for them by means of CAMEL, a Python tool² described in Epinat et al. (2012). We then performed tilted-ring models of WHISP data cubes using the package 3D-Based Analysis of Rotating Object via Line Observations (3D-Barolo, Di Teodoro & Fraternali 2015). This software creates a model data cube from input centre of the rings, inclination and position angle, systemic velocity, scale height of the disc, and velocity dispersion, which is fitted to the observations. To determine the kinematic centre, we first kept the inclination and position angle fixed, with the centre of radial rings and the systemic velocity as free parameters. Then, we fixed the centre and let the inclination and position angles vary to get the warp parameters.

The final sample of galaxies is presented in Table 1 which lists the coordinates, distance, systemic velocity, the mean inclination and position angle of the discs, and the properties of the synthesised beam. With mean source distances of 8 and 9 Mpc for THINGS and WHISP respectively, the mean beam B_maj of 11″ and 18″ correspond to average linear resolutions of 470 pc and 80 pc, respectively. The THINGS sample is the main reference and the WHISP sample is used to check the observed trends are not linked to a particular telescope. Section 4 shows that the WHISP data fully support the trends found with THINGS data only.

3. Methodology

Observations yield the integrated flux of gas emission lines, the line-of-sight velocity and velocity dispersion. This observed dispersion cannot be modelled directly as it results from the convolution of the random motion in the plane of galaxies (σ) with the Line Spread Function of the instrumental device (σ_LSF) and the broadening from thermal processes occurring in the gas (σ_T). The observed dispersion is then expressed by:

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

Additionally, due to the finite angular resolution of observations, σ is affected by a smearing effect. This comes from unresolved velocity gradients in the line-of-sight velocity, which broadens emission lines.

Below we describe how we modelled the BS effect, and how a corrected velocity dispersion has been inferred. Then, we describe how Fourier transforms constrain the amplitudes and phases of asymmetries in the velocity dispersion maps.

3.1. Modelling the effect of BS in a velocity dispersion map

Beam smearing impacts the flux, velocity field, and velocity dispersion maps extracted from datacubes in a different way. Its main effect on velocity fields is to weaken the inner gradient of the rotation curve by mixing data from adjacent regions in the disc. de Blok et al. (2008) showed that BS does not have a significant impact on THINGS galaxies rotation velocities. They built mock datacubes, smoothed the cubes to the THINGS resolution, and derived the velocity fields and rotation curves from the mock datacubes. Compared to the input used to make the datacube, deviations were less than ∼1 km s⁻¹.

For velocity dispersion maps, one expects BS to be more prominent where large variations of radial velocities are observed locally, for instance due to rotation curve inner gradients and to variations of cos(θ), θ being the azimuth angle (see Eq. (C.1)). We need to model the impact of BS on velocity dispersion maps by taking into account the variations of velocities in the two spatial dimensions. We use the tool MocKinG³ based on an analytical formula presented in Epinat et al. (2010) who derived the observed velocity dispersion in a galaxy following:

$$\begin{array}{c}\hfill {\sigma }^2={\sigma }_{\mathrm{corr}}^2+{\sigma }_{\mathrm{bs}}^2,\end{array}$$(2)

with

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{corr}}^2=\frac{\int {\int }_{\mathrm{pix}}\left[{\sigma }_{\mathrm{loc}}^{2}M\right]\otimes \mathrm{PSF}\mathrm{d}s}{\int {\int }_{\mathrm{pix}}M\otimes \mathrm{PSF}\mathrm{d}s},\end{array}$$(3)

and

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{bs}}^2=\frac{\int {\int }_{\mathrm{pix}}\left[V_{\mathrm{los}}^{2}M\right]\otimes \mathrm{PSF}\mathrm{d}s}{\int {\int }_{\mathrm{pix}}M\otimes \mathrm{PSF}\mathrm{d}s}-{\left(\frac{\int {\int }_{\mathrm{pix}}\left[V_{\mathrm{los}}M\right]\otimes \mathrm{PSF}\mathrm{d}s}{\int {\int }_{\mathrm{pix}}M\otimes \mathrm{PSF}\mathrm{d}s}\right)}^2,\end{array}$$(4)

where σ_loc is the local velocity dispersion map, V_los is the line-of-sight velocity field, M is the line flux map, ⊗PSF represents the two-dimension convolution by the PSF, and ∬_pix ds integrates over the surface of the pixel. Equation (3) shows that σ_corr is impacted by BS on the local velocity dispersion σ_loc whereas Eq. (4) shows that unresolved velocity shears in the first moment map create an artificial line broadening. Our goal is to study σ_corr by correcting the observed velocity dispersion σ for the term σ_bs of Eq. (4). The maps of first two moments (M and V_los) in this equation should be at high-resolution and free from BS to properly account for velocity variations inside the beam and pixels. In practice, such high-resolution maps do not exist so the observed flux is usually used for M, and either a high-resolution model (see e.g. Epinat et al. 2010) or an observed velocity field are used for V_los. We used MocKinG with the observed H I flux maps and velocity fields to compute the BS correction. It is thus asymmetric by construction and accounts for both circular and non-circular motions from the kinematics. This leads to a fair correction, although not perfect, and is discussed in Sect. 5.2.3, 5.2.4, 5.3.1, and Appendix B. These BS dispersion maps were modelled assuming a bi-dimensional Gaussian synthesised beam (see Table 1). The effect of the observed dirty beams, which are not entirely elliptical, is addressed in Sect. 5.2.2.

An example of BS modelling with the galaxy NGC 2841 is shown in Fig. 1. This fast-rotating galaxy (second panel) has an inclination of ∼74°. The fourth panel shows the corresponding σ_bs, which, once subtracted to the observed dispersion (middle panel), yields the corrected velocity dispersion map (right panel) from which asymmetries are measured. The general behaviour of the smearing effect is thus that σ_bs is larger at low radius and exhibits a typical X-shape pattern. On average, σ_bs decreases with radius, but non-negligible values are also observed near non-axisymmetric density and velocity features (see Sect. 5.2.1), underlining the benefits of modelling the smearing from 2D data.

Fig. 1. Density and velocity maps of NGC 2841. From left to right we show: the observed flux density map, the observed velocity field, the observed velocity dispersion, the BS model, and the velocity dispersion map corrected from the BS effect.

3.2. Fourier series modelling of velocity dispersion maps

In order to quantify asymmetries, we performed FFTs of the maps of H I random motions. Discrete FFTs of the corrected ${\sigma }_{\text{los}}^2$ maps were performed using the FFTPACK package of the SciPy libraries in Python (see Appendix A). In cylindrical coordinates, the discrete FFT of the velocity dispersion is:

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{asym}}^2={\displaystyle \sum _{k=0}^{N-1}{\sigma }_k^{2}cos(k(\theta -{\varphi }_k)),}\end{array}$$(5)

where σ_k and ϕ_k are the amplitude and phase angle of an asymmetry of order k, N is the number of orders in the decomposition, equivalent to the number of elements along the considered ring, and θ is the azimuthal angle in the galaxy plane measured from the semi-major axis of the receding disc half. With this FFT formalism, we make the choice of studying the azimuthal asymmetry, that is we cannot study the asymmetry with respect to the disc plane in the vertical direction. Nevertheless, the present study includes radial, azimuthal and vertical components projected on the sky plane, assuming either that contributions along the direction orthogonal to the disc plane are averaged across the line of sight, which is in practice more valid for face-on than for edge-on galaxies, or that discs are infinitely thin. We use the squared velocity dispersion because ${\sigma }_{\text{los}}^2$ is a simple sum of quadratic terms (Eqs. (1) and (2), see also Appendix C, Eq. (C.12)). Consistent results are obtained when FFTs of the linear dispersion are calculated.

In practice, we decomposed a galactic disc into a series of concentric rings whose geometry was defined by the tilted-ring models (Sect. 2). Given the different angular sizes of the galaxies, the adopted ring width of one full width at half beam power leads to differing numbers of rings per galaxy, from 7 for NGC 3627 to more than 140 for NGC 2403. Such a non-uniformity in the number of radial bins among galaxies has no visible consequence on the results (see Sect. 4).

For each dispersion map, we first subtracted ${\sigma }_{\text{LSF}}^2$ and ${\sigma }_{\text{T}}^2$ from the squared observed dispersion, as well as the corresponding 2D map of ${\sigma }_{\text{bs}}^2$ described in Sect. 3.1. All pixels with resulting negative quadratic velocity dispersion after these subtractions were discarded from the maps at this stage of the process, because such values are unphysical. The instrumental dispersion is σ_LSF = Δ_V/2.35 and σ_LSF = 1.2Δ_V/2.35 for the THINGS and WHISP galaxies, respectively, where Δ_V is the velocity channel width listed in Table 1. H I is a mixture of cool (∼100 K) and warm (∼5000 − 8000 K) gas. In M 33, the H I velocity dispersion was sometimes narrower than σ_T if the warm gas was assumed to dominate (Chemin et al. 2020). Hence, we do not consider σ_T here. We note that this has no consequence hereafter, because considering gas as a warm neutral medium is equivalent to subtracting quadratically σ_T ∼ 6 km s⁻¹ from the axisymmetric term σ₀. This latter term, measured as the mean dispersion of a given ring, is then subtracted quadratically from all pixel values inside the considered ring. Therefore, the observed ${\sigma }_{\text{asym}}^2-{\sigma }_0^2$ are centred on 0, and can be negative (see Eq. (5)). We also point out that σ_T could be asymmetric and vary over the disc. In such a case, fluctuations of σ_T as a function of the position have been accounted in the measurements of σ_asym, though it is not possible to disentangle these effects from those arising from other local motions of, for example, gravitational origins, without being able to measure locally the gas temperature.

We then sorted the squared velocity dispersions with increasing values of azimuth and apply the FFT, leading to a harmonic decomposition with N/2 terms, where N is the number of pixels in the considered radial bin. Incomplete coverage of azimuths as caused by missing pixels (not-a-number values) could seen by the FFT as artificial perturbations. In the THINGS sample, about 77% of the 883 available tilted rings show less than 1% of missing values of all available pixels, and 95% less than 4% of missing azimuthal angles (Fig. 2). The missing factor is thus low, and we verified that it has no impact in the analysis. Within WHISP data, we rejected 30% of the initial 203 rings because they had more than 40% of missing values. Within the 144 remaining rings, 15% of the pixels have missing values. Even though no impact was detected, we replaced the missing values by the azimuthally averaged dispersion before the modelling.

Fig. 2. Distribution of the fraction of missing values for the radial bins in the THINGS galaxies.

3.3. Generating toy models to study the impact of unresolved ordered velocity variations

As discussed in Sect. 3.1, BS is present in our data and is corrected in this study using observed line flux maps and velocity fields, which limits the accuracy of the correction. In order to study the impact of residual BS effects, we perform toy models for cases with barely resolved motions due to large-scale axisymmetric rotation of various strengths, with and without additional asymmetric velocity perturbations on local scales, under both isotropic and anisotropic hypotheses. We build toy models to produce data cubes, velocity, and σ_los fields, and calculate the FFT of σ_los. This enables us to generate different configurations and identify the conditions for which systematic values in the distributions of σ_k and/or ϕ_k could occur. Our toy models have no dynamical basis, yet they are very useful to assess the effects of BS, anisotropy, and streamings in σ_los. The full description of the toy models is presented in Appendix B.

We produced mock datacubes of 400 × 400 pixels (scale of 1″, or ∼50 pc at a typical distance of 10 Mpc), with 200 spectral elements with a 3 km s⁻¹ velocity sampling using 5 × 10⁶ uniformly distributed points. To first order, the galaxy is assumed to be an axisymmetric rotating disc to which velocity perturbations can be added. Two rotation curves, one with a weak velocity gradient and a moderate velocity plateau, the other with a steep inner gradient quickly reaching a high velocity plateau, were used. Sharp planar velocity perturbations are produced by a bisymmetric spiral pattern, with five possible inner angles. No lopsidedness (k = 1 mode) is introduced for simplicity. Therefore, the perturbed models have intrinsically anisotropic velocities on large scales, as illustrated by the ellipsoid elongated in the radial dimension in Appendix B. We also synthesise velocity anisotropy in the toy models by modifying the shape of the velocity ellipsoid locally, that is by choosing a uniform radial bias σ_θ = 0.7σ_R with null covariance, with σ_R = 8 km s⁻¹ and σ_z = 5 km s⁻¹, in opposition to isotropic velocity distributions (σ_R = σ_θ = σ_z = 8 km s⁻¹). The choice of σ_θ = 0.7σ_R comes from the radial bias seen of young stellar populations in the disc of the Milky Way (Gaia Collaboration 2023). This is an additional effect to the streaming-driven velocity anisotropy, and switching it off and on in the toy models enables to assess other specific sources of anisotropy not accounted for by the streaming perturbations in the planar velocity components.

For each particle, velocities are drawn randomly with a mean velocity and a velocity dispersion for each component. The mean in radial, azimuthal and vertical velocities are computed from the rotation curve and the perturbation. The velocities are then projected along the line-of-sight for three inclinations (45°, 60°, and 75°). For each inclination, cubes were created for 12 orientations of the spiral perturbation. Gaussian smoothing was then applied to mimic BS (8 pixels FWHM, corresponding to ∼400 pc), before extracting moment maps. Beam-smearing corrections were applied as for observed data. Because of projection effects, a perturbation along the radial (azimuthal) direction is mainly seen along the minor (major) axis in the line-of-sight velocity fields and velocity dispersion maps (see Fig. B.2). FFTs of ${\sigma }_{\text{los}}^2$ were calculated for 35 independent rings, yielding profiles and distributions of σ_k and ϕ_k, as in Sect. 4 using the 36 models (3 inclinations, 12 orientations of the perturbation). Figure 3 shows a velocity perturbation (case where the velocity perturbation vector is orientated towards the closest position of the spiral) for the steep rotation curve. From the maps, we can see residual BS effects due to both large-scale rotation and to the unresolved perturbation. The velocity perturbation is more orientated azimuthally in the centre than in the outer parts. The effect of the radial bias is also clearly seen along the minor axis in both the velocity dispersion map and in the histograms of the second order phases. Other cases are shown in Appendix B. An analysis of BS residuals induced by large-scale rotation is provided in Sect. 5.2.4 and the analysis of projection effects from asymmetric perturbations in Sect. 5.3.1.

Fig. 3. Maps and histograms of toy models with a steep rotation curve (v0 = 250 km s−1, rs = 20″ and γ = 2) with a spiral perturbation. We show (top), from left to right, the density map, the velocity field, and the BS corrected dispersion maps in the isotropic and anisotropic cases, obtained with an inclination of 60° and with ϕsp = 135°. The bottom panels show amplitude and phase histograms of orders k = 2 and k = 4. Blue and green histograms correspond to the cases with local uniform isotropy with large-scale asymmetries, and local uniform anisotropy with large-scale asymmetries, respectively.

4. Asymmetries in velocity dispersion maps

Coefficients of the FFTs were calculated for the 883 and 144 radial rings of THINGS and WHISP sources, respectively. Of the first 20 orders, the dominant asymmetries are those up to the fourth harmonic, in agreement with Chemin et al. (2020), so we limit the analysis to k ≤ 4 hereafter.

4.1. A case study: NGC 2841

The results of the Fourier analysis for all galaxies from the THINGS sample are presented in Appendix D. To illustrate examples of results, Fig. 4 shows the process and results for NGC 2841, from ${\sigma }_{\text{corr}}^2$ to the FFT and the individual FFT components. The normalised phase angles ϕ_k/T_k, where T_k = 2π/k, are shown, with values of 0.5 and 1 corresponding to half and a full period. The normalised phases ϕ₁/T₁ = 0, 0.5, 1 and ϕ_k/T_k = 0 for other orders are aligned along the major axis of the galaxy, and ϕ₂/T₂ = 0.5 along the minor axis. In NGC 2841, ${\sigma }_{\text{corr}}^2$ is high along a cross shaped pattern, with strong k = 2 and k = 4 perturbations and the k = 2 perturbation grows stronger with radius. The residual ${\sigma }_{\text{res}}^2$ is weak, showing that the first 4 orders reproduce the structure in ${\sigma }_{\text{corr}}^2$. The galaxy is lopsided at small radii (e.g. Baldwin et al. 1980) and this is detected by the k = 1 mode of the FFT. The phase angle of the bisymmetry does not vary much beyond R ∼ 170″, at a value of ∼0.7 times the period of the k = 2 asymmetry.

Fig. 4. Examples of FFT results with the galaxy NGC 2841. Top: observed squared velocity dispersion and its modelling through FFT up to the order 4, and their residuals. Middle panel: individual orders of the FFT projected in the plane of the galaxy. Bottom: amplitudes (left) and phase angles (right) of FFT coefficients as a function of the radius.

4.2. Census of asymmetries in H I random motions

4.2.1. General trends

Figure 5 presents the distributions of the amplitudes and normalised phase angles of the Fourier modes for our reference THINGS sample (blue histograms), and the WHISP sample appended to the THINGS sources (green histograms). Table 2 lists the mean, median and standard deviation of the amplitudes for the 883 rings of the THINGS sample. Globally, the distributions of amplitudes have comparable shapes. With an average of 11 km s⁻¹, the strongest Fourier mode is of second order, and the weakest is of third order (8 km s⁻¹). Most k = 2 amplitudes are between 2.5 and 7.5 km s⁻¹ with a few rings below 2.5 km s⁻¹, possibly due to noise.

Fig. 5. Results of the FFT measurements of the H I velocity dispersion maps. The histograms show the number of rings as a function of amplitude (top row) and normalised phase angle (bottom row) of the FFT harmonics of the THINGS sample (883 rings, in blue), the interpolated distribution with 20 rings per galaxy for THINGS (280 rings, in red, normalised to the maximum values of the blue histograms), and the WHISP sample (144 rings, in green). Phase angles are normalised to the period Tk of each order.

For k = 1, we see in Fig. 5 (bottom) that the phase angle peaks at ϕ₁/T₁ = 0.45 and 0.95 with dips at 0.15 and 0.75. The bisymmetric mode is maximal around ϕ₂/T₂ = 0.65, and is minimum at 0. The k = 3 distribution is relatively flat with a maximum located at ϕ₃/T₃ = 0.15. The k = 4 histogram is strongly peaked at ϕ₄/T₄ = 0.5. The uncertainties are approximately the bin size of the distributions (0.1 T_k).

The distributions of Fig. 5 may be biased by a few galaxies with more radial rings than others. To assess the impact of different numbers of rings, we measured new profiles of amplitudes and phase angles using 20 equally spaced rings per galaxy as all galaxies but NGC 3627 have at least 20 independent rings. The initial profiles were interpolated at the 20 new radii. This yields new distributions built on 280 radii for 14 galaxies (NGC 3627 was excluded for this assessment). These histograms are shown as red lines in Fig. 5. The agreement between the blue and red histograms shows that the unequal number of radial bins has little impact on the results.

We can estimate the significance of the dips and peaks of the phase distributions. Assuming that the gas velocity ellipsoids are isotropic, the velocity dispersion maps should exhibit asymmetries randomly distributed over the plane of the sky. The 883 rings observed for the THINGS sample should thus yield, on average, N = 88 counts in each of the 10 bins of the histograms. If we quantify a confidence level as $\zeta =\sqrt{N}$, then the uncertainty is ζ ∼ 9. Keeping only the bins showing an excess of at least 3ζ, the peaks observed at 0.45 and 0.95 for k = 1 are detected at a level of 6.3ζ and 3.3ζ (147 and 119 counts, respectively), 5.5ζ at 0.65 for k = 2 (140 counts), 4.8ζ at 0.15 for k = 3 (133 counts), and up to 11ζ for the bin at 0.45 of k = 4 (194 counts). Now keeping only the bins showing a deficit of at least 3ζ, the minima observed at 0.15 and 0.75 for k = 1 correspond to a confidence level of 4.2ζ (49 counts), 3.7ζ at 0.95 for k = 2 (53 counts), and 4.5 to 5.3ζ at 0.05, 0.15, 0.25 and 0.95 for k = 4 (45, 42, 46 and 38 counts, respectively). Therefore, the systematic phase angles are not consistent with a random fluctuation of the orientations of asymmetries in the galaxies. We note that in the case of k = 3 no trend is seen around the unique bin exceeding 3ζ. It may be that this peak is caused by chance due to the limited statistics, and not by a systematic effect, unlike the peaks or dips seen in the other asymmetries. With the interpolated profiles at 20 bins for each galaxy (red histograms of Fig. 5), we find less significant peaks and dips than with the 883 rings due to lower statistics, but the trends are preserved.

4.2.2. Correlations between the Fourier modes

The upper row of Fig. 6 compares the Fourier amplitudes as a function of the phase angle ϕ₂/T₂ for the strongest orders (k = 1, 2, 4) for THINGS galaxies. The colour code is the number of tilted rings within bin widths of 0.05 for ϕ₂/T₂, and 1.75 for the amplitude difference. The density of tilted rings is highest around a difference of 0 km s⁻¹ (in agreement with values given in Table 2), irrespective of the value of ϕ₂/T₂, which implies correlated amplitudes. We measure a Pearson correlation coefficient of 0.8 ± 0.1 between the amplitudes of orders 2 and 4.

Fig. 6. Comparisons of amplitudes and phase angles between the Fourier k = 1, 2 and 4 orders, shown as differences of amplitudes versus ϕ2/T2 (upper row) and histograms of phase angles difference normalised to π (bottom row). The comparisons between the orders k = 1 and k = 2, k = 1 and k = 4, and k = 2 and k = 4 are shown in the left, middle, and right columns, respectively. In the upper row, the amplitude difference is colour coded with the logarithmic number of the tilted rings.

The bottom row of Fig. 6 compares the phase angle differences between orders Δϕ_m, n = ϕ_m − ϕ_n, and shows the differences within one period of the highest order, that is within ±T_m/2, with m > n. The distributions of Δϕ_2, 1 and Δϕ_4, 2 are peaked and symmetric around zero, implying a correlation between the k = 1 and k = 2 modes and between k = 2 and k = 4. These correlations suggest that k = 4 is a harmonic of the second order perturbation.

4.3. Correlations with the H I flux density

Several processes may induce correlations or anti-correlations between gas velocity dispersion and gas density. A large velocity dispersion can arise from unresolved motions, such as insufficient spatial resolution or the presence of unresolved multiple peak profiles. These unresolved velocity gradients may result from various physical processes, including steep density gradients, gas compression in density waves, instabilities associated with spiral-like features, or starburst outflows. The faintest regions may also display a high velocity dispersion because the S/N affects the profile widths, noisier profiles appearing broader. Leaving aside the aforementioned observational caveats (low resolution, low S/N) and focussing on the resolved regions, correlations between dispersion and density can have different origins. Large-scale star formation is clearly linked to H I content but this is not true for small scales (e.g. Zhou et al. 2018). As long as the H I density is below the density necessary to gravitationally collapse and form molecular hydrogen, the H I gas clouds have no reason to cool and the velocity dispersion to decrease. In some galaxies (e.g. NGC 4214, Wilcots & Thurow 2001), the largest Hα velocity dispersion is observed in the diffuse ionised gas regions which often has a low H I column density. Around massive star forming regions, the large velocity gradients might be explained by outflows. In low H I column density regions far from massive stars (e.g. NGC 2366, Hunter et al. 2001), in the absence of a local source, the velocity width might be induced by long-range turbulent pressure.

Using NGC 2841 as an example, Fig. 1 shows that the gas density (Σ_los) is principally distributed in bright rings and outer spiral arms, whereas larger velocity dispersion (σ_los) is along a centred cross-shaped structure. In other words, the gas density is either anti-correlated or correlated with velocity dispersion, sometimes fainter in large velocity dispersion regions (e.g. in the galaxy centre), brighter at low velocity dispersion (e.g. in the NE and SW quadrants), but also brighter at larger velocity dispersion (e.g. NW and SE quadrants). Across the sample, similar features are observed, and velocity dispersion seems correlated to spiral arms. Therefore, due to the wide diversity of processes, a galaxy-by-galaxy or a pixel-by-pixel comparison between the flux density and the velocity dispersion is beyond the scope of this paper. We rather aim at identifying global trends with respect to our analysis of the velocity dispersion.

The correlation between gas density and velocity dispersion has therefore been investigated for the whole sample using FFTs of the density maps, in the same way as for the velocity dispersion. The phase angle differences Δϕ_k ≡ ϕ_k(Σ)−ϕ_k(σ) are shown in Fig. 7. As in NGC 2841, we observe that all orders except k = 4 exhibit correlated (peaks) and uncorrelated (troughs) phases. An anti-correlation is observed for k = 1 (dip at Δϕ_k ∼ 0), and a correlation for k = 2 (peak at Δϕ_k ∼ 0). These findings suggest that the asymmetries observed in the velocity dispersion are related to those of the gas distribution that are mainly induced by spiral arms, bars, warps and lopsidedness.

Fig. 7. Histograms of phase angles (normalised by π) difference between gas density and velocity dispersion ϕk(Σ)−ϕk(σ) for k = 1, 2 and 4 (left, middle, and right panels, respectively).

4.4. FFT results versus galaxy properties

Our sample is somewhat biased by the fact that the smallest galaxies, which are usually the latest type and faintest ones, are also the closest⁴. To look for an effect, we divided the sample into two sub-samples of seven galaxies each around the median value of each parameters: morphological type, optical radius, absolute magnitude, and metallicity. Class I (II) is the sub-sample with the more distant (nearest), or the brightest (faintest), or the earliest (latest) type, or the largest (smallest) galaxies.

In Fig. 8, we show histograms of σ₂, ϕ₂, and ϕ₄ for classes I (red) and II (blue) of absolute magnitude, which is the parameter providing the largest difference between classes I and II. As expected, σ₂ is lower for small galaxies. For k = 2, the incidence of ϕ₂/T₂ ∼ 0.5 is low (high, respectively) for class I galaxies (class II). ϕ₂/T₂ in class II galaxies agrees with the observations of M 33 (Chemin et al. 2020), which is a class II galaxy. For orders k = 1 and 3, the phase does not depend significantly on the class. ϕ₄/T₄ has a bigger peak for class I than for class II but the difference is smaller for other parameters than the absolute magnitudes. The ϕ₄/T₄ distribution is almost flat for faint galaxies and peaked for bright ones. This could be due to residual BS (see Sect. 5.2.4) or because higher gas density contrasts are expected in earlier type spirals, inducing larger perturbations (see Sect. 5.3.1).

Fig. 8. Histograms of the velocity dispersion amplitudes σk (k = 2, left panel) and of the normalised phase angles ϕk/Tk (k = 2, middle panel and 4, right panel) for the 7 galaxies brighter and fainter than M = −20.7 mag (red and blue colours, respectively) for the THINGS sample. The sub-samples are the same as when splitting between high (red) and low (blue) projected maximum velocity (see Sect. 5.2.4). The grey histograms are the total distributions from Fig. 5.

5. Discussion

Our FFT analysis of the velocity dispersion of the H I gas shows that common values for ϕ_k/T_k are (1) ϕ₁/T₁ ≈ 0.45 or 0.95, (2) ϕ₂/T₂ ≈ 0.35 and 0.65 and (3) ϕ₄/T₄ ≈ 0.45, with no major structure in the k = 3 mode. These peaks correspond to major axis of the 1st order asymmetry and half period of the 4th order asymmetry. Interestingly, finding k = 2 asymmetries projected preferentially near the major axis is significantly less likely.

In the absence of instrumental effects such as BS at any scale, and assuming gas is isotropic, we expect the principal axes of the galaxies to be randomly orientated with respect to the spiral pattern or bar, so the distributions of ϕ_k/T_k should be uniform. In this section, we discuss the origin of the alignment of phases around particular angles.

5.1. Systematic effects in the FFT measurements and the tilted-ring models

5.1.1. Significance of FFT coefficients

In Appendix A, we describe the method used to infer the accuracy $\text{d}({\sigma }_k^2)$ of a measured amplitude ${\sigma }_k^2$ (Eq. (A.5)), and its significance $s_k={\sigma }_k^2/\text{d}({\sigma }_k^2)$. For this exercise, we assume a reasonable uncertainty dσ = 1 km s⁻¹ on the standard deviation of the line spread function. By keeping only the rings for which the significance is greater than 3, we find that amongst the 883 rings of the THINGS sample, 482 exceed this threshold for the order k = 1, 587 for k = 2, 431 for k = 3, and 536 for k = 4.

We further divided the rings into high and low significance samples using the median significance of each order. The median values of the significance s_k are 3.3, 5.2, 2.9, and 3.9 for k = 1, 2, 3 and 4, respectively. The resulting histograms are given in Fig. 9, with the low(high)-significance samples shown in red (blue, respectively). A first result is that apart larger tails in the distributions of high-significance samples, the amplitudes are not much affected by the division into two sub-samples. Since some high amplitudes rings are located in the innermost part of galaxies where the number of pixels is much lower than in the outskirt, we observed that some high amplitudes rings (red tail) are less significant than some of the low amplitude ones (blue distribution at low amplitudes). On the other hand, the low-significance samples show essentially random phases, except for order k = 4, while the high-significance show the trends identified before (ϕ₁/T₁ ≈ 0.45, ϕ₂/T₂ ≈ 0.65 and ϕ₄/T₄ ≈ 0.45) but more clearly. In other words, the rings of lower significance are consistent with having asymmetries randomly distributed in the sky plane (except for k = 4), while those of higher significance seem to be more representative of the systematically orientated asymmetries in the velocity dispersion maps. Moreover, the fact that both sub-samples have ϕ₄ centrally peaked around the half period is a hint that low amplitude residuals due to BS still affect the data, despite the correction applied before deriving the FFTs. We address this in Sect. 5.2.

Fig. 9. Histograms of amplitudes (upper row) and normalised phase angles (bottom row) resulting from the FFT analysis for the THINGS sample when splitting the sample in two, using the median significance of each order (see Sect. 5.1.1). The blue and red distributions are for the rings with a significance respectively greater and smaller than the median value.

5.1.2. Deviations in the tilted-ring model projection parameters

We now test the robustness of our results with respect to variations of kinematic centre, inclination and position angle. On the one hand, we changed the position of the kinematic centre of the rings by adding a random value selected from a Gaussian distribution of standard deviation of 1″ and 3″, respectively for the right ascension and declination, as expected from their average uncertainties (Trachternach et al. 2008). On the other hand, we ignored warps by using the mean inclination and position angle for each galaxy, as given in de Blok et al. (2008).

The new distributions of phases are then compared to those of Fig. 5, by averaging the positive differences of histogram heights of all bins. Overall, the trends described in Sect. 4 are still present. Varying the kinematic centre leads to differences of 10%, 7.2%, 5.7% and 7.7%, for k = 1, 2, 3 and 4, respectively. Not taking into account the disc warps has little impact as well, although the variations are more important, of 8%, 15.6%, 10% and 22.8% for k = 1, 2, 3 and 4, respectively. We can thus exclude the possibility that the systematic phase angles of the asymmetries in the H I velocity dispersion are caused by incorrect kinematic parameters.

5.2. Systematic signatures associated to BS

Beam smearing might provide residuals in the velocity dispersion beyond our correction process if locally the velocity gradient is not resolved as shown from Eq. (4), or if the radio beam is not well represented by the 2D Gaussian model described in Sect. 3.1. An incomplete BS correction could contribute to or create velocity dispersion asymmetries.

5.2.1. Properties of the BS pattern

To numerically explore the effect of BS, we first created a mock galaxy having both exponential flux and rotation curve models. In order to work with the same amplitude of the projected velocity field and thus with a similar BS induced velocity dispersion at zeroth order, we use a maximum projected velocity of 100 km s⁻¹ independent from the inclination. The Gaussian beam width is set to 10″ × 15″, and in order to sample correctly the exponential disc, we adopt a pixel size of 1″ and a scale-length of 50″. Figure 10 shows (left) the BS induced velocity dispersion maps at 30° (top) and 70° (bottom) inclinations. The central and right panels show the amplitudes and phase angles of the FFT decompositions. In our mock galaxy and PSF models, the BS pattern can be described using only orders k = 2 and k = 4. The ratio between the amplitude of these two orders vary with the inclination, with the k = 4 strength being predominant at large inclination. This latter is maximum around the radius where the maximum rotation curve is reached. The phases ϕ₂ and ϕ₄ are rather constant with radius, with their value depending on inclination. At fixed inclination, this constant value is related to the position angle of the galactic disc relative to that of the beam. In case both position angles match, or if the beam is circular, we have exactly ϕ₂ = 0.5 T₂ and ϕ₄ = 0.5 T₄. For an inclination of 70°, ϕ₂ and ϕ₄ are always in the range of 0.4 T₂ − 0.6 T₂ and 0.4 T₄ − 0.6 T₄ respectively, with ϕ₄ being close to its half period for most position angles.

Fig. 10. Results of the FFTs applied to the BS dispersion maps σbs, as built from disc models projected at 30° (top row) and 70° (bottom row). Left panel: map of σbs. The synthesised beam adopted in the toy model is shown on the bottom left corner. Middle panel: amplitudes of the Fourier modes for orders k = 0 to k = 4. Right panel: phase angles of the Fourier modes. A horizontal dashed line shows the half period of the Fourier modes. The disc scale-length is shown as a white contour (left) and as a vertical dashed line (middle and right).

Figure 11 shows that asymmetries in σ_bs are centrally peaked and lower than the observed asymmetries. The phase angles of σ_bs are similar to those seen in the galaxies. Subtle differences are observed, like the peak of probability at ϕ₂/T₂ ∼ 0.5 not observed in galaxies, or the stronger and narrower peak at ϕ₄/T₄ = 0.5. Finally, we also compared the results of the FFTs obtained with and without BS correction, and found negligible differences. These findings indicate either that our data are sparsely affected by BS, or that the BS correction is underestimated. In Sects. 5.2.3 and 5.2.4, we further investigate the impact of using higher resolution data to infer BS correction and the amplitude of BS residuals respectively.

Fig. 11. Histograms of amplitudes (upper row) and normalised phase angles (bottom row) resulting from the FFT analysis for the THINGS sample (in blue) and for the BS dispersion maps σbs (in red).

5.2.2. Uncertainty on the beam shape

The observed dirty beam of interferometric observations is complex. In this section, we measure σ_bs using the observed dirty beam for 12 of the galaxies from THINGS. Maps of the dirty beam are mainly composed by a combination of two Gaussian functions, a narrow one to describe the PSF core, similar to the one of our model, and a second one with larger full width at half maximum to describe the PSF wings. The extent of dirty beam maps is often larger than those of the galaxies but the beam power at large angular distance from its centre is very low. We apply a circular mask to the dirty beam map to deal mostly with positive PSF values, and get rid of small variations at large radius from the beam centre that could induce spurious effects in our analysis. We vary the radius of the masks from 5 to 50 pixels, in order to be consistent with our 2D Gaussian model which has an extent of 5 to 10 pixels within the THINGS sample. Using a dirty beam map truncated at the same size and scale as our 2D Gaussian model, the resulting σ_bs maps are very comparable to those obtained with the Gaussian beam. For example, within a cut of five pixels around the beam centroid, the amplitude of σ_bs is ∼10% lower on average than for the Gaussian shape throughout the THINGS sample. However, when the spatial extent of the dirty beam is larger, the patterns in σ_bs are more prominent, leading to larger amplitudes. For instance, with a mask cut at a radius of 20 pixels, the amplitudes are twice larger than our model, on average. As for the phase angles, the minor difference between the two beam shapes is that the peaks are narrower and slightly larger for k = 2 and 4 in the case dirty beam.

5.2.3. Kinematics at higher angular resolution

Daigle et al. (2006) and Dicaire et al. (2008) published Hα data cubes for respectively 28 and 37 galaxies from the SINGS sample (Kennicutt et al. 2003). Those observations were led with a Fabry-Perot interferometer at the 1.6 m telescope of the Observatoire du Mont Megantic (OMM, Canada), and the 3.6 m telescopes at the Canada-France-Hawaii Telescope and the European Southern Observatory (ESO, La Silla, Chile). In our sample, all galaxies were observed at OMM except NGC 3621 and NGC 7793, observed at ESO. In this section, we benefit from the higher angular resolution of the Hα velocity fields to estimate, and correct from, the BS effects on H I THINGS velocity fields. The angular resolution of the Hα data is seeing limited, ∼1″ at ESO and ≲3″ at OMM, which represents a gain of ∼2 to 13 with respect to THINGS data. The underlying hypothesis is that the Hα kinematics have a similar behaviour as the H I one. Nevertheless, the flux distribution of the ionised gas is more peaked due to the star forming regions. To model the Hα BS map, we used the Hα velocity field weighted by the flux distribution of the neutral gas. Since the kinematical maps available from Daigle et al. (2006) and Dicaire et al. (2008) were obtained using adaptative smoothing, we compute new maps from the sky-removed data cubes. We first apply a 2-pixels Gaussian smoothing to increase the Hα S/N and then extract kinematical maps using the same barycentre method used in the two original papers. We also reproject the THINGS flux density maps on the same grid as Hα maps in order to match the positions and pixel size of both datasets. To compute the BS maps, we apply the same methodology as defined in Sect. 3.1, using the reprojected H I flux distributions, the 2-pixels Gaussian smoothed Hα velocity field and the single 2D Gaussian PSF to describe the new H I angular resolution.

In Fig. 12, we show the resulting Hασ_bs map (left panel), compared to the H Iσ_bs map (middle panel), as well as the ratio of the two maps, for an example galaxy (NGC 3521). The BS pattern is similar in both cases, with the typical cross-shape in the centre. While the BS correction inferred from H I data is less peaked in the centre than for Hα, the latter also decreases more slowly with radius. This explains the large values in the normalised ratio found far from the centre. Beam-smearing amplitudes are higher when modelling with the ionised gas. As the amplitude may be an important concern in our modelling, we study galaxy by galaxy the amplitude ratio for the two different cases.

Fig. 12. Comparison of the BS models using optical and cm data for the galaxy NGC 3521. Left panel: Hα data based BS model. Middle panel: H I data based BS model, trimmed at the same extent as the optical model. Right panel: ratio of the BS models, defined as BS(Hα)/BS(H I).

In Fig. 13, a comparison of the BS contributions is made for the THINGS sample. For each galaxy, a 2D histogram with the amplitude of the H I BS correction on the x-axis and the ratio of the BS corrections in Hα to H I on the y-axis, is made. Then, the histograms were normalised by the number of pixels in the maps, to avoid being biased by objects with a larger number of pixels. Figure 13 is then obtained by summing all the 2D histograms. It shows that the ratio is > 1, which is not surprising given that the Hα BS benefits from a higher angular resolution. On average, the H I BS is 2.3 times smaller than the Hα BS. We also observe a trend at low H I BS with high ratio induced either by noise or by discontinuities between regions in the Hα data. Indeed, a continuously increasing ratio with decreasing H I BS can be modelled by assuming that the Hα BS is equal to the quadratic sum of the H I BS and a constant term of ∼3 to 10 km s⁻¹. This constant term can be related to pixel-to-pixel variations either induced by noise in the Hα data, or by true local velocity variations or discontinuities unseen in H I. Amplitudes vary within the sample, due to various intrinsic velocity gradients in each galaxy, as well as its distance and geometry.

Fig. 13. Pixel-by-pixel comparison of the BS contribution using Hα SINGS and H I THINGS data. The y-axis shows the ratio of the BS amplitudes Hα/H I and in the x-axis is the amplitude of the H I model.

Despite these findings, we decided to continue working with the H Iσ_bs for several reasons. First, the spatial extent of H I gas is larger, making it possible to evaluate the local BS contribution of spiral arms or other features that are not located in the inner disc regions. A second reason is that the Hα gas does not trace the same interstellar component, therefore the geometry of the discs might differ, and the flux distribution are noticeably different. Furthermore, it has been shown that the imprint of the systematic phase along the minor axis of the galaxy M 33 did not depend on the angular resolution (Chemin et al. 2020), while a range of 7 in the size of the synthesised beam was probed. Using H I maps to derive σ_bs should thus not be an issue in the analysis (see also Sect. 5.2.4).

To assess the impact of BS under more conservative conditions, we then subtracted respectively twice and thrice the H I BS from the THINGS velocity dispersions. The result of this process is that a large number of rings are removed in this case. Indeed, in order to avoid sampling issues in the FFT, we removed rings having more than 5% of pixels with undefined value after the BS correction, due to this quadratic correction being larger than the actual measured dispersion. The median fraction of such pixels per ring is 0.7%, 4.0% and 16.6% when subtracting once, twice and thrice the BS contribution, respectively. Only 621 and 258 tilted rings remain in the two latter cases, respectively, with less than 5% of pixels with an undefined value, decreasing drastically the statistics. With these particularly conservative BS subtractions, some galaxies totally disappear from the remaining rings, which indicates that the velocity dispersion would be totally dominated by BS, which should not be the case owing to the (high) H I resolution.

In Fig. 14, we show the phase angles for k = 2 and k = 4 resulting from the subtraction of once or twice the contribution σ_bs (red and blue histograms, respectively). We note that the amplitude distributions are not strongly affected. As mentioned before, more than 200 rings were removed, which explains why the average number is lower for the red histograms. For k = 2, the peak seen in the blue distribution is no more present in the red one. This behaviour is also observed for k = 1 (not shown here). For k = 4, in both cases we observe the strong peak at ϕ₄/T₄ ∼ 0.5. The interpretation is not straightforward, as the observed phase angle can be explained either by a residual BS signature or by the correlation of the asymmetry with the BS. Nevertheless, the fact that a strong peak remains in ϕ₄ tends to indicate that BS alone cannot explain the observed trends.

Fig. 14. Comparison of the normalised phase angles histograms of orders k = 2 and k = 4 (left and right panels, respectively) in the velocity dispersion by subtracting once and twice the contribution from the BS effect (blue and red colours, respectively).

We also smoothed the Hα data at the average resolution of the H I, 11″ for all the sample with the goal to verify that the same level of smoothing of the gradients implies a similar BS effect. The results are as expected, being that the normalised residuals for each galaxy show a distribution peaked on 0, indicating that both modelling are consistent with each other.

Similarly, to verify the consistency of the findings described above, we also performed the Fourier analysis of an example galaxy (NGC 2841) using the 6″ angular resolution THINGS velocity and density maps, as obtained from robust weightings. The amplitude of the BS is larger than for natural weightings, which is due to the combined effects from the higher angular resolution with the noisier velocity field for robust weightings than for natural weightings. This finding matches perfectly the result found with the Hα kinematics. Nevertheless, no significant difference in the strength and phase angles of the Fourier modes is observed with respect to the natural weighting for this galaxy.

5.2.4. Amplitude of BS residuals

The correlation between BS and observed signatures found in Sect. 5.2.1 may indicate that BS residuals remain after BS correction, due to the use of observed flux maps and velocity fields rather than high resolution data to infer the correction (see Sect. 3.1). To quantify residuals, we use the toy models described in Sect. 3.3 and Appendix B. We focus on the isotropic case with no perturbation and a steep inner gradient of the rotation curve, for which the impact of BS is the largest, and added realistic noise to dispersion maps (cf. Appendix A). Correcting for BS as described in Sect. 3.1 reduces the strengths of orders k = 2 and k = 4 by a factor between 2 and 4 and is more efficient for less inclined galaxies and for order k = 2. We also show that the maximum strength of orders k = 2 and k = 4 of BS residuals is reached at a lower radius than without BS correction, and that it decreases for decreasing inclinations. For inclinations below 60°, residual BS cannot account for the k = 2 signature at radii greater than a few beam FWHM and at radii > 10 beam FWHM for k = 4. At 75°; however, the BS signature remains well above the uncertainties because (i) the projected velocity gets higher than at lower inclination, and (ii) the spatial resolution along the minor axis gets poor, inducing more hidden velocity gradients within the beam and pixel size. We also performed an FFT analysis separately for the 50% largest and smallest radii and found, as expected, that the amplitude of orders k = 2 and k = 4 is reduced (by a factor larger than two) and that the significance of the peak in phase histograms of the fourth order is much lower (by a factor around 4) at large radii. For the model with a weaker velocity gradient, when noise is added, all the systematic phase angles disappear and the distributions become flat. The impact of BS induced by large-scale rotation depends on the physical resolution, the rotation curve shape, the projected maximum velocity, galaxy inclination and high resolution flux distribution. A complete study of BS on velocity dispersion maps depending on all these parameters is beyond the scope of the present study.

About half galaxies in our sample have projected maximum rotation velocities above 150 km s⁻¹, as inferred from the profile width at 20% of the peak intensity given in Walter et al. (2008), and about 60% have an inclination larger than 60°, which may therefore present BS residuals signatures in the fourth order. We first split the sample between inner and outer rings. The FFT amplitudes of the innermost rings are 27% to 43% larger than for the outermost rings, depending on the order. The phase angle distributions are shown in Fig. 15, with the innermost and outermost rings shown in blue and red, respectively. These distributions are different. For k = 2, the distribution for the outermost rings is similar to that of the global distribution whereas the distribution does not peak on 0.5 ϕ₂/T₂ for the inner rings, as would be expect if the dispersion maps were dominated by BS residuals. The inner and outer regions both show a prominent peak k = 4 close to the half period, the peak being slightly offset towards lower phase angles in the outer parts. We also split the sample in low versus high projected maximum velocities, leading to the same sub-samples as for the absolute magnitude (see Sect. 4.4 and Fig. 8). As already discussed, the peak observed in the phase angle histograms of the fourth order is much more pronounced when projected velocities are high. This may be a proof that BS residuals are affecting our results at least for the fourth order. Nevertheless, given the large values of order four strengths σ₄ compared to what is expected from residuals in toy models, it may also be that such a signature is intrinsic to massive galaxies with high gas density contrasts. On the other hand, the peak in the second order around half period of galaxies with low projected velocities cannot be attributed to BS.

Fig. 15. Histograms of normalised phase angles resulting from the FFT analysis for the THINGS sample splitting the rings as a function of radius (inner part in blue, outer part in red).

5.2.5. A possible link between velocity gradients and elevated velocity dispersion

We have shown that BS may explain the signatures we observe if not properly taken into account, as there are strong correlations between the modelled BS pattern and the observed velocity dispersion. Nevertheless, we have also shown that the amplitude of BS induced by large-scale motions is expected to be too low to have a significant impact on the data. The fact that BS seems to be coupled to observed signatures in the phase angles of asymmetries may indicate that the regions with the strongest velocity gradients or discontinuities, especially in the outer parts of galaxies, are also regions with an intrinsic large velocity dispersion at a higher level than what would be induced by BS from the observed gradients. These regions with large velocity gradients are often related to perturbations like bars, spiral arms, inter-arms and warps, especially at large radius, which means that there might be locally induced turbulence responsible for perturbations in both the velocity field and the velocity dispersion. This may also be that BS occurs on much smaller scale than that reached within our observations and that what we observe is not related to large-scale motions. A complementary way to investigate the effects of unresolved gradients of asymmetric velocities on the velocity dispersion is addressed in Sect. 5.3.1.

5.3. Projection effects in random motions

Finding systematic phase angles of asymmetries in velocity dispersion maps is a surprising result. The projection of asymmetric patterns along the line-of-sight should be randomly distributed. Flat distributions of phase angles were then expected in Fig. 5 from a theoretical perspective (see discussion in Appendix C). Privileged phases thus imply a dependency of the patterns on the orientation of the disc principal axes with respect to the observer, because the major axis of discs has been chosen at the origin of the azimuthal angles (θ = 0 along the major axis). In other words, this suggests the presence of projection effects in the velocity dispersion maps. In this section, we thus want to study the possible origin(s) of projection effects in the velocity dispersion. This is achieved by investigating the effect of asymmetries in the ordered motions on the random motions (Sect. 5.3.1). We also address the effect of correlated velocity components (Sect. 5.3.2).

5.3.1. Effects from the asymmetric ordered motions

It is usual to assume that the velocity ellipsoid of gas is isotropic in galactic discs. Under this assumption, no deprojection of data is required, so that the observed velocity dispersion is a direct proxy of the gaseous random motions. However, this assumption is very idealised, and can only apply to axisymmetric kinematics, which rarely occurs in reality. Indeed, the velocity fields (V_los) always exhibit disturbances, as caused by the projection of asymmetric ordered motions (V_R, V_θ, V_z) due to, for example, bars, rings, spiral arms, or warps. The effect of any of such disturbances in V_los should also propagate to σ_los because of the implied gradients in one or several directions in the plane. We test this hypothesis by studying the toy models with asymmetric kinematics perturbations described in Sect. 3.3 and Appendix B in both isotropic and anisotropic cases. We mainly focus hereafter on cases with the weakest velocity gradient (shown in Fig. B.2) to reduce possible residual BS signatures related to unresolved rotation (see Sect. 5.2.4) and emphasise on local perturbations.

In the axisymmetric and isotropic model (used as reference), the strength of the Fourier coefficients of the various orders are comparable (∼0.5 km s⁻¹ on average). In all other models, σ₁ and σ₃ barely vary from this value, while the k = 2 and k = 4 components are significantly larger, being mostly > 2 km s⁻¹, except in the axisymmetric case with uniform local anisotropy, where σ₄ is very comparable to the odd amplitudes. The distributions of the phase angles for the negligible modes k = 1 and k = 3 do not show systematic peaks. The even orders thus always dominate the dispersion asymmetries, and the bisymmetry is the strongest perturbation, on average.

Any model with a velocity ellipsoid showing a radial bias σ_θ < σ_R creates a bisymmetry theoretically aligned with the minor axis. The strength of k = 2 asymmetries in σ_los within mock data obtained from models combining effects from the uniform radial bias and V_R and/or V_θ streamings is between 20 and 40% larger than those without the radial bias, on average, while the amplitude of the mode k = 4 is not affected. The distribution of ϕ₂ is sharply peaked at π/2 for the axisymmetric case of anisotropy with the radial bias, and the peak is slightly enlarged when there are asymmetries in V_R and/or V_θ, regardless of their strength. Moreover, in this latter case with the velocity perturbation exclusively set along the azimuthal direction (ΔV_R = 0, ΔV_θ ≠ 0), the incidence of ϕ₂ = 0 or 1 remains extremely low. These findings can be explained with the help of Eq. (3) from Chemin et al. (2020), where having k = 2 asymmetries in σ_los aligned with the minor axis is a genuine imprint of a radially biased velocity ellipsoid (see also Eq. (C.11), but with null cross terms). Thus, a local and uniform anisotropy of the velocity ellipsoid with a radial bias drives the phases of asymmetries in σ_los more strongly than the effect of the velocity anistotropy induced at larger scale by the projection of V_R and V_θ gradients. This is because the radial bias is applied to all velocities in the mock data, while the V_R and V_θ streamings are only restricted to the locations of the perturbations. If a tangential bias σ_θ > σ_R had been assumed, the bisymmetry would have been exclusively aligned along the disc major axis, to the detriment of the minor axis. Because THINGS and WHISP galaxies do not show a ϕ₂ distribution strongly concentrated around π/2, anisotropic gaseous velocity ellipsoids consistent with a radial bias uniformly distributed everywhere in discs, like the one presented here, seem very unlikely. A solution to avoid the ϕ₂/T₂ concentration at 0.5 is to have either an axis ratio σ_θ/σ_R which varies with radius, or tilted velocity ellipsoids. Section 5.3.2 provides a hint of how this goal could be achieved.

The impact of velocity anisotropy arising from asymmetries in V_R (or V_θ, respectively) is to greatly decrease the likelihoods of finding ϕ₂ close to zero (π/2), in the case the velocity perturbation is along the radial axis, that is ΔV_θ = 0 and ΔV_R ≠ 0 (azimuthal axis, ΔV_R = 0, ΔV_θ ≠ 0). This is because there is little projection of V_θ (V_R, respectively) onto the minor (major) axis. These results remain valid regardless of the uniform isotropy/radial bias assumptions. If there were an equal number of tilted rings dominated by radial and azimuthal velocity streamings, this would lead to distributions of ϕ₂ with minima around 0 and π/2, and conversely with maxima ϕ₂ = 0.25π and 0.75π, which values are not far from those observed in the sample within the quoted uncertainties.

Then, still for models with asymmetries in V_R and/or V_θ, the distributions of σ₄ and ϕ₄ are very similar in both cases of uniform isotropy or radial bias. Thus, the uniform velocity anisotropy has negligible effect on the fourth order perturbation of σ_los, and this latter is an harmonics of the second order perturbation in our toy models. The information about the orientation of the bisymmetry (i.e. along the major or minor axis) is nevertheless lost modulo π/2 in the fourth order⁵. Moreover, it is obvious that asymmetries have destroyed the systematic peak at ϕ₄ = 0.5 T₄ seen in the reference axisymmetric models. This is because this peak could be caused by small residuals from BS correction, as mentioned previously. Unless our velocity asymmetries in the toy models are significantly overestimated with respect to the reality, the occurrence of a systematic peak at ϕ₄ = 0.5 T₄ in the observations indicates that the fourth order asymmetry is very likely an harmonic of a second order which would have preferential phase angles ϕ₂ close to either 0.25π and 0.75π, which values remains close to the observed ϕ₂ ∼ 0.35π and 0.65π within the quoted uncertainty of 0.1 T₂.

We note that combining a steep velocity gradient top the perturbations (see Fig. B.3) does not change qualitatively the previous results for the order k = 2 though the contrast between high and low probability is a bit reduced, especially in the isotropic case with the velocity perturbation along V_R, due to residual BS inducing a larger probability for null phase angles whatever the direction of the perturbation. Anisotropy always dominates order k = 2 when present. For order k = 4, a peak related to residual BS clearly remains around ϕ₄/T₄ = 0.5, which increases the trend already observed with weak rotation-induced BS.

Finally, although not shown here, we observe that when the BS correction is made with flux density and velocity maps at the highest angular resolution of the mock data (native resolution of ∼50 pc, before spatial smoothing), no significant residuals due to BS remain in the distribution of phase angles. We note here that this value of 50 pc is only applicable to our toy models, and not to real observations, because we do not know precisely the angular scale of velocity streamings in galaxies. This indicates that the velocity anisotropy from the streamings can only affect the phase angles when the streamings are not properly resolved. Therefore, our method consisting in searching for systematic effects in the phase angles of σ_los from observations of ‘limited’ resolutions like THINGS is a powerful tool to predict the directions in the plane towards which the perturbations and anisotropy are more significant. In particular, it implies that the lopsidedness perturbation of H I kinematics in THINGS is more likely to be stronger along the tangential dimension than along the radial direction, hence a preferentially large-scale tangential bias in the velocity anisotropy, in order to explain the higher incidence of ϕ₁ along the major axis of the galaxies. Similarly, it is less probable to find a bisymmetry perturbation stronger in the tangential direction, to explain the lower incidence of ϕ₂ towards the major axis of the galaxies. The observations of ϕ₂ aligned near the disc minor axis in lower mass galaxies (Sect. 4.4), similarly to the outer part of M 33, would thus likely be explained by dominant radial streamings, hence a stronger large-scale radial bias in the velocity anisotropy (Chemin et al. 2020). Another consequence of the analysis is that velocity asymmetries in the first moment maps of galaxies must be observed on a scale of about the PSF width, though of reduced amplitude. We thus expect that systematic effects are present in the phase angles of V_los asymmetries, by analogy with those evidenced in σ_los, tracing the direction of dominant bias of the velocity anisotropy from the perturbations. Comparing the properties of asymmetries of V_los and σ_los is beyond the scope of this article, however.

5.3.2. Effects from anisotropic velocity ellipsoids with tilted axes

Until now, the tests on the shape of the velocity ellipsoid has only considered components of the velocity tensor which are independent from each other. However, the asymmetries in V_R, V_θ and V_z all result from the same perturbed gravitational potential. By consequence, the components must show a certain degree of correlation, at least at the vicinity of perturbations. Correlated velocities are evidenced in the kinematics of stellar populations of the Milky Way, through tilted axes of the stellar velocity ellipsoid. The observed tilt angle of the axes of the radial and vertical velocity distribution of Galactic halo stars can inform on the 3D shape of the mass distribution and the dark matter halo (Siebert et al. 2008; Smith et al. 2009; Wegg et al. 2019), while that of the radial and tangential velocity distribution (the vertex deviation) is caused by a bisymmetric perturbation in the disc, such as a bar, a spiral structure, a warp, or more simply an elliptical disc (Kuijken & Tremaine 1991; Debattista et al. 2019).

We can test the effect of correlations between H I velocities by means of the uniform velocity anisotropy presented above, by considering the possibility of tilted ellipsoids, as we do for stars. This has a consequence on the projection of the ellipsoid on the line-of-sight, and following the complete mathematical description from Appendix C, the squared velocity dispersion is given by:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\sigma }_{\mathrm{los}}^2=& (\frac{{\sigma }_{\theta }^2-{\sigma }_{\mathrm{R}}^2}{2}cos2\theta +{\sigma }_{\mathrm{R}\theta }^{2}sin2\theta ){sin}^{2}i\hfill \\ \hfill & +({\sigma }_{\mathrm{R}\mathrm{z}}^{2}cos\theta +{\sigma }_{\theta z}^{2}sin\theta )sin2i\hfill \\ \hfill & +{\sigma }_z^2{cos}^{2}i+\frac{{\sigma }_{\theta }^2+{\sigma }_{\mathrm{R}}^2}{2}{sin}^{2}i,\hfill \end{array}\end{array}$$(6)

where ${\sigma }_{\text{R}\theta }^2$, ${\sigma }_{\theta z}^2$ and ${\sigma }_{\text{Rz}}^2$ are the covariance terms marking the correlation between the components of the velocity tensor in the cylindrical frame of a disc. This equation can then be expressed as a sum of three terms of orders 0, 1 and 2, with an analogue formalism as the one we used to study asymmetries:

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{los}}^2=k_0+k_{1}cos(\theta -{\alpha }_1)+k_{2}cos(2(\theta -{\alpha }_2)),\end{array}$$(7)

where k_i are the amplitudes of these terms and α_i are their phases, that all depend on the various velocity dispersion terms. More specifically, we have

$$\begin{array}{c}\hfill tan2{\alpha }_2=2{\sigma }_{\mathrm{R}\theta }^2/({\sigma }_{\theta }^2-{\sigma }_{\mathrm{R}}^2).\end{array}$$(8)

Assuming no tilted axes for the uniform velocity anisotropy in the previous section left no choice than having a dispersion bisymmetry systematically aligned with the minor or major axis (α₂ = π/2 or 0, respectively), depending on the sign of ${\sigma }_{\theta }^2-{\sigma }_{\text{R}}^2$. This is at odds with the observation that the second order phase angles ϕ₂ = 0.35T₂ and ϕ₂ = 0.65T₂ are particularly more likely. Equations (7) and (8) thus imply that it is possible to get various ϕ₂ orientations if the velocity anisotropy is not as uniform through the disc, and more importantly if the covariance varies. In the Milky Way disc, the measurement of the axis ratio of the tangential-to-radial velocity dispersions of young stars as measured with Gaia spectro-astrometric kinematics, is σ_θ/σ_R = 0.66 (Gaia Collaboration 2023). Since young stars kinematics are expected to be comparable to that of the gas, we further assume for this exercise that the velocity anisotropy of the planar components of the interstellar gas for any galaxies in our sample is comparable to the value found for young stars in the MW disc. Reminding now that ${\sigma }_{\text{R}\theta }^2={\rho }_{\text{R}\theta }{\sigma }_{\text{R}}{\sigma }_{\theta }$, with ρ_Rθ the correlation coefficient, an appropriate choice of covariance and correlation coefficient which explain the higher incidence of ϕ₂ ≡ α₂ = 0.35π and 0.65π for gas is ρ_Rθ = ±0.6. Furthermore, Eq. (7) implies that correlations between the vertical and the two planar components should also be considered to get systematic phases of first order phases. In order to explain those seen in Fig. 5, we would thus need to have α₁ ∼ 0 or π to match the observations, that is, $tan{\alpha }_1={\sigma }_{\theta z}^2/{\sigma }_{Rz}^2~0$, following Eq. (C.13). This is possible if ρ_θzσ_θ ≪ ρ_Rzσ_R. Here again, studies of the kinematics of young stars in the Galaxy should help to test these predictions, which is relevant only if the shape of gaseous velocity ellipsoids are comparable to those of young stars in nearby discs. It would then be interesting to investigate how these predictions on ρ_Rθ compare with the observed planar correlations of young stars of the Galactic disc, and how ρ_θzσ_θ compare with ρ_Rzσ_R for the same stellar populations. To our knowledge, none of such distributions have ever been measured yet for the Milky Way.

6. Summary and concluding remarks

For this work, we have performed the first systematic search and characterisation of asymmetries in velocity dispersion maps of H I gas in nearby galaxies. We used the best H I data available in the archives, namely 32 galaxies from THINGS and WHISP samples with a spatial resolution better than 1.5 kpc and a neutral gas distribution extended by at least seven times the resolution, to carry out this pioneering work. We performed a Fourier analysis of the observed second moment maps previously corrected from both the line spread function and the two-dimensional BS contribution tracing the velocity gradients occurring at a large angular scale. This allowed us to measure the strength and phase angle for each of the first four Fourier harmonics, which are the most important perturbations in the random motions of H I gas.

We find a wide range of strength of asymmetries, up to ∼35 km s⁻¹. The shape of the likelihood distributions of the amplitude of the Fourier modes are similar, with a peak around 5 − 10 km s⁻¹, and an extended tail towards larger values. Overall, the strongest Fourier mode is the bisymmetry (k = 2), with a mean amplitude of ∼11 km s⁻¹, followed by the lopsidedness k = 1 asymmetry. On the other hand, the shape of the likelihood distributions of phase angles strongly depends on the Fourier mode, and the distributions present significant variations. For the k = 1 asymmetry, the probability to find a phase angle near the major axis of the discs is larger, while finding an asymmetry orientated near 0.15 and 0.75 × 2π is less probable. For k = 2, the likelihood is larger at phase angles of ∼0.35 and 0.65 × π, and lower near the disc major axis, whereas that of the k = 4 mode is concentrated at a value of 0.5 × π/2, implying an orientation mainly lying at an equivalent separation from the disc major and minor axes. These systematic phase angles are robust against the number of titled rings and against the significance of the detected signal. Using a uniform number of rings for all galaxies lowers the significance of the result, but the trends for the shape of the probability distributions are preserved. More uniform likelihood distributions were expected for randomly distributed perturbations in the galaxies and with isotropic velocity ellipsoids for the H I gas. This is evidence that strong projection effects imprint on the velocity dispersion maps of H I gas in nearby galaxies.

Finding the origin(s) of the asymmetries and its systematics is a difficult task because of the nature of the velocity dispersion itself, which stems from instrumental effects and various processes inherent to the galactic random motions (dynamical, hydrodynamical, and thermal). In that aspect, we have shown that large-scale-induced BS might produce non-negligible residuals, especially in the fourth order for highly inclined galaxies with high projected rotation velocities. Typical signatures of BS on phase angles are orientated near the major axis (k = 1), the minor axis (k = 2), and at 0.5 × π/2 (k = 4). Nevertheless, these orientations are not fully consistent with the observations, especially for the second order, and the strength of Fourier modes in the BS contribution is much smaller than the observed asymmetries. This suggests that unresolved velocity streamings could be responsible for some variations in the likelihood distributions of the phase angles. The correlations found between Fourier decomposition of H I gas density and velocity dispersion maps suggest that the observed bisymmetry, and thus streaming motions, are induced by spiral arms, bars, or warps.

We thus further addressed the impact on the random motions from such velocity streamings present in the ordered motions V_R and V_θ. To do so, toy models mocking discs with a barely resolved spiral perturbation in both the density and kinematics were made, from which mock moment maps of discs seen under various projection angles, strengths, and orientations, of the perturbations were derived, and Fourier analyses were performed, by analogy to the observations. When streamings in the ordered motions are stronger in the radial dimension, hence in the presence of a radially biased velocity anisotropy on large scales, they propagate to the random motions through the BS effect, but on small scales, and it is more likely to find phase angles of the k = 2 mode orientated along the minor axis. The same effect applies to the major axis in the presence of a tangential large-scale anisotropy as caused by stronger streamings in the tangential dimension. When both V_R and V_θ streaming motions are present, the phase angles of the bisymmetry of velocity dispersion reflect the competition between the two directions. Therefore, producing distributions of k = 2 phase angles similar to the observations would probably require specific mixing between the radial and tangential perturbations in discs, such as having different regions of dominance of streamings in one direction over the other. Also, no k = 4 phase angles centred around π/4, as in the observations, could be produced within the toy models. Our results suggest that the fourth order asymmetry is likely a combination between BS residuals and harmonics of the second order asymmetry in the velocity dispersion.

This analysis thus shows that it is possible to constrain which asymmetry between V_R or V_ϕ dominates the streamings by studying the systematic orientations of the k = 1 and k = 2 phase angles of σ_los. In particular, to explain the systematic orientation of k = 1 phase angles along the major axis observed in the samples, this work suggests that lopsidedness in V_los must be dominated by tangential perturbations, while the systematic orientations of k = 2 phase angles outside the major axis require elliptical streamings in V_los affected by a stronger contribution from radial perturbations. By consequence, the observation of velocity asymmetries more significant in one direction in the plane than the other implies that H I velocities are highly anisotropic, at least at the angular scales probed by the Fourier analysis and the resolution of the observations. Yet it is possible that the velocity ellipsoid of gas remains isotropic at angular scales much smaller than those probed by the data, thus very locally, perhaps at the scale of individual clouds. Moreover, we predict that the signature of velocity anisotropy on large scales is also present in velocity fields of galaxies, again in the shape of systematic orientations of asymmetries near the principal axes of the discs. This prediction should greatly benefit from analyses of residual velocity fields of galaxies.

We further show that an anisotropic velocity ellipsoid with correlations between velocity components can generate asymmetries in the velocity dispersion with any orientation. Within this formalism, we predict that the probability to find coefficients of a correlation between V_R and V_θ must be larger towards |ρ_Rθ|∼0.6, in order to explain the larger occurrence of gaseous ϕ₂ at ∼0.35π and ∼0.65π, assuming σ_θ/σ_R ∼ 0.7 as for the youngest populations of stars in the Milky Way. It would mean that the correlations between velocities in nearby gaseous discs are not random values, and thus that perturbations of the gravitational potential share common properties among the galaxies. This prediction could be compared with the velocity correlations measured for the youngest stellar populations of the Milky Way, to which the kinematics of interstellar gas should be comparable in many aspects.

This pioneering work has highlighted the importance of velocity anisotropy in shaping the asymmetries of H I velocity dispersions, through velocity streamings in the H I gas which are dominant in particular directions in the plane and/or tilted velocity ellipsoids. The present study is intended to serve as a stepping stone towards the design and the analysis of future large-scale surveys of galaxies, observed at higher sensitivity and resolution than those samples studied here, such as those upcoming with the Square Kilometre Array and its precursors for the H I gas, or those already observed in the molecular interstellar medium of galaxies by ALMA, as well as in ionised gas and stellar velocity fields. In that aspect, comparisons of the shapes of asymmetries in σ_los from both the stellar and gaseous components will be crucial to support the proposal that gas velocities are anisotropic. Harmonic decompositions of V_los will also be helpful. This will be the subject of future papers following this work.

Acknowledgments

We are grateful to F. Walter for sharing the maps of the dirty beams of THINGS observations. P. Adamczyk and L. Chemin warmly acknowledge the Comité Mixto European Southern Observatory-Gobierno de Chile. The research of L. Chemin is funded by the Chilean Agencia Nacional de Investigación y Desarrollo (ANID) through the Fondo Nacional de Desarrollo Científico y Tecnológico (FONDECYT) Regular Project 1210992.

References

Appendix A: Significance of FFT results with noisy data

In inner regions of galaxies, rings are smaller and contain less pixels than in the outermost parts. Moreover, the average velocity dispersion tends to decrease with radius, and can thus vary significantly from one ring to another, so that the properties of rings are not necessarily uniform as a function of radius. The results of a FFT, and in particular the precision that can be obtained on Fourier amplitudes, is impacted these parameters. It is essential to understand how the noise in the data affects the results of the FFT and to estimate to which limits we can recover the real signal, or in other words the reliability of the recovered signal. In our study, we apply the FFT implemented in the Python package FFTPACK from the Scipy library on the squared velocity dispersion σ². To have a realistic estimate of the noise σ², we express the uncertainties on σ² using the differential of the squared velocity dispersion as:

$$\begin{array}{c}\hfill d({\sigma }^2)=2{\sigma }_{0}d\sigma ,\end{array}$$(A.1)

where ${\sigma }_0^2=\langle {\sigma }^2\rangle$ is the value of the zeroth order of the FFT, because all higher orders cancel out on average over a ring, and where dσ is the uncertainty on the observed velocity dispersion and that roughly corresponds to the uncertainty on the instrumental resolution at a given signal-to-noise ratio. Since we do not have uncertainty maps due to the processing techniques (moments method) on THINGS data, dσ is not known precisely. Since the spectral resolution of both THINGS and WHISP data corresponds to about σ_LSF ∼ 2 − 5 km s⁻¹, and because the actual uncertainty on the line width depends on the S/N, we assume that the uncertainty on the velocity dispersion is of the order dσ ∼ 1 km s⁻¹.

The way the uncertainties propagate from the data to the FFT coefficients depends on how the FFT is implemented. A discrete Fourier transform is used in FFTPACK to compute Fourier coefficients:

$$\begin{array}{c}\hfill C_k={\displaystyle \sum _{m=0}^{N-1}{c}_{m}cos(-\frac{2\pi km}{N})+j\sum _{m=0}^{N-1}{c}_{m}sin(-\frac{2\pi km}{N}),}\end{array}$$(A.2)

where j is the unit imaginary number, C_k is the Fourier coefficient of order k as computed by FFTPACK, and c_m are the N elements in the discrete array. In our study, c_m are the squared velocity dispersion values along rings, σ²(m), and we deduce the amplitude of each order k ≠ 0 as ${\sigma }_k^2=2\left|C_k\right|/N$, and ${\sigma }_0^2=\left|C_0\right|/N$ as the average of the squared velocity dispersion values along the ring. From this formalism, it can be shown that for k ≠ 0:

$$\begin{array}{c}\hfill \begin{array}{cc}& \mathrm{Var}[{\sigma }_k^2]=\frac{4}{N^2}{\displaystyle \sum _{m=0}^{N-1}\mathrm{Var}[\sigma {(m)}^2]}\hfill \\ \hfill & +\frac{8}{N^2}{\displaystyle \sum _{\mathrm{\Delta }m=1}^{N-1}\sum _{m=0}^{N-1-\mathrm{\Delta }m}\mathrm{Cov}[\sigma {(m)}^2,\sigma {(m+\mathrm{\Delta }m)}^2]cos\left(\frac{2\pi k\mathrm{\Delta }m}{N}\right),}\hfill \end{array}\end{array}$$(A.3)

with Var and Cov being respectively the variance and covariance. Assuming that the variance is constant across the ring (Var[σ(m)²]=Var[σ²]) and that covariance is null between different pixels, we can get at first approximation that:

$$\begin{array}{c}\hfill \mathrm{Var}[{\sigma }_k^2]=\frac{4}{N}\mathrm{Var}[{\sigma }^2],\end{array}$$(A.4)

leading to the uncertainty on coefficients:

$$\begin{array}{c}\hfill d({\sigma }_k^2)=\frac{4{\sigma }_{0}d\sigma }{\sqrt{N}}·\end{array}$$(A.5)

Interestingly, this equation shows that the uncertainty is the same for all orders, that it is proportional to the zeroth order, and that it decreases with radius since the number of pixels increases linearly with radius and because σ₀ rarely tends to decrease with radius.

We have checked numerically the consistency between this expression and the actual results from FFT using a toy model. We define the noise δ on ${\sigma }_{\mathit{obs}}^2$ as a Gaussian distribution centred on 0 and with a standard deviation of σ_Gauss = d(σ²). We create a toy model of the observed velocity dispersion with a constant velocity dispersion to which we added noise:

$$\begin{array}{c}\hfill {\sigma }_{\mathit{obs}}^2={\sigma }_0^2+\delta ,\end{array}$$(A.6)

with δ a random value in the noise distribution. Doing so, the distribution of each observed point is taken independently. We used several values for dσ, σ₀, and N the number of points in the FFT, since N varies from one ring to another and from one galaxy to another. We do 1000 iterations of the FFT and analyse the amplitudes for orders 1 to 4, which are null in our model, as shown in Eq. A.6. Fig. A.1 shows the measured uncertainties of the toy models and confirms that Eq. A.5 is a good description of the uncertainties. Equation A.5 can therefore be used to infer at first order the significance of Fourier coefficients, $s_k={\sigma }_k^2/d({\sigma }_k^2)$, depending on the number of rings and on the mean velocity dispersion value.

Fig. A.1. Square root of the noise on the squared velocity dispersion amplitude for orders k = 1 to 4 of the FFT as a function of σ0 for two distinct values of dσ and various number of points N in the FFT. The curves correspond to Eq. A.5 with dσ = 1 km s−1 (dashed) and dσ = 4 km s−1 (dotted), each corresponding to a different value of N (from top to bottom: 1000, 2000, 4000 and 8000). Coloured symbols correspond to different orders measured on toy models. The similar behaviour of orders leads to superposed points.

Appendix B: Toy model description to study projection and resolution effects on random motions

Here, we describe the detailed generation of the toy models presented in Sect. 3.3 and used in Sections 5.2.1 and 5.3.1 to investigate the impact of projection effects on the study of random motions in both axisymmetric and asymmetric cases as well as with and without local anisotropy. We first built disc-like kinematics with 5 × 10⁶ points representing individual particles or gas clouds uniformly distributed over a disc of 200″ radius. The adopted axisymmetric contribution of azimuthal velocities V_θ follows V_θ(R) = v₀/(1 + (r_s/R)^γ)^1/γ for two different sets of parameters to describe either (i) a rotation curve slowly rising to a moderate maximum rotation velocity with v₀ = 200 km s⁻¹, r_s = 100″ and γ = 1, or (ii) a rotation curve with a steep inner gradient and a sharp transition to a high velocity plateau reached at small radius with v₀ = 250 km s⁻¹, r_s = 20″ and γ = 2. We considered those two rotation curves to have either a low or a high impact of large-scale rotation on BS. The disc was not assumed contracting or expanding, that is, V_R(R) = 0, and a negligible average vertical component, V_z(R) = 0. A few solutions exist to describe the departure of V_θ and V_R from axisymmetry at a given R, θ position in the plane. For example, planar asymmetries can follow Fourier harmonics of second order to mimic elliptical orbits in a bar or spiral potential (e.g. Spekkens & Sellwood 2007; Gaia Collaboration 2023). This is useful to describe large-scale perturbations of V_θ and V_R, but not necessarily appropriate to mock sharp velocity gradients across and along spiral arms (Chemin et al. 2016). Therefore, we adopted a different prescription than the simple cylindrical harmonics, and defined a spiral pattern by:

$$\begin{array}{c}\hfill \theta ={\varphi }_{\mathit{sp}}+2\pi R/r_0,\end{array}$$(B.1)

where R and θ are radial and azimuthal coordinates in the plane of the galaxy, ϕ_sp is the phase of the pattern at the centre, and r₀ is the radius interval necessary for the spiral to complete 2π. We used r₀/2π = 50″ in all models, and since the interesting parameter to probe here is the difference between the orientation of the velocity perturbation and a reference axis in the disc, we vary the initial angle of the spiral perturbation with a sampling of 15°, from 0° to 165° with respect to the major axis of the galaxy, leading to 12 possible orientations for the pattern. For each R, θ position, the closest point along the spiral in the galaxy frame was found to infer the distance D_sp to the spiral and the unitary vector u between those points. A density perturbation following the spiral pattern was introduced as a constant over-density for all points having their distance D_sp lower than a maximum distance ΔD where the spiral has an impact on the axisymmetric model. We note that the distribution of points is not important in this exercise, as it does not affect the kinematics. It can only impact the derivation of the BS model, though with negligible overall effect. We did not model the case of asymmetric vertical motions, and defined the amplitude of the velocity perturbation by ΔV = ΔV₀ × (1 − D_sp/ΔD) for D_sp ≤ ΔD, ΔV₀ being the maximum amplitude of the velocity perturbation. We used ΔV₀ = 25 km s⁻¹ and ΔD = 25″. The value of 25 km s⁻¹ was chosen to match the maximum variation of velocity through the spiral arms of the grand-design spiral Messier 99 with respect to the axisymmetric circular velocity (Fig. 8 of Chemin et al. 2016), or the maximum strength of the elliptical motions of RGB stars in the Galactic bar (Fig. 19 of Gaia Collaboration 2023). We assumed the six following configurations:

• (1) Axisymmetric models with no perturbation (ΔV = ΔV_R = ΔV_θ = 0). These models are used as references from which comparisons can be made, with respect to previous asymmetric models.

• Asymmetric velocity models with a spiral velocity perturbation along either (2) the azimuthal direction (ΔV_θ = ΔV × u.u_θ/|u.u_θ|, ΔV_R = 0), u_θ being the unitary azimuthal vector, or (3) the radial direction (ΔV_R = ΔV × u.u_R/|u.u_R|, ΔV_θ = 0), u_R being the unitary radial vector. These cases are necessary to identify the impact of velocity gradients along one planar direction, independently from the other direction. This thus represents a departure from V_θ(R) of at least ∼13%, and significantly more from V_R(R).

• Asymmetric velocity models with a spiral velocity perturbation along both azimuthal and radial directions for the cases (4) ΔV_θ = 2ΔV_R = 2ΔV × u.u_θ/|u.u_θ|, (5) ΔV_R = 2ΔV_θ = 2ΔV × u.u_R/|u.u_R|, and (6) with the spiral velocity perturbation ΔV along u, that is ΔV_R = ΔV × u.u_R and ΔV_θ = ΔV × u.u_θ. In the latter case, the orientation of the velocity perturbation changes with radius, being more azimuthal in the inner parts and more radial in the outer parts. These three models are useful to identified the combined effect of planar asymmetries with different strengths.

The azimuthal velocity at the coordinate R, θ was thus obtained by V_θ(R, θ) = V_θ(R)±|ΔV_θ(R, θ)|, and similarly for the radial component. This creates a sharp kinematic discontinuity through the spiral feature, and the asymmetries in V_θ or V_R make the velocity ellipsoid intrinsically anisotropic on scales larger than that of single gas clouds. To illustrate this, Fig. B.1 shows the velocity ellipsoid in the case of local velocity isotropy without large-scale anisotropy due to a spiral streaming (top panel), and in the case of local velocity isotropy but with anisotropy induced at large scale by the V_R streaming described above (bottom panel). The signature of the radial bias of the anisotropy is clearly seen as contours of density elongated in the radial dimension of the velocity tensor.

Fig. B.1. Velocity ellipsoid of the mock disc with local velocity isotropy for the cases without large-scale anisotropy (upper panel), and with a large-scale anisotropy induced by the spiral VR streamings (bottom panel), assuming radial motions are null on average. The white contours represent densities of 100, 1000, 5000, and 10000 toy-model particles.

Each particle was then assigned V_θ, V_R and V_z components, which were drawn randomly following a Gaussian probability distribution centred on the local values described above, V_θ(R, θ),V_R(R, θ) and V_z(R), and with standard deviation σ_θ, σ_R and σ_z, respectively. The two following cases have been explored: σ_θ = σ_R = σ_z = 8 km s⁻¹, corresponding to traditional isotropic velocity ellipsoids in the absence of asymmetries, and σ_R = 8 km s⁻¹, σ_θ = 0.7σ_R = 5.6 km s⁻¹, and σ_z = 5 km s⁻¹, which is a direct way to also generate a velocity anisotropy, uniformly, and in addition to the velocity streamings. To keep the modelling simple, we thus did not allow direct variations of σ_θ and σ_R within the disc. Combined to the six previously enumerated configurations, this leads to a total of 12 cases for each rotation curve. We did not explore the possibility σ_θ > σ_R, as this situation only occasionally occurs in the stellar populations of the Milky Way or the Large Magellanic Cloud (Gaia Collaboration 2021b, 2023).

Then, the mock cubes of data were built, giving each particle a V_los following Eq. C.1. Each data cube contains 200 spectral channels of 3 km s⁻¹ width, and 400 × 400 squared pixels, with a pixel scale of 1″. The assigned channel map of a particle corresponds to the one closest to V_los, and the line profile in a spatial pixel was obtained by summing the counts of all individual particles from this pixel. Constant position angle of the disc major axis and disc inclination were assumed as a function radius. The position angle was fixed at 45° in all models, and inclinations of 45°, 60°, and 75° were used. For each inclination, 12 mock data cubes were generated, corresponding to the 12 orientations of the spiral perturbation described above. To mimic an observational BS effect, the data cubes were then smoothed by a Gaussian function with a FWHM of 8 pixels (corresponding to ≈465 pc, assuming a distance of 12 Mpc). These are idealised data cubes to which no further smoothing by synthesised light and point spread functions were applied, and no further random noise was added. The first and second moments of the data cubes were then derived, the BS correction applied, using the low-resolution flux and velocity maps, to remain consistent with the methodology applied to the observed data.

In Figures B.2 (weak velocity gradient, moderate maximum rotation velocity) and B.3 (steep inner gradient, high maximum rotation velocity reached at small radius), we show both the effects of models of uniform isotropy and anisotropy for (1) the axisymmetric configuration, and for the spiral velocity perturbation along either (2) radial or (3) azimuthal directions. An example with the velocity perturbation along u is shown in Fig. 3 for the rotation curve with the steep gradient. We show for each case the maps obtained with an inclination of 60° and with ϕ_sp = 135°. Whereas no obvious BS residuals due to rotation are visible in Fig. B.2, stronger residuals are observed in Fig. B.3, especially in the inner parts. This is due to both the shape and amplitude of the rotation curve. Nevertheless, we checked that using the moment maps of the high-resolution mock data cubes, that is those obtained before the spatial smoothing, leads, as expected, to an accurate correction of the BS effect (see Sect. 3.1), except at the very centre for the steeper rotation curve. We also point out that when the spiral perturbation is set along the radial (azimuthal) direction, the impact is minimum along the major (minor) axis in both the velocity field and the velocity dispersion map. Fig. B.2 also shows the histograms of amplitudes and phases for orders two and four resulting from the Fourier analysis of the corresponding models, described and discussed in Sections 5.2.1 and 5.3.1.

Fig. B.2. Kinematics maps and velocity dispersion Fourier analysis histograms for the toy models with a weak velocity gradient rotation curve (v0 = 200 km s−1, rs = 100″ and γ = 1) without any perturbation (top), and with a spiral velocity perturbation towards either radial (middle) or azimuthal (bottom) components. We show, from left to right the density map, the velocity field and the BS corrected dispersion maps in the isotropic and anisotropic cases, obtained with an inclination of 60° and with ϕsp = 135°. We show from left to right amplitude histograms of orders k = 2 and k = 4 and normalised phase angle histograms of these orders, obtained from the FFT analysis of toy models dispersion maps. Blue and green histograms correspond to isotropic and anisotropic cases respectively.

Fig. B.3. Same as Fig. B.2, but for the toy models with a steep velocity gradient rotation curve (v0 = 250 km s−1, rs = 20″ and γ = 2).

Appendix C: Line-of-sight velocity dispersion for a disc in cylindrical coordinates

In this Appendix, we introduce and discuss the geometrical and mathematical framework used to study asymmetries and anisotropies in the neutral gas spatially resolved velocity dispersion fields. Velocity dispersion is a second order moment, and therefore is a more complex quantity to study than flux and radial velocity which are respectively the zero and first moment order.

Due to its hydrodynamical properties, the gas likely lies in a plane. The velocity vector in the frame of the galactic plane is described by two components both lying in this plane. Thus, in the cylindrical frame, the observed velocity along the line-of-sight is expressed as:

$$\begin{array}{c}\hfill V_{\mathrm{los}}=V_{\mathrm{sys}}+V_{\theta }cos\theta sini+V_{R}sin\theta sini+V_{z}cosi,\end{array}$$(C.1)

where V_θ is the azimuthal velocity, that is the rotation motions; V_R is the radial velocity, that is the inward or outward motions; and V_z is the vertical velocity, that is motions perpendicular to the galactic plane.

The observed velocity V_los (los for line-of-sight) is linked to the projection of V_θ, V_R and V_z along the line-of-sight through five additional parameters: PA, the position angle of the major axis of the galaxy (measured anticlockwise from the north to the direction of receding side of the galaxy); i, the inclination of the galactic disc with respect to the sky plane; V_sys, the systemic velocity of the galaxy; and x_c and y_c, the coordinates of the rotation centre in Cartesian coordinates (sky projection).

Both the radial, azimuthal, and vertical components can vary with R and θ, which are the polar coordinates in the plane of the galaxy with respect to the centre, choosing the major axis as reference θ = 0 (receding side). The azimuth in the plane of the galaxy, θ, is linked to the position angle PA, the inclination i, the position x (east-west), y (north-south) and centre x_c, y_c in the sky by the set of equations C.2 to C.7:

$$\begin{array}{c}\hfill Rcos\theta =rcos\psi ,\end{array}$$(C.2)

$$\begin{array}{c}\hfill Rsin\theta =r\frac{sin\psi }{cosi},\end{array}$$(C.3)

$$\begin{array}{c}\hfill cos\psi =\frac{(y-y_c)cos\mathit{PA}-(x-x_c)sin\mathit{PA}}{r},\end{array}$$(C.4)

$$\begin{array}{c}\hfill sin\psi =-\frac{(x-x_c)cos\mathit{PA}+(y-y_c)sin\mathit{PA}}{r},\end{array}$$(C.5)

$$\begin{array}{c}\hfill r=\sqrt{{(x-x_c)}^2+{(y-y_c)}^2},\end{array}$$(C.6)

$$\begin{array}{c}\hfill R=r\sqrt{{cos}^2\psi +\frac{{sin}^2\psi }{{cos}^{2}i}},\end{array}$$(C.7)

with ψ being the anticlockwise angle from the major axis, and r being the distance to the centre, both in the plane of the sky.

It is possible to define velocity dispersion components:

$$\begin{array}{c}\hfill {\sigma }_{\mathit{ij}}^2=\overline{V_i{V}_j}-\overline{V_i}\overline{V_j},\end{array}$$(C.8)

where the subscripts i and j denote the different coordinate directions⁶: radial (R), azimuthal (θ), and vertical (z), whereas V_i and V_j give the corresponding velocities. When i ≠ j, the term is the covariance between V_i and V_j, whereas in the specific case of i = j, the term is the variance of V_i and we note ${\sigma }_{ii}^2={\sigma }_i^2$. By definition, the line-of-sight velocity dispersion is thus defined as

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{los}}^2=\overline{V_{\mathrm{los}}^2}-{\overline{V_{\mathrm{los}}}}^2,\end{array}$$(C.9)

which leads, by replacing the terms by their corresponding expression using Eq. C.1, by using Eq. C.8, and because V_sys is constant, to:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\sigma }_{\mathrm{los}}^2=& {\sigma }_{\theta }^2{cos}^2\theta {sin}^{2}i+{\sigma }_R^2{sin}^2\theta {sin}^{2}i+{\sigma }_z^2{cos}^{2}i\hfill \\ \hfill & +2({\sigma }_{R\theta }^{2}cos\theta sin\theta {sin}^{2}i+{\sigma }_{\theta z}^{2}sin\theta cosisini\hfill \\ \hfill & +{\sigma }_{\mathit{Rz}}^{2}cos\theta cosisini),\hfill \end{array}\end{array}$$(C.10)

where ${\sigma }_{R\theta }^2,{\sigma }_{Rz}^2$ and ${\sigma }_{\theta z}^2$ are the squared cross terms of the velocity dispersion tensor. Those terms can be expressed as ${\sigma }_{R\theta }^2={\rho }_{R\theta }{\sigma }_R{\sigma }_{\theta }$, ${\sigma }_{Rz}^2={\rho }_{Rz}{\sigma }_R{\sigma }_z$, ${\sigma }_{\theta z}^2={\rho }_{\theta z}{\sigma }_{\theta }{\sigma }_z$, owing to the definition of the correlation coefficients between the components ρ_Rθ, ρ_Rz and ρ_θz that are bounded between -1 and 1. The term ${\sigma }_{R\theta }^2$ can thus be negative, depending on the sign of the correlation coefficient.

Equation C.10 can be recast in a sum of trigonometric polynomials of second degree:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\sigma }_{\mathrm{los}}^2=& (\frac{{\sigma }_{\theta }^2-{\sigma }_R^2}{2}cos2\theta +{\sigma }_{R\theta }^{2}sin2\theta ){sin}^{2}i\hfill \\ \hfill & +({\sigma }_{\mathit{Rz}}^{2}cos\theta +{\sigma }_{\theta z}^{2}sin\theta )sin2i\hfill \\ \hfill & +{\sigma }_z^2{cos}^{2}i+\frac{{\sigma }_{\theta }^2+{\sigma }_R^2}{2}{sin}^{2}i.\hfill \end{array}\end{array}$$(C.11)

We can formulate differently this equation to match better the formalism introduced with Fourier Transforms (see Sect. 3.2):

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{los}}^2=k_0+k_{1}cos(\theta -{\alpha }_1)+k_{2}cos(2(\theta -{\alpha }_2)),\end{array}$$(C.12)

with:

$$\begin{array}{c}\hfill k_0={\sigma }_z^2{cos}^{2}i+\frac{{\sigma }_{\theta }^2+{\sigma }_R^2}{2}{sin}^{2}i,\end{array}$$(C.13)

$$\begin{array}{c}\hfill k_1=\sqrt{{\sigma }_{\mathit{Rz}}^4+{\sigma }_{\theta z}^4}\times sin2i,\end{array}$$(C.14)

$$\begin{array}{c}\hfill cos{\alpha }_1=\frac{{\sigma }_{\mathit{Rz}}^2}{k_1}sin2i,\end{array}$$(C.15)

$$\begin{array}{c}\hfill sin{\alpha }_1=\frac{{\sigma }_{\theta z}^2}{k_1}sin2i,\end{array}$$(C.16)

$$\begin{array}{c}\hfill k_2=\sqrt{{\left(\frac{{\sigma }_{\theta }^2-{\sigma }_R^2}{2}\right)}^2+{\sigma }_{R\theta }^4}\times {sin}^{2}i,\end{array}$$(C.17)

$$\begin{array}{c}\hfill cos2{\alpha }_2=\frac{({\sigma }_{\theta }^2-{\sigma }_R^2)}{2k_2}{sin}^{2}i,\end{array}$$(C.18)

$$\begin{array}{c}\hfill sin2{\alpha }_2=\frac{{\sigma }_{R\theta }^2}{k_2}{sin}^{2}i.\end{array}$$(C.19)

This formalism further motivated our choice of performing Fourier Transforms on squared velocity dispersion. With these definitions, we find

$$\begin{array}{c}\hfill tan2{\alpha }_2=2{\sigma }_{R\theta }^2/({\sigma }_{\theta }^2-{\sigma }_R^2),\end{array}$$(C.20)

so that α₂ is the tilt angle of the velocity ellipsoid in (R, θ) (see e.g. Smith et al. 2009), that is, the vertex deviation (Kuijken & Tremaine 1991, 1994; Binney & Merrifield 1998), analogue to the tilt angle of the velocity ellipsoid in (R, z) plane used, for example, in Hagen et al. (2019).

In case there are asymmetries, those are embedded in the terms k_i and α_i. For the sake of simplicity, let’s assume that cross terms are null, which leads to

$$k_1=0,k_2=\frac{{\sigma }_{\theta }^2-{\sigma }_R^2}{2}\times {sin}^{2}i,\mathrm{and}{\alpha }_2=0.$$

If we only have an asymmetry of order k (sinusoidal) in both ${\sigma }_R^2$ and ${\sigma }_{\theta }^2$ with similar amplitudes σ_R, k = σ_θ, k = σ_k and with phases ϕ_R, k and ϕ_θ, k − ϕ_R, k = Δϕ_k, we can therefore express k₀ and k₂ as:

$$\begin{array}{cc}& k_0={\sigma }_z^2{cos}^{2}i+{\sigma }_k^{2}cos(\mathrm{\Delta }{\varphi }_k/2)cos(k\theta -{\varphi }_{R,k}-\mathrm{\Delta }{\varphi }_k/2){sin}^{2}i,\hfill \end{array}$$(C.21)

$$\begin{array}{cc}& k_2={\sigma }_k^{2}sin(\mathrm{\Delta }{\varphi }_k/2)sin(k\theta -{\varphi }_{R,k}-\mathrm{\Delta }{\varphi }_k/2){sin}^{2}i,\hfill \end{array}$$(C.22)

and Eq. C.12 becomes:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\sigma }_{\mathrm{los}}^2=& {\sigma }_z^2{cos}^{2}i+{\sigma }_k^{2}cos(\mathrm{\Delta }{\varphi }_k/2)cos(k\theta -{\varphi }_{R,k}-\mathrm{\Delta }{\varphi }_k/2){sin}^{2}i\hfill \\ \hfill & +{\sigma }_k^2/2sin(\mathrm{\Delta }{\varphi }_k/2)sin((k+2)\theta -{\varphi }_{R,k}-\mathrm{\Delta }{\varphi }_k/2){sin}^{2}i\hfill \\ \hfill & +{\sigma }_k^2/2sin(\mathrm{\Delta }{\varphi }_k/2)sin((k-2)\theta -{\varphi }_{R,k}-\mathrm{\Delta }{\varphi }_k/2){sin}^{2}i.\hfill \end{array}\end{array}$$(C.23)

Orders k, k − 2 and k + 2 appear. Nevertheless, even if Δϕ_k is constant, the phase ϕ_R, k is supposed to remain random, so no peculiar direction should be favoured.

The general case without the assumption that σ_R, k = σ_θ, k leads to:

$$\begin{array}{c}\hfill k_2=c_{k}cos(k\theta -{\varphi }_{R,k})+s_{k}sin(k\theta -{\varphi }_{R,k}),\end{array}$$(C.24)

with

$$\begin{array}{cc}& c_k=\frac{{\sigma }_{\theta ,k}^{2}cos\mathrm{\Delta }{\varphi }_k-{\sigma }_{R,k}^2}{2},\mathrm{and}\hfill \end{array}$$(C.25)

$$\begin{array}{cc}& s_k=\frac{{\sigma }_{\theta ,k}^{2}sin\mathrm{\Delta }{\varphi }_k}{2},\hfill \end{array}$$(C.26)

which can also be written as

$$\begin{array}{c}\hfill k_2=x_{k}cos(k\theta -{\varphi }_{R,k}-{\varphi }_k),\end{array}$$(C.27)

with

$$x_k=\sqrt{c_k^2+s_k^2},cos{\varphi }_k=c_k/x_k,\mathrm{and}sin{\varphi }_k=s_k/x_k.$$

The general case would therefore also lead to orders k − 2 and k + 2.

If σ_R, k = σ_θ, k = σ_k, we obtain that $c_k=-{\sigma }_k^2{sin}^2(\mathrm{\Delta }{\varphi }_k/2)$ and $s_k={\sigma }_k^{2}sin(\mathrm{\Delta }{\varphi }_k/2)cos(\mathrm{\Delta }{\varphi }_k/2)$, and therefore $x_k={\sigma }_k^{2}sin(\mathrm{\Delta }{\varphi }_k/2)$ and ϕ_k = Δϕ_k/2 + π/2, which enables us to recover Eq. C.22.

Appendix D: Results of the FFTs for individual galaxies in THINGS

In this section we show the various kinematic product for each THINGS galaxy as presented in Figure 1 and Figure 4. The upper row plots represents from left to right the observed flux density map, the observed velocity field, the observed velocity dispersion, the BS model, and the velocity dispersion map corrected from the BS effect. The second row shows the observed squared velocity dispersion and its modelling through FFT up to k = 4, and their residuals. The third row shows the individual orders of the FFT in the galaxy. The bottom row presents the amplitudes (left) and phase angles (right) of the FFT coefficients as a function of radius.

Fig. D.1. H I density and velocity maps and FFT results of NGC925. (Top panel, from left to right) Observed flux density map, observed velocity field, observed velocity dispersion, beam smearing model, and velocity dispersion map corrected from the beam smearing effect. (Second panel, from left to right) Squared observed velocity dispersion and its corresponding modelling through FFT up to the order 4, and the residuals between those two maps. (Third panel) Individual squared orders of the FFT modelling projected in the plane of the galaxy. (Bottom panel) Radial variation of the FFT amplitudes (left) and phase angles (right) for orders k = 0 to k = 4.

Fig. D.2. Same as Fig. D.1 for NGC2403.

Fig. D.3. Same as Fig. D.1 for NGC2903.

Fig. D.4. Same as Fig. D.1 for NGC2976.

Fig. D.5. Same as Fig. D.1 for NGC3031.

Fig. D.6. Same as Fig. D.1 for NGC3198.

Fig. D.7. Same as Fig. D.1 for NGC3521.

Fig. D.8. Same as Fig. D.1 for NGC3621.

Fig. D.9. Same as Fig. D.1 for NGC3627.

Fig. D.10. Same as Fig. D.1 for NGC4736.

Fig. D.11. Same as Fig. D.1 for NGC5055.

Fig. D.12. Same as Fig. D.1 for NGC6946.

Fig. D.13. Same as Fig. D.1 for NGC7331.

Fig. D.14. Same as Fig. D.1 for NGC7793.

All Tables

All Figures

	Fig. 1. Density and velocity maps of NGC 2841. From left to right we show: the observed flux density map, the observed velocity field, the observed velocity dispersion, the BS model, and the velocity dispersion map corrected from the BS effect.
In the text

	Fig. 2. Distribution of the fraction of missing values for the radial bins in the THINGS galaxies.
In the text

Fig. 3. Maps and histograms of toy models with a steep rotation curve (v0 = 250 km s−1, rs = 20″ and γ = 2) with a spiral perturbation. We show (top), from left to right, the density map, the velocity field, and the BS corrected dispersion maps in the isotropic and anisotropic cases, obtained with an inclination of 60° and with ϕsp = 135°. The bottom panels show amplitude and phase histograms of orders k = 2 and k = 4. Blue and green histograms correspond to the cases with local uniform isotropy with large-scale asymmetries, and local uniform anisotropy with large-scale asymmetries, respectively.

In the text

	Fig. 4. Examples of FFT results with the galaxy NGC 2841. Top: observed squared velocity dispersion and its modelling through FFT up to the order 4, and their residuals. Middle panel: individual orders of the FFT projected in the plane of the galaxy. Bottom: amplitudes (left) and phase angles (right) of FFT coefficients as a function of the radius.
In the text

Fig. 5. Results of the FFT measurements of the H I velocity dispersion maps. The histograms show the number of rings as a function of amplitude (top row) and normalised phase angle (bottom row) of the FFT harmonics of the THINGS sample (883 rings, in blue), the interpolated distribution with 20 rings per galaxy for THINGS (280 rings, in red, normalised to the maximum values of the blue histograms), and the WHISP sample (144 rings, in green). Phase angles are normalised to the period Tk of each order.

In the text

Fig. 6. Comparisons of amplitudes and phase angles between the Fourier k = 1, 2 and 4 orders, shown as differences of amplitudes versus ϕ2/T2 (upper row) and histograms of phase angles difference normalised to π (bottom row). The comparisons between the orders k = 1 and k = 2, k = 1 and k = 4, and k = 2 and k = 4 are shown in the left, middle, and right columns, respectively. In the upper row, the amplitude difference is colour coded with the logarithmic number of the tilted rings.

In the text

	Fig. 7. Histograms of phase angles (normalised by π) difference between gas density and velocity dispersion ϕk(Σ)−ϕk(σ) for k = 1, 2 and 4 (left, middle, and right panels, respectively).
In the text

Fig. 8. Histograms of the velocity dispersion amplitudes σk (k = 2, left panel) and of the normalised phase angles ϕk/Tk (k = 2, middle panel and 4, right panel) for the 7 galaxies brighter and fainter than M = −20.7 mag (red and blue colours, respectively) for the THINGS sample. The sub-samples are the same as when splitting between high (red) and low (blue) projected maximum velocity (see Sect. 5.2.4). The grey histograms are the total distributions from Fig. 5.

In the text

	Fig. 9. Histograms of amplitudes (upper row) and normalised phase angles (bottom row) resulting from the FFT analysis for the THINGS sample when splitting the sample in two, using the median significance of each order (see Sect. 5.1.1). The blue and red distributions are for the rings with a significance respectively greater and smaller than the median value.
In the text

Fig. 10. Results of the FFTs applied to the BS dispersion maps σbs, as built from disc models projected at 30° (top row) and 70° (bottom row). Left panel: map of σbs. The synthesised beam adopted in the toy model is shown on the bottom left corner. Middle panel: amplitudes of the Fourier modes for orders k = 0 to k = 4. Right panel: phase angles of the Fourier modes. A horizontal dashed line shows the half period of the Fourier modes. The disc scale-length is shown as a white contour (left) and as a vertical dashed line (middle and right).

In the text

	Fig. 11. Histograms of amplitudes (upper row) and normalised phase angles (bottom row) resulting from the FFT analysis for the THINGS sample (in blue) and for the BS dispersion maps σbs (in red).
In the text

	Fig. 12. Comparison of the BS models using optical and cm data for the galaxy NGC 3521. Left panel: Hα data based BS model. Middle panel: H I data based BS model, trimmed at the same extent as the optical model. Right panel: ratio of the BS models, defined as BS(Hα)/BS(H I).
In the text

	Fig. 13. Pixel-by-pixel comparison of the BS contribution using Hα SINGS and H I THINGS data. The y-axis shows the ratio of the BS amplitudes Hα/H I and in the x-axis is the amplitude of the H I model.
In the text

	Fig. 14. Comparison of the normalised phase angles histograms of orders k = 2 and k = 4 (left and right panels, respectively) in the velocity dispersion by subtracting once and twice the contribution from the BS effect (blue and red colours, respectively).
In the text

	Fig. 15. Histograms of normalised phase angles resulting from the FFT analysis for the THINGS sample splitting the rings as a function of radius (inner part in blue, outer part in red).
In the text

Fig. A.1. Square root of the noise on the squared velocity dispersion amplitude for orders k = 1 to 4 of the FFT as a function of σ0 for two distinct values of dσ and various number of points N in the FFT. The curves correspond to Eq. A.5 with dσ = 1 km s−1 (dashed) and dσ = 4 km s−1 (dotted), each corresponding to a different value of N (from top to bottom: 1000, 2000, 4000 and 8000). Coloured symbols correspond to different orders measured on toy models. The similar behaviour of orders leads to superposed points.

In the text

	Fig. B.1. Velocity ellipsoid of the mock disc with local velocity isotropy for the cases without large-scale anisotropy (upper panel), and with a large-scale anisotropy induced by the spiral VR streamings (bottom panel), assuming radial motions are null on average. The white contours represent densities of 100, 1000, 5000, and 10000 toy-model particles.
In the text

Fig. B.2. Kinematics maps and velocity dispersion Fourier analysis histograms for the toy models with a weak velocity gradient rotation curve (v0 = 200 km s−1, rs = 100″ and γ = 1) without any perturbation (top), and with a spiral velocity perturbation towards either radial (middle) or azimuthal (bottom) components. We show, from left to right the density map, the velocity field and the BS corrected dispersion maps in the isotropic and anisotropic cases, obtained with an inclination of 60° and with ϕsp = 135°. We show from left to right amplitude histograms of orders k = 2 and k = 4 and normalised phase angle histograms of these orders, obtained from the FFT analysis of toy models dispersion maps. Blue and green histograms correspond to isotropic and anisotropic cases respectively.

In the text

	Fig. B.3. Same as Fig. B.2, but for the toy models with a steep velocity gradient rotation curve (v0 = 250 km s−1, rs = 20″ and γ = 2).
In the text

Fig. D.1. H I density and velocity maps and FFT results of NGC925. (Top panel, from left to right) Observed flux density map, observed velocity field, observed velocity dispersion, beam smearing model, and velocity dispersion map corrected from the beam smearing effect. (Second panel, from left to right) Squared observed velocity dispersion and its corresponding modelling through FFT up to the order 4, and the residuals between those two maps. (Third panel) Individual squared orders of the FFT modelling projected in the plane of the galaxy. (Bottom panel) Radial variation of the FFT amplitudes (left) and phase angles (right) for orders k = 0 to k = 4.

In the text

	Fig. D.2. Same as Fig. D.1 for NGC2403.
In the text

	Fig. D.3. Same as Fig. D.1 for NGC2903.
In the text

	Fig. D.4. Same as Fig. D.1 for NGC2976.
In the text

	Fig. D.5. Same as Fig. D.1 for NGC3031.
In the text

	Fig. D.6. Same as Fig. D.1 for NGC3198.
In the text

	Fig. D.7. Same as Fig. D.1 for NGC3521.
In the text

	Fig. D.8. Same as Fig. D.1 for NGC3621.
In the text

	Fig. D.9. Same as Fig. D.1 for NGC3627.
In the text

	Fig. D.10. Same as Fig. D.1 for NGC4736.
In the text

	Fig. D.11. Same as Fig. D.1 for NGC5055.
In the text

	Fig. D.12. Same as Fig. D.1 for NGC6946.
In the text

	Fig. D.13. Same as Fig. D.1 for NGC7331.
In the text

	Fig. D.14. Same as Fig. D.1 for NGC7793.
In the text