H i Kinematics and Mass Distribution of Messier 33

Figure 2 presents selected channel maps of the combined DRAO+Arecibo datacube. The Galactic H i emission does not appear in these channel maps. The foreground Galactic H i is only detected at ${V}_{{\rm{hel}}}\geqslant -50\,\mathrm{km}\,{{\rm{s}}}^{-1}$ (see also Chemin et al. 2009) and does not contaminate the M33 gas emission. The variation in orientation of the contours illustrates the perturbed H i disk of M33. In particular, the contours of the lower flux density do not draw the V-shape typical of an unperturbed disk, but are elongated and twisted at their edges. In addition, the orientation of gas at both ends of the minor axis (heliocentric velocity −178 $\mathrm{km}\,{{\rm{s}}}^{-1}$) is almost perpendicular to that of the inner distribution. These are the signatures of the H i warp of M33.

Zoom In Zoom Out Reset image size

The H i integrated profile shown in Figure 3 is asymmetric, with slightly more gas in the receding southern half (3% relative to the approaching northern half). The intensity-weighted systemic velocity is −180.3 ± 2.3 $\mathrm{km}\,{{\rm{s}}}^{-1}$, as derived from any channels with velocities $\lt \,-50\,\mathrm{km}\,{{\rm{s}}}^{-1}$. We define the maximum flux of the global H i profile as the average value of the two maxima. The width of the H i profile at 50% of the maximum flux is ${W}_{50}=183\,\mathrm{km}\,{{\rm{s}}}^{-1}$, with a corresponding mid-point velocity of −180.4 $\mathrm{km}\,{{\rm{s}}}^{-1}$, similar to the intensity-weighted mean value. At the 20% level, the velocity width is ${W}_{20}=200\,\mathrm{km}\,{{\rm{s}}}^{-1}$. The total H i mass is $1.95\pm 0.36\,{10}^{9}\,{M}_{\odot }$. This value is in agreement with the mass found by Corbelli et al. (2014) from VLA–GBT data, and unsurprisingly with the mass derived by Putman et al. (2009). The total neutral gas mass is about six times higher than the molecular gas mass of $\sim 3.3\times {10}^{8}$ ${M}_{\odot }$ given in Gratier et al. (2010).

Zoom In Zoom Out Reset image size

Figure 4 presents position–velocity (PV) diagrams of the datacube, one made along the major axis of the inner disk (position angle PA of 202°, see Section 4.2), one along the major axis of the outer warped disk (PA = 165°), both of them centered on the photometric center, and another one along the direction PA = 175°, but slightly off-centered from the photometric center (bottom panel). The center of this last PV diagram is ${(\alpha ,\delta )}_{{\rm{J}}2000}=({01}^{{\rm{h}}}{33}^{{\rm{m}}}45\mathop{.}\limits^{{\rm{s}}}7,+30^\circ 42^{\prime} 51^{\prime\prime} )$. This slice orientation is chosen because it goes through two regions of higher velocity dispersion at the NNW and SSE (see Section 4.1). The width of each PV slice is equivalent to two pixels ($\sim 45^{\prime\prime}$). The yellow line in the diagrams traces a cut in the velocity field and highlights nicely to which extent the high-brightness emission traces the rotation of the disk.

Zoom In Zoom Out Reset image size

Along the direction PA = 202°, the high-brightness gas is typical of a disk with a slightly rising/constant rotation curve (offsets $\leqslant 30^{\prime} \mbox{--}35^{\prime}$). The brightness contours then exhibit decreasing velocities at absolute offsets $\gt 35^{\prime}$. This apparent decrease occurs in the same region in which the twist of the iso-brightness contours is observed in the channel maps, as caused by the M33 warp. Along the direction PA = 165°, the high-brightness emission again follows a slowly rising/flat rotational pattern.

The most important result here is the detection of a low-brightness H i component that has an extremely anomalous velocity behavior with respect to the high-density gas. The contours of low-density gas display "beard-like" features in the diagrams, in reference to earlier works (Fraternali et al. 2001; Sancisi et al. 2001, and references therein), which means that the contours are systematically stretched toward a line-of-sight velocity that deviates from those of contours of gas with the highest density.

First, the low-brightness emission extends toward the systemic velocity of the system, as seen in the top left and bottom right quadrants of all diagrams (open arrows). Second, this low-brightness gas extends to higher velocities than the high-brightness disk contours (filled arrows). This is more obvious for PA = 165° and 175° than for PA = 202°. In this PV diagram, only an absolute offest $\sim 5^{\prime}$ shows a high-velocity bump. Third, the low-brightness emission leaks in the forbidden velocity zones, which causes it to be in apparent counterrotation with respect to the high-density gas. In our diagrams, the forbidden zones are the bottom left and top right quadrants. A H i component in the approaching disk half is in the forbidden velocity zone when its radial velocity is higher than the systemic velocity ("receding component on the approaching side"). Conversely, a H i component in the receding disk half is in the forbidden velocity zone when its radial velocity is lower than the systemic velocity ("approaching component on the receding side"). The forbidden velocity gas is detected at low offset and seems to make a link between the low- and high-velocity gas. Note also the isolated component at 30' for PA = 202°. We show in Figure 5 the gas component with forbidden velocity by extracting H i profiles at several positions along the PA = 165° and 202° directions. Only channels with ${T}_{{\rm{B}}}\gt 5\sigma$ are shown for clarity. The equatorial coordinates of these spectra are given in Table 3. The H i gas in the forbidden zone is the asymmetric tail of the profiles, with velocities differing by up to 80–100 $\mathrm{km}\,{{\rm{s}}}^{-1}$ from the main peak. H i gas in the PA = 202° forbidden zone is also observed as distinct H i peaks with velocities differing by 100 $\mathrm{km}\,{{\rm{s}}}^{-1}$ from the main peak (offsets = −30', +10').

Zoom In Zoom Out Reset image size

Table 3.  Coordinates of Five Representative Spectra Exhibiting a H i Velocity Component in the Forbidden Velocity Zones from the PV Diagrams of Figure 4

PV Diagram	Offset	Coordinates (R.A., Decl.) (J2000)
PA = 165°	$+7^{\prime}$	${01}^{{\rm{h}}}{33}^{{\rm{m}}}40\mathop{.}\limits^{{\rm{s}}}7$, +30°31'56''
PA = 165°	$-11^{\prime}$	${01}^{{\rm{h}}}{33}^{{\rm{m}}}22\mathop{.}\limits^{{\rm{s}}}0$, +30°49'46''
PA = 202°	$+10^{\prime}$	${01}^{{\rm{h}}}{33}^{{\rm{m}}}17\mathop{.}\limits^{{\rm{s}}}1$, +30°29'43''
PA = 202°	$-12^{\prime}$	${01}^{{\rm{h}}}{33}^{{\rm{m}}}52\mathop{.}\limits^{{\rm{s}}}6$, +30°50'09''
PA = 202°	$-30^{\prime}$	${01}^{{\rm{h}}}{34}^{{\rm{m}}}26\mathop{.}\limits^{{\rm{s}}}7$, +31°07'49''

Download table as:  ASCIITypeset image

