# Mass models of disk galaxies from the DiskMass Survey in MOND

## Abstract

This article explores the agreement between the predictions of Modified Newtonian Dynamics (MOND) and the rotation curves and stellar velocity dispersion profiles measured by the DiskMass Survey. A bulge-disk decomposition was made for each of the thirty published galaxies, and a MOND Poisson solver was used to simultaneously compute, from the baryonic mass distributions, model rotation curves and vertical velocity dispersion profiles, which were compared to the measured values. The two main free parameters, the stellar disk’s mass-to-light ratio () and its exponential scale-height (), were estimated by Markov Chain Monte Carlo modelling. The average best-fit K-band stellar mass-to-light ratio was . However, to match the DiskMass Survey data, the vertical scale-heights would have to be in the range to pc which is a factor of two lower than those derived from observations of edge-on galaxies with a similar scale-length. The reason is that modified gravity versions of MOND characteristically require a larger to fit the rotation curve in the absence of dark matter and therefore predict a stronger vertical gravitational field than Newtonian models. It was found that changing the MOND acceleration parameter, the shape of the velocity dispersion ellipsoid, the adopted vertical distribution of stars, as well as the galaxy inclination, within any realistic range, all had little impact on these results.

###### keywords:

cosmology: dark matter; galaxies: kinematics and dynamics; methods: numerical## 1 Introduction

Understanding the dynamics of disk galaxies is essential to the vetting process of theories of galaxy formation and cosmology (Flores & Primack 1994; de Blok & McGaugh 1998; de Blok et al. 2001; van den Bosch & Swaters 2001; Swaters et al. 2003; Gentile et al. 2004; Gilmore et al. 2007; de Blok 2010). Disk galaxies moderately inclined to the line of sight (50-80) can provide an HI rotation curve from which the dark matter (DM) content and distribution can be deduced in the context of Newtonian dynamics (Bosma 1978; Rubin et al. 1978; Bosma 1981a, b; Sofue & Rubin 2001). However, there exist degeneracies between the DM halo density profile, the stellar mass-to-light ratios of the stellar components and the scale-height of the stellar disk (van Albada et al. 1985; Kuijken & Gilmore 1989, 1991; Angus et al. 2012).

The vertical stellar distribution is assumed to be a declining exponential on both sides of the mid-plane. Thus, the scale-height of the stellar disk is simply the exponential rate at which the stellar luminosity density drops, with increasing height above or below the mid-plane of the disk. It is assumed to be constant with radius, which is supported by observations (e.g. van der Kruit & Searle 1981, 1982; Bizyaev & Mitronova 2002).

The stellar mass-to-light ratio () is the ratio between the mass of a stellar population and the luminosity, as observed through a particular bandpass. It cannot be observed directly, so it is ordinarily inferred by numerical stellar population synthesis models. This modelling crucially depends on the star-formation and chemical enrichment history considered and on the initial mass function (IMF; Bell & de Jong 2001; Chabrier 2003; Kroupa 2001). The IMF is the spectrum of stellar types formed given a molecular cloud of certain initial mass, metallicity and other relevant properties. It therefore has the scope to vary from galaxy to galaxy. In addition to the IMF, stellar population synthesis models have other sources of uncertainty (Conroy et al. 2009, 2010; Conroy & Gunn 2010).

In standard fitting of rotation curves of highly inclined disk galaxies, it is common to invoke the “maximum disk hypothesis”, (van Albada & Sancisi 1986; Sackett 1997; Courteau & Rix 1999) i.e. that the stellar disk contributes maximally to the rotation curve. This hypothesis is supported by observations that deduce the microlensing optical depth in the Milky Way (e.g., Bissantz & Gerhard 2002), the baryonic Tully-Fisher relation (McGaugh & Schombert 2015) and measurements that place the co-rotation radius of barred galaxies just beyond the end of the bar (e.g., Sellwood & Debattista 2014). It yields values for the that are typically in accordance with the predictions of stellar population synthesis models.

This, however, does not demonstrate that the hypothesis is correct (see e.g. Herrmann & Ciardullo 2009; Dutton et al. 2011) and it would be ideal to have a robust, independent measurement that breaks the disk mass degeneracy. This is theoretically possible because the can be determined dynamically through the vertical velocity dispersions of stars over the full projected area of the galactic disk, if the scale-height is known (Bahcall 1984). For close to edge-on disk galaxies, however, one cannot measure the vertical velocity dispersion and thus the technique is limited to only those galaxies with moderate inclinations to the line of sight (), where is seen as optimal (Bershady et al. 2010b; hereafter DMSii).

In order to break the degeneracy between the DM halo and , the DiskMass Survey (DMS; Bershady et al. 2010a; hereafter DMSi) made observations of the line of sight velocity dispersion profiles of 46 nearly face-on disk galaxies (DMSi): 30 of which have been published and a further 100 are part of the larger survey. They also measured their surface brightness profiles and rotation curves. There were two further, essential ingredients in the analysis that come from scaling relations. The first is the luminous Tully-Fisher (TF) relation between the absolute magnitude of a galaxy and a measure of its outer rotation speed. This helps to isolate the inclination of the galaxy, which for nearly face-on galaxies can otherwise be obtained using the tilted disk method of Andersen & Bershady (2013).

The second is the correlation between the disk stellar scale-height, a measure of the thickness of the disk, and scale-length, a measure of the radial extent. This relationship is explored in detail by DMSii (their §2.2), and is derived from the studies of Kregel et al. (2002) (hereafter K02) and Pohlen et al. (2000); Schwarzkopf & Dettmar (2000); Xilouris et al. (1997, 1999). This allowed Bershady et al. (2011) and Martinsson et al. (2013a, b) (hereafter DMSvi and DMSvii) to infer the of the disk. The data imply that the stellar disks are “sub-maximal” (K-band or lower, see Swaters et al. 2014) which means they do not contribute maximally to the rotation curve in the central regions, contrary to the value found from population synthesis models that assume a Kroupa IMF (; McGaugh & Schombert 2014). This leaves more room for DM in the central regions.

