The degeneracy of M33 mass modelling and its physical implications

The Local Group galaxy M33 exhibits a regular spiral structure and is close enough to permit high resolution analysis of its kinematics, making it an ideal candidate for rotation curve studies of its inner regions. Previous studies have claimed the galaxy has a dark matter halo with an NFW profile, based on statistical comparisons with a small number of other profiles. We apply a Bayesian method from our previous paper to place the dark matter density profile in the context of a continuous, and more general, parameter space. For a wide range of initial assumptions we find that models with inner log slope $\gamma_{\rm in}<0.9$ are strongly excluded by the kinematics of the galaxy unless the mass-to-light ratio of the stellar components in the $3.6\mu$m band satisfies $\Upsilon_{3.6}\geq2$. Such a high $\Upsilon_{3.6}$ is inconsistent with current modelling of the stellar population of M33. This suggests that M33 is a galaxy whose dark matter halo has not been significantly modified by feedback. We discuss possible explanations of this result, including ram pressure stripping during earlier interactions with M31.


INTRODUCTION
Cosmological models of the formation of dark matter haloes predict cusped density profiles (Dubinski & Carlberg 1991;Navarro et al. 1996), which do not appear to match the dark matter density profiles inferred from observations of rotation curves of disk galaxies (Gentile et al. 2004).
Observations of galaxies from THINGS, The HI Nearby Galaxy Survey ) have provided good observational constraints on the rotation curves (and thus density profiles) of nearby galaxies, as explored in de Blok et al. (2008) and Hague & Wilkinson (2013) (hereafter HW13). These improved velocity data have allowed more precise constraints on halo density profiles than were possible in previous papers that addressed the cusp-core problem.
In this paper we apply the Bayesian method presented in HW13 to the galaxy M33, which has been previously claimed to have a density profile described by a single power law ρ ∝ r −1.3 , compatible with the NFW profile (Corbelli & Salucci 2000). We also examine the more recent claim of Seigar (2011) that the NFW profile itself best represents the dark matter halo of M33. This result is at odds with the observational claim (e.g. de Blok et al. 2008) that smaller galaxies are best described with cored halo profiles. If confirmed, it might have implications for the role of baryonic feedback on the evolution of halo profiles.
In Section 2 we describe how we reproduce the baryonic mass modelling and rotation curve of M33, and the MCMC technique we use on this model. In Section 3 we peter.hague@le.ac.uk miw6@le.ac.uk analyse the output of the MCMC chains and in Section 4 we discuss our result in the context of previous papers and the current paradigm of galaxy formation.

MODELLING OF M33
We decompose the rotation curves into four components: a stellar disk, a bulge, a gas disk and a dark matter halo. The circular velocity contribution of each component is added in quadrature to produce a proposed rotation curve, to be compared with observations.

