HI Kinematics and Mass Distribution of Messier 33

TL;DR

This study presents a detailed HI survey of Messier 33, revealing complex gas kinematics, a warped disk, and a cuspy dark matter halo, challenging some predictions of the LambdaCDM model.

Contribution

It provides new deep HI observations of M33, identifying a low brightness anomalous HI component and modeling its mass distribution with implications for dark matter halo profiles.

Findings

Detection of a perturbed outer gas distribution and kinematic twist.

Identification of a low brightness HI component with anomalous velocities.

Mass modeling favors a cuspy dark matter halo with a high concentration.

Abstract

A new deep HI survey of the galaxy Messier 33 is presented, based on observations obtained at the Dominion Radio Astrophysical Observatory. We observe a perturbed outer gas distribution and kinematics in M33, and confirm the disk warping as a significant kinematical twist of the major axis of the velocity field, though no strong tilt is measured, in agreement with previous work. Evidence for a new low brightness HI component with anomalous velocity is reported. It harbours a large velocity scatter, as its kinematics both exceeds and lags the rotation of the disk, and leaks in the forbidden velocity zone of apparent counter-rotation. The observations also reveal wide and multiple peak HI profiles which can be partly explained by crowded orbits in the framework of the warp model. Asymmetric motions are identified in the velocity field, as possible signatures of a lopsided potential and the warp. The mass distribution modeling of the hybrid Halpha-HI rotation curve favours a cuspy dark matter halo with a concentration in disagreement with the LambdaCDM dark halo mass-concentration relationship. The total mass enclosed in 23 kpc is 8 10^10 Msol, of which 11% are stars and gas. At the virial radius of the cuspy halo, the resulting total mass is 5 10^11 Msol, but with a baryonic mass fraction of 2% only. This strongly suggests a more realistic radius encompassing the total mass of M33 well smaller than the virial radius of the halo, maybe comparable to the size of the HI disk.

Tables

Table 1. Table 1: Parameters of M 33.

Parameters	Value	Source
Morphological type	SA(s)cd	RC3
R.A. (2000)	01h 33m 33.1s	RC3
Dec. (2000)	+30° 39${}^{{}^{\prime }}$ 18′′	RC3
Systemic Velocity (km s-1)	$-179\pm 3$	RC3
Distance (Mpc)	0.84
Scale (pc/arcmin)	244
Disk Scale length (kpc)	1.6 (@ 3.6 $\mu$m)	Kam15
Optical radius, $R_{25}$	35$\stackrel{\prime }{.}$4 $\pm$ 1$\stackrel{\prime }{.}$0	RC3
Inclination, $i$	52°  $\pm$ 3°	WWB
Position angle (major axis)	22.5° $\pm$ 1°	WWB
Apparent magnitude, mV	5.28	RC3
Absolute magnitude, MV	$-19.34$
Total Hi mass (M⊙)	${\mathrm{1.95\hspace{0.17em}10}}^9$	Sec. 2
Systemic Velocity (km s-1)	$-180\pm 3$	Sec. 2
$V_{\mathrm{rot}}$ maximum (km s-1)	125	Sec. 4
Stellar mass (M⊙), mass models	${\mathrm{5.5\hspace{0.17em}10}}^9$	Sec. 5
Dynamical mass (M⊙), inside $R=23$ kpc	${\mathrm{7.9\hspace{0.17em}10}}^{10}$	Sec. 5

Table 2. Table 2: Summary of the six 21 cm H i   line synthesis fields centred on and surrounding M33, carried out with the DRAO Synthesis Telescope.

Observ.	Field Centre	Beam Parameters
Date	(RA, DEC) (J2000.0)	${\theta }_{\mathrm{maj}}{(}^{\prime })\times {\theta }_{\min }{(}^{\prime })$, CCWE
09/29/08	01${}^{\text{h}}$33${}^{\text{m}}$50.9${}^{\text{s}}$, +30$°$39$\prime$36$″$	$1.90\times 0.97\prime$,$-$89.69$°$
09/29/08	01${}^{\text{h}}$36${}^{\text{m}}$10.2${}^{\text{s}}$, +31$°$50$\prime$34$″$	$1.85\times 0.97$, $-$89.82$°$
11/05/08	01${}^{\text{h}}$31${}^{\text{m}}$38.2${}^{\text{s}}$, +29$°$28$\prime$12$″$	$1.98\times 0.97$, $-$89.91$°$
11/05/08	01${}^{\text{h}}$34${}^{\text{m}}$45.8${}^{\text{s}}$, +31$°$08$\prime$13$″$	$1.86\times 0.97$, $-$90.11$°$
12/04/08	01${}^{\text{h}}$32${}^{\text{m}}$56.4${}^{\text{s}}$, +30$°$11$\prime$01$″$	$1.94\times 0.97$, $-$90.30$°$
12/04/08	01${}^{\text{h}}$33${}^{\text{m}}$50.9${}^{\text{s}}$, +30$°$39$\prime$36$″$	1$\prime$.91$\times$0$\prime$.97, $-$90.49$°$

Table 3. Table 3: Coordinates of five representative spectra exhibiting a H i  velocity component in the forbidden velocity zones from the PV diagrams of Fig.  4 .

PV diagram	Offset	Coordinates (RA, DEC) (J2000)
PA=165°	$+7\prime$	01${}^{\text{h}}$33${}^{\text{m}}$40.7${}^{\text{s}}$, +30$°$31$\prime$56$″$
PA=165°	$-11\prime$	01${}^{\text{h}}$33${}^{\text{m}}$22.0${}^{\text{s}}$, +30$°$49$\prime$46$″$
PA=202°	$+10\prime$	01${}^{\text{h}}$33${}^{\text{m}}$17.1${}^{\text{s}}$, +30$°$29$\prime$43$″$
PA=202°	$-12\prime$	01${}^{\text{h}}$33${}^{\text{m}}$52.6${}^{\text{s}}$, +30$°$50$\prime$09$″$
PA=202°	$-30\prime$	01${}^{\text{h}}$34${}^{\text{m}}$26.7${}^{\text{s}}$, +31$°$07$\prime$49$″$

Table 4. Table 4: Results of the tilted-ring model of the H i  velocity field of M33.

