© ESO, 2014

1. Introduction

Rotation curves of spiral galaxies are fundamental tools to study the visible mass distributions in galaxies and to infer the properties of any associated dark matter halo. Knowledge of the radial halo density profile from the center to the outskirts of galaxies is essential for solving crucial issues at the heart of galaxy formation theories, including the nature of the dark matter itself. Numerical simulations of structure formation in the flat cold dark matter cosmological scenario (hereafter ΛCDM) predict a well-defined radial density profile with a central cusp for the collisionless particles in virialized structures: the Navarro-Frenk-White (hereafter NFW) profile, originally derived by Navarro et al. (1996, 1997). The two parameters of the “universal” NFW density profile are the halo overdensity and the scale radius, or (in a more useful parameterization) the halo concentration and its virial mass. For hierarchical structure formation, small galaxies should show the highest halo concentrations and massive galaxies the lowest concentrations (Macciò et al. 2007). The relation between these parameters is still under investigation since numerical simulations improve steadily by considering more stringent constraints on the cosmological parameters, by using a higher numerical resolution, and by considering phenomena that can affect the growth or time evolution of dark matter halos. The ΛCDM dark halo models have been shown to be appropriate for explaning the rotation curve of some spiral galaxies (e.g., Martinsson et al. 2013), including that of the nearest spiral, M 31, which has been traced out to large galactocentric distances by the motion of its satellites (Corbelli et al. 2010). Dwarf galaxies with extended rotation curves have instead contradicted cold dark matter scenarios since central regions show constant density cores (Gentile et al. 2005, 2007). Shallow central density distributions and very low dark matter concentrations have been supported also by high resolution analysis of rotation curves of low surface brightness galaxies (de Blok et al. 2001; Gentile et al. 2004). This has cast doubt on the nature of dark matter on one hand, but has also given new insights into the possibility that the halo structure might have been altered by galaxy evolution (e.g., Gnedin & Zhao 2002). Recent hydrodynamical simulations in ΛCDM framework have also shown that strong outflows in dwarf galaxies inhibit the formation of bulges and decrease the dark matter density, thus reconciling dwarf galaxies with ΛCDM theoretical predictions (Governato et al. 2009, 2010). Also, shallower density profiles than true ones might be recovered in the presence of small bulges or bars (Rhee et al. 2004). M 33 provides useful tests for the ΛCDM cosmological scenario since it is higher in mass than a dwarf, but it hosts no bulge nor prominent bars (Corbelli & Walterbos 2007). The absence of prominent baryonic mass concentrations alleviates the uncertainties in the inner circular velocities and makes it unlikely that the baryonic collapse and disk assembly alter the distribution of dark matter in the halo.

M 33 is a low-luminosity spiral galaxy, the most common type of spiral in the Universe (Marinoni et al. 1999) and the third most luminous member of the Local Group. It yields a unique opportunity to study the distribution of dark matter because it is richer in gas, and more dark matter dominated than our brighter neighborhood M 31. Because of its proximity, the rotation curve can be traced with unprecedented spatial resolution. M 33 is in fact one of the very few objects for which it has been possible to combine a high quality and high resolution extended rotation curve (Corbelli & Salucci 2000) with a well determined distance (assumed in this paper to be D = 840 kpc, Freedman et al. 1991; Gieren et al. 2013), and with an overwhelming quantity of data across the electromagnetic spectrum. Previous analyses have shown that the HI and CO velocity fields are very regular and cannot be explained by Newtonian dynamic without including a massive dark matter halo (Corbelli 2003). On the other hand, the dark matter density distribution inside the halo is still debated because of the unconstrained value of the mass-to-light ratio of the stellar disk. These uncertainties leave open the possibility to fit the rotation curve using non-Newtonian dynamic in the absence of dark matter (Corbelli & Salucci 2007).

The extent of the HI disk of M 33 has been determined through several deep 21 cm observations in the past (e.g., Huchtmeier 1978; Corbelli et al. 1989; Corbelli & Schneider 1997; Putman et al. 2009). Beyond the HI disk, Grossi et al. (2008) found some discrete HI clouds but there is no evidence of a ubiquitous faint HI distribution. This is likely due to a sharp HI-HII transition as the HI column density approaches N_HI ~ 10¹⁹ cm^-2: hydrogen becomes mostly ionized by the UV extragalactic background radiation (Corbelli & Salpeter 1993) and the HI column density often drops below the detection limits. Previous observations of the outer disk have been carried out at low spatial resolution, mostly with single dish radio telescopes (FWHM> 3 arcmin) that needed sidelobe corrections. The new combined Very Large Array (VLA) and Green Bank Telescope (GBT) HI survey of M 33, presented in this paper, has an unprecedented sensitivity and spatial resolution. It is one of the most detailed databases ever made of a full galaxy disk, including our own. Our first aim is to use the combined VLA+GBT HI survey to better constrain the amplitude and orientation of the warp, and the rotation curve. The final goal will be to constrain the baryonic content of M 33 disk and the distribution of dark matter in its halo through the dynamical analysis of the rotation curve.

There are also crucial, but now questionable, assumptions that were made in the prior analyses that leads us to further investigate the baryonic and dark matter content of M 33 once again. Firstly, a constant mass-to-light ratio has been used throughout the disk and determined via a dynamical analysis of the rotation curve. Futhermore, the light distribution has also been fitted using one value of the exponential scalelength in the K-band. This might not be consistent with the negative radial metallicity gradient that supports an inside-out formation scenario (Magrini et al. 2007, 2010; Bresolin 2011). Portinari & Salucci (2010) investigated the radial variations of the stellar mass-to-light ratio using chemo-photometric models and galaxy colors. They indeed found a radially decreasing value that affects mass modeling. Evidence of a radially varying mass-to-light ratio has also been found through spectral synthesis techniques in large samples of disk dominated galaxies (Zibetti et al. 2009; González Delgado et al. 2014). In order to overcome this simplified constant M/L ratio model, and given the proximity of M 33, in this paper we shall build a surface density map for the bright star forming disk of M 33. Deep optical surveys of the outskirts of M 33 have shown that the outer disk is not devoid of stars and that several episodes of star formation have occurred (e.g., Davidge et al. 2011; Grossi et al. 2011; Sharma et al. 2011, and references therein). It is therefore important to account for the stellar density and radial scalelength in the outer disk of M 33 before investigating the dark matter contribution to the outer rotation curve.

In Sect. 2, we describe the data used to retrieve the velocities and gas surface densities along the line of sight. The details of the reduction and analysis of the optical images used to build the stellar surface density map are given in Sect. 3. In Sect. 4, we infer the disk orientation and the rotation curve. Some details of the warp and disk thickness modeling procedures are given in the Appendix. In Sect. 5, we derive the surface density of the baryons. The rotation curve dynamical mass models and the dark halo parameters are discussed in Sect. 6, also in the framework of the ΛCDM scenario. Section 7 summarizes the main results on the baryonic and dark matter distribution in this nearby spiral.

2. The gaseous disk

The HI data products utilized for our study originate from a coordinated observing campaign using the VLA and the GBT. Multi-configuration VLA imaging was obtained for programs AT206 (B, CS) and AT268 (D-array), conducted in 1997–1998 and 2001, respectively. During the AT206 (B, CS-array) observations, we employed a six-point mosaicking strategy, with the mosaic spacing and grid orientation designed to provide relatively uniform sensitivity over the high-brightness portion of the M 33 disk (r< 9.5 kpc, i.e., 39 arcmin). The D configuration observations (AT268) were obtained with the VLA antennas continuously scanning over a larger, rectangular region centered on M 33. A grid of 99 discrete phase tracking centers was utilized for this on-the-fly (OTF) observing mode. The aggregate visibility data from all three (B, CS, and D) VLA configurations is sensitive to HI structures at angular scales of ~4–900 arcsec. Early results based on the B and CS configuration data were presented in Thilker et al. (2002), while an independent reduction of all three VLA configurations has been presented by Gratier et al. (2010). The datacubes analyzed herein benefit from a combination of VLA observations with GBT total power data, and hence should be considered superior to all datasets published earlier, with the sole exception of Braun (2012) for which the dataset is identical.

Our single dish GBT observations were obtained during 2002 October, using the same setup described by Thilker et al. (2004) for M 31 observations. The telescope was scanned over a 5 × 5 deg region centered on M 33 while HI spectra were measured at a velocity resolution of 1.25 km s^-1. At M 33’s distance of 840 kpc, the GBT beam (9.́1 FWHM) projects to 2.2 kpc, and our total power survey covers a 74 × 74 kpc^-2 region. At a resolution of 3 kpc and 18 km s^-1, the achieved rms flux sensitivity was 7.2 mJy/beam, corresponding to an HI column density of 2.5 × 10¹⁷ cm^-2 in a single channel averaged over the beam. Joint deconvolution of all of the VLA interferometric data was carried out after inclusion of the appropriately scaled GBT total power images following the method employed for the reduction of the M 31 data described in detail by Braun et al. (2009).

The angular resolution of the HI products used for our study was substantially less than the limit of our VLA+GBT dataset. In practice, we used datacubes with 20 arcsec (81 pc) and 130 arcsec (0.53 kpc) restored beam size. Masking of these cubes was completed to limit our analysis to regions with sufficient signal-to-noise and mitigate contributions from foreground Milky Way emission. With the support from masked regions outside of M 33, we interpolated Galactic emission structures across M 33 and subtracted this component to account crudely for its presence.

From the masked 20′′ and 130′′ datacubes we produced moment-0, moment-1, peak intensity, and velocity at peak intensity images. These image products were generated in the usual manner with Miriad’s moment task. The peak intensity and velocity at peak intensity images are derived from a three-point quadratic fit. Finally, Miriad’s imbin task was used to generate rebinned cubes with only one pixel per resolution element, such that adjacent pixels were uncorrelated. These binned cubes were the products used to constrain the galaxy model.

Fig. 1 Upper panels: integrated 21 cm line intensity and the intensity weighted mean velocity over the whole galaxy using data at 130 arcsec resolution. Lower panels: peak intensity and the velocity at peak maps using VLA+GBT data at a resolution of 20 arcsec.

2.1. The molecular gas map

Rotational velocities in the innermost regions are complemented by the independent CO J = 1–0 line dataset described by Corbelli (2003) at a spatial resolution of 45 arcsec (FCRAO survey). For the H₂ surface mass density radial scalelength, we use the mean value between that determined by the Five College Radio Astronomy Observatory (FCRAO) survey and that determined via the IRAM-30 m CO J = 2–1 line map of Gratier et al. (2010) (1.9 and 2.5 kpc, respectively). We use the following expression to fit the molecular gas surface density: $\begin{array}{ccc}{\mathrm{\Sigma }}_{{\mathrm{H}}_{\mathrm{2}}}\mathrm{=}\mathrm{10}\mathrm{\times }{\mathrm{e}}^{\mathrm{-}\mathit{R}\mathit{/}\mathrm{2.2}} {\mathit{M}}_{\mathrm{\odot }}{\mathrm{pc}}^{-2}& & \end{array}$(1)where R is in kpc. The total gas mass in molecular form according to the above expression is 3.0 × 10⁸ M_⊙, if integrated through the whole HI extent. This is in agreement with the molecular mass inferred by Corbelli (2003), Gratier et al. (2010), and more recently by Druard et al. (2014).

