<figure><img src="image_1.png" /></figure>

# Dark Matter in Ultra-diffuse Galaxies in the Virgo Cluster from Their Globular Cluster Populations

Elisa Toloba¹, Sungsoon Lim²,³, Eric Peng²,³, Laura V. Sales⁴, Puragra Guhathakurta⁵,⁶,

J. Christopher Mihos⁶, Patrick Côté⁷, Alessandro Boselli⁸, Jean-Charles Cuillandre⁹, Laura Ferrarese⁷, Stephen Gwyn⁷, Ariane Lançon¹⁰, Roberto Muñoz¹¹, and Thomas Puzia¹¹

¹Department of Physics, University of the Pacific, 3601 Pacific Avenue, Stockton, CA 95211, USA; etoloba@pacific.edu

²Department of Astronomy, Peking University, Beijing 100871, People's Republic of China

³Kavli Institute for Astronomy and Astrophysics, Peking University, Beijing 100871, People's Republic of China

⁴Department of Physics and Astronomy, 900 University Avenue, Riverside, CA 92521, USA

⁵UCO/Lick Observatory, University of California, Santa Cruz, 1156 High Street, Santa Cruz, CA 95064, USA

⁶Department of Astronomy, Case Western Reserve University, Cleveland, OH 44106, USA

⁷National Research Council of Canada, Herzberg Astronomy and Astrophysics Program, Victoria, BC V9E 2E7, Canada

⁸Aix Marseille University, CNRS, LAM, Laboratoire d'Astrophysique de Marseille, Marseille, France

⁹CEA/IRFU/SAp, Laboratoire AIM Paris-Saclay, CNRS/INSU, Université Paris Diderot, Observatoire de Paris, PSL Research University, F-91191 Gif-sur-Yvette Cedex, France

¹⁰Observatoire Astronomique de Strasbourg, Université de Strasbourg, CNRS, UMR 7550, 11 rue de l'Université, F-67000 Strasbourg, France

¹¹Institute of Astrophysics, Pontificia Universidad Católica de Chile, Av. Vicuña Mackenna 4860, 7820436 Macul, Santiago, Chile
Received 2017 October 27; revised 2018 February 23; accepted 2018 March 10; published 2018 March 29

## Abstract

We present Keck/DEIMOS spectroscopy of globular clusters (GCs) around the ultra-diffuse galaxies (UDGs) VLSB-B, VLSB-D, and VCC615 located in the central regions of the Virgo cluster. We spectroscopically identify 4, 12, and 7 GC satellites of these UDGs, respectively. We find that the three UDGs have systemic velocities ($V_{sys}$) consistent with being in the Virgo cluster, and that they span a wide range of velocity dispersions, from ~16 to ~47 km s⁻¹, and high dynamical mass-to-light ratios within the radius that contains half the number of GCs (407±916, 21±15, 60±65, respectively). VLSB-D shows possible evidence for rotation along the stellar major axis and its $V_{sys}$ is consistent with that of the massive galaxy M84 and the center of the Virgo cluster itself. These findings, in addition to having a dynamically and spatially (~1 kpc) off-centered nucleus and being extremely elongated, suggest that VLSB-D could be tidally perturbed. On the contrary, VLSB-B and VCC615 show no signs of tidal deformation. Whereas the dynamics of VLSB-D suggest that it has a less massive dark matter halo than expected for its stellar mass, VLSB-B and VCC615 are consistent with a ~10¹² M☉ dark matter halo. Although our samples of galaxies and GCs are small, these results suggest that UDGs may be a diverse population, with their low surface brightnesses being the result of very early formation, tidal disruption, or a combination of the two.

**Key words:** galaxies: clusters: individual (Virgo) – galaxies: evolution – galaxies: formation – galaxies: individual (VLSB-B, VLSB-D, VCC615) – galaxies: kinematics and dynamics

## 1. Introduction

Ultra-diffuse galaxies (UDGs) are extremely low surface brightness galaxies (central surface brightness $\mu_{g,0} \gtrsim 24$ mag arcsec⁻²) with luminosities in the dwarf galaxies regime ($M_V \gtrsim -16$), and sizes in the massive galaxies regime (half-light radius $R_e \gtrsim 1.5$ kpc). UDGs are characterized by spheroidal shapes, nearly exponential surface brightness profiles, and quenched stellar populations (Mihos et al. 2015, 2017; van Dokkum et al. 2015).

Large low surface brightness galaxies were found for the first time in the Virgo cluster photographic plates by Sandage & Binggeli (1984) and Binggeli et al. (1987). Later on, other studies found a few more of these diffuse galaxies (Impey et al. 1988; Dalcanton et al. 1997; Caldwell 2006). However, with the new deep imaging surveys, a plethora of these systems are being found mainly in cluster environments (Koda et al. 2015; Mihos et al. 2015; Muñoz et al. 2015; van Dokkum et al. 2015; Martínez-Delgado et al. 2016; Toloba et al. 2016b; van der Burg et al. 2016; Janssens et al. 2017; Mihos et al. 2017; Román & Trujillo 2017; Venhola et al. 2017).