Modified Newtonian Dynamics (MOND Milgrom 1983, see Famaey & McGaugh 2012 for a recent review) is a theory which proposes a modification of dynamics whose impact is most apparent in regions of low acceleration. Most current working versions of MOND consist of an actual modification of gravity, i.e., at the classical level, a modification of the Newtonian Poisson equation (but see also Milgrom 2011). This modification occurs due to the hypothesised existence of a new constant of physics with dimensions of acceleration, . For accelerations much stronger than this threshold, , there is no discerned deviation from Newtonian gravity. However, far below the threshold the true acceleration perceived by a test mass is found from - where is the expected Newtonian gravitational field.

Nipoti et al. (2007) and Bienaymé et al. (2009) made studies of the vertical dynamics of the Galaxy in MOND, showing that it could be possible to distinguish between MOND and the DM paradigm with data from the Milky Way. The extra constraint on the dynamics from vertical velocity dispersions in nearly face-on disks of external galaxies provides a new test of this hypothesis. This article addresses whether MOND can simultaneously account for the measured vertical velocity dispersions and rotation curves, while keeping in line with galaxy scaling relations.

In section 2 the framework is presented for the joint modelling of galaxy rotation curves and stellar vertical velocity dispersions in the MOND context. In section 3 the methods are discussed and this includes a discussion of the bulge-disk decomposition, the accuracy of the Poisson solver, and the observational error budget for the main data. In section 4 the primary results are presented, this includes a discussion of the fits to the vertical dynamics and rotation curves, the confidence ranges of the fitted parameters, and how well the fitted parameters mesh with other observations. In section 5 possible scenarios that could alter the results are discussed. In section 6 conclusions are drawn and their implications are explored.

## 2 Dynamical analysis of the DiskMass survey

### 2.1 Rotation curve fitting

The following reviews how to fit the measured rotation curve of a disk galaxy in an idealised case. This is done in order to expose the free parameters and it is generalised for both DM and MOND.

The total model rotation speed is required from the total model potential in order to compare with the measured rotation curve. This can be found from

(1) |

Here, is the cylindrical radius in the disk mid-plane and is the total gravitational potential. Next, the total potential is required from the total mass distribution. This is computed via the Poisson equation, which in Newtonian dynamics is

(2) |

where is Newton’s gravitational constant, is the total mass density from all sources (see the following and §2.1.1) and is the total Newtonian potential. In Newtonian dynamics, is fully equivalent to , but consists of and . In MOND, is fully equivalent to , but a second step is made to find from , which is (Milgrom 2010)

(3) |

where is the acceleration threshold of MOND and is an interpolating function, chosen here to depend on its argument as

(4) |

where is the simple -function and is the standard -function (Famaey & McGaugh 2012 Eqs. 51 and 53), or likewise by

(5) |

#### Baryonic density

The last question here is the distribution of . The derivation of the contribution from atomic and molecular gas is described in DMSvii. Both are assumed to have non-smooth, axisymmetric, radial surface densities. Both gas components are included in the modelling, but they are considered fixed in mass and distribution. They are given nominal scale-heights of 200 pc, to which reasonable variations are inconsequential.

There are typically two stellar components, a bulge and a disk. Below, the bulge and disk surface brightnesses are given to expose the free parameters in their fitting. The bulge is assumed to be spherical and to follow a Sérsic profile

(6) |

where the effective surface brightness (), the effective projected radius () and the Sérsic index () can be fitted to the observed surface brightness distribution.^{1}

(7) |

where is the total luminosity of the disk, is the scale-length and is the scale-height. The vertical distribution is characterised by the exponential function, but the commonly used function (van der Kruit & Searle 1981; Bottema 1993) would be equally appropriate. Here is the scale-height which corresponds to at large . Using the vertical distribution does not change the conclusions (see §5.4).

In general, the surface density can be found by integrating along the line of sight. For a face-on galaxy, the luminosity density of Eq 7 can be projected to give the surface brightness

(8) |

Here the density profiles are always assumed to be smooth. To find the mass density of the bulge and disk, the luminosity densities of the bulge and the disk must be multiplied by their respective so that

(9) |

where and are the stellar values of the bulge and disk respectively. The total density of baryons is then , where is the atomic and molecular gas density. As stated previously, in MOND, is equivalent to the total mass density, (Eq 2), because there is no DM in MOND galaxies. For Newtonian gravity .

For the baryonic mass models, there are 8 parameters: . Of those 9, the surface brightness parameters are either directly observed or unambiguously fitted to the surface brightness profile. This leaves only . Since the DMS sample is chosen so that the total bulge luminosity to total disk luminosity is low, is relatively insignificant and is never independently varied in the modelling performed here (). Therefore, and are the only two free parameters that are relevant to a theory like MOND. In principle, the inclination of the galaxy also has a small amount of freedom but it is strongly curtailed by the luminous TF relation. The aforementioned free parameters are fitted for through a simultaneous comparison of the model vertical velocity dispersions and rotation curves with the observed ones, as is described in §3.3.

#### Inclination

The derivation of the rotation curve of a moderate or high-inclination disk galaxy from the measured 2D velocity field, permits the fitting of tilted rings (see Begeman 1989; van der Hulst et al. 1992). These tilted rings allow us to model the variation in the inclination and position angle of numerous concentric annuli at different galactocentric radii. The inclinations of the various rings, as a function of radius, can vary by 10 (e.g. de Blok et al. 2008) depending on the quality of the data, the regularity of the velocity field and characteristic inclination of the disk.

The DMS galaxies have low inclinations (they are close to face-on), thus the inclination has a lot of leverage on the inferred rotation speed because of the shape of the sine function. It is possible to derive accurate kinematic inclinations for nearly face-on disks using the tilted disk (as opposed to tilted ring) technique of Andersen & Bershady (2013). It is also possible to infer the inclination using the luminous TF relation (Verheijen 2001). This relates the absolute K-band magnitude of the galaxy, , to a measure of the outer rotation speed, such that

(10) |

Anderson et al. (2015, in prep) have determined that the kinematic and TF inclinations for the DiskMass galaxies generally agree well, although there are some outliers.