Radius	Radius	$V_{\mathrm{rot}}$	$\Delta V_{\mathrm{rot}}$	$i$	PA	Radius	Radius	$V_{\mathrm{rot}}$	$\Delta V_{\mathrm{rot}}$	$i$	PA
(′)	(kpc)	(km s-1)	(km s-1)	(°)	(°)	(′)	(kpc)	(km s-1)	(km 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

Table 5. Table 5: Results of the mass models.

Model	Parameter	${\rm Y}=0.52$	${\rm Y}=0.72$	Free ${\rm Y}$
ISO	${\rho }_0$	$29.8\pm 1.1$	$10.6\pm 0.3$	$6.5\pm 0.6$
	$R_c$	$3.1\pm 0.1$	$6.1\pm 0.2$	$8.9\pm 0.7$
	${\rm Y}$	0.52	0.72	$0.80\pm 0.01$
NFW	$V_{200}$	$114.4\pm 0.7$	$152.5\pm 0.4$	$139.9\pm 0.5$
	$c$	$6.05\pm 0.05$	$3.44\pm 0.01$	$4.55\pm 0.01$
	${\rm Y}$	0.52	0.72	$0.56\pm 0.01$
MOND	$a_0$	$1.22\pm 0.01$	$0.78\pm 0.01$	$1.76\pm 0.11$
Standard	${\rm Y}$	0.52	0.72	$0.50\pm 0.03$
MOND	$a_0$	$0.68\pm 0.01$	$0.33\pm 0.01$	$1.84\pm 0.11$
Simple	${\rm Y}$	0.52	0.72	$0.35\pm 0.02$

Table 6. Table 6: Comparison between mass models.

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

Table 7. Table 7: Comparison between Dark Matter and MOND mass models.

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

Table 8. Table 8: H i  dispersion and mass surface density of M33

Radius	${\sigma }_{\mathrm{los}}$	${\Sigma }_{H\text{i}}$	Radius	${\sigma }_{\mathrm{los}}$	${\Sigma }_{H\text{i}}$	Radius	${\sigma }_{\mathrm{los}}$	${\Sigma }_{H\text{i}}$	Radius	${\sigma }_{\mathrm{los}}$	${\Sigma }_{H\text{i}}$
2	0.7	7.89	26	7.6	7.54	50	8.5	0.49	74	8.2	0.03
4	9.5	7.73	28	7.6	6.50	52	8.2	0.36	76	6.5	0.04
6	9.3	8.26	30	8.0	4.99	54	8.3	0.23	78	6.5	0.04
8	9.6	8.52	32	8.9	3.50	56	7.6	0.18	80	7.4	0.04
10	0.0	8.23	34	9.8	2.51	58	7.7	0.12	82	6.5	0.08
12	0.1	7.66	36	0.2	1.75	60	6.8	0.08	84	5.6	0.09
14	9.1	8.00	38	0.3	1.34	62	7.4	0.04	86	6.1	$\cdots$
16	8.2	8.61	40	0.5	1.00	64	7.5	0.03	88	5.7	$\cdots$
18	8.0	8.14	42	9.5	0.75	66	6.7	0.03	90	5.7	0.09
20	8.4	8.19	44	8.8	0.67	68	6.6	0.00	92	5.6	$\cdots$
22	7.9	7.61	46	9.1	0.56	70	6.4	0.03	94	6.4	$\cdots$
24	8.0	7.57	48	9.0	0.63	72	8.0	0.07	96	7.5	$\cdots$

Table 9. Table 9: Hybrid H α 𝛼 \alpha - H i  rotation curve of M33

$\ln \left(R\right)$	$V_{\mathrm{rot}}$	$\Delta V_{\mathrm{rot}}$	$\ln \left(R\right)$	$V_{\mathrm{rot}}$	$\Delta V_{\mathrm{rot}}$	$\ln \left(R\right)$	$V_{\mathrm{rot}}$	$\Delta V_{\mathrm{rot}}$	$\ln \left(R\right)$	$V_{\mathrm{rot}}$	$\Delta V_{\mathrm{rot}}$	$\ln \left(R\right)$	$V_{\mathrm{rot}}$	$\Delta V_{\mathrm{rot}}$
1.61	6.20	2.32	5.87	71.95	4.15	6.56	85.63	0.99	6.96	97.13	0.56	7.25	108.25	6.96
2.30	9.30	2.71	5.89	72.19	4.56	6.57	84.64	0.62	6.97	100.87	1.88	7.25	107.50	7.89
2.71	13.00	0.55	5.90	72.40	5.20	6.57	87.03	0.66	6.97	100.32	2.41	7.25	108.45	5.87
3.00	18.25	3.79	5.91	71.80	5.68	6.58	89.26	1.30	6.98	97.17	0.57	7.26	108.88	5.48
3.22	20.73	3.31	5.93	71.14	4.69	6.59	90.95	1.33	6.98	98.41	1.45	7.26	107.92	5.56
3.40	22.69	3.50	5.94	71.81	4.74	6.59	92.36	1.26	6.98	98.67	0.55	7.27	105.22	6.43
3.56	24.54	1.37	5.95	71.89	7.37	6.60	92.75	0.66	6.99	101.05	0.62	7.27	105.77	6.42
3.69	26.23	0.68	5.97	72.47	7.82	6.61	91.76	2.36	6.99	101.51	0.55	7.27	110.49	4.32
3.81	30.62	1.15	5.98	73.25	7.98	6.61	91.09	3.59	7.00	102.93	1.20	7.28	110.51	4.05
3.91	31.84	0.94	5.99	73.46	8.12	6.62	91.89	2.46	7.00	103.97	1.17	7.28	105.47	5.25
4.01	34.51	1.18	6.00	74.46	8.06	6.63	91.47	2.86	7.01	103.50	1.60	7.28	105.33	5.86
4.09	35.42	1.46	6.02	75.18	5.68	6.63	92.49	1.66	7.01	104.44	1.15	7.29	109.38	1.96
4.17	37.37	2.00	6.03	76.20	4.64	6.64	92.69	1.02	7.02	104.81	1.10	7.29	113.28	4.31
4.25	38.38	2.49	6.04	77.74	3.39	6.65	93.88	0.66	7.02	104.98	1.67	7.29	115.69	2.75
4.32	39.20	1.20	6.05	78.30	2.56	6.65	95.00	0.86	7.03	105.49	1.90	7.30	116.42	2.97
4.38	41.82	2.47	6.06	78.07	2.50	6.66	94.41	0.56	7.03	106.07	1.72	7.30	115.31	3.52
4.44	41.08	2.43	6.08	77.82	3.33	6.67	93.85	2.13	7.03	106.00	0.98	7.30	112.56	3.18
4.50	41.74	1.71	6.09	78.52	2.93	6.67	93.11	3.81	7.04	107.46	0.76	7.31	111.32	2.71
4.55	42.57	2.18	6.10	78.85	2.04	6.68	94.65	2.99	7.04	109.03	1.80	7.31	112.61	2.99
4.61	43.40	4.11	6.11	78.78	1.12	6.68	96.50	2.17	7.05	109.25	2.05	7.31	113.70	2.34
4.65	43.85	3.64	6.12	78.14	0.57	6.69	96.79	1.76	7.05	108.48	0.71	7.32	114.91	0.89
4.70	45.95	3.03	6.13	79.82	1.28	6.70	96.32	0.60	7.06	107.45	1.20	7.32	115.79	5.15
4.74	48.40	0.96	6.14	79.36	1.60	6.70	96.00	1.00	7.06	107.63	1.50	7.32	112.92	9.94
4.79	49.04	2.33	6.15	79.58	1.10	6.71	95.99	0.55	7.06	107.73	1.95	7.33	119.07	5.09
4.83	49.37	4.79	6.16	79.11	0.55	6.72	94.58	0.82	7.07	108.36	0.88	7.33	111.82	13.33
4.87	49.58	5.19	6.17	79.62	0.88	6.72	95.24	0.56	7.07	108.79	0.61	7.33	113.09	13.94
4.91	50.83	3.94	6.18	82.16	2.31	6.73	95.08	0.57	7.08	108.59	0.70	7.34	110.64	13.95
4.94	52.27	4.79	6.19	82.16	2.24	6.73	96.30	2.29	7.08	108.83	0.71	7.34	108.71	15.32
4.98	52.64	6.30	6.20	82.59	3.54	6.74	96.56	1.87	7.09	104.87	3.58	7.34	106.90	18.46
5.01	52.61	7.27	6.21	82.91	4.06	6.75	96.82	0.55	7.09	103.69	3.97	7.35	106.47	17.74
5.04	52.97	6.74	6.22	82.02	3.83	6.75	98.17	0.70	7.09	106.19	1.31	7.35	107.07	17.69
5.08	54.45	5.94	6.23	81.00	4.64	6.76	97.69	1.36	7.10	107.36	0.62	7.35	107.73	19.17
5.11	54.75	6.23	6.24	81.58	3.35	6.76	96.95	0.65	7.10	107.99	1.15	7.36	110.16	17.09
5.14	53.44	5.79	6.25	82.93	4.11	6.77	95.99	0.64	7.11	108.06	0.70	7.36	109.82	17.98
5.16	55.63	4.15	6.26	83.81	4.91	6.77	95.19	0.90	7.11	107.26	1.15	7.36	111.36	19.37
5.19	55.99	1.39	6.27	83.69	5.27	6.78	94.76	2.86	7.11	107.30	1.04	7.37	108.01	16.48
5.22	54.21	2.04	6.28	82.93	4.98	6.79	94.75	1.55	7.12	105.96	1.11	7.37	104.27	15.46
5.25	52.71	3.09	6.29	81.05	2.57	6.79	94.29	2.72	7.12	106.36	0.91	7.37	104.57	16.63
5.27	52.96	2.03	6.30	80.03	0.96	6.80	95.88	4.38	7.13	107.76	1.47	7.37	107.86	17.65
5.30	55.06	0.83	6.31	79.84	1.28	6.80	98.05	6.72	7.13	106.64	0.70	7.43	106.76	2.16
5.32	55.61	1.61	6.32	78.11	0.64	6.81	98.67	6.77	7.13	105.83	1.00	7.50	107.33	3.02
5.35	56.01	1.60	6.33	78.30	0.80	6.81	99.73	8.12	7.14	107.77	1.97	7.56	108.29	4.01
5.37	56.37	1.47	6.34	79.95	1.11	6.82	100.12	8.06	7.14	109.46	0.58	7.62	109.72	4.01
5.39	57.87	0.86	6.35	81.21	2.43	6.82	101.28	8.81	7.15	111.04	0.96	7.68	111.98	4.78
5.42	58.99	0.80	6.35	82.36	4.01	6.83	100.50	8.65	7.15	105.97	2.53	7.73	116.06	2.18
5.44	59.17	1.33	6.36	81.44	3.53	6.84	100.15	8.55	7.15	104.96	5.49	7.78	117.23	2.45
5.46	60.13	2.08	6.37	80.95	4.15	6.84	102.06	8.80	7.16	111.44	8.30	7.83	116.46	6.48
5.48	60.68	1.47	6.38	81.85	3.98	6.85	98.96	5.91	7.16	112.71	7.21	7.88	115.68	8.07
5.50	61.23	2.09	6.39	81.83	2.81	6.85	96.91	3.30	7.17	112.66	7.25	7.92	117.40	8.23
5.52	61.59	1.58	6.40	81.39	3.26	6.86	97.79	4.86	7.17	99.94	4.97	7.97	116.84	8.93
5.54	61.26	1.64	6.41	81.60	3.75	6.86	97.54	4.49	7.17	95.60	9.12	8.01	115.70	9.64
5.56	62.38	2.57	6.41	83.01	1.77	6.87	96.80	3.84	7.18	97.16	9.68	8.05	115.09	7.69
5.58	63.79	3.62	6.42	82.73	1.34	6.87	96.34	3.32	7.18	97.93	4.86	8.08	117.09	5.11
5.60	64.71	4.29	6.43	82.19	2.18	6.88	96.61	2.97	7.19	102.58	5.33	8.12	118.20	3.15
5.62	64.44	4.18	6.44	83.79	0.84	6.88	96.58	2.36	7.19	103.68	4.01	8.15	118.42	1.42
5.63	64.76	5.09	6.45	84.23	2.02	6.89	97.27	3.13	7.19	104.75	4.99	8.19	118.20	1.78
5.65	66.67	4.51	6.45	85.03	4.17	6.89	96.27	3.97	7.20	106.04	7.02	8.22	117.45	2.38
5.67	68.62	3.97	6.46	86.46	2.47	6.90	97.57	4.77	7.20	106.78	5.93	8.25	119.56	0.79
5.69	68.09	3.94	6.47	86.55	1.43	6.90	99.50	4.53	7.20	105.56	6.01	8.28	118.60	1.53
5.70	65.49	4.83	6.48	86.27	0.82	6.91	103.83	7.82	7.21	104.72	4.71	8.31	122.63	0.50
5.72	65.61	5.44	6.48	86.51	0.77	6.91	103.79	7.01	7.21	106.75	4.59	8.34	124.10	2.86
5.74	66.85	4.35	6.49	86.51	1.19	6.92	103.69	5.59	7.22	107.06	5.78	8.37	125.03	2.19
5.75	67.49	5.55	6.50	86.73	2.25	6.92	102.16	4.12	7.22	104.89	6.04	8.40	125.49	2.54
5.77	67.81	5.14	6.51	85.88	1.65	6.93	99.43	3.12	7.22	106.56	4.29	8.42	125.24	8.08
5.78	68.85	4.03	6.51	85.95	0.60	6.93	99.70	2.12	7.23	107.95	4.16	8.45	121.96	9.76
5.80	68.92	4.30	6.52	85.53	0.55	6.94	103.82	6.87	7.23	109.25	3.79	8.48	120.39	8.47
5.81	70.75	3.98	6.53	85.95	2.42	6.94	98.75	0.60	7.23	110.55	3.30	8.50	114.04	26.64
5.83	71.41	3.79	6.54	85.40	1.30	6.95	99.30	1.32	7.24	111.55	3.30	8.52	110.03	34.64
5.84	71.66	4.62	6.54	82.53	0.55	6.95	99.55	1.55	7.24	110.80	1.73	8.55	98.69	27.45
5.86	72.38	4.01	6.55	85.17	0.67	6.96	104.52	5.50	7.24	110.85	2.90	8.57	100.07	33.42
8.59	104.32	35.17	8.62	101.19	27.35	8.64	123.49	39.13

Equations

$${\mathrm{V}}_{\mathrm{obs}}={\mathrm{V}}_{\mathrm{sys}}+{\mathrm{V}}_{\mathrm{rot}}\cos \theta \sin i$$

$${\mathrm{V}}_{\mathrm{obs}}=c_0+\sum _{\mathrm{k}=1}\left(c_k\cos k\theta +s_k\sin k\theta \right)\sin i$$

$$\Sigma [{\mathrm{M}}_{\odot }\mathrm{p}{\mathrm{c}}^{-2}]={\rm Y}\times 10^{-0.4({\mu }_{3.6}-C_{3.6})}$$

$${\mathrm{V}}_{\mathrm{rot}}^2={\mathrm{V}}_{\star }^2+{\mathrm{V}}_{\mathrm{gas}}^2+{\mathrm{V}}_{\mathrm{DM}}^2$$

$${\rho }_{\mathrm{ISO}}(R)=\frac{{\rho }_0}{1+(R/R_c{)}^2}$$

$${\mathrm{V}}_{\mathrm{ISO}}(R)=\sqrt{4\pi G{\rho }_0{R}_c^2\left(1-R/R_c\arctan (R/R_c)\right)}$$

$${\rho }_{\mathrm{NFW}}(R)=\frac{{\rho }_i}{R/R_s(1+R/R_s{)}^2}$$

$${\mathrm{V}}_{\mathrm{NFW}}(R)={\mathrm{V}}_{200}\sqrt{\frac{\ln (1+\mathrm{cx})-cx/(1+\mathrm{cx})}{\mathrm{x}(\ln (1+\mathrm{c})-c/(1+\mathrm{c}))}}$$

$$\mu (x=g/a_0)g=g_N$$

$$\mu (x)=\frac{x}{\sqrt{1+x^2}}$$

$${\mathrm{V}}_{\mathrm{STD}}={\mathrm{V}}_{\mathrm{lum}}\sqrt{\sqrt{0.5\left(1+\sqrt{1+(2{\mathrm{a}}_{0}R/{\mathrm{V}}_{\mathrm{lum}}^2{)}^2}\right)}}$$

$$\mu (x)=\frac{x}{1+x}$$

$${\mathrm{V}}_{\mathrm{SIM}}={\mathrm{V}}_{\mathrm{lum}}\sqrt{1+0.5\left(\sqrt{1+4{\mathrm{a}}_{0}R/{\mathrm{V}}_{\mathrm{lum}}^2}-1\right)}$$

Full text

Local Group – techniques: interferometric – galaxies: kinematics and dynamics – galaxies: individual: M33

††slugcomment: Accepted for publication in the Astronomical Journal

1 Introduction

With the Milky Way and Andromeda, the Triangulum galaxy (Messier 33, the third most massive disk galaxy of the Local Group) has long been among the most studied nearby galaxies to scrutinize the chemical, dynamical and structural properties of the stellar populations and of the interstellar medium. In particular, as it is a prototype of gas-rich spirals of moderate inclination (Table 1), it is very appropriate to study the relationships of the atomic, molecular and ionized gas content with star formation inside the disk. Since it is now well admitted that M33 has undergone a tidal encounter with his massive companion M31, as shown by the perturbation of its close environment (e.g. Braun & Thilker, 2004; McConnachie et al., 2010; Wolfe et al., 2013), it is expected that gas expelled during that interaction is currently returning into the M33 disk, fueling the active star formation (Putman et al., 2009).

The implied important population of Hii regions (e.g. Boulesteix et al., 1974; Zaritsky et al., 1989; Relaño et al., 2013) has thus been a motivation for us to present the first large-scale, arcsec-resolution, 3D spectroscopy survey of the disk of M33 in the H$\alpha$ emission line (Kam et al., 2015, hereafter Kam15). On one hand, the objectives of this survey were to measure the internal kinematics of star forming regions. For instance, Kam15 presented detailed velocity fields of compact and extremely large Hii regions, like NGC 604, or the relation between the velocity dispersion and the integrated intensity of the H$\alpha$ line, underlying the physical processes occurring in these regions and in the diffuse interstellar medium (stellar winds, expansion, etc.). The catalog of Hii regions to be provided from this survey will also be useful to study the relationships with other tracers of star formation in M33 at an unprecedented level of details. On the other hand, the H$\alpha$ mapping has been a unique opportunity to determine for the first time the most extended H$\alpha$ velocity field of M33. The modeling of the velocity field led Kam15 to conclude that the kinematical parameters of the ionized gas disk are very consistent with those of the stellar disk. The most important result of Kam15 is the determination of the H$\alpha$ rotation curve out to 8 kpc sampled every 20 pc, which is, as of date, the most resolved rotation curve obtained for any massive spiral galaxies other than the Milky Way. That curve perfectly traces the velocity gradient in the inner disk, and presents many wiggles characteristic of spiral arms perturbations, barely seen in previous Hi and CO observations of M33. These irregularities however do not prevent the H$\alpha$, Hi or CO rotation velocities to remain in good agreement.

The derivation of the H$\alpha$ rotation curve actually constitutes the first pillar of a broader project devoted to revisit the modeling of the mass distribution of M33. The other pillar inherent to such modeling is to get the Hi rotation curve to cover as far as possible the outer disk, in order to obtain accurate fits of the dark matter (DM) distribution, or of alternate gravity models such as Modified Newtonian Dynamics (MOND, Milgrom, 1983a, b).

Existing Hi studies of M33 attest to which extent the tidal interaction with Andromeda has not been without consequences on the gas distribution in the disk outskirts, revealing perturbed features like the strong warp, arc-like structures or diffuse and discrete gas around the disk (Corbelli & Schneider, 1997; Putman et al., 2009; Lockman et al., 2012). Recently, Corbelli et al. (2014) combined VLA and GBT observations to model the Hi warp and rotation curve. They have presented a new rotation curve that extends about 5 kpc further out than previous studies. They have also shown that the distribution of dark matter is consistent with a model for which the mass density steeply decreases in the centre of the halo (the cosmological cusp à la Navarro-Frenk-White, see Navarro et al., 1997), and whose concentration agrees well with the halo mass-concentration relation from $\Lambda$ Cold Dark Matter simulations (Ludlow et al., 2014). Comparable results have been obtained by Hague & Wilkinson (2015), but from a lower resolution Hi  curve derived earlier by Corbelli & Salucci (2000). They concluded that models with density profiles with an inner slope shallower than 0.9 (as measured at $R\sim 0.5$ kpc, the first velocity point of their rotation curve) are only compatible with 3.6-$\rm\mu m$ mass-to-light ratios $\Upsilon>2$, while cuspier densities can coexist with $\Upsilon<2$. However, a halo with a constant density core in M33 would be difficult to reconcile with stellar populations models, or with other observational large scale dynamical studies (e.g. Lelli, McGaugh & Schombert, 2016). This is something we want to revisit with a new set of data.

In this context, we have performed an Hi survey of M33 at the Dominion Radio Astrophysical Observatory (DRAO). This arcminute-resolution survey is a good intermediate between VLA and Arecibo or GBT measurements. The objectives of this article are first to present the survey, the gas content and distribution of M33, a new tilted-ring model of the Hi velocity field and the Hi rotation curve. We also want to examine the shape and amplitude of the outer rotation curve, at radii beyond $R=17$ kpc where Corbelli et al. (2014) presented new velocities. The second objective of the article is to benefit from the high-resolution survey of Kam15 and perform mass distribution models from an hybrid resolution H$\alpha$-Hi rotation curve. In particular we want to determine the most appropriate stellar mass, density profile of dark matter, and infer the total mass of M33. In addition, we want to compare dark matter mass models with Modified Newtonian Dynamics.

Throughout the article, we adopt a Hubble constant of 68 km s*-1* Mpc*-1* (Planck Collaboration, 2016) and a distance to M33 of 0.84 Mpc (see Kam15 and references therein). The basic parameters of M33 are summarized in Table 1. The observations and reduction of the new DRAO data are presented in Section 2, which also gives general characteristics of a combined DRAO+Arecibo Hi datacube. Section 4 presents the analysis of the Hi distribution, the tilted-ring and Fourier models of the Hi velocity field. It also determines a hybrid H$\alpha$-Hi rotation curve to be used for the modeling of the mass distribution performed in Section 5.

2 21-cm observations and data reduction

The primary observations for this study were made with the Synthesis Telescope (ST) at the Dominion Radio Astrophysical Observatory (DRAO). This telescope is an East-West interferometer consisting of seven ($\sim$9 m diameter) dishes spaced variously across a baseline range of 13 to 617 m. At 1420 MHz the longest baseline achieves a synthesized half-power beamwidth of 49$\arcsec$(EW) by 49$\arcsec$/sin $\delta$(NS) with uniform weighting, although in the Hi  line we use a Gaussian taper in the $u,v$ plane to increase the sensitivity of each velocity channel at the slight expense of resolution (58$\arcsec\times$58$\arcsec/\textrm{sin}\delta$). As the goal of this study is to trace the extended rotation curve of the galaxy, the DRAO instrument was chosen because of two inherent strengths: its wide field (3.1$\arcdeg$), combined with its deep integration (144 h per field). The DRAO instrument is also unmatched in its absolutely calibrated polarization capability (see Landecker et al., 2000; Kothes et al., 2010, for specifications of the DRAO telescope). Observations of M33 consist of Stokes I,Q,U & V made in a 30 MHz continuum band centred at 1420 MHz ($\lambda$21 cm), Stokes I in a 2 MHz band at 408 MHz ($\lambda$74 cm), and 256 channels within a 2 MHz band centred on the Hi  line. For this study we use only the Hi  data products; future studies will present the total power and polarized radio emission from M33. The datacube velocity resolution is 2.64 km s*-1* at 21-cm and each channel is $\rm\Delta V=1.65$ km s*-1* wide. The band was centred on a heliocentric velocity of $\rm V_{\textsc{hel}}=-180$ km s*-1*. Data processing and mosaic-making steps proceeded as per those in Chemin et al. (2009) for the M31 observations with the DRAO telescope. The measured noise per channel is $\sim 12-13$ mJy beam*-1* or $\Delta T_{B}\sim 1.1$ K in the elliptical synthesized beam (58″ $\times$ 114″). To increase sensitivity to the faint outermost Hi  disk of M33 a total of 6 full-synthesis pointings on and around the galaxy were observed and mosaiced together (see Table 2).

To recover the large-scale Hi  structures in M33, we merge single-dish (aka “short-spacing”) data obtained from the Turn-On GALFA Survey (TOGS) portion of the GALFA-Hi  survey at Arecibo. These data, previously published by Putman et al. (2009), have an angular resolution of 3$\farcm$4 and a velocity resolution of 5.15 km s*-1*, and provide exceptional spatial frequency overlap with the DRAO synthesis data (which are missing structures larger than about 45$\arcmin$). Most importantly, TOGS data are fully corrected for stray-radiation entering the sidelobes, thus preventing contamination of particularly the M33 outskirts and allowing the extended gas distribution and kinematics to be traced. TOGS Hi  data are added to the calibrated interferometer-only mosaic in the same manner as in Chemin et al. (2009) after converting them to the same spatial and velocity-resolution and grid. Figure 1 compares the full resolution Hi gas distribution to the WISE W1 (3.4 $\mu$m) and W3 (12 $\mu$m) images in which the foreground stars have been identified and removed (Jarrett et al., in prep).

To gain sensitivity, we made two final datacubes smoothed to velocity resolutions of 5.3 and 10.6 km s*-1*, and spatially smoothed to a circular 120$\arcsec\times$120$\arcsec$ Gaussian beam. The measured 1$\sigma$ noise at the pointing centre of the final short-spacings-added 2$\arcmin$ spatial resolution and 10 km s*-1* velocity resolution mosaic is 2 mJy beam*-1*, or $\sim 80$ mK channel*-1*.

While our final Hi  cube is as sensitive on a per-beam basis as the VLA BCD data presented in Gratier et al. (2010) and in Corbelli et al. (2014), the 120$\arcsec$ and 10 km s*-1* resolution mosaic certainly lacks the fine detail and velocity resolution of the current state-of-the-art mapping of M33. Rather, it is intended to trace the extended Hi  disk out to very large radii, and large Hi structures. Indeed, since the typical cloud-cloud velocity dispersion in the halo of the Milky Way is 25 km s*-1* (de Heij et al., 2002), our 10 km s*-1* velocity resolution cube is more than sufficient to fully sample similar Hi clouds in the M33 disk outskirts.

3 General properties of the atomic neutral gas in M33

3.1 HI channel maps, profile and mass

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

The Hi integrated profile shown in Fig. 3 is asymmetric, with slightly more gas in the receding Southern half (3% relatively to the approaching Northern half). The intensity weighted systemic velocity is $-180.3\pm 2.3$ km s*-1*, as derived from any channels with velocities $<-50$ km s*-1*. We define the maximum flux of the global Hi profile as the average value of the two maxima. The width of the Hi profile at 50% of the maximum flux is $W_{50}=183$ km s*-1*, with a corresponding mid-point velocity of $-180.4$ km s*-1*, similar to the intensity weighted mean value. At the 20% level, the velocity width is $W_{20}=200$ km s*-1*. The total Hi mass is $1.95\pm 0.36\,10^{9}$ M⊙. This value is in agreement with the mass found by Corbelli et al. (2014) from VLA-GBT data, and unsurprisingly with the one derived by Putman et al. (2009). The total neutral gas mass is about six times larger than the molecular gas mass of $\sim 3.3\times 10^{8}$ M⊙ given in (Gratier et al., 2010).

3.2 Anomalous velocity gas in M33

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 §4.2), one along the major axis of the outer warped disk (PA = 165°), both of them being centered on the photometric centre, and another one along the direction PA=175°, but slightly off-centered from the photometric centre (bottom panel). The center of this last PV diagram is $(\alpha,\delta)_{\rm J2000}=(01{\rm h}33{\rm m}45.7{\rm s},+30\arcdeg 42^{\prime}51^{\prime\prime})$. This slice orientation is chosen because it goes through two regions of larger velocity dispersion at the NNW and SSE (see §4.1). The width of each PV slice is equivalent to two pixels ($\sim 45\arcsec$). 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.