There are three main possible mechanisms that could explain the observed properties of the UDGs. (1) They could be extended dwarf galaxies. Some simulations predict that they are rapidly rotating (Amorisco & Loeb 2016), while others suggest that their extended sizes are the result of strong gas outflows (Di Cintio et al. 2017). (2) They could be tidal galaxies formed from the debris of harassed and ram pressure stripped galaxies that lost large fractions of stars. In these two scenarios, the UDGs are expected to have shallow potential wells, which makes them vulnerable to the cluster environment (e.g., Moore et al. 1996). However, UDGs are found in extremely dense regions, such as the cores of the Virgo and Fornax clusters (Mihos et al. 2015, 2017; Muñoz et al. 2015) and the Coma cluster (van Dokkum et al. 2015). They could be falling in the cluster for the first time. (3) They could be “failed” massive galaxies, where the environment and/or internal feedback stopped the star formation and, as a result, the number of stars is smaller than expected for that size. In this scenario, UDGs have a massive dark matter halo that makes them less prone to disruption.

The large number of globular clusters (GCs) found in some UDGs (van Dokkum et al. 2017) points to the third scenario
given that these numbers are more typical of massive galaxies than dwarfs (Peng et al. 2006). However, an analysis of a larger sample of UDGs suggests that they do not have a statistically significant excess of GCs compared to normal dwarf galaxies of the same stellar mass (Amorisco et al. 2018).

Measuring the dark matter halo would help to distinguish between formation scenarios. A massive dark halo can explain their survival in high density environments and their origin as “failed” massive galaxies gets stronger. A low-mass dark halo would suggest that UDGs are puffed up dwarf galaxies that are likely on the verge of disruption. However, if disruption is currently happening, it is hard to interpret dark matter halo mass estimates based on observed velocities.

There are three UDGs with kinematic measurements in the literature. All three seem to have massive dark matter halos ($10^{10}-10^{11} M_{\odot}$, Beasley et al. 2016; van Dokkum et al. 2016, 2017). We analyze here the internal dynamics of three UDGs in the central regions of the Virgo cluster, doubling the current statistics. We assume that the distance to the Virgo cluster is 16.5 Mpc (Mei et al. 2007; Blakeslee et al. 2009).

## 2. Data

### 2.1. Sample Selection

We target GC candidates in the Virgo UDGs VLSB-B, VLSB-D, and VCC615 (Binggeli et al. 1987; Mihos et al. 2015, 2017). The GCs are selected from the Next Generation Virgo cluster Survey (NGVS; Ferrarese et al. 2012). Point-like sources are split into three categories, attending to their probability of being foreground stars, GCs in the Virgo cluster, and background galaxies. These probabilities are obtained combining the position of all point-like sources in different color-color diagrams based on $u^*$, g, i, z photometry (and Ks only available for VLSB-B, Muñoz et al. 2014) with the inverse concentration parameter (ic), which measures how point-like or extended the object is (see Figure 1; Powalka et al. 2016, E. W. Peng et al. 2018, in preparation).