The first anomalous component is reminiscent of the slow-rotation extraplanar H i layer seen in other galaxies (e.g., NGC 2403, NGC 891, and NGC 253, Fraternali et al. 2002; Oosterloo et al. 2007; Lucero et al. 2015). The second component recalls the population of low-mass high-velocity clouds in the halo of our Galaxy (Wakker & van Woerden 1997) or M31 (Westmeier et al. 2008). Putman et al. (2009) have reported high-velocity gas around M33, but at larger distances than the new gas detected here, well outside the field of view of our observtions. These authors proposed a scenario with extraplanar gas falling onto M33 after a tidal event with M31. As for the third component, it would not be the first time that forbidden velocity gas is reported in nearby galaxies (NGC 2403, Fraternali et al. 2001, 2002). Interestingly, Fraternali and collaborators reported that the low angular momentum and the forbidden velocity components in NGC 2403 seem to form a coherent structure. M33 is thus similar to NGC 2403 in this aspect.

The detailed derivation of the distribution, kinematics, and mass of the three anomalous components is beyond the scope of this article and will be presented in a future paper (L. Chemin et al. 2017, in preparation). This work will require a more appropriate processing of the datacube into several components than the single-moment maps analysis we did here, possibly like in Fraternali et al. (2002) with NGC 2403, who fit to and subtracted from the datacube a Gaussian profile centered on the highest H i peak, or like in Chemin et al. (2012), who fit multiple Gaussian peaks to the current M33 data set.

The task moments in gipsy has been used to compute the moment maps using the datacube smoothed at $120^{\prime\prime} \times 120^{\prime\prime}$. Figures 6 and 7 show the maps of the zeroth moment (integrated emission), first moment (velocity field), and second moment (velocity dispersion).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The H i emission map shows that the high surface density gas is contained within the stellar disk ($R\lesssim 30\mbox{--}35^{\prime}$, or ≲8 kpc). It displays a multiple spiral arms pattern that coincides well with the spiral structure observed in the Hα disk (Kam15) and molecular gas disk (Druard et al. 2014). Beyond the stellar disk, the gas column density sharply decreases, reaching $\sim 5\times {10}^{18}$ cm⁻² in the outermost regions. Gaseous tails and arcs are observed to the northwest and southeast of the disk as part of the H i warp. The SE feature is less extended than the NW arc-like structure, which roughly points in the direction of M31. These perturbations could be signatures of the past interaction between the two galaxies. Note also that M33 is surrounded by a prominent stellar structure that provides additional evidence of an encounter with M31 (McConnachie et al. 2010), or at least of recent gravitational interaction. This stellar structure extends ∼2° (30 kpc in projection) to the northwest toward M31, nearly three times farther out than the size of the M33 stellar disk, and thus farther than the H i gas in the warp.

Our H i map is in very good agreement with atomic gas distributions seen elsewhere (Corbelli et al. 2014). The H i mass surface density profile derived from the column density map with the adopted kinematical parameters given below is shown in Figure 8. The ${{\rm{H}}}_{2}$ mass surface density from Druard et al. (2014) is also shown. Naturally, ${{\rm{H}}}_{2}$ is much more concentrated than the atomic gas. Both gas components reach surface densities ∼10 ${M}_{\odot }$ pc⁻² (8 for H i and 12 for ${{\rm{H}}}_{2}$) in the inner disk regions.

Zoom In Zoom Out Reset image size

The velocity field is quite regular in the inner regions ($R\lesssim 8$ kpc). Beyond this radius, in the regions where the gas distribution is perturbed, the velocity field is irregular. This is apparent in the prominent twist of the velocity contours, whose feature was also visible in the channel maps (Section 3). Note also an outer southwestern cloud or extension, highlighted in yellow and orange in the velocity map, that rotates more slowly than gas inside the disk at similar radius. This extension has been discussed by Putman et al. (2009) (see also Grossi et al. 2008).

In the velocity dispersion map, an incomplete ring of larger dispersion is observed in the transition between the high column density inner disk and the low column density outer disk. In this structure, profiles exhibit multiple H i peaks or wider single H i peaks (see also Chemin et al. 2012). These features are likely due to the crowding of gas orbits caused by the disk warping (see Section 4.4 and Figure 11).

The task rotcur in gipsy was used to derive the kinematical parameters and the rotation curve. Assuming negligible radial motions, the tilted-ring model fits the expression

$$\begin{eqnarray}&&{V}_{\mathrm{obs}}={V}_{\mathrm{sys}}+{V}_{\mathrm{rot}}\cos \theta \sin i\end{eqnarray} \tag{ 1 }$$

to the H i velocity field, where θ is the azimuthal angle in the plane of the galaxy and i the inclination. The position angle of the kinematical major axis is defined as the counterclockwise angle in the plane of the sky from the north to the receding-side semimajor axis. The angle θ is measured relatively to the semimajor axis. The southwestern cloud has been masked for the derivation of the parameters and rotation curve.

In a first run, the position of the rotation center and the systemic velocity ${V}_{\mathrm{sys}}$ were let free to vary, using fixed kinematical inclination and major-axis PA at the optical values (52° and 202°, respectively). We find a systemic velocity of ${V}_{\mathrm{sys}}\sim -183\,\mathrm{km}\,{{\rm{s}}}^{-1}$, and a center at ${(\alpha ,\delta )}_{{\rm{J}}2000}=({01}^{{\rm{h}}}{33}^{{\rm{m}}}50\mathop{.}\limits^{{\rm{s}}}9,+30^\circ 40^{\prime} 20^{\prime\prime} )$. Since the positional difference with respect to the photometric center is much smaller than the beam size, we decided to fix the kinematical center at the position of the photometric center. Moreover, because the inferred systemic velocity is in agreement with the mean value determined from the integrated profile (Section 4.1), we adopted ${V}_{\mathrm{sys}}=-180$ $\mathrm{km}\,{{\rm{s}}}^{-1}$ as given by the profile (both intensity-weighted mean and mid-point velocity). A second run allowed us to fit the inclination and major-axis PA using the fixed center and systemic velocity given by the adopted values. The results of the tilted-ring model are shown in Figure 9, and the adopted inclination and position angle are listed in Table 4.

Zoom In Zoom Out Reset image size

Table 4.  Results of the Tilted-ring Model of the H i Velocity Field of M33