The atomic and molecular hydrogen gas surface densities have been multiplied by 1.33 to account for the helium mass when deriving the total gas mass surface density.

3. The stellar mass surface density

In this section, we describe how the stellar mass surface density maps are derived: in Sect. 3.1, we present the method, then in Sects. 3.2 and 3.3, we describe the data and their processing and, finally, in Sect. 3.4, we present the results. The extension of the stellar mass distribution to distances larger than ~5 kpc is described in Sect. 3.5.

3.1. Method: from light to stellar mass

Maps of stellar mass surface densities can be derived from multiband optical imaging, as shown in the seminal work of Zibetti et al. (2009, ZCR09 hereafter), thanks to the relatively tight correlation existing between colors and the apparent stellar mass-to-light ratio (M/L, e.g., Bell & de Jong 2001). Here we implement an extension of the ZCR09 method, based on a fully Bayesian approach, similar to that adopted and validated by a number of works in the literature of the last decade (e.g., Kauffmann et al. 2003; Salim et al. 2005; Gallazzi et al. 2005; da Cunha et al. 2008; Gallazzi & Bell 2009). Our imaging dataset is composed of mosaic maps in the B, V, and I bands from the Local Group Survey (Massey et al. 2006) and g and i bands from the Sloan Digital Sky Survey (SDSS, York et al. 2000). After processing the images as described in the following subsections, we perform a pixel by pixel stellar population analysis by comparing the measured magnitudes in the five bands with the corresponding magnitudes derived from a comprehensive model library of 150 000 synthetic spectral energy distributions (SED). Each model SED is computed with a stellar population synthesis (SPS) technique based on the most recent revision of the Bruzual & Charlot (2003, hereafter BC03; Charlot & Bruzual, in prep.) simple stellar populations (SSPs), which use the “Padova 1994” stellar evolutionary library (see BC03 for details) and Chabrier (2003) stellar initial mass function (IMF). The SSPs are linearly combined according to a variety of star formation histories (SFH), parameterized and distributed as in Gallazzi et al. (2005, see their Sect. 2.3 for details), where random (in terms of time of occurrence, intensity, and duration) bursts of star formation are superimposed onto continuous exponentially declining SFHs, with random initial time t_form and e-folding timescales 1/γ. The metallicity of the stellar populations is fixed in each model along the SFH and is randomly chosen between 0.02 and 2.5 Z_⊙ from a uniform linear variate. Dust attenuation is also implemented following the simplified two-component parameterization as in à la Charlot & Fall (2000) for the diffuse ISM and the stellar birth clouds: we adopt the same prescriptions as in da Cunha et al. (2008, see their Sect. 2.1 for details). Each model is originally normalized so that the total mass in living stars and stellar remnants (white dwarfs, neutron stars, and black holes) at the end of the SFH is 1 M_⊙. At each pixel, for each model i we define the chi-square as follows: $\begin{array}{ccc}{\mathit{\chi }}_{\mathit{i}}^{\mathrm{2}}\mathrm{=}\sum _{\mathit{j}\mathrm{=}\mathrm{1}}^{\mathit{N}\mathrm{bands}}\frac{{\mathrm{(}{\mathit{m}}_{\mathrm{obs}\mathit{,j}}\mathrm{-}\mathrm{(}{\mathit{m}}_{\mathrm{model}\mathit{,i,j}}\mathrm{+}\mathit{Z}{\mathit{P}}_{\mathit{i}}{\mathrm{)}}^{\mathrm{)}}}^{\mathrm{2}}}{{\mathit{\sigma }}_{\mathit{j}}^{\mathrm{2}}}& & \end{array}$(2)where j runs over all bands used for the mass estimation, m_obs,j is the observed magnitude of the pixel in the jth band and σ_j the corresponding error, m_model,i,j is the magnitude of the model i in the jth band plus distance modulus, and ZP_i is the zeropoint offset to apply to the model to match the observations. We then compute the ZP_i that minimizes ${\mathit{\chi }}_{\mathit{i}}^{\mathrm{2}}$. Since the models are originally normalized to 1 M_⊙ in stars, the best-fitting stellar mass ${\mathit{M}}_{\mathit{i}}^{\mathrm{\ast }}$ associated with the model i is simply given by: $\begin{array}{ccc}{\mathit{M}}_{\mathit{i}}^{\mathrm{\ast }}\mathrm{=}{\mathrm{10}}^{\mathrm{(}\mathrm{-}\mathrm{0.4}\mathit{Z}{\mathit{P}}_{\mathit{i}}\mathrm{)}}{\mathit{M}}_{\mathrm{\odot }}\mathit{.}& & \end{array}$(3)To each model i we further associate a likelihood $\begin{array}{ccc}{\mathrm{ℒ}}_{\mathit{i}}\mathrm{\propto }\exp \mathrm{(}\mathrm{-}{\mathit{\chi }}_{\mathit{i,}\min }^{\mathrm{2}}\mathit{/}{\mathrm{2}}^{\mathrm{)}}\mathit{.}& & \end{array}$(4)We can now build the unity-normalized probability distribution function (PDF) for the stellar mass at a given pixel: we consider the median-likelihood value (i.e., the value of stellar mass corresponding to the cumulative normalized PDF value of 0.5) as our fiducial estimate. The mass range between the 16th and 84th percentile of the PDF divided by two provides the error on the estimate. The so-computed errors account for both photometric uncertainties and, importantly, also for intrinsic degeneracies between model colors and M/L.

The adoption of this Bayesian median-likelihood estimator, instead of the simple conversion of light to mass via the M/L derived from a relation with one or two colors as done in ZCR09, allows to take full advantage of our multiband dataset and reduce the impact of systematic effects in the individual bands (see next subsection).

We note that the BC03 revised models used here produce significantly higher stellar masses than those adopted in ZCR09, by some 0.1 dex. The main difference concerns the implementation of the TP-AGB evolutionary phase, which in the so-called “CB07” version predicts much higher luminosities than in the original BC03 models at fixed stellar mass. Recent works (e.g., Kriek et al. 2010; Conroy & Gunn 2010; Zibetti et al. 2013), however, have shown that the BC03 models are indeed closer to real galaxies than the CB07 models and this motivates our choice to use the BC03 models.

3.2. The SDSS data

The Sloan Digital Sky Survey (SDSS, York et al. 2000) has covered M 33 in several imaging scans, which must be properly background-subtracted and combined to obtain a suitable, science-ready map. We only use the g and i bands from SDSS, which deliver the best signal-to-noise ratio out to the optical radius. We discard the r band because of its contamination by Hα line emission (see discussion in ZCR09). We reconstruct the separate scan stripes that cover M 33 by stitching together the individual fields that we retrieved from the DR8 (Aihara & et al. 2011) data archive server. The SDSS pipeline automatically subtracts the background by fitting low order polynomials along each column in the scan direction. The process, however, is optimized to work with galaxies with sizes of a few arcminutes at most. In the case of M 33, which is more than 1 degree across, the standard background subtraction algorithm fails by subtracting substantial fractions of the galaxy light. Therefore, we have added back the SDSS background and obtained the original scans, where only the flatfield correction has been applied.

For each column along the scan direction, we determine the background as the best-fitting 5th degree Legendre polynomial, after masking out the elliptical region that roughly corresponds to the optical extent of M 33 (80 arcmin × 60 arcmin). This is achieved using the background task in the twodspec IRAF package, with a 2-pass sigma-clipping rejection at 3σ. This procedure leaves large-scale residuals in the backround regions, which we estimate at the level of ≈26.5 mag arcsec^-2 in g and ≈25.5–26 mag arcsec^-2 in i (r.m.s.). Residuals are visible in the final mosaic maps as structures along the scan direction. We also checked that the extrapolation of the fit over the masked central area does not introduce spurious effects. In fact, the amplitude of the fitted polynomial along the scan direction never exceeds values corresponding to ≈24.5 mag arcsec^-2 in the central regions, thus affecting their surface brightness by a few percent at most.

The scans are then astrometrically registered and combined in a widefield mosaic using SWARP (Bertin et al. 2002). Bright Milky Way stars must be removed before surface photometry can be obtained. Stars are identified using SExtractor (Bertin & Arnouts 1996) on an unsharp-masked version of the mosaic. Bright young stellar associations and compact star forming regions are manually removed, as our intent is to model the underlying (more massive) stellar disk, and such bright young features could bias the inferred M/L ratio downward. The original mosaic is then interpolated over the detected stars with a constant background value averaged around each star, using the IRAF task imedit. The size of the apertures for interpolation are optimized iteratively to minimize the residuals, and the brightest stars are edited manually for best results.

The edited SDSS mosaic is then smoothed with a Gaussian of ≈100 pc (25 arcsec) FWHM to remove the stochastic fluctuations due to bright stars occupying the disk of M 33, and finally resampled on a pixel scale of 8.3 arcsec/pix (i.e., ≈3 pixels per resolution element). The g- and i-band surface brightness maps are finally corrected for the foreground MW extinction, by 0.138 and 0.071 mag, respectively. Thanks to the massive smoothing and resampling, photon noise is largely negligible. The photometric errors are dominated by systematics in the photometric zero point and in the background subtraction. We conservatively estimate these sources of error as 0.05 mag for the zero-point (also including uncertainties related to the unknown variations of foreground Galactic extinction) and as 30% of the surface brightness for the background level uncertainty at μ_g = 25 and μ_i = 24.5 mag arcsec^-2, which takes the aforementioned fluctuations in the residual background into account. The two contributions then are added in quadrature pixel by pixel. The zero-point uncertainty dominates in the inner and brighter parts of the galaxy, while the background uncertainty dominates at low surface brightness.

Fig. 2 Panel a): color composite RGB image of the B − V − I mosaics from the Massey et al. (2006) survey. Panel b): stellar mass map derived from the median likelihood estimator (described in Sect. 3.1) based on LGS BVI bands. Panel c): stellar mass map obtained from the Bayesian median likelihood estimator, based on the full LGC and SDSS BVIgi dataset. Panel d): apparent stellar M/L derived from median likelihood stellar mass map based on LGS+SDSS BVIgi bands and the observed LGS V-band surface brightness.

3.3. The Local Group Survey

The Local Group Survey (Massey et al. 2006) provides a set of deep widefield images of M 33 in the UBVRI passbands (only B, V, and I will be used to estimate the stellar mass in the following), together with a photometric catalog of 146 622 stars. Because of restrictions in the field of view, the galaxy was covered by acquiring three partially overlapping fields (C,N,S) in each band, a total of 15 images. All these images, together with those secured in narrow bands, are publicly accessible¹ but no photometric calibration is provided because of the composite nature of the mosaicking CCD detector. In order to calibrate the images, we select from the photometric catalog a set of stars brighter than V = 16.00: 195 stars in the C field, 114 in the N, and 76 in the S. First, we measure the size of the PSF in the images and then degraded all of them to match the PSF of the image with the worst seeing (the N field in the I band with a Gaussian σ = 0.35 arcsec). Next, we perform multi-aperture photometry of the selected stars and derive the calibration by comparing fluxes with the values reported in the catalog. The multi-aperture photometry also provides the background offsets between the 3 different fields in each band. After readjusting the background and the scale of the N and S fields to C, the fields are combined to produce an image of the entire galaxy, which is later sky subtracted by using the most external regions. After registering the positions to the V-band reference, the images are again calibrated with the same procedure outlined above but including color terms. A full resolution RGB color composite image of the B − V − I bands is shown in Fig. 2, top left panel.