Baryonic Mass Models
The Seigar (2011) model consists of a gas component taken from Corbelli & Salucci (2000), along with two stellar components -a bulge and a disk. These components are distinguished photometrically, rather than by velocity structure. The disk is assumed to be exponential, whilst the bulge is taken to be a Sérsic profile (although in this case the best fit was found to be n = 1, making it into an exponential profile.) Table 1 shows the values used to generate the two components. Values for the bulge appear to differ from those quoted in Seigar (2011) as we have converted them to the equivalents for an exponential disk. We explore models with both fixed and freely varying baryonic mass-to-light ratios. We do not include a mass-to-light gradient in the disk as the estimated gradient in Seigar (2011) would have negligible impact compared to the rotation curve errors, but would introduce additional complications into the parameter space in the case where we use a free baryonic mass-to-light ratio (calculated in the Spitzer 3.6µm band) Υ 3.6 .
Our first model for the stellar mass allows the massto-light ratio of the stellar disk to vary freely over a 1.14 × 10 9 Gas disk scale height zgas (kpc) 0.5 Gas disk mass Mgas (M ) 3 × 10 9 Note. -The two stellar components are modelled as exponential disks using parameters from Seigar (2011), and the gas component is modelled using radially binned surface density data provided by Corbelli & Salucci (2000). large range. We consider a second, fixed stellar model using the mass modelling of Oh et al. (2008), along with the J − K values for M33 taken from the 2MASS Large Galaxy Atlas (Jarrett et al. 2003). These values are measured within a 20 mag arcsec −2 isophote, which in M33 corresponds to a radius r ≈ 6kpc. For M33, [J − K] = 0.891 which gives a mass-to-light ratio Υ 3.6 = 0.67. This is more consistent with current estimates of mass-to-light ratios of similar nearby galaxies (e.g. Meidt et al. 2014) than the value from Seigar (2011).
We have been provided with the radial surface density of neutral atomic gas and molecular gas used in Corbelli & Salucci (2000) by the authors. We processed this using the ROTMOD task in GIPSY 1 , which employs the method described in Casertano (1983) to generate a rotation curve contribution. Using a sech 2 vertical scaling law, and the value for z gas given in Corbelli & Salucci (2000), we are able to produce a contribution consistent with that shown in Figure 6 of that paper. The gas contribution used in Seigar (2011) is identical to this.

Dark Matter Models
To model the dark matter halo we used an α − β − γ profile, but as in HW13 transformed to remove the degeneracy between ρ s and r s ρ(r) = Σ max G v 2 max ( r rs ) γ (1 + ( r rs ) 1/α ) α(β−γ) (1) where r s is the scale radius, v max is the peak velocity of the dark matter rotation curve, α, β and γ are shaping parameters, and Σ max is the normalised surface density at r max ≡ r(v max ), given by Contrary to the statements in Adams et al. (2014), it does not matter that the parameters of this halo profile are still degenerate to some extent as we focus on physical properties of each halo model, within the data range, rather than parameters such as γ. Our approach has the advantage that we do not impose as strong a prior link between the log slope of the halo at small and large radii (due to the larger parameter space) and we do not extrapolate beyond the data range. We have confirmed the utility of this approach through extensive testing in HW13.
1 http://www.astro.rug.nl/∼gipsy/ Note. -The number of radial bins in models in the second set of models is reduced by removing the feature at the outermost part of the rotation curve in Corbelli & Salucci (2000). The priors for Υ 3.6 are (1) a freely varying Υ 3.6 in the range [0.1, 5] (models A, D); (2) the value from Seigar et al. (2008), using a central mass-tolight ratio in the Spitzer 3.6µm band of Υ 3.6 = 1.25 ± 0.10 taken from Seigar et al. (2008) (model B); (3) a value for Υ 3.6 derived more recently with assumptions from Oh et al. (2008), as described in the text (model C).

MCMC analysis
Following the method in HW13 we ran 8 parallel MCMC chains, with a total of ∼ 4 × 10 7 models, with each model being evaluated on the basis of χ 2 , and model selection using the Metropolis-Hastings algorithm (Hastings 1970) implemented through CosmoMC (Lewis & Bridle 2002).
We present the results for eight runs, shown in Table  2. For the A runs we allowed a free mass-to-light ratio, Υ 3.6 , with a range [0.1, 5] to generously cover possible stellar contributions from no disk contribution through to a super-maximal disk. For the B runs we calculated the baryonic components as in Seigar (2011); for the C runs we use the mass-to-light ratio calculated above in Section 2.1, and the D runs use two independent values for Υ 3.6 for the disk and bulge components, using the same ranges as the A runs.
The parameter space is mirrored around γ = 0, using a range [-2,2] for model selection but taking the absolute value for likelihood testing, so that potentially viable cored profiles are not located at the boundary of the parameter space (see HW13 for details).
We have discarded 10,000 models from the beginning of each chain to allow for burn-in, which we find is sufficient for all the chains to move to areas of high likelihood (χ 2 red < 2.5). However, after this point some chains explore secondary peaks before finding the main peak. These secondary peaks are a genuine part of the distribution, as can be verified by the fact that chains sometimes leave the main peak to explore them for an extended period. Our choice to include these early models (after the first 10,000) gives us a more complete picture of the multi-modal probability distribution. However we should note that relative likelihood comparisons between the peaks cannot be done, because exploration of secondary peaks is influenced by the starting point of the chain.

