Andromeda XXI -- a dwarf galaxy in a low density dark matter halo
Michelle L. M. Collins, Justin I. Read, Rodrigo A. Ibata, R. Michael Rich, Nicolas F. Martin, Jorge Peñarrubia, Scott C. Chapman, Erik J. Tollerud, Daniel R. Weisz
MMNRAS , 1–16 (2020) Preprint 25 February 2021 Compiled using MNRAS L A TEX style file v3.0
Andromeda XXI – a dwarf galaxy in a low density dark matter halo
Michelle L. M. Collins ★ , Justin I. Read , Rodrigo A. Ibata , R. Michael Rich ,Nicolas F. Martin , Jorge Peñarrubia , Scott C. Chapman , , Erik J. Tollerud ,Daniel R. Weisz Physics Department, University of Surrey, Guildford, GU2 7XH, UK Observatoire de Strasbourg,11, rue de l’Université, F-67000, Strasbourg, France Department of Physics and Astronomy, University of California, Los Angeles, CA 90095-1547 Institute for Astronomy, University of Edinburgh, Royal Observatory, Blackford Hill, Edinburgh, EH9 3HJ, UK NRC Herzberg Institute for Astrophysics, 5071 West Saanich Road, Victoria, V9E 2E7, British Columbia, Canada Dalhousie University Dept. of Physics and Atmospheric Science Coburg Road Halifax, B3H1A6, Canada Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA Department of Astronomy, University of California Berkeley, Berkeley, CA 94720, USA
Accepted XXX. Received YYY; in original form ZZZ
ABSTRACT
Andromeda XXI (And XXI) has been proposed as a dwarf spheroidal galaxy with a centraldark matter density that is lower than expected in the Standard Λ Cold Dark Matter ( Λ CDM)cosmology. In this work, we present dynamical observations for 77 member stars in this system,more than doubling previous studies to determine whether this galaxy is truly a low densityoutlier. We measure a systemic velocity of 𝑣 𝑟 = − . ± . − and a velocity dispersionof 𝜎 𝑣 = . + . − . km s − , consistent with previous work and within 1 𝜎 of predictions madewithin the modified Newtonian dynamics framework. We also measure the metallicity of ourmember stars from their spectra, finding a mean value of [ Fe / H ] = − . ± . 𝜌 DM ( ) = . + . − . × M (cid:12) kpc − at 68% confidence, and a density attwo half light radii of 𝜌 DM ( . ) = . + . − . × M (cid:12) kpc − at 68% confidence. These areboth a factor ∼ − Λ CDM.We show that this cannot be explained by ‘dark matter heating’ since And XXI had too littlestar formation to significantly lower its inner dark matter density, while dark matter heatingonly acts on the profile inside the half light radius. However, And XXI’s low density can beaccommodated within Λ CDM if it experienced extreme tidal stripping (losing >
95% of itsmass), or if it inhabits a low concentration halo on a plunging orbit that experienced repeatedtidal shocks.
Key words: keyword1 – keyword2 – keyword3
In recent years, a number of galaxies with surprisingly low centraldark matter densities have been discovered in the Local Group. Ithas long been observed that low-surface brightness galaxies havelower central densities than predicted by pure cold dark matter(CDM) simulations. Early work showed that their density profilesare more consistent with centrally flat cores rather than the steepcusps predicted by CDM (e.g. Flores & Primack 1994; de Bloket al. 2001; Read et al. 2017). More recently, comparisons of thedwarf spheroidal (dSph) satellites of the Milky Way (MW) andAndromeda (M31) with dark matter only simulations show that theobserved dSphs also possess a lower density within in their half- ★ E-mail: [email protected] (MLMC) light radii ( 𝑟 half ) as a population than expected, and this is referredto as the “Too Big to Fail” (TBTF) problem (e.g. Boylan-Kolchinet al. 2013; Tollerud et al. 2012; Collins et al. 2014). For mostof these systems, it is likely that repeated gravitational potentialfluctuations, driven by feedback from star formation, has played arole in softening their central cusps over time into flatter cores (e.g.Navarro et al. 1996a; Read & Gilmore 2005; Pontzen & Governato2012; Zolotov et al. 2012; Brooks & Zolotov 2014; Oñorbe et al.2015; Read et al. 2016, 2019). This process – that has become knownas ‘dark matter heating’ – is only effective in galaxies above somestellar mass-to halo mass ratio threshold ( 𝑀 ∗ / 𝑀 ∼ × − ; DiCintio et al. 2012), and for systems with extended star formation(Read et al. 2016, 2019).However, the low density problem seems particularly acutein some of the more diffuse companions of the Milky Way and © a r X i v : . [ a s t r o - ph . GA ] F e b M. L. M. Collins et al.
Andromeda. Around our Galaxy, two extremely diffuse satellites –Antlia 2 and Crater II (Torrealba et al. 2016, 2019) – have been un-covered, whose half light radii far exceed that expected from theirsimilarly bright counterparts. Both objects have a surface bright-ness of 𝜇 >
30 mag arcsec − , and effective radii in excess 1 kpc( 𝑟 half = . . 𝜇 = . ± . − and 𝑟 half = . + . − . kpc (Martinet al. 2016). All three of these “feeble giants” also appear to re-side in surprisingly low mass halos. Through measured velocitydispersions, their mass within 𝑟 half are far lower than expected forsystems of their size or brightness (Caldwell et al. 2017; Fu et al.2019; Torrealba et al. 2019; Collins et al. 2020), raising questionsabout whether they can be understood in the context of the Λ ColdDark Matter ( Λ CDM) framework. For these systems, it is unlikelythat star formation feedback alone can explain their low densities.Thus far, detailed star formation histories are not available for theseobjects, but even assuming highly efficient dark matter heating, theircurrent sizes, surface brightness and velocity dispersions cannot bereproduced (Torrealba et al. 2019). Instead, it is assumed they mustalso have undergone extreme tidal interactions while orbiting theirhost galaxy.There are two main tidal effects that can act to lower the den-sity of orbiting satellites: tidal stripping and tidal shocking. Tidalstripping occurs when the tidal force from the host galaxy exceedsthat from the satellite, causing dark matter and/or stars to becomeunbound (e.g. von Hoerner 1957; Read et al. 2006b). Since thispeels away stars and dark matter from the outside in, it only lowersthe inner density after extreme stripping has occurred. For cuspydark matter profiles, this means losing >99% of the satellite’s initialmass (Peñarrubia et al. 2008, 2010a; Errani & Peñarrubia 2020).Shallower or cored dark matter profiles require less extreme massloss (Read et al. 2006a; Peñarrubia et al. 2010a; Brooks & Zolo-tov 2014). Tidal shocking occurs for satellites moving on eccentricorbits if the external gravitational field changes more rapidly thanthe internal dynamical time of the stars and/or dark matter. Thismeans that tidal shocks are maximised at pericentre, where the ex-ternal field changes most rapidly (e.g. Spitzer 1958; Gnedin et al.1999), and for low density satellites, since the internal orbit timegoes as 𝜌 − / , where 𝜌 is the total density. For this reason, in thecontext of Λ CDM, satellites on plunging orbits can only efficientlylower their density through tidal shocks if they have a central darkmatter core (e.g. Read et al. 2006a; Errani et al. 2017; Errani &Peñarrubia 2020; van den Bosch & Ogiya 2018), or if they inhabitlow-concentration dark matter halos before infall (Amorisco 2019).Unlike dark matter heating, the combination of tidal stripping andshocking will lower the satellite’s density at all radii, not just in thecentre (e.g. Kazantzidis et al. 2004; Read et al. 2006a). However,it can be challenging to detect, observationally, when this has oc-curred. The tell-tale signatures of tides – stellar streams, distortedouter stellar isophotes and/or tangential anisotropy – all manifest inthe low surface brightness stellar outskirts (e.g. Read et al. 2006a;Ural et al. 2015; Amorisco 2019; Genina et al. 2020).For Antlia 2 and Crater II, proper motions can be measuredfrom Gaia DR2 (Fritz et al. 2018; Fu et al. 2019; Torrealba et al.2019), and both are consistent with being on radial orbits that bringthem within a few 10s of kpc from the Galactic centre – the regime inwhich tidal shocks are likely important (Read et al. 2006a). As such,it has been argued that tidal processes govern the evolution of thesegalaxies (Sanders et al. 2018; Amorisco 2019). Indeed, a plungingorbit for Antlia 2 could also explain observed disturbances in the outer HI disc of the Milky Way (Chakrabarti et al. 2019). A preciseorbit for Andromeda XIX has not been determined, but from itsstructural and dynamical properties, it is possible that this galaxy isalso on a similarly plunging orbit (McConnachie et al. 2008; Martinet al. 2016; Collins et al. 2013, 2020). To determine whether all lowmass extended galaxies can be explained in this manner in Λ CDM,however, we need to study more systems.In this work, we turn our attention to another significantlow-mass outlier, Andromeda XXI (And XXI). This M31 satel-lite has a luminosity of 𝐿 = . × L (cid:12) and a half-light radiusof 𝑟 half = + − pc (Martin et al. 2016; Weisz et al. 2019a),approximately 3 times more extended than other galaxies of com-parable luminosity. And, based on an earlier dynamical study, it hasa low central velocity dispersion of 𝜎 𝑣 = . + . − . km s − , consis-tent with residing in a very low density halo (Collins et al. 2013,2014). However, these findings were based on dynamics of only32 member stars, meaning that detailed dynamical modelling of itshalo could not be carried out. In this work, we reanalyse the massfor And XXI using dynamics for 77 member stars, modelling itsradial dark matter density profile for the first time with an updatedversion of the GravSphere Jeans code (Read et al. 2017; Read &Steger 2017; Read et al. 2018; Gregory et al. 2019; Genina et al.2019). We combine our dynamical analysis with the measured starformation history for this galaxy from Weisz et al. (2019b) to de-termine whether its central density can be explained by dark matterheating driven by stellar feedback (e.g. Navarro et al. 1996b; Read& Gilmore 2005; Pontzen & Governato 2012; Read et al. 2016,2019), and to assess what role – if any – tidal interactions may playin explaining its properties.This paper is organised as follows. In §2, we start by dis-cussing our observations. We then present the observed dynamicsand metallicity of our member stars in §3. In §4, we describe severalimprovements we have made to the GravSphere code to improve itsperformance for modelling systems with small numbers of stars, andwe present our dynamical modelling of And XXI (tests of our newmethod on mock data are included for completeness in AppendixA). In §5, we discuss our results and set them in the context of priorwork in the literature. Finally, in §6 we present our conclusions. Spectroscopic observations of And XXI were undertaken using theDeep Extragalactic Imaging Multi-Object Spectrograph (DEIMOS,Faber et al. 2003; Cooper et al. 2012), which is mounted on the KeckII telescope in Mauna Kea. The multi-object mode of DEIMOS al-lows us to simultaneously observe ∼
150 stellar targets in a singlepointing, spread across the significant field of view of DEIMOS(16 (cid:48) × (cid:48) , which translates to ∼ . × . + − pc (Martin et al. 2016;Weisz et al. 2019a) fits comfortably within a single mask. The ob-servations presented in this work are taken from 3 DEIMOS masks,which were observed in September 2011 (previously presented inCollins et al. 2013) and October 2013 as part of the Z-PAndASspectroscopic survey. The instrumental set-up for each mask wasidentical, and used the 1200 line mm − grating (resolution of 1.4ÅFWHM). A central wavelength of 7800Å was used, allowing usspectral coverage from ∼ − 𝜆 ∼ MNRAS000
150 stellar targets in a singlepointing, spread across the significant field of view of DEIMOS(16 (cid:48) × (cid:48) , which translates to ∼ . × . + − pc (Martin et al. 2016;Weisz et al. 2019a) fits comfortably within a single mask. The ob-servations presented in this work are taken from 3 DEIMOS masks,which were observed in September 2011 (previously presented inCollins et al. 2013) and October 2013 as part of the Z-PAndASspectroscopic survey. The instrumental set-up for each mask wasidentical, and used the 1200 line mm − grating (resolution of 1.4ÅFWHM). A central wavelength of 7800Å was used, allowing usspectral coverage from ∼ − 𝜆 ∼ MNRAS000 , 1–16 (2020) ndromeda XXI – a detailed study Table 1.
Details of And XXI spectroscopic observations. A total of 88 And XXI velocities were measured, for 77 independent stars (11 repeat measurements).Mask name Date RA Dec Position angle (deg) Exposure time ( 𝑠 ) No. targets No. members7And21 26 Sep 2011 23:54:47.70 42:28:33.6 180 3600 157 32A21maj 01 Oct 2013 23:54:47.70 42:28:15.0 147 7200 112 26A21min 01 Oct 2013 23:54:47.70 42:28:15.0 57 7200 110 30 strong absorption feature that we use to determine both velocitiesand metallicities for our observed stars. The exposure time for themask observed in 2011 was 3600s (in 3 × × (cid:48)(cid:48) , 0.8 (cid:48)(cid:48) and0.7 (cid:48)(cid:48) , resulting in typical S:N values of ∼ − . Finally, we also correctthese velocities to the heliocentric frame.As we aim to combine kinematic data from three masks, ob-served at different times and in different conditions, it is imperativethat we correct our final velocities for any systematic shifts caused bymis-alignments of stars within the slits themselves. Such misalign-ments can be caused by astrometry errors or from a slight offset ofthe position angle of the mask on the sky, and can introduce velocityshifts in our spectra of up to ∼
15 km s − (Collins et al. 2013). Suchshifts can be corrected for using atmospheric telluric absorptionlines, which are imprinted on each of our spectra. As they originatefrom the atmosphere, they should always be observed at their rest-frame wavelength. By cross-correlating the science spectra with amodel telluric spectrum, one can measure and correct for any ve-locity shifts introduced by the types of misalignments mentionedabove.Within our catalogue, we have 23 stars with more than one ve-locity measurement. We use these to check our velocity calibrationsdescribed above are reasonable. We find a mean (median) offsetvelocity for these stars of 𝑣 offset = . ( . ) km s − , both of whichare within our typical velocity uncertainties. We provide all mea-sured velocities, metallicities and other reduced properties for ourmember and non-member stars in an electronic file, available onlineat the journal website. The reduced spectra are also available uponrequest, while the raw data are available through the Keck archive. Subaru Suprime-Cam imaging of And XXI were undertaken on21-22 August 2009 in Cousins 𝑉 and 𝑖 𝑐 filters, in photometricconditions with an average seeing of ∼ (cid:48)(cid:48) . A single field was ( g - i ) i P m e m b e r Figure 1.
A colour magnitude diagram for And XXI, constructed from Sub-aru Suprime-cam data. Likely members of And XXI with spectroscopic dataare highlighted as large points, colour-coded by probability of membership(derived from imaging and spectroscopic data). The red giant branch canbe clearly seen. The solid line is an isochrone from the Dartmouth databasewith [Fe/H] = − . [ 𝛼 / Fe] = .
4, age = 12 Gyr, shifted to the distancemodulus of And XXI, 𝑚 − 𝑀 = . observed, with 5 × 𝑠 in the 𝑉 − band, and 20 × 𝑠 in the 𝑖 𝑐 − band. The data were processed with the CASU pipeline (Irwin& Lewis 2001), which debiased, flat-fielded, trimmed and gaincorrected the images. A catalogue was generated, and each sourcewas morphologically classified as either stellar, non-stellar or noise-like. Finally, as we wished to combine these with PAndAS imagingof And XXI (Martin et al. 2016), we transformed these into CFHT-MegaCam 𝑔 and 𝑖 − band magnitudes, which also allowed for a fullcalibration of our catalogue, using the following calibration, 𝑔 = 𝑉 − . × ( 𝑉 − 𝑖 𝑐 ) + .
060 (1) 𝑖 = 𝑖 𝑐 + . × ( 𝑉 − 𝑖 𝑐 ) + . . (2)Finally, we extinction correct all our data using the dust mapsof Schlegel et al. (1998). Our final colour magnitude data for thesedata are shown in fig. 1. In this section we present our measurements of the basic chemo-dynamic properties of And XXI. All numerical results are reportedin table 2. Throughout our analysis we make use of various python
MNRAS , 1–16 (2020)
M. L. M. Collins et al.
Table 2.
The properties of And XXI.Property 𝛼, 𝛿 (J2000) 𝑚 𝑉 , 𝑎 . ± . 𝑀 𝑉 , 𝑎 − . ± . 𝑏 + − 𝑟 half (arcmin) 𝑎 . + . − . 𝑟 half (pc) 𝑎 ± 𝜇 (mag per sq. arcsec) 𝑎 . ± . 𝐿 ( 𝐿 (cid:12) ) 𝑎 . + . − . × 𝑣 𝑟 ( km s − ) 𝑐 − . ± . − 𝜎 𝑣 ( km s − ) 𝑐 . + . − . kms − 𝑀 ( 𝑟 < . 𝑟 half ( 𝑀 (cid:12) ) 𝑐 . + . − . × [ 𝑀 / 𝐿 ] half ( 𝑀 (cid:12) / 𝐿 (cid:12) ) 𝑐 + − [ Fe / H ] (dex) 𝑐 − . ± . 𝑎 Martin et al. (2016), 𝑏 Conn et al. (2012), 𝑐 This work packages, specifically NumPy, SciPy, Astropy and Matplotlib(Oliphant 06 ; Jones et al. 2001; Astropy Collaboration et al. 2018;Hunter 2007).
To investigate the mass profile of And XXI, we first identify whichstars observed with DEIMOS are probable members following theprocedure outlined in Collins et al. (2013) and Collins et al. (2020).Briefly, this method assigns probability of membership based onthree criteria: (1) the stars position on the colour magnitude diagramof the dwarf galaxy, 𝑃 CMD (2) the distance of the star from the centreof the dwarf galaxy 𝑃 dist and (3) the velocity of the star, 𝑃 vel . Theprobability of membership can then be expressed as a multiplicationof these three criteria: 𝑃 member ∝ 𝑃 CMD × 𝑃 dist × 𝑃 vel (3) 𝑃 CMD is determined using the colour magnitude diagram(CMD) of And XXI. We implement a method based on that ofTollerud et al. (2012), using an isochrone to isolate those stars mostlikely to be associated with And XXI (see figure 1). We use anold, metal poor isochrone from the Dartmouth stellar evolution-ary models ([Fe/H] = − . [ 𝛼 / Fe] = .
4, age = 12 Gyr, shifted to adistance modulus of 𝑚 − 𝑀 = .
65 Dotter et al. 2008; Weisz et al.2019a) that well represents the RGB of the dwarf galaxy. We thenmeasure the minimum distance of a star from this isochrone ( 𝑑 min ),and assign a probability using the following equation: 𝑃 CMD = exp (cid:32) − 𝑑 𝜎 (cid:33) (4)where 𝜎 CMD = . 𝑃 dist is determined using the known radialsurface brightness profile of the dwarf, modelled as a Plummer pro-file, using the half-light radius and ellipticity parameters for AndXXI as determined from PAndAS data (Martin et al. 2016). 𝑃 vel isdetermined by simultaneously fitting the velocities of all observedstars assuming that 3 dynamically distinct, Gaussian componentsare present: the MW foreground contamination ( 𝑃 MW , with sys-temic velocity 𝑣 MW and velocity dispersion of 𝜎 𝑣, MW ), the M31halo contamination ( 𝑃 M31 , with systemic velocity 𝑣 M31 and veloc-ity dispersion of 𝜎 𝑣, M31 ), and And XXI, 𝑃 A21 , with an arbitrary
Table 3.
Prior values used in our emcee analysisComponent Prior 𝑣 𝑟 ( km s − ) prior 𝜎 𝑣 ( km s − ) And XXI − < 𝑣 𝑟 < −
340 0 < 𝜎 𝑣 < − < 𝑣 𝑟 < < 𝜎 𝑣 < − < 𝑣 𝑟 < −
220 0 < 𝜎 𝑣 < systemic velocity, 𝑣 𝑟 and velocity dispersion 𝜎 𝑣 . As such, the prob-ability that a given star belongs to component 𝑦 (where 𝑦 is the MW,M31 or And XXI) is 𝑃 𝑦,𝑖 = √︃ 𝜋 ( 𝜎 𝑣,𝑦 + 𝛿 𝑣𝑟,𝑖 ) × exp − (cid:169)(cid:173)(cid:173)(cid:171) 𝑣 𝑦 − 𝑣 𝑟,𝑖 √︃ 𝜎 𝑣,𝑦 + 𝛿 𝑣𝑟,𝑖 (cid:170)(cid:174)(cid:174)(cid:172) , (5)where 𝑣 𝑟,𝑖 and 𝛿 𝑣𝑟,𝑖 are the velocity and uncertainty for the 𝑖 − thstar. In reality, the Milky Way distribution is more complex thana single Gaussian assumes (see, e.g. Gilbert et al. 2006; Collinset al. 2013), however as And XXI is well-separated kinematicallyfrom the Milky Way population, this assumption does not affectour membership determination. The overall likelihood function canthen be simply written as a combination of these three components, L 𝑖 ( 𝑣 𝑟,𝑖 , 𝛿 𝑣𝑟,𝑖 |P) = ( − 𝜂 MW − 𝜂 M31 ) × 𝑃 A21 + 𝜂 MW × 𝑃 MW + 𝜂 M31 × 𝑃 M31 , (6)where 𝜂 MW and 𝜂 M31 are the fraction of our sample found withinthe Milky Way and M31 halo components of the model. We usethe MCMC emcee code (Foreman-Mackey et al. 2013) to explorea broad parameter space for these components. We set uniformpriors for each of our parameters (the velocities and dispersionsfor each population, see table 3 for details). In addition, we set0 < 𝜂 MW / M31 <
1. The dynamical values we measure using thistechnique will likely resemble the final values derived from a prob-ability weighted analysis, but without including prior informationabout CMD and spatial positions of stars, these may be biased by theinclusion of contaminant M31 halo stars. This is thus the first stepin determining the kinematic parameters. This first-pass kinematicanalysis gives 𝑣 𝑟 − . ± . − and 𝜎 𝑣 = . + . − . km s − forAnd XXI.We then combine these three probabilities ( 𝑃 CMD , 𝑃 dist , 𝑃
A21 )and calculate the likelihood of a given star being a member of AndXXI. The results of this analysis are shown in fig. 2, where wedisplay the kinematic distribution of all stars in the left hand panel,with the most likely stars highlighted in red. The right panel showsthe distance of the stars from the centre of And XXI as a functionof their velocity, colour coded by their probability of membership, 𝑃 member . The dashed line represents 𝑟 half . Already an interestingtrend can be seen here, as the stars within the half-light radius showa far broader dispersion than those located further out. While only10 stars contribute to this outer population, it is certainly striking.In total, we identify 77 stars with 𝑃 member > . 𝑃 member < .
02, significantly lower than our lowest probabilitymember.We next use these probabilities as weights in our analysis. Todetermine the systemic velocity ( 𝑣 𝑟 ) and velocity dispersion ( 𝜎 𝑣 ) MNRAS000
02, significantly lower than our lowest probabilitymember.We next use these probabilities as weights in our analysis. Todetermine the systemic velocity ( 𝑣 𝑟 ) and velocity dispersion ( 𝜎 𝑣 ) MNRAS000 , 1–16 (2020) ndromeda XXI – a detailed study
500 400 300 200 100 0 v (kms ) N *
450 400 350 300 250 200 150 100 50 v (km/s) D i s t a n ce ( a r c m i n ) P m e m b e r Figure 2. Left:
Velocity histogram of all stars observed in our 3 DEIMOS fields, highlighting And XXI members ( 𝑃 member > .
1) in red.
Right:
The distancefrom the centre of And XXI as a function of velocity. All probable And XXI members are colour coded by probability of membership. The half-light radius isindicated as a dashed line, with the grey shaded region indicating the 1 𝜎 uncertainties on this value (Martin et al. 2016). Interestingly, a clear difference canalready be seen at this boundary, where the stars within 𝑟 half show a significantly larger dispersion than those beyond 𝑟 half . of And XXI, we define a likelihood function, L , that describes asingle Gaussian population:log L 𝐴 (cid:18) 𝑣 𝑟,𝑖 |(cid:104) 𝑣 𝑟 (cid:105) , 𝜃, √︃ 𝜎 𝑣𝑟 + 𝛿 𝑣𝑟,𝑖 (cid:19) = − 𝑁 ∑︁ 𝑖 = log ( 𝜎 )+ (cid:32) ( 𝑣 𝑟 − 𝑣 𝑟,𝑖 ) 𝜎 (cid:33) + log ( 𝜋 ) + log ( 𝑃 member ,𝑖 ) (7)where 𝜎 = √︃ 𝜎 𝑣 + 𝛿 𝑣,𝑖 is the combination of the underlying veloc-ity dispersion of And XXI and the velocity uncertainty of individualstars ( 𝛿 𝑣,𝑖 ). 𝑃 member ,𝑖 is the probability of membership of the 𝑖 − thstar. We again use emcee to investigate the parameter space, us-ing the results from our first pass as the initial starting guess forAnd XXI’s parameters, and the same uniform priors. The result-ing posterior distribution can be seen in fig. 3. Both the systemicvelocity and velocity dispersion are well resolved, giving medianvalues of 𝑣 𝑟 = − . ± . − and 𝜎 𝑣 = . + . − . km s − ,where the uncertainties are the 68% percentiles of the posteriordistributions. These values are nearly identical to our 3 compo-nent, kinematic-only approach above. They are also consistent withC13 values derived from 32 stars of 𝑣 𝑟,𝐶 = − . ± . − and 𝜎 𝑣,𝐶 = . + . − . km s − . Given this similarity, it is likely thatAnd XXI still resides in a low mass halo, as found by C13. Wereturn to this issue in § 4.To investigate the visually striking change in dispersion at 𝑟 half seen in fig. 2, we run the same MCMC process on stars insideand outside this boundary. The results for 𝑣 𝑟 are fully consistentfor both samples, but 𝜎 𝑣 changes dramatically. For the 67 centralstars, we measure 𝜎 𝑣 ( 𝑟 < 𝑟 half ) = . + . − . km s − , consistent withour main analysis. However, for the 10 stars at larger radii we find 𝜎 𝑣 ( 𝑟 > 𝑟 half ) = . + . − . km s − . This is a puzzling difference. As v A ( km / s ) = − . + . − . − − − − − v A ( km / s ) σ A ( k m / s ) σ A ( km / s ) σ A ( km / s ) = 6 . + . − . Figure 3.
A corner plot showing the results of our kinematic analysis usingemcee. Both velocity and velocity dispersion are well resolved, giving 𝑣 𝑟 = − . ± . 𝜎 𝑣 = . + . − . . we only have 10 stars in the outer sample, it is hard to draw a strongconclusion about this significant drop at this time. With our sample of 77 stars in And XXI, we can measure how theradial velocity and velocity dispersion behave as a function of radius,and attempt to map its mass profile. Typically, dSphs are dispersion
MNRAS , 1–16 (2020)
M. L. M. Collins et al. v e l ( k m / s ) And XXI - radial profile
Radius (arcmin) V e l o c it y d i s p e r s i on ( k m / s ) This workC13
Figure 4.
Radial velocity profile (top) and velocity dispersion profile (bot-tom) for And XXI. Individual stellar velocities are also shown in the toppanel as small grey points, with those from our original survey highlightedas encircled points. The blue dashed lines and shaded regions represent thesystemic velocity and dispersion for the whole dataset, plus their 1 𝜎 confi-dence intervals as calculated using our MCMC approach. In the lower panel,we also show the average dispersion derived from our previous work (C13)as a dot-dashed line and grey shaded region. Our new value is slightly higher,but perfectly consistent within the uncertainties with C13. From the velocityprofile, it appears that this system is not in equilibrium, as variations aroundthe systemic are often significant ( > − 𝜎 ). The dispersion profile is alsosurprising. As all our data are within ∼ × 𝑟 half , one expects the profile toremain constant as a function of radius owing to the extended dark matterhalo. Instead, we see that the majority of the profile is flat, with a dispersionof ∼ − − , with the exception of a far hotter bin at ∼ (cid:48) . This bin hasa dispersion that is inconsistent with the others at a level of ∼ 𝜎 . supported, with little to no rotation, and we expect them to have areasonably constant systemic velocity (and velocity dispersion) asa function of radius. In fig. 4 we show how the systemic velocity(top) and velocity dispersion (bottom) vary as a function of radiusin the system. We construct this figure by arranging our stars into5 equal-sized bins, each with 19 stars. We only include stars with 𝑃 member > .
1. In each bin, we compute 𝑣 𝑟 and 𝜎 𝑣 using the sametechnique as before (equation 6). What is immediately striking isthat there is a non-symmetric variation of both these quantities withrespect to the global values (dashed lines) derived in § 3.1. From thevelocity profile, we note that the 𝑣 𝑟 determined for the first and thirdbins (at ∼ . > 𝜎 . Such fluctuationsmay indicate that this is not a system in dynamical equilibrium.Equally striking is the behaviour of the velocity dispersion asa function of radius. As dSph galaxies are thought to be deeply em-bedded within extended dark matter halos, one expects the velocitydispersion of the galaxy to trace this component. As such, it shouldremain constant as a function of radius out to many half-light radii.However, we have already seen that the dispersion in the outskirtsof And XXI is significantly lower than the central value, and this isfurther demonstrated here. What we see is that the majority of thestellar population has a dispersion consistent with ∼ − − ,with the exception of one bin at ∼ (cid:48) ( ∼
480 pc), which appears sig-nificantly dynamically hotter. This can even be seen in the raw data shown in fig. 2. Is this some signature of substructure, or merely anartefact of small number statistics?Such bumps and wiggles in the dynamical profiles of dSphshave previously been observed (e.g. Andromeda II and XIX, Amor-isco et al. 2014; Collins et al. 2020). These have been interpretedas evidence of mergers, tidal effects, substructure, or just randomfluctuations in the data. To investigate how often one expects to seesuch variation about the mean, we run tests on mock data.First, we measure the chi-square goodness of fit our velocityand dispersion data (where we compare each bin to the mean of thefull sample). We measure this to be 𝜒 𝑣 = . 𝜒 𝑠 = . 𝑣 𝑟 and 𝜎 𝑣 for each realisation with our observations. We findthat 38/1000 realisations have 𝜒 𝑣 ≥ .
0, 35/1000 realisations have 𝜒 𝑠 > .
4, and 13/1000 cases satisfy both criteria. We show examplesof a random selection of ‘normal’ (better chi-square) and ‘outlier’(higher chi-square) realisations in fig. 5. Statistically, our observeddistribution can be drawn from a dwarf galaxy with a flat, Gaussianprofile ∼ .
3% of the time (a 2.5 𝜎 event). As such, it is an interestinganomaly, but not strong evidence for substructure. More data wouldresolve whether this is a true feature, or merely noise in the data.Finally, we construct major and minor axis profiles of And XXI(shown in fig. 6) to see whether there are any signs of rotation aboutthese axes. Both appear flat, suggesting there is no (significant)rotation in And XXI. As a final check, we use the same MCMCtechnique from Collins et al. (2020) to search for a signature ofrotation along any arbitrary axis, however no such signal is found. We measure the metallicity of our And XXI member stars usingthe equivalent widths of the three Ca II triplet absorption lines.These features are well-known to allow a proxy for iron abundancemeasurements, [Fe/H], in RGB stars (e.g. Armandroff & Da Costa1991). We perform this for all stars with 𝑆 / 𝑁 > 𝑆 / 𝑁 cut, and that have metallicity uncertainties of less than 0.8 dexin fig. 7. We see that the distribution is centred around [ Fe / H ] = − . ± .
1. We cannot resolve a metallicity spread given the largeuncertainties in our measurements, measuring 𝜎 [ Fe / H ] = . ± . 𝜎 [ Fe / H ] < . MNRAS000
1. We cannot resolve a metallicity spread given the largeuncertainties in our measurements, measuring 𝜎 [ Fe / H ] = . ± . 𝜎 [ Fe / H ] < . MNRAS000 , 1–16 (2020) ndromeda XXI – a detailed study v e l ( k m s ) And XXI - mock profile, normal
Radius (arcmin) V e l o c i t y d i s p e r s i o n ( k m / s ) v e l ( k m s ) And XXI - mock profile, outlier
Radius (arcmin) V e l o c i t y d i s p e r s i o n ( k m / s ) Figure 5. Left:
Radial velocity profile (top) and velocity dispersion profile (bottom) for 3 random realisations of our mock data (shown with red circles, greensquares and blue triangles). These profiles all have chi-square values lower than our observed data, and so show ‘normal’ flat velocity and dispersion profiles.
Right:
As left, but for mock realisations where the chi-square values for both systemic velocity and velocity dispersion are similar to our observed data. Wesee very similar behaviour, where most of the galaxy has a dispersion of 3 − − , but 1-2 bins show dispersions of order 10 km s − . This is seen in ∼ . v e l ( k m / s ) And XXI - major axis profile
Radius (arcmin) V e l o c it y d i s p e r s i on ( k m / s ) v e l ( k m / s ) And XXI - minor axis profile
Radius (arcmin) V e l o c it y d i s p e r s i on ( k m / s ) Figure 6.
Velocity profiles and dispersion profiles along the photometric major (left) and minor (right) axis. Neither the major nor minor axis shows evidencefor rotation or a velocity gradient, and there is significant scatter around both the systemic velocity and average dispersion. the Milky Way as blue triangles, Kirby et al. 2013; Tollerud et al.2012; Collins et al. 2013, 2017, 2020; Kirby et al. 2013; Martinet al. 2013; Ho et al. 2015; Wojno et al. 2020). The best fit to thisrelation from Kirby et al. (2013) is shown as a dashed black line,with the grey band representing the 1 𝜎 scatter. Using our derived velocity dispersion, we can measure the centralmass, density and mass-to-light ratio of And XXI, and compare itwith similar dwarf galaxies in the Local Group. It has been well-established that one can constrain the mass within a given radius, 𝑅 , for dwarf spheroidals using their velocity dispersions. Typically, 𝑅 is similar or equal to the projected half light radius of the galaxy MNRAS , 1–16 (2020)
M. L. M. Collins et al. [Fe/H] N * [Fe/H] = +0.100.09 . . . . [Fe/H] . . . . [ F e / H ] . . . . [Fe/H] [Fe/H] = +0.120.08 Figure 7.
The metallicity distribution function for And XXI member stars with 𝑆 / 𝑁 > = − . ± . 𝜎 [ Fe / H ] < . Luminosity ( L ) [ F e / H ] Figure 8.
The luminosity-metallicity relation for Local Group dSphs.And XXI is shown as a red star. The dSphs of M31 are represented with darkgrey circles, while those of the Milky Way are light grey triangles. The blackdashed line is the best fit relation for Milky Way dSphs from Kirby et al.(2013), with the grey band showing the 1 𝜎 scatter. And XXI is perfectlyconsistent with this relation. (Walker et al. 2009; Wolf et al. 2010). Recent work by Errani et al.(2018) use the mass within 𝑅 = . 𝑟 half , and no assumption on theshape of the mass (dispersion) profile, such that: 𝑀 ( 𝑅 < . 𝑟 half ) = . 𝑟 half 𝜎 𝑣 𝐺 . (8) Using this, we measure a central mass of 𝑀 ( 𝑅 < . 𝑟 half ) = . + . − . × M (cid:12) . This gives a central mass-to-light ratio of [ 𝑀 / 𝐿 ]( 𝑅 < . 𝑟 half ) = + − M (cid:12) / L (cid:12) , significantly dark mat-ter dominated. We measure a central dark matter density of 𝜌 ( 𝑅 < . 𝑟 half ) = . ± . × − M (cid:12) pc − . In our previouswork, we found that And XXI had a lower central mass and densitythan dwarf galaxies of a comparable size or brightness. With ourupdated measurements, we find this remains true. In fig. 9, we showthe mass and density within 1 . 𝑟 half for Milky Way (blue trian-gles) and M31 (red circles) dSphs, compared to the best fit NFWmass(/density) profiles for this population from Collins et al. (2014)in grey (Walker et al. 2007, 2009; Simon & Geha 2007; Simon et al.2011, 2015; Martin et al. 2007, 2013, 2014; Ho et al. 2012; Tollerudet al. 2012; Collins et al. 2013, 2020; Kirby et al. 2015, 2017). Inboth panels, we see that And XXI is of significantly lower massand density than the best fit profile would predict. Could it be lowerdensity because it harbours a central core, as is supposed for somelow-density Milky Way dSphs (e.g. Fornax, Crater 2 and Antlia 2(Goerdt et al. 2006; Walker & Peñarrubia 2011; Amorisco & Evans2012; Read et al. 2019; Caldwell et al. 2017; Torrealba et al. 2019)?And is this growing number of low density dwarf galaxies a chal-lenge for the cold dark matter paradigm? To test this, we use Jeansmodelling to measure the density profile of And XXI, and compareto Λ CDM expectations.
GravSphereWe use an updated version of the GravSphere Jeans modellingcode to model our dynamical and photometric data for And XXI, A version of the GravSphere code with the free form massmodel, amongst others (Genina et al. 2019), is available to down-load from https://github.com/AnnaGenina/pyGravSphere. The updatedMNRAS000
GravSphereWe use an updated version of the GravSphere Jeans modellingcode to model our dynamical and photometric data for And XXI, A version of the GravSphere code with the free form massmodel, amongst others (Genina et al. 2019), is available to down-load from https://github.com/AnnaGenina/pyGravSphere. The updatedMNRAS000 , 1–16 (2020) ndromeda XXI – a detailed study
100 1000 r (pc) M ( r < . × r h a l f ) ( M )
100 1000 r (pc) ( r < . × r h a l f ) ( M / p c ) Figure 9. Left:
Mass enclosed within 1 . × 𝑟 half for the Local Group dSphs. The MW subsystem are shown with light grey triangles, while the M31 systemare shown as dark grey circles. And XXI is highlighted as a red star. The black line shows the best fit NFW halo profile for the Local Group dSphs from Collinset al. (2014), with the gray bad representing the 1 𝜎 scatter, And XXI is significantly below this fit, indicating it has a lower mass than expected for its size. Right:
Enclosed density for the Local Group dSphs, plus the best fit density profile, as in the left-hand panel. And XXI is lower density than typically expected. with the end goal of measuring its dark matter density profile, 𝜌 ( 𝑟 ) .GravSphere is described in detail in Read & Steger (2017) andRead et al. (2018); here we briefly summarise its implementation.GravSphere solves the spherical Jean equation Jeans (1922) for aset of ‘tracers’ (i.e. our member stars with radial velocity measure-ments), to determine both 𝜌 ( 𝑟 ) , and the velocity anisotropy profile, 𝛽 ( 𝑟 ) . The Jeans equation is given by (van der Marel 1994; Mamon& Łokas 2005): 𝜎 ( 𝑅 ) = Σ ( 𝑅 ) ∫ ∞ 𝑅 (cid:18) − 𝛽 ( 𝑟 ) 𝑅 𝑟 (cid:19) 𝜈 ( 𝑟 ) 𝜎 𝑟 𝑟 d 𝑟 √ 𝑟 − 𝑅 , (9)where Σ ( 𝑅 ) is the surface brightness profile at projected radius 𝑅 , 𝜈 ( 𝑟 ) is the spherically averaged tracer density as a function ofspherical radius, 𝑟 , and 𝛽 ( 𝑟 ) is the velocity anisotropy: 𝛽 = − 𝜎 𝑡 𝜎 𝑟 . (10)where 𝜎 𝑟 and 𝜎 𝑡 are the radial and tangential velocity dispersionprofiles, respectively, and 𝜎 𝑟 is given by: 𝜎 𝑟 ( 𝑟 ) = 𝜈 ( 𝑟 ) 𝑔 ( 𝑟 ) ∫ ∞ 𝑟 𝐺 𝑀 ( ˜ 𝑟 ) 𝜈 ( ˜ 𝑟 ) ˜ 𝑟 𝑔 ( ˜ 𝑟 ) d 𝑟, (11)where: 𝑔 ( 𝑟 ) = exp (cid:18) ∫ 𝛽 ( 𝑟 ) 𝑟 d 𝑟 (cid:19) , (12)and 𝑀 ( 𝑟 ) is the cumulative mass profile of the system. GravSphere code described in this paper, along with the new bin-ulator binning method (§4.1.1), is available to download fromhttps://github.com/justinread/gravsphere.
The light profile, Σ ( 𝑅 ) , is modelled as a superposition of Plum-mer spheres (Plummer 1911; Rojas-Niño et al. 2016). The cumula-tive mass profile is given by: 𝑀 ( 𝑟 ) = 𝑀 ∗ ( 𝑟 ) + 𝑀 cNFWt ( 𝑟 ) , (13)where the cumulative stellar light profile, 𝑀 ∗ ( 𝑟 ) , is normalised toasymptote to the total stellar mass at infinity, and allowed to varywithin some flat prior range (see §4.1.2), and 𝑀 cNFWt ( < 𝑟 ) is thecoreNFWtides profile from Read et al. (2018) that describes thecumulative dark matter mass profile: 𝑀 cNFWt ( 𝑟 ) = 𝑀 cNFW ( < 𝑟 ) 𝑟 < 𝑟 𝑡 𝑀 cNFW ( 𝑟 𝑡 ) + 𝜋𝜌 cNFW ( 𝑟 𝑡 ) 𝑟 𝑡 − 𝛿 (cid:20)(cid:16) 𝑟𝑟 𝑡 (cid:17) − 𝛿 − (cid:21) 𝑟 > 𝑟 𝑡 (14)where: 𝑀 cNFW ( < 𝑟 ) = 𝑀 NFW ( < 𝑟 ) 𝑓 𝑛 (15)and: 𝑀 NFW ( 𝑟 ) = 𝑀 𝑔 𝑐 (cid:34) ln (cid:18) + 𝑟𝑟 𝑠 (cid:19) − 𝑟𝑟 𝑠 (cid:18) + 𝑟𝑟 𝑠 (cid:19) − (cid:35) (16)with: 𝑓 𝑛 = (cid:20) tanh (cid:18) 𝑟𝑟 𝑐 (cid:19)(cid:21) 𝑛 (17)and 𝜌 cNFW is given by: 𝜌 cNFW ( 𝑟 ) = 𝑓 𝑛 𝜌 NFW + 𝑛 𝑓 𝑛 − ( − 𝑓 ) 𝜋𝑟 𝑟 𝑐 𝑀 NFW (18)
MNRAS , 1–16 (2020) M. L. M. Collins et al. where: 𝜌 NFW ( 𝑟 ) = 𝜌 (cid:18) 𝑟𝑟 𝑠 (cid:19) − (cid:18) + 𝑟𝑟 𝑠 (cid:19) − (19)and the central density 𝜌 and scale length 𝑟 𝑠 are given by: 𝜌 = 𝜌 crit Δ 𝑐 𝑔 𝑐 / 𝑟 𝑠 = 𝑟 / 𝑐 (20) 𝑔 𝑐 = ( + 𝑐 ) − 𝑐 + 𝑐 (21)and 𝑟 = (cid:20) 𝑀 𝜋 Δ 𝜌 crit (cid:21) / (22)where 𝑐 is the dimensionless concentration parameter ; Δ = 𝜌 crit = .
05 M (cid:12) kpc − is the crit-ical density of the Universe at redshift 𝑧 = 𝑟 is the ‘virial’radius at which the mean enclosed density is Δ × 𝜌 crit ; and 𝑀 isthe ‘virial’ mass – the mass within 𝑟 . The coreNFWtides profileadds to these a parameter 𝑛 that determines how centrally cuspedthe density is ( 𝑛 = 𝑛 = 𝜌 ∝ 𝑟 − cusp), 𝑟 𝑐 that sets the size of this core, and 𝑟 𝑡 thatdetermines an outer ‘tidal radius’ beyond which the density fall-offsteepens as 𝜌 ∝ 𝑟 − 𝛿 .The velocity anisotropy profile is given by: 𝛽 ( 𝑟 ) = 𝛽 + ( 𝛽 ∞ − 𝛽 ) + (cid:0) 𝑟 𝑟 (cid:1) 𝑞 (23)which makes the solution to equation 12 analytic.GravSphere fits the surface brightness profile, Σ ( 𝑅 ) , and ra-dial velocity dispersion profile, 𝜎 los , using the Markov chain MonteCarlo code Emcee (Foreman-Mackey et al. 2013). A symmeterisedversion of 𝛽 ( 𝑟 ) is used to avoid infinities, defined as:˜ 𝛽 = 𝜎 𝑟 − 𝜎 𝑡 𝜎 𝑟 + 𝜎 𝑡 = 𝛽 − 𝛽 , (24)where ˜ 𝛽 = 𝛽 = − 𝛽 = 𝜌 and 𝛽 (Merrifield & Kent 1990;Richardson & Fairbairn 2014; Read & Steger 2017).GravSphere has been extensively tested on mock data forspherical and triaxial systems (Read et al. 2017, 2021), for tidallydisrupting mocks (Read et al. 2018) and for realistic mocks drawnfrom a cosmological simulation (Genina et al. 2019). In most cases,GravSphere is able to recover the dark matter density profile withinits 95% confidence intervals over the range 0 . < 𝑅 / 𝑅 / < 𝑅 / is the half light radius. However, the code has a tendencyto underestimate the density at large radii 𝑅 (cid:38) 𝑅 / (Read et al.2017, 2021). Furthermore, for small numbers of stars, or where themeasurement error on each star is large, the binning method withinthe code can become biased (Gregory et al. 2019; Zoutendijk et al.2021).To address the above concerns, for this paper we have exten-sively updated and improved GravSphere. The primary change isa new data binning module, binulator, that we describe in §4.1.1.Smaller additional changes are: 1) a switch to using the coreN-FWtides profile as the default dark matter mass model, as in Readet al. (2018); Read & Erkal (2019), rather than using the ‘non-parametric’ series of power laws centred on a set of radial bins, described in Read & Steger (2017); and 2) using tighter priors on 𝛽 ( 𝑟 ) by default. We describe all of our priors in §4.1.2 and presentmock data tests of our new methodology in §A.The reasons for shifting away from the ‘non-parametric’ darkmatter mass profile from Read & Steger (2017) is twofold. Firstly, for ∼ + tracers, the results from using the coreNFWtides profileand the ‘non-parametric’ profile agree within their respective 68%confidence intervals Alvarez et al. (2020), yet the coreNFWtidesmodel fit parameters are more cosmologically useful (Read et al.2018; Read & Erkal 2019). Secondly, in the absence of data, thedefault priors on the non-parametric mass profile favour a highdensity in the centre that falls very steeply outwards, inconsistentwith expectations in most popular cosmological models. By con-trast, the coreNFWtides model naturally defaults to cosmologicalexpectations at large radii in the absence of data.The reason for the tighter priors on ˜ 𝛽 ( 𝑟 ) is that, theoretically,we expect that dynamical systems in pseudo-equilibrium shouldbe close to isotropic in the centre, with radial anisotropy, or weaktangential anisotropy, at large radii (e.g. Pontzen et al. 2015; Geninaet al. 2019; Alvey et al. 2021). For this reason, we allow for onlymild tangential anisotropy ( ˜ 𝛽 ( 𝑟 ) > − .
1) and demand ˜ 𝛽 ( 𝑟 ) → 𝑟 → binulatorThe main improvement we make to GravSphere is a completereworking of its data binning routines into a separate code: the bin-ulator. This improves the binning by fitting a generalised Gaussianprobability distribution function (PDF) to each bin to estimate itsmean, variance and kurtosis. This has the advantages that: 1) thedistribution function can be readily convolved with the error PDF ofeach star, survey selection functions, binary star velocity PDFs, andsimilar; and 2) the method returns a robust estimate of the mean,variance and kurtosis, and their uncertainties, even in the limit of avery small number of tracer stars.The full method proceeds, as follows. Firstly, the data for Σ ( 𝑅 ) are fit using Emcee (Foreman-Mackey et al. 2013) to obtain thebest-fit multi-Plummer model for the light profile (see above). Thisprovides an initial guess for the later GravSphere fits, and will beimportant for calculating the virial shape parameters (see below). Toprovide maximum flexibility, these fits allow for individual Plummercomponents to have ‘negative mass’ while still ensuring that the totalsurface density is positive definite, as in Rojas-Niño et al. (2016).Next, the discreet stellar velocity data are sorted into equal numberbins in radius, weighted by membership probability. We use 25stars per bin by default; our results are not sensitive to this choice.A velocity PDF is then fit to the stars in each bin to determinethe mean, variance and kurtosis of the bin, similarly to the methoddescribed in Sanders & Evans (2020). However, for our velocityPDF, we assume a generalised Gaussian: 𝑝 𝑖 = 𝛽 𝑣 𝛼 𝑣 Γ ( . / 𝛽 𝑣 ) exp (cid:16) −(| 𝑣 los ,𝑖 − 𝜇 𝑣 |/ 𝛼 𝑣 ) 𝛽𝑣 ) (cid:17) (25)where 𝑣 los ,𝑖 is the line of sight velocity of a star, 𝑖 , Γ ( 𝑥 ) is the GammaFunction, and 𝜇 𝑣 , 𝛼 𝑣 and 𝛽 𝑣 are parameters fit to each bin that relateto moments of the velocity distribution function. The mean is givenby 𝜇 𝑣 ; the variance by 𝜎 = 𝛼 𝑣 Γ ( . / 𝛽 𝑣 )/ Γ ( . / 𝛽 𝑣 ) ; and thekurtosis by 𝜅 = Γ ( . / 𝛽 𝑣 ) Γ ( . / 𝛽 𝑣 )/[ Γ ( . / 𝛽 𝑣 )] .To account for errors on the velocity of each star, the abovegeneralised Gaussian should each be convolved with the error prob-ability distribution function (PDF) for each star. This convolution MNRAS000
1) and demand ˜ 𝛽 ( 𝑟 ) → 𝑟 → binulatorThe main improvement we make to GravSphere is a completereworking of its data binning routines into a separate code: the bin-ulator. This improves the binning by fitting a generalised Gaussianprobability distribution function (PDF) to each bin to estimate itsmean, variance and kurtosis. This has the advantages that: 1) thedistribution function can be readily convolved with the error PDF ofeach star, survey selection functions, binary star velocity PDFs, andsimilar; and 2) the method returns a robust estimate of the mean,variance and kurtosis, and their uncertainties, even in the limit of avery small number of tracer stars.The full method proceeds, as follows. Firstly, the data for Σ ( 𝑅 ) are fit using Emcee (Foreman-Mackey et al. 2013) to obtain thebest-fit multi-Plummer model for the light profile (see above). Thisprovides an initial guess for the later GravSphere fits, and will beimportant for calculating the virial shape parameters (see below). Toprovide maximum flexibility, these fits allow for individual Plummercomponents to have ‘negative mass’ while still ensuring that the totalsurface density is positive definite, as in Rojas-Niño et al. (2016).Next, the discreet stellar velocity data are sorted into equal numberbins in radius, weighted by membership probability. We use 25stars per bin by default; our results are not sensitive to this choice.A velocity PDF is then fit to the stars in each bin to determinethe mean, variance and kurtosis of the bin, similarly to the methoddescribed in Sanders & Evans (2020). However, for our velocityPDF, we assume a generalised Gaussian: 𝑝 𝑖 = 𝛽 𝑣 𝛼 𝑣 Γ ( . / 𝛽 𝑣 ) exp (cid:16) −(| 𝑣 los ,𝑖 − 𝜇 𝑣 |/ 𝛼 𝑣 ) 𝛽𝑣 ) (cid:17) (25)where 𝑣 los ,𝑖 is the line of sight velocity of a star, 𝑖 , Γ ( 𝑥 ) is the GammaFunction, and 𝜇 𝑣 , 𝛼 𝑣 and 𝛽 𝑣 are parameters fit to each bin that relateto moments of the velocity distribution function. The mean is givenby 𝜇 𝑣 ; the variance by 𝜎 = 𝛼 𝑣 Γ ( . / 𝛽 𝑣 )/ Γ ( . / 𝛽 𝑣 ) ; and thekurtosis by 𝜅 = Γ ( . / 𝛽 𝑣 ) Γ ( . / 𝛽 𝑣 )/[ Γ ( . / 𝛽 𝑣 )] .To account for errors on the velocity of each star, the abovegeneralised Gaussian should each be convolved with the error prob-ability distribution function (PDF) for each star. This convolution MNRAS000 , 1–16 (2020) ndromeda XXI – a detailed study − −
20 0 20 400 . . f r e q u e n c y β v , κ = 1 . , . β v , κ = 2 . , . β v , κ = 5 . , . − −
20 0 20 40v los (km/s) − − R e s i du a l ( % ) Figure 10.
The generalised Gaussian velocity PDF (equation 25), convolvedwith a Gaussian error PDF of width 𝜎 e = 𝛼 𝑣 =
15 km/s, with varying 𝛽 𝑣 , as marked. The correspondingkurtosis, 𝜅 , for each model is marked in the legend. Notice that 𝛽 𝑣 = ( PDFtrue − PDFfast )/ max [ PDFtrue ] × 𝜌 − 𝛽 degeneracy. integral, however, is expensive to compute. To speed up the cal-culation, we employ an analytic approximation to this convolutionintegral that assumes Gaussian errors for the stellar velocities andis exact in the limits that: 1) the normalised Gaussian approaches aGaussian ( 𝛽 𝑣 = 𝛼 → ˜ 𝛼 = 𝛼 + 𝜎 ,𝑖 Γ ( . / 𝛽 𝑣 )/ Γ ( . / 𝛽 𝑣 ) (26)where 𝜎 e ,𝑖 is the Gaussian width of the error PDF of star, 𝑖 . Thequality of this approximation is shown in Figure 10, along with someexample generalised Gaussian PDFs. Notice that for this example(typical of current data for nearby dwarf spheroidals), the error inthis approximation is typically less than 5%, and everywhere lessthan 10%.We fit the above normalised Gaussian to each bin with Emceeusing the likelihood: L = 𝑁 (cid:214) 𝑖 𝑝 𝑖 (27)where 𝑁 is the (weighted) number of stars in each bin.We assume flat priors on the parameters: − . < 𝜇 𝑣 / kms − < .
1; 1 < 𝛼 𝑣 / kms − <
25; and 1 < 𝛽 𝑣 <
5. This allows for akurtosis in the range 2 < 𝜅 < (cid:104) 𝑣 los (cid:105) , 𝜎 los and 𝜅 for each bin. From these, we can also calculate the marginalised PDF of fourth velocity moments (cid:104) 𝑣 (cid:105) = 𝜅𝜎 . We use this, alongwith the fit to Σ ( 𝑅 ) , to determine the marginalised PDFs of the firstand second virial shape parameters: 𝑣 𝑠 = ∫ ∞ Σ (cid:104) 𝑣 (cid:105) 𝑅 d 𝑅 (28)and 𝑣 𝑠 = ∫ ∞ Σ (cid:104) 𝑣 (cid:105) 𝑅 d 𝑅 . (29)where the above integrals are calculated numerically for 2,500 ran-dom draws from the marginalised distribution of (cid:104) 𝑣 (cid:105) in each bin,assuming that (cid:104) 𝑣 (cid:105) is constant beyond the outermost bin. We findthat the marginalised PDFs for 𝜎 los for each bin are close to Gaus-sian. For this reason, GravSphere assumes Gaussian uncertaintieson 𝜎 los when performing its fit to 𝜎 los ( 𝑅 ) . However, the PDFs for 𝑣 𝑠 and 𝑣 𝑠 are typically non-Gaussian and so a final improvementis that GravSphere now incorporates these non-Gaussian PDFs for 𝑣 𝑠 and 𝑣 𝑠 self-consistently into its likelihood function.Since we now fit a velocity PDF to each bin, the binulatorwill work equally well with very low numbers of tracers stars and/orwhen the velocity errors are large, resolving the issues with earlierversions of GravSphere, outlined above. We present tests of thisupdated binulator+GravSphere code on mock data in AppendixA. We use priors on the coreNFWtides model parameters of: 7 . < log ( 𝑀 / M (cid:12) ) < .
5; 7 < 𝑐 < − < log ( 𝑟 𝑐 / kpc ) <
1; 1 < log ( 𝑟 𝑡 / kpc ) <
20; 3 < 𝛿 <
5; and 0 < 𝑛 <
1. Weuse priors on the symmetrised velocity anisotropy of: ˜ 𝛽 ( ) = − . < ˜ 𝛽 ∞ < − < log ( 𝑟 / kpc ) <
0; and 1 < 𝑞 <
3. Finally,we use a flat prior on the stellar mass of 5 . < 𝑀 ∗ /( 𝑀 (cid:12) ) < . GravSphere modelling of And XXI
We now apply binulator+GravSphere to And XXI, with the goalof constraining its dark matter density profile. We take our surfacebrightness profile from our Subaru Suprime-cam imaging, whichcovers And XXI out to 4 effective radii (this agrees within quoteduncertainties with that derived from PAndAS photometry (Martinet al. 2016)), and our dispersion profile is constructed from our77 spectroscopically identified members. To generate this profile,the probabilities of our likely members are summed to give an“effective” number of members, 𝑁 eff = 𝑁 mem ∑︁ 𝑖 = 𝑃 mem , i . (30)For our sample, 𝑁 eff =
50. We bin the data radially from thecentre of And XXI in two bins, using an effective 25 stars per bin. Weshow the Σ ∗ ( 𝑅 ) and 𝜎 LOS ( 𝑅 ) profiles in fig. 11. The vertical blueline represents the half-light radius determined from the PAndASimaging. We model And XXI under the assumption that it is aspherical, non-rotating system, and we use the stellar mass fromMcConnachie (2012) of 𝑀 ∗ = . × M (cid:12) , assuming an error of25%.Using the above, we model the density profile of And XXIand show the result in the lower left panel of fig. 11. We measure a MNRAS , 1–16 (2020) M. L. M. Collins et al. − R [kpc]10 − − − − − Σ ∗ [ N k p c − ] − . − . . . . [ R/ kpc]02468101214 σ L O S [ k m s − ] − . − . . . . [ r/ kpc] − . − . − . − . . . . . . ˜ β − r [kpc]10 ρ [ M (cid:12) k p c − ] FitAbund. match (SFR)Abund. match ( M ∗ ) Figure 11. Top left:
Surface brightness profile for And XXI. The blue points show photometric data from deep Subaru imaging. The best-fit from GravSPhereis shown as a solid black line, with the gray shaded regions showing the 68% and and 95% confidence intervals.
Top right:
The velocity dispersion profilefor And XXI. Each bin contains an effective ∼
25 stars weighted by membership probability. Again, the best fit and confidence intervals from GravSphere areshown in black/grey.
Bottom left:
The symmetrised radial anisotropy profile for And XXI.
Bottom right:
The final dark matter density profile inferred byGravSphere. Notice that it is consistent with both a cusp and core within GravSphere’s 95% confidence intervals (light grey). However, it is less dense at allradii than expected from abundance matching in Λ CDM, whether abundace matching with the stellar mass (blue) or the mean star formation rate (red). In allpanels, the blue vertical line shows the half-light radius for And XXI determined from PAndAS photometry (Martin et al. 2016). central dark matter density within 150 pc of 𝜌 DM ( ) = . ± . × M (cid:12) kpc − at 68% confidence. The data are consistentwith a cusp or core within 2 𝜎 (light gray shading). Overlaid onthe Figure are two density profiles determined from abundancematching in Λ CDM, using the stellar mass (blue) and the mean starformation rate (red), as in Read & Erkal (2019) (a more detaileddescription of how these are derived is presented in §5.1, below).For the latter, we use the star formation history for And XXI fromWeisz et al. (2019b). Notice that in both cases, our determinationof the dark matter density profile for And XXI is lower at all radii than expectations from abundance matching in Λ CDM by a factorof ∼ − Λ CDM context, and investigate whether dark matterheating from star formation or tidal processes could explain thisresult.
Star formation can lower the central density of a dark matter halo,gradually turning a dark cusp into a core within ∼ the half lightradius of the stars, 𝑅 / (Navarro et al. 1996b; Read & Gilmore2005; Pontzen & Governato 2012; Read et al. 2016, 2019). Thestar formation history of And XXI has been measured from shallow MNRAS000
Star formation can lower the central density of a dark matter halo,gradually turning a dark cusp into a core within ∼ the half lightradius of the stars, 𝑅 / (Navarro et al. 1996b; Read & Gilmore2005; Pontzen & Governato 2012; Read et al. 2016, 2019). Thestar formation history of And XXI has been measured from shallow MNRAS000 , 1–16 (2020) ndromeda XXI – a detailed study HST data (reaching ∼ ∼
50% of its stars prior to 8.3 Gyr ago, and90% prior to 5.8 Gyr ago. There is no sign of any star formation inthe last 2 − 𝑀 for And XXI. For this we use the abundancematching machinery from Read & Erkal (2019). Classic abundancematching estimates 𝑀 statistically from the stellar mass, 𝑀 ∗ .For And XXI, using the abundance matching relation from Readet al. (2017) (consistent with Behroozi et al. 2013), this yields 𝑀 = . ± . × M (cid:12) . However, for satellite galaxies likeAnd XXI that quench on infall to a larger host, 𝑀 ∗ becomesa poor proxy for the pre-infall 𝑀 (e.g. Ural et al. 2015; To-mozeiu et al. 2016; Read & Erkal 2019). Read & Erkal (2019)argue that using instead the mean star formation rate, (cid:104) SFR (cid:105) ,evaluated over the period during which the satellite was formingstars, solves this ‘quenching problem’. They show empirically forMilky Way satellites that (cid:104)
SFR (cid:105) correlates much better with 𝑀 than does 𝑀 ∗ . If we abundance match And XXI using instead its (cid:104) SFR (cid:105) = . ± . × − M (cid:12) / yr (calculated from And XXI’s starformation history as in Read & Erkal 2019), we obtain a higher 𝑀 = . ± . × M (cid:12) .Next, we must calculate whether or not there has been sufficientstar formation within And XXI to cause its central density to be low-ered significantly. For this, we use the coreNFW profile calibratedon the simulations from Read et al. 2016. Read et al. (2016) deter-mine the amount of core formation based on the total ‘star formationtime’, defined as in Read & Erkal (2019) as 𝑡 SF = 𝑀 ∗ /(cid:104) SFR (cid:105) . ForAnd XXI, this gives 𝑡 SF = . 𝑛 = tanh ( . 𝑡 SF / 𝑡 dyn ) , where 𝑡 dyn is the dynamical timeat the scale radius, 𝑟 𝑠 (equation 16; and see Read et al. 2016). (Re-call that 𝑛 = 𝑛 = 𝑟 − density cusp.)Assuming the Dutton & Macciò (2014) 𝑀 − 𝑐 relationfor Λ CDM, and its scatter, we use the above coreNFW profile toestimate the range of theoretical expectations for And XXI’s densityprofile in Λ CDM. This is shown in Figure 11, bottom right panel.The blue band shows the expected range assuming 𝑀 derivedfrom abundance matching with 𝑀 ∗ ; the red band shows the sameusing abundance matching with (cid:104) SFR (cid:105) . In both cases, the width ofthe band incorporates uncertainties in 𝑀 and the expected 2 − 𝜎 scatter in 𝑐 . Notice that in all models, the central density hasbeen slightly lowered by star formation, but not enough to producea flat core (we find that our models span the range 𝑛 = . − . 𝑅 > 𝑅 / .The end result of this exercise is that And XXI’s density is lowerthan expected for isolated halos in Λ CDM by a factor of ∼ − 𝑅 / = .
875 kpc) where the density is best-constrained. And, itremains true when marginalising over the expected 2 − 𝜎 scatterin 𝑐 in Λ CDM, and when accounting for dark matter heatinglowering the central dark matter density. We conclude, therefore,that And XXI’s density at all radii is lower than expected for isolatedhalos in Λ CDM. We discuss next whether this can be explained bytides.
In the discussion above, we have assumed that And XXI is unaffectedby tides. Tidal stripping, and even more so tidal shocking, will actto lower And XXI’s density over time (e.g. Gnedin et al. 1999;Read et al. 2006b,a; Amorisco 2019). If we assume that significanttidal stripping (losing 90-99% of its original mass) has taken place,we can use the relations for tidal evolution from Peñarrubia et al.(2010b) to assess the impact on And XXI’s density. These tidaltracks can be implemented for either cuspy or cored dark matterhalos, and the effects for these are very different. If And XXI possessan NFW cusp ( 𝛾 = 𝛾 = .
5, and removal of 95% of theoriginal halo, we can find central density values that are consistentwith our measurements for And XXI within their 68% confidenceintervals.However, tidal shocking can be much more efficient than tidalstripping. The effect is maximised if And XXI is on a plunging orbitand starts out with low density, either due to some inner cusp-coretransformation, or due to it inhabiting a low concentration halo (seediscussion in §1 and Amorisco 2019).The above results imply that we can explain And XXI’s lowdensity through a combination of dark matter heating and tidal pro-cesses, with the dominant effect coming from tides. But, how likelyis it that And XXI is on an orbit that has allowed it to experienceextreme tidal stripping or shocking? Both tidal scenarios requirea highly radial orbit, with a small pericentre (<20 kpc; Read et al.2006a). And XXI is currently far from its host, at a 3D distance of 𝐷 M31 = + − kpc (Weisz et al. 2019a). Without proper motions,it is difficult for us to place meaningful constraints on the currentorbit of And XXI to determine whether it has recently passed closeto M31. Even with proper motions, its current orbit is not necessar-ily a robust indicator of its past close interactions (e.g. Lux et al.2010; Genina et al. 2020).The current light profile for And XXI shows no obvious signsof tidal stripping in the outskirts, but this does not necessarily implythat tidal stripping of the stellar component has not taken place (e.g.Read et al. 2006a; Peñarrubia et al. 2009; Ural et al. 2015; Geninaet al. 2020). Furthermore, significant tidal shocking can lower thedensity of And XXI if it is on a sufficiently plunging orbit, withoutany tidal stripping of stars taking place (Read et al. 2006a; Amorisco2019). The unusual dynamics of And XXI reported in § 3.2 couldimply that the system is not in dynamical equilibrium, but this 2 𝜎 finding is not adequate to provide unambiguous evidence of tidalstripping. With future data and modelling, better constraints can beplaced on the orbit of And XXI and on its faint stellar outskirts.This will help us understand whether tidal effects have acted tolower the central density of this dwarf galaxy, or whether its lowdensity points to physics beyond Λ CDM.
Finally, we consider other theories which may also explain AndXXI’s low central density and other properties. In McGaugh & Mil-grom (2013) , they calculate the velocity dispersions for AndromedadSphs in the modified Newtonian dynamics (MOND) framework.As dwarf galaxies would be free from dark matter in this framework,one may expect to see a lower central mass or density than predicted
MNRAS , 1–16 (2020) M. L. M. Collins et al. by Λ CDM. Indeed, MOND has been shown to nicely reproduce thelow density cores seen in isolated low surface brightness galaxies(e.g. Sanders & McGaugh 2002; Famaey & McGaugh 2012). Ingeneral, the MOND predictions of McGaugh & Milgrom (2013)for M31 dwarf spheroidals show a good consistency with the mea-sured velocity dispersions presented in both Tollerud et al. (2012)and Collins et al. (2013) (though note that some of the Milky Waydwarfs are more problematic for MOND; Angus et al. e.g. 2014;Read et al. e.g. 2019). Here, we compare our updated dispersion forAnd XXI with these prior predictions.McGaugh & Milgrom (2013) give two estimates for the ve-locity dispersion of M31 dwarfs. The first is the isolated case,and the second is for dwarf galaxies embedded in an externalfield (the external-field effect, EFE). For satellite galaxies, it isnot clear that the isolated approximation is valid, and for And XXIin particular, they recommend using the EFE calculations. Theypredict a velocity dispersion of 𝜎 EFE = . + . − . km s − , wherethe presented value uses an assumption of the stellar mass tolight ratio of [ 𝑀 / 𝐿 ∗ ] = (cid:12) / L (cid:12) , and the uncertainties use [ 𝑀 / 𝐿 ∗ ] = (cid:12) / L (cid:12) . The uncertainty on the predicted dis-persion reflects this range of M/L they considered to be plausible.This value is lower than our measured value of 𝜎 𝑣 = . + . − . km s − ,but it is within 1 𝜎 . As such, the MOND prediction still agrees withour new data. We have presented results from our chemodynamical study of thelow mass M31 satellite, And XXI. We have identified 77 probablestellar members in this galaxy, and our key findings are as follows: • We measure a systemic velocity for And XXI of 𝑣 𝑟 = − . ± . − and a velocity dispersion of 𝜎 𝑣 = . + . − . km s − , con-sistent with the findings of Collins et al. (2013) from a much smallersample. • We find that the stars in the outskirts of And XXI have a farlower velocity dispersion than those within the half-light radius. Wealso see that the systemic velocity and velocity dispersion vary withradius, which could indicate that it is not in dynamical equilibrium. • From those stars with high enough 𝑆 / 𝑁 , we additionally mea-sure the metallicity distribution function of And XXI. We measurea mean metallicity of [ Fe / H ] = − . ± . • We model the dark matter density profile of And XXI usingan updated version of GravSphere, with the main improvementbeing a complete reworking of its data binning routines into a sepa-rate code: the binulator (§4.1.1). This fits a generalised GaussianPDF to each bin to estimate its mean, variance and kurtosis, andtheir uncertainties. This has the advantages that: 1) the distributionfunction can be readily convolved with the error PDF of each star,survey selection functions, binary star velocity PDFs, and similar;and 2) the method returns a robust estimate of the mean, varianceand kurtosis even in the limit of a very small number of tracer stars.We tested this new method on mock data in Appendix A, show-ing that it provides an unbiased estimate of the underlying densityprofile even for ∼
100 tracer stars. • Using binulator+GravSphere, we find a central dark matterdensity for And XXI of 𝜌 DM ( ) = . + . − . × M (cid:12) kpc − at 68% confidence, and a density at two half light radii of 𝜌 DM ( . ) = . + . − . × M (cid:12) kpc − at 68% confidence. Wecannot distinguish between a cusped or cored profile, however thedensity at all radii is a factor ∼ − Λ CDM. We show that this cannot beexplained by ‘dark matter heating’ since And XXI had too little starformation to significantly lower its inner dark matter density, whiledark matter heating only acts on the profile inside the half light ra-dius. However, And XXI’s low density can be accommodated within Λ CDM if it experienced extreme tidal stripping (losing >
95% of itsmass), or if it inhabits a low concentration halo on a plunging orbitthat experienced repeated tidal shocks. Future work establishing theorbit of And XXI would help us understand the origin of its lowinner density and thereby improve our understanding of the natureof dark matter. • When comparing our measured velocity dispersion forAnd XXI with expectations from MOND, we find it to be con-sistent with predictions by McGaugh & Milgrom (2013) within 1 𝜎 .Given this object sits within the low acceleration regime, it willcertainly be of interest for future work. DATA AVAILABILITY
All raw DEIMOS spectra are available via the Keck archive. Anelectronic table with reduced properties (coordinates, magnitudes, 𝑆 / 𝑁 , velocities and metallicities) for all stars will be provided on thejournal website. Fully reduced 1D spectra will be made available onreasonable request to the lead author as these are not hosted on theKeck archive. The data used in the GravSphere modelling are avail-able (with the code) at https://github.com/justinread/gravsphere. ACKNOWLEDGMENTS
RI and NM acknowledge funding from the European ResearchCouncil (ERC) under the European Unions Horizon 2020 researchand innovation programme (grant agreement No. 834148).
REFERENCES
Alvarez A., Calore F., Genina A., Read J., Serpico P. D., Zaldivar B., 2020,JCAP, 2020, 004Alvey J., et al., 2021, MNRAS, 501, 1188Amorisco N. C., 2019, arXiv e-prints, p. arXiv:1901.05460Amorisco N. C., Evans N. W., 2012, MNRAS, 419, 184Amorisco N. C., Evans N. W., van de Ven G., 2014, Nature, 507, 335Angus G. W., Gentile G., Diaferio A., Famaey B., van der Heyden K. J.,2014, MNRAS, 440, 746Armandroff T. E., Da Costa G. S., 1991, AJ, 101, 1329Astropy Collaboration et al., 2018, AJ, 156, 123Behroozi P. S., Wechsler R. H., Conroy C., 2013, ApJ, 770, 57Boylan-Kolchin M., Bullock J. S., Sohn S. T., Besla G., van der Marel R. P.,2013, ApJ, 768, 140Brooks A. M., Zolotov A., 2014, ApJ, 786, 87Caldwell N., et al., 2017, ApJ, 839, 20Chakrabarti S., Chang P., Price-Whelan A. M., Read J., Blitz L., HernquistL., 2019, ApJ, 886, 67Collins M. L. M., et al., 2013, ApJ, 768, 172Collins M. L. M., et al., 2014, ApJ, 783, 7Collins M. L. M., Tollerud E. J., Sand D. J., Bonaca A., Willman B., StraderJ., 2017, MNRAS, 467, 573Collins M. L. M., Tollerud E. J., Rich R. M., Ibata R. A., Martin N. F.,Chapman S. C., Gilbert K. M., Preston J., 2020, MNRAS, 491, 3496Conn A. R., et al., 2012, ApJ, 758, 11 MNRAS000
Alvarez A., Calore F., Genina A., Read J., Serpico P. D., Zaldivar B., 2020,JCAP, 2020, 004Alvey J., et al., 2021, MNRAS, 501, 1188Amorisco N. C., 2019, arXiv e-prints, p. arXiv:1901.05460Amorisco N. C., Evans N. W., 2012, MNRAS, 419, 184Amorisco N. C., Evans N. W., van de Ven G., 2014, Nature, 507, 335Angus G. W., Gentile G., Diaferio A., Famaey B., van der Heyden K. J.,2014, MNRAS, 440, 746Armandroff T. E., Da Costa G. S., 1991, AJ, 101, 1329Astropy Collaboration et al., 2018, AJ, 156, 123Behroozi P. S., Wechsler R. H., Conroy C., 2013, ApJ, 770, 57Boylan-Kolchin M., Bullock J. S., Sohn S. T., Besla G., van der Marel R. P.,2013, ApJ, 768, 140Brooks A. M., Zolotov A., 2014, ApJ, 786, 87Caldwell N., et al., 2017, ApJ, 839, 20Chakrabarti S., Chang P., Price-Whelan A. M., Read J., Blitz L., HernquistL., 2019, ApJ, 886, 67Collins M. L. M., et al., 2013, ApJ, 768, 172Collins M. L. M., et al., 2014, ApJ, 783, 7Collins M. L. M., Tollerud E. J., Sand D. J., Bonaca A., Willman B., StraderJ., 2017, MNRAS, 467, 573Collins M. L. M., Tollerud E. J., Rich R. M., Ibata R. A., Martin N. F.,Chapman S. C., Gilbert K. M., Preston J., 2020, MNRAS, 491, 3496Conn A. R., et al., 2012, ApJ, 758, 11 MNRAS000 , 1–16 (2020) ndromeda XXI – a detailed study Conn A. R., et al., 2013, ApJ, 766, 120Cooper M. C., Newman J. A., Davis M., Finkbeiner D. P., Gerke B. F., 2012,spec2d: DEEP2 DEIMOS Spectral Pipeline, Astrophysics Source CodeLibrary (ascl:1203.003)Di Cintio A., Knebe A., Libeskind N. I., Brook C., Yepes G., Gottloeber S.,Hoffman Y., 2012, preprint, ( arXiv:1204.0515 )Dotter A., Chaboyer B., Jevremović D., Kostov V., Baron E., Ferguson J. W.,2008, ApJS, 178, 89Dutton A. A., Macciò A. V., 2014, MNRAS, 441, 3359Errani R., Peñarrubia J., 2020, MNRAS, 491, 4591Errani R., Peñarrubia J., Laporte C. F. P., Gómez F. A., 2017, MNRAS, 465,L59Errani R., Peñarrubia J., Walker M. G., 2018, MNRAS, 481, 5073Faber S. M., et al., 2003, in Iye M., Moorwood A. F. M., eds, Presented at theSociety of Photo-Optical Instrumentation Engineers (SPIE) ConferenceVol. 4841, Society of Photo-Optical Instrumentation Engineers (SPIE)Conference Series. pp 1657–1669Famaey B., McGaugh S. S., 2012, Living Reviews in Relativity, 15, 10Flores R. A., Primack J. R., 1994, ApJ, 427, L1Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, PASP, 125,306Fritz T. K., Battaglia G., Pawlowski M. S., Kallivayalil N., van der MarelR., Sohn S. T., Brook C., Besla G., 2018, A&A, 619, A103Fu S. W., Simon J. D., Alarcón Jara A. G., 2019, ApJ, 883, 11Genina A., et al., 2019, arXiv e-prints, p. arXiv:1911.09124Genina A., Read J. I., Fattahi A., Frenk C. S., 2020, arXiv e-prints, p.arXiv:2011.09482Gilbert K. M., et al., 2006, ApJ, 652, 1188Gnedin O. Y., Hernquist L., Ostriker J. P., 1999, ApJ, 514, 109Goerdt T., Moore B., Read J. I., Stadel J., Zemp M., 2006, MNRAS, 368,1073Gregory A. L., Collins M. L. M., Read J. I., Irwin M. J., Ibata R. A., MartinN. F., McConnachie A. W., Weisz D. R., 2019, MNRAS, 485, 2010Ho N., et al., 2012, ApJ, 758, 124Ho N., Geha M., Tollerud E. J., Zinn R., Guhathakurta P., Vargas L. C.,2015, ApJ, 798, 77Hunter J. D., 2007, Computing In Science & Engineering, 9, 90Ibata R., Sollima A., Nipoti C., Bellazzini M., Chapman S. C., DalessandroE., 2011, ApJ, 738, 186Irwin M., Lewis J., 2001, New Astronomy Review, 45, 105Jeans J. H., 1922, MNRAS, 82, 122Jones E., Oliphant T., Peterson P., et al., 2001, SciPy: Open source scientifictools for Python,
Kazantzidis S., Mayer L., Mastropietro C., Diemand J., Stadel J., Moore B.,2004, ApJ, 608, 663Kirby E. N., Cohen J. G., Guhathakurta P., Cheng L., Bullock J. S., GallazziA., 2013, ApJ, 779, 102Kirby E. N., Simon J. D., Cohen J. G., 2015, ApJ, 810, 56Kirby E. N., Cohen J. G., Simon J. D., Guhathakurta P., Thygesen A. O.,Duggan G. E., 2017, ApJ, 838, 83Lux H., Read J. I., Lake G., 2010, MNRAS, 406, 2312Mamon G. A., Łokas E. L., 2005, MNRAS, 362, 95Martin N. F., Ibata R. A., Chapman S. C., Irwin M., Lewis G. F., 2007,MNRAS, 380, 281Martin N. F., et al., 2013, ApJ, 772, 15Martin N. F., et al., 2014, ApJ, 793, L14Martin N. F., et al., 2016, ApJ, 833, 167McConnachie A. W., 2012, AJ, 144, 4McConnachie A. W., et al., 2008, ApJ, 688, 1009McGaugh S., Milgrom M., 2013, ApJ, 766, 22Merrifield M. R., Kent S. M., 1990, AJ, 99, 1548Navarro J. F., Eke V. R., Frenk C. S., 1996a, MNRAS, 283, L72Navarro J. F., Eke V. R., Frenk C. S., 1996b, MNRAS, 283, L72Oñorbe J., Boylan-Kolchin M., Bullock J. S., Hopkins P. F., Kereš D.,Faucher-Giguère C.-A., Quataert E., Murray N., 2015, MNRAS, 454,2092Oliphant T., 2006–, NumPy: A guide to NumPy, USA: Trelgol Publishing,
Orkney M. D. A., et al., 2021, arXiv e-prints, p. arXiv:2101.02688Peñarrubia J., Navarro J. F., McConnachie A. W., 2008, ApJ, 673, 226Peñarrubia J., Navarro J. F., McConnachie A. W., Martin N. F., 2009, ApJ,698, 222Peñarrubia J., Benson A. J., Walker M. G., Gilmore G., McConnachie A. W.,Mayer L., 2010a, MNRAS, 406, 1290Peñarrubia J., Belokurov V., Evans N. W., Martínez-Delgado D., GilmoreG., Irwin M., Niederste-Ostholt M., Zucker D. B., 2010b, MNRAS, 408,L26Plummer H. C., 1911, MNRAS, 71, 460Pontzen A., Governato F., 2012, MNRAS, 421, 3464Pontzen A., Read J. I., Teyssier R., Governato F., Gualandris A., Roth N.,Devriendt J., 2015, MNRAS, 451, 1366Read J. I., Erkal D., 2019, MNRAS, 487, 5799Read J. I., Gilmore G., 2005, MNRAS, 356, 107Read J. I., Steger P., 2017, MNRAS, 471, 4541Read J. I., Wilkinson M. I., Evans N. W., Gilmore G., Kleyna J. T., 2006a,MNRAS, 367, 387Read J. I., Pontzen A. P., Viel M., 2006b, MNRAS, 371, 885Read J. I., Agertz O., Collins M. L. M., 2016, MNRAS, 459, 2573Read J. I., Iorio G., Agertz O., Fraternali F., 2017, MNRAS, 467, 2019Read J. I., Walker M. G., Steger P., 2018, MNRAS, 481, 860Read J. I., Walker M. G., Steger P., 2019, MNRAS, 484, 1401Read J. I., et al., 2021, MNRAS, 501, 978Richardson T., Fairbairn M., 2014, MNRAS, 441, 1584Rojas-Niño A., Read J. I., Aguilar L., Delorme M., 2016, MNRAS, 459,3349Sanders J. L., Evans N. W., 2020, MNRAS, 499, 5806Sanders R. H., McGaugh S. S., 2002, ARA&A, 40, 263Sanders J. L., Evans N. W., Dehnen W., 2018, MNRAS, 478, 3879Schlegel D. J., Finkbeiner D. P., Davis M., 1998, ApJ, 500, 525Simon J. D., Geha M., 2007, ApJ, 670, 313Simon J. D., et al., 2011, ApJ, 733, 46Simon J. D., et al., 2015, ApJ, 808, 95Spitzer Lyman J., 1958, ApJ, 127, 17Starkenburg E., et al., 2010, A&A, 513, A34+Tollerud E. J., et al., 2012, ApJ, 752, 45Tomozeiu M., Mayer L., Quinn T., 2016, ApJ, 827, L15Torrealba G., Koposov S. E., Belokurov V., Irwin M., 2016, MNRAS, 459,2370Torrealba G., et al., 2019, MNRAS, 488, 2743Ural U., Wilkinson M. I., Read J. I., Walker M. G., 2015, Nature Communi-cations, 6, 7599Walker M. G., Peñarrubia J., 2011, ApJ, 742, 20Walker M. G., Mateo M., Olszewski E. W., Gnedin O. Y., Wang X., Sen B.,Woodroofe M., 2007, ApJ, 667, L53Walker M. G., Mateo M., Olszewski E. W., Peñarrubia J., Wyn Evans N.,Gilmore G., 2009, ApJ, 704, 1274Weisz D. R., et al., 2019a, MNRAS, 489, 763Weisz D. R., et al., 2019b, ApJ, 885, L8Wojno J., Gilbert K. M., Kirby E. N., Escala I., Beaton R. M., TollerudE. J., Majewski S. R., Guhathakurta P., 2020, arXiv e-prints, p.arXiv:2004.03425Wolf J., Martinez G. D., Bullock J. S., Kaplinghat M., Geha M., MuñozR. R., Simon J. D., Avedo F. F., 2010, MNRAS, 406, 1220Zolotov A., et al., 2012, ApJ, 761, 71Zoutendijk S. L., Brinchmann J., Bouché N. F., den Brok M., KrajnovićD., Kuijken K., Maseda M. V., Schaye J., 2021, arXiv e-prints, p.arXiv:2101.00253de Blok W. J. G., McGaugh S. S., Bosma A., Rubin V. C., 2001, ApJ, 552,L23van den Bosch F. C., Ogiya G., 2018, MNRAS, 475, 4066van der Marel R. P., 1994, MNRAS, 270, 271von Hoerner S., 1957, ApJ, 125, 451MNRAS , 1–16 (2020) M. L. M. Collins et al. ρ D M ( M (cid:12) k p c − ) True100 tracers1,000 tracers10,000 tracers − r (kpc)0 . . . . . ˜ β ( r ) − r (kpc)0 . . . . . Figure A1.
Testing the updated GravSphere code and its new binulatormodule on spherical mock data from the Gaia Challenge suite. The leftpanels show results for the PlumCoreOm mock; the right panels for thePlumCuspOm mock. The top panels show the recovery of the density profilefor 100 tracers (green), 1,000 tracers (blue) and 10,000 tracers (black), wherethe width of the bands mark the 95% confidence intervals in each case. Thered dashed lines show the true solutions. The bottom panels show the samefor the recovery of the symmetrised velocity anisotropy profile. Notice that,even with 100 tracers, we obtain an unbiased recovery of the radial densityprofile. For the PlumCuspOm mock, we are also able to weakly detect theradial anisotropy at large radii. As the sampling increases, the confidenceintervals narrow around the correct solution.
APPENDIX A: TESTING
BINULATOR + GRAVSPHERE
ONMOCK DATA
In §4.1, we introduced an updated version of the GravSpherecode (Read & Steger 2017; Read et al. 2018) with a new reworkedbinning method, the binulator. In this Appendix, we present testsof the new code on mock data drawn from the Gaia Challenge spherical and triaxial suite (Read et al. 2021). We focus on twomocks that were particularly challenging for the previous versionof GravSphere: PlumCuspOm and PlumCoreOm. The former is aPlummer sphere embedded in a spherical, cuspy dark matter halo;the latter is the same embedded in a dark matter halo with a centralconstant density core. In both cases, the velocity dispersion profileis isotropic at the centre and maximally radially anisotropic at largeradii. Both mocks are described in detail in Read et al. (2021).Here, we model both assuming Gaussian velocity errors on each starof 2 km/s for 100, 1,000 and 10,000 randomly sampled kinematictracers. We assume that the photometric light profile is in all caseswell-sampled with 10,000 stars, similarly to the situation with realdata for nearby dwarf spheroidal galaxies (e.g. Read et al. 2019).Previously, GravSphere struggled on these mocks. For bothPlumCoreOm and PlumCuspOm with 1,000 tracers, it returned adark matter density that was systematically low for 𝑅 > 𝑅 / ,shifting, however, to 𝑅 > 𝑅 / for 10,000 tracers. The velocityanisotropy was well-recovered in all cases, except PlumCoreOm for http://astrowiki.ph.surrey.ac.uk/dokuwiki/ which 1,000 tracers was insufficient to detect the radial anisotropy atlarge radii (Read et al. 2021). In Gregory et al. (2019), GravSpherewas tested on mock data for fewer tracers (20, 100 and 500 tracerstars). There, a similar bias towards low densities at large radii wasreported.We show results applying the updated GravSphere with itsnew binulator binning module to the PlumCoreOm and Plum-CuspOm mocks in Figure A1. In all cases, we obtain an unbiasedrecovery of both the dark matter density profile, 𝜌 DM ( 𝑟 ) , and thesymmetrised stellar velocity anisotropy, ˜ 𝛽 ( 𝑟 ) . The code is no longerbiased to low density at large radii. This owes to a combination ofthe improved data binning, the proper inclusion of the full errorPDF for the virial shape parameters, and to our switch to usingthe coreNFW profile (recall that the non-parametric power-law-in-bins profile used previously by default in GravSphere biasesmodels towards falling very steeply at large radii in the absence ofdata). Notice that now, even with just 100 tracers, we are able todetect that the PlumCuspOm is denser and falls more steeply thanPlumCoreOm. For PlumCuspOm, we are also able to weakly detectthat the velocity distribution is radially anisotropic at large radii.As we move to 1,000 and 10,000 tracers, we obtain an increasinglyhigh-fidelity recovery of both the density profile and the velocityanisotropy. This paper has been typeset from a TEX/L A TEX file prepared by the author.MNRAS000