The B, V, and I band calibrated maps are then edited, degraded in resolution and resampled to match the SDSS images described in the previous subsection. From the calibration procedure of these final maps, again by comparison with the photometry of the stars with V ≤ 16.00 in the Local Group Survey catalog, we determine uncertainties in the zero-point of 0.02. 0.10, and 0.18 mag in the B, V, and I band, respectively. These final maps are also used for determining the background subtraction. Twenty fields containing no stars, each 20 × 20 arcsec in extent, were selected in the far NE and SW regions at the borders of the image; these fields show minimal emission in all three bands. In each bandpass, we assume the background to be the mean of the 20 median values, and the background uncertainty to be the 1σ of the sample (not of the mean); this uncertainty amounts to 30% of the surface brightness at μ_B = 25.0, μ_V = 24.3, μ_I = 23.5.

3.4. Stellar mass map(s)

The left bottom panel of Fig. 2 shows the map of stellar mass derived from the pixel-by-pixel analysis described in Sect. 3.1, using both the BVI maps from LGS and the gi maps from SDSS. The smoothness and azimuthal symmetry of the map is remarkable, although, especially at low surface density values, one can clearly see background residuals appearing as large-scale “spot” fluctuations (coming from the LGS maps) and fluctuations along stripes (coming from the SDSS scans). Random errors on this map (as derived from the interpercentile range of the PDF described in Sect. 3.1) range from 0.06 to 0.1 dex (approximately 15 to 25%), including the contribution of both photometric uncertainties and intrinsic model degeneracies between colors and M/L, but not systematics (see end of this section).

In the bottom right panel of Fig. 2, we display the map of the ratio between the estimated stellar mass and the V-band surface brightness (as measured) in solar units. This figure makes it obvious that a simple rescaling of the surface brightness by a constant mass-to-light ratio (as is usually done in the dynamical analysis) is a very crude approximation to the real stellar mass distribution. In particular, the M/L gradients apparent from this map substantially affect the slope of the mass profile, as discussed more in detail in Sect. 6.

Figure 2, with its many apparent artifacts, also highlights the limits of our photometry, mainly due to the difficulties of obtaining an accurate background subtraction in such a large-scale mosaic. On the other hand, the redundancy of our dataset in terms of wavelength coverage reveals its importance to limit the impact of systematic effects when we compare the mass surface density obtained from the full BVIgi dataset (bottom left panel of Fig. 2) to that obtained from the LGS dataset alone (upper right panel). In particular, the comparison shows that i) the anomalously high surface density that appear in the outskirts of the map based on LGS alone are corrected by also including SDSS imaging, most likely thanks to a more accurate background subtraction in the SDSS; ii) the north-south asymmetry evident in the LGS map is also corrected by a more stable and accurate zero-point determination in the SDSS. Moreover, we note that the stripe structure introduced by the background subtraction in the SDSS images is partly removed by the inclusion of LGS images.

Despite some differences in the overall normalization, the mass surface density variations in the two maps are consistent. The stellar mass map derived from the full dataset (LGS+SDSS) will be our reference mass distribution but we shall carry out the dynamical analysis using both maps, limiting the usage of the mass map from the LGS dataset to R< 4 kpc. Considering the uncertainties derived from our Bayesian marginalization approach (thus combining model degeneracies and random photometric errors) and comparing the two mass maps obtained with the two photometric datasets, we estimate a typical stellar mass uncertainty of ~30% (0.11 dex) for μ_B< 25, which we apply to the dynamical analysis. We stress once again that the uncertainties derived from the interpercentile range of the PDF over a large and comprehensive library of model star-formation histories, metallicity, and dust distributions, de facto realistically include the contribution of the systematic uncertainties related to the choice of models. This is one of the key advantages of using this approach with respect to other maximum-likelihood fitting algorithms, either parametric or non-parametric.

The only possible systematic uncertainties that are not included in the error budget are those related to the IMF, which is kept fixed in all models (to the Chabrier IMF), and those related to the basic ingredients of stellar population synthesis, i.e., the base SSPs. One can, of course, include models with different IMFs in the library, but broadband optical colors are essentially insensitive to the IMF, hence they do not provide any constraint in that respect. Possible different IMFs should thus be treated as an extra freedom in the stellar mass normalization. Testing different SSPs is clearly out of the scope of this paper. However, it is worth mentioning that after many years of debate about the role of TP-AGB stars in stellar populations, the community is reaching a broad consensus on their limited role (e.g., Kriek et al. 2010; Conroy & Gunn 2010; Zibetti et al. 2013), thus leaving little freedom on the intrinsic M/L value of simple stellar populations at a fixed SED. Therefore, this contribution to the mass error budget appears negligible in comparison with the variations induced by different star formation histories, metallicity and dust distributions, already accounted for by our method of error estimates.

3.5. The outermost stellar disk

As pointed out by several surveys (e.g., McConnachie et al. 2010), the stellar disk of M 33 extends out to the edge of the HI disk following a similar warped structure. To account for this fainter, but extended disk, we first extrapolate the stellar mass surface density of the inner stellar map outward to 10 kpc, using the radial scalelength inferred by the 3.6 μm map (Verley et al. 2007, 2009). This scalelength, 1.8 ± 0.1 kpc, is consistent with that inferred in our stellar surface density maps beyond R ~ 2 kpc. Stellar counts from deep observations of several fields further out, in the outer disk of M 33, have been carried out by Grossi et al. (2011) using the Subaru telescope. The typical face-on radial scalelength of the stellar mass density in the outer disk inferred from these observations is 25 arcmin (~6 kpc), much larger than that of the brighter inner disk and very similar to the HI scalelength in the same region. We shall use this radial scaling to extrapolate exponentially the stellar surface density further out and truncate the profile with a sharp falloff beyond 23 kpc.

4. The rotation curve

In this section, we derive the rotation curve of M 33 starting from an analysis of the spatial orientation of the disk.

4.1. The warped disk

To make a dynamical mass model of a disk galaxy it is necessary to reconstruct the three-dimensional velocity field from the velocities observed along the line of sight. If velocities are circular and confined to a disk, one needs to establish the disk orientation, i.e., the position angle of the major axis (PA), and the inclination (i) of the disk with respect to the line of sight. If the disk exhibits a warp these parameters vary with galactocentric radius. This is often the case for gaseous disks that extend outside the optical radius and show different orientations than the inner disks. Tilted ring models are used to infer the rotation curve of warped galactic disks. Corbelli & Schneider (1997) have fitted a tilted ring model to the 21 cm line data over the full extent of the M 33-HI disk. However, the disk was not fully sampled by the data since the observations (carried out with the Arecibo flat feed, FWHM = 3.9 arcmin) only sampled the disk with 4.5 arcmin spacing. To determine the disk orientation using the new 21 cm all-disk survey described in this paper, we follow a tilted ring model procedure that departs from the usual schemes and has been used by Corbelli & Schneider (1997) and later implemented by Józsa et al. (2007) (TiRiFiC, a Tilted Ring Fitting Code, available to the public). A set of free rings is considered, each ring being characterized by its radius R and by seven additional parameters: the HI surface density S_HI, the circular velocity V_c, the inclination i and the position angle PA, the systemic velocity V_sys, and the position and velocity shifts of the orbital centers with respect to the galaxy center (Δx_c,Δy_c). Importantly, rather than fitting the moments of the flux distribution, we infer the best-fitting parameters by comparing the synthetic spectra of the tilted ring model to the full spectral database, as explained in details in Appendix A. In Appendix A we described also the two minimization methods used, the “shape” and “v-mean” methods, which give consistent results for the M 33 ring inclinations and position angles.

We display the best-fitting values of i and PA in Fig. 3 with their relative uncertainties. As we can see from Fig. 3 the orientation of the outermost rings has high uncertainties. This is especially remarkable when using the “shape” method, hence, we fix the outermost ring parameters to be equal to those of second-last ring for deconvolving the data. The intermediate regions solutions appear to be more robust. In the next subsection the best-fitting values of systemic velocity shifts and ring center displacements (shown in Appendix A) will be used together with the rings orientation angles, i and PA, to derive the rotational velocities, V_r, and the face-on surface brightness, Σ_HI, from high and low resolution data. Notice that the values of V_r and Σ_HI we infer are consistent with the value of the free parameters V_c and S_HI of the adopted ring model but will be sampled at a higher resolution along the radial direction.

Fig. 3 Inclination and position angle of the 11 free tilted rings used for deconvolving the 21 cm data. The open triangles (in black in the online version), connected by a dashed line, shows the value using the shape minimization method, the open squares (in red in the online version), connected by a continuous line, are relative to the v-mean minimization method (see Appendix A for details). The uncertainties are computed by varying one parameter of one free ring at a time around the minimal solution until the χ2 of the tilted ring model increases by 2%.

4.2. The 21 cm velocity indicators

To derive the rotation curve, we use the 21 cm datacube at a spatial resolution of 20 and 130 arcsec. Emission in the high spatial resolution dataset is visible only out to R = 10 kpc while in the lower resolution dataset the HI is clearly visible over a much extended area. To trace the disk rotation, we consider both the peak and the flux weighted mean velocities (moment-1) along the line of sight at the original spectral resolution. The velocity at the peak of the line is in general a better indicator of the disk rotation than the mean velocity when the signal-to-noise is high and the rotation curve is rising (e.g., de Blok et al. 2008), and hence we will use this velocity indicator inside the optical disk. In the outer disk, where the signal-to-noise is lower and the rotation curve is flatter we use the flux weighted mean velocities (as we show below the two velocity indicators in this region give quite consistent rotation curves).

The rotation curve of M 33 is extracted from the adopted line-of-sight velocity indicator using as a deconvolution model the best-fit disk geometrical parameters derived via the shape and v-mean methods (see Appendix A). We shall refer to these as model-shape and model-mean parameters, they are shown in Figs. 3 and A.1. If 21 cm emission is present in an area located at a larger radius than the outermost free ring, we deconvolve the observed velocities assuming a galaxy disk orientation as the outermost free ring.

When using the mean instead of the peak velocity, the same deconvolution model results in a rotation curve with lower rotational velocities. In the optical disk, the difference between the rotation curve traced by the peak or the mean velocities can be as high as 10 km s^-1, but the scatter between the results relative to the two deconvolution models is small (it is always less than 2 km s^-1 and for most bins less than 1 km s^-1). Further out there are no systematic differences between the curve traced by the peak and the mean velocity, but the outermost part is very sensitive to the orientation parameters adopted. These parameters are not tightly constrained in the outermost regions because of the partial coverage of the modeled rings with detectable 21 cm emission, and as a result of some multiple component spectra. Here it will be important to consider the deconvolution model uncertainties in the rotation curve.

Fig. 4 Bottom panel: rotational velocities for each pixel in the peak map (small dots, red dots for the southern half and blue dots for the northern half in the online version) and average values in each radial bin (large symbols). The model-shape has been used to deconvolve the data. Triangles (red in the online version) are for the southern half while squares (blue in the online version) are for the northern half. Top panel: average rotation curve for peak velocities. Bins are 2 arcmin wide in the radial direction.