RESULTS
A typical rotation curve is shown in Figure 1. This is taken from the peak of the distribution produced by the A1 chain and shows a cusped halo.
For A1, we measured the log slope of the dark matter profile at the innermost radial bin, γ in , and found it to be degenerate with the mass-to-light ratio as shown in The dark blue curve is the gas (atomic and molecular) contribution, the orange curve is the central bulge-like stellar component, the red curve is the stellar disk, the green line is a proposed dark matter halo (taken from the most occupied bin from the parameter space of the A1 run) and the light blue line is the expected rotation curve. Observed data are in black. The mass-to-light ratio of the stellar components has been found by fitting the rotation curve rather than by modelling the stellar population in this case. Figure 2. The distribution is binned on a 128 × 128 grid and then contours placed that enclose 68%, 95%, and 99% of all models. Our analysis favours models with steep inner cusps and high mass-to-light ratios, with a tail in the distribution moving towards flat haloes with even higher values of Υ 3.6 . The two areas are connected by a bridge of models which is not shown here as the density is below the 3σ level. The exact combination of stellar disk mass and dark matter halo favoured in Seigar (2011) (Υ 3.6 = 1.25 and γ in 1) is disfavoured at over 3σ here.
To check that all 8 chains are converged on the same distribution, we calculated the ratio of the variance of the means of the chains σ(x), to the mean of the variances of the chainsσ(x), for each parameter. For γ this gave a value of 147.7 as the chains did not converge to the same distribution in this parameter, but after that the highest value was 0.1 for β, indicating the rest of the parameters are well converged. Inspection of the distribution for each chain showed that in some cases the chain had found only the main peak shown in Figure 2 whereas in some cases the chain only found the left hand side of the tail. As some chains managed to integrate the entire distribution as shown in Figure 2, it is reasonable to assume that given enough time any chain will converge on the same result. Combining the distributions may produce an incorrect relative model density between the density of the main peak and the tail, but our analysis does not make use of this.
At the low γ in end of the plot is the maximal disk case, where the baryonic component of the galaxy contributes almost all of the rotation curve. Note that despite being -Contour plot of the mass-to-light Υ 3.6 versus inner log slope γ in . The red contour contains 68% of all the models in the MCMC chains, the green contour contains 95% and the blue contour contains 99%.
flat, this region is not flush with the upper boundary of the Υ 3.6 range -higher values are excluded by the rotation curve.
The rotation curve presented in Corbelli & Salucci (2000) shows an apparent feature in the outermost part of the rotation curve where the rotation velocity begins to increase, after having levelled off (see Figure 1). This feature may be modelled by a flat dark matter density profile extending throughout the radial range, coupled with the maximal baryonic component to model the shape of the rotation curve at low r. We confirm this by calculating a measure of the core radius, r 1 , which is the radius at which the log slope of the dark matter halo reaches -1. Figure 3 shows that this value is high (on the order of the radial extent of the data) for high values of Υ 3.6 .
If these outermost data points represented a genuine feature of the density profile and the disk were not maximal, it would require a sharp increase in the dark matter density at this point as the baryonic component is marginal here. As there is no obvious mechanism to form such a shell of dark matter, we cannot take this model to be correct at high r.
We consider an artefact of the tilted ring method used to generate this rotation curve to be more a likely explanation of this rotation curve feature. In Corbelli & Salucci (2000), an initial set of radial bins in this rotation curve are generated by fitting a parameterised ring to the HI velocity field of the galaxy, under the assumption of entirely circular motion, and then additional radial bins are calculated from interpolating between the neighbouring rings. If the assumptions of the model do not hold e.g. in the presence of significant radial motion, then the parameters of a particular ring may be invalid. An overestimation of inclination would lead to an overestimation in rotation velocity, and due to the interpolation, a single incorrect inclination can account for the apparent feature q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q seen in the last two radial bins. Considered without the final two radial bins, it is not clear the feature exists at all. Figure 4 shows that for run B1 (where Υ 3.6 is taken from Seigar (2011)), the distribution of γ in favours a cusped density profile, steeper than the best fitting NFW profile which we take to be representative of cusped profiles in general for this purpose. The distribution is bimodal, but removing the last two data points removes the second peak. This peak is then purely dependent on a feature which may be an artefact.
Runs D1 and D2, with an additional free parameter for the Υ 3.6 value of the inner stellar component, produced a similar degeneracy to the A runs. Figure 5 shows the relation between γ in and both values of Υ 3.6,bulge . There is a weak, secondary peak of models featuring a maximal bulge and a flat (γ in < 0.5) inner slope, but the main peak in Υ 3.6,disk versus γ in is unaffected.