Relating the measured, inclined outer rotation velocity with the expected outer velocity from the TF relation (Eq 10) allows the expected inclination to be deduced. This inclination is only that expected for the outer parts of the rotation curve, and thus the inner parts can vary somewhat due to a warp.

In addition to the Luminous TF relation, there is a baryonic TF relation (McGaugh et al. 2000; McGaugh 2005) which relates the total baryonic mass of a galaxy to its outer, flat rotation speed, where .

This relation is fundamental to MOND and the exponent and the constant of proportionality are predictions which agrees well with the observed relation (McGaugh 2005). Thus, the MOND baryonic TF relation can be written in a similar form to Eq 10 as

(11) |

Here, and are the fractions of the total luminosity contributed by the bulge and disk respectively. The absolute magnitude of the Sun in the K-band is (Blanton et al. 2003) and is the gas mass.

It was found that the inclinations from Eq 11 are typically between 5-15% larger than those found with Eq 10, depending on the used (here the was taken to be between 0.6 and 1): a smaller implies a smaller corresponding rotation velocity, hence a larger inclination.

For a large enough sample there should, in principle, be no systematic deviation from either TF relation. When modelling rotation curves in general, it is not always clear when the rotation curve has reached the terminal velocity, so some margin of error must be granted. Since Eq 10 has no dependence on , and to make it easier to compare with the DMS results, the luminous TF relation inclinations (Eq 10) are used in this article.

### 2.2 Stellar vertical velocity dispersions

#### Choice of stellar velocity dispersion ellipsoid parameters

In addition to the measured rotation curve, the DMS also measured the line of sight velocity dispersion profile of the stars over the full projected area of the disk. This is then azimuthally averaged to give a 1D line of sight velocity dispersion, which can be converted to a vertical velocity dispersion () through the equation (Westfall et al. 2011)

(12) |

Here, the inclination of the galaxy to the line of sight is again, , and & provide information about the stellar velocity dispersion ellipsoid. Generally, and are expected to take on certain values from measurements in the Solar neighbourhood (Binney & Merrifield 1998; Gerssen & Shapiro Griffin 2012), and is presumed to take on specific values from the epicycle approximation. However, beyond the Milky Way their variation is not empirically well known (see discussion in sect 2.1 of DMSii, and also Westfall et al. 2011; Westfall 2015; Gentile et al. 2015). By choosing galaxies that are nearly face-on, the DMS reduces their importance (cf. Fig 5). The statistical variation of these parameters is discussed in DMSii §2.1, and the DMS analysis establishes and . The mean values used by the DMS are chosen for the default values in this article.

#### Model Stellar vertical velocity dispersions

In order to compare with the observations, the model vertical velocity dispersions of the galaxies must be computed. The vertical velocity dispersion at a height above the mid-plane, at a radial distance from the centre of the disk galaxy is found from (see Nipoti et al. 2007)

(13) |

and the equivalent of the observed vertical velocity dispersion at any radius, R, weighted by the local stellar surface density is given by

(14) |

These equations effectively reduce to

(15) |

In the Newtonian gravity framework, Eq 15 depends mainly on 2 fitted parameters: and . As stated previously, inclination could also be varied but only in a tight range around the TF relation values. Thus, combining simultaneous fits to the observed rotation curves and vertical velocity dispersions is a strong test of the MOND paradigm.

It is worth noting that the infinite potential well of isolated galaxies in MOND is irrelevant here since the vertical gravitational field in Eq 13 is convolved with the exponentially declining stellar density, and thus the MOND gravity is only relevant where there are stars.

#### The DiskMass Survey Method

In the analysis of DMSvi, the measured rotation curve of each galaxy is used to fit the DM halo parameters and then the vertical velocity dispersion is used to directly give the mass surface density using , where is assumed to be 1.5 to describe an exponential vertical stellar distribution. From this surface density, the gas disk was subtracted. This left the stellar disk surface density, . This also includes an unknown component of DM. The stellar as a function of radius was given by . An average of this out to a given radius then defines the quoted disk . This approach, although straight-forward, is only accurate when the derivative of the rotation curve is small. It is also not transferable to MOND because of the non-linearity of the theory and that the in MOND affects both the rotation curve and vertical velocity dispersion. It is more secure to make the reverse calculation and go from observed surface brightness, sample a to give surface density, then use Eq 15 to give a model vertical velocity dispersion which can be compared with the observed vertical velocity dispersions. However, this requires calculation of the full three dimensional potential from a model galaxy - which is performed here.

## 3 Preliminary Modelling

### 3.1 Bulge-Disk decomposition

The inclination corrected surface brightness profiles presented by DMSvi were analysed with a Markov Chain Monte Carlo (MCMC) approach to fit bulge plus disk surface brightness models using Eqs 6 and 8. The seeing was simultaneously accounted for with a Gaussian convolution of appropriate radius (see tables 3 and 4 of DMSvii). All 5 surface brightness parameters (2 disk and 3 bulge) were fitted, but the bulge surface brightness parameters were found to be degenerate with each other due to seeing effects and the lack of data points at low radii. Thus, the bulge parameters were fixed at the maximum likelihood values, which allowed better sampling of the more important disk parameters. In Figs 11-13 the posterior probability distributions are presented for the two disk surface brightness parameters of each galaxy. For the disk scale-length, the range found by DMSvi (red curve) is also plotted.

The best fits to the surface brightness profiles are given in Figs 14-16 and are generally good. Both linear radius and log radius are plotted on the -axis to expose the quality of the fits to the bulge and the outer disk. The fitted scale-lengths found here (see Figs 11-13) are generally consistent with those fitted by the DMS, although theirs have smaller error ranges.

Once the bulge and disk surface brightness profiles were established, an Abel transform was used to de-project the bulge surface brightness, which was stored numerically. The disk has cylindrical symmetry, so it does not require de-projection. Following this, rejection sampling was used to generate N-particle representations of the galaxies, with half of all particles representing the stellar disk and the other half split evenly between the stellar bulge and the atomic and molecular gas disks. The fitted scale-lengths are used for all the following analysis, but they are not used to generate the N-particle representations of the stellar disk. Instead, the “observed” surface brightness of the disk is sampled after subtracting the fitted bulge - since the fit to the surface brightness can be poor at large radii.