Fig. 5 Bottom panel: rotational velocities for each pixel in the first moment-1 map (small dots), the average values in each ring (large symbols, red symbols are for the southern half and blue symbols for the northern half in the online version). The model-mean has been used to deconvolve the data. Triangles are for the southern half while squares are for the northern half. The average rotation curve using moment-1 (mean) velocities is shown in the top panel. Bins are 2 arcmin wide in the radial direction.

In Figs. 4 and 5, we show the rotation curve from the peak and the mean velocities adopting the tilted ring model-shape and the model-mean, respectively. For the outer disk, we use the low resolution dataset with 2 arcmin radial bins (0.5 kpc wide). For the inner disk the high resolution data are binned in 40 arcsec radial bins (163 pc wide). In Fig. 6 we show the curve traced by peak velocities after averaging the rotational velocities relative to the deconvolution model-shape and deconvolution model-mean. In the same figure, the filled dots show the curve from the low resolution dataset, which perfectly matches the high resolution curve. To not overweight the inner rotation curve in the dynamical analysis, we will use the low resolution 21 cm dataset in addition to the velocities of the CO-J = 1–0 line peaks. The average values of PA and i between deconvolution model-shape and model-mean are used to trace the rotation curve with CO J = 1–0 line (star symbols in Fig. 6). These deconvolution parameters have also been used to trace the inner curve with moment-1 velocities (shown for comparison in Fig. 6 with open square symbols). There is not much difference between the rotational velocities retrieved from model-shape or model-mean in the inner region, but the moment-1 velocities clearly yield a lower rotation curve than the peak velocities.

Fig. 6 Inner rotation curve as traced by the peak velocities at high resolution (filled squares in blue in the online version), and at low resolution (filled dots in red in the online version). The CO data are shown with star symbols (in green in the online version). The open square symbols (in red in the online version) are for rotational velocities as traced by the moment-1 map. For each curve the weighted mean velocities between that relative to deconvolution model-shape and model-mean are used.

Fig. 7 Tilted ring position angles and inclinations (open square symbols and filled dot symbols respectively for the ROTCUR and KINEMETRY fitting routines).

Fig. 8 Some of the tilted ring parameters fitted by the routines ROTCUR and KINEMETRY. The top panel shows the term k5, which is a higher order term in the KINEMETRY outputs. In the bottom panel we show the shifts of the systemic velocity, and, in the middle panel, the orbital center shifts along the y-axis in the south-north direction (filled square symbols) and along the x-axis oriented from east-to-west (open square symbols).

Fig. 9 Adopted rotation curve (open black circles) for the inner regions (bottom panel) using the peak velocities and for the outer regions (top panel) using the moment-1. The dispersions take into account variation of Vr between different deconvolution models, the 2% variations around the minimum χ2, and the standard deviation of the mean velocity in each bin. In both panels, the open square and asterisk symbols (in magenta and green in the online version) show Vr according to the ROTCUR and KINEMETRY deconvolution model, respectively. In the top panel the open triangles are for the moment-1 velocities deconvolved according to model-shape and model-mean (in red and blue colors, respectively in the online version).

4.3. Rotation curve: a comparative approach

In this subsection, we compare the rotation curve derived through our technique of fitting the 21 cm line data with a tilted ring model, to those resulting from other approaches and methods. These complementary approaches are useful to strengthen the robustness of the rotation curve and to better define the uncertainties. We consider two additional methods devised for deriving galaxy rotation curves from two-dimensional moment maps. The first method is the standard least-square fitting technique developed by Begeman (1987) as implemented in the ROTCUR task within the analysis package NEMO (Teuben 1995). The other method is based on the harmonic decomposition of the velocity field along ellipses and for this purpose we use the software KINEMETRY developed and provided by Krajnović et al. (2006). Both methods work on 2D momentum maps, i.e., on one velocity per pixel, and minimize the free parameters of one ring at a time, without accounting for possible overlaps of ring pieces in the beam. We run ROTCUR in two steps: at first we let the ring centers and V_sys vary. The orbital center shifts are larger in the outermost regions and roughly in the direction of M 31, as found by Corbelli & Schneider (1997), and are shown in Fig. 8. In the second iteration, we fix the ring centers and V_sys to the average values found in the first iteration, x_c = −33 arcsec y_c = 105 arcsec V_sys = −178 km s^-1. We run a second iteration because by fixing the orbital centers the routines converge further out, at radii as large as R = 23.5 kpc. In this second attempt we run ROTCUR with 51 free rings, uniform weight, excluding data within a 20° around the minor axis. The resulting PA and i are shown in the top panel of Fig. 7. With the KINEMETRY routines, we fit with a smaller number of free rings, 10, and convergence is achieved over a smaller radial range, out to R = 17.5 kpc. Higher order terms of the harmonic decomposition are useful tools when looking for non-circular motion in the disk. In the bottom panel of Fig. 7, we show the usual PA and i but in addition we plot the coefficient k₅ in Fig. 8. This is a higher order term that represents deviations from simple rotation due to a separate kinematic component like radial infall or non-circular motion. The amplitude of k₅ is negligible in the optical disk but it reaches 8 km s^-1 in the outer disk. If the anomalous velocities are in the radial direction, this gas can fuel star formation in the inner disk.

Figure 9 shows that the rotation curves obtained with the NEMO and KINEMETRY routines are very similar to that derived from the deconvolution model-shape and -mean inside the optical disk. There are some differences in the outer disk because of differences in the resulting warp orientation. These differences become significant beyond 15 kpc, and we take them into account.

4.4. The adopted rotation curve and its uncertanties

For the final rotation curve of M 33 we use the following recipe: inside the optical disk we use 21 cm peak velocities and compute the circular velocities as the weighted mean of velocities originating from our two deconvolution choices (model-shape and model-mean). To these we append the CO dataset. Beyond 9 kpc, in each radial bin we use the weighted mean circular velocity computed by averaging the 21 cm moment-1 velocities resulting from: (1) the deconvolution model-shape; (2) the deconvolution model-mean; (3) the package KINEMETRY; and (4) the task ROTCUR in NEMO package. By extrapolating the orientation of the outermost fitted ring for a few kpc outward, we trace the rotation curve out to 23 kpc.

The rotation curve uncertainties not only take the standard deviations around the mean of the 2 or 4 velocities averaged into account, but also the data dispersion in each radial bin and the uncertainties relative to deconvolution models (ΔV_r corresponding to 2%χ² variations of the tilted ring model-shape and -mean through PA and i displacements). We apply to the rotation curve the small finite disk thickness corrections, described in Appendix B. The final rotation curve is given in Table 1.

5. The surface density of the baryons

In this section, we use the radial profile of the stellar mass surface density from the BVIgi maps, together with the gaseous surface density, to compute the total baryonic face-on surface mass density of M 33 and the stellar mass-to-light ratio. The surface density of stars perpendicular to the disk is shown in Fig. 10: it drops by more than 3 orders of magnitudes from the center to the outskirts of the M 33 disk. According to the BVIgi-mass map, the total stellar mass is 4.9 × 10⁹ M_⊙ (5.5 × 10⁹ M_⊙ for the BVI map) of which about 12% resides in the outer disk. At the edge of the optical disk Barker et al. (2011) find a surface density of stars ~1.7 M_⊙/pc² for a Chabrier IMF down to 0.1 solar masses, assuming the inclination of our tilted ring model. This is consistent to what our modeled radial stellar distribution predicts at 9 kpc: 1.7 ± 0.5M_⊙/pc². The outermost field (S2) of Barker et al. (2011) has been placed outside the warped outer disk and hence traces only a possible stellar halo.

Using the best-fitting tilted ring model we derive the radial distribution of neutral atomic gas, perpendicular to the galactic plane, in the optically thin approximation. This is shown in Fig. 10. The total HI mass computed by integrating the surface density distribution in Fig. 10 out to 23 kpc is about 20% higher than the true HI mass of the galaxy, which is 1.53 × 10⁹ M_⊙. This is because we average the flux of all pixels with nonzero flux in each ring and the HI emission in the outermost rings does not cover the whole ring surface. We do this because it is likely that undetected pixels at 21 cm are not empty areas but host ionized gas. In fact, a sharp HI edge has been detected in this galaxy as the gas column density approaches 2 × 10¹⁹ cm^-2 and interpreted as an HI–>HII transition (Corbelli & Salpeter 1993). If we consider the irregular outer HI contours of M 33 as due to ionization effects, the outermost part of the atomic radial profile is in reality an HI+HII profile since there is HII at a similar column density where HI is lacking. This assumption has a negligible effect on the dynamical analysis of the rotation curve. The continuous line in Fig. 10 is the log of the function used to compute the dynamical contribution of the hydrogen mass to the rotation curve. The total baryonic surface density is computed adding to the hydrogen gas the stellar mass surface density, the molecular and the helium mass surface density (as given in Sect. 2.2). Stars dominate the potential in the star-forming disk (R< 7 kpc), beyond this radius stars and gas give a similar contribution to the baryonic mass surface density and decline radially in the outer disk with a similar scalelength.

Fig. 10 Neutral atomic gas surface density perpendicular to the galactic plane of M 33 (small filled triangles) and the function that fits the data (continuous line,red in the online version) after the 21 cm line intensity has been deconvolved according to tilted ring model-shape. Asterisk symbols indicate the stellar mass surface density using the BVIgi stellar surface density map. The dashed line is the fit to the stellar surface density and the extrapolation to larger radii. Open squares(in blue in the online version) show for comparison the surface density using the BVI mass map. The heavy weighted line is the total baryonic surface density, the sum of atomic and molecular hydrogen, helium and stellar mass surface density.

Fig. 11 Mass-to-light ratio as a function of galactocentric radius according to the BVIgi stellar surface density map (filled triangles, in black in the online version). This is computed by averaging the stellar mass map and the V-band surface brightness along ellipses corresponding to the tilted ring model. We then compute M∗/LV as the ratio between these two quantities in units of M⊙/L⊙. We show for comparison the same ratio relative to the BVI stellar surface density map (open triangles, in red in the online version). The error bars are the standard deviations relative to azimuthal averages and do not take into account uncertainties in the mass determination.

The mass-to-light ratio as a function of galactocentric radius is computed by averaging the BVIgi stellar surface density map and the V-band surface brightness along ellipses corresponding to the tilted ring model. The filled triangles in Fig. 11 show the ratio between these quantities, M_∗/L_V, in units of M_⊙/L_⊙. We show for comparison the same ratio relative to the BVI stellar surface density map (open triangles). Despite some differences in the two maps, the shape of the azimuthally averaged radial profiles of the mass surface density are consistent. There is some difference in the radial decline only in the innermost 1 kpc. The error bars are standard deviations relative to azimuthal averages and do not take uncertainties in the mass determination into account. We clearly see radial variations of the mass-to light ratio, especially in the innermost 1.5 kpc relative to areas at larger galactocentric radii. This is an effect of the inside-out disk formation and evolution, which adds to more localized variations present in the map such as in arm-interarm contrast.

Average values or radial dependencies of the mass-to-light ratios, shown in Fig. 11, can be compared with those derived previously using different methods. The dynamical analysis of planetary nebulae in M 33 (Ciardullo et al. 2004) yields mass-to-light ratios that increase radially, being 0.2 M_⊙/L_V⊙ at R = 2 kpc and 0.8 M_⊙/L_V⊙ at 5 kpc. However, this result is based on the assumption of radially constant vertical scaleheight and of a much longer radial scalelength of the stellar surface density than suggested by near-infrared photometry. These assumptions also imply a strong radial decrease of the velocity dispersion if the disk is close to marginal stability, in disagreement with what is observed for the atomic gas. A more recent revised model for the marginally stable disk of M 33 (Saburova & Zasov 2012) implies a flatter mass-to-light ratio, of order 2 M_⊙/L_V⊙ in closer agreement with what we find. An increase of the mass-to-light ratio in the central regions of M 33 and a stellar mass surface density with very similar values to what we show in Fig. 10 has been found by Williams et al. (2009) using resolved stellar photometry and modeling of the color-magnitude diagrams.