Comparison with Previous Work
We have found a result more in agreement with the value of γ in = 1.3 implied by Corbelli & Salucci (2000) than the assertion by Seigar (2011) that M33 is best described by an NFW halo. Whilst the NFW halo does fit the rotation curve decomposition with χ 2 red = 1.18, there are many other haloes (mostly steeper) that also fit the same data equally well.
It is important to note that, as discussed in HW13, with multiple parameter models, reduced chi-squared (χ 2 red ) is not reliable across the entire parameter space because it is calculated assuming that the degrees of freedom are constant across the parameter space, which is not always the case. For instance, if a model includes a high stellar mass, and a reduced contribution of the dark matter halo to the rotation curve at small r, then any shape parameters of the halo are going to become less relevant to the quality of the fit. This was described, in an extreme case, in HW13 for the case of constructed high surface brightness rotation curves where the dark matter contribution to the rotation speed was smaller than the error bars. Even in less extreme cases, the parameter β is often not fully utilised, if the scale radius of the halo is large enough that it does not become the dominant shape parameter within the data range.
The probability of adding a model, a point in param- eter space, to one of our MCMC chains is not based on the absolute value of its χ 2 but on the gradient of the goodness of fit, i.e. the relative goodness of fit compared to some other point the chain may arrive from. The final result is essentially an integral of this value over all possible starting points -but with a substantial weighting towards nearby points due to the Gaussian shape of the model selection function in the Metropolis-Hasting algorithm. This means that the MCMC result does not rely on the goodness of fit being a globally correct representation of likelihood. Given that in a tilted ring model, errors are computed from azimuthal variations in the inferred circular velocity, it is not immediately clear they are strictly robust, and an analysis (such as ours) that does not rely on their robustness is preferable.
In HW13 we showed that the MCMC method was able to recover the correct halo model from synthetic data with artificially inflated error bars more tightly than would be naively expected. We attribute this not only to the method being less dependent on the actual values of the goodness of fit, but also to the physically reasonable prior assumption of a smooth dark matter density profile being able to guide the fit more than would be possible with the errors alone.
In Seigar (2011) it is claimed that fitting of the NFW profile to rotation curve bins outside r = 7kpc is evidence that the this profile "best represents these data". Only the prior assumption of an NFW profile makes the fit to the outer rotation curve relevant to the determination of the inner density profile. Without assuming a strong link between inner and outer data points as a prior, a steeper inner density profile is favoured when using the Υ 3.6 value used in that paper.
The second claim of Seigar (2011) is that the NFW concentration parameter c vir is related to the spiral arm pitch angle P . Taking c vir purely as a density profile parameter without any implications for inner profile slope, correlations between it and other physical parameters are not necessarily in conflict with our conclusions.