Radius	Radius	${V}_{\mathrm{rot}}$	${\rm{\Delta }}{V}_{\mathrm{rot}}$	i	PA	Radius	Radius	${V}_{\mathrm{rot}}$	${\rm{\Delta }}{V}_{\mathrm{rot}}$	i	PA
(')	(kpc)	($\mathrm{km}\,{{\rm{s}}}^{-1}$)	($\mathrm{km}\,{{\rm{s}}}^{-1}$)	(°)	(°)	(')	(kpc)	($\mathrm{km}\,{{\rm{s}}}^{-1}$)	($\mathrm{km}\,{{\rm{s}}}^{-1}$)	(°)	(°)
2	0.5	42.0	2.4	53.2	201.3	50	12.2	115.7	9.6	55.6	173.8
4	1.0	58.8	1.5	53.3	201.3	52	12.7	115.1	7.7	55.7	173.4
6	1.5	69.4	0.4	53.4	201.3	54	13.2	117.1	5.1	55.8	172.9
8	2.0	79.3	4.0	53.5	201.3	56	13.7	118.2	3.2	55.9	172.5
10	2.4	86.7	1.8	53.6	201.3	58	14.2	118.4	1.4	56.0	172.1
12	2.9	91.4	3.1	53.7	201.3	60	14.7	118.2	1.8	56.1	171.6
14	3.4	94.2	4.8	53.8	201.3	62	15.1	117.5	2.4	56.2	171.2
16	3.9	96.5	5.5	53.9	201.3	64	15.6	119.6	0.8	56.3	170.8
18	4.4	99.8	3.9	54.0	201.3	66	16.1	118.6	1.5	56.4	170.3
20	4.9	102.1	1.7	54.1	201.3	68	16.6	122.6	0.5	56.5	169.9
22	5.4	103.6	0.4	54.2	201.3	70	17.1	124.1	2.9	56.6	169.5
24	5.9	105.9	0.7	54.3	201.3	72	17.6	125.0	2.2	56.7	169.0
26	6.4	105.7	1.7	54.4	201.3	74	18.1	125.5	2.5	56.8	168.6
28	6.8	106.8	2.2	54.5	201.3	76	18.6	125.2	8.1	56.9	168.2
30	7.3	107.3	3.0	54.6	201.3	78	19.1	122.0	9.8	57.0	167.7
32	7.8	108.3	4.0	54.7	198.8	80	19.5	120.4	8.5	57.1	167.3
34	8.3	109.7	4.0	54.8	195.4	82	20.0	114.0	26.6	57.2	166.8
36	8.8	112.0	4.8	54.9	192.0	84	20.5	110.0	34.6	57.3	166.4
38	9.3	116.1	2.2	55.0	188.7	86	21.0	98.7	27.4	57.5	166.0
40	9.8	117.2	2.5	55.1	185.3	88	21.5	100.1	33.4	57.6	165.5
42	10.3	116.5	6.5	55.2	181.9	90	22.0	104.3	35.2	57.7	165.1
44	10.8	115.7	8.1	55.3	178.5	92	22.5	101.2	27.4	57.8	164.7
46	11.2	117.4	8.2	55.4	175.2	94	23.0	123.5	39.1	57.8	164.3
48	11.7	116.8	8.9	55.5	174.2	96	23.5	115.3	26.7	57.9	163.8

Note. i and PA are the adopted H i inclination and major-axis position angle and ${V}_{\mathrm{rot}}$ the resulting H i rotation curve with associated velocity uncertainties (see text for details).

Download table as:  ASCIITypeset image

The H i rotation curve (Figure 9) is then derived with fixed values for the kinematical center and systemic velocity, and the adopted models of inclination and position angle (solid lines in the middle and bottom panels). Also shown are the obtained H i rotation curves for the approaching (${V}_{{\rm{a}}}$) and receding sides (${V}_{{\rm{r}}}$) fit separately with similar kinematical parameters. Differences between ${V}_{{\rm{a}}}$ and ${V}_{{\rm{r}}}$ are mostly smaller than 10 $\mathrm{km}\,{{\rm{s}}}^{-1}$. The largest differences, up to 60 $\mathrm{km}\,{{\rm{s}}}^{-1}$, occur at $R\gt 80^{\prime}$ in the perturbed NW and SE structures from the warp.

The H i rotation curve is reported in Table 4. The total velocity uncertainty ${\rm{\Delta }}{V}_{\mathrm{rot}}$ is defined by ${\rm{\Delta }}{V}_{\mathrm{rot}}^{2}={\epsilon }^{2}+| ({V}_{{\rm{a}}}-{V}_{{\rm{r}}})/2{| }^{2}$, with being the formal rms error for both sides of the model (solid line in the top panel of Figure 9). The dominant error of the total uncertainty is the velocity difference between the two disk sides. Because of this definition, the errors appear larger than in many other studies. However, we prefer being conservative because we model the rotation curve with axisymmetric components (Section 5) and consider that our definition is more representative of the true asymmetry in the observed kinematics and of uncertainties when doing mass models.

Figure 10 shows the line-of-sight dispersion profile of M33, derived by azimuthally averaging the dispersion field using the adopted inclination and position angle profiles. Table 8 of Appendix A lists the mean dispersion. On average, the dispersion is ∼9 $\mathrm{km}\,{{\rm{s}}}^{-1}$ in the central regions ($R\lt 30^{\prime}$). The signature of the ring of higher dispersion is observed as a rise of 3 $\mathrm{km}\,{{\rm{s}}}^{-1}$ ($R=35^{\prime} \mbox{--}40^{\prime}$). The dispersion then continuously decreases at larger radius. Cuts made along the semimajor axis of the approaching and receding sides in the dispersion map show significant differences inside $R=30^{\prime}$ and beyond $R=50^{\prime}$. Moreover, the mean dispersion is higher than the values along the semimajor axes. Such differences demonstrate the intrinsic asymmetry of the dispersion field, similarly to the disturbed velocity field, as caused by any of the spiral, warp, and lopsidedness perturbations (see Sections 4.4 and 4.5).

Zoom In Zoom Out Reset image size

The warping of M33 is mainly observed as a significant kinematical twist of the major-axis position angle, starting by a sharp decrease from ∼202° to ∼170° within $7\lt R\lt 11\,\mathrm{kpc}$, and then a smooth variation down to 165° at large radius. A kinematical tilt is also observed, although less spectacular, as a small rise of inclination (∼5° from the center to the disk outskirts).

The geometry of the H i orbits implied by the axisymmetric tilted-ring model is shown in Figure 11. Within the transition region of significant PA twist, the modeled rings are very close to each other. Gas clouds are likely colliding at these radii. This statement is confirmed by the observation of double H i peaks at, e.g., locations C and D, or more generally by the incomplete ring-like structure of higher velocity dispersion (Figure 6). The change in inclination is so insignificant here that a configuration in which the line-of-sight crosses the disk more than once is hardly possible. Moreover, we verified that the wide and double H i peaks are not caused by the "low" resolution of the DRAO data. Indeed, first the velocity gradient is not important in this region because the rotation curve is almost flat. Then, such peculiar H i profiles are also observed in the 12'' VLA datacube of Gratier et al. (2010), as seen at position C (the location D is at the periphery of the VLA field of view, the corresponding profile is only made of noise).

Zoom In Zoom Out Reset image size

It is worthwhile to mention that not all wider profiles can be explained by our idealized model in the transition region because some of them are also observed elsewhere than at the location of crowded rings. This is seen at angular offsets ∼30'–35' in the PA = 175° PV diagram of Figure 4, whose offsets correspond to ${(\alpha ,\delta )}_{{\rm{J}}2000}=({01}^{{\rm{h}}}{33}^{{\rm{m}}}57\mathop{.}\limits^{{\rm{s}}}6,+30^\circ 05^{\prime} 02^{\prime\prime} )$ and $({01}^{{\rm{h}}}{33}^{{\rm{m}}}20\mathop{.}\limits^{{\rm{s}}}2,+31^\circ 12^{\prime} 30^{\prime\prime} )$. Here, these wider profiles coincide with holes in the gas distribution and trace high- and low-velocity faint gas (Section 3.2).

Except in the transition region, the inferred model geometry is regular and no wide profiles are observed (positions A and B in the inner disk, positions E and F in the outer disk).

Interestingly, the PA twist model implies that the major axes of the outermost ellipses are aligned with the direction M33–M31. This could be a consequence of the gravitational interaction between the two galaxies. In this scenario, the interaction stretched the outer H i disk in the direction of M31, and this perturbation is acting down to a radius of $R\sim 7\,\mathrm{kpc}$. Circular orbits are excluded in these outer regions as circularity cannot create precessed and closely grouped rings such as those implied by the projected model. Gas orbits are likely elongated and lopsided beyond this radius.