We select objects with $g < 24.5$ and with higher probability of being GCs than being foreground stars or background galaxies. Due to the large field of view of the DEIMOS spectrograph ($16.3 \times 5'$), we also include some foreground stars that, due to their position on the sky, are candidates for Virgo Overdensity and Sagittarius Streams (Figure 1). Their analysis will be presented in a future paper.

### 2.2. Observations and Data Reduction

The observations were carried out with the DEIMOS spectrograph (Faber et al. 2003) located at the Keck II 10 m telescope (Mauna Kea Observatory). We designed one mask per UDG and the 600 lines/mm grating centered at 7200 Å with slit widths of 1" and the GG455 blocking filter. The wavelength coverage is 4700-9200 Å with a pixel scale of 0.52 Å/pixel, and a spectral resolution of 2.8 Å (FWHM).

The three DEIMOS slitmasks had position angles (P.A.) of 105°, -152°, and -140°, respectively, for VLSB-B, VLSB-D, and VCC615. All the slits were aligned with the slitmasks, except for VCC615, for which the slits had a 10° offset, resulting in the slits having P.A. = -130°.

The slitmasks were observed on 2017 March 04, with exposure times of 83 minutes for VLSB-B, 78 minutes for VLSB-D, and 87 minutes for VCC615. The average seeing was 0.6 (FWHM).

We reduced the data with the SPEC2D pipeline (Cooper et al. 2012; Newman et al. 2013) with improvements described by Kirby et al. (2015a, 2015b). The wavelength solution is improved by tracing the sky lines along the slit and improving the extraction of the one-dimensional spectra by accounting for the differential atmospheric refraction along the slit. The main steps in the reduction process consisted of flat-field corrections, wavelength calibration, sky subtraction, and cosmic-ray cleaning. In Figure 2, we show an example of the fully reduced spectra for objects with high, intermediate, and low S/N.

## 3. Results

### 3.1. Radial Velocity Measurements and Membership Criteria

Line-of-sight radial velocities are measured following the same steps as those described in Toloba et al. (2016a). In short, we feed the penalized pixel-fitting software (pPXF; Cappellari & Emsellem 2004) with 17 high signal-to-noise ($100 < S/N < 800 \AA^{-1}$) stellar templates observed with the same instrumental setup as the science data. To reduce the mismatch fitting problem, the stellar templates include stellar types from B1 to M0 and luminosity classes from supergiants to dwarfs. The radial velocity uncertainties are calculated running 1000 Monte Carlo simulations, where the flux of each spectrum is perturbed within the flux uncertainty obtained during the reduction process, assuming that it is Gaussian.

The final radial velocities are corrected by small offcentering effects across the slits. This affects unresolved sources and is quantified using the atmospheric B and A bands at 6850-7020 Å and 7580-7690 Å. The resulting radial velocity uncertainties are the quadrature sum of the uncertainty in the observed radial velocity and the A and B bands.

The membership criteria is described in Toloba et al. (2016a) and summarized in Figure 1. Those GCs that are within a box of $\Delta R/R_e < 10$ and $|\Delta V| < 150 \text{ km s}^{-1}$, approximately three times the typical velocity dispersion of dwarf galaxies, are considered GC satellites. The expected contamination for all the GC candidates combined is $1.1 \pm 0.2$ within this box. The contaminants would be intracluster GCs, GC satellites of other galaxies in Virgo, and Milky Way (MW) stars. Background galaxies are spectroscopically identified for their emission lines and removed from the sample. The GCs classified as satellites have median photometric probabilities of being GCs of 91% and photometric probabilities of being stars smaller than 5%.

### 3.2. Velocity Dispersion and Velocity Gradient

The dynamical properties of the three UDGs are analyzed using the Markov Chain Monte Carlo (MCMC) method (Foreman-Mackey et al. 2013). We make two implementations to avoid having more than two free parameters at a time. Assuming that the line-of-sight radial velocities (v) come from a Gaussian distribution, the logarithmic probability of the observed velocities for a certain systemic velocity ($V_{\text{sys}}$) and velocity dispersion ($\sigma$) is

$$ \mathcal{L}(V_{\text{sys}}, \sigma) = -\frac{1}{2} \sum_{n=1}^{N} \log(2\pi(\sigma^2 + \delta v_n^2)) - \sum_{n=1}^{N} \frac{(v_n - V_{\text{sys}})^2}{2(\sigma^2 + \delta v_n^2)}. \quad (1) $$

where N is the number of GC satellites and $\delta v$ are the radial velocity uncertainties, which contribute to increase the width of the Gaussian distribution.
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Upper panel: g-band image for the central 3 × 2 deg of the Virgo cluster from Mihos et al. (2017). Second row of panels: NGVS g-band images of the three UDGs. Different colors indicate objects of different nature based on their spectrophotometric properties (symbols are the same as in the lowest panels). Note that VLSB-B is not nucleated. Third row of panels: color-color diagrams. Symbols are split into high inverse concentration (ic; i.e., extended) and low ic (point-like) sources. The K band, only available for VLSB-B, clearly separates the GC and stellar locus. The blue ($u-g \sim 0.2$) vertical band is the locus of background galaxies. Lower panels: membership diagrams. The orange triangle in VLSB-B indicates an object whose probability of being a GC based on multi-wavelength photometry, that includes the K band, and the extreme deconvolution technique is just slightly higher than its probability of being an MW star (60% vs. 38%). To be cautious, we will not consider this object as a GC satellite.</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Examples of three GC spectra with different S/Ns put into the rest frame using their radial velocities. The panels are arranged in order of decreasing S/N and luminosity from top to bottom. In each panel, the S/N, heliocentric velocity, and g-band magnitude are shown. Three wavelength regions are shown for each spectrum: the region that includes the Hβ and the Mg triplet lines, the region that includes the Hα line, and the region that includes the Ca triplet lines. These lines are indicated with vertical dashed black lines.</figcaption></figure>

We run the implementation based on Equation (1) twice. The second time $V_{sys}$ is fixed to the heliocentric velocity of the nucleus. This can only be done for VCC615 and VLSB-D, the two nucleated UDGs.

Given the low numbers of GC satellites, we perform simulations to test the statistical significance and possible biases in the calculations. We simulate Gaussian distributions with input dispersions from 10 to 100 km s⁻¹ in steps of 10 km s⁻¹. For each one of these distributions, we randomly select 4, 7, and 12 velocities with uncertainties that are the average uncertainty of our observed radial velocities. For each randomly selected sample, we apply the same MCMC method as described above. In general, the input velocity dispersion is always recovered with possibly a small bias of ~5 km s⁻¹ for large input velocity dispersions: >80 km s⁻¹ for samples of 12 GC satellites; >40 km s⁻¹ for samples of seven GC satellites. For samples of four GC satellites and input dispersions >40 km s⁻¹, the overestimation can be as high as 10 km s⁻¹. However, this small bias is always within the measured error bars. On the contrary, if $V_{sys}$ is fixed and $\sigma$ is the only free parameter, the difference between the input and output dispersion is always $\le 5$ km s⁻¹ for input dispersions $\le 60$ km s⁻¹.

We estimate whether VLSB-D, the UDG with the largest number of spectroscopically confirmed GC satellites, shows internal rotation using the MCMC implementation described in Martin & Jin (2010). The logarithmic probability in this case is

$$ \mathcal{L}(V_{sys}, \sigma, dv/dr, \phi) = -\frac{1}{2} \sum_{n=1}^{N} \log(2\pi(\sigma^2 + \delta v_n^2)) - \sum_{n=1}^{N} \frac{(v_n - V_{sys} - \frac{dv}{dr}r_n)^2}{2(\sigma^2 + \delta v_n^2)}, \quad (2) $$

where $V_{sys}$ is fixed to the value obtained running the MCMC in Equation (1) and $dv/dr$ is the velocity gradient along the projected distance $r$ with P.A. = $\phi$:

$$ r = (\text{R.A.} - \text{R.A.}_0)\cos(\text{decl.}_0)\sin(\phi) + (\text{decl.} - \text{decl.}_0)\cos(\phi). \quad (3) $$

$R.A.0$ and $decl.0$ are the coordinated of the photometric galaxy center.

The first time we run the MCMC following Equation (2), we include $\phi$ as a free parameter. From that analysis, we find the angle that maximizes $dv/dr$ (shown in Figure 3). In the final run, we fix $\phi$ to this suggested position angle.
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Upper and middle panels: two-dimensional and marginalized posterior probability density functions for the systemic velocity ($V_{\text{sys}}$), velocity dispersion ($\sigma$), and velocity gradient ($dv/dr$). The orange lines represent the 16th, 50th, and 84th quartiles from left to right. These are considered the best value and $1\sigma$ uncertainties. The three galaxies have systemic velocities consistent with being in the Virgo cluster and dispersions within the typical values for low luminosity galaxies ($\lesssim 50$ km s$^{-1}$ for $M^* < 10^{9} M_{\odot}$; Geha et al. 2003; Toloba et al. 2014). Lower panels: velocity gradient for VLSB-D.</figcaption></figure>

In both MCMC implementations, we use flat priors within plausible physical ranges: $V_{\text{sys}}$ is within the typical values for Virgo cluster galaxies ($-500 < V_{\text{sys}} < 3000$ km s$^{-1}$); dispersions are within $0 < \sigma < 200$ km s$^{-1}$; and velocity gradients are within $-30 < dv/dr < 30$ km s$^{-1}$ arcmin$^{-1}$.

The upper panels of Figure 3 show the MCMC results for Equation (1). All three UDGs have $V_{\text{sys}}$ consistent with being galaxies in the Virgo cluster and show a wide range of low velocity dispersions ($<50$ km s$^{-1}$). The $V_{\text{sys}}$ and location in the sky suggest that VLSB-B and VCC615 are members of the Virgo subcluster A within $\sim 1.1\sigma$ of its velocity distribution (see Boselli et al. 2014, for a description of the spectrophotometric parameters of Virgo substructures). The $V_{\text{sys}}$ of VLSB-D is smaller than the value measured for its nucleus (see Figure 1 and Table 1). This suggests that the nucleus is not at the center of the gravitational potential, which is also supported by this nucleus being $\sim 1$ kpc spatially off-centered. The $V_{\text{sys}}$ of VCC615 coincides with the velocity and position of its nucleus, which suggests that it is at the center of the gravitational potential. The systemic velocity of VLSB-B is consistent with zero, which makes the available Ks-band photometry for these sources essential (see Figure 1 and Muñoz et al. 2014). Our four VLSB-B GC satellites have probabilities of being GCs >86%, while their probabilities of being MW stars are <6%.

We perform two sets of simulations to investigate the effect that having one MW star in our sample of four GC satellites in VLSB-B would have in our measured velocity dispersion. In the first set of simulations, we randomly select samples of three and four objects within a Gaussian distribution with widths from 10 to 100 km s$^{-1}$ in steps of 10 km s$^{-1}$. For each randomly selected sample, we calculate the velocity dispersion following Equation (1). The velocity dispersions always agree within the error bars, although the uncertainties for calculations done with three objects are 17% larger. In the second set of simulations, we select three objects from a Gaussian distribution with a width of 45 km s$^{-1}$, which represent GC satellites, and one object from a Gaussian distribution with a width of 100 km s$^{-1}$, assuming that the halo of the MW has the same dispersion as that of M31 (Gilbert et al. 2014). We also include that the probability of this object being an MW star is <6%, as obtained from their photometric information. We calculate the velocity dispersion of the four objects following Equation (1). The results of these simulations suggest that we can reject with 90% confidence the hypothesis of having an MW star in our sample. In summary, all these simulations indicate that the probability of having an MW star in our sample of GC satellites is very low but, if it is there, it does not affect the measured velocity dispersion, only increases its uncertainty.

The lower panels of Figure 3 show the measured velocity gradient for VLSB-D. We use our simulations presented in
<table><caption>**Table 1**

Properties of the UDGs</caption><thead><tr><th></th><th>VLSB-B</th><th>VLSB-D</th><th>VCC615</th></tr></thead><tbody><tr><td>R.A. (hh:mm:ss)</td><td>12:28:10.6</td><td>12:24:42.1</td><td>12:23:04.7</td></tr><tr><td>Decl. (dd:mm:ss)</td><td>+12:43:28</td><td>+13:31:02</td><td>12:00:56</td></tr><tr><td>M<sub>V</sub> (mag)</td><td>-13.5 ± 0.2</td><td>-16.2 ± 0.4</td><td>-14.7 ± 0.1</td></tr><tr><td>R<sub>e</sub> (kpc)</td><td>2.9 ± 0.2</td><td>13.4 ± 2.0</td><td>2.4 ± 0.1</td></tr><tr><td>ε</td><td>0.17 ± 0.15</td><td>0.55 ± 0.10</td><td>0.05 ± 0.05</td></tr><tr><td>⟨μ<sub>V</sub>⟩<sub>e</sub> (mag arcsec<sup>-2</sup>)</td><td>27.5 ± 0.1</td><td>27.6 ± 0.2</td><td>25.8 ± 0.1</td></tr><tr><td>M<sup>*</sup> (×10<sup>7</sup> M<sub>☉</sub>)</td><td>0.6 ± 0.1</td><td>7.9 ± 0.1</td><td>2.1 ± 0.1</td></tr><tr><td>R<sub>h</sub> (kpc)</td><td>1.8<sup>+0.8</sup><sub>-0.6</sub></td><td>8.4<sup>+8.7</sup><sub>-2.8</sub></td><td>1.9<sup>+0.7</sup><sub>-0.5</sub></td></tr><tr><td>N<sub>GC,tot</sub></td><td>12<sup>+7</sup><sub>-5</sub></td><td>36<sup>+47</sup><sub>-17</sub></td><td>14<sup>+6</sup><sub>-5</sub></td></tr><tr><td>N<sub>GC,spec</sub></td><td>4</td><td>12</td><td>7</td></tr><tr><td>V (km s<sup>-1</sup>)</td><td>24.9<sup>+22.3</sup><sub>-36.2</sub></td><td>1033.8<sup>+5.9</sup><sub>-5.5</sub></td><td>2094.0<sup>+14.9</sup><sub>-13.0</sub></td></tr><tr><td>V<sub>nuc</sub> (km s<sup>-1</sup>)</td><td>—</td><td>1040.1 ± 1.4</td><td>2094.1 ± 2.7</td></tr><tr><td>σ (km s<sup>-1</sup>)</td><td>47<sup>+53</sup><sub>-29</sub></td><td>16<sup>+6</sup><sub>-4</sub></td><td>32<sup>+17</sup><sub>-10</sub></td></tr><tr><td>dv/dr (km s<sup>-1</sup> arcmin<sup>-1</sup>)</td><td>—</td><td>5.9<sup>+11.7</sup><sub>-11.9</sub></td><td>—</td></tr><tr><td>V<sub>rot</sub> (km s<sup>-1</sup>)</td><td>—</td><td>17.2<sup>+33.9</sup><sub>-34.7</sub></td><td>—</td></tr><tr><td>M<sub>1/2</sub> (×10<sup>9</sup> M<sub>☉</sub>)</td><td>4.9<sup>+11.1</sup><sub>-4.9</sub></td><td>3.2<sup>+2.4</sup><sub>-1.7</sub></td><td>2.5<sup>+2.7</sup><sub>-1.6</sub></td></tr><tr><td>M/L<sub>V</sub> (M<sub>☉</sub>/L<sub>☉</sub>)</td><td>407<sup>+916</sup><sub>-407</sub></td><td>21<sup>+15</sup><sub>-11</sub></td><td>60<sup>+65</sup><sub>-38</sub></td></tr><tr><td>f<sub>DM</sub> (%)</td><td>>99 ± 1</td><td>99 ± 1</td><td>>99 ± 1</td></tr></tbody></table>

**Note.** Rows 1–8: photometric parameters. Rows 9–18: spectroscopic measurements. The central coordinates are in J2000. The magnitudes are in the Vega system. $R_e$ is the stellar half-light radius. $\epsilon$ is the ellipticity. $\langle\mu_V\rangle_e$ is the average surface brightness within the $R_e$. These three parameters are taken from Mihos et al. (2015, 2017). $M^*$ is the total stellar mass. $R_h$ is the radius that contains half the total number of GCs. $N_{GC,tot}$ is the total number of GCs. $N_{GC,spec}$ is the number of spectroscopically confirmed GCs. The remaining candidates were not observed or had too low S/Ns to estimate reliable velocities. $V$ is the heliocentric systemic velocity. $V_{nuc}$ is the heliocentric velocity for the nucleus. $\sigma$ is the velocity dispersion of the GC system (it includes rotation if present). $\frac{dv}{dr}$ is the velocity gradient along the P.A. indicated in Figure 3. $V_{rot}$ is the rotation speed derived from the velocity gradient. $M_{1/2}$ is the dynamical mass within the $R_h$. $M/L_V$ is the dynamical mass-to-light ratio within the $R_h$. $f_{DM}$ is the dark matter fraction within the $R_h$ and its error bar refers to the 99% confidence interval.

Toloba et al. (2016a) to address the reliability of this velocity gradient given the low number statistics. These simulations show that for samples smaller than 10 GCs and velocity uncertainties $\delta v \gtrsim 10$ km s$^{-1}$ (or 15%–30% relative velocity uncertainties for low-mass galaxies with $V_{rot}/\sigma = 0–2$), the velocity gradient measured for a galaxy that is not rotating and for a galaxy rotating with $V_{rot}/\sigma \sim 1$ is undistinguishable. This means that any rotation measured under these conditions can be purely by chance. These conditions are met for VLSB-B and VCC615. However, if the number of GCs is >10, the average velocity uncertainty is $\delta v < 10$ km s$^{-1}$, and the galaxy is rotating with $V_{rot}/\sigma \sim 2$, the recovered dispersion and rotation coincide with the input values within the error bars. This suggests that VLSB-D could be rotating along its major axis; however, due to our sample consisting only of 12 GC satellites, more data are needed to confirm this result.

### 3.3. Total Mass and Dark Matter Content

We derive the total mass of the UDGs using the estimator for dynamically hot systems in equilibrium by Wolf et al. (2010):

$$M_{1/2} = 930 \frac{\sigma^2}{\text{km}^2 \text{s}^{-2}} \frac{R_h}{\text{pc}} M_{\odot}, \quad (4)$$

where $\sigma$ is the velocity dispersion measured from Equation (1) and $R_h$ is the radius that contains half the population of the dynamical tracers. In this case, it is the radius that contains half the number of the GCs (see Table 1). The diffuse nature of the UDGs makes it challenging to decide where the GC population ends, and as a result, $R_h$ is very uncertain and it is usually assumed that $R_h = R_e$ (e.g., Beasley et al. 2016). Using NGVS images, we use MCMC to fit the GC number density profile with a Sérsic function with index $n = 1$ assuming circular GC distribution. We obtain $R_h = 0.85_{-0.21}^{+0.26}$ arcmin for VCC1287, which is $1.12R_e$ (assuming $R_e = 45''5$). Estimating $M_{1/2}$ using $R_h$ results in a slightly larger stellar mass and mass-to-light ratio than that obtained by Beasley et al. (2016), $M_{1/2} = 4.1_{-2.7}^{+4.2} \times 10^9 M_{\odot}$ and $M/L_g = 179_{-117}^{+182}$. For our UDGs, $R_h < R_e$ (see Table 1), although they are consistent within the uncertainties. We use $R_h$ in our calculations, but if we used $R_e$ instead, the derived $M_{1/2}$ would be ~50%–60% larger.

The total masses found for the three UDGs are much higher than the expected values for their stellar masses (see Figure 4). However, their $N_{GC}$ are consistent with the number expected for galaxies with that $M_{1/2}$, although VLSB-B appears to be on the low side of the relation.

We estimate the fraction of dark matter within the $R_h$, assuming that these galaxies do not have gas. We use $g-i$ to estimate the total stellar mass (Taylor et al. 2011) and assume that within the $R_h$ the stellar mass is half, although $R_h < R_e$, which makes the stellar mass within the $R_h$ less than half. The results suggest these UDGs are heavily dark matter dominated (see Table 1).

## 4. Discussion and Conclusions

We spectroscopically confirm 4 GC satellites in VLSB-B, 12 in VLSB-D, and 7 in VCC615. We use them to measure $V_{sys}$ and $\sigma$ of the three UDGs and confirm their dynamical association with the Virgo cluster. We estimate their total $M/L$ within the $R_h$ and find that these galaxies have extremely large values for their stellar mass. Assuming that they follow an NFW profile (Navarro et al. 1997), where the stellar mass is negligible as suggested by their high $M/L$, we find that VLSB-B and VCC615 very likely have dark matter halos of ~$10^{12} M_{\odot}$ (Figure 5). These are typical values for galaxies that have stellar masses two orders of magnitude higher than that of these UDGs.

The interpretation of the dark matter halo of VLSB-D is uncertain given that it may not be in equilibrium. The tidal features, the spatially and dynamically off-center nucleus, and the velocity gradient suggest that VLSB-D is being tidally stripped as it orbits through Virgo. VLSB-D could have recently interacted with M84, given their similar $V_{sys}$ (1017 km s$^{-1}$, Cappellari et al. 2011), and the fact that VLSB-D’s tidal tails align along the direction of M84.

VLSB-B and VCC615 show smooth and round stellar distributions (Mihos et al. 2015). If they are in dynamical equilibrium, these could be within the most dark matter dominated galaxies known, only comparable to other UDGs and Local Group dSphs. However, more GCs should be observed to confirm the estimated $\sigma$. Such high $M/L$ (Figure 4) can only be explained with massive halos and relatively high concentrations, at least for VLSB-B (Figure 5), which might suggest an early collapse and early infall into the cluster
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Relations with the dynamical mass. Left panel: mass-to-light ratio vs. the total mass within the $R_h$. Our data is shown with red symbols. The blue square and diamond are Dragonfly 44 and DX1 from van Dokkum et al. (2017). The blue dot is VCC1287 from Beasley et al. (2016), this value is calculated using the half-light radius of the galaxy instead of $R_h$. VCC1287 has $R_h > R_{h,}$ thus, a more fair comparison would imply moving the blue dot downward. In gray, we show the sample for Local Group dwarf spheroidals (LG dSphs Wolf et al. 2010), Virgo cluster dwarf early types (dEs Toloba et al. 2014), ATLAS$^{3D}$ early-type galaxies (Cappellari et al. 2013), and clusters of galaxies (Zaritsky et al. 2006). Right panel: total number of GCs vs. the total mass within the $R_h$. Comparison sample by Jordán et al. (2009) and Harris et al. (2013).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Circular velocity (left) and total mass (right) estimates at the half-light radius, symbols are the same as in Figure 4. The black lines show velocity and mass profiles corresponding to NFW halos with fiducial masses $M_{200} = 10^{11}$ and $10^{12} M_{\odot}$. We choose the average concentration $c = 8.3, 10$, respectively, following Dutton & Macciò (2014). The shaded region indicates the scatter expected by allowing the concentration to change by 25%. Whereas VLSB-D seems to have a lower dark matter content than expected given its stellar mass, VLSB-B and VCC615 are consistent with a relatively massive MW-like halo.</figcaption></figure>

(Navarro et al. 1997). This scenario suggests that VLSB-B and VCC615 could be “failed” galaxies that formed fewer stars than expected for their likely massive dark matter halos. This could be due to an extremely low star formation efficiency or an abrupt truncation of their star formation due to the early interaction with the hot intracluster medium.

Our data suggests a structurally and dynamically diverse population of galaxies, where round and extremely low surface brightness galaxies could be rapidly rotating. It is important to statistically quantify the significance of such rotation not only in

VLSB-B but also in other UDGs by increasing the number of GCs observed and the number of UDGs studied dynamically.

E.T. acknowledges the support from the Eberhardt Fellowship awarded by the University of the Pacific. E.T. and P.G. acknowledge the NSF grants AST-1010039 and AST-1412504. S.L. and E.W.P. acknowledge support from NSFC grant 11573002. L.V.S. acknowledges support from HST-AR-14583 and the Hellman Foundation. The authors thank the referee for useful suggestions that have helped to improve this manuscript.
ORCID iDs

Elisa Toloba https://orcid.org/0000-0001-6443-5570

Eric Peng https://orcid.org/0000-0002-2073-2781

Puragra Guhathakurta https://orcid.org/0000-0001-8867-4234

J. Christopher Mihos https://orcid.org/0000-0002-7089-8616

Patrick Côté https://orcid.org/0000-0003-1184-8114

Alessandro Boselli https://orcid.org/0000-0002-9795-6433

Laura Ferrarese https://orcid.org/0000-0002-8224-1128

Roberto Muñoz https://orcid.org/0000-0003-1743-0456

Thomas Puzia https://orcid.org/0000-0003-0350-7061

References

Amorisco, N. C., & Loeb, A. 2016, MNRAS, 459, L51

Amorisco, N. C., Monachesi, A., & White, S. D. M. 2018, MNRAS, 475, 4235

Beasley, M. A., Romanowsky, A. J., Pota, V., et al. 2016, ApJL, 819, L20

Binggeli, B., Tammann, G. A., & Sandage, A. 1987, AJ, 94, 251

Blakeslee, J. P., Jordán, A., Mei, S., et al. 2009, ApJ, 694, 556

Boselli, A., Voyer, E., Boissier, S., et al. 2014, A&A, 570, A69

Caldwell, N. 2006, ApJ, 651, 822

Cappellari, M., & Emsellem, E. 2004, PASP, 116, 138

Cappellari, M., Emsellem, E., Krajnović, D., et al. 2011, MNRAS, 413, 813

Cappellari, M., McDermid, R. M., Alatalo, K., et al. 2013, MNRAS, 432, 1862

Cooper, M. C., Newman, J. A., Davis, M., Finkbeiner, D. P., & Gerke, B. F. 2012, spec2d: DEEP2 DEIMOS Spectral Pipeline, Astrophysics Source Code Library, ascl:1203.003

Dalcanton, J. J., Spergel, D. N., Gunn, J. E., Schmidt, M., & Schneider, D. P. 1997, AJ, 114, 635

Di Cintio, A., Brook, C. B., Dutton, A. A., et al. 2017, MNRAS, 466, L1

Dutton, A. A., & Macciò, A. V. 2014, MNRAS, 441, 3359

Faber, S. M., Phillips, A. C., Kibrick, R. I., et al. 2003, Proc. SPIE, 4841, 1657

Ferrarese, L., Côté, P., Cuillandre, J.-C., et al. 2012, ApJS, 200, 4

Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306

Geha, M., Guhathakurta, P., & van der Marel, R. P. 2003, AJ, 126, 1794

Gilbert, K. M., Kalirai, J. S., Guhathakurta, P., et al. 2014, ApJ, 796, 76

Harris, W. E., Harris, G. L. H., & Alessi, M. 2013, ApJ, 772, 82

Impey, C., Bothun, G., & Malin, D. 1988, ApJ, 330, 634

Janssens, S., Abraham, R., Brodie, J., et al. 2017, ApJL, 839, L17

Jordán, A., Peng, E. W., Blakeslee, J. P., et al. 2009, ApJS, 180, 54

Kirby, E. N., Guo, M., Zhang, A. J., et al. 2015a, ApJ, 801, 125

Kirby, E. N., Simon, J. D., & Cohen, J. G. 2015b, ApJ, 810, 56

Koda, J., Yagi, M., Yamanoi, H., & Komiyama, Y. 2015, ApJL, 807, L2

Martin, N. F., & Jin, S. 2010, ApJ, 721, 1333

Martínez-Delgado, D., Läsker, R., Sharina, M., et al. 2016, AJ, 151, 96

Mei, S., Blakeslee, J. P., Côté, P., et al. 2007, ApJ, 655, 144

Mihos, J. C., Durrell, P. R., Ferrarese, L., et al. 2015, ApJL, 809, L21

Mihos, J. C., Harding, P., Feldmeier, J. J., et al. 2017, ApJ, 834, 16

Moore, B., Katz, N., Lake, G., Dressler, A., & Oemler, A. 1996, Natur, 379, 613

Muñoz, R. P., Eigenthaler, P., Puzia, T. H., et al. 2015, ApJL, 813, L15

Muñoz, R. P., Puzia, T. H., Lançon, A., et al. 2014, ApJS, 210, 4

Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493

Newman, J. A., Cooper, M. C., Davis, M., et al. 2013, ApJS, 208, 5

Peng, E. W., Jordán, A., Côté, P., et al. 2006, ApJ, 639, 95

Powalka, M., Lançon, A., Puzia, T. H., et al. 2016, ApJS, 227, 12

Román, J., & Trujillo, I. 2017, MNRAS, 468, 703

Sandage, A., & Binggeli, B. 1984, AJ, 89, 919

Taylor, E. N., Hopkins, A. M., Baldry, I. K., et al. 2011, MNRAS, 418, 1587

Toloba, E., Guhathakurta, P., Peletier, R. F., et al. 2014, ApJS, 215, 17

Toloba, E., Li, B., Guhathakurta, P., et al. 2016a, ApJ, 822, 51

Toloba, E., Sand, D. J., Spekkens, K., et al. 2016b, ApJL, 816, L5

van der Burg, R. F. J., Muzzin, A., & Hoekstra, H. 2016, A&A, 590, A20

van Dokkum, P., Abraham, R., Brodie, J., et al. 2016, ApJL, 828, L6

van Dokkum, P., Abraham, R., Romanowsky, A. J., et al. 2017, ApJL, 844, L11

van Dokkum, P. G., Abraham, R., Merritt, A., et al. 2015, ApJL, 798, L45

Venhola, A., Peletier, R., Laurikainen, E., et al. 2017, A&A, 608, A142

Wolf, J., Martinez, G. D., Bullock, J. S., et al. 2010, MNRAS, 406, 1220

Zaritsky, D., Gonzalez, A. H., & Zabludoff, A. I. 2006, ApJ, 638, 725