In this analysis the static N-particle representations are only used to find the potential and gravitational field, there is no N-body evolution of the simulated galaxies. The masses recovered after generating the N-particle representations of each component of each galaxy generally agreed very well with the masses reported by the DMS.

### 3.2 The MOND Poisson solver

The MOND Poisson solver described in Angus et al. (2012) is used to compute the radial and vertical gravitational fields of the N-particle galaxy models. The code solves the modified Poisson equation of the Quasi-linear version of MOND (QUMOND) given by Eq 3. The code uses a 3D grid and the cloud-in-cell technique to numerically discretise the 3D density of an N-particle distribution. It then employs finite differencing and multigrid methods to iterate from a test potential to the final potential which accurately reflects the density.

#### Comparison of theoretical and numerical vertical gravity

The Poisson solver is used to compute the gravitational field of the baryons in the radial direction, , as a function of radius, to find the model rotation speed (Eq 1). Simultaneously, the gravitational field of the baryons is solved for in the vertical direction, , as a function of height above the disk, at several discrete radii: to 9.5 kpc in steps of 1 kpc and then , 14.5, 18.5 and 25.5 kpc. In MOND, both gravitational fields noted above are equivalent to the gradient of the total potential (i.e. ) used in Eqs 1 and 13 and thus the rotation speed can be straight-forwardly calculated and Eqs 13 and 14 can be integrated to find at different radii, .

The most important, and difficult to compute, quantity is the vertical gravity profile at large radii. For an isolated, double exponential disk, like that introduced in Eq 7, the Newtonian vertical gravity profile at a given radius can be calculated numerically - as described in detail by Kuijken & Gilmore (1989).

In the top panel of Fig 1 the Newtonian vertical gravity profile is plotted for three discrete cylindrical radii , 8.5 and 18.5 kpc. The left and right hand panels correspond to two different double exponential disk galaxies (as per Eq 7). The left hand panel has scale-length kpc, scale-height kpc and total mass . The right hand panel has scale-length kpc, scale-height kpc and total mass . The computation from the Poisson solver using Newtonian dynamics is the thick black line and the theoretical value using Eq 27 of Kuijken & Gilmore (1989) is the red line. The agreement is generally very good, except for small heights above the disk. This is due to a combination of the limit of spatial resolution at large radii, owing to the centrally refining mesh, and limited particles at large radii. The average percentage error for these discrete radii (3.5, 8.5, 11.5, 14.5, 18.5 and 25.5 kpc) between the numerical calculation of (using Eq 15) and the theoretical values are 0.05, 0.8, 0.7, 5.8, 5.2 and 4.4%. Given that only 6 of the 30 galaxies have vertical velocity dispersion data points beyond 11 kpc and the errors on those data points are typically more than 10%, this is more than adequate.

In the bottom two panels of Fig 1 the Newtonian (black lines) and MOND (; turquoise lines) vertical gravity profiles are compared, both of which were computed with the Poisson solver. The MOND profiles are strongly boosted relative to the Newtonian profiles, however, the higher surface density galaxy (left hand panel) receives less of a boost than the other.

### 3.3 Error budget

The extra information that allows the DMS to close their set of equations is that, in general, disk scale-heights are observed to correlate with their scale-lengths. These scale-heights and scale-lengths have been fitted to edge-on galaxies and certain assumptions are made in their modelling (like constant inclination, that the disks are well described by exponential or distributions). The DMS made an analysis of literature measurements and derived a simple relationship between and such that in units of kpc (see DMSii)

(16) |

This relation has a 1 scatter of roughly 25% (DMSii).

In the analysis presented here, the scale-height and disk are fitted to the observed vertical velocity dispersions and rotation curves of the sample of 30 galaxies. After the distribution of fitted scale-height and scale-lengths is known, they will be compared for consistency with the direct observations of and for the sample of edge on galaxies compiled by the DMS. Therefore, a simultaneous fit must be made to the observed vertical velocity dispersions and rotation curves. A critical concern is the contribution of each of the two datasets to the overall likelihood. Since the vertical velocity dispersions are the primary data set, the errors on each data point as taken as computed by the DMS. Given that there is the potential for some systematic errors from disk warping, the rotation curve should not be given too large a weighting, especially the inner parts. The decision was made to increase the error bars on the rotation curve data points to 10 . Lastly, the two separate reduced s from fitting the vertical velocity dispersion profile and the rotation curve are combined.

Therefore, the likelihood is given by

(17) | |||||

where and are the number of relevant data points in the rotation curve and vertical velocity dispersion profile respectively. The prior then multiplies to give the un-normalised posterior probability.

## 4 Primary results

### 4.1 Parameter posterior probabilities

MCMC sampling was used to find the posterior probability distribution for each of the free parameters: disk scale-height (), disk () and inclination (). The scale-height was varied by considering the ratio of the fitted scale-height to the one derived from observations (Eq 16). A broad Gaussian prior of width 1.5 was placed on this ratio. The prior placed on the disk was centred on and had a width 0.5 . The ratio between the fitted inclination and the luminous TF relation inclination (from Eq 10) received a fairly tight Gaussian prior of 0.15 (or 15%), since there should be little deviation from the TF relation. The main reason for including inclination as a free parameter is just in case there is a sharp increase in posterior probability for a small change in inclination. The full expression for the likelihood is given by Eq 17. This defines how the goodness of fit for the two data sets is combined: rotation curve and vertical velocity dispersion.

The re-normalised posterior probability distributions of each of the three parameters (scale-height, inclination and ) is plotted in Figs 17-19 for each galaxy individually. There are two main lines in each of the three columns: MOND with and 2 (see Eq 4) black lines, solid and dashed respectively. The second column (inclination) also has the curve of the prior (turquoise). The third column (M/L) has a green line which represents the Newtonian gravity (with DM halo) found by the DMS (DMSvi), shown only for reference. All the lines are re-normalised to have the same maximum posterior probability. The two MOND fits vary little in terms of goodness of fit.