The model velocity map based on the H i RC allows us to make a residual velocity map, defined here as the observation minus the model velocity field (Figure 12). No systematic velocity asymmetry is found, which shows the goodness of the adopted mass center and systemic velocity. Significant residuals are tightly linked to the warp perturbation described above. The lowest residuals occur in the limit of the unperturbed inner H i disk, at the level of a few $\mathrm{km}\,{{\rm{s}}}^{-1}$. The largest residuals are observed in the outer disk, starting from the transition zone in which the major-axis PA varies abruptly, up to ∼30 $\mathrm{km}\,{{\rm{s}}}^{-1}$ (absolute value). The ring of higher velocity dispersion (Figure 6) thus coincides with large residuals. Note also the significant residual for the SW clump ($\gt 30\,\mathrm{km}\,{{\rm{s}}}^{-1}$, absolute value), confirming that it does not rotate in the same way as the disk. We finally note that the larger residuals we observe in the NW and SE H i extensions confirm the asymmetric gas orbits at the disk periphery.

Zoom In Zoom Out Reset image size

We can assess the perturbations in the H i velocity field of M33 further by expanding the standard model of Equation (1) to higher-order cosine and sine terms. Franx et al. (1994) and Schoenmakers et al. (1997) were among the first to propose Fourier analyses of velocity fields to estimate asymmetries caused by, e.g., oval distortions, lopsided potentials, spiral arms, or warps. These authors argued that a perturbing mode of the gravitational potential of order m generates, to the first order, $k=m-1$ and $k=m+1$ Fourier components in velocity fields.

The following Fourier series model was thus fit to the H i velocity field:

$$\begin{eqnarray}&&{V}_{\mathrm{obs}}={c}_{0}+\displaystyle \sum _{{\rm{k}}=1}({c}_{k}\cos k\theta +{{\rm{s}}}_{k}\sin k\theta )\sin i,\end{eqnarray} \tag{ 2 }$$

where c_k and s_k are the velocity coefficients of harmonic order k (k is an integer). The c₀ coefficient is the systemic velocity ${V}_{\mathrm{sys}}$ of Equation (1), the c₁ coefficient is the rotation curve, the s₁ term is the noncircular radial velocity, while higher-order terms constrain deviations from axisymmetry of c₁ and s₁.

We restrict the model to k = 4, implying that a perturbation of the potential on the order of up to m = 3, possibly m = 5, is likely to be detected, if it exists. For the harmonic model the inclination and position angle were fixed at the values of the tilted-ring model of Section 4.2. The angular sampling is ∼490 pc (120''). We derived the amplitudes of the noncircular and asymmetric motions, given by ${v}_{1}=| {{\rm{s}}}_{1}|$ and ${v}_{k}={({c}_{k}^{2}+{{\rm{s}}}_{k}^{2})}^{1/2}$, respectively, for $k\gt 1$, and ${v}_{0}=\sqrt{{({c}_{0}-{V}_{\mathrm{sys}})}^{2}}$. These amplitudes are shown in Figure 13. The uncertainties are the $1\sigma$ formal errors from the fit.

Zoom In Zoom Out Reset image size

One first observes increasing amplitudes for every k terms at $R\gtrsim 15\,\mathrm{kpc}$. These may be other signatures of the interaction with M31 on the outer M33 velocity field. Interestingly, a similar trend was observed for v₂, v_3, and v₄ in the outer velocity field of another grand-design spiral galaxy, Messier 99, which is also perturbed by its environment (Chemin et al. 2016).

Then, the variations of v₀ and v₂ are similar. Larger amplitudes are observed in the inner 5 kpc, as well as at $R\sim 8$ and 12 kpc. This is evidence of a kinematical lopsidedness, i.e., an m = 1 perturbation of the M33 gravitational potential. Larger k = 0 and k = 2 terms in the innermost disk region are other similarities with M99. The particularity of M33, however, is that the difference between the rotation curves for the approaching and receding disk sides ($R\lt 5$ kpc) is not as prominent as in M99. Chemin et al. (2016) showed from an asymmetric 3D mass model that the central peak of these even terms cannot be modeled solely by a lopsided gravitational potential of luminous matter, and proposed that DM may also be lopsided in the innermost regions of M99. We refer to a future work a similar modeling of asymmetric mass distribution of M33 to investigate whether lopsided gas and stellar potentials are enough to explain the kinematical asymmetry.

Second, the variations of v₁ and v₃ are similar. The k = 1 and k = 3 terms are small inside R = 7 kpc and larger at $R\sim 11\,\mathrm{kpc}$, within the zone of strong major-axis PA variation. Lower values in the inner disk perfectly reflect the observation that no bisymmetric structure dominates in the inner density map (Figure 6). On the other hand, higher values at large radius reflect the dynamical impact of the warp of M33.

Third, a bump is observed for the k = 4 term at R ∼ 8 kpc that coincindes with the bumps seen for the k = 0 and k = 2 terms. Whether this feature is caused by the same m = 1 perturbation remains unclear. It is unlikely to be caused by higher-order perturbations, as no m = 3 or m = 5 modes are observed in the gaseous and stellar density maps.

The average amplitude of the noncircular motion is $\langle {v}_{1}\rangle \,=3.8\pm 0.6\,\mathrm{km}\,{{\rm{s}}}^{-1}$. The average amplitudes of the asymmetries are $\langle {v}_{0}\rangle =3.3\pm 0.3\,\mathrm{km}\,{{\rm{s}}}^{-1}$, $\langle {v}_{2}\rangle =4.9\pm 0.6\,\mathrm{km}\,{{\rm{s}}}^{-1}$, $\langle {v}_{3}\rangle \,=4.9\pm 0.7\,\mathrm{km}\,{{\rm{s}}}^{-1}$, and $\langle {v}_{4}\rangle =3.7\pm 0.4\,\mathrm{km}\,{{\rm{s}}}^{-1}$, again showing the more important impact of the m = 1 and m = 2 perturbations on the gravitational potential of M33. The level of these asymmetries is consistent with the asymmetry measured between the rotation curves of the approaching and receding disk sides.

4.6.1. DRAO-Arecibo versus VLA–GBT

The significant twist of the major axis of the H i velocity field we measure at large radius from our DRAO-Arecibo data set is very consistent with results obtained from VLA–GBT data by Corbelli et al. (2014), or from earlier Arecibo 21 cm measurements by Corbelli & Schneider (1997). However, the comparison with the kinematical tilt found in Corbelli et al. (2014) is made difficult because these authors found two different trends of inclination variation at large radii. On the one hand, Corbelli et al. (2014) found a disk inclination that decreases, as based on a three-dimensional modeling of the H i datacube. On the other hand, they found a slightly increasing inclination from more traditional tilted-ring models of the velocity field. Our warp model is therefore in agreement with this part of their modeling, but not entirely with their 3D datacube modeling.

Comparing our rotation curve with the curve from Corbelli et al. (2014) is not straightforward either because their curve was derived from these two different models of a tilted disk. The DRAO and VLA rotation curves are shown in Figure 14 and agree well within the uncertainties for $R\leqslant 20\,\mathrm{kpc}$. Then the shapes unsurprinsingly differ beyond R = 20 kpc. A more detailed investigation shows that when we compare our result with the curve they obtained from a rotcur-like model that is similar to ours, the shapes become comparable at these radii, showing a drop out to $R\sim 22\,\mathrm{kpc}$, followed by a rise. This means that the action of the 3D modeling they made of the datacube causes the shape difference at the largest radii.