The former M 33 rotation curve (Corbelli 2003) was best fitted using CDM dark halo models by assuming an average stellar mass-to-light ratio in the range 0.5–0.9 M_⊙/L_V⊙, a somewhat smaller value than derived here. However, Corbelli (2003) considered an additional bulge component whose presence has not been supported by a subsequent analysis (Corbelli & Walterbos 2007). From the stellar mass model presented here it seems more likely that the M 33 disk has a larger mass-to-light ratio in the innermost 1.5 kpc rather than a genuine bulge. A constrained, radially varying, mass-to-light ratio without a bulge component, as given by our mass map, is more consistent with CDM halo models for M 33 than with other dark matter distributions (see next section).

6. Tracing dark matter via dynamical analysis of the rotation curve

The dynamical analysis of a high resolution rotation curve, such as that presented in this paper for M 33, together with detailed maps of the baryonic mass components, provides a unique test for the dark matter halo density and theoretical models of structure formation and evolution. In this section, we shall consider a spherical halo with an NFW dark matter density profile (Navarro et al. 1996, 1997) for galaxies forming in a cold dark matter scenario. We also consider the Burkert dark matter density profile (or core model) (Burkert 1995) since this successfully fitted the rotation curve of dark matter dominated dwarf galaxies (e.g., Gentile et al. 2007,and references therein). Both models describe the dark matter halo density profile using two parameters that we determine through the best fit to the rotation curve. Our last attempt will be to fit the rotation curve using MOdified Newtonian Dynamics (MOND).

We perform a dynamical analysis of the rotation curve in the radial range: 0.4 ≤ R ≤ 23 kpc. For the gas and the stellar surface density distribution, we consider the azimuthal averages shown in Fig. 10. Given the 30% uncertainty for the stellar mass surface densities we shall consider total M_∗ in the following intervals: 9.57 ≤ log M_∗ ≤ 9.81 M_⊙ when using the BVIgi mass map and 9.62 ≤ log M_∗ ≤ 9.86 M_⊙ when using the BVI mass map. We assume a vertically uniform gaseous disk with half thickness of 0.5 kpc and a flaring disk for the stellar disk with a half thickness of 100 pc at the center, reaching 1 kpc at the outer disk edge. The contribution of the disk mass components to the rotation curve is computed according to Casertano (1983), who generalizes the formula for the radial force to thick disks. We use the reduced chi-square statistic, χ², to judge the goodness of a model fit.

Fig. 12 Rotation curve of M 33 (filled triangles) and best-fitting dynamical models (solid line) for NFW dark matter halo profiles (dot-dashed line). The small and large dashed lines (red and blue in the online version) show the gas and stellar disk contributions to the rotation curve, respectively. In the bottom panel, the highest stellar contribution and lowest dark halo curve are for the best-fitting dynamical model using the BVIgi mass map, the other curves are for the BVI mass map. The top panel shows the most likely dynamical models for the two mass maps when the likelihood of the dynamical fit is combined with that of the stellar mass surface density and of the C − Mh relation resulting from simulations of structure formation in a ΛCDM cosmology.

6.1. Collisionless dark matter: a comparison with LCDM simulations

The NFW density profile is usually written as: $\begin{array}{ccc}\mathit{\rho }\mathrm{\left(}\mathit{R}\mathrm{\right)}\mathrm{=}\frac{{\mathit{\rho }}_{\mathrm{NFW}}}{\frac{\mathit{R}}{{\mathit{R}}_{\mathrm{NFW}}}(\mathrm{1}\mathrm{+}\frac{\mathit{R}}{{\mathit{R}}_{\mathrm{NFW}}}{)}^{\mathrm{2}}}& & \end{array}$(5)where ρ_NFW and R_NFW are the characteristic density and scale radius, respectively. Numerical simulations of galaxy formation find a correlation between ρ_NFW and R_NFW, which depends on the cosmological model (e.g., Navarro et al. 1997; Avila-Reese et al. 2001; Eke et al. 2001; Bullock et al. 2001). Often this correlation is expressed using the concentration parameter C ≡ R_vir/R_NFW, and V_h or M_h, the characteristic velocity or mass at R_vir. The virial radius, R_vir, is the radius of a sphere containing a mean density Δ times the cosmological critical density. The value of Δ varies between 93 and 97 depending upon the adopted cosmology (Macciò et al. 2008). This corresponds to a mean halo overdensity of about 360 with respect to the cosmic matter density at z = 0. In this paper we compare our best-fitted parameters C and M_h with the results of N-body simulations in a flat ΛCDM cosmology using relaxed halos for WMAP5 cosmological parameters (Macciò et al. 2008). In particular, we shall refer to the relation between the mean concentration and virial mass of dark halos resulting from the numerical simulation of Macciò et al. (2008) whose dispersion is ±0.11 around the mean of log  C. The resulting C − M_h relation is similar to that found by Macciò et al. (2007) and more recently by Prada et al. (2012).

We now fit the M 33 rotation curve using the dark halo concentration, C, and mass, M_h as free parameters. We fix the stellar mass surface density distribution to that given by our mass maps extrapolated outward, as explained in the previous section. We allow a scaling factor to account for the model and data uncertainties, i.e., we consider total stellar masses log  M_∗/M_⊙ = 9.69 ± 0.12 and log  M_∗/M_⊙ = 9.74 ± 0.12 for the BVIgi and BVI mass distributions, respectively. The best fits using the two mass maps are very similar as shown in the bottom panel of Fig. 12. The reduced χ² are 0.96 and 1.08 for the BVI and BVIgi stellar mass distributions, respectively, with log M_∗/M_⊙ = 9.8, close to the upper limit of our considered range. The dark matter halo for the best fits has concentration C = 6.7 and mass M_h = 5 × 10¹¹h^-1 M_⊙. A close inspection of the 1,2,3-σ confidence areas in the log C − log M_h/h^-1 plane and in the log  M_∗ − log  M_h/h^-1 plane, shown in Fig. 13, reveals that indeed there is still some degeneracy in the C − M_h plane: a lighter stellar disk can provide good fits to the rotation curve if a less massive dark halo with a higher concentration is in place. Confidence areas are traced in Fig. 13 only for the allowed stellar mass range. In Fig. 13, we also show the most likely value of the stellar mass computed by the synthesis models (dot-dashed line in the upper panels) and the log C − log M_h/h^-1 relation as from ΛCDM numerical simulations (continuous line Macciò et al. 2008) and its dispersion around the mean (dashed lines). The values of C and M_h suggested by the dynamical analysis of the M 33 rotation curve are in good agreement with the ΛCDM predictions.

We now compute the most likely values of the free parameters, C, M_h, and M_∗ by considering the composite probability of three events: i) the dynamical fit to the rotation curve; ii) the stellar mass determined via synthesis models; and iii) the log C − log M_h/h^-1 relation found by numerical simulations of structure formation in a ΛCDM cosmology. The composite probability gives smaller confidence areas in the free parameter space, shown in Fig. 14, which are very similar for the two mass maps. The model with the highest probability has the following parameters: ${\mathit{M}}_{\mathrm{\ast }}\mathrm{=}\mathrm{4.}{\mathrm{9}}_{-0.7}^{\mathrm{+}\mathrm{0.5}}\mathrm{\times }{\mathrm{10}}^{\mathrm{9}}{\mathit{M}}_{\mathrm{\odot }}$, ${\mathit{M}}_{\mathrm{h}}\mathrm{=}\mathrm{3.}{\mathrm{9}}_{-0.6}^{\mathrm{+}\mathrm{1.0}}\mathrm{\times }{\mathrm{10}}^{\mathrm{11}}{\mathit{M}}_{\mathrm{\odot }}$, C = 10 ± 1, for the BVIgi mass map and dimensionless Hubble parameter h = 0.72. The fit to the rotation curve, shown in the upper panel of Fig. 12, has a χ² = 1.3. For the BVI map, we get ${\mathit{M}}_{\mathrm{\ast }}\mathrm{=}\mathrm{4.}{\mathrm{8}}_{-0.4}^{\mathrm{+}\mathrm{0.4}}\mathrm{\times }{\mathrm{10}}^{\mathrm{9}}{\mathit{M}}_{\mathrm{\odot }}$, ${\mathit{M}}_{\mathrm{h}}\mathrm{=}\mathrm{4.}{\mathrm{6}}_{-0.6}^{\mathrm{+}\mathrm{0.7}}\mathrm{\times }{\mathrm{10}}^{\mathrm{11}}{\mathit{M}}_{\mathrm{\odot }}$, C = 9 ± 1 with a χ² = 1.1. The dynamical fit of this model is shown in the upper panel of Fig. 12. Given the marginal differences in the two sets of free parameters we can summarize the best-fit ΛCDM dark halo model for M 33 as follows: $\begin{array}{ccc}{\mathit{M}}_{\mathrm{h}}\mathrm{=}\mathrm{(}\mathrm{4.3}\mathrm{\pm }\mathrm{1.0}\mathrm{)}\mathrm{\times }{\mathrm{10}}^{\mathrm{11}}{\mathit{M}}_{\mathrm{\odot }} \mathit{C}\mathrm{=}\mathrm{9.5}\mathrm{\pm }\mathrm{1.5}& & \end{array}$(6)and the total stellar disk mass estimate is M_∗ = (4.8 ± 0.6) × 10⁹ M_⊙. The resulting dynamical model implies that the contribution of the dark matter density to the gravitational potential is never negligible, although it becomes dominant outside the star-forming disk (R> 7 kpc) where the stellar and gaseous disk gives a small, but similar contribution to the rotation curve.

Fig. 13 Confidence areas in the log C − log Mh/h-1 plane and in the M∗ − Mh/h-1 plane for the dynamical fit to the rotation curve using a NFW dark matter halo profile: 68.3% (darker regions, blue in the online version), 95.4% (magenta) and 99.7% (lighter regions, cyan). The mass distribution is that given by the BVIgi colors for the areas shown in a) and b), and by the BVI colors for the areas in c) and d). In panels a) and c) we show the C–Mh relation (solid line) and its dispersion (dashed line) for ΛCDM cosmology. In panels b) and d) the dash-dot line indicates the most likely stellar mass according the synthesis models, which rely on the colors shown next to panel labels.

Fig. 14 68.3% (darker regions, blue in the online version), 95.4% (magenta) and 99.7% (lighter regions, cyan) confidence areas in the log C − log Mh/h-1 plane and in the M∗ − Mh/h-1 plane when the probability of the dynamical fit is combined with the probability of the stellar mass surface density and of the C − Mh distribution from simulations of structure formation in a ΛCDM cosmology. The mass distribution is that given by the BVIgi map for the areas shown in a) and b), and by the BVI map for the areas in c) and d).