From the left hand column it is clear that the majority of the fits require substantially lower scale-heights than those derived from observations (Eq 16), which in those panels are unity. There is a preference with most galaxies to have a higher inclination because that decreases the amplitude of the rotation curve. A lower rotation curve requires a lower . A lower allows a larger scale-height - which allows better agreement with the observations of edge-on galaxies (Eq 16). However, the fairly tight prior on the inclination prevents it from changing substantially. Usually it changes by less than 10%.

In Table 1 various 1 confidence ranges for the fitted scale-heights and stellar mass-to-light ratios are presented using MOND with varying constraints (from interpolating functions and which parameters are left free). The preferred scale-heights, those from Eq 16 and used by the DMS, are given for each galaxy (column 2). For reference (column 3) the found by the DMS when using Newtonian gravity and DM halos is shown, followed by the best fit MOND scale-heights (as a ratio between the fitted scale-height and the one derived from observations of edge-on galaxies - Eq 16) and the for three different interpolating functions (the two used above and one other). These fits are for the scenario where the scale-height and are barely constrained, but the inclination is tightly constrained by its prior.

#### Comparison of the MOND fits with the DMS data

In Figs 23-25 the best fits to the vertical velocity dispersion (first and third columns) and rotation curve (second and fourth columns) are plotted in the MOND case with . The two interpolating functions are not plotted together because different best fit inclinations mean the data points vary in each case. There is little difference between the two MOND cases in terms of quality of fit. Most fits to the vertical velocity dispersion are good, but the quality of fits to the rotation curves range from poor to good. Rotation curves with good fits include the galaxies UGC 448, 1081, 1087, 1908, 3091, 4036, 4368, 4380, 7244, 6918, 8196, 9965, 11318, 12391. Quite good fits are found for UGC 463, 1529, 1635, 1862, 3701, 4256; whereas the fits to UGC 3140, 3997, 4107, 4458, 4555, 4622, 6903, 7917, 9177, 9837 are poor.

MOND fits to the rotation curves of nearly face-on galaxies are very sensitive to the inclination (see e.g. de Blok & McGaugh 1998). In the outer parts, they depend on the fourth power of . Warping of the disk might be invoked to improve certain rotation curve fits, however, in some cases the DMS galaxies have already been corrected for warps.

To compare with models that use the scale-height from Eq 16, the vertical velocity dispersion (red line) is over-plotted for each galaxy using that larger scale-height, but the and inclination of the best fit are kept. In many cases, the red line has a far greater amplitude than the data points and best fit (black line). There are, however, a few cases where the red and black lines are close together, such as UGC 1862, UGC 4368 and UGC 4458 - although the latter only has a single relevant data point.

### 4.2 Comparison of the fitted vs with observations

From fitting the observed surface brightness profiles of the DMS galaxies, a value for the radial scale-length was estimated for each galaxy disk. As can be seen in Figs 11-13 (second and fourth columns) the radial scale-length fitted by DMSvi is generally slightly larger than the best fit found here, but the majority are consistent within the errors. This fact makes the analysis presented here slightly more favourable to MOND. The difference in fitted scale-lengths is in part due to the fact that in DMSvi a finite central region was used to fit the scale-length, whereas the full extent of the galaxy is used here. Thus, for any galaxies whose disk surface brightness is not well fit by a single exponential there is a difference. For example, UGC 4036 has a shallow inner surface brightness profile and then a steep outer decline, thus the scale-length found here naturally has a smaller scale-length. Recall that the N-particle realisations (§3.1) use the observed surface brightness, not the fitted exponential disk.

Similarly, by fitting the vertical velocity dispersion profile and the rotation curve of each DMS galaxy, a confidence range for the scale-height of each galaxy has been derived. In principle, this combination of and parameters can be compared to the measurements of scale-lengths and scale-heights of a sample of edge-on galaxies. In Fig 2 the vs diagram for the MOND fits (, ) is plotted along with the observations of these parameters for edge-on galaxies made by K02.

In order to visualise this more plainly, contours in the vs plane were added to Fig 2. Here, the re-normalised posterior probability density of all individual data points are co-added on a fine grid, with each galaxy receiving equal weighting. A point with larger error bars will spread its probability density over a larger range. The blue contours represent the co-added MOND probability densities and the black triangles give the location of each individual galaxy point. The error bars are not included to minimise clutter. The measurements by K02 have red shaded contours with green circles. The MOND contours are significantly inclined relative to the K02 contours. The equivalent plot for the interpolating function is almost identical.

### 4.3 Two dimensional Kolmogorov-Smirnov test

To demonstrate statistically the the large offset in the vs plane between the K02 data points in Fig 2 and the fitted MOND points, a 2D Kolmogorov-Smirnov (KS) test (see Press et al. 1992) chapter 14.7) was employed. In this test, the first K02 data point was taken and made the origin of four quadrants. Then the fraction of K02 points in each quadrant and the fraction of MOND points in the same quadrants. The largest difference between the two fractions in a single quadrant is stored and then the rest of the K02 points are cycled through, making each one the origin in turn. The largest of the largest differences is stored and then the procedure is repeated centring on the MOND points. The largest difference in fractions is ordinarily the same regardless of whether the K02 or MOND points are central. This number is the statistic of merit for the 2D KS test and provides a significance level that the two samples originate from the same parent population.

Since the data points have fairly large error bars, if the mean of a particular data point lies in the first quadrant, it is not ideal to add the whole point to that quadrant. Instead, the fraction of the point inside each quadrant determined by the point’s error bars is added.

According to the 2D KS test, the significance levels that the MOND fits with and 2 come from the same parent distribution as the K02 sample are and respectively. Using the interpolating function the significance level is . Comparing with DM models, it is easy to obtain significance levels of order unity, so the method seems robust.

As a separate test, a straight line was fitted to each data set of vs with a fixed intercept. The resulting gradients differ at the 2 level, regardless of interpolating function. Using a variable intercept only improves the agreement slightly.

### 4.4 Superposition of distributions