Zoom In Zoom Out Reset image size

It is clear, however, that both H i rotation curves are perturbed at these radii, harboring larger scatter and uncertainties. We conclude that the H i rotation curve of M33 is only reliable out to ∼20 kpc, regardless of the method used to derive it. Fortunately, this outermost region has negligible impact on the mass distribution models described in Section 5 because we used the inverse of the squared errors as weighting function of the rotation velocities.

4.6.2. The Hybrid Hα–H i Rotation Curve of M33

Kam15 showed that the M33 stellar bulge has a negligible impact on the mass distribution (see also Corbelli & Walterbos 2007; Corbelli et al. 2014). Therefore only the velocity contribution from the stellar disk, ${V}_{\star }$, is needed in the present study. It has been derived with the task rotmod of gipsy (van der Hulst et al. 1992) from near-infrared surface photometry, which is well known to give the best representation of the old stellar disk population that contributes most to the stellar mass. Kam15 derived the surface brightness profile of M33 at 3.6 μm from Spitzer/IRAC data, where contaminating bright stars have been removed. More complete details on the derivation of the surface brightness profile are given in Kam15. The disk mass-to-light ratio at 3.6 μm, ϒ, has been estimated from infrared color and stellar population synthesis (SPS) models following the prescriptions given by Oh et al. (2008) and de Blok et al. (2008):

$$\begin{eqnarray}&&{\rm{\Sigma }}[{M}_{\odot }\ {\mathrm{pc}}^{-2}]={\rm{\Upsilon }}\times {10}^{-0.4({\mu }_{3.6}-{C}_{3.6})},\end{eqnarray} \tag{ 3 }$$

where ${\mu }_{3.6}$ is the surface brightness and ${C}_{3.6}=24.8$ is a correction value as given in Oh et al. (2008). Using the J − K color index from Jarrett et al. (2003), Kam15 derived ${\rm{\Upsilon }}=0.72\pm 0.1$. This corresponds to a stellar mass of $(7.6\pm 1.1)\ {10}^{9}\,{M}_{\odot }$, which is 38% (58%) higher than the fixed (best-fit) stellar mass given in Corbelli et al. (2014). Note that ${\rm{\Upsilon }}\sim 0.7$ is the upper limit expected from mass distribution modeling of the large galaxy sample of Lelli et al. (2016).

This discrepancy is the reason why we have also performed models with ${\rm{\Upsilon }}=0.52$ to facilitate comparisons with the results of Corbelli et al. (2014). This value of ϒ was obtained by scaling the maximum of the stellar velocity curve to ∼70 $\mathrm{km}\,{{\rm{s}}}^{-1}$, which corresponds to the highest contribution from stars in the fixed stellar mass fits of Corbelli et al. (2014). ${\rm{\Upsilon }}=0.52$ is valid here only in the Spitzer/IRAC 3.6 μm band, and it corresponds to a stellar mass of $5.5\ {10}^{9}$ ${M}_{\odot }$, which thus agrees perfectly with the value from the SPS model of Corbelli et al. (2014). Assuming ${\rm{\Upsilon }}=0.52$ that is constant with radius in our models is a good way to account for a stellar velocity contribution that is fairly similar in every aspects to Corbelli et al. (2014). Note also that it corresponds to the mean value found for the sample of Lelli et al. (2016), and appears more in agreement with most galaxies of their sample as bright as M33 than ${\rm{\Upsilon }}=0.72$.

The velocity contribution from the disk of neutral gas, ${V}_{\mathrm{atom}}$, is derived from our deep DRAO+Arecibo data set by integrating the total atomic gas mass surface densities (Figure 8) multiplied by a factor of ∼1.3 to take into account the helium contribution.

The molecular gas that is traced by the CO line is mainly concentrated in the innermost kpcs. Molecular gas in M33 has been observed by Tosaki et al. (2011) with the Nobeyama Radio Observatory at 19'' resolution, and by Gratier et al. (2010) and Druard et al. (2014) with the IRAM 30 m dish at a resolution of $\sim 12^{\prime\prime} \mbox{--}15^{\prime\prime}$. Tosaki et al. (2011) reported that the H i and CO peaks are not always correlated and that the density of the atomic gas is higher than that of the molecular gas in the inner parts, while Gratier et al. (2010) showed that the H i density is lower than that of the molecular gas. Using a conversion factor of $N({{\rm{H}}}_{2})/{I}_{\mathrm{CO}(1\to 0)}=4\,\times \,{10}^{20}\,{\mathrm{cm}}^{-2}/$(K $\mathrm{km}\,{{\rm{s}}}^{-1}$), twice the value found for the Milky Way (M33 has half the solar metallicity), they measured an average density of ${{\rm{\Sigma }}}_{{{\rm{H}}}_{2}}=8.5\pm 0.2$ ${M}_{\odot }$ pc⁻² for the central kpc, and a total molecular gas mass of $\sim 3.3\times {10}^{8}$ ${M}_{\odot }$ for the entire M33 disk. We used the ${{\rm{H}}}_{2}$ density profile derived from Druard et al. (2014), scaled by a factor of ∼1.3 to infer the velocity contribution of the total molecular gas component, ${V}_{\mathrm{mol}}$.

The total rotation velocity ${V}_{\mathrm{rot}}$ in the models with a DM component is defined by

$$\begin{eqnarray}&&{V}_{\mathrm{rot}}^{2}={V}_{\star }^{2}+{V}_{\mathrm{gas}}^{2}+{V}_{\mathrm{DM}}^{2},\end{eqnarray} \tag{ 4 }$$

where ${V}_{\star }$ is the contribution from the stellar disk and ${V}_{\mathrm{gas}}={({V}_{\mathrm{atom}}^{2}+{V}_{\mathrm{mol}}^{2})}^{1/2}$ the total contribution from the gaseous disk, as deduced in Section 5.2, and ${V}_{\mathrm{DM}}$ from a DM halo that is assumed to be spherical.

5.3.1. The Pseudo-isothermal Sphere Model

Here, the density profile of DM is given by

$$\begin{eqnarray}&&{\rho }_{\mathrm{ISO}}(R)=\displaystyle \frac{{\rho }_{0}}{1+{(R/{R}_{c})}^{2}}.\end{eqnarray} \tag{ 5 }$$

The corresponding circular velocities are

$$\begin{eqnarray}&&{V}_{\mathrm{ISO}}(R)=\sqrt{4\pi G{\rho }_{0}{R}_{c}^{2}(1-R/{R}_{c}\arctan (R/{R}_{c}))},\end{eqnarray} \tag{ 6 }$$

where ${\rho }_{0}$ and R_c are the central density and the core radius of the halo, respectively. We can describe the steepness of the inner mass density profile by a power law $\rho \sim {R}^{\alpha }$. In the case of the ISO halo, $\alpha \to 0$ (the halo core has a density that is almost constant).

5.3.2. The Navarro–Frenk–White Model

The NFW model—the so-called "universal halo"—is deduced from cold dark matter simulations (Navarro et al. 1997). The density profile is cuspy, following a $\rho \propto {R}^{-1}$ law in the center, and is given by