In Table 1 we display the rotational velocities, V_r and their dispersions, the azimuthal averages of the atomic gas surface mass density, Σ_HI, and of the modelled stellar mass density, Σ_∗. The values of Σ_∗ given in Table 1 correspond to the most likely value of the stellar mass distribution according to BVIgi maps and to rotation curve fit for ΛCDM halo models.

6.2. Core models of dark matter halos

The dark matter halo density profile proposed by Burkert (1995) is a profile commonly used to represent the family of cored density distributions (Donato et al. 2009) and it is given by: $\begin{array}{ccc}\mathit{\rho }\mathrm{\left(}\mathit{R}\mathrm{\right)}\mathrm{=}\frac{{\mathit{\rho }}_{\mathit{B}}}{\mathrm{(}\mathrm{1}\mathrm{+}\frac{\mathit{R}}{{\mathit{R}}_{\mathit{B}}}\mathrm{)}(\mathrm{1}\mathrm{+}\mathrm{(}\frac{\mathit{R}}{{\mathit{R}}_{\mathit{B}}}{\mathrm{)}}^{\mathrm{2}})}& & \end{array}$(7)where ρ_B is the dark matter density of the core that extends out to R_B (core radius). The baryonic contribution to the M 33 curve is declining between 5 and 10 kpc. Hence, to have an independent estimate of the dark matter density distribution according to Salucci et al. (2010) we must probe regions that are beyond 10 kpc. The extended rotation curve of M 33, which increases by about 10–20% between 5 kpc and 23 kpc, the outermost probed radius, is therefore appropriate for this purpose. We have searched the best-fitting parameters to the rotation curve, ρ_B, R_B, and M_∗ (with M_∗ within the range allowed by our stellar surface density maps). There are no acceptable fits if the BVIgi mass map is used, since the value of the minimum χ² is 3.2. We find a minimum χ² of 1.34 for the BVI mass map with M_∗ = 7.2 × 10⁹ M_⊙ and the following dark halo parameters: $\begin{array}{ccc}{\mathit{R}}_{\mathit{B}}\mathrm{=}\mathrm{7.5}\mathrm{kpc} {\mathit{\rho }}_{\mathit{B}}\mathrm{=}\mathrm{0.018}{\mathit{M}}_{\mathrm{\odot }}{\mathrm{pc}}^{-3}\mathit{.}& & \end{array}$(8)In this galaxy we are able to probe the dark matter out to three times R_B. The best fit is shown in the bottom panel of Fig. 15 and the 95.4% and 99.7% confidence areas in Fig. 16. Noticeably, the core model provides acceptable fits when the stellar mass is at the upper boundary of the interval compatible with the stellar surface density map. A good quality fit for the BVIgi stellar mass profile with a halo core model requires lower core densities and about 10¹⁰M_⊙ of stars, which is extreme for a galaxy like M 33. This massive stellar disk implies a factor 2 higher stellar surface density than that computed here via synthesis models, unless the IMF has a Salpeter slope over the entire stellar mass range. However, despite the heavy stellar disk predicted by the cored dark matter density distribution, the central surface density of the halo core is compatible with the value found by Donato et al. (2009), log  (ρ_BR_B/(M_⊙ pc²)) = 2.15 ± 0.2, by fitting a very large sample of rotation curves for galaxies of any luminosity and Hubble type. We plan to investigate cored dark matter halos in more detail in a subsequent, dedicated paper.

Fig. 15 Rotation curve of M 33 (filled triangles) and the best-fitting model (solid line) for a dark matter halo with a constant density core (bottom panel) and according to MOND (top panel). The small and large dashed lines (red and blue in the online version) show the gas and stellar disk Newtonian contributions to the rotation curves, respectively; the dark halo contribution in the bottom panel is shown with a long dashed line.

6.3. Non-Newtonian dynamics without dark matter

An alternative explanation for the mass discrepancy in galaxies has been proposed by Milgrom by means of the modified Newtonian dynamics or MOND (Milgrom 1983). Outside the bulk of the mass distribution, MOND predicts a much slower decrease of the (effective) gravitational potential, with respect to the Newtonian case. This is often sufficient to explain the observed non-Keplerian behavior of RCs (Sanders & McGaugh 2002). According to this theory, the dynamics becomes non-Newtonian below a limiting acceleration value, a₀, where the effective gravitational acceleration takes the value g = g_n/μ(x), with g_n the acceleration in Newtonian dynamics, x = g/a₀, and μ(x) is an interpolating function between the Newtonian regime and the case g ≪ a₀. Here, we shall use the critical acceleration value a₀ derived from the analysis of a sample of rotation curves a₀ = 1.2 ± 0.3 × 10^-8 cm s^-2 (Sanders & McGaugh 2002; Gentile et al. 2011). We have tested MOND for two choices of the interpolating function μ(x) (see Famaey & Binney 2005,for details). In particular, we have used the “standard” and the “simple” interpolation function and found that the former provides better fits to the M 33 data. Using the “standard” interpolating function, $\mathit{\mu }\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{=}\mathit{x}\mathit{/}\sqrt{\mathrm{1}\mathrm{+}{\mathit{x}}^{\mathrm{2}}}$, and the BVIgi stellar surface density map, we find a minimum χ² = 1.77 just outside the 3σ confidence limits. For the best fit, MOND requires M_∗ = 6.5 × 10⁹ M_⊙ and a₀ = 1.4 × 10^-8 cm s^-2. A slightly lower χ² (1.74) is found using the BVI mass map with M_∗ = 5 × 10⁹ M_⊙ and a₀ = 1.5 × 10^-8 cm s^-2. The fit provided by MOND does not improve considerably by increasing the stellar mass or the value of the critical acceleration a₀.

Fig. 16 The 95.4% (magenta in the online version) and 99.7% (lighter regions, cyan) confidence areas in the log r0 − log ρ0 plane and in the log r0 − log M∗ plane relative to the rotation curve fit using a core model and the BVI mass map. In the top panel, the dash-dot line indicates the most likely stellar mass according the synthesis models, which rely on the BVI colors. There are no core models in the 68.3% confidence area for the BVI mass map and no acceptable core models for the BVIgi mass map down to the 99.7% level. The solid line in the bottom panel indicates the inverse log ρ0 − log r0 correlation found by Donato et al. (2009).

6.4. The baryonic fraction

The rotation curve of M 33 is well fitted by a dark matter halo with a NFW density profile and a total mass of (4.4 ± 1.0) × 10¹¹ M_⊙. This is much larger than the halo mass, which in M 33 gives a baryonic fraction equal to the cosmic value, f_c = 0.16 (Spergel et al. 2003). Taking the best-fitting value of stellar mass, we compute a baryonic fraction of about 0.022. A baryonic fraction lower than the cosmic value is of no surprise since this is a common results in low luminosity galaxies. Feedback from star formation, such as supernovae driven outflows, is likely responsible for such a low baryonic fraction. Intergalactic filaments are in fact enriched of metals thanks to the exchange of matter between galaxies and their environment. Since the loss of gas from the galaxy will depend on the depth of the potential well, we expect that low-mass halos will be more devoid of baryons. The relationship between the stellar mass of galaxies and the mass of the dark matter halos has been derived by a statistical approach matching N-body simulated halo abundances as a function of the halo mass, to the observed abundance of galaxies as a function of their stellar mass (Moster et al. 2010). The resulting stellar mass fraction is an increasing function of the halo mass up to 10¹²M_⊙. The analysis of the rotation curves is a different way of testing this scenario. The dark halo and stellar masses resulting from the dynamical analysis of the M 33 rotation curve are compatible with the statistical relationship found by Moster et al. (2010, see their Fig. 6) when scatter is taken into account. A related question is whether the presence of outflows in the early evolutionary phases of M 33 might have affected the dark matter NFW profile. The recent work by Di Cintio et al. (2014b,a) has shown that the M 33 halo mass is just at the edge of where its inner density profiles is expected to be modified by baryonic feedback, with a cuspy-like preferred inner slope.

7. Summary and conclusions

The advantage of studying a nearby galaxy such as M 33 is the possibility of combining high resolution 21 cm datasets with a overwhelming amount of multifrequency data. In this work we took advantage of existing widefield optical images in various bands to construct a map of stellar mass surface density. Two different sets of images, namely the Local Group Survey (Massey et al. 2006) and the SDSS (York et al. 2000), have been combined in the innermost 5 kpc. Further out, these optical images are not sensitive enough to constrain the stellar surface density of M 33 against variations in the background light. Thus, we make use of the Spitzer 3.6 μm map, and of deep images of outer disk fields made with large optical telescopes, to estimate the stellar surface density out to the edge of the HI map. Using several methods for estimating the radially varying spatial orientation of the M 33 disk, we trace the rotation curve out to R = 23 kpc. The stellar and atomic gas maps, together with the available informations on the molecular gas distribution, have been used to derive the baryonic surface density perpendicular to the galactic plane with small uncertainties. The knowledge of the potential well due to the baryons, constrains the dynamical analysis of the rotation curve and the dark matter halo models in a more stringent way than previously possible. The radial distribution of the stellar mass surface density inferred from the maps, and the comparison with the light distribution, emphasizes the importance of combining the dynamical analysis with synthesis models. In fact, the mass-to-light ratio has non-negligible radial variations in the mapped region; additionally, local variations are present in the disk such as between the arm and interarm regions.

Numerical simulations of hierarchical growth of structure in a ΛCDM cosmological model give detailed predictions of the dark matter density distribution inside the halos. The universal NFW radial profile provides an excellent fit to the M 33 rotation curve. The free parameters, halo concentration and halo mass, are found to be C = 9.5 ± 1.5 and M_h = 4.3 ± 1.0 × 10¹¹ M_⊙, when we take the C − M_h relation resulting from numerical simulations (Macciò et al. 2008) and the stellar mass surface density via synthesis models into account. The best estimate of the stellar mass of M 33 is M_∗ = 4.8 × 10⁹ M_⊙, with 12% residing in the outer disk. When added to the gas this gives a baryonic fraction of order of 0.02. A comparison of this baryonic fraction with the cosmic inferred value suggests an evolutionary history that should account for a loss of a large fraction of the original baryonic content.

A naive view of the distribution of the baryons inside the NFW dark matter halo of M 33, as modeled in this paper, and the rotation curve are shown in Fig. 17.

The baryonic matter distribution in the framework of the modified Newtonian dynamic (MOND) does not provide good fits to the M 33 data once the stellar content is constrained. A dark halo with a constant density core is marginally compatible with the stellar mass distribution and with the dynamical analysis of the M 33 rotation curve. It requires a heavy stellar disk at the limit of the range allowed by our mass maps. The presence of a dark cuspy core in M 33, as predicted by structure formation in a ΛCDM hierarchical universe, is in agreement with numerical simulation of baryonic feedback effects on the density profile of dark matter (Di Cintio et al. 2014b). For our best fit halo mass, energy from stellar feedback is insufficient to significantly alter the inner dark matter density, and the galaxy retains a cuspy profile. Given the uncertainty in the M 33 stellar-to-dark matter ratio, however, we cannot exclude a slight modification of the original profile.

Fig. 17 A naive view of M 33: inner stellar surface density (from red to magenta colors in the online version) superimposed on the HI gas distribution in log-scale (light pink) are shown together with the stellar+gas contribution to the rotation curve (dashed cyan line) and with observed rotation curve (filled dots, in blue). The analysis presented in this paper has pinned down the average dark matter density distribution in the halo (light gray) which explains the observed rotation curve. The figure extrapolates the halo beyond the region analyzed in this paper in a naive smooth fashion. The total mass model fit to the rotational data, through the region analyzed in this paper, is shown by the continuous blue line.

