Spatially-resolved stellar kinematics of the ultra diffuse galaxy Dragonfly 44. I. Observations, kinematics, and cold dark matter halo fits

TL;DR

This study provides detailed stellar kinematics of the ultra diffuse galaxy Dragonfly 44, revealing a high dark matter content and insights into its dark matter halo profile, challenging some galaxy formation models.

Contribution

First spatially-resolved kinematic analysis of Dragonfly 44, showing a rising velocity dispersion profile and supporting a cored dark matter halo model.

Findings

No evidence of rotation in Dragonfly 44.

High mass-to-light ratio indicating dark matter dominance.

Dark matter halo possibly has a core rather than a cusp.

Abstract

We present spatially-resolved stellar kinematics of the well-studied ultra diffuse galaxy (UDG) Dragonfly 44, as determined from 25.3 hrs of observations with the Keck Cosmic Web Imager. The luminosity-weighted dispersion within the half-light radius is $\sigma_{1/2}=33^{+3}_{-3}$ km/s. There is no evidence for rotation, with $V/\sigma<0.12$ (90% confidence) along the major axis, in apparent conflict with models where UDGs are the high-spin tail of the normal dwarf galaxy distribution. The spatially-averaged line profile is more peaked than a Gaussian, with Gauss-Hermite coefficient $h_4=0.13\pm 0.05$. The mass-to-light ratio within the effective radius is $M/L=26^{+7}_{-6}$, similar to other UDGs and higher by a factor of six than normal galaxies of the same luminosity. This difference between UDGs and other galaxies is, however, sensitive to the aperture that is used, and is much reduced when the $M/L$ ratios are measured within a fixed radius of 10 kpc. Dragonfly 44 has a rising velocity dispersion profile, from $\sigma=26^{+4}_{-4}$ km/s at R=0.2 kpc to $\sigma=41^{+8}_{-8}$ km/s at R=5.1 kpc. The profile can only be fit with a cuspy NFW profile if the orbital distribution has strong tangential anisotropy, with $\beta=-0.8^{+0.4}_{-0.5}$. An alternative explanation is that the dark matter profile has a core: a Di Cintio et al. (2014) density profile with a mass-dependent core provides a very good fit to the kinematics for a halo mass of $\log (M_{200}/{\rm M}_{\odot})=11.2^{+0.6}_{-0.6}$ and $\beta=-0.1^{+0.2}_{-0.3}$, i.e., isotropic orbits. This model predicts a slight positive kurtosis, in qualitative agreement with the measured $h_4$ parameter. UDGs such as Dragonfly 44 are dark matter dominated even in their centers, and can constrain the properties of dark matter in a regime where baryons usually dominate the kinematics: small spatial scales in massive halos.

Tables

Table 1. Table 1 Exposure times

Date	Science	Sky
	[sec]	[sec]
January 22	7,200	2,400
February 11	10,800	4,800
February 12	7,200	4,800
February 13	10,800	4,800
February 17	3,600	2,400
February 18	5,400	2,400
April 13	5,400	2,400
April 17	1,800	2,400
May 10	9,000	3,600

Table 2. Table 2 Velocity Dispersion Profile a a The data in Fig.  14 are σ eff = ( σ 2 + v 2 ) 0.5 subscript 𝜎 eff superscript superscript 𝜎 2 superscript 𝑣 2 0.5 \sigma_{\rm eff}=(\sigma^{2}+v^{2})^{0.5} .

$R$bbLuminosity-weighted average radius of elliptical aperture; $R_{\mathrm{maj}}\approx 1.20R$.	$R$	$v$	$\sigma$
[arcsec]	[kpc]	[km s-1]	[km s-1]
$0\stackrel{\prime \prime }{.}5$	0.23	${0.2}_{-2.6}^{+2.2}$	${26.1}_{-3.5}^{+4.4}$
$1\stackrel{\prime \prime }{.}0$	0.49	$-{3.4}_{-2.8}^{+2.6}$	${26.7}_{-3.4}^{+4.1}$
$1\stackrel{\prime \prime }{.}6$	0.79	$-{3.1}_{-2.5}^{+2.2}$	${26.5}_{-2.9}^{+4.4}$
$2\stackrel{\prime \prime }{.}3$	1.13	$-{0.7}_{-1.8}^{+2.0}$	${31.8}_{-2.9}^{+3.3}$
$3\stackrel{\prime \prime }{.}2$	1.53	${0.5}_{-2.6}^{+2.0}$	${29.1}_{-2.4}^{+3.4}$
$4\stackrel{\prime \prime }{.}0$	1.94	${0.7}_{-2.0}^{+2.0}$	${29.5}_{-2.6}^{+3.0}$
$5\stackrel{\prime \prime }{.}3$	2.55	$-{0.4}_{-2.0}^{+2.3}$	${29.3}_{-2.5}^{+3.1}$
$7\stackrel{\prime \prime }{.}4$	3.62	${0.3}_{-2.5}^{+2.3}$	${34.4}_{-3.6}^{+3.8}$
$10\stackrel{\prime \prime }{.}6$	5.13	${5.9}_{-4.1}^{+3.8}$	${40.2}_{-7.6}^{+7.9}$

Table 3. Table 3 Radial Photometry

$R$	$R$	$V_{606}$	$V_{606}-I_{814}$
[arcsec]	[kpc]	[mag arcsec-2]	[mag]
$0\stackrel{\prime \prime }{.}5$	0.23	$24.75\pm 0.01$	$0.53\pm 0.02$
$1\stackrel{\prime \prime }{.}0$	0.49	$24.77\pm 0.01$	$0.52\pm 0.02$
$1\stackrel{\prime \prime }{.}6$	0.79	$24.84\pm 0.01$	$0.49\pm 0.03$
$2\stackrel{\prime \prime }{.}3$	1.13	$24.93\pm 0.02$	$0.47\pm 0.03$
$3\stackrel{\prime \prime }{.}2$	1.53	$25.09\pm 0.02$	$0.48\pm 0.03$
$4\stackrel{\prime \prime }{.}0$	1.94	$25.27\pm 0.02$	$0.50\pm 0.04$
$5\stackrel{\prime \prime }{.}3$	2.55	$25.54\pm 0.03$	$0.47\pm 0.05$
$7\stackrel{\prime \prime }{.}4$	3.62	$26.07\pm 0.04$	$0.51\pm 0.08$
$10\stackrel{\prime \prime }{.}6$	5.13	$26.67\pm 0.08$	$0.55\pm 0.13$

Equations

$$\ln p=-\frac{1}{2}\sum _{\lambda =4700}^{5500}\left[\frac{{\left(F_{\lambda }-M_{\lambda }\right)}^2}{e_{\lambda }^2}+\ln e_{\lambda }^2\right],$$

$$M(\lambda )={\alpha }_0{S}_{\mathrm{avg}}(\lambda )+\sum _{i=1}^{10}{\alpha }_i{\mathrm{PCA}}_i(\lambda )+{\alpha }_{11}{T}_{\mathrm{gal}}(\lambda ).$$

$$M_{\mathrm{sky}}(\lambda )=M_{\mathrm{full}}(\lambda )-{\alpha }_{11}{T}_{\mathrm{gal}}(\lambda ),$$

$$d=\sum _r{\left[{\mu }_{\mathrm{KCWI}}(r)-(a{\mu }_{\mathrm{HST}}(r)+b)\right]}^2,$$

$$F_{\mathrm{sub}}(\lambda )=F(\lambda )-\left({\alpha }_0{S}_{\mathrm{avg}}(\lambda )+\sum _{i=1}^{10}{\alpha }_i{\mathrm{PCA}}_i(\lambda )+{\alpha }_{11}-0.85\right).$$

$$F_{\mathrm{avg}}(\lambda )=\frac{{\sum }_{i=1}^{34}{n}_i{w}_i{F}_{\mathrm{sub},i}(\lambda )}{{\sum }_{i=1}^{34}{n}_i{w}_i},$$

$$w_{i}=\frac{a_{i}^{2}}{m_{i}}\Big{\langle}\frac{a_{i}^{2}}{m_{i}}\Big{\rangle}^{-1}_{i},$$

$$e(F_{\mathrm{avg}})(\lambda )=\frac{s_{\mathrm{bi}}}{m_{\mathrm{avg}}}\sqrt{M_{\mathrm{avg}}(\lambda )}.$$

$$F_{\mathrm{opt}}(\lambda )=\frac{{\sum }_{i=1}^8{w}_i{F}_{\mathrm{avg},\mathrm{aper}i}(\lambda )}{{\sum }_{i=1}^8{w}_i},$$

$$w_i=\frac{f_{\mathrm{avg},\mathrm{aper}i}}{e_{\mathrm{avg},\mathrm{aper}i}}.$$

$${\sigma }_{\mathrm{instr}}(\lambda )=0.377-5.79\times 10^{-5}{\lambda }_{5000}-1.144\times 10^{-7}{\lambda }_{5000}^2,$$

$$L(y)=\frac{\exp (-y^2/2)}{\sigma \sqrt{2\pi }}\left[1+h_3\frac{y(2y^2-3)}{\sqrt{3}}+h_4\frac{4(y^2-3)y^2+3}{\sqrt{24}}\right],$$