Along the direction $\rm PA=202\arcdeg$, the high brightness gas is typical of a disk having slightly rising/constant rotation curve (offsets $\leq 30-35\arcmin$). The brightness contours then exhibit decreasing velocities at absolute offsets $>35\arcmin$. This apparent decrease occurs in the same region as the one where the twist of the iso-brightness contours are observed in the channel maps, as caused by the M33 warp. Along the direction PA=165°, the high brightness emission follows again a slowly rising/flat rotational pattern.

The most important result here is the detection of a low brightness Hi component that has an extremely anomalous velocity behaviour 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 towards a line-of-sight velocity that deviates from those of contours of gas with the highest density.

First, the low brightness emission extends towards 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\arcmin$ shows a high-velocity bump. Third, the low brightness emission leaks in the forbidden velocity zones, which makes it being in apparent counter-rotation with respect to the high-density gas. In our diagrams, the forbidden zones are the bottom-left and top-right quadrants. A Hi component in the approaching disk half is in the forbidden velocity zone when its radial velocity is larger than the systemic velocity (“receding component on the approaching side”). Conversely, a Hi component in the receding disk half is in the forbidden velocity zone when its radial velocity is smaller 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 Fig. 5 the gas component with forbidden velocity by extracting Hi profiles at several positions along the PA=165° and 202° directions. Only channels with $\rm T_{B}>5\sigma$ are shown for clarity. The equatorial coordinates of these spectra are given in Tab. 3. The Hi gas in the forbidden zone is the asymmetric tail of the profiles, with velocities differing by up to 80-100km s*-1* from the main peak. Hi gas in the PA=202° forbidden zone is also observed as distinct Hi peaks with velocities differing by 100 km s*-1* from the main peak (offsets$=-30\arcmin,+10\arcmin$).