$$\begin{eqnarray}&&{\rho }_{\mathrm{NFW}}(R)=\displaystyle \frac{{\rho }_{i}}{R/{R}_{{\rm{s}}}{(1+R/{R}_{{\rm{s}}})}^{2}},\end{eqnarray} \tag{ 7 }$$

where ${\rho }_{i}\approx 3{H}_{0}^{2}/(8\pi G)$ is the critical density for closure of the universe and R_s is a scale radius. The velocity contribution corresponding to this halo is given by

$$\begin{eqnarray}&&{V}_{\mathrm{NFW}}(R)={V}_{200}\,\sqrt{\displaystyle \frac{\mathrm{ln}(1+{cx})-{cx}/(1+{cx})}{x(\mathrm{ln}(1+c)-c/(1+c))}},\end{eqnarray} \tag{ 8 }$$

with ${V}_{200}$ that is the velocity at a radius ${R}_{200}$ at which the density is 200 times that for closure of the universe, $c={R}_{200}/{R}_{{\rm{s}}}$ gives the concentration parameter of the halo, and $x=R/{R}_{{\rm{s}}}$.

An alternative to DM to explain the missing mass problem is MOND (Milgrom 1983a, 1983b). MOND has been successful in correctly reproducing many galaxy rotation curves (e.g., Sanders & Verheijen 1998; Gentile et al. 2010). It postulates that in a regime of acceleration much smaller than a universal constant acceleration, a₀, the classical Newtonian dynamics is no longer valid and the law of gravity is modified.

In the MOND framework, the gravitational acceleration of a test particle is given by

$$\begin{eqnarray}&&\mu (x=g/{a}_{0})g={g}_{N},\end{eqnarray} \tag{ 9 }$$

where g is the acceleration, g_N is the Newtonian acceleration, and $\mu (x)$ is an interpolating function that must satisfy $\mu (x)=x$ for $x\ll 1$ and $\mu (x)=1$ for $x\gg 1$. The MOND velocity profile thus depends on $\mu (x)$. For the "standard" μ-function proposed by Milgrom (1983a),

$$\begin{eqnarray}&&\mu (x)=\displaystyle \frac{x}{\sqrt{1+{x}^{2}}},\end{eqnarray} \tag{ 10 }$$

the MOND rotation velocity is

$$\begin{eqnarray}&&{V}_{\mathrm{STD}}={V}_{\mathrm{lum}}\sqrt{\sqrt{0.5\ (1+\sqrt{1+{(2{a}_{0}R/{V}_{\mathrm{lum}}^{2})}^{2}})}},\end{eqnarray} \tag{ 11 }$$

with ${V}_{\mathrm{lum}}^{2}={V}_{\star }^{2}+{V}_{\mathrm{gas}}^{2}$, ${V}_{\star }$ and ${V}_{\mathrm{gas}}$ are the same as in Equation (4). The standard scale acceleration is ${a}_{0}\,={({H}_{0}/75)}^{2}\times 1.2\ {10}^{-8}=0.99\,\mathrm{cm}\,{{\rm{s}}}^{-2}$ (${H}_{0}=68$ $\mathrm{km}\,{{\rm{s}}}^{-1}$ Mpc⁻¹).

For the simple function

$$\begin{eqnarray}&&\mu (x)=\displaystyle \frac{x}{1+x}\end{eqnarray} \tag{ 12 }$$

that was shown to apply better to galaxy rotation curves (Famaey & Binney 2005), the velocity is given by

$$\begin{eqnarray}&&{V}_{\mathrm{SIM}}={V}_{\mathrm{lum}}\sqrt{1+0.5(\sqrt{1+4{a}_{0}R/{{\rm{V}}}_{\mathrm{lum}}^{2}}-1)}.\end{eqnarray} \tag{ 13 }$$

The ${V}_{\mathrm{STD}}$ and ${V}_{\mathrm{SIM}}$ models were both fit to the rotation of M33, with free ${a}_{0}$ and fixed or free ϒ.

We performed nonlinear Levenberg–Marquardt least-squares fits of the ISO, NFW, and MOND models. The ISO and NFW fits have two free parameters at fixed ϒ, three when ϒ is left free. The MOND fits have one free parameter at fixed ϒ, and two at free mass-to-light ratio. A normal weighting function set to ${\rm{\Delta }}{V}_{\mathrm{rot}}^{-2}$ is used. Figures 15 and 16 show the mass models and Table 5 lists the fitted parameters. A drawback of using different data sets to make a hybrid rotation curve is that the resulting distribution of velocity uncertainties is rarely Gaussian. We measure for instance a median Hα velocity uncertainty that is ∼60% that of the median H i velocity uncertainty. A consequence of the error non-homogeneity, combined with the nonuniform sampling of the rotation curve, is to yield high ${\chi }^{2}$ (reduced values ≳6), which, taken individually, means that a fit has no statistical significance despite the ∼350 degrees of freedom. We thus estimated the differences of ${\chi }^{2}$, relatively to a given model, that is, ${\rm{\Delta }}{\chi }^{2}=({\chi }_{i}^{2}-{\chi }_{j}^{2})/{\chi }_{i}^{2}$, to compare the models i and j with N free parameters, and find those that are more likely (see also Martinsson et al. 2013). The ${\rm{\Delta }}{\chi }^{2}$ values are reported in Table 6.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 5.  Results of the Mass Models

Model	Parameter	${\rm{\Upsilon }}=0.52$	${\rm{\Upsilon }}=0.72$	Free ϒ
ISO	${\rho }_{0}$	29.8 ± 1.1	10.6 ± 0.3	6.5 ± 0.6
	Rc	3.1 ± 0.1	6.1 ± 0.2	8.9 ± 0.7
	ϒ	0.52	0.72	0.80 ± 0.01
NFW	V200	114.4 ± 0.7	152.5 ± 0.4	139.9 ± 0.5
	c	6.05 ± 0.05	3.44 ± 0.01	4.55 ± 0.01
	ϒ	0.52	0.72	0.56 ± 0.01
MOND	a0	1.22 ± 0.01	0.78 ± 0.01	1.76 ± 0.11
Standard	ϒ	0.52	0.72	0.50 ± 0.03
MOND	a0	0.68 ± 0.01	0.33 ± 0.01	1.84 ± 0.11
Simple	ϒ	0.52	0.72	0.35 ± 0.02

Note. ϒ is in ${M}_{\odot }$/${L}_{\odot }$, ${\rho }_{0}$ and ${\rho }_{-2}$ in 10⁻³ ${M}_{\odot }$ pc⁻³, R_c and ${R}_{-2}$ in kpc, V₂₀₀ in $\mathrm{km}\,{{\rm{s}}}^{-1}$. For the MOND model, a₀ is in units of 10⁻⁸ cm s⁻². Uncertainties are statistical (formal) 1σ errors from the fits.

Download table as:  ASCIITypeset image

Table 6.  Comparison between Mass Models