In the approach presented here to modelling the DMS data the is fitted for as well as the scale-height. At near infrared wavelengths a smaller variation of the is expected from one galaxy to the next than at optical wavelengths (Bell & de Jong 2001; Bell et al. 2007). To find the posterior probability distribution of across the full sample of galaxies, the summation of the re-normalised posterior probability density for the of all 30 galaxies is plotted in Fig 3. This is done for three MOND interpolating functions and the DMS computation of in the DM scenario. This means the same weighting is given to each galaxy’s posterior probability density in the summation. Recall that in the fitting process a Gaussian prior is used with width 0.5 centred on .

The MOND interpolating function (solid black line) has an average best fit , whereas the MOND interpolating function (dashed black line) requires significantly larger values, . Bear in mind that the produces less of a boost to the gravity than the interpolating function at intermediate gravities around . The interpolating function (which uses Eq 5) has a similar distribution as . There is no impact on the derived values from the weak prior. The green line, given only for reference, is the derived from the DMS analysis with Newtonian gravity and DM.

### 4.5 Why does MOND need such thin disks?

The results presented here suggest that MOND requires disks that have roughly half the vertical scale-height as those inferred from observations of edge on galaxies (DMS II or Eq 16). There are two key differences between the MOND vertical velocity dispersions and those using DM. The first effect is related to the required of the galaxies. Comparing a disk with the same surface brightness in MOND and Newtonian gravity, the Newtonian disk is surrounded by a dark matter halo. This dark matter halo can account for the radial gravitation field that produces the observed rotation speed. This means the of the Newtonian disk is only weakly constrained by the rotation curve (van Albada et al. 1985). In MOND, the is essentially fixed by the rotation speed because it is the only parameter that can be varied to allow a match between the observed and expected rotation speeds. Based on the results presented here, the required of the disk is a factor of between 1.5 and 2.5 higher in MOND (cf. Fig 3), depending on the interpolating function used. Since the is typically much larger in MOND, the resulting vertical gravitational field is larger than the dark matter equivalent. Moreover, even if the disks in MOND and Newtonian gravity had the same , the MOND gravitational field is still boosted relative to the Newtonian one as seen in Fig 1. This follows because the nearly spherical dark matter halo typically has a weak influence on the vertical gravitational field within the disk.

The influence of MOND on the vertical gravitational field, relative to Newton, for the same disk is of course greater when the surface density is low. This is demonstrated by comparing the vertical gravitational field of a high surface brightness galaxy (Fig 1 bottom left panel) with a lower surface brightness one (Fig 1 bottom right panel).

In relation to the equations of §2, it can be seen that increasing or increases the amplitude of the model . However, in MOND only increasing increases the model rotation curve at intermediate to large radii. Thus, once is fixed by the rotation curve, only can be varied to fit the observed vertical velocity dispersion. Since the observed vertical velocity dispersions are typically much lower than predicted by MOND, must be decreased to fit them. It is for these reasons, the fitted scale-heights in MOND are significantly lower than those expected from observations (Eq 16).

### 4.6 Central dynamical surface density versus central surface luminosity density

Swaters et al. (2014) used the DMS data to present a correlation between the extrapolated disk central surface luminosity density, from the photometry, and the extrapolated disk dynamical central surface mass density (see §2.2.3). The majority of galaxies are consistent with the relation , but at the low surface luminosity density end a clump of 6 galaxies lie above the relation. In the Newtonian context, this can be attributed to greater relative DM dominance for low surface luminosity density galaxies (Swaters et al. 2014), but this is also reminiscent of the MOND transition in rotation curves of spiral galaxies, where the effect becomes pronounced at low surface densities. Thus, a prediction for dynamical surface density versus central surface luminosity density in MOND would be valuable.

To this end, a representative galaxy from the DMS was selected - UGC 4036. The stellar bulge and all gas was removed from the N-particle realisation, thus it has the same stellar disk surface brightness profile as UGC 4036, and the same vertical profile. For this model, the central vertical velocity dispersion, from Eq 15, can be found for a chosen . From here, the surface brightness of the stellar disk can be scaled to compute at a broad range of central surface luminosity densities which allows a comparison with the data from Swaters et al. (2014), which is done in Fig 4.

The data points come from Swaters et al. (2014) Fig 1 (left hand panel) and both coordinates are extrapolated from exponential fits to the observed projected disk surface luminosity density and observed vertical velocity dispersion. The red line is the relation . The black solid line is the MOND prediction with and the scale-height expected from observations of edge-on galaxies (DMSii; Eq 16 in this paper). The black dashed line has the same , but used disks that are twice as thin, which MOND has been shown here to require in order to mesh with the DMS data. The two blue dashed lines also use disks twice as thin as observations of edge-on galaxies, with the upper line using and the lower line using . Thus, MOND provides a natural explanation for the deviation of low surface luminosity density galaxies from the central dynamical surface density versus central surface luminosity density relation.

## 5 Secondary results: Varying additional parameters

In this section the impact of other variables on the results from §4 is considered.

### 5.1 Variation of the Stellar Velocity dispersion Ellipsoid

As discussed in §2.2, the stellar velocity dispersion ellipsoid (SVE) parameters are fixed to those used by the DMS i.e. and . To demonstrate the degree to which a reasonable variation in these parameters might affect the results, in Fig 5 the vertical velocity dispersion ratio between three combinations of parameters and the default values are plotted against inclination.

Using different values for in Eq 12 can make a sizeable difference to the derived , especially for large inclinations. The dotted line of Fig 5 shows that increasing to 0.9, which is 2 above the default value used by DMSvi, would increase the derived by around 10% for standard DMS inclinations of . On the other hand, decreasing to 0.48 and using (dot-dashed line in Fig 5), which was the value derived by Westfall et al. (2011) for UGC 463, would decrease the derived by around 15% for .

To test the impact of these SVE parameters on the fits to the vertical velocity dispersions and rotation curves, MCMC modelling was used to simultaneously fit , and , with inclination fixed by the TF relation and set by observations of edge-on galaxies (Eq 16). Here, only flat priors are used such that and . The has a Gaussian prior of width 0.5 centred on 0.3. The most relevant parameter is , due to its impact on Eq 12.

In the top panel of Fig 6 the co-addition (for all 30 galaxies) of the re-normalised posterior probability is plotted as a function of . Each of the 30 galaxies are given equal weighting, regardless of their individual posterior probability. This shows the impact of increasing from the nominal value of 0.6 used here and by the DMS. The posterior probability increases roughly a factor of two and a half from to 1.