$$L^{\prime }(y)=\{\begin{array}{ll}L(y)& \text{if L(y)>0}\\ 0& \text{otherwise}\end{array}.$$

$${\lambda }_z^{\prime }={\lambda }_z+c_1({\lambda }_z-⟨{\lambda }_z⟩)+c_2({\lambda }_z-⟨{\lambda }_z⟩{)}^2,$$

$$\sigma =(24.9\pm 1.0)+(2.9\pm 0.4)R,$$

$$M(r<r_{1/2})\approx 9.3\times 10^5{\sigma }_{\mathrm{e}}^2{R}_{\mathrm{e},\mathrm{c}},$$

$$\log M/L_I(<r_{1/2})=0.69-0.194L_8+0.0838(L_8{)}^2,$$

$$M(r<10\mathrm{kpc})=\frac{10}{r_{1/2}}M(r<r_{1/2}),$$

$$n=\{\begin{array}{ll}1& \text{if LI<1010L⊙}\\ 1+2.5\left[\log (L_I)-10\right]& \text{otherwise}\end{array}.$$

$$\log M/L_I(<10\mathrm{kpc})=1.52-0.427L_8+0.0682(L_8{)}^2,$$

$$\rho (r)={\rho }_s{\left(\frac{r}{r_s}\right)}^{-\gamma }{\left[1+{\left(\frac{r}{r_s}\right)}^{\alpha }\right]}^{\frac{\gamma -\beta }{\alpha }},$$

$$\gamma =-0.06+\log \left[{\left(10^{X+2.56}\right)}^{-0.68}+10^{X+2.56}\right].$$

$$M_{200}=200{\rho }_{\mathrm{crit}}\frac{4\pi r_{200}^3}{3},$$

$$P(\log (M_{\mathrm{stars}}/L_V))$$

$$P(\log (M_{200}))$$

$$P(-\log (1-\beta ))$$

Full text

††slugcomment: The Astrophysical Journal, in press

1. Introduction

Over the past several years it has been found that large, quiescent galaxies with very low central surface brightness are surprisingly common (van Dokkum et al. 2015; Koda et al. 2015; van der Burg, Muzzin, & Hoekstra 2016). Ultra diffuse galaxies (UDGs), with half-light radii $R_{\rm e}\gtrsim 1.5$ kpc and central surface brightness $\mu(g,0)\gtrsim 24$ mag arcsec*-2*, dominate the population of large galaxies in rich clusters (Danieli & van Dokkum 2019) and have also been found in groups and the general field (Merritt et al. 2016; Martínez-Delgado et al. 2016; Román & Trujillo 2017; van der Burg et al. 2017).

UDGs exhibit a wide variety of properties, as might perhaps be expected given their broad selection criteria: many are smooth and round, resembling very large dwarf spheroidals (van Dokkum et al. 2015), some are clearly tidally-disrupted (such as the spectacular boomerang-shaped galaxy M101-DF4; Merritt et al. 2016), and others are gas rich with widely distributed low level star formation (e.g., Leisman et al. 2017). An intriguing aspects of UDGs is that they often have many globular clusters. The number of clusters varies strongly from galaxy to galaxy, but on average it is 5–7 times higher than in other galaxies of the same luminosity (Beasley et al. 2016; Peng & Lim 2016; van Dokkum et al. 2017a; Amorisco et al. 2018; Lim et al. 2018; Forbes et al. 2018). In at least some UDGs the clusters have similar colors to the smooth galaxy light, and in those UDGs both the clusters and the diffuse light appear to be old, metal poor, and $\alpha-$enhanced, similar to many globular clusters in the Milky Way (e.g., Beasley & Trujillo 2016; Gu et al. 2018; van Dokkum et al. 2018c). Other UDGs appear to be younger, and may have more complex histories (Ferré-Mateu et al. 2018; Fensch et al. 2019; Martín-Navarro et al. 2019).

From a galaxy formation perspective UDGs pose an interesting challenge, as their existence was not explicitly predicted. There are ways to puff up galaxies after their initial formation, for example through external tides (Hayashi et al. 2003; Yozin & Bekki 2015; Ogiya 2018; Jiang et al. 2019; Carleton et al. 2019; Martin et al. 2019) or strong supernova feedback (Agertz & Kravtsov 2016; Di Cintio et al. 2017; Chan et al. 2018). Such “processing” scenarios likely play an important role, although they do not easily account for the high globular cluster numbers (Lim et al. 2018) or the apparent structural integrity (Mowla et al. 2017) of many UDGs in the Coma cluster. Other models seek the origin of UDGs in a combination of low mass, high spin, and late formation (Amorisco & Loeb 2016; Rong et al. 2017, Liao et al. 2019). These models have the benefit of explaining their ubiquity but generically predict that UDGs are young disks, requiring additional processing to turn them into old spheroidal objects (see Liao et al. 2019). In van Dokkum et al. (2018c) we suggested that globular cluster-rich UDGs had extremely high gas densities at the time of their formation, and that feedback from an intense, compact star burst that created the globular clusters caused both the cessation of star formation and the expansion of the galaxies. However, this idea is little more than speculation at this point (see Katz & Ricotti 2013, for related ideas).

Somewhat irrespective of their structural evolution and star formation history, UDGs may provide constraints on the nature and spatial distribution of dark matter. As has long been recognized for dwarf spheroidals (Lin & Faber 1983; Walker et al. 2007; Battaglia et al. 2008; Walker & Peñarrubia 2011) and low surface brightness gas-rich dwarfs and spirals (de Blok et al. 2001; Swaters et al. 2003; Hayashi et al. 2004), galaxies with a low baryon density offer relatively unambiguous information on the dark matter profile. These profiles are often found to be shallower than the cuspy NFW form (Navarro, Frenk, & White 1997). The origin of these shallow profiles (cores) is not well understood; proposed explanations include tidal effects (Read & Gilmore 2005), baryonic processes such as supernova feedback (Governato et al. 2010; Pontzen & Governato 2012; Di Cintio et al. 2014b), warm or mixed dark matter (see Macciò et al. 2012), and “fuzzy” dark matter (e.g., Marsh & Silk 2014). The latter idea postulates that the dark matter is an ultra-light axion with a de Broglie wavelength of hundreds of parsecs.

In this context UDGs such as VCC 1287 (Beasley et al. 2016), Dragonfly 17 (Peng & Lim 2016), and Dragonfly 44 (van Dokkum et al. 2016) occupy an interesting region of parameter space, as their stellar masses are a factor of $\sim 100$ higher than those of dwarf spheroidal galaxies. If they lie on or above the canonical halo mass – stellar mass relation (Moster, Naab, & White 2013) their halos masses are also much higher than those of dwarf spheroidals, and they could help disentangle the processes shaping the distribution of dark matter on kpc scales. An example is the formation of cores, as tidal and baryonic explanations for their presence are most effective at particular mass scales, possibly around $M_{\rm halo}\sim 10^{11}$ M*⊙* (Di Cintio et al. 2014b). Similarly, in fuzzy dark matter models the size of the central “soliton” is expected to scale with halo mass as $r_{\rm sol}\propto M_{\rm halo}^{-1/3}$ (Schive et al. 2014; Wasserman et al. 2019), which means that it is expected to be a more distinct feature in the kinematic profiles of more massive halos (see, e.g., Hui et al. 2017).

So far, only galaxy-integrated measurements of UDG kinematics have been made, and only for a handful of galaxies. They paint a confusing picture, and suggest a remarkable range of dark matter properties (see Spekkens & Karunakaran 2018; Toloba et al. 2018; Alabi et al. 2018; Ferré-Mateu et al. 2018). At one extreme is the large Coma UDG Dragonfly 44, with a stellar velocity dispersion of $\sigma=47^{+8}_{-6}$ km s*-1*,111At least according to van Dokkum et al. (2016) – the true value is almost certainly lower, as shown in this paper. $74\pm 18$ globular clusters, and an estimated halo mass of $M_{200}=10^{11}-10^{12}$ M*⊙* (van Dokkum et al. 2016; Di Cintio et al. 2017). At the other are the galaxies NGC 1052-DF2 and NGC 1052-DF4, with dispersions of $\sigma=8.5^{+2.2}_{-3.1}$ km s*-1* (Danieli et al. 2019) and $\sigma=4.2^{+4.4}_{-2.2}$ km s*-1* (van Dokkum et al. 2019) respectively. These velocity dispersions are consistent with those expected from the stellar mass alone (van Dokkum et al. 2018a; Wasserman et al. 2018a). This large apparent difference in dark matter content is surprising as Dragonfly 44 and the NGC 1052 galaxies have very similar stellar masses and morphologies, and are all rich in globular clusters.

Here we present constraints on the mass profile of a UDG, as derived from spatially-resolved kinematics. This measurement has recently become possible thanks to the arrival of the Keck Cosmic Web Imager (KCWI; Morrissey et al. 2012, 2018) on the Keck II telescope. KCWI is a low surface brightness-optimized integral field unit (IFU) spectrograph, and the combination of its relatively high spectral resolution, low read noise, and blue wavelength coverage make it a near-perfect instrument for UDG spectroscopy. We chose the well-studied UDG Dragonfly 44 for our program (see van Dokkum et al. 2017a). This paper presents the observations, data analysis, kinematics, and dark matter halo fits. Two companion papers discuss constraints on the axion mass in fuzzy dark matter models (Wasserman et al. 2019) and the stellar population of Dragonfly 44 (A. Villaume et al., in preparation). We assume that Dragonfly 44 is at the distance of the Coma cluster, and for convenience we take 100 Mpc for that distance (see Carter et al. 2008). All wavelengths are in air, not vacuum.

2. Observations

2.1. Expected integrated-light spectra for diffuse galaxies

We begin by discussing the expected integrated-light spectra of quiescent low surface brightness galaxies, as this is a relatively new topic (see van Dokkum et al. 2016; Martín-Navarro et al. 2019; Emsellem et al. 2019; Danieli et al. 2019; Chilingarian et al. 2019). Dynamical studies of low luminosity galaxies in the Local Group are typically based on velocity measurements of individual stars. As one example out of many, Geha et al. (2010) obtained the velocities of 520 stars in NGC 147 and 442 stars in NGC 185 to measure the velocity and velocity dispersion profiles of these two Andromeda satellites out to $\sim 8$ effective radii. They used the DEIMOS spectrograph on Keck, which offers excellent multiplexing capability and a resolution of $\sigma_{\rm instr}\approx 20$ km s*-1* in the Ca triplet region. Velocities are routinely determined to an accuracy of 1–2 km s*-1*, enabling measurements of velocity dispersions down to $\lesssim 5$ km s*-1* (e.g., Ibata et al. 2006; Martin et al. 2007; Geha et al. 2010; Collins et al. 2010).

For low mass quiescent galaxies beyond a few Mpc distance only integrated-light measurements can be obtained, and velocity dispersions can only be determined from the broadening of absorption lines of the entire stellar population. This is difficult: the surface brightness is low, which means that the signal-to-noise (S/N) ratio per pixel is low; the metallicity is low, which means that the metal lines are weak; and the velocity dispersion is low, which means that the instrumental resolution needs to be relatively high. The observed velocity dispersion is related to the intrinsic velocity dispersion as $\sigma_{\rm obs}^{2}=\sigma_{\rm stars}^{2}+\sigma_{\rm instr}^{2}$: if the stellar dispersion is, say, 20 % of the instrumental resolution, the observed dispersion is only 2 % larger than the instrumental broadening. If $\sigma_{\rm stars}<\sigma_{\rm instr}$ systematic effects such as template mismatch, small errors in the wavelength calibration, and uncertainties in the (wavelength-dependent) spectral resolution often dominate the error budget (see, e.g., Kelson et al. 2000; Chilingarian et al. 2008; Strader et al. 2009; Emsellem et al. 2019).222When features are not well resolved, fitting codes effectively match the total absorption (the product of the width and depth of the line) rather than (just) the width. This is why template mismatch is a particularly onerous problem when $\sigma_{\rm stars}<\sigma_{\rm instr}$: codes appear to fit velocity widths, but what they are actually fitting is equivalent widths.

Fortunately, if the instrumental resolution is sufficiently high, the low metallicity is compensated to some degree by the low intrinsic velocity dispersion: they conspire to yield absorption features whose observed strength is fairly independent of the surface brightness of the galaxy. This is illustrated in Fig. 1, where we show the expected integrated-light spectra of galaxies with a range of central surface brightness. The spectra are synthetic stellar population synthesis models (Conroy, Gunn, & White 2009), generated at a resolution of $R=10,000$ ($\sigma_{\rm temp}=13$ km s*-1*) using the MIST isochrones (Choi et al. 2016). We assume that the velocity dispersion is related to the surface brightness as $\mu=35-7.5\log\sigma$. This relation is consistent with the central surface brightnesses and velocity dispersions of elliptical galaxies (which have $\mu\sim 17.5$ mag arcsec*-2* and $\sigma\sim 200$ km s*-1*; Franx, Illingworth, & Heckman 1989) and those of dwarf galaxies in the Local Group ($\mu\sim 27$ mag arcsec*-2* and $\sigma\sim 10$ km s*-1*; McConnachie 2012). The relation between velocity dispersion and metallicity is obtained from an approximate fit to the data in Fig. 12 of Gu et al. (2018): [Fe/H] $=1.5\log\sigma-3.15$. The age is assumed to be 10 Gyr for all objects.

Black spectra are at the intrinsic resolution of the galaxies,333The templates were smoothed by a Gaussian of width $\sigma_{\rm sm}^{2}=\sigma_{\rm gal}^{2}-\sigma_{\rm temp}^{2}$. that is, for a hypothetical instrumental resolution of $\sigma_{\rm instr}=0$ km s*-1*. It is clear that even the faintest galaxies with the lowest metallicity have many strong spectral features in the optical, reaching a continuum absorption strength of $\sim 50$ %. However, this is not the case at low spectral resolution. The orange lines show the same spectra for $\sigma_{\rm instr}=100$ km s*-1*. At this instrumental resolution the features become gradually weaker for fainter galaxies, with the Balmer lines and the Ca triplet the only reasonably strong lines redward of $\lambda=4000$ Å for galaxies in the UDG regime ($\mu\gtrsim 24$ mag arcsec*-2*).

The models in Fig. 1 demonstrate qualitatively that it is possible to measure velocity dispersions from metal lines in very low surface brightness galaxies, despite their low metallicity. The fact that the lines are weak is compensated by the fact that they are not blended and have maximum absorption depths that are roughly independent of surface brightness. Most studies of Local Group dwarfs measure stellar velocities from H$\alpha$ and the Ca $\lambda\lambda 8662.1,8542.1,8498.0$ triplet lines, largely because the DEIMOS spectrograph on Keck has its highest sensitivity and spectral resolution in the red. However, at sufficiently high instrumental resolution the integrated-light spectrum blueward of $6000$ Å actually has the highest information content, as is evident in Fig. 1. The strongest features are the Ca $\lambda\lambda 3968.5,3933.7$ H+K lines. These intrinsically-broad lines cannot be used for velocity dispersion measurements but are probably the best features to target for redshift measurements of faint low surface brightness objects.

2.2. Keck Cosmic Web Imager spectroscopy

IFU spectroscopy of Dragonfly 44 was obtained with KCWI on Keck II in the Spring of 2018, following initial observations in 2017 June during commissioning of the instrument. The commissioning data informed the observing strategy in 2018 but are not used in the analysis: conditions were variable, the observing strategy was not yet optimized, and aspects of the instrument and data processing were still being finalized. A list of 2018 dates and exposure times is provided in Table 1. The list does not include nights that had to be discarded due to cirrus or clouds. Maunakea was plagued by bad weather during the entire winter and spring of 2018, and the amount of useable time was about one third of the total allocated time for this program.

The medium slicer was used with the medium resolution BM grating, for a field of view of $16\arcsec\times 20\arcsec$ and an approximate spectral resolution $R\sim 4000$ (see § 4.1 for a measurement of the instrumental resolution as a function of wavelength). The data were taken with $2\times 2$ binning to reduce read noise. A sky position angle of $-32^{\arcdeg}$ was used, as this places the major axis of the galaxy along the long ($20\arcsec$) axis of the IFU. The field of view of KCWI is shown in Fig. 2, along with an HST WFC3/UVIS image of Dragonfly 44 (from van Dokkum et al. 2017a). The central wavelength is $\lambda_{\rm cen}\approx 5050$ Å, with small (10 Å – 20 Å) variations between observing nights so the same wavelengths do not always fall on the same part of the detector. The effective wavelength range is approximately $4650$ Å – 5450 Å.

The observing strategy typically constituted of the following steps. First a nearby star (see the left panel of Fig. 2) was acquired, using the slit viewing camera to center the star in the KCWI field. Then a pre-determined offset was applied to place Dragonfly 44 close to the center of the field. This offset was varied by a few arcseconds for each exposure; this is helpful for diagnosing flat fielding and sky subtraction issues, and yields data over a slightly larger area in the final combined frames than the instantenous KCWI field of view. We then obtained a science exposure of 1,800 s. After the first science exposure we moved the telescope to a relatively empty area about $1\farcm 5$ away, and obtained a 1,200 s “sky” exposure. These are crucial for accurate sky subtraction, as Dragonfly 44 overfills the KCWI field of view (see Fig. 2). We then moved back to Dragonfly 44 and obtained another 1,800 s science exposure. A typical nightly sequence was science – sky – science – science – sky – science – science – sky, but this varied somewhat during the runs due to changing conditions, telescope/instrument problems, and other issues.

The total exposure time of frames that went into our final stack is $61,200$ s, or 17 hrs. The total exposure time that went into our sky analysis is $30,000$ s, and the total science + sky time that is used in the analysis is 25.3 hrs. In addition to these science and blank field data we obtained standard sets of daytime darks, flat fields, and arc lamp exposures.

3. Data reduction

3.1. Pipeline processing

The KCWI Data Extraction and Reduction Pipeline (KDERP) is maintained in a public github repository.444https://github.com/Keck-DataReductionPipelines/KcwiDRP We used this pipeline, with default settings, to turn individual science and sky frames into wavelength-calibrated, flat-fielded, and cosmic ray-cleaned data cubes. What follows is a brief summary of the pipeline processing steps; we refer to § 4 of Morrissey et al. (2018) and the documentation in the github repository for more detailed information.

The pipeline is modular, with eight stages (nine when including a “book keeping” preparation step). In the first stage bias and overscan are subtracted, the data are converted to electrons, and cosmic rays are removed using the L.A.Cosmic (van Dokkum 2001) algorithm. The second stage subtracts the dark and removes scattered light. In the third stage analytic functions are found that describe both geometric distortions and the wavelength calibration. These are based on cross-correlations of the data with arc lamp spectra and a pattern of continuum bars. The output of this stage includes a map that provides the 3D data cube position for each pixel in the 2D image. As discussed in § 4.1, errors in the wavelength calibration are $\approx 0.08$ Å, corresponding to $\approx 5$ km s*-1*. The fourth stage applies a pixel-to-pixel flat field correction, as well as a correction for vignetting and the overall illumination pattern. Stage five is a sky subtraction step, which in our analysis is carried out at a later stage. In stage six the data cubes are generated, based on the functions that were derived in stage three. In stage seven the data are corrected for differential atmospheric refraction. Stage eight is flux calibration, using a standard star; this stage is skipped in our analysis as all our measurements are insensitive to the overall continuum calibration.

The pipeline products that are used in the subsequent steps are the “ocubed” files: the rectified, but un-skysubtracted cubes. The data are sampled on a three-dimensional grid with pixel size $0\farcs 68\times 0\farcs 29\times 0.5$ Å.

3.2. Correction for residual spatial variation

Inspection of the wavelength-collapsed 2D science and sky frames shows a small gradient in the background at a level of $\approx 3$ %. Within each night the gradient has a similar amplitude (as a fraction of the background) in the sky frames and the science frames, and as these have different exposure times we conclude that it is likely a multiplicative rather than an additive effect. For each night a correction flat is created by fitting low order 2D surfaces to all the wavelength-collapsed sky frames and then averaging these fits. The fits are done iteratively with aggressive outlier rejection so that serendipitous objects in the sky frames are ignored. The science frames are then divided by the correction flat. Although the amplitude varies somewhat from night to night the correction flats always show the same pattern, a negative gradient in the $x$-direction (the “short”, $16\arcsec$, axis, which corresponds to systematic slice-to-slice variation over the whole detector rather than within slices). We tested that the gradient is not driven by a particular wavelength region, and that treating the variation as an additive rather than a multiplicative effect does not change the results.

3.3. Alignment

Spatially aligning the data cubes is not straightforward, as there is no compact object that is bright enough to determine an accurate position for every exposure. Object C1 (Fig. 2) is usually near the edge of the field, and in about half the exposures outside of it. Instead of using a single object we fit each collapsed data cube to a 2D model of the flux in the entire KCWI field of view. This model is created from the $V_{606}$ HST WFC3/UVIS image of Dragonfly 44, shown in the left panel of Fig. 2. This image is convolved by a Gaussian with a FWHM of $1\farcs 0$ and projected onto the same spatial grid (with a position angle of $-32^{\arcdeg}$ and $0\farcs 68\times 0\farcs 29$ pixels) as the KCWI data. The resulting model for the KCWI spatial flux distribution is shown in the second panel of Fig. 2.

Before performing the fit an approximate background is removed by subtracting the average of the flux in the outer 1-pixel wide perimeter of the collapsed science frame. Both the collapsed frame and the model are normalized so the total flux in each is 1. The best fitting shift with respect to the model is found by a simple grid search, subtracting the model from the shifted science frame at each step and minimizing the square of the residuals. All science exposures yield clear minima and stable solutions. The data cubes are shifted to the common reference frame of the model using linear interpolation.

The final spatial resolution of the combined datacube is a reflection of the seeing during the observations, guiding errors, and the uncertainties in the alignment of individual exposures. We assess the spatial resolution using object C1 in the summed, wavelength-collapsed data cube. After removing the flux from the galaxy by fitting linear functions in $x$ and $y$, we fit one-dimensional Gaussian profiles in both directions. The FWHM in the $x$-direction is 1.98 pixels, corresponding to $1\farcs 3$. The FWHM in the better-sampled $y$ direction is 3.96 pixels, or $1\farcs 1$. At the distance of Coma these values correspond to 0.63 kpc and 0.53 kpc respectively.

3.4. Sky subtraction

The sky subtraction is the most critical step in the data reduction. The median galaxy signal ranges from $4.5$ e*-* pix*-1* in the center to 0.8 e*-* pix*-1* in the outer annulus, whereas the sky continuum is $\approx 40$ e*-* pix*-1*. Typical metal line absorption depths of 10 % therefore correspond to 0.002 $\times$ the sky brightness in the outer annulus, which means that both emission and absorption features in the sky spectrum need to be modeled and subtracted with great precision.

3.4.1 Principal component analysis

A particular complication in our dataset is that Dragonfly 44 is larger than the KCWI field of view, which means that we cannot determine the sky spectrum from empty areas in the science cubes. We opted not to use the nod-and-shuffle technique, as this comes with severe penalties of a factor of two in S/N ratio (a factor of $\sqrt{2}$ due to the fact that only 50 % of the exposure time is spent on-target, and another factor of $\sqrt{2}$ because the noise from the offset field is added to that of the science field) and a factor of four in spectral coverage. Instead, as explained in § 2.2 we interspersed the 1800 s science exposures with 1200 s blank sky frames. We cannot use an average of all the spatially-collapsed blank sky spectra in a particular night to subtract the background from the science cubes of that same night, as the sky continuum, emission, and absorption all vary too much within and between exposures. Instead, we use all the spatially-collapsed individual sky spectra over multiple nights to capture the wavelength-dependent time variation in the sky spectrum. This variation is parameterized with a principal component analysis (PCA), and the PCA components are then fitted to the background in each of the science cubes. We note that this PCA analysis is different from that used in the “Zurich Atmosphere Purge” (ZAP; Soto et al. 2016) algorithm; ZAP captures any spatial variation in the sky spectrum caused by flat fielding errors or other systematic effects, in regions away from objects of interest. In what follows we describe our methodology in more detail.

We split the data in three overlapping sets, as we find that the sky variation in winter is somewhat different than in spring. This could be a seasonal effect, but it is perhaps more likely that it is due to the fact that the time of night when Dragonfly 44 is accessible changes during the Coma season (the end of the night in January and the start of the night in May). The first set consists of 14 sky exposures and the second and third consist of 12 each. A 1D spectrum is created from each sky cube by collapsing both spatial axes, after carefully masking all objects in the field. Next, a PCA is performed on the 1D spectra within each set, using the scikit-learn python implementation of singular value decomposition.555https://scikit-learn.org/ Ten components are used; we verified that the results are nearly identical for eight, nine, or eleven.

The top left panel of Fig. 3 shows the average sky spectrum from the first set of 14 collapsed blank sky cubes.666Although they differ in the details the other two sets produce very similar-looking results. The spectrum is complex, with many absorption and emission features reaching $\gtrsim 10$ % of the continuum. The extremely strong [O i] $\lambda 5577.3$ line was interpolated over, as it falls outside of the wavelength range of interest (indicated by the vertical dotted lines) and would otherwise dominate many of the PCA components. The ten eigenspectra are labeled PC1 – PC10. They disentangle several distinct causes of the variation in the night sky. PC1 is an excellent match to the solar spectrum, as will be demonstrated in § 4.1. The contribution of the Sun is of course higher when data are taken closer to morning and evening twilight and also when the moon is above the horizon.

PC2 mostly reflects the variation in lines from OH radicals; these are produced in the upper atmosphere from recombination of atomic oxygen and are strongest near dawn and dusk. The presence of these lines may seem surprising, as the OH “forest” is usually considered to only be present at wavelengths $\lambda>6100$ Å. At bluer wavelengths the lines are weaker as they have small transition probabilities, but as shown in Fig. 3 the 8–1 and 9–2 Meinel bands (see Osterbrock et al. 1996, 2000) are important contributors to the sky variation. Other important varying lines are H$\beta$ (from geocoronal atomic hydrogen; Burrage, Yee, & Abreu 1989), the [N i] $\lambda\lambda 5198.2,5200.5$ doublet (Sharpee et al. 2005), and Hg i $\lambda$5460.7 (Osterbrock et al. 2000). These lines do not change in lockstep, and the PC3 – PC10 eigenspectra mostly capture different combinations of positive and negative variations of the OH bands and the individual lines.

3.4.2 Fitting and subtracting the sky in the science cubes

The sky in each of the 34 science cubes is fitted with a linear combination of templates. For each science cube, an average “sky + galaxy” spectrum was created by averaging all pixels after masking most of the light of Dragonfly 44 and other objects in the field of view (see § 3.5). Besides sky emission this spectrum also contains flux from the galaxy, as Dragonfly 44 extends beyond the KCWI field of view. In order to model this spectrum we maximize the likelihood

[TABLE]

with $F$ the extracted sky + galaxy spectrum, $e$ the errors in the data, and $M$ a model of the form

[TABLE]

The model is a linear combination of the twelve templates shown in Fig. 3: the average sky spectrum, the ten eigenspectra describing the variation in the sky, and a template for the contribution of Dragonfly 44. This galaxy template, $T_{\rm gal}$, is created by redshifting a 10 Gyr, [Fe/H] = $-1.0$ stellar population synthesis model (see § 2.1) to the velocity of Dragonfly 44, convolving it to a resolution of 40 km s*-1*, and multiplying it by a low order polynomial so that the continuum roughly corresponds to that of the observed galaxy spectrum. This last step accounts for the fact that no flux calibration was done in the analysis. These details do not influence the fit very much, as the main function of $T_{\rm gal}$ is to allow for an additive component that is not part of the sky spectrum. As we show later Dragonfly 44 is actually not the only contribution of this kind; there is an additional unidentified additive component in the science data which is “absorbed” by $\alpha_{11}$, the coefficient for the galaxy template. In this context it is reassuring that we find nearly indistinguishable results for the best-fitting sky models if $T_{\rm gal}$ is replaced by a featureless spectrum.

The fit is performed with the emcee Markov chain Monte Carlo (MCMC) sampler (Foreman-Mackey et al. 2013), finding the best fitting twelve free parameters $\alpha_{0},...,\alpha_{11}$ for each of the 34 science cubes. One hundred walkers are used and 1000 samples are generated. Burn-in is assumed to occur after 800 samples. Broad, uniform priors are used. The walkers always converge quickly and produce well-defined best fit values for all coefficients. The uncertainties in $\alpha_{1},...,\alpha_{10}$ are uncorrelated, as expected. The fits are generally excellent, with residuals consistent with the expected photon noise. An example is shown in Fig. 4, for a science exposure from the night of February 13 2018. The full sky + galaxy model, as given by Eq. 2, is an excellent fit to the data. For clarity only the wavelength region around the redshifted H$\beta$ line is shown; the fit is of similar quality elsewhere.

For each of the 34 science exposures the best-fitting sky model for that exposure is subtracted from each spatial pixel. The sky model is given by

[TABLE]

that is, all the sky components of the full model but not its galaxy template. The residual $F-M_{\rm sky}$ is shown in the bottom panel of Fig. 4, along with the galaxy template. The H$\beta$ line is detected in the outer regions of the galaxy in this individual 1800 s exposure, illustrating the power of KCWI for studying faint, spatially-extended emission.

3.4.3 Residual offset

Our sky subtraction methodology is insensitive to background signal that shows no variation with wavelength, and also to signal that is not present in the blank sky exposures but only in the science cubes. There is evidence for such signal, as the value of $\alpha_{11}$, which accounts for all contributions to the science data that are not accounted for by the sky templates, varies between $\alpha_{11}=0.6$ and $\alpha_{11}=3.1$. Inspection of the H$\beta$ absorption line in the sky-subtracted spectra (see Fig. 4) indicates that it is not the galaxy flux that varies,777Furthermore, it would be difficult to come up with an explanation why the galaxy flux would vary by a factor of $\sim 5$, as the 34 science exposures that were retained were taken under (nearly) photometric conditions. as the absolute absorption (in e*-* pix*-1*) is independent of the value of $\alpha_{11}$.

The remaining offsets in the sky-subtracted science cubes are measured in the following way. We measure the average flux per pixel in elliptical apertures within the wavelength-collapsed data cubes (the third panel from left in Fig. 2), carefully masking contaminating objects. For a perfect sky subtraction this measurement should correspond to the surface brightness profile of Dragonfly 44. We compare these profiles to the actual surface brightness profile, determined from the (degraded) HST $V_{606}$ image of Dragonfly 44 (second panel in Fig. 2). We minimize the difference

[TABLE]

with $a$ and $b$ free parameters. The fit is done for radii $r>3\arcsec$ so that variations in the seeing and centering do not influence the results. The values of $b$ should correspond to the remaining background levels in the 34 sky-subtracted science frames. The process is illustrated for one exposure (the same one as shown in Fig. 4) in the inset of Fig. 5. For $b=0$ the scaled HST surface brightness profile (dotted line) is not a good fit to the profile measured from the collapsed data cube (points). The best fit is obtained for $b=0.52$ (solid line).

In the main panel of Fig. 5 the offset that is derived from the spectral fit ($\alpha_{11}$) is compared to the offset as derived from the surface brightness profile fit ($b$). There is excellent agreement with an offset, demonstrating that $\alpha_{11}$ indeed represents both residual background and galaxy flux. We determine the contribution from the galaxy to $\alpha_{11}$ by calculating the intersection of the relation between $\alpha_{11}$ and $b$ (red line) with the line $b=0$ (dashed line). We find $\alpha_{11}(0)=0.85\pm 0.03$, that is, the galaxy flux within the aperture that is used for the sky background fit is 0.85 e*-1* pix*-1* and the residual offset that needs to be subtracted from the science data cubes is $\alpha_{11}-0.85$. We note that this process for determining the absolute background level ignores a possible color gradient between the effective wavelength of the $V_{606}$ filter ($0.59\,\mu$m) and the central wavelength of the spectra ($0.51\,\mu$m). As shown in Appendix A, the $V_{606}-I_{814}$ color of Dragonfly 44 is constant within the measurement errors at $r>3\arcsec$.

Taking all the results from this section together, the final, sky- and background-subtracted data cubes are given by

[TABLE]

The values of $\alpha_{0}-\alpha_{11}$ are different for each of the 34 cubes, but we apply the same sky and background correction for every spatial pixel. We find no evidence for systematic variations within the KCWI field of view but this cannot be ruled out.

3.5. Extraction and combination of spectra

3.5.1 Definition of apertures

One-dimensional spectra are extracted from each of the 34 sky- and background-subtracted science cubes, using a variety of spatial apertures. The apertures that are used in the kinematic modeling are eight elliptical annuli that follow the isophotes of the galaxy (see Fig. 6). An approximate square root spacing is used; given the surface brightness profile of the galaxy this yields an approximately equal S/N ratio for each of the extracted spectra. The radii, measured along the major axis of the ellipse, are listed in Table 2. The axis ratio of the ellipses is 0.69, and the radius along the major axis can be converted to the circularized radius through $R_{\rm circ}=0.83R_{\rm maj}$.

We also define apertures along the major and minor axis, as well as masks for contaminating objects and special apertures for particular objects in the field. For each aperture and each science exposure a spectrum is obtained by summing the spectra of all unmasked pixels and dividing the summed spectrum by the number of unmasked pixels. The number of unmasked pixels in the outer apertures is not the same for each exposure, due to the spatial dithering.

3.5.2 Combination of spectra

Before combining the 34 individual spectra for each aperture the barycentric velocity correction needs to be applied, as the data were taken over a long time period. The correction ranges from $22$ km s*-1* for the January data to $-18$ km s*-1* in April, which means the variation is of the same order as the instrumental resolution and the central velocity dispersion of the galaxy. To minimize interpolation-induced smoothing the spectra are resampled onto a $2\times$ finer grid with d$\lambda=0.25$ Å in this step. In this resampling step a small wavelength calibration correction is applied, as derived from the fit of the Solar spectrum to PC1 (see § 4.1). Omitting this correction does not lead to discernable changes in the results.

The individual shifted spectra are then combined:

[TABLE]

with $n_{i}$ the number of unmasked pixels of exposure $i$ and $w_{i}$ the weight. These weights are defined as

[TABLE]

with $a$ the galaxy scaling parameter from Eq. 4, $m$ the wavelength-averaged flux of model $M$ (Eq. 2), and angular brackets indicating the mean. The weights vary between 0.79 and 1.29 with an rms of 0.12.

The uncertainties in the averaged spectra are determined from a combination of the averaged sky spectrum and the residuals from a fit to a stellar population synthesis model:

[TABLE]

Here $M_{\rm avg}$ is calculated in an analogous way to $F_{\rm avg}$, $m_{\rm avg}$ is the wavelength-averaged mean of $M_{\rm avg}$, and $s_{\rm bi}$ is the biweight scatter (Beers, Flynn, & Gebhardt 1990) in the residuals from a fit to a stellar population synthesis model (see § 4). This approach ensures that the fits in § 4 have acceptable $\chi^{2}$ and that the ensuing uncertainties in the fit parameters are properly normalized. Finally, we mask pixels with residuals that exceed $2.5\times$ the expected error. These are invariably associated with the strongest sky lines, which are imperfectly modeled with the PCA analysis. In all apertures the fraction of masked pixels in the wavelength range $4850$ Å $<\lambda<5450$ Å is 2–3 %.

3.6. Optimal extraction

In addition to the aperture spectra described above we create a combined spectrum of the entire galaxy that maximizes the S/N ratio:

[TABLE]

with the sums over the 8 elliptical apertures and the weight given by

[TABLE]

Here $f_{{\rm avg,aper}\,\,i}$ is the wavelength-averaged signal and $e_{{\rm avg,aper}\,\,i}$ is the wavelength-averaged error in aperture $i$.

This spectrum provides the best constraints on the average stellar population of Dragonfly 44 (A. Villaume et al., in preparation). In the present study it is used to select the template that is fit to the data in our kinematic modeling (§ 4.2) and to constrain the kinematic line profile (§ 5.3). It is shown in Fig. 7, along with the best-fitting model from § 4 for reference. The median S/N ratio is 48 per 0.25 Å pixel, equivalent to 96 Å*-1*.

4. Measurements of Kinematics

4.1. Modeling of the spectral resolution

As discussed below we use a stellar population synthesis model as a template to measure the kinematics of Dragonfly 44. This requires that the template and the data have the same resolution. We therefore need to accurately characterize the (wavelength-dependent) resolution that is delivered by the instrument.

4.1.1 Fitting the solar spectrum

The instrumental resolution is typically determined from the widths of emission lines in arc lamp exposures; however, both the light path and the data handling of the calibration lamps are different from the science data. The data reduction process of the science data involves the combination of long exposures over many nights and this is likely to impact the effective spectral resolution. Furthermore, the kinematics are measured from template fits to absorption line spectra rather than from the fits of Gaussians to individual emission lines.

Ideally, the instrumental resolution is measured directly from the science data, for example by using higher resolution observations of the same objects as templates (see van Dokkum et al. 2017b). In our case, we make use of the fact that one of the eigenspectra of the sky variation (PC1; see Fig. 3) comprises scattered and reflected sunlight. We fit PC1 with a high resolution solar spectrum obtained from the BAse de données Solaire Sol (BASS2000888http://bass2000.obspm.fr/solar_spect.php), in small wavelength intervals. Both the model and the data were divided by a polynomial of order $(\Delta\lambda/100)+1$. Free parameters in the fit are the radial velocity, the velocity dispersion, and an additive constant. The instrumental line profile is held fixed, using $h_{3}^{+}=-0.005$ and $h_{4}^{+}=-0.094$ (see below). The fit is done using the emcee-based code described in van Dokkum et al. (2016).

The best fits are shown in Fig. 8. The correspondence between PC1 and the solar spectrum is remarkably good, as illustrated by the insets. The resulting instrumental resolution is shown by the solid points in the top left panel of Fig. 9. The line is the best-fitting relation, of the form

[TABLE]

with $\lambda_{5000}=\lambda-5000$ and $\sigma_{\rm instr}$ the second moment of the instrumental resolution in units of Å. This corresponds to $\sigma_{\rm instr}=24.0$ km s*-1* at $\lambda=4800$ Å and $\sigma_{\rm instr}=18.6$ km s*-1* at $\lambda=5400$ Å, and $R\approx 5,600$ at $\lambda=5000$ Å. Forcing the instrumental profile to be Gaussian ($h_{3}^{+}=h_{4}^{+}=0$) produces very similar results, as shown by the open symbols and dashed line.

The relative velocity shift with wavelength, expressed in Å, is shown in the bottom left panel of Fig. 8. We find a small but systematic wavelength calibration error, with peak-to-peak variation of $\pm 0.1$ Å ($\pm 6$ km s*-1*). The rms variation is $0.06$ Å (4 km s*-1*). A fourth-order polynomial fit to this variation (indicated by the solid line) is included in the resampling of the science data (see § 3.5.2).

4.1.2 Instrumental line profile

A major uncertainty in measuring mass profiles from kinematic data is the degree of anisotropy in the velocity distribution, and this can, in principle, be constrained by deviations from a Gaussian profile: flat-topped profiles indicate tangential anisotropy, peaked profiles radial anisotropy (see, e.g., Bender, Saglia, & Gerhard 1994; Thomas et al. 2007; Amorisco & Evans 2012). However, this relies on an excellent characterization of the instrumental line profile, as well as a very high S/N ratio and adequate control of systematics such as the wavelength calibration.

To characterize the instrumental line profile we follow common practice and parameterize deviations from a Gaussian profile with the (asymmetric) $h_{3}$ and (symmetric) $h_{4}$ components of a Gauss-Hermite expansion (van der Marel & Franx 1993; Cappellari et al. 2007). The line of sight velocity distribution is then parameterized by

[TABLE]

with $y=(v-V)/\sigma$ (see, e.g., Cappellari 2017). As the line profile cannot be negative, we impose the modification

[TABLE]

We use $h_{3}^{+}$, $h_{4}^{+}$ to specify the components of this modified profile, and include these in the emcee-based dispersion fitting code.

We find that the line profile is flat-topped, with negative $h_{4}^{+}$, as expected for the slitwidth-limited resolution provided by the medium slicer (see, e.g., Casini & de Wijn 2014). There is no clear trend with wavelength, and the average values for the Gauss-Hermite components in three wavelength regions are $\langle h_{3}^{+}\rangle=-0.005\pm 0.005$ and $\langle h_{4}^{+}\rangle=-0.094\pm 0.014$. The corresponding profile is shown in the right panel of Fig. 9. We use this line profile in the template construction in § 4.2.

4.2. Template construction

The kinematics are measured by fitting a template spectrum to the observed spectra. The template is one of a set of synthetic stellar population synthesis models that have a native resolution of $R=10,000$ and are based on the same set of libraries as discussed in § 2.1. These model spectra are convolved to the resolution of the Dragonfly 44 spectra. This convolution takes the native resolution of the templates into account, as well as the line profile and the wavelength dependence of the KCWI resolution (as determined in § 4.1). We note that the line profile of the synthetic templates is (exactly) Gaussian.

The best-fitting template is determined by fitting models with a range of discrete ages and metallicities to the optimally-extracted spectrum shown in Fig. 7 and determining the relative likelihood. The results are shown in Fig. 10, with the grey level indicating the relative log likelihood. The distribution of the grey points follows the well-known age-metallicity degeneracy, and the best fit is obtained for an age of 10 Gyr and a metallicity [Fe/H]$=-1.25$. This result is in good agreement with the previous measurement by Gu et al. (2018), who found an age of $8.9^{+4.3}_{-3.3}$ Gyr and [Fe/H]$=-1.3^{+0.4}_{-0.4}$ for Dragonfly 44 from deep MaNGA spectroscopy. Quantitative constraints on the (spatially-resolved) stellar population of Dragonfly 44, derived from more flexible lower resolution model fits to the KCWI data, will be presented in A. Villaume et al., in preparation.

The high resolution models we use have solar abundance ratios, whereas stellar populations with ages of $\sim 10$ Gyr are often $\alpha$-enhanced. We model an enhanced Mg abundance by artificially increasing the depth of the Mg triplet lines, using high resolution integrated KCWI spectra of the old and metal poor Milky Way globular clusters M3 and M13. These data were obtained as templates to measure the velocity dispersion of the UDG NGC 1052-DF2 (see Danieli et al. 2019). Here they are used to slightly increase (by $\approx 5$ %) the depth of the Mg lines in the template spectrum, matching them to the observed lines in M3 and M13. We tested that this enhancement has a negligible effect on the final dispersions. The final template is, then, a synthetic stellar population synthesis model with an age of 10 Gyr, [Fe/H]$=-1.25$, and [Mg/Fe]$\approx 0.3$, convolved to match the wavelength-dependent non-Gaussian line profile of KCWI.

4.3. Velocity and velocity dispersion measurements

The kinematic fits are done with the MCMC methodology described briefly above and more extensively in van Dokkum et al. (2016). The template and data were continuum filtered with a polynomial of order $\Delta\lambda/100+1$. The template normalization, velocity, $\sigma$, $h_{3}$, and $h_{4}$ are fit parameters. We also include two optional parameters $c_{1}$, $c_{2}$ to describe any remaining wavelength calibration errors. The template is sampled onto a wavelength grid defined as

[TABLE]

with $\lambda_{z}=\lambda\times(1+z)$. For our default measurements $c_{1}=c_{2}=0$. As discussed above, in the following “velocity dispersions” are actually second moments of the velocity distribution $L(y)$. We use 100 walkers and 1000 samples, with burn-in assumed to occur after 800. The fits are well-behaved, and provide stable and converged minima. We fit all of the apertures described in § 3.5.

In Fig. 11 we show the best fits for the nine elliptical apertures from Fig. 6, with a focus on the region around the H$\beta$ line and the Mg triplet. Qualitatively, the fits are excellent, with the red model generally within the grey band around the data. This is quantified by calculating the residuals from the model fits and comparing these to the expected errors. For all apertures we find that the biweight scatter in the residuals corresponds to the median error within $10$ %, as expected.

We assess the importance of systematic errors by varying the analysis. Neglecting to correct the Mg strength to match globular clusters leads to a change in derived dispersions of $\sim 1$ km s*-1*. Allowing an additive offset or linear combinations of multiple templates leads to $\sim 5$ km s*-1* changes in the derived dispersions. We note, however, that these should not be considered free parameters; the surface brightness profile of the HST image of Dragonfly 44 sets the overall background level, and the template corresponds to the best stellar population synthesis fit (with small errors) to a smoothed version of the data. Changing the wavelength region of the fit or forcing $h_{3}=h_{4}=0$ has a 1–2 km s*-1effect. Splitting the spectra in three equal-length wavelength regions also leads to 1–2 km s-1* effects. Not applying the wavelength calibration correction function (the fit in the bottom left panel of Fig. 9), or fitting for $c_{1}$ and $c_{2}$, has a negligible effect on the dispersions and also on the $h_{3}$ and $h_{4}$ parameters. Allowing $c_{1}$ and $c_{2}$ to be free does have an effect on the derived velocities; these show larger scatter with larger uncertainties, as expected.

5. Kinematics of Dragonfly 44

5.1. Major and minor axis kinematics

We first consider whether the galaxy is supported by rotation, as might be expected in some UDG formation models (e.g., Amorisco & Loeb 2016). The kinematics along the major and minor axis are shown in Fig. 12. They are measured in wedges that grow from $2\arcsec$ width in the center to $\approx 10\arcsec$ at the largest distances from the center. The rotation velocity, with respect to the mean, is shown in the top panels, and the velocity dispersion is shown in the bottom panels. The major and minor axis profiles are shown separately. There is no evidence for rotation. We determine the maximum rotation speed by fitting the normalization of a model rotation curve to the velocity data. The model is the best fitting Jeans model to the rotationally-supported dE galaxy NGC 147, as derived by Geha et al. (2010). NGC 147 is a satellite of M31 with a similar stellar mass as Dragonfly 44. This is an ad hoc way of generating a plausible rotation curve shape; our results are not sensitive to the precise form of the model.

The best-fitting maximum rotation velocities are $V_{\rm max}=1$ km s*-1* and $V_{\rm max}=3$ km s*-1* for the major and minor axis. Both values are consistent with zero, and the 90 % upper limits are $V_{\rm max}<3.4$ km s*-1* and $V_{\rm max}<5.5$ km s*-1* respectively. The mean velocity dispersion is 27 km s*-1* for the major axis and 32 km s*-1* for the minor axis, which implies $V_{\rm max}/\langle\sigma\rangle<0.12$ (major axis) and $V_{\rm max}/\langle\sigma\rangle<0.17$ (minor axis), with 90 % confidence. We note that a physically better-motivated measure of rotational support is $\left(\langle V^{2}\rangle/\langle\sigma^{2}\rangle\right)^{0.5}$, that is, the rms of the locally-measured ratios (see Binney 2005; Cappellari et al. 2007). However, the S/N ratio in individual spatial bins is not sufficiently high to measure this quantity reliably from our data.

The limit on $V_{\rm max}/\langle\sigma\rangle$ along the major axis of Dragonfly 44, combined with its axis ratio of $b/a=0.69$, means that the galaxy is not rotationally-supported. In Fig. 13 Dragonfly 44 is placed on the well-known Binney (1978) diagram of $V_{\rm max}/\langle\sigma\rangle$ versus observed ellipticity ($\epsilon=1-b/a$). The solid line is for edge-on oblate spheroids with no anisotropy (see the discussion in § 6.1 of Cappellari et al. 2007). In such models $V_{\rm max}/\langle\sigma\rangle=0.6$ for $\epsilon=0.31$, an order of magnitude higher than the upper limit for Dragonfly 44. Dotted lines are for increasing anisotropy, here parameterized with $\delta=(2\beta-\gamma)/(2-\gamma)$ (see, e.g., Eq. 4–7 in Cappellari et al. 2007). The other data points in Fig. 13 are dwarf galaxies in and near the Local Group, taken from the compilation by Wheeler et al. (2017). The most straightforward interpretation of the distribution of the Local Group dwarfs is that they are supported by random motions rather than rotation (see Wheeler et al. 2017), and that the observed shapes are due to anisotropy in the velocity dispersion tensor. We note that these galaxies typically have low Sersic indices, similar to Dragonfly 44. The grey level and size of the symbol indicate the stellar masses of the galaxies, which range from $10^{4}-10^{8}$ M*⊙. With $M_{\rm stars}\approx 3\times 10^{8}$ M⊙* (van Dokkum et al. 2016) the stellar mass of Dragonfly 44 is just above the highest mass galaxy in the Wheeler et al. (2017) sample. Its $V_{\rm max}/\langle\sigma\rangle$ is at the low end of the distribution of dwarf galaxies, particularly when compared to more massive dwarfs which tend to have slightly more rotational support (as also noted in Wheeler et al. 2017).

Finally, we note several caveats. First, the $V/\sigma$ diagnostic diagram has mostly been applied to luminous early-type galaxies, which are baryon-dominated within the effective radius. Second, the dynamics of galaxies can be quite complex. A specific example is the early-type galaxy NGC 4550, whose relatively flat rotation curve and radially-increasing velocity dispersion profile (Rix et al. 1992) are similar to what is seen in Dragonfly 44 (see § 5.2). In this case these observed features are due to two counter-rotating disks. Although this particular scenario is unlikely, as the line profile would be flat-topped rather than peaked (§ 5.3), this illustrates the difficulty in interpreting incomplete data in a unique way.

5.2. Radial velocity dispersion profile

Having established that the galaxy does not show appreciable rotation, we assume that the kinematics can be meaningfully characterized in the elliptical apertures shown in Figs. 6 and 11. The radial velocity dispersion profile is shown in Fig. 14, in linear units (top panel) and logarithmic units (bottom panel). We find that the velocity dispersion profile gradually increases with radius, and shows no ”second order” features such as a bump on small ($0.5-1$ kpc) scales. A simple linear fit gives

[TABLE]

with $R$ in kpc and $\sigma$ in km s*-1*. Note that only random errors are taken into account in the uncertainties.

We include the UDG NGC 1052–DF2 in Fig. 14, from Danieli et al. (2019). The new data for Dragonfly 44 confirm that the two galaxies have very different kinematics, even though they have a similar luminosity, morphology, and stellar population, and both have a relatively high number of globular clusters. These two objects highlight the large object-to-object scatter that appears to exist within the UDG population (see also Spekkens & Karunakaran 2018; Toloba et al. 2018).

We also compare the Dragonfly 44 measurements to radial velocity dispersion profiles of Local Group dwarf spheroidal galaxies (dSphs). UDGs resemble dSphs in terms of their visual morphology, mean Sersic index, axis ratio distribution, surface brightness, and high dark matter fraction within the effective radius; the main (baryonic) differences are their $\sim 10\times$ larger sizes, corresponding $\sim 100\times$ larger luminosities and stellar masses, and their higher average globular cluster specific frequency (see, e.g., van Dokkum et al. 2015; Lim et al. 2018). Thin curves in the bottom panel of Fig. 14 show the radial dispersion profiles of seven classical dSphs, obtained from kinematic data of Walker et al. (2007). The velocities and velocity dispersions were added in quadrature, and rebinned to logarithmic bins with a size of 0.1 dex. The radial profile of Dragonfly 44 is not very different from scaled-up versions of those of dSph galaxies. In particular, Draco has a positive dispersion gradient that is similar to that of Dragonfly 44.

5.3. Line profile

The line profile of Dragonfly 44 is non-Gaussian, with reasonably high significance. In Fig. 15 the best-fitting $h_{3}$ and $h_{4}$ parameters are shown, for the fit to the optimally-extracted spectrum. We find $h_{3}=-0.03\pm 0.04$ and $h_{4}=0.13\pm 0.05$. The $h_{3}$ parameter is consistent with zero, indicating that the line profile is close to symmetric. This is necessarily the case if the galaxy is axisymmetric and in dynamical equilibrium, as any asymmetric deviations on one side of the galaxy should be inverted on the opposite side (see, e.g., Fig. 15 in Bender et al. 1994).

The $h_{4}$ parameter, which measures symmetric deviations from a Gaussian, is positive with a formal significance of $3\sigma$. We performed many tests to assess whether this result is reliable or driven by some subtle systematic error in the analysis. The result persists when the spectrum is split in separate wavelength intervals (see Fig. 15) or when $h_{3}$ is forced to be zero; it is seen in almost all radial bins (albeit with low significance for each individual bin); it persists when the instrumental line profile is assumed to be Gaussian instead of flat-topped; it is not sensitive to the wavelength calibration corrections described in § 3.5.2 and § 4.3; and it is insensitive to the exact template that is used.

A positive $h_{4}$ can indicate radial orbital anisotropy (see, e.g., Bender et al. 1994; Thomas et al. 2007; Amorisco & Evans 2012) but can also be due to other effects, such as flatteningf of the galaxy along the line of sight (e.g., Magorrian & Ballantyne 2001). Although there is no unique interpretation of the $h_{4}$ measurement it is difficult to reconcile with tangential anisotropy, which is required to fit the radial velocity dispersion profile for certain halo models (as we show later).

5.4. Robustness of the radial dispersion profile

We end this section with an empirical assessment of the accuracy of the error bars on the dispersion measurements. The residuals from the template fits to the spectra are consistent with the errors in the data, but as discussed in § 4.3 this does not address possible systematic errors in the dispersion measurements. Furthermore, the interpretation of the uncertainties is complicated by the fact that the spectra were resampled to a finer grid.

The major and minor axis kinematics of § 5.1 provide a test that relies only on the assumption that the dispersion profiles are symmetric with respect to the center of the galaxy. The major axis profile includes seven dispersion measurements on the “left” side and seven measurements on the “right”. Similarly, the minor axis profile includes six measurements on either side of the galaxy’s center. Assuming that the profile is symmetric, we therefore have thirteen pairs of independent measurements of the same quantity, determined from regions that span the entire detector. The scatter in the differences between these paired dispersions is $5.7\pm 1.1$ km s*-1*. The expected variation is $\langle(e_{1,i}^{2}+e_{2,i}^{2})^{0.5}\rangle_{i}=6.2\pm 0.4$ km s*-1*, where $e_{1}$ and $e_{2}$ are the uncertainties in the measurements on the two sides of the galaxy. The observed variation is consistent with the expected variation, and we conclude that the quoted uncertainties accurately reflect the observed variation between independent measurements.

6. Dynamical Mass and M/L Ratio

6.1. Is Dragonfly 44 in dynamical and structural equilibrium?

Before interpreting the kinematics we first ask whether the galaxy is in equilibrium. Some UDGs clearly are not, with the “boomerang galaxy” M101-DF4 (Merritt et al. 2016) a case in point, and it has been suggested that many Coma UDGs are in the process of tidal disruption by the cluster potential (Yozin & Bekki 2015). The rising velocity dispersion profile of Dragonfly 44 could be interpreted as evidence for such disruption, as unbound material may inflate the observed velocity dispersion at large radii. A demonstration of this was given by Muñoz, Majewski, & Johnston (2008), who reproduced the rising velocity dispersion profile of the Carina dSph in the context of such models.

However, rising velocity dispersion profiles such as that of Carina, Draco and Dragonfly 44 do not necessarily imply tidal disruption. In fact, the kinematics of Draco have been reproduced with relatively simple mass models (Kleyna et al. 2002), and the galaxy has been described as “flawless” because of its lack of detected tidal features in deep imaging data (Ségall et al. 2007). Furthermore, as noted in § 1, tidal disruption scenarios cannot be easily reconciled with the high globular cluster counts of UDGs in Coma (Lim et al. 2018) and their lack of obvious tidal features (Mowla et al. 2017).

In the specific case of Dragonfly 44 tidal heating by the Coma cluster is unlikely because it appears to live in a dynamically-cold environment. Three other low surface brightness galaxies in the vicinity of Dragonfly 44 (Dragonfly 42, DFX1, and DFX2) have redshifts, all from the DEIMOS multi-slit spectroscopy described in van Dokkum et al. (2016, 2017a), and Alabi et al. (2018). Of this sample of four faint galaxies, three (Dragonfly 42, Dragonfly 44, and DFX2) are within $\approx 100$ km s*-1* of each other (see van Dokkum et al. 2017a; Alabi et al. 2018). The velocity dispersion of the Coma cluster is $\sim 1000$ km s*-1* (Colless & Dunn 1996), and the probability that three out of four randomly selected galaxies have a velocity range $\Delta v<100$ km s*-1* is $3\times 10^{-3}$.

It is unclear whether Dragonfly 44 is in a cold clump that is falling into the cluster, a filament, or a structure that is unrelated to Coma; this can be constrained by measureing redshifts for more galaxies in the vicinity of Dragonfly 44. Irrespective of the precise interpretation, the small redshift range strongly suggests that the galaxy has not been part of the Coma cluster for a long time, and is therefore unlikely to be significantly affected by tidal heating or other cluster-driven processes. Although it is difficult to completely rule out tidal effects, in the following we assume that the galaxy is in equilibrium and that its dynamics reflect the galaxy’s gravitational potential.

6.2. Dynamical mass within the effective radius

The velocity dispersion profile extends slightly beyond the projected half-light radius of Dragonfly 44, $R_{\rm e,maj}=4.7$ kpc. As shown in Wolf et al. (2010), the luminosity-weighted velocity dispersion within the projected circularized half-light radius $R_{\rm e,c}=R_{\rm e,maj}(b/a)^{0.5}$ provides a robust estimate of the dynamical mass within the 3D half-light radius $r_{1/2}$ that is insensitive to anisotropy or the form of the density profile. We measure the luminosity-weighted dispersion directly from a luminosity-weighted extracted spectrum within the half-light radius. The resulting dispersion is $\sigma_{\rm e}=33^{+3}_{-3}$ km s*-1*. This value is lower than that reported in van Dokkum et al. (2016), who found $\sigma=47_{-6}^{+8}$ km s*-1* based on a DEIMOS spectrum in the H$\alpha$ region. A re-assessment of the 2016 analysis uncovered an error; the revised DEIMOS dispersion is $\sigma=42^{+7}_{-7}$ km s*-1*, closer to the KCWI value.999In the 2016 analysis the spectroscopic data were combined without applying barycentric velocity corrections to each individual dataset. However, as the DEIMOS data were taken over a period of several months the peak-to-peak velocity corrections are $\approx 30$ km s*-1*. As a result, the combined spectrum was slightly broadened, leading to a dispersion measurement that was biased high. After applying the required corrections we derive $\sigma=42$ km s*-1* instead of 47 km s*-1*. Despite this better agreement, the probability that the difference can be attributed to chance is only 3.3 %. We do not have an explanation for the discrepancy but speculate that it may be caused by the large weight of the Balmer H$\alpha$ line in the van Dokkum et al. (2016) analysis or systematic errors introduced by the crosstalk corrections that were needed.

The Wolf et al. (2010) estimator,

[TABLE]

gives $M(r<r_{1/2})=3.9_{-0.5}^{+0.5}\times 10^{9}$ M*⊙. The total $I_{814}$ magnitude of Dragonfly 44 is $M_{I}=-16.7$ (van Dokkum et al. 2017a), or $L_{I}=3.0^{+0.6}_{-0.6}\times 10^{8}$ L⊙, assuming a 20 % error in the total luminosity. Therefore, the $M/L_{I}$ ratio within the 3D half light radius is $M/L_{I}(r<r_{1/2})=26^{+7}_{-6}$ M⊙/L⊙*.

The expected $M/L_{I}$ ratio from the stellar population is $M/L_{I}\approx 1-1.5$, and the implication is that Dragonfly 44 is extremely dark matter dominated within its half radius. This is generally not the case for galaxies in this mass and luminosity range, as demonstrated in the left panels of Fig. 16. Here we show the relation between the dynamical $M/L_{I}$ ratio within the half-light radius as a function of the total luminosity. Grey points show samples of “normal” galaxies from Zaritsky, Gonzalez, & Zabludoff (2006), Wolf et al. (2010), and Toloba et al. (2015). The line is a fit to these samples for $L_{I}<10^{11}$ L*⊙*, of the form

[TABLE]

with $L_{8}\equiv\log L_{I}-8$, $M/L_{I}(<r_{1/2})$ the dynamical mass-to-light ratio within the half-light radius, and $L_{I}$ the total $I$ band luminosity. The rms scatter around this line is only 0.22 dex, independent of luminosity.

Large black and colored points are UDGs from Beasley et al. (2016), Toloba et al. (2018), Martín-Navarro et al. (2019), Danieli et al. (2019), van Dokkum et al. (2019), and Chilingarian et al. (2019).101010We only show objects with $R_{\rm e}>1.5$ kpc from Chilingarian et al. (2019). This is an update to Fig. 3 in van Dokkum et al. (2016) and Fig. 4 in Toloba et al. (2018), which showed the relation between $M/L$ ratio and dynamical mass. Dragonfly 44, and other UDGs with measured kinematics, are in a “no man’s land” in this parameter space, with $M/L$ ratios that are similar to much fainter and much brighter galaxies. Phrased differently, the small scatter in the well-known U-shaped distribution of galaxies in this plane is likely partially due to selection effects.

6.3. Dynamical mass within a fixed radius

It is tempting to interpret the vertical axis of the left panels of Fig. 16 in terms of the ratio between total halo mass and total stellar mass. This is often done implicitly; e.g., Martin et al. (2018) use the $M/L$ ratio within the effective radius as a probe of the total halo mass when assessing whether the UDG NGC 1052-DF2 is lacking in dark matter. However, a second parameter, the effective radius, plays an important role. The effective radius always contains 50 % of the light, but it does not contain a fixed fraction of the dark matter. At fixed halo mass, virial radius, and concentration, the enclosed dark matter mass within the half-light radius (and therefore the dynamical $M/L$ ratio) is expected to scale with that radius.111111The same is true if the mass were measured within, say, $2R_{\rm e}$ or some other measure of the optical extent of the galaxies.

We assess this effect by estimating the $M/L_{I}$ ratio within a fixed 3D radius of $r=10$ kpc for all the galaxies in the samples quoted above. The mass is extrapolated by assuming a flat rotation curve, that is,

[TABLE]

with $r_{1/2}\approx 4/3R_{\rm e,c}$. The luminosity is extrapolated by numerically integrating the Sersic (1968) profile out to $r=10$ kpc, where the Sersic index is assumed to be

[TABLE]

The results are shown in the right panels of Fig. 16. Low luminosity galaxies tend to have small effective radii, and their $M/L$ ratios within 10 kpc are much higher than those within $r_{1/2}$. High luminosity galaxies have $r_{1/2}\sim 10$ kpc, and their $M/L$ ratios within 10 kpc are similar to those within $r_{1/2}$. Because of this correlation of the effective radius with luminosity, the distribution of galaxies in the right panels is very different than in the left panel. The best fitting relation,

[TABLE]

has a factor of seven higher normalization at $L=10^{8}$ L*⊙* than the relation between $M/L(<r_{1/2})$ and luminosity.

However, Dragonfly 44, as well as other UDGs, stay at nearly the same location. As a result, Dragonfly 44 is now consistent with the relation defined by normal galaxies, whereas NGC 1052-DF2 now falls far below it: in the left panels its $M/L$ ratio is similar to that of other galaxies, but given its large effective radius its $M/L$ ratio should have been much higher if it had a normal dark matter halo. We conclude that it is hazardous to interpret the $M/L$ ratio within the effective radius in terms of halo masses. The $M/L$ ratio within a fixed large aperture should provide a better indication, but for most galaxies in Fig. 16 this represents a significant extrapolation beyond the regime where the kinematics are measured.

7. Dark matter halo fits

Here we seek to interpret the measured kinematics in the context of parameterized models for the mass distribution. In particular, we ask what classes of models can reproduce the rising velocity dispersion profile and the positive $h_{4}$ parameter, and what the implications are for the halo mass of Dragonfly 44.

7.1. Procedure

We use the methodology that is outlined in Wasserman et al. (2018b) and Wasserman et al. (2018a). Briefly, spherical mass models with a given density profile are fit to the observed velocity dispersion profile using a Bayesian Jeans modeling formalism. The mass distribution is modeled as the sum of the stellar distribution and a parameterized dark matter halo profile. For the halo we use two descriptions that are both instances of the general $(\alpha,\beta,\gamma)$ profile,

[TABLE]

with $r_{s}$ the scale radius and $\rho_{s}$ the scale density (Hernquist 1990). These profiles have a powerlaw slope $-\gamma$ on small scales and $-\beta$ on large scales, with the form of the profile near the transition controlled by $\alpha$.

The first model is a standard cuspy Navarro-Frenk-White (NFW) profile, with $\alpha=1$, $\gamma=1$, and $\beta=3$ (Navarro et al. 1997). The second is a cored model with a flatter inner density profile. This can be achieved by reducing the value of $\gamma$ in Eq. 21 while retaining $\alpha=1$ and $\beta=3$ (e.g., Zhao 1996; Wyithe, Turner, & Spergel 2001), or by using physically-motivated fitting functions such as the “CORENFW” profile of Read, Agertz, & Collins (2016). Here we use the parameterization of Di Cintio et al. (2014a), who use high resolution hydrodynamical simulations to derive “empirical” relations between the stellar-to-halo mass $X=\log(M_{\rm stars}/M_{\rm halo})$ and $(\alpha,\beta,\gamma)$. These relations are specified in Eq. 3 of Di Cintio et al. (2014a). The relation for the inner density profile is given by

[TABLE]

The Di Cintio profile is similar to NFW for very low mass galaxies and for $L_{*}$ galaxies, and has a core that is maximal for stellar masses of $\sim 10^{8.5}$ M*⊙*, that is, the stellar mass of Dragonfly 44 and of most other UDGs that have been studied in detail so far. For this mass $X\sim-2.5$ (for UDGs that lie on the stellar mass – halo mass relation) and $\gamma\sim 0.3$. We note that, as the stellar mass of Dragonfly 44 is fairly well constrained, this model does not have significantly more freedom than the standard NFW profile.

Halo masses are expressed in terms of $M_{200}$ and the concentration $c_{200}$, with

[TABLE]

with $c_{200}=r_{200}/r_{s}$. We assume the median concentration from the halo mass – concentration relation determined by Diemer & Kravtsov (2015). In addition to the halo parameters the anisotropy $\beta=1-\sigma_{\rm tan}^{2}/\sigma_{\rm rad}^{2}$ is a fit parameter.121212We use two different parameters in this paper that are both denoted “$\beta$”: the second coefficient in the ($\alpha$, $\beta$, $\gamma$) profile (Eq. 21), and Binney’s anisotropy parameter. These are both conventional expressions, and we believe changing either of them would be confusion. Hopefully it is always clear from context which $\beta$ the text is referring to. For simplicity the anisotropy is assumed to be constant with radius; models with varying anisotropy give qualitatively similar results, albeit with more freedom in the mass profiles.

The top panels of Fig. 17 illustrate the behavior of the models and what parameters can be constrained by the data. The curves are based on Jeans modeling with fixed model parameters; that is, they are not fits to the data but they do take the observed surface brightness profile of Dragonfly 44 into account. More massive halos obviously produce higher velocity dispersions, but the effect is relatively small: about a factor of $\approx 1.5$ change in velocity dispersion for a factor of 10 change in halo mass. At fixed halo mass and anisotropy the predicted dispersions are higher for NFW halos than for the cored Di Cintio halos, although the difference vanishes for low halo masses. The Di Cintio halos readily predict rising velocity dispersion profiles, particularly for halo masses $M_{200}\sim 10^{11}$ M*⊙* where the cores are maximal. However, the shape of the velocity dispersion profile is degenerate with the anisotropy parameter $\beta$: generically, radial anisotropy produces falling profiles whereas tangential anisotropy produces rising profiles. In principle, this degeneracy can be resolved by including the form of the absorption line profile in the analysis (see § 5.3 and also, e.g., Amorisco & Evans 2012). We will return to this below.

7.2. Results

The models are fit to the data using the emcee MCMC sampler (Foreman-Mackey et al. 2013), as described in Wasserman et al. (2018b). The following priors are used:

[TABLE]

where $N(\mu,\sigma^{2})$ is the normal distribution and $U({\rm min},{\rm max})$ is a uniform distribution. The only informative prior is Eq. 24. The mean $M/L_{V}$ ratio comes from stellar population synthesis modeling (e.g., van Dokkum et al. 2016), with a standard deviation obtained from Taylor et al. (2011).

Good fits are obtained for both standard NFW halos and Di Cintio cored halos, as shown in the bottom panels of Fig. 17. The distributions of MCMC samples for $M_{200}$ and $\beta$ are shown in Fig. 18. NFW halos require strong tangential anisotropy whereas the Di Cintio profiles do not. For an NFW profile the best fitting halo mass is $\log(M_{200}/{\rm M}_{\odot})=10.6^{+0.4}_{-0.3}$ and $\beta=-0.8^{+0.4}_{-0.5}$, whereas these values are $11.2^{+0.6}_{-0.6}$ and $-0.1^{+0.2}_{-0.3}$ respectively for the Di Cintio profile. One way to view these results is that the Di Cintio models “naturally” predict rising velocity dispersion profiles for the stellar mass regime of Dragonfly 44, whereas NFW profiles predict decreasing profiles unless strong tangential anisotropy is invoked.

As shown in Fig. 18 Dragonfly 44 is consistent with the stellar mass – halo mass relation of Moster et al. (2010) within $1\sigma$; this relation is very similar to that of Behroozi et al. (2013a) and others in this regime. However, the total halo mass is not particularly well constrained in either model, as the data sample the halo only out to a small fraction of the virial radius. The uncertainty in the halo mass is much larger than the probable scatter of 0.2 dex in the relation (see, e.g., Behroozi, Wechsler, & Conroy 2013b; Gu, Conroy, & Behroozi 2016). The halo mass is higher for a Di Cintio profile than for an NFW profile; in particular, a halo mass of $10^{12}$ M*⊙*, as is indicated by the globular cluster counts of Dragonfly 44 (van Dokkum et al. 2016, 2017a; Harris, Blakeslee, & Harris 2017; Forbes et al. 2018), is within the $1\sigma$ contour for the Di Cintio model whereas it falls outside of the $2\sigma$ contour for the NFW model. We note that the stellar mass is not constrained by the kinematics; the distribution of the MCMC samples simply follows the prior.

The integrated mass profiles for both types of halos are shown in Fig. 19, along with the stellar mass profile. The galaxy is dominated by dark matter at all radii even in the cored Di Cintio model. The dark matter fraction within 1 kpc is somewhat model dependent, but it is well-constrained at $f_{\rm dm}\sim 95$ % on scales of a few kpc where we have the most information. These results extend the analysis of § 6, where we showed that Dragonfly 44 has the expected $M/L$ ratio within its effective radius for a “normal” dark matter halo. UDGs should have very high $M/L$ ratios, because they are so large: galaxies such as NGC 1052-DF2, with a $M/L$ ratio within the effective radius that is not very different from other galaxies of the same luminosity (Fig. 13, and Danieli et al. 2019), are the outliers.

7.3. Anisotropy or core?

We have shown that the velocity dispersion profile of Dragonfly 44 can be reproduced with two classes of models: a standard cuspy NFW profile combined with strong tangential anisotropy, or a cored profile that is close to isotropic. As discussed in § 5.3 the shape of the absorption line profile can, in principle, constrain the degree of anisotropy and therefore help decide which of these two options is more likely. Generically, a positive $h_{4}$ parameter indicates radial rather than tangential anisotropy, as radial orbits create excesses at both zero velocity and in the wings of the velocity distribution (see, e.g., Fig. 2 in van der Marel & Franx 1993). Qualitatively, the positive $h_{4}=0.13\pm 0.05$ is inconsistent with both models, as neither model has significant radial anisotropy (see the right panel of Fig. 18). It is more inconsistent with the NFW profile, as this model requires strong tangential anisotropy.

However, non-zero kurtosis can have other causes than anisotropy in the orbital distribution, such as deviations from spherical symmetry (e.g., Read & Steger 2017) and the presence of a core in the density profile. As shown by Łokas (2002), a cored profile leads to positive kurtosis, and this may be the reason for the positive $h_{4}$ parameter that we measure. We quantify this by determining the kurtosis ($\xi_{2}=\mu_{4}/\mu_{2}^{2}$, where $\mu_{i}$ is the $i^{\rm th}$ moment) and the excess kurtosis ($\kappa=\xi_{2}-3$) from the posteriors of the model fits. Figure 20 shows the radial dependence of $\kappa$ for both choices of the density profile. The NFW model produces slightly negative kurtosis, as expected from the tangential anisotropy. The cored model produces positive kurtosis despite its nearly isotropic orbital distribution, as was also found by Łokas (2002).

The data point in Fig. 20 is the $h_{4}$ measurement from the optimally-extracted spectrum. The weighted radius of this extraction is 1.3 kpc (as determined from the mean flux and weight of each spectrum that contributes to it). The Gauss-Hermite coefficient $h_{4}$ was converted to excess kurtosis using $\kappa\approx 8\sqrt{6}h_{4}$ (van der Marel & Franx 1993). The observed kurtosis is higher than in either of the models, but closer to the cored model than to the NFW one: the distance to the Di Cintio model is $2.0\sigma$ and the distance to the NFW model is $2.4\sigma$. This discrepancy may indicate that the central density profile is even flatter than $\gamma\sim 0.3$, which it is in the Di Cintio model. It could also reflect the limitations of the assumption of spherical symmetry in the Jeans modeling (which is known to be incorrect, as Dragonfly 44 has $b/a=0.69$; see also Burkert 2017), or even the assumption that the galaxy is in equilibrium. Finally, we cannot exclude undiagnosed systematic errors in the line profile measurement. With these caveats, we cautiously conclude that the observed line profile is more consistent with a cored profile than with an NFW profile.

8. Discussion and Conclusions

In this paper we present spatially-resolved kinematics of the UDG Dragonfly 44, obtained with KCWI on the Keck II telescope. We find no evidence for rotation, which is significant as Dragonfly 44 is one of the more flattened UDGs: its axis ratio is $b/a=0.69$, whereas the median for Coma UDGs is 0.74 (van Dokkum et al. 2015). The limit that we derive is more stringent than for many other low luminosity galaxies (see Fig. 13). This result is difficult to reconcile with models in which UDGs are the high spin tail of the distribution of normal dwarf galaxies, as was proposed by Amorisco & Loeb (2016). Amorisco et al. suggest processing by the cluster environment may decrease $V/\sigma$, but as we discuss in § 6.1 Dragonfly 44 appears to be in a dynamically-cold environment.

The velocity dispersion within the effective radius is lower than what we reported in van Dokkum et al. (2016), and as discussed in § 6.2 this is partly due to an error in our earlier analysis. The corrected value is marginally consistent with our new measurement ($\sigma=42^{+7}_{-7}$ km s*-1* from DEIMOS and $\sigma=33^{+3}_{-3}$ km s*-1* from KCWI), but we cannot exclude other systematic effects. It may be that the large weight of a Balmer line (H$\alpha$) in the analysis, or the cross-talk corrections we had to apply, influenced the earlier result. The $M/L$ ratio of Dragonfly 44 is $M_{\rm dyn}/L_{I}=26^{+7}_{-6}$ M${}_{\odot}/$L*⊙* within the effective radius, and the galaxy is dominated by dark matter even in the center. This does not necessarily mean that the galaxy has an “overmassive” halo; as discussed in Fig. 6.3 and shown in Fig. 16 UDGs are expected to have very high $M/L$ within the effective radius, simply by virtue of being large. UDGs with “normal” $M/L(<R_{\rm e})$ ratios for their luminosity, such as NGC 1052-DF2, are the ones that deviate from the expectations from the stellar mass – halo mass relation.

We find that the velocity dispersion profile gradually increases with radius. The profile cannot be fit with a standard NFW halo and an isotropic velocity distribution: Dragonfly 44 either has a relatively flat density profile (a core) or strong tangential anisotropy. The Di Cintio et al. (2014a) model, with a mass-dependent core, fits the data remarkably well. It reproduces the observed velocity dispersion profile with isotropic orbits, has a halo mass that is in very good agreement with the stellar mass – halo mass relation, and is in qualitative agreement with the positive $h_{4}$ parameter. Another way to phrase this result is that the kinematics of Dragonfly 44 are similar to other galaxies in this stellar mass range, which also show evidence for cores (see Di Cintio et al. 2014b, 2014a, and references therein). Our results lend support to the model of Carleton et al. (2019), who show that cores may lead to ultra diffuse galaxy formation in clusters as a result of tidal stripping. This model is certainly consistent with the kinematics of Dragonfly 44, but perhaps not with its dynamically-cold local environment. It also remains to be seen whether such tidal models can explain the high globular cluster counts of Dragonfly 44 and other UDGs.

Irrespective of the detailed mass distribution it is clear that Dragonfly 44 has a gravitationally-dominant dark matter halo, similar to many other UDGs (Beasley et al. 2016; Toloba et al. 2018; Martín-Navarro et al. 2019), and in apparent contrast to the UDGs NGC 1052-DF2 and NGC 1052-DF4 (see, e.g., van Dokkum et al. 2018b, 2019; Martin et al. 2018; Famaey, McGaugh, & Milgrom 2018; Emsellem et al. 2019; Danieli et al. 2019). With a robust velocity dispersion measurement for NGC 1052-DF2 from stellar kinematics (Danieli et al. 2019), the identification of a second galaxy in the same class (NGC 1052-DF4; van Dokkum et al. 2019), and the results presented in this study, there can be little doubt that large, diffuse, spheroidal galaxies with stellar masses of a few $\times 10^{8}$ M*⊙* have a remarkable range in their kinematics and, hence, dark matter properties on kpc scales (see Fig. 16). This qualitatively addresses a point raised by Kroupa (2012), who noted that the then-observed low scatter in the SMHM relation is difficult to explain in the standard cosmological model. Similar arguments have been made by McGaugh (2012) on the basis of the Tully-Fisher relation.

The converse of this argument is that a high scatter is difficult to explain in alternatives to dark matter, such as Modified Newtonian Dynamics (MOND; Milgrom 1983) and Emergent Gravity (Verlinde 2017). The observed dispersion of Dragonfly 44 is higher than the MOND prediction whereas the dispersions of NGC 1052-DF2 and NGC 1052-DF4 are lower. Specifically, for a MOND acceleration scale of $a_{0}=3.7\times 10^{3}$ km2 s*-2* kpc*-1* the predicted velocity dispersion of Dragonfly 44 is $\sigma_{\rm M}\approx(0.05GM_{\rm stars}a_{0})^{1/4}\approx 23$ km s*-1*, lower than the luminosity-weighted dispersion within the effective radius of $\sigma_{\rm e}=34^{+3}_{-3}$ km s*-1*. The predicted dispersions for NGC 1052-DF2 and NGC 1052-DF4 are of the same order, whereas the observed dispersions are $<10$ km s*-1* (Danieli et al. 2019; van Dokkum et al. 2019). The “external field effect” (Famaey et al. 2018) mitigates this tension, but it is difficult to explain dispersions as low as $5-10$ km s*-1* for these galaxies even when this effect is maximal (see Müller, Famaey, & Zhao 2019). A possible way to reconcile alternative dark matter models with UDG kinematics is to invoke strong variations in the stellar initial mass function (see, e.g., Conroy & van Dokkum 2012; Geha et al. 2013), such that Dragonfly 44 has a bottom-heavy mass function and the NGC 1052 dwarfs are bottom-light. There is no prior motivation for these specific variations, but they may be testable.

The fact that it is dark matter dominated on all scales makes Dragonfly 44, and other UDGs like it, important in the quest to understand the distribution of dark matter on $\lesssim 1$ kpc scales. As noted above, there is a long history of using galaxies with a low baryonic density to constrain the density profiles of dark matter halos (see, e.g., Aaronson 1983; de Blok et al. 2001; Kleyna et al. 2002; Swaters et al. 2003, and many others). As UDGs such as Dragonfly 44 are dark matter dominated on all scales, as shown explicitly in Fig. 19, they offer a “pristine” view of their dark matter even on small spatial scales, which is unusual for galaxies with stellar masses of $M_{\rm stars}=10^{8}-10^{9}$ M*⊙* (and higher). In particular, dwarf elliptical galaxies, “classical” low surface brightness disk galaxies, and gas-rich dwarf galaxies all have typical $M/L$ ratios in the range $5-10$ within the optical extent, and even lower values in the center (Geha, Guhathakurta, & van der Marel 2002; Swaters et al. 2003; Zaritsky et al. 2006; Wolf et al. 2010).

A particularly interesting deviation in the dark matter profile occurs in “fuzzy” dark matter models, where the dark matter particle is an ultra-light axion with a de Broglie wavelength of 100s of parsecs (e.g., Marsh & Silk 2014). At low halo masses the soliton core in these models is difficult to distinguish, but at higher masses ($10^{11-12}$ M*⊙) the predicted density profile can display a characteristic bump on small scales. As pointed out by Hui et al. (2017), massive UDGs may be able to constrain such models. In Fig. 21 we show predicted velocity dispersion profiles for a galaxy with the surface brightness profile of Dragonfly 44, using an NFW + soliton model of the form proposed by Marsh & Pop (2015). We should not expect a bump on 500 pc scales in Dragonfly 44, as the soliton core becomes a distinct feature only for halo masses of $\gtrsim 10^{12}$ M⊙, or velocity dispersions of $\sim 60$ km s-1*. In that regime the effects of the soliton can, with sufficiently accurate data, be distinguished from those of anisotropy and variations in the halo mass (see also Robles, Bullock, & Boylan-Kolchin 2019). In this context the recent announcement of a UDG with an apparent stellar dispersion of $\sigma=56\pm 10$ km s*-1* (Martín-Navarro et al. 2019) is exciting, as it suggests that Dragonfly 44 does not define the upper end of the halo mass range of UDGs. We also note that, even though we cannot determine whether soliton models provide a better fit than standard NFW models, we can place constraints on the particle mass in the context of fuzzy dark matter models. These quantitative constraints are given in a companion paper (Wasserman et al. 2019).

This paper is dedicated to Mariella Silvia, a young woman with a passion for astronomy who has inspired us with her fortitude. Support from HST grant HST-GO-14643 and NSF grants AST-1312376, AST-1616710, AST-1518294, and AST-1613582 is gratefully acknowledged. This work was partially supported by a NASA Keck PI Data Award, administered by the NASA Exoplanet Science Institute. The data presented herein were obtained at the W. M. Keck Observatory, which is operated as a scientific partnership among the California Institute of Technology, the University of California and the National Aeronautics and Space Administration. The Observatory was made possible by the generous financial support of the W. M. Keck Foundation. We are grateful to the staff of Keck Observatory, and in particular Luca Rizzi, for their excellent support and help. The authors wish to recognize and acknowledge the very significant cultural role and reverence that the summit of Maunakea has always had within the indigenous Hawaiian community. We are most fortunate to have the opportunity to conduct observations from this mountain. AJR is a Research Corporation for Science Advancement Cottrell Scholar.

Appendix A A. Radial photometry

In the sky subtraction process the surface brightness profile of the galaxy is used to determine the absolute background level within the KCWI field of view (see § 3.4.3). This profile is derived from the HST WFC3/UVIS $V_{606}$ image that was described in van Dokkum et al. (2017a). The central wavelength of the $V_{606}$ filter is $0.59\,\mu$m, slightly larger than the central wavelength of the KCWI spectrum ($0.51\,\mu$m), and the sky subtraction could be in error if the galaxy has a strong color gradient. Here we assess whether such a gradient exists, making use of the fact that Dragonfly 44 was also observed in the $I_{814}$ filter (see van Dokkum et al. 2017a).

The $V_{606}$ image and surface brightness profile are shown in the top panels of Fig. 22. A color image created from the $V_{606}$ and $I_{814}$ data is shown in the bottom left panel. To measure the color gradient the $I_{814}$ image is smoothed to the ground-based seeing and resampled to the KCWI pixel scale, mimicking the procedure that is used for the $V_{606}$ image. Next, $I_{814}$ fluxes are measured in the same elliptical apertures (with the same masking of stars, globular clusters, and background galaxies) as are used for the extraction of the KCWI spectra and the measurements of the $V_{606}$ fluxes. Finally, colors are measured using the WFC3 zeropoints,131313The WFC3/UVIS $V_{606}$ and $I_{814}$ AB zeropoints are 26.103 and 25.139 respectively. with a $-0.01$ mag correction to account for Galactic reddening. Uncertainties are dominated by the uncertainty in the background level in the HST images, and determined from the variation in the mean flux in empty apertures placed near Dragonfly 44. The results are shown in the bottom right panel of Fig. 22 and listed in Table 3. There is no evidence for a gradient at $r>3\arcsec$, the part of the profile that is used to scale the profile derived from the KCWI data.