Modeling with Dark Matter	Value	Diagnosis
$({\chi }_{\mathrm{ISO},{\rm{\Upsilon }}=0.52}^{2}-{\chi }_{\mathrm{NFW},{\rm{\Upsilon }}=0.52}^{2})/{\chi }_{\mathrm{NFW},{\rm{\Upsilon }}=0.52}^{2}$	0.52	NFW, ${\rm{\Upsilon }}=0.52$ more likely than ISO, ${\rm{\Upsilon }}=0.52$
$({\chi }_{\mathrm{ISO},{\rm{\Upsilon }}=0.72}^{2}-{\chi }_{\mathrm{NFW},{\rm{\Upsilon }}=0.72}^{2})/{\chi }_{\mathrm{NFW},{\rm{\Upsilon }}=0.72}^{2}$	−0.68	ISO, ${\rm{\Upsilon }}=0.72$ more likely than NFW, ${\rm{\Upsilon }}=0.72$
$({\chi }_{\mathrm{ISO},\mathrm{free}\ {\rm{\Upsilon }}}^{2}-{\chi }_{\mathrm{NFW},\mathrm{free}\ {\rm{\Upsilon }}}^{2})/{\chi }_{\mathrm{NFW},\mathrm{free}\ {\rm{\Upsilon }}}^{2}$	0.13	NFW, free ϒ more likely than ISO, free ϒ
$({\chi }_{\mathrm{NFW},{\rm{\Upsilon }}=0.72}^{2}-{\chi }_{\mathrm{NFW},{\rm{\Upsilon }}=0.52}^{2})/{\chi }_{\mathrm{NFW},{\rm{\Upsilon }}=0.52}^{2}$	0.52	NFW, ${\rm{\Upsilon }}=0.52$ more likely than NFW, ${\rm{\Upsilon }}=0.72$
$({\chi }_{\mathrm{ISO},{\rm{\Upsilon }}=0.72}^{2}-{\chi }_{\mathrm{ISO},{\rm{\Upsilon }}=0.52}^{2})/{\chi }_{\mathrm{ISO},{\rm{\Upsilon }}=0.52}^{2}$	−0.67	ISO, ${\rm{\Upsilon }}=0.72$ more likely than ISO, ${\rm{\Upsilon }}=0.52$
Modeling with MOND	Value	Diagnosis
$({\chi }_{\mathrm{SIM},{\rm{\Upsilon }}=0.52}^{2}-{\chi }_{\mathrm{STD},{\rm{\Upsilon }}=0.52}^{2})/{\chi }_{\mathrm{STD},{\rm{\Upsilon }}=0.52}^{2}$	0.43	STD, ${\rm{\Upsilon }}=0.52$ more likely than SIM, ${\rm{\Upsilon }}=0.52$
$({\chi }_{\mathrm{SIM},{\rm{\Upsilon }}=0.72}^{2}-{\chi }_{\mathrm{STD},{\rm{\Upsilon }}=0.72}^{2})/{\chi }_{\mathrm{STD},{\rm{\Upsilon }}=0.72}^{2}$	0.52	STD, ${\rm{\Upsilon }}=0.72$ more likely than SIM, ${\rm{\Upsilon }}=0.72$
$({\chi }_{\mathrm{SIM},\mathrm{free}\ {\rm{\Upsilon }}}^{2}-{\chi }_{\mathrm{STD},\mathrm{free}\ {\rm{\Upsilon }}}^{2})/{\chi }_{\mathrm{STD},\mathrm{free}\ {\rm{\Upsilon }}}^{2}$	0.11	STD, free ϒ more likely than SIM, free ϒ
$({\chi }_{\mathrm{STD},{\rm{\Upsilon }}=0.72}^{2}-{\chi }_{\mathrm{STD},{\rm{\Upsilon }}=0.52}^{2})/{\chi }_{\mathrm{STD},{\rm{\Upsilon }}=0.52}^{2}$	0.33	STD, ${\rm{\Upsilon }}=0.52$ more likely than STD, ${\rm{\Upsilon }}=0.72$
$({\chi }_{\mathrm{SIM},{\rm{\Upsilon }}=0.72}^{2}-{\chi }_{\mathrm{SIM},{\rm{\Upsilon }}=0.52}^{2})/{\chi }_{\mathrm{SIM},{\rm{\Upsilon }}=0.72}^{2}$	0.44	SIM, ${\rm{\Upsilon }}=0.52$ more likely than SIM, ${\rm{\Upsilon }}=0.72$

Note. STD and SIM are for MOND models with standard and simple (respectively) interpolation functions.

Download table as:  ASCIITypeset image

As in many studies of galaxy mass distribution from fitting of high-resolution rotation curves, the M33 mass models were not successful in reproducing all the irregularities of the rotation curve. This is explained by the impossibility of our axisymmetric modeling and inputs (spherical DM densities, planar gas and surface density profiles) to mimic the asymmetries in the disks (spiral arms, warp, lopsidedness) that reflect in wiggles in the axisymmetric rotation curve.

5.5.1. Results for Dark Matter Models

5.5.2. Results for MOND

All MOND models strongly favor the original standard interpolation function over the more simple one (Table 6). Figure 16 thus only presents results obtained with the standard interpolation function. Moreover, configurations with ${\rm{\Upsilon }}=0.52$ are preferred over ${\rm{\Upsilon }}=0.72$. This latter model obviously fails at reproducing the rotation curve from R = 9 kpc. This result is reflected in a best-fit mass-to-light ratio of 0.5 for the preferred standard μ-function, which corresponds to a stellar mass of $5.3\ {10}^{9}$ ${M}_{\odot }$. The scale acceleration ${a}_{0}=1.76\ {10}^{-8}$ cm s⁻² is by 78% higher than the standard value of ${a}_{0}=1.2\ {10}^{-8}\times {({H}_{0}/75)}^{2}\,\mathrm{cm}$ s⁻², for ${H}_{0}=68$ $\mathrm{km}\,{{\rm{s}}}^{-1}$ Mpc⁻¹.

The comparison between MOND and the ISO or NFW models cannot be made directly from ${\chi }^{2}$ residuals because these models have not the same number of free parameters. Instead, we made use of Akaike information criterion (AIC, Akaike 1974) because it is a simple linear combination of both the number of free parameters, N, and the fits quality, i.e., $\mathrm{AIC}=2N+{\chi }^{2}$. Consequently, comparing AIC residuals ${\rm{\Delta }}\mathrm{AIC}={\mathrm{AIC}}_{\mathrm{MOND}}\mbox{--}{\mathrm{AIC}}_{\mathrm{DM}}$ is well appropriate to determine which of MOND and ISO or NFW is more likely than the other (see Chemin et al. 2011, 2016). The AIC residuals are reported in Table 7 for the standard μ-function. The meaning of ${\rm{\Delta }}\mathrm{AIC}$ is similar to that of ${\rm{\Delta }}{\chi }^{2}$. Negative residuals mean that MOND is more likely than DM models, while positive residuals imply that MOND is less likely than DM models. It is found that at fixed mass-to-light ratio ${\rm{\Upsilon }}=0.52$, the NFW model is more appropriate than MOND, which is more likely than the ISO model. At ${\rm{\Upsilon }}=0.72$, the opposite result is found, the ISO model is more appropriate than MOND, which is more likely than the cusp. Another important point is that at free ϒ, both DM haloes are more appropriate than MOND. In other words, a more elaborate form than MOND, with a density law of hidden mass and an additional free parameter, has had a significant impact on the modeling of the mass distribution of M33. This result is even more significant for NFW than for the pseudo-isothermal sphere. This does not imply that MOND has to be rejected, however. It simply illustrates the preference for DM-based models to explain the mass distribution underlying the rotation curve of M33.

Table 7.  Comparison between Dark Matter and MOND Mass Models