Mass-to-light Degeneracy
It is clear from our result that there is a substantial degeneracy between the baryonic mass-to-light ratio Υ 3.6 and the inner log slope γ in . The scale length of the stellar disk is approximately equal to the radius at which the dark matter halo transitions from rising to flat, so the most obvious cause of the degeneracy is the degree to which the rising part of the stellar contribution to the circular velocity is used to model the rising part of the overall rotation curve. To investigate this, we assume that the degeneracy of the fit near the peak of the disk rotation curve approximates that of the entire rotation curve, and that the inner dark matter halo can be taken to be a single power law (which in our results it can; over 98% of models in A1 have a scale radius outside the data range). We can then estimate the degeneracy using v 2 max,star + v DM (R max ) 2 = constant (3) where v max,star is the maximum velocity contribution of the stellar component and v DM (R max ) is the velocity contribution of the dark matter halo at that radius. Given that v 2 max,star ∝ Υ 3.6 and v DM (R max ) 2 ∝ (R max /r s ) 2−γin (as γ in is assumed to be representative of the log slope at this radius, due to our assumption of a single power law) we can write where a, b, and c are constants. We binned γ in values of a subset of models from A1 (34846 models chosen from all chains with a probability 10 −5 , ignoring the first 10,000 models for burn-in) to produce a set of 2000 bins with a uniform number of points. We fit the relationship to these data and find a = 2.6, b = 6.3, and c = 0.089 as shown in Figure 6. The baryonic mass-to-light ratio in the maximal disk case, Υ 3.6,max = a − bc 2 and in the no disk case, γ in,no = 2 − (log a − log b)/log c. q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq -Fit to the degeneracy between log slope at the inner most bin (γ in ) and mass-to-light ratio in the Spitzer 3.6µm band (Υ 3.6 ). Points are 34846 models chosen from the chain, after burnin, with a probability 10 −5 . The line shows the relationship Υ 3.6 = 2.6−6.3×0.089 2−γ in , fitted to 2000 bins, and should not be taken as valid outside the range of the models shown. Green and blue points are constraints on γ in found by models B1 and C1 respectively. See text for discussion.
In runs D1 and D2, we tested whether this degeneracy would be changed by having independent Υ 3.6 values for both stellar components. The scale radius of the inner component is too short to meaningfully contribute to the circular velocity at the point in the rotation curve which we have determined is the strongest driver of the degeneracy between γ in and Υ 3.6,disk . The distribution of Υ 3.6,bulge is unsurprisingly wide, due to its dependence on a small number of the innermost data points.

Alternate Mass Models
It is possible the apparently large stellar mass in our cored models could be accounted for by molecular gas. The data we obtained from Corbelli & Salucci (2000) included the molecular gas fraction, calculated from CO emission, along with the atomic gas component of the rotation curve. The factor used to calculate total molecular gas mass from CO emission, X CO , is estimated based on observations of the Milky Way and may not be correct for M33. In the mass modelling we use, the stellar disk has a mass 3.8 × 10 9 M and the molecular gas disk component has a mass 3 × 10 8 M , so the latter would represent an increase in Υ 3.6 of ∼ 0.05. To account for the difference between the modelled Υ 3.6 = 0.67 and the Υ 3.6 = 2.5 required for a flat halo would require X CO to be 37 times larger. In Dame et al. (2001) the 1σ relative error for this ratio was found to be less than 0.17 for nearby clouds in the Milky Way and thus a factor 37 increase seems unlikely.
The fact that the NFW profile is able to fit the rotation velocities in previous work does not itself convey the properties of the NFW profile (i.e. a log slope of -1 at r = 0) upon the galaxy, given that we have shown there are many profiles with differing properties that also provide good fits. Asserting a specific density profile imposes a relation between the log slopes at small and large radii as a prior, which means that the parameters that produce the highest likelihood might not correspond to a halo that fits all radii optimally. Thus, it cannot be concluded that a low reduced χ 2 value for an NFW fit gives a high posterior probability for specific analytic properties of the NFW profile, e.g. the central ρ ∝ r −1 cusp. Allowing more freedom in the profile, and fully exploring the parameter space with MCMC, resolves these issues and provides a more robust description of the dark matter density profile across the entire radial range of the rotation curve data. 4.4. Ram pressure stripping Bekki (2008) found that M33 encountered M31 with a periapsis of ∼ 100kpc at 4-8Gyr before present. Ram pressure stripping of M33 by an outflow from M31 could have had an impact on the progression of feedback in M33 (e.g. Nayakshin & Wilkinson 2013). The prompt removal of low density gas by stripping would cause the dark matter halo to relax into a shallow, less concentrated state, as modelled in HW14, so it might be reasonably assumed that such an event would lead to a cored dark matter halo, but this is not necessarily the case.
It is reasonable to assume the the process of contraction would mean that after initial baryon infall, but before feedback begins, the dark matter halo would have an inner log slope steeper than γ in = 1. Feedback models such as Read & Gilmore (2005) and Ogiya & Mori (2012) require multiple outflow and inflow events to account for the transition from such steep initial haloes to the presently observed flatter inner haloes. If this process were interrupted early on, it could prevent a sufficiently large amount of feedback that, even though the event itself would flatten the halo slightly, it would still retain an inner halo that is steep relative to those of similar galaxies.