This increase does not necessarily mean that the fit qualities are good. To demonstrate this, it is worth comparing the fits where is a free parameter (but scale-height is fixed) to the fits where the scale-height is free. These fits where scale-height is free represent, for the majority of galaxies, a good quality fit to the vertical velocity dispersion. For the remainder of this section, this fit is referred to as the benchmark fit. does not influence the fit to the rotation curve.

In the bottom panel of Fig 6 the -axis is the ratio between the posterior probability of the fit with fixed scale-height (Eq 16) and inclination (luminous TF relation), and best fit and to the benchmark fit. This is plotted against the inclination of the galaxy derived from the luminous TF relation to show that increasing may help with some (but not all) of the high inclination galaxies, but the majority of the sample would have extremely poor fits to the DMS data even with .

### 5.2 Fixed scale-height, unconstrained inclination

A relevant question is how far the fitted inclination would have to deviate from the luminous TF relation before a good fit to the DMS data could be achieved, whilst remaining consistent with the values derived from fits to the photometry of edge-on galaxies. To do this, the scale-height was fixed to the aforementioned value (Eq 16) for each galaxy and fitted only for the inclination and . The parameters and were used. For the sake of brevity, in Fig 7 only the co-added, posterior probability densities for the two parameters are plotted.

The top panel shows the , which is shifted to values roughly three times smaller (cf. dashed line of Fig 3) than the original models with scale-height free. In the last column of Table 1 the range is given for each galaxy in this scenario individually. The values are much smaller than the other models where inclination is constrained. In the bottom panel of Fig 7 the inclination, which is given as a ratio of the fitted value () to the luminous TF relation value (), is plotted. With an unchanged inclination (), the quality of the fits to the majority of the DMS galaxies are poor, but by increasing the fitted inclination to the fits are of a high quality. This is partly owed to the resulting rotation speeds generally having lower amplitudes. Such significant deviations from V01’s luminous TF relation seems unlikely.

To demonstrate this last point, in Fig 8 the luminous TF relation from V01 (Eq 10) is plotted, which is used for the DMS and this study. The solid black line is the best fit TF relation and the black data points are the values of the DMS galaxies. The blue, green and red points are the positions of those galaxies if their inclinations were increased by a factor 1.1, 1.2 or 1.4 respectively. The dashed and dotted lines are the 1 and 2 scatter in the observed TF relation. Clearly even changing the inclinations by a factor of 1.1 would be in stark conflict with the TF relation.

### 5.3 Variable MOND acceleration parameter

The DMS data seem to be at odds with MOND when using the default acceleration parameter . To test this parameter’s influence, models were tested where the scale-height and inclination are fixed, but the MOND acceleration parameter is variable - along with the . Although varying the MOND acceleration parameter from galaxy to galaxy is not an acceptable solution, this could identify a unique value for that is more consistent with the DMS data. In the top panel of Fig 9 another co-added, posterior probability density is plotted, this time against the MOND acceleration parameter. This is done for three different interpolating functions (, 4 and 8 - black, red and green lines respectively). The parameters and 8 are trialled to see if the sharper transition from the MOND regime to the Newtonian regime would affect the ability to fit the DMS data. It appears that little is gained overall by varying and the same is true on a case by case basis. This can be demonstrated by comparing the fit with free (and scale-height fixed) to the benchmark fit (similar to what was done in §5.1 and Fig 6 bottom panel). In the bottom panel of Fig 9, for each galaxy individually, the ratio is plotted between the posterior probability of the best fit using a variable (with fixed scale-height) and the benchmark fit. The vast majority of the galaxies would have poor fits with a variable and fixed scale-height.

Varying suffers from the same drawbacks as varying the . Increasing can allow a fit to the rotation curve with a lower , but this simultaneously increases the vertical gravitational field, which in turn requires the scale-height to be reduced.

In the second (with ) and third last (with ) columns of Table 1 the range is displayed for each galaxy individually when the MOND acceleration parameter is variable.

### 5.4 vertical distributions

Regarding whether the poor results were caused by using an exponential function to describe the vertical scale-height, the DMS data was modelled using a function (van der Kruit & Searle 1981) to describe the vertical distribution of the stellar disk in the N-particle models. The scale-height used in the function is . In Fig 10 the posterior probability densities are plotted for the ratios of the fitted to for the and 2 interpolating functions. These are compared with data from Mosenkov et al. (2010) and Bizyaev & Mitronova (2002). From Fig 5d of Mosenkov et al. (2010) the indicative range of observed is between 12 and 22. A solid vertical red line is placed in Fig 10 to mark the value where their histogram begins to decline and a dashed vertical red line to mark where their data effectively terminates.

The 153 galaxies of Bizyaev & Mitronova (2002) are also used to produce a probability density (green curve in Fig 10) to compare with the posterior probability density found here. The distributions (solid and dashed black lines in Fig 10) fitted here peak beyond where Mosenkov et al. (2010) and Bizyaev & Mitronova (2002) have any evidence of such values. The distribution found here also extends to much larger values. To demonstrate the mismatch statistically, a 2D KS test was performed to compare the and distributions found here with the data of Bizyaev & Mitronova (2002). The confidence levels found for the and 2 interpolating functions were and . Evidently, the vertical distribution does little to aid the agreement between MOND and the DMS data.

## 6 Conclusions and Discussion

By assuming the scale-heights of face-on stellar disks are similar to the measured scale-heights of edge-on galaxies (DMSii), the DMS found that stellar disks are sub-maximal, with K-band mass-to-light ratios of the order of (DMSvi).

In this paper the measured vertical velocity dispersions and rotation curves from DMSvi were analysed in the context of MOND. The problem is that the velocity dispersions measured by the DiskMass Survey can only allow for thick enough disks if the vertical gravitational field is weak. In Newtonian gravity this can be accommodated by lowering the . In MOND, the must be larger to account for the rotation curve in the absence of DM. This means the vertical gravitational field is strong. Moreover, it is amplified by MOND relative to the vertical gravitational field of a disk with the same in Newtonian gravity.