MOND versus DM	Value	Diagnosis
${\mathrm{AIC}}_{\mathrm{STD},{\rm{\Upsilon }}=0.52}-{\mathrm{AIC}}_{\mathrm{NFW},{\rm{\Upsilon }}=0.52}$	7.3	NFW, ${\rm{\Upsilon }}=0.52$ more likely than STD, ${\rm{\Upsilon }}=0.52$
${\mathrm{AIC}}_{\mathrm{STD},{\rm{\Upsilon }}=0.72}-{\mathrm{AIC}}_{\mathrm{NFW},{\rm{\Upsilon }}=0.72}$	−2.2	STD, ${\rm{\Upsilon }}=0.72$ more likely than NFW, ${\rm{\Upsilon }}=0.72$
${\mathrm{AIC}}_{\mathrm{STD},\mathrm{free}{\rm{\Upsilon }}}-{\mathrm{AIC}}_{\mathrm{NFW},\mathrm{free}{\rm{\Upsilon }}}$	5.1	NFW, free ϒ more likely than STD, free ϒ
${\mathrm{AIC}}_{\mathrm{STD},{\rm{\Upsilon }}=0.52}-{\mathrm{AIC}}_{\mathrm{ISO},{\rm{\Upsilon }}=0.52}$	−16.5	STD, ${\rm{\Upsilon }}=0.52$ more likely than ISO, ${\rm{\Upsilon }}=0.52$
${\mathrm{AIC}}_{\mathrm{STD},{\rm{\Upsilon }}=0.72}-{\mathrm{AIC}}_{\mathrm{ISO},{\rm{\Upsilon }}=0.72}$	16.5	ISO, ${\rm{\Upsilon }}=0.72$ more likely than STD, ${\rm{\Upsilon }}=0.72$
${\mathrm{AIC}}_{\mathrm{STD},\mathrm{free}{\rm{\Upsilon }}}-{\mathrm{AIC}}_{\mathrm{ISO},\mathrm{free}{\rm{\Upsilon }}}$	1.9	ISO, free ϒ more likely than STD, free ϒ
${\mathrm{AIC}}_{\mathrm{STD},{\rm{\Upsilon }}=0.52}-{\mathrm{AIC}}_{\mathrm{NFW},{\rm{\Upsilon }}=0.72}$	−16.9	STD, ${\rm{\Upsilon }}=0.52$ more likely than NFW, ${\rm{\Upsilon }}=0.72$
${\mathrm{AIC}}_{\mathrm{STD},{\rm{\Upsilon }}=0.72}-{\mathrm{AIC}}_{\mathrm{NFW},{\rm{\Upsilon }}=0.52}$	21.9	NFW, ${\rm{\Upsilon }}=0.52$ more likely than STD, ${\rm{\Upsilon }}=0.72$
${\mathrm{AIC}}_{\mathrm{STD},{\rm{\Upsilon }}=0.52}-{\mathrm{AIC}}_{\mathrm{ISO},{\rm{\Upsilon }}=0.72}$	1.9	ISO, ${\rm{\Upsilon }}=0.72$ more likely than STD, ${\rm{\Upsilon }}=0.52$
${\mathrm{AIC}}_{\mathrm{STD},{\rm{\Upsilon }}=0.72}-{\mathrm{AIC}}_{\mathrm{ISO},{\rm{\Upsilon }}=0.52}$	−1.9	STD, ${\rm{\Upsilon }}=0.72$ more likely than ISO, ${\rm{\Upsilon }}=0.52$

Note. $\mathrm{AIC}=2N+{\chi }^{2}$ is the Akaike information criterion (Akaike 1974), where N is the number of parameters to fit. AIC residuals have been normalized to 100 for clarity. Results are for the standard (STD) interpolation function only, which is more likely than the simple μ-function.

Download table as:  ASCIITypeset image

The finding of a dependency of the inner shape of the DM halo on the mass-to-light ratio is in agreement with the analysis of Hague & Wilkinson (2015). These authors used a Bayesian approach to constrain the mass distribution of M33 from the kinematics of Corbelli & Salucci (2000). They excluded a combination of inner DM density slopes shallower than 0.9 with ${\rm{\Upsilon }}\lt 2$. Their modeling differs from ours, however, since they used an additional low-mass bulge component, and found that models with density slopes steeper than the NFW cusp, and with ${\rm{\Upsilon }}\sim 1.5$ are more likely. Such large ϒ can only be consistent with shallow DM density profiles with our higher-resolution data, which is ruled out by the present analysis, or by analysis of larger galaxy samples (e.g., Lelli et al. 2016).

At fixed stellar mass, finding a combination ${\rm{\Upsilon }}=0.52$-NFW halo as the most likely result is in very good agreement with the model of Corbelli et al. (2014). The corresponding halo concentration (c = 6.1) also agrees with the value given by these authors (c = 6.7). This concentration is nevertheless not consistent with c = 9, which is the value inferred from the Millenium Simulation halo mass–concentration relationship (Ludlow et al. 2014) at a redshift of z = 0 and for our virial mass ${M}_{200}$.

At free stellar mass, the inferred stellar mass is $5.9\ {10}^{9}$ ${M}_{\odot }$ (corresponding to a maximum velocity of 73 $\mathrm{km}\,{{\rm{s}}}^{-1}$), which is 23% higher than the most likely mass of Corbelli et al. (2014) (corresponding to a maximum velocity of 60 $\mathrm{km}\,{{\rm{s}}}^{-1}$). This difference is not significant, however, owing to the range of stellar mass predicted by their SPS models. The halo concentration c = 4.6 still disagrees with the halo-mass concentration relation, but also with the most likely concentration given in Corbelli et al. (2014), $c=9.5\pm 1.5$. The most likely values of Corbelli et al. (2014) were obtained by combining the probability density functions that best fit their rotation curve, a stellar mass compatible with SPS models, and a concentration compatible with the halo mass–concentration relation.

We verified that the concentration discrepancy is not caused by the choice of the adopted higher-resolution rotation curve by fitting NFW models at fixed or free stellar mass, and to various rotation curves. Such curves either combined our outer ($R\gt 6.5$ kpc) DRAO velocities with the inner ($R\lt 6.5$ kpc) curve from Corbelli et al. (2014), or our inner Hα velocities with the outer GBT points from Corbelli et al. (2014). All the fits yielded $c\lesssim 7.3$ and a stellar disk mass comparable to our estimate. We also performed a model at fixed concentration c = 9.5 and free mass-to-light ratio and found ${V}_{200}=99$ $\mathrm{km}\,{{\rm{s}}}^{-1}$ and ${\rm{\Upsilon }}=0.32$. This corresponds to a maximum velocity of 55 $\mathrm{km}\,{{\rm{s}}}^{-1}$ for the stellar contribution. It is another way to illustrate the halo concentration-stellar mass degeneracy shown in Section 5.5.1 and in Hague & Wilkinson (2015).

To summarize, the origin of the concentration difference between both studies can only be the composite likelihood assumption made in Corbelli et al. (2014). Our higher-resolution data set and the adopted cusp model not tied to the halo mass–concentration relation strongly favors a concentration in disagreement with ΛCDM simulations. Since the stellar mass is the key parameter that allows the concentration to match the NFW cusp with the central mass distribution imposed by the rotation curve, the cusp concentration contradiction in M33 can be solved only when the stellar mass is known with better accuracy.