Having determined the stellar mass of M 33 and considering a star formation rate of 0.45 M_⊙ yr^-1 (Verley et al. 2009) at z = 0, we estimate a specific star formation rate (SSFR) of about 10^-10 yr^-1, in agreement with the value inferred for large galaxy samples at z = 0 and by multi-epoch abundance matching models (e.g., Moster et al. 2013). A future analysis of the map of the SSFR of M 33, in the framework of a chemical evolution model, which takes into account gas outflow and inflow and reproduces the observed metallicity gradients (Magrini et al. 2010), will provide useful insights on the star formation history of M 33 and, in general, on galaxy evolutionary models.

Online material

Appendix A: The tilted ring model

To better constrain the disk orientation using the VLA+ GBT 21 cm datacube we first spatially smooth the data described in the previous section at 130 arcsec resolution to gain sensitivity. At this spatial resolution, we reach a brightness sensitivity of 0.25 K. Considering a typical signal width of 20 km s^-1 our sensitivity should be appropriate for detecting HI gas at column densities as low as 10¹⁹ cm^-2. To determine the disk orientation, we compare the tilted ring models directly against the full spectral database considering channels 12.88 km s^-1 wide (rather than fitting the moments of the flux distribution). The cube consists of 2475 positions (i.e., pixels 40 × 40 arcsec² wide) in which 21 cm emission has been detected, and for each position we have 25 velocity channels covering from –20 to –342 km s^-1 heliocentric velocities. We summarize below the main features of our method.

We use 110 tilted concentric rings in circular rotation around the center to represent the overall distribution of HI. Each ring is characterized by its radius R and by seven additional parameters: the H I surface density S_HI, the circular velocity V_c, the inclination i and the position angle PA with respect to the line of sight, the systemic velocity V_sys and the position shifts of the orbital centers with respect to the galaxy center (Δx_c,Δy_c). These last three parameters allow the rings to be nonconcentric and to have velocity shifts with respect to the systemic due to local perturbations (such as gas outflowing or infalling into the ring or M 31 tidal pull). Of these large set of rings, we allow only the parameters of eleven equally spaced rings, called the “free” rings, to vary independently. We set the properties of the innermost free ring to be the same as those of the next adiacent free ring (because of its small size it turns out to be highly unconstrained) and we keep the free rings to be the 11,22,33 ... 110th ring. The properties of rings between each of the free ring radii were then linearly interpolated. Each of the seven parameters of the ten free rings were allowed to vary. We assume that the emission is characterized by a Gaussian line of width w, which is an additional free parameter of the model, centered at V_c. We compute the 21 cm emission along each ring as viewed from our line of sight, and the synthetic spectrum at each pixel by convolving the emission from various ring pieces with the beam pattern. We then test how well the synthetic and observed spectra match by comparing the flux densities in 25 velocity channels. With this method, we naturally account for the possibility that the line flux in a pixel might result from the superposition along the line of sight of emission from different rings. As an initial guess for the free parameters of the tilted ring model we follow the results of Corbelli & Schneider (1997).

The assignment of a measure of goodness of the fit is done following two methods: the “shape” and the “v-mean” method. In the shape method, we minimize a χ² given by the sum of two terms: the flux and the shape term. The flux term is set by the difference between the observed and modeled fluxes in each pixel. The shape term retains information about the line shape only, which are lost when just the first few moments are examined. This is essential in the regions of M 33 where the emission is non-Gaussian, for example when the velocity distribution of the gas is bimodal (this is indeed the case for some regions in the outer disk of M 33, see CS). The shape term is given by the difference between the observed and the normalized modeled fluxes in each pixel and for every spectral channel. The normalized model spectrum is the flux predicted by the model in a given channel multiplied by the ratio of the observed to model integrated emission. In doing so the shape term is no longer dependent on the total flux. The shape term is computed only for pixels with flux larger than 0.2 Jy km s^-1/beam i.e. N_HI ≥ 1.3 × 10¹⁹ cm^-2. The error for the shape term is the rms in the baselined spectra, σ_n,i, which is the experimental uncertainty in the flux in the i-channel at the n-pixel. The noise in each channel of the datacube, σ_ch, is uniform and estimated to be 0.0117 Jy/beam/channel. As suggested by Corbelli & Schneider (1997) we estimate σ_n,i as: $\begin{array}{ccc}{\mathit{\sigma }}_{\mathit{n,i}}\mathrm{=}\mathrm{2}\frac{{\mathit{w}}_{\mathit{i}}\mathrm{\times }{\mathit{\sigma }}_{\mathrm{ch}}}{\sqrt{{\mathit{w}}_{\mathit{i}}\mathit{/}\mathrm{2.57}}}& & \end{array}$(A.1)where w_i = 12.88 km s^-1 is the width of the channel used to compare the data with the model and 2.57 km s^-1 is the database channel width. The flux term is affected by the uncertainty in the integrated flux and by calibration uncertainties, proportional to the flux. The calibration error forces the minimization to be sensitive to weak-line regions and is 5%. The resultant reduced χ² formula is the following: $\begin{array}{ccc}{\mathit{\chi }}_{\mathrm{shape}}^{\mathrm{2}}& \mathrm{=}& \frac{\mathrm{1}}{\mathit{N}\mathrm{-}{\mathit{N}}_{\mathrm{f}}}\sum _{\mathit{n}\mathrm{=}\mathrm{1}}^{\mathit{N}}\{\frac{\mathrm{(}{\mathit{I}}_{\mathit{n}}^{\mathrm{mod}}\mathrm{-}{\mathit{I}}_{\mathit{n}}^{\mathrm{obs}}{\mathrm{)}}^{\mathrm{2}}}{\mathrm{25}{\mathit{\sigma }}_{\mathit{n,i}}^{\mathrm{2}}\mathrm{+}\mathrm{(}\mathrm{0.05}{\mathit{I}}_{\mathit{n}}^{\mathrm{obs}}{\mathrm{)}}^{\mathrm{2}}}\\ & & \mathrm{+}\frac{\mathrm{1}}{\mathrm{25}}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathrm{25}}\frac{\mathrm{\left[}\mathrm{\right(}{\mathit{I}}_{\mathit{n}}^{\mathrm{obs}}\mathit{/}{\mathit{I}}_{\mathit{n}}^{\mathrm{mod}}\mathrm{)}{\mathit{I}}_{\mathit{n,i}}^{\mathrm{mod}}\mathrm{-}{\mathit{I}}_{\mathit{n,i}}^{\mathrm{obs}}{\mathrm{]}}^{\mathrm{2}}}{{\mathit{\sigma }}_{\mathit{n,i}}^{\mathrm{2}}}\begin{array}{c}\mathrm{⎫}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎬}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎭}\end{array}\end{array}$(A.2)where N is the number of pixels (2475) and N_f is the number of free parameters (71 for our basic model), leaving a large number of degrees of freedom. Given the difficulty in finding a unique minimum, we use a two-step method to converge toward the minimal solution. Since some of the parameters might be correlated, we begin by searching for minima over a grid of the parameters surrounding our first guess. We carried out several optimization attempts under a variety of initial conditions and with different orderings for adjusting the parameters. After iterating to smaller ranges of variation, we choose the parameter values that yield the minimum χ². In the second step, we check the minimal solution by applying a technique of partial minima. We evaluate the χ² by varying each parameter separately. We checked our solution by surrounding the galaxy with zero-flux observations for stabilizing the outermost ring. It is important to notice that the flux and the shape term yield a similar contribution to the minimum χ² value.

In the second method we determine a solution using the deviations of the integrated flux and of the intensity-weighted mean velocity along the line of sight at each pixel. We carried this out with another two-step procedure, allowing all 71 parameters to be varied. We start by keeping the ring centers and systemic velocity fixed; then the rings centers and velocities are considered as free parameters in the minimization as well. In order to keep the model sensitive to variations of parameters of the outermost rings, each pixel in the map is assigned equal weight. Pixels with higher or lower 21 cm surface brightness contribute equally to determine the global goodness of the model fit. Since the original data has a velocity resolution of 1.25 km s^-1, we arbitrarily set σ_m, the uncertainty in the mean velocity, equal to the width of about five channels (~6 km s^-1). This is simply a scaling factor that gives similar weight to the two terms in the χ² formula. The equation below defines the reduced χ² of the v-mean method: $\begin{array}{ccc}{\mathit{\chi }}_{\mathrm{mean}}^{\mathrm{2}}& \mathrm{=}& \frac{\mathrm{1}}{\mathit{N}\mathrm{-}{\mathit{N}}_{\mathrm{f}}}\sum _{\mathit{n}\mathrm{=}\mathrm{1}}^{\mathit{N}}\{\frac{\mathrm{(}{\mathit{I}}_{\mathit{n}}^{\mathrm{mod}}\mathrm{-}{\mathit{I}}_{\mathit{n}}^{\mathrm{obs}}{\mathrm{)}}^{\mathrm{2}}}{\mathrm{25}{\mathit{\sigma }}_{\mathit{n,i}}^{\mathrm{2}}\mathrm{+}\mathrm{(}\mathrm{0.05}{\mathit{I}}_{\mathit{n}}^{\mathrm{obs}}{\mathrm{)}}^{\mathrm{2}}}\\ & & \mathrm{+}\mathrm{(}{\mathit{V}}_{\mathit{n}}^{\mathrm{mod}}\mathrm{-}{\mathit{V}}_{\mathit{n}}^{\mathrm{obs}}{\mathrm{)}}^{\mathrm{2}}\mathit{/}{\mathit{\sigma }}_{\mathit{m}}^{\mathrm{2}}\}\end{array}$(A.3)where V_mod is the mean velocity predicted by the tilted ring model at the pixel n.

Given the large number of degrees of freedom, the increase of χ² corresponding to 1σ probability interval for Poisson statistics would be very small. The χ² standard deviation is of order 0.03 for the flux term and of order 0.006 for the shape term. Since the presence of local perturbations does not allow the model χ² to approach unity, we consider χ² fractional variations corresponding to the mean value of the two terms (i.e., 2%). By testing the χ² variations for each variable independently, we should have an indication on which ring and parameter is well constrained by the fitting procedure. Hence, we first arbitrarily collect all possible sets of tilted ring models which give local minima in the χ² distribution with values within 2% of the lowest χ² (which is 7.3 and 6.8 for the v-mean and shape method respectively). In Fig. A.1, we show V_c, Δx_c, Δy_c and ΔV_sys corresponding to an assortment of models whose χ² is within 2% of the absolute minimum value found. The adopted systemic velocity is V_sys = −179.2 km s^-1. The displacements of the ring centers and systemic velocities are not very large and increase going radially outward, as does the scatter between solutions corresponding to partial minima. The value of the velocity dispersion we find from the minimization is of order 10 km s^-1.

Fig. A.1 Parameters Vc, Δxc, Δyc, and ΔVsys corresponding to an assortment of tilted ring models whose χ2 is within 2% of the absolute minimum value found. In the panels to the left the parameters for the v-mean minimization are shown, and in the panels to the right those for the shape minimization. The continuous line connects the parameters of the best fit tilted ring models used for deriving the rotation curve.

For each minimization method we then select a tilted ring model between those with acceptable χ² using the maximum