Thus, the most straight-forward way to reduce the model vertical velocity dispersions in MOND is to decrease the stellar scale-height. If these were reduced from the values derived from observations of edge-on galaxies by roughly a factor of 2 then the DMS vertical velocity dispersions and rotation curves would be compatible with MOND. Regardless of stellar scale-length, these disks would have scale-heights between 200 and 400 pc. According to a two-dimensional K-S tests, such thin disks are strongly at odds with observations.

All other key parameters involved in analysing the DMS data in MOND were modelled, such as the stellar velocity dispersion ellipsoid parameters, inclination, the choice of stellar vertical distribution profile, the MOND acceleration parameter and interpolating function. No clear way was found to reconcile MOND with the data.

If the derived rotation velocities were lowered, by imposing 30-40% higher inclinations, it became possible to simultaneously fit the low vertical velocity dispersions and low rotation speeds with low values for the stellar . By design, since the outer rotation speeds are fitted, this is in accordance with the MONDian baryonic TF relation (with logarithmic slope fixed to 4, Eq 11). Despite this, these circular velocities are inconsistent with the luminous TF relation which is quite well established for high-inclination galaxies. To this end, it would be interesting to make a fully self-consistent analysis that models the effect of the varying inclination on the two dimensional velocity fields and the surface brightness profiles.

A prediction was made for the scaling of the central dynamical surface density versus central surface luminosity density in MOND and compared with the DMS galaxies (see Swaters et al. 2014). MOND correctly predicts that the low surface luminosity density galaxies deviate from the simple relation .

In addition to viewing MOND as a modification of the law of gravity, Milgrom (2011) has maintained that a modification of inertia might be better positioned to explain all available galactic dynamics data. In such a theory, the orbital trajectories of particles would have different inertias meaning the interpolating function in the radial direction would differ from the vertical directions. This is an intriguing possibility since the problem with MOND outlined here is that for the radial force defined by a galaxy, the corresponding vertical force is too large. It is interesting therefore that the vertical velocity dispersions of disk galaxies must increase to agree with MOND, but the observed radial velocity dispersions of some dwarf spheroidal galaxies (Angus 2008; Angus et al. 2014) must decrease. This anecdotal evidence should be followed up with more rigorous study.

Something worth further investigation is the impact from super-thin disks, much thinner than the observed scale-heights represented here (Schechtman-Rook et al. 2012; Schechtman-Rook & Bershady 2014), on the measurement of vertical velocity dispersions by the DMS survey. Specifically, how prominent are super-thin disks in all disk galaxies, and are the vertical velocity dispersions that the DMS measures meaningfully influenced by them? The super-thin disk found by Schechtman-Rook & Bershady (2014) was truncated at around 3 kpc, suggesting the influence is minimal. As discussed in DMSi, it is expected that stars in or near the mid-plane diffuse over time. The youngest stars may contribute to the light because of the OB stars, but they contribute little to the absorption lines, because those OB stars have weak spectral features (except for H and He).

It is possible to make even stronger claims about the status of MOND-like theories, as well as conventional dynamics, with certain desiderata. From a theoretical point of view, this would be a detailed study of the impact of the neglected cross-term (tilt of the velocity ellipsoid) in Eq 15. From an observational point of view, a larger sample of edge-on galaxies (like K02) with near-infrared photometry, to more precisely confirm the correlation between scale-length and scale-height. It would also be beneficial to have an increased sample of nearly face-on galaxies, with measured stellar velocity dispersions and rotation curves and with greater sensitivity to be able to measure the stellar velocity dispersions to larger radii.

## 7 acknowledgements

GWA is a postdoctoral fellow of the FWO Vlaanderen (Belgium). AD acknowledges partial support from the INFN grant Indark, the grant Progetti di Ateneo TO Call 2012 0011 âMarco Poloâ of the University of Torino and the grant PRIN 2012 “Fisica Astroparticellare Teorica” of the Italian Ministry of University and Research. The authors thank M. Bershady, K. Westfall and F. Lelli for valuable comments and are deeply indebted to the DiskMass Survey team, and Thomas Martinsson in particular, for taking and sharing their data. The reviewer also suggested several important improvements to the article.

DMS | MOND | ||||||||||

and free; constrained | and free | and free | |||||||||

Galaxy | (kpc) | ||||||||||

U00448 | 0.460.10 | 0.310.21 | |||||||||

U00463 | 0.450.10 | 0.290.17 | |||||||||

U01081 | 0.400.09 | 0.400.16 | |||||||||

U01087 | 0.410.09 | 0.320.13 | |||||||||

U01529 | 0.440.10 | 0.290.16 | |||||||||

U01635 | 0.390.09 | 0.260.12 | |||||||||

U01862 | 0.240.06 | 0.560.29 | |||||||||

U01908 | 0.530.12 | 0.260.18 | |||||||||

U03091 | 0.440.10 | 0.230.10 | |||||||||

U03140 | 0.430.10 | 0.330.14 | |||||||||

U03701 | 0.440.11 | 0.850.40 | |||||||||

U03997 | 0.580.14 | 0.450.19 | |||||||||

U04036 | 0.490.12 | 0.290.11 | |||||||||

U04107 | 0.410.09 | 0.280.14 | |||||||||

U04256 | 0.520.12 | 0.150.19 | |||||||||

U04368 | 0.410.10 | 0.570.34 | |||||||||

U04380 | 0.540.12 | 0.270.09 | |||||||||

U04458 | 0.790.17 | 0.420.24 | |||||||||

U04555 | 0.480.11 | 0.350.20 | |||||||||

U04622 | 0.700.16 | 0.280.17 | |||||||||

U06903 | 0.490.11 | 0.240.17 | |||||||||

U06918 | 0.210.05 | 0.060.24 | |||||||||

U07244 | 0.460.11 | 0.410.18 | |||||||||

U07917 | 0.760.17 | 0.290.33 | |||||||||

U08196 | 0.530.12 | 0.940.41 | |||||||||

U09177 | 0.670.15 | 0.360.26 | |||||||||

U09837 | 0.600.14 | 0.320.14 |