CONCLUSION
We have modelled the rotation curve of M33 using the MCMC-based approach we presented in HW13. We have quantified and understood the degeneracy between baryonic mass-to-light ratio Υ 3.6 and, and the log slope at the inner bin γ in . We cannot resolve the conflict between observations of similar galaxies (Kuzio de Naray et al. (2006) and Kuzio de Naray et al. (2008)) and the MCMC analysis of the M33 rotation curve without assuming Υ 3.6 > 2, which is difficult to reconcile with stellar population modelling. We find that with a lower fixed Υ 3.6 = 0.67, an NFW halo is compatible with the data, but that this part of parameter space is not strongly favoured when we relax the constraint on Υ 3.6 . We strongly exclude the combination of Υ 3.6 < 2 and a halo profile inner log slope γ in < 0.9, for a comprehensive range of assumptions.
The relation we find between Υ 3.6 and γ in admits at least the three following scenarios: 1. there is a great deal more mass in the disk of M33 than is accounted for by standard modelling of stellar populations and molecular gas clouds.
2. the halo of M33 deviates significantly from spherical symmetry, being flattened at small disk radius and less so in the outer part of the galaxy.
3. the dark matter halo has a much steeper inner profile than would be expected from hydrodynamical simulations of galaxy formation (e.g. Governato et al. 2010;Maccio et al. 2011). This could occur if M33 were dominated by the process of contraction. The above simulations show both contraction of the halo steepening the inner profile, and feedback flattening it. In the absence of any obvious source of significant additional disk mass, and assuming no fundamental error in the view of baryon-dark matter interaction in galaxy formation, we propose that the history of the dark matter halo in M33 is dominated by contraction. Ram pressure stripping by M31 before feedback flattening the halo is a possible physical mechanism by which this could have happened.
We will investigate the last two possibilities by applying our modelling scheme to cosmological simulations in a future work.

ACKNOWLEDGEMENTS
This work is based [in part] on observations made with the Spitzer Space Telescope, which is operated by the Jet Propulsion Laboratory, California Institute of Technology under a contract with NASA. The CosmoMC code was written by Anthony Lewis (Lewis & Bridle 2002). We thank Edvige Corbelli for supplying the initial rotation curve decompositions. We acknowledge Wyn Evans for valuable discussions during a visit by MIW to the Aspen Center for Theoretical Physics, and Simon White for bringing work on M33 to our attention. MIW acknowledges the Royal Society for support through a University Research Fellowship. PRH acknowledges STFC for financial support. This research used the Complexity HPC cluster at Leicester which is part of the DiRAC2 national facility, jointly funded by STFC and the Large Facilities Capital Fund of BIS.