north-south symmetry criterion for rotation curves relative to the two galaxy halves. The corresponding values of i and PA are shown in Fig. 3 of Sect. 4 with the relative uncertainties. The uncertainties are computed by varying one parameter of each free ring at a time, around the minimal solution until the χ² increases by 2%. Simultaneous parameter variations within the given uncertainties gives χ² variations larger than 2%. In deriving the rotation curve, we take into account the uncertainties considering deconvolution models in which the inclination or PA of all the rings vary simultaneously. In this case we consider only 35% of the uncertainties displayed in Fig. 3, in either PA or i, in order to have a χ² within 2% of the minimal solution.

Appendix B: Finite disk thickness corrections

The “rotation curve” is the azimuthal component of the velocity in the equatorial plane of the disk at a given galactocentric distance. However, the 21 cm spectrum observed at a certain position in an inclined disk depends not only on the azimuthal velocity in the plane but also on two additional effects: the smearing due to the extent of the telescope beam and the vertical extension of the disk; both become more severe with increasing inclination. The tilted ring model fit runs over a smooothed database, spatially and in frequency, whose final geometrical parameters are then used to derive the rotation curve from a higher resolution dataset. Therefore, instead of including disk finite thickness effects in the tilted ring models we prefer to account for this and for the beam smearing in the 21 cm spectral cube at the original resolution using a a set of numerical simulations.

We assume the gas to be in an azimuthally symmetric disk inclined with respect to the line of sight according to the tilted ring model fit. The gas radial distribution is set equal to that given by the integrated spectral profile, while the vertical profile is modeled by an exponential with a folding length of 0.3 kpc. Only a disk component is considered with no allowance for a halo component (Oosterloo et al. 2007). For the beam we used a Gaussian with FWHM = 20 arcsec truncated at a radius of 24 arcsec; the channel spacing in the spectrum is 1.25 km s^-1. The input rotation curve is that obtained by the observed peak and mean velocities, and the random velocity is assumed to be isotropic and spatially constant with FWHM = 8.0 km s^-1, as observed in most of the disk. We ran simulations with and without a vertical rotation lag according to Oosterloo et al. (2007).

Using the above parameters, we simulate a synthetic HI data cube. The cube is used to derive the corresponding velocity for each position either the peak or the flux-averaged mean, and then build a simulated rotation curve. The differences between the simulated and the input rotation give the average corrections to the rotation curve as function of radius. We shall refer to these corrections as finite disk thickness corrections. They are of order 2–3 km s^-1 and reach values of 5 km s^-1 only within the innermost 200 pc. We apply these corrections to the rotational velocities.

Acknowledgments

We are grateful to Stéphane Charlot and Gustavo Bruzual for kindly providing us with the latest revision of the BC03 stellar population synthesis models ahead of publication. We aknowledge financial support from PRIN MIUR 2010-2011, project “The Chemical and Dynamical Evolution of the Milky Way and Local Group Galaxies”, prot. 2010LY5N2T.

References

All Tables

All Figures

	Fig. 1 Upper panels: integrated 21 cm line intensity and the intensity weighted mean velocity over the whole galaxy using data at 130 arcsec resolution. Lower panels: peak intensity and the velocity at peak maps using VLA+GBT data at a resolution of 20 arcsec.
In the text

Fig. 2 Panel a): color composite RGB image of the B − V − I mosaics from the Massey et al. (2006) survey. Panel b): stellar mass map derived from the median likelihood estimator (described in Sect. 3.1) based on LGS BVI bands. Panel c): stellar mass map obtained from the Bayesian median likelihood estimator, based on the full LGC and SDSS BVIgi dataset. Panel d): apparent stellar M/L derived from median likelihood stellar mass map based on LGS+SDSS BVIgi bands and the observed LGS V-band surface brightness.

In the text

Fig. 3 Inclination and position angle of the 11 free tilted rings used for deconvolving the 21 cm data. The open triangles (in black in the online version), connected by a dashed line, shows the value using the shape minimization method, the open squares (in red in the online version), connected by a continuous line, are relative to the v-mean minimization method (see Appendix A for details). The uncertainties are computed by varying one parameter of one free ring at a time around the minimal solution until the χ2 of the tilted ring model increases by 2%.

In the text

Fig. 4 Bottom panel: rotational velocities for each pixel in the peak map (small dots, red dots for the southern half and blue dots for the northern half in the online version) and average values in each radial bin (large symbols). The model-shape has been used to deconvolve the data. Triangles (red in the online version) are for the southern half while squares (blue in the online version) are for the northern half. Top panel: average rotation curve for peak velocities. Bins are 2 arcmin wide in the radial direction.

In the text

Fig. 5 Bottom panel: rotational velocities for each pixel in the first moment-1 map (small dots), the average values in each ring (large symbols, red symbols are for the southern half and blue symbols for the northern half in the online version). The model-mean has been used to deconvolve the data. Triangles are for the southern half while squares are for the northern half. The average rotation curve using moment-1 (mean) velocities is shown in the top panel. Bins are 2 arcmin wide in the radial direction.

In the text

Fig. 6 Inner rotation curve as traced by the peak velocities at high resolution (filled squares in blue in the online version), and at low resolution (filled dots in red in the online version). The CO data are shown with star symbols (in green in the online version). The open square symbols (in red in the online version) are for rotational velocities as traced by the moment-1 map. For each curve the weighted mean velocities between that relative to deconvolution model-shape and model-mean are used.

In the text

	Fig. 7 Tilted ring position angles and inclinations (open square symbols and filled dot symbols respectively for the ROTCUR and KINEMETRY fitting routines).
In the text

Fig. 8 Some of the tilted ring parameters fitted by the routines ROTCUR and KINEMETRY. The top panel shows the term k5, which is a higher order term in the KINEMETRY outputs. In the bottom panel we show the shifts of the systemic velocity, and, in the middle panel, the orbital center shifts along the y-axis in the south-north direction (filled square symbols) and along the x-axis oriented from east-to-west (open square symbols).

In the text

Fig. 9 Adopted rotation curve (open black circles) for the inner regions (bottom panel) using the peak velocities and for the outer regions (top panel) using the moment-1. The dispersions take into account variation of Vr between different deconvolution models, the 2% variations around the minimum χ2, and the standard deviation of the mean velocity in each bin. In both panels, the open square and asterisk symbols (in magenta and green in the online version) show Vr according to the ROTCUR and KINEMETRY deconvolution model, respectively. In the top panel the open triangles are for the moment-1 velocities deconvolved according to model-shape and model-mean (in red and blue colors, respectively in the online version).

In the text

Fig. 10 Neutral atomic gas surface density perpendicular to the galactic plane of M 33 (small filled triangles) and the function that fits the data (continuous line,red in the online version) after the 21 cm line intensity has been deconvolved according to tilted ring model-shape. Asterisk symbols indicate the stellar mass surface density using the BVIgi stellar surface density map. The dashed line is the fit to the stellar surface density and the extrapolation to larger radii. Open squares(in blue in the online version) show for comparison the surface density using the BVI mass map. The heavy weighted line is the total baryonic surface density, the sum of atomic and molecular hydrogen, helium and stellar mass surface density.

In the text

Fig. 11 Mass-to-light ratio as a function of galactocentric radius according to the BVIgi stellar surface density map (filled triangles, in black in the online version). This is computed by averaging the stellar mass map and the V-band surface brightness along ellipses corresponding to the tilted ring model. We then compute M∗/LV as the ratio between these two quantities in units of M⊙/L⊙. We show for comparison the same ratio relative to the BVI stellar surface density map (open triangles, in red in the online version). The error bars are the standard deviations relative to azimuthal averages and do not take into account uncertainties in the mass determination.

In the text

Fig. 12 Rotation curve of M 33 (filled triangles) and best-fitting dynamical models (solid line) for NFW dark matter halo profiles (dot-dashed line). The small and large dashed lines (red and blue in the online version) show the gas and stellar disk contributions to the rotation curve, respectively. In the bottom panel, the highest stellar contribution and lowest dark halo curve are for the best-fitting dynamical model using the BVIgi mass map, the other curves are for the BVI mass map. The top panel shows the most likely dynamical models for the two mass maps when the likelihood of the dynamical fit is combined with that of the stellar mass surface density and of the C − Mh relation resulting from simulations of structure formation in a ΛCDM cosmology.

In the text

Fig. 13 Confidence areas in the log C − log Mh/h-1 plane and in the M∗ − Mh/h-1 plane for the dynamical fit to the rotation curve using a NFW dark matter halo profile: 68.3% (darker regions, blue in the online version), 95.4% (magenta) and 99.7% (lighter regions, cyan). The mass distribution is that given by the BVIgi colors for the areas shown in a) and b), and by the BVI colors for the areas in c) and d). In panels a) and c) we show the C–Mh relation (solid line) and its dispersion (dashed line) for ΛCDM cosmology. In panels b) and d) the dash-dot line indicates the most likely stellar mass according the synthesis models, which rely on the colors shown next to panel labels.

In the text

Fig. 14 68.3% (darker regions, blue in the online version), 95.4% (magenta) and 99.7% (lighter regions, cyan) confidence areas in the log C − log Mh/h-1 plane and in the M∗ − Mh/h-1 plane when the probability of the dynamical fit is combined with the probability of the stellar mass surface density and of the C − Mh distribution from simulations of structure formation in a ΛCDM cosmology. The mass distribution is that given by the BVIgi map for the areas shown in a) and b), and by the BVI map for the areas in c) and d).

In the text

Fig. 15 Rotation curve of M 33 (filled triangles) and the best-fitting model (solid line) for a dark matter halo with a constant density core (bottom panel) and according to MOND (top panel). The small and large dashed lines (red and blue in the online version) show the gas and stellar disk Newtonian contributions to the rotation curves, respectively; the dark halo contribution in the bottom panel is shown with a long dashed line.

In the text

Fig. 16 The 95.4% (magenta in the online version) and 99.7% (lighter regions, cyan) confidence areas in the log r0 − log ρ0 plane and in the log r0 − log M∗ plane relative to the rotation curve fit using a core model and the BVI mass map. In the top panel, the dash-dot line indicates the most likely stellar mass according the synthesis models, which rely on the BVI colors. There are no core models in the 68.3% confidence area for the BVI mass map and no acceptable core models for the BVIgi mass map down to the 99.7% level. The solid line in the bottom panel indicates the inverse log ρ0 − log r0 correlation found by Donato et al. (2009).

In the text

Fig. 17 A naive view of M 33: inner stellar surface density (from red to magenta colors in the online version) superimposed on the HI gas distribution in log-scale (light pink) are shown together with the stellar+gas contribution to the rotation curve (dashed cyan line) and with observed rotation curve (filled dots, in blue). The analysis presented in this paper has pinned down the average dark matter density distribution in the halo (light gray) which explains the observed rotation curve. The figure extrapolates the halo beyond the region analyzed in this paper in a naive smooth fashion. The total mass model fit to the rotational data, through the region analyzed in this paper, is shown by the continuous blue line.

In the text

Fig. A.1 Parameters Vc, Δxc, Δyc, and ΔVsys corresponding to an assortment of tilted ring models whose χ2 is within 2% of the absolute minimum value found. In the panels to the left the parameters for the v-mean minimization are shown, and in the panels to the right those for the shape minimization. The continuous line connects the parameters of the best fit tilted ring models used for deriving the rotation curve.

In the text