The first anomalous component is reminiscent of the slow rotation extraplanar Hi layer seen in other galaxies (e.g. NGC2403, NGC891, NGC253, Fraternali et al., 2002; Oosterloo et al., 2007; Lucero et al., 2015). The second component reminds 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) already evidenced high-velocity gas around M33, but at larger distances than the new one 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 (NGC2403, Fraternali et al., 2001, 2002). Interestingly, these authors reported that in NGC2403 the low angular momentum and the forbidden velocity components seem to form a coherent structure. M33 is thus similar to NGC2403 in that 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 (Chemin et al., 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, maybe like in Fraternali et al. (2002) with NGC2403, who fitted and subtracted to the datacube a Gaussian profile centred on the highest Hi peak, or in Chemin et al. (2012), who fitted multiple Gaussians peaks to the current M33 dataset.

4 Hi distribution and kinematics

4.1 Moment maps

The task moments in gipsy has been used to compute the moment maps using the data cube smoothed at $120\arcsec\times 120\arcsec$. Figures 6 and 7 show the maps of the $0^{\rm th}$ moment (integrated emission), $1^{\rm st}$ (velocity field) and $2^{\rm nd}$ (velocity dispersion).

The Hi emission map shows that the high surface density gas is contained within the stellar disk ($R\lesssim 30-35\arcmin$, or $\lesssim 8$ kpc). It displays a multiple spiral arms pattern that coincides well with the spiral structure evidenced in the H$\alpha$ 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*-2* in the outermost regions. Gaseous tails and arcs are observed to the North-West and South-East of the disk as part of the Hi 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 $\sim 2$° (30 kpc in projection) to the North-West towards M31, nearly three times farther out than the size of the M33 stellar disk, and thus farther than the Hi gas in the warp.

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

The velocity field is quite regular in the inner regions ($R\lesssim 8$ kpc). Beyond that 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 evidenced in the channel maps (Sect. 3). Note also an outer South-Western cloud or extension, illustrated by yellow and orange colours in the velocity map, that rotates slower than gas inside the disk at similar radius. This extension has already 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 Hi peaks or wider single Hi peaks (see also Chemin et al., 2012). These features are likely due to the crowding of gas orbits because of the disk warping (see §4.4 and Fig. 11).

4.2 Derivation of the kinematical parameters

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

[TABLE]

to the Hi velocity field, where $\theta$ 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 semi-major axis. The angle $\theta$ is measured relatively to the semi-major axis. The South-Western Cloud has been masked for the derivation of the parameters and rotation curve.

In a first run, the position of the rotation centre and the systemic velocity $\rm V_{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 $\rm V_{sys}\sim-183$ km s*-1*, and a centre at $(\alpha,\delta)_{\rm J2000}=(01{\rm h}33{\rm m}50.9{\rm s},+30\arcdeg 40^{\prime}20^{\prime\prime})$. Since the positional difference with respect to the photometric centre is much less than the beam size, we decided to fix the kinematical centre at the position of the photometric centre. Moreover, because the inferred systemic velocity is in agreement with the mean value determined from the integrated profile (§4.1), we adopted $\rm V_{sys}=-180$ km 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 fixed centre and systemic velocity given by the adopted values. The results of the tilted ring model are shown in Fig. 9 and the adopted inclination and position angle are listed in Tab. 4.

4.3 The HI rotation curve and velocity dispersion profile

The Hi rotation curve (Fig. 9) is then derived with fixed values for the kinematical centre 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 Hi rotation curves for the approaching ($\rm V_{a}$) and receding sides ($\rm V_{r}$) fitted separately with similar kinematical parameters. Differences between $\rm V_{a}$ and $\rm V_{r}$ are mostly smaller than 10 km s*-1*. The largest differences, up to 60 km s*-1*, occur at $R>80\arcmin$ in the perturbed NW and SE structures from the warp.

The Hi rotation curve is reported in Tab. 4. The total velocity uncertainty $\rm\Delta V_{rot}$ is defined by $\rm\Delta V^{2}_{rot}=\epsilon^{2}+|(V_{a}-V_{r})/2|^{2}$ with $\epsilon$ being the formal RMS error for the model with both sides fitted simultaneously (solid line in the top panel of Fig. 9). The dominant error of the total uncertainty is the velocity difference between the two disk sides. Because of that definition, the errors appear larger than in many other studies. However, we prefer being conservative because we model the rotation curve with axisymmetric components (Sect. 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-sigth 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 $\sim 9$ km s*-1* in the central regions ($R<30\arcmin$). The signature of the ring of higher dispersion is observed as a rise of 3 km s*-1* ($R=35-40\arcmin$). The dispersion then continuously decreases at larger radius. Cuts made along the semi-major axis of the approaching and receding sides in the dispersion map show significant differences inside $R=30\arcmin$ and beyond $R=50\arcmin$. Moreover, the mean dispersion is larger than values along the semi-major 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 §4.4 and 4.5).

4.4 The HI warp of M33

The warping of M33 is mainly observed as a significant kinematical twist of the major axis position angle, starting by a sharp decrease from $\sim 202\arcdeg$ to $\sim 170\arcdeg$ within $7<R<11$ kpc, and then a smooth variation down to $165\arcdeg$ at large radius. A kinematical tilt is also evidenced, though less spectacular, as a small rise of inclination ($\sim 5\arcdeg$ from the centre to the disk outskirts).

The geometry of the Hi orbits implied by the axisymmetric tilted-ring model is shown in Fig. 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 Hi peaks at e.g. locations C and D, or more generally by the incomplete ring-like structure of larger velocity dispersion (Fig 6). The change of inclination is so insignificant here that a configuration where the line-of-sight crosses the disk more than once is hardly possible. Moreover, we verified that the wide and double Hi 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 Hi 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).

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 $\sim 30-35\arcmin$ in the PA=175° PV diagram of Fig. 4, which offsets correspond to $(\alpha,\delta)_{\rm J2000}=(01{\rm h}33{\rm m}57.6{\rm s},+30\arcdeg 05^{\prime}02^{\prime\prime})$ and $(01{\rm h}33{\rm m}20.2{\rm s},+31\arcdeg 12^{\prime}30^{\prime\prime})$. Here, these wider profiles coincide with holes in the gas distribution and trace high- and low-velocity faint gas (§3.2).

Elsewhere than in the transition region, the inferred model geometry is regular and wide profiles are not 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 outermost ellipses have their major axis 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 Hi disk in the direction of M31, and that perturbation is acting down to a radius of $R\sim 7$ 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 that radius.

The model velocity map based on the Hi RC allows us to make a residual velocity map, defined here as the observation minus the model velocity field (Fig. 12). No systematic velocity asymmetry is found, implying the goodness of the adopted mass centre 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 Hi disk, at the level of a few km 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 $\sim 30$ km s*-1* (absolute value). The ring of higher velocity dispersion (Fig. 6) thus coincides with large residuals. Note also the significant residual for the SW clump ($>30$ km s*-1*, absolute value), confirming that it does not rotate similarly as the disk. One can notice finally that observing larger residuals in the NW and SE Hi extensions confirms the existence of asymmetric gas orbits at the disk periphery.

4.5 Asymmetric and non-circular motions in M33

We can assess further the perturbations in the Hi velocity field of M33 by expanding the standard model of Eq. 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 fitted to the Hi velocity field:

[TABLE]

where $c_{k}$ and $s_{k}$ are the velocity coefficients of harmonic order $k$ ($k$ is an integer). The $c_{0}$ coefficient is the systemic velocity $\rm V_{sys}$ of Eq. 1, the $c_{1}$ coefficient is the rotation curve, the $s_{1}$ term is the noncircular radial velocity, while higher order terms constrain deviations from axisymmetry of $c_{1}$ and $s_{1}$.

We restrict the model to $k=4$, implying that a perturbation of the potential of order of up to $m=3$, maybe $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 Sect. 4.2. The angular sampling is $\sim$490 pc (120″). We derived the amplitudes of the noncircular and asymmetric motions, given respectively by ${\rm v_{1}}=|s_{1}|$ and ${\rm v_{k}}=(c_{k}^{2}+s_{k}^{2})^{1/2}$ for $k>1$, and ${\rm v_{0}}=\sqrt{(c_{0}-{\rm V_{sys}})^{2}}$. These amplitudes are shown in Fig. 13. The uncertainties are the $1\sigma$ formal errors from the fit.

One first observes increasing amplitudes for every $k$ terms at $R\gtrsim 15$ kpc. These may be other signatures of the interaction with M31 on the outer M33 velocity field. Interestingly, a similar trend was observed for v2, v3 and v4 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 v0 and v2 are similar. Larger amplitudes are observed in the inner 5 kpc, as well as at $R\sim 8$ and 12 kpc. This is evidence for a kinematical lopsidedness, i.e. a $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<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 dark matter may be lopsided as well 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 v1 and v3 are similar. The $k=1$ and $k=3$ terms are small inside $R=7$ kpc and larger at $R\sim 11$ kpc, within the zone of strong major axis PA variation. Smaller values in the inner disk reflects perfectly the observation that no bisymmetric structure dominates in the inner density map (Fig. 6). On another hand, larger values at large radius reflects the dynamical impact of the M33 warp.

Thirdly, a bump is evidenced for the $k=4$ term at $R\sim$ 8 kpc, coincinding with those seen for the $k=0$ and $k=2$ terms. Whether this feature is caused by the same $m=1$ perturbation is something that 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{\rm v}_{1}\rangle=3.8\pm 0.6$ km s*-1*. The average amplitude of asymmetries are $\langle{\rm v}_{0}\rangle=3.3\pm 0.3$ km s*-1*, $\langle{\rm v}_{2}\rangle=4.9\pm 0.6$ km s*-1*, $\langle{\rm v}_{3}\rangle=4.9\pm 0.7$ km s*-1*, and $\langle{\rm v}_{4}\rangle=3.7\pm 0.4$ km s*-1*, again showing the more important impact of $m=1$ and $m=2$ perturbations on the gravitational potential of M33. The level of such asymmetries is consistent with the asymmetry measured between the rotation curves of the approaching and receding disk sides.

4.6 Comparisons with previous works

4.6.1 DRAO-Arecibo versus VLA-GBT

The significant twist of the major axis of the Hi velocity field we measure at large radius from our DRAO-Arecibo dataset is very consistent with results obtained from VLA-GBT data by Corbelli et al. (2014), or from older 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 Hi 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 one from Corbelli et al. (2014) is not straightforward as well because their curve was derived from these two different models of tilted disk. The DRAO and VLA rotation curves are shown in Fig. 14 and agree well within the uncertainties for $R\leq 20$ kpc. Then the shapes unsurprinsingly differ beyond $R=20$ kpc. A more detailed investigation shows that when we compare with the curve they obtained from a rotcur-like model similar to ours, the shapes become comparable at these radii, showing a drop out to $R\sim$22 kpc, followed by a rise. Therefore it is the action of the 3D modeling they made of the datacube that causes the shape difference at the largest radii.

What is clear, however, is that both Hi rotation curves are perturbed at these radii, harbouring larger scatter and uncertainties. We conclude that the Hi rotation curve of M33 is only reliable out to $\sim$20 kpc, irrespective of the method used to derive it. Fortunately, this outermost region has negligible impact on the mass distribution models described in Sect.5 because we used as weighting function of the rotation velocities the inverse of the squared errors.

4.6.2 The hybrid H$\alpha$ - HI rotation curve of M33

Figure 14 also shows the very good agreement between the DRAO Hi rotation curve with that of the ionized gas disk from Kam15. The H$\alpha$ rotation curve naturally harbours more wiggles than the 21-cm data because of its higher resolution.

The analysis of the mass distribution (Sect. 5) is based on a hybrid rotation curve which is obtained by combining the H$\alpha$ velocities from Kam15 for $R\leq 6.5$ kpc with the Hi velocities for $R>6.5$ kpc. This radius was chosen sligthly before the location from where the Hi warping starts. The advantages of using a hybrid curve for mass models are to benefit from the high sampling of the H$\alpha$ data (20 pc, 5″) to constrain the central velocity gradient more accurately than with our radio interferometry only (490 pc, 120″ sampling), and from the large extent of the neutral gas disk. The H$\alpha$ rotation curve within $R=6.5$ kpc has $\sim 10$ times more velocity points than the Hi curve for $R>6.5$ kpc, yielding a final rotation curve with 353 data points at a hybrid resolution of 20 pc-490 pc. The hybrid curve is listed in Tab. 9 of Appendix B.

5 Mass distribution models of Messier 33

The modeling of the mass distribution fits a model velocity profile to the hybrid H$\alpha$-Hi rotation curve of Messier 33. The contributions from the luminous matter are those of the gaseous and stellar disks, inferred from mass surface density profiles (§5.1 and §5.2). We have considered models with dark matter (DM) using two different forms for the DM halo: the pseudo-isothermal sphere (ISO) and the Navarro-Frenk-White (NFW) models (§5.3). We also did the modeling within the framework of Modified Newtonian Dynamics (§5.4).

5.1 Stellar component

In Kam15, the M33 stellar bulge was shown to have a negligible impact in the mass distribution (see also Corbelli & Walterbos, 2007; Corbelli et al., 2014). Therefore only the velocity contribution from the stellar disk, $\rm 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 the most to the stellar mass. Kam15 derived the surface brightness profile of M33 at 3.6 $\mu$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 $\mu$m, $\Upsilon$, has been estimated from infrared colours and stellar population synthesis (SPS) models following the prescriptions given by Oh et al. (2008) and de Blok et al. (2008):

[TABLE]

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$ colour index from Jarrett et al. (2003), Kam15 derived $\Upsilon=0.72\pm 0.1$. This corresponds to a stellar mass of $(7.6\pm 1.1)\ 10^{9}$ M⊙, which is 38% (58%, respectively) higher than the fixed (best-fit) stellar mass given in Corbelli et al. (2014). Note that $\Upsilon\sim 0.7$ is the upper limit expected from mass distribution modeling of the large galaxy sample of Lelli, McGaugh & Schombert (2016).

This discrepancy is the reason why we have also performed models with $\Upsilon=0.52$ to facilitate comparisons with the results of Corbelli et al. (2014). This value of $\Upsilon$ was obtained by scaling the maximum of the stellar velocity curve to $\sim 70$ km s*-1*, which corresponds to the highest contribution from stars in the fixed stellar mass fits of Corbelli et al. (2014). $\Upsilon=0.52$ is valid here in the Spitzer/IRAC 3.6 $\mu$m band only, and corresponds to a stellar mass of $5.5\ 10^{9}$ M⊙, which thus agrees perfectly with the value from the SPS model of Corbelli et al. (2014). Assuming $\Upsilon=0.52$ that is constant with radius in our models is a good way to account for a stellar velocity contribution 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, McGaugh & Schombert (2016), and appears more in agreement with most of galaxies of their sample as bright as M33 than $\Upsilon=0.72$.

5.2 Atomic and molecular gas components

The velocity contribution from the disk of neutral gas, $\rm V_{atom}$, is derived from our deep DRAO+Arecibo dataset by integrating the total atomic gas mass surface densities (Fig. 8) multiplied by a factor of $\sim$1.3 to take into account the contribution from Helium.

The molecular gas, 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\arcsec-15\arcsec$. Tosaki et al. (2011) report that the Hi 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 Hi density is lower than that of the molecular gas. Using a conversion factor of $\rm N(H_{2})/I_{CO(1\rightarrow 0)}=4\times 10^{20}cm^{-2}/$(K km s*-1*), twice the value found for the Milky Way (M33 is half the solar metallicity), they measured an average density of $\Sigma_{H_{2}}=8.5\pm 0.2$ M⊙ $\rm pc^{-2}$ for the central kpc, and a total molecular gas mass of $\sim 3.3\times 10^{8}$ M⊙ for the entire M33 disk. We used the H2 density profile derived from Druard et al. (2014), scaled by a factor of $\sim 1.3$ to infer the velocity contribution of the total molecular gas component, $\rm V_{mol}$.

5.3 Dark Matter halo component

The total rotation velocity $\rm V_{rot}$ in the models with a DM component is defined by:

[TABLE]

where $\rm V_{\star}$ is the contribution from the stellar disk and $\rm V_{\rm gas}=(V_{atom}^{2}+V_{mol}^{2})^{1/2}$ the total contribution from the gaseous disk, as deduced in §5.2, and $\rm V_{DM}$ from a DM halo assumed spherical.

5.3.1 The pseudo-isothermal sphere model

Here, the density profile of DM is given by:

[TABLE]

The corresponding circular velocities are:

[TABLE]

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\rightarrow 0$ (the halo core has a density 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:

[TABLE]

where $\rho_{i}\thickapprox 3{\rm 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:

[TABLE]

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

5.4 Modified Newtonian Dynamics mass models

An alternative to dark matter to explain the missing mass problem is MOdified Newtonian Dynamics (Milgrom, 1983a, b). MOND has been successful to reproduce correctly 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_{0}$, the classical Newtonian dynamics is no more valid and the law of gravity is modified.

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

[TABLE]

where $g$ is the acceleration, $g_{N}$ the Newtonian acceleration, and $\mu(x)$ is an interpolating function that must satisfy: $\mu(x)=x$ for $x<<1$ and $\mu(x)=1$ for $x>>1$. The MOND velocity profile thus depends on $\mu(x)$. For the “standard” $\mu$-function proposed by Milgrom (1983a)

[TABLE]

the MOND rotation velocity is:

[TABLE]

with $\rm V_{lum}^{2}=V_{\star}^{2}+V_{gas}^{2}$, $\rm V_{\star}$ and $\rm V_{gas}$ being as in Eq. 4. The standard scale acceleration is $a_{0}=\rm(H_{0}/75)^{2}\times 1.2\ 10^{-8}=0.99$ cm s*-2* ($\rm H_{0}=68$ km s*-1* Mpc*-1*).

For the simple function

[TABLE]

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

[TABLE]

Both $\rm V_{STD}$ and $\rm V_{SIM}$ models were fitted to the rotation of M33, with free $a_{0}$ and fixed or free $\Upsilon$.

5.5 Results and analysis

We performed non-linear Levenberg-Marquardt least-square fits of the ISO, NFW and MOND models. The ISO and NFW fits have two free parameters at fixed $\Upsilon$, three when $\Upsilon$ is left free. The MOND fits have one free parameter at fixed $\Upsilon$, and two at free mass-to-light ratio. A normal weighting function set to $\rm\Delta V_{rot}^{-2}$ is used. Figures 15 and 16 show the mass models and Tables 5 list the fitted parameters. A drawback of using different datasets to make an hybrid rotation curve is that the resulting distribution of velocity uncertainties is rarely Gaussian. We measure for instance a median H$\alpha$ velocity uncertainty that is $\sim 60\%$ that of the median Hi velocity uncertainty. A consequence of the error non-homogeneity, combined with the non-uniform sampling of the rotation curve, is to yield large $\chi^{2}$ (reduced values $\gtrsim 6$), which, taken individually, means that a fit has no statistical significance despite of the $\sim 350$ degrees-of-freedom. We thus estimated the differences of $\chi^{2}$, relatively to a given model, that is, $\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 $\Delta\chi^{2}$ values are reported in Tabs 6.

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

The analysis of the residual rotation velocities (bottom insert in Fig. 15) shows that NFW halo reproduces more correctly the inner $R=2$ kpc of the rotation curve than the core-dominated model, which later implies a total velocity model always smaller than the observation, irrespective of the value of $\Upsilon$. An opposite trend is observed within $R=9-16$ kpc for $\Upsilon=0.52$ and $R=9-18$ kpc for $\Upsilon=0.72$, where the ISO halo is more appropriate than the cusp, which later implies a total model always smaller than the rotation curve. At intermediate radii ($2<R<6.5$ kpc) and in the outer disk ($R\gtrsim 19$ kpc), the rotation curve is reproduced equivalently by the cusp and core models. None of the models can reproduce the rotation curve for $R=6.5-8$ kpc, where the warp starts.

In a more global framework, the analysis of $\Delta\chi^{2}$ first shows that the density shape which is more approriate at $\Upsilon=0.52$ is that of the NFW halo, while at $\Upsilon=0.72$ it is that of the core-dominated halo. Then, going from $\Upsilon=0.52$ to $\Upsilon=0.72$ improves the modeling for the ISO halo, while it is the opposite way from $\Upsilon=0.72$ to $\Upsilon=0.52$ that the modeling has been improved for the NFW halo. Therefore, the NFW halo can only coexist with $\Upsilon=0.52$ (lower stellar mass) whereas the ISO model can only coexist with $\Upsilon=0.72$ (higher stellar mass). The models compensate by a cuspier density from 0.72 to 0.52 to get a good fit in the central regions. This trend is confirmed by the free $\Upsilon$ (best-fit) models, as $\Upsilon=0.56$ is found for NFW (stellar mass of $5.9\ 10^{9}$ M⊙) and $\Upsilon=0.8$ is found for ISO (stellar mass of $8.5\ 10^{9}$ M⊙). Note also the anti-correlation of the cusp concentration with the stellar mass (the larger the mass, the smaller the concentration).

We also deduce from $\Delta\chi^{2}$ that the combination $\Upsilon=0.52$/NFW is more likely than the combination $\Upsilon=0.72$/ISO. The NFW halo as more appropriate model is also apparent at free mass-to-light ratio. Consequently, our mass models are more consistent with the SPS modeling of Corbelli et al. (2014) than with the one used by Kam15. A stellar mass-to-light ratio as large as $\Upsilon=0.72$ is therefore excluded for M33.

We therefore adopt the cusp obtained with $\Upsilon=0.52$ as the most likely dark matter halo for M33, since this mass-to-light ratio has a physical meaning, as based on SPS models described in Corbelli et al. (2014). The results obtained at free $\Upsilon$ can then be seen as a way to confirm the $\Upsilon=0.52$ results. With $\Upsilon=0.52$, $\rm V_{200}=114$ km s*-1*, $c=6.1$, the inferred mass of Messier 33 is $\log(M_{R\leq 168}/\rm M_{\odot})=11.72$ within the virial radius of the cusp, $R=R_{200}=168$ kpc. This total mass estimate is roughly half those of the Milky Way or the Andromeda galaxy (Chemin et al., 2009; Bland-Hawthorn & Gerhard, 2016, and references therein). The implied baryonic fraction is 2% at that radius, which strongly differs from the cosmic value, $\Omega_{b}/\Omega_{m}=15.7\%$ (Planck Collaboration, 2016). This estimate does not include the unknown mass of warm and hot gas. It thus points out that the M33 virial radius cannot be as large as $R_{200}$, unless M33 violates the cosmic value.

Within a radius of $R=23$ kpc, at the last point of the rotation curve, the total mass is $\log\left(M_{R\leq 23}/\rm M_{\odot}\right)=10.90$ (or $\sim 13$ times less massive than the Galaxy or M31). The mass fraction of baryons is 11%, which is more consistent with the cosmic value. This implies that if one assumes that the baryonic fraction of M33 and the cosmic value must be similar, then $R\sim 23$ kpc is close to the real location enclosing the real total mass of M33. By integrating the density of the adopted NFW halo, we estimate that $R=17-18$ kpc is the radius where the M33 baryonic mass fraction equals the cosmic value. Interestingly, this location is very close to the radius where the rotation curve starts to drop ($R\sim 19$ kpc), whose characteristic is expected to occur beyond the radius that encompasses the total galaxy mass.

5.5.2 Results for MOND

All MOND models strongly favour the original, standard interpolation function over the more simple one (Tab. 6). Figure 16 thus only presents results obtained with the standard interpolation function. Moreover, configurations with $\Upsilon=0.52$ are preferred over $\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 $\mu$-function, which corresponds to a stellar mass of $5.3\ 10^{9}$ M⊙. The scale acceleration $a_{0}=1.76\ 10^{-8}$ cm s*-2* is by 78% larger than the standard value of $a_{0}=1.2\ 10^{-8}\times(\rm H_{0}/75)^{2}$ cm s*-2*, for $\rm H_{0}=68$ km s*-1* Mpc*-1*.

The comparison between MOND and the ISO or NFW models can not be done 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. ${\rm AIC}=2N+\chi^{2}$. Consequently, comparing AIC residuals $\rm\Delta AIC=AIC_{MOND}-AIC_{DM}$ is well appropriate to find which of MOND and ISO or NFW is more likely than the other (see Chemin et al., 2011; Chemin et al., 2016). The AIC residuals are reported in Tab. 7 for the standard $\mu$-function. The meaning of $\rm\Delta AIC$ is similar to that of $\Delta\chi^{2}$. Negative residuals mean MOND more likely than dark matter models, while positive residuals imply MOND less likely than dark matter models. It is found that at fixed mass-to-light ratio $\Upsilon=0.52$, the NFW model is more appropriate than MOND, which is more likely than the ISO model. At $\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 $\Upsilon$, both DM haloes are more appropriate than MOND. In other words, a more elaborated form than MOND, with a density law of hidden mass and an additional free parameter have 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. This simply illustrates the preference for DM-based models to explain the mass distribution underlying the rotation curve of M33.

5.6 Comparison with previous works

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 $\Upsilon<2$. Their modeling however differs from ours since they used an additional low mass bulge component, and found more likely models having density slopes steeper than the NFW cusp, and with $\Upsilon\sim 1.5$. Such large $\Upsilon$ 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, McGaugh & Schombert, 2016).

At fixed stellar mass, finding a combination $\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_{\rm 200}$.

At free stellar mass, the inferred stellar mass is $5.9\ 10^{9}$ M⊙ (corresponding to a maximum velocity of 73 km s*-1*), which is 23% larger than the most likely mass of Corbelli et al. (2014) (corresponding to a maximum velocity of 60 km s*-1*). Such a 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 of best-fitting of 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>6.5$ kpc) DRAO velocities with the inner ($R<6.5$ kpc) curve from Corbelli et al. (2014), or our inner H$\alpha$ 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 $\rm V_{200}=99$ km s*-1* and $\Upsilon=0.32$. This corresponds to a maximum velocity of 55 km s*-1* for the stellar contribution. It is another way to illustrate the halo concentration-stellar mass degeneracy shown in Sect. 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 dataset and adopted cusp model not tied to the halo mass-concentration relation strongly favours a concentration in disagreement with $\Lambda$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, it is only once it is known with more accuracy that the cusp concentration conflict will be alleviated in M33.

6 Conclusion

New high sensitivity Hi observations of M33 obtained with the DRAO interferometer have been presented. Combined with the single dish Arecibo data from Putman et al. (2009), the dataset reach column densities as low as $\sim 5~{}\times~{}10^{18}$ cm*-2* in the outer Hi disk of M33.

The main results on the Hi distribution and kinematics of M33 are:

• •

While the bulk of the Hi gas is found within the stellar disk ($\leq R_{25}\sim 8$ kpc), the Hi distribution is traced out to $\geq 2.7R_{25}$. It is irregular in the outer region, in the form of tails or arc-like features or isolated clumps. More gas is detected in the Southern receding disk side, and the Hi emission is more extended to the North-West than to the South-East of the disk, implying a lopsided Hi distribution. The surface density is nearly constant out to the edge of the stellar disk, and then drops abruptly. At the adopted distance of 840 kpc, the Hi mass is $\sim 2\times 10^{9}$ M⊙  for a $\rm M_{HI}/L_{V}\simeq 0.2$.

• •

Position-velocity diagrams make it possible to evidence “beard-like” contours from a low brightness component with an important velocity scatter, as it is observed lagging and exceeding the disk rotation, as well as leaking in the forbidden velocity zone of apparent counter-rotation.

• •

The rotation curve is in very good agreement with the H$\alpha$ curve of Kam et al. (2015) in the inner disk, and consistent with the Hi rotation curve of Corbelli et al. (2014) inside $R=20$ kpc. Beyond that radius, the Hi rotation curves of the approaching and receding disk sides differ by up to 60 km s*-1*.

• •

The warp of M33 consists mainly in a strong twist of the position angle of the kinematical major axis ($\sim 40\arcdeg$) beyond $R=7$ kpc. This result is in perfect consistency with previous Hi studies. Only a minor increase of inclination is detected throughout the entire disk ($\sim 5\arcdeg$).

• •

Wider and double-peaked Hi profiles are evidenced in a large-scale, incomplete ring-like structure of larger dispersion. They coincide with the transition zone of twist of the major axis position angle between the inner and outer regions. The crowding of rings inferred by our warp model naturally explains part of the larger dispersion and multiple peaks. Collisions of gas clouds are expected in this region, and the gas orbits are likely elongated in the direction to the companion Messier 31. Other wider Hi profiles that are not in the crowded rings zone are associated to holes in the Hi distribution.

• •

A Fourier series analysis of the velocity field reveal non-circular and asymmetric motions, suggesting perturbations of the first and second order of the gravitational potential of M33. The asymmetric motions are all observed to increase in the disk outskirts.

• •

The past tidal interaction with Messier 31 already evidenced by large-scale gas and stellar surveys of the M33 enviromnent (e.g. Braun & Thilker, 2004; Putman et al., 2009; Ibata et al., 2014) is likely at the origin of most of the perturbed Hi kinematics and morphology presented here.

The main results of the M33 mass distribution modeling from the hybrid H$\alpha$-Hi rotation curve are:

• •

The most likely density shape of dark matter is cuspy, to the detriment of a pseudo-isothermal sphere. The concentration of the most likely NFW halo disagrees with that expected by the halo mass-concentration from CDM numerical models, or from previous Hi studies. Modified Newtonian Dynamics is less likely than models with a dark matter halo.

• •

The mass enclosed within the virial radius of the best-fit NFW halo (168 kpc) is $\sim 5.2\,10^{11}$ M⊙, implying a very low baryonic mass fraction (2%), in conflict with the cosmic value of 15.7%. That result suggests a more plausible M33 virial radius well smaller than that of the adopted cusp.

• •

The most likely mass of the stellar disk is $5.5\,10^{9}$ M⊙, only about 3 times larger than the mass of disk of neutral Hydrogen. The enclosed mass within $R=23$ kpc at the last point of the rotation curve is $\sim 7.9\,10^{10}$ M⊙. Luminous matter represents about 11% of that mass, in better agreement with the cosmic value. A radius as low as the radius of the Hi disk could thus be very nearby the true location encompassing the total mass of M33.

Acknowledgements

The Dominion Radio Astrophysical Observatory is operated as a national facility by the National Research Council of Canada. ZKS’s work was supported by CC’s Discovery grant of the Natural Sciences and Engineering Research Council of Canada. The work of CC and THJ is based upon research supported by the South African Research Chairs Initiative (SARChI) of the Department of Science and Technology (DST), the Square Kilometer Array South Africa (SKA SA) and the National Research Foundation (NRF). The research of EE has been supported by an SKA SA fellowship. LC acknowledges financial supports from CNPq (PCI/INPE) project 301176/2016-7, Comité Mixto ESO-Chile, and DGI University of Antofagasta. TF’s work has partially been supported by a grant from the Brandon University Research Committee (BURC). We are very grateful to Kevin Douglas and Mary Putman for providing us the Arecibo TOGS data, to Clément Druard for having provided the emission map of the molecular gas, to Pierre Gratier for sharing the high-resolution VLA datacube, and to Jonathan Braine for his comments that greatly helped improved the analysis and the manuscript.

Appendix A Hi dispersion and mass surface density of M33

Appendix B Hybrid H$\alpha$-Hi rotation curve of M33