Mon.Not.R.Astron.Soc.000,1–11(2014) Printed18January2017 (MNLATEXstylefilev2.2) A test for skewed distributions of dark matter, and a detection in galaxy cluster Abell 3827 Peter Taylor1,2(cid:63), Richard Massey1,3, Mathilde Jauzac1,3,4, Fr´ed´eric Courbin5, David Harvey5, R´emy Joseph5 & Andrew Robertson3 1 Centre for Extragalactic Astronomy, Durham University, South Road, Durham DH1 3LE, UK 2 Mullard Space Science Laboratory, University College London, Holmbury St Mary, Dorking, Surrey RH5 6NT, UK 3 Institute for Computational Cosmology, Durham University, South Road, Durham DH1 3LE, UK 4 Astrophysics and Cosmology Research Unit, School of Mathematical Sciences, University of KwaZulu-Natal, Durban 4041, South Africa 7 5 Laboratoire d’astrophysique, Ecole Polytechnique F´ed´erale de Lausanne, Observatoire de Sauverny, CH-1290 Versoix, Switzerland 1 0 2 Accepted—.Received—;inoriginalform18January2017 n a J ABSTRACT Simulations of self-interacting dark matter (SIDM) predict that dark matter should 6 lag behind galaxies during a collision. If the interaction is mediated by a high-mass 1 forcecarrier,thedistributionofdarkmattercanalsodevelopasymmetricdarkmatter ] tails.Tosearchforthisasymmetry,wecomputethegravitationallensingpropertiesof O a mass distribution with a free skewness parameter. We apply this to the dark matter C aroundthefourcentralgalaxiesinclusterAbell3827.Inthegalaxywhosedarkmatter h. peakhaspreviouslybeenfoundtobeoffset,wemeasureskewnesss=0.23+−00..0252 inthe same direction as the peak offset. Our method may be useful in future gravitational p lensing analyses of colliding galaxy clusters and merging galaxies. - o Key words: dark matter — astroparticle physics — galaxies: clusters: individual: r t Abell 3827 — gravitational lensing: strong s a [ 1 v 1 INTRODUCTION 2016).Thesescatteringscantemporarilyseparatedarkmat- 2 ter from its associated galaxies. Such dark matter lags be- MostofthemassintheUniverseisdarkmatter(e.g.Planck 1 hind the galaxies, toward the position of diffuse gas that is 4 Collaboration2016).Darkmatterappearsinvisible,because slowed by ram pressure (Clowe et al. 2004; Lage & Farrar 4 itdoesnotinteract(orinteractsveryweakly)withStandard 2014; Harvey et al. 2015). Furthermore, if the particle in- 0 Model particles including photons. teractionismediatedbyahigh-massforcecarrier,thescat- . As the nature of dark matter remains unknown, there 1 tered dark matter will typically be ejected from the halo isnoreasontoaprioriassumeaparticulartheoryofitsori- 0 into asymmetric tails in the distribution of dark matter 7 gin.Thewiderangeofproposeddarkmattermodelspredict (Kahlhoefer et al. 2014). 1 different spatial distributions, particularly on small scales. Gravitationallensingoffersthemostdirectwaytomap : Darkmatterparticlesthatinteractwitheachother(SIDM) v the spatial distribution of dark matter, and hence to infer were proposed in (Spergel & Steinhardt 2000) to explain i itsparticleproperties.Gravitationallensingreferstothede- X small scale discrepancies between observations and simula- flection of light rays passing near any mass, including dark tions of collisionless dark matter. In the SIDM paradigm, r matter.Thankstothisdeflection,(unrelated)objectsbehind a energytransferbetweenparticlesmakesthecentreofgalax- darkmatterappearcharacteristicallydistorted,orevenvis- ies (Vogelsberger et al. 2012) and galaxy clusters (Rocha iblealongmorethanone(curved)line-of-sight.Eventhough etal.2013)morecircularandlessdense,potentiallyresolv- dark matter is invisible, it is possible to invert this process ing the core/cusp problem. Small substructures can also be andinferwhereitmustbe,byundistortingtheobservedim- erased–leadingtotheobservedunderabundanceofgalaxies ages, or ray-tracing multiple images back onto each other. in the Local Group, relative to simulations. During merg- Galaxy cluster Abell 3827 (22h01(cid:48)49(cid:48).(cid:48)1 −59◦57(cid:48)15(cid:48)(cid:48), ers between galaxies or galaxy clusters, dark matter in- redshift z=0.099) is particularly well suited for this kind of teractions transfer momentum between the colliding dark study. It gravitationally lenses a z=1.24 galaxy with spi- matter haloes (Randall et al. 2008; Kahlhoefer et al. 2014; ral arms and several knots of star formation that can be Robertson et al. 2017; Kim et al. 2016; Robertson et al. treated as independent background sources (Massey et al. 2015, hereafer M15). While most clusters contain only one (cid:63) e-mail:[email protected] brightest central galaxy, Abell 3827 contains four equally- (cid:13)c 2014RAS 2 P. Taylor et al. bright galaxies within its central 10kpc (Carrasco et al. N 2u0ra1t0i;onWmilleiaamnss t&haStashoam2e01o1f).thTehigsalhaixgihelsyaupnpuesauralclcoosnefitgo- -59.942 E gravitationallylensedimages.Thus,thedistributionoftheir 2’’ dtearr’ksrmelaatttievrecparnoxbimemityea(siunrteedrmprseocfisgerlya.vBitaectiaounsaelolefntshiengc)lu,sit- -59.943 Ae1 Af1 ispossibletoresolvesmallspatialoffsetsbetweenthedistri- Ad1 bution of dark matter and stars in the foreground galaxies. Ao1 Ae3 proaIcnhtthoissepaarpcehrfwore tphreesepnrtedaicnteedwapsayrmammeettrriyc liennstihnegdaips-- -59.944 Ab1Aa1 Ac1 AcA3a3 Ao3 Af3 Ab3 tribution of dark matter during mergers. A previous search N3 Ad3 lnoeonkte(dHfaorrvreeysiedtuaall.s2a0ft1e7r),subbuttrtahcatitnmgathyebseylmesmsesterniscitciovmepboe-- -59.945 N2 Ao7 N6 cause a tail of scattered particles shifts the best-fit position oofftthheesyremsimduetarli.cWcoemipnosnteeandtbcaocnkswtraurcdts,athsuinsgrleemhoavliongmsoomdeel -59.946 N1 Ac4Ao4 Ao8 N4 Af2 with a free skewness parameter that captures the asymme- Aa5 Aa4 Ae4 Ad2 Ae2 try found in simulations. We implement and distribute this Aa6 Ab4 Af4 Ao2 metoadle.l20in07t)h.eWpeubtelisctlyitaovnaiblaobtlhemLoecnkstdoaotal,wsohftewrearteh1e(sJkuelwlo- -59.947 Ad4 Ab2 Ac2 ness of the lens is known a priori, and on Abell 3827. Sec- Aa2 ttiioonn 32indtersocdriubceesseoxuirstninegwopbasrearmvaettiroincsleonfsAmboedlell.38S2e7ct.ioSnec4- -59.948 contains an analysis of Abell 3827. To be consistent with M15 we assume throughout this paper a flat ΛCDM cos- AmtoltohgeyrwedisthhifΩtMof=A0b.3e,llΩ3Λ8=270,.71(cid:48)(cid:48)ancodrrHes0p=o7n0dksmtos−1.18M28pkcp−c1.. -59.949 330.475 330.470 33 2 DATA Figure 1. Hubble Space Telescope image of galaxy cluster Abell 3827 in F160W (red), F814W (green) and F330W (blue) Broad-band imaging of Abell 3827 has been obtained from bands, after using MuSCADeT (Joseph et al. 2016) to fit and theGeminitelescopeatopticalwavelengths(Carrascoetal. subtractforegroundemission.ResidualemissionfromtwoMilky 2010) and from the Hubble Space Telescope (HST) in the Waystarshasbeenmaskedout,andremainsvisibleatlowlevel UV, optical and near-infrared (M15). This revealed four aroundthefourbrightcentralgalaxiesN1–N4.Circlesshowmul- similarly-bright elliptical galaxies within 10 kpc of each tipleimageidentifications,withtheradiusofthecirclereflecting other, plus a background spiral galaxy, whose multiply- uncertaintyontheirpositions(Ao8hasonlybeendetectedfrom theground). lensedimagesarethreadedthroughouttheclustercore(fig- ure 1). IntegralFieldUnitspectroscopyhasbeenobtainedfrom the VLT. An initial 1 hour exposure with the Multi-Unit perturbation around the Pseudo Isothermal Elliptical Mass SpectroscopicExplorer(MUSE)identifiedfourmaingroups Distribution (PIEMD; Kassiola & Kovner 1993). of lensed images, and suggested two low S/N peaks as can- didates for a demagnified central image (M15). A subse- 3.1 Pseudo-Isothermal Elliptical Mass quent additional 4 hour exposure (programme 295.A-5018; Distribution Massey et al. in prep.) confirms both candidates (Ao7 at RA: 330.47047, Dec: −59.945183, Ao8 at RA: 330.47079, The 2D surface mass density Σ of a circularly symmetric Dec: −59.946112). Indeed, Ao7 is also visible in HST imag- pseudo-isothermal mass distribution, projected along a line ing, after using the MuSCADeT multiwavelength method of sight, is: (Joseph et al. 2016) to estimate and subtract bright fore- (cid:32) (cid:33) σ2r 1 1 groundemission(Figure1).Wethereforeusealltheimages Σ(r)≡ 0 cut √ − , (1) identified by M15, plus the two new ones. 2G(rcut−rcore) rc2ore+r2 (cid:112)rc2ut+r2 whereσ isthe1Dvelocitydispersion,andwherer (r ) 0 core cut isaninner(outer)radius.ToconvertthisintoaPIEMDwith 3 LENS MODELLING METHOD ellipticity (cid:15) = a−b (cid:62) 0, where a and b are the semi-major a+b andsemi-minoraxesrespectively,Kassiola&Kovner(1993) We shall model the distribution of mass in the galaxy clus- apply their coordinate transformation (2.3.6): teras asum ofcluster-scale plusgalaxy-scale halos (follow- x y ing e.g. Limousin et al. 2007; Jauzac et al. 2014), each a x→x = , y→y = . (2) em 1+(cid:15) em 1−(cid:15) This maps a circle onto an ellipse centered at the origin, 1 http://projets.lam.fr/projects/lenstool/wiki with its major axis along the x axis. In general, including (cid:13)c 2014RAS,MNRAS000,1–11 Dark matter in galaxy cluster Abell 3827 3 a rotation to set the major axis at angle φ , this can be (cid:15) expressed in polar coordinates as: r2 →r2 = r2 (cid:2)1+(cid:15)2−2(cid:15)cos(cid:0)2(θ−φ )(cid:1)(cid:3). (3) em (1−(cid:15)2)2 (cid:15) Theangleαbywhichalightrayisdeflectedasitpasses nearthelens,andtheequivalent2Dgravitationalpotential ψ canbecomputedbyintegratingthedensitydistribution: α(r)=4GDlDls (cid:90)Σ(cid:0)r(cid:48)(cid:1) r−r(cid:48) d2r(cid:48) c2 Ds |r−r(cid:48)|2 (4) ψ(r)=4GDlDls (cid:90)Σ(cid:0)r(cid:48)(cid:1)log(cid:12)(cid:12)r−r(cid:48)(cid:12)(cid:12) d2r(cid:48) , c2 D s where D, D and D are the angular diameter distance l s ls fromtheobservertothelens,observertothesource,andlens to the source respectively. For general mass distributions, these integrals are difficult to solve – but closed forms have been found for the PIEMD, using techniques from complex analysis that exploit its elliptical symmetry (Bourassa & Kantowski 1975). A related halo model is the Pseudo Isothermal Ellipti- cal Potential (PIEP; Kassiola & Kovner 1993). In this, the coordinate transformation is applied to a circular potential Figure2.IsodensitycontoursofPISMDmassdistributions.The ψ(ratherthanthedensity).Itisthenmathematicallyeasier core radius is about the same as the innermost density contour. toobtainthedeflectionangleanddensityviadifferentiation: Thick, black lines show an ordinary PIEMD; the bottom-left is circular, and the others have ellipticity (cid:15)=0.15, with the major axisatvariousangles.Thinner,greylinesshowthesamedensity α(r)=∇ψ(r) profileswithskews=0.1,0.2,0.3totheright. c2 D (5) Σ(r)= s ∇2ψ(r). 8πGDD l ls In detail, the PIEP potential ψ is transformed so that: potential corresponding to the PIEMD.2 Analytic (albeit cumbersome) expressions for deflection angle and density ψ(x,y)→ψ(cid:48)(x,y)≡ψ(cid:0)x(cid:48),y(cid:48)(cid:1). (6) can be readily calculated via differentiation (equation 5). Thefirstandsecondderivativescanthenbecomputedwith We denote this the Pseudo Isothermal Skewed Mass Dis- applications of the chain rule. For example, the first x- tribution (PISMD); its isodensity contours are shown for derivative of the potential is: various values of (cid:15) and s in Figure 2. Foranyskew,thepeakofthemassdistributionalways ψx(cid:48) =(cid:0)ψx(cid:48)(cid:0)x(cid:48),y(cid:48)(cid:1)x(cid:48)x+ψy(cid:48)(cid:0)x(cid:48),y(cid:48)(cid:1)yx(cid:48)(cid:1)(cid:12)(cid:12)(x,y), (7) lies at the same position, so it will be possible to use the samemetricasM15tomeasureanyoffsetbetweenthemost where the subscript denotes partial differentiation. The re- gravitationallyboundstarsanddarkmatter.Thetotalmass sulting mass distribution is not the same as a PIEMD, be- changesslightlywithincreasingskew,butthiscanberecal- cause of the way the coordinate transformation propagates culated after a fit. through the chain rule (or back up the integrals in equa- Like the PIEP, the density distribution of the PISMD tion 4). For large (cid:15), the mass distribution corresponding to exhibits undesired behaviour with large skews, because aPIEPhasundesirablefeaturesincludingconcave(peanut- the coordinate transformation was applied to ψ, not Σ. shaped) isodensity contours (Kassiola & Kovner 1993). Isophotes of the density distribution become concave, and the skew ellipticity can overwhelm the underlying elliptic- 3.2 Pseudo-Isothermal Skewed Mass Distribution ity. We avoid these effects by restricting |s| < 0.3.3 Since the PISMD is invariant under transformations s → −s To perturb the mass distribution in a way that resembles and φ → φ + π, we fit parameters within the domain the behaviour of SIDM in numerical simulations (see figure s s s ∈ [−0.3,0.3] and φ in some interval of length π. This 5 of Kahlhoefer et al. 2014), we apply a further coordinate s ensures that the parameter space can be explored symmet- transformation that maps a circle onto an ellipse with its focus (rather than centre) at the origin: r2(cid:0)1−s2(cid:1)3/2 2 WewouldideallyapplythistransformationtothePIEMDmass r2 →r(cid:48)2 = (8) distribution, but the relevant integrals (equation 4) do not con- (1+scos[θ−φs])2 tracttoasimpleform.Askewedmassdistributioncouldalsobe derived from the potential corresponding to a PIEP. We choose with s being the third eccentricity such that s = (cid:112) toperturbthePIEMDsothatwerecoverthiswidely-usedmass 1−b2/a2, and the power 3/2 being introduced to pre- distributioninthes→0limit,andtominimiseundesiredconvex serve area. Note the asymmetric cos(θ) terms rather than curvatureindensityisophotes. thecos(2θ)termsinthemappingdescribedbyequation(3). 3 SeetheAppendixforanalternativemodelthatdoesnotsuffer We apply this transformation to the 2D gravitational fromthiseffect,buthasotherdisadvantages. (cid:13)c 2014RAS,MNRAS000,1–11 4 P. Taylor et al. rically about s = 0, allowing the case of zero skew to be implemented the polar option, so that Lenstool’s MCMC recovered without bias. We set the edge of the domain of optimiser can explore a circularly symmetric region, with φ wellawayfromanypreferreddirection(inpractice,hav- no preferred direction that could bias the inferred skew. s ing explored parameter space via a quick search), to make Lenstool also defines ellipticities in this way, for the same sure an MCMC sampler operates efficiently near regions of reason. interest. Nonetheless,itmayoftenbedesirabletoknowthepos- teriorprobabilitydistributionofskewnessalonge.g.adirec- tionofmotion,andperpendiculartothat(i.e.theCartesian 3.3 Testing an implementation in Lenstool form). The posterior probability distributions of skew and We have implemented the PISMD as potential 813 in the skew angle are returned by Lenstool (in runmode=3) by publicly-available software Lenstool (Jullo et al. 2007). thesamplingdensityoftheMCMCchain.Thiscanbecon- Given a parameterised mass distribution, and the location verted to the posterior of the skew in some direction φ by ofbackgroundsources,Lenstoolcancomputetheposition projecting and then weighting each sample by: ofobservedmultipleimages.Giventhepositionofobserved |s| images,itcanalsouseMarkov-ChainMonteCarlo(MCMC) w= (cid:112) . (9) 0.32−s2cos2(φ −φ) optimisation to fit parameters of the lensing mass distribu- s tion. The numerator is the Jacobian to convert the area of pa- To test whether Lenstool can accurately recover a rameterspacefrompolartoCartesiancoordinates.Thede- known input skew, we run two sets of tests. We first con- nominator corrects for prior bias, because the restriction sider an isolated lens, with three background sources at |s| ∈ [−0.3,0.3] leads to a (semi-)circular prior on the pro- different redshifts: the example with images configuration jected skew. that is packaged with Lenstool. As a null test, we adopt the input mass distribution with skew s = 0. From true the observed positions of multiple images, Lenstool suc- 4 STRONG LENS ANALYSIS OF ABELL 3827 cessfully recovered a best fit (maximum likelihood) value s=−0.0008+0.02. We use the observed positions of lensed multiple images to −0.02 We then set skew s = 0.2 and φ = 1.6. We set fitamassmodelofthecluster.Ourchoiceofmodelparam- true strue sourcepositionsbyprojectingoneimageofeachsourceback eters and their priors is based on those of M15, with some toitssourceplane,thencreateamocksetofmultiplylensed additional degrees of freedom. We assume 0.8(cid:48)(cid:48) uncertainty imagesbyre-projectingthissourceforwardthroughthelens. on the position of lensed image Ao8, which has only been When fitting this mock data, Lenstool successfully recov- detected from the ground. We assume 0.2(cid:48)(cid:48) uncertainty on ers best-fit values s=0.2+0.001 and φ =1.6+0.04. the position of all other lensed images, which are identified −0.001 s −0.05 Second, we test the recovery of input skews in a com- by HST. plex cluster lens with a mass distribution based on the Abell 3827. Choosing one of the quadruply lensed back- 4.1 Fiducial mass model ground galaxy images, we repeat the procedure outlined above: projecting the light backwards and then forwards The cluster’s large-scale mass distribution is modelled as a through a cluster lens with known mass distribution. The single PIEMD. Based on a comprehensive (but slow) initial clusterisgiventhesameparametersasourfiducialmodelfor exploration of parameter space, its position is given by a Abell3827(see§4.1),withtheexceptionoftheskewparam- broad Gaussian prior with σ = 2(cid:48)(cid:48) = 3.66kpc, centred on eters. In this test, the dark matter associated with galaxy the position of galaxy N2. Flat priors are imposed on its N1isgivenskewnessstrue =0.25andφstrue =1.6.Asanull ellipticity ((cid:15) < 0.75), core size (rcore < 40(cid:48)(cid:48)) and velocity test, galaxies N2–N4 are set to have no skew, strue =0. dispersion (300<σv<1000km/s). Its cut radius is fixed at We run Lenstool with the same free parameters and r =1000(cid:48)(cid:48),welloutsidethestronglensingregion,i.e.away cut priors as in §4.1). Within such a highly dimensional pa- from any multiple image constraints. rameter space, we find that the best-fit values are some- Central galaxies N1–N4 are each modelled as a stellar timesnoisy,forparametersthatmakeonlyasmalldifference component (which was not included in the fiducial model to the overall goodness of fit. However, the full posterior of M15), plus a dark matter one. Following Giocoli et al. probabilitydistributionfunction(PDF)issmoothandwell- (2012),thestellarcomponentsaremodelledwithHernquist sampled. Hence, for the rest of this paper, we shall quote (1990) profiles: themodalpeakand68%widthoftheposteriorPDF,which ρ Lenstoolalsoreturns.Thismakesnodifferenceforthesim- ρ (r)= s , (10) plemodelabove,andsuccessfullyrecoverss=0.24+0.04 and star (r/rs)(1+r/rs)3 −0.31 φs = 1.60−.09.299 for galaxy N1, and s = 0.01+−00..1143, 0.07+−00..1105, where the scale radius rs is related to the half mass radius 0.11+0.11 for galaxies N2, N3, N4 (with very weakly con- R , such that R = r /0.551, and the scale density ρ = −0.16 e e s s strained φs). Mtotal/(cid:0)2πrs3(cid:1). We fix the mass of the stellar component, and its half-mass radius, using the optical magnitudes and profiles measured by M15. These parameters are listed in 3.4 Prior bias for polar parameters Table 1. Askewisatwo-componentvector,andcanbeexpressedin The four central galaxies’ dark matter components are polarformasamagnitude|s|anddirectionφ ,orinCarte- now modelled as PISMDs. We impose flat priors on their s sianformasanamountinorthogonaldirections{s ,s }.We positions,in4(cid:48)(cid:48)×4(cid:48)(cid:48)boxescentredontheirluminositypeaks, x y (cid:13)c 2014RAS,MNRAS000,1–11 Dark matter in galaxy cluster Abell 3827 5 Table 1. ParametersofthefiducialmassmodelfittedbyLenstool.Quantitiesinsquarebracketsarefixed.Errorsonotherquantities show68%statisticalconfidencelimits,marginalisingoveruncertaintyinallotherparameters.Stellarmasscomponentsaremodelledas Hernquistprofiles,withamass(computedfromfluxintheF606Wband),scaleradiusandellipticity(fittedusingGalfit;galaxyN4is contaminated by a nearby star). Dark matter components are modelled as PISMDs, with a 1D velocity dispersion, core and cut radii, ellipticityandskewness.Positionsaregiveninarcsecondsrelativeto(R.A.:4330.47515,Dec.:−59.945996),exceptgalaxies’darkmatter components,whicharerelativetothepositionoftheirstars.AnglesareanticlockwisefromEast. ∆xx[(cid:48)[(cid:48)(cid:48)](cid:48)] ∆yy[(cid:48)[(cid:48)(cid:48)](cid:48)] Mσvas[skm[M/(cid:12)s]] rrcoscre[(cid:48)[(cid:48)(cid:48)](cid:48)] rcut[(cid:48)(cid:48)] (cid:15) φ(cid:15)[◦] s φs[◦] N1 stars [−0.06] [0.04] [1.00×1011] [0.53] [0.12] [61] darkmatter −0.29+0.25 −0.71+0.30 149+8 [0.1] [40] 0.02+0.33 151+19 0.21+0.06 86+44 −0.14 −0.16 −12 −0.01 −116 −0.22 −44 N2 stars [5.07] [2.05] [2.46×1011] [0.79] [0.17] [39] darkmatter −0.23+0.30 0.00+0.30 182+29 [0.1] [40] 0.42+0.05 23+32 0.03+0.11 117+41 −0.16 −0.30 −22 −0.22 −12 −0.14 −80 N3 stars [9.69] [3.98] [2.77×1011] [0.33] [0.05] [31] darkmatter −0.05+0.25 −0.06+0.18 213+8 [0.1] [40] 0.49+0.01 15+14 −0.02+0.08 169+7 −0.25 −0.29 −10 −0.16 −8 −0.11 −109 N4 stars [9.26] [−1.08] [2.08×1011] [1.37] [0.39] [127] darkmatter −1.35+0.39 0.51+0.35 255+8 [0.1] [40] 0.02+0.25 136+17 0.08+0.08 147+21 −0.34 −0.27 −10 −0.01 −28 −0.09 −80 N6 stars [18.54] [2.47] [0] darkmatter [0] [0] 38+26 [0.1] [40] [0] [0] [0] [0] −25 Cluster dm 5.53+1.46 2.33+1.97 683+139 30.12+9.23 [1000] 0.56+0.13 63+2 [0] [0] −1.61 −1.59 −75 −6.43 −0.10 −3 plus flat priors on their ellipticity ((cid:15) < 0.5) and velocity dispersion (v < 600km/s). We fix r = 40(cid:48)(cid:48) = 73kpc disp cut (Limousin et al. 2007). Galaxy N6 is much fainter than the others, so we ap- proximate its total mass distribution as a single PIEMD. This has a fixed position and ellipticity to match the light N3 distribution, and only its velocity dispersion is optimised (with a flat prior v <500km/s). disp N2 We optimise the free parameters using Lenstool, with runmode=3 and BayesRate=0.1 (Jullo et al. 2007). (Modal) maximum likelihood parameters are shown in Ta- ble 1, and the corresponding mass distribution is shown in Figure3.ThebestfitmodelachievesaRMS offsetbetween N1 the observed and predicted positions of multiple images of N4 (cid:104)rms(cid:105) =0.26(cid:48)(cid:48), which is consistent with 0.2(cid:48)(cid:48) observational i uncertainty. The modal χ2/dof=67.1/19 with a log likeli- hood of 8.18. The full posterior probability distribution for thedarkmatterassociatedwithgalaxiesN1–N4isshownin Figures 4 and 5. Colours: mass in stars Contours: total mass (white), dark matter belonging to galaxies (black) Figure 3.Thebestfittingmassdistributioninthegravitational 4.2 Sensitivity to model choices lensAbell3827,integratedalongourlineofsight.Forreference, 4.2.1 Flat priors on position thebackgroundcolourscaleshowsthemodelledstellarmassden- sity. Red spots indicate the position of the luminosity peak in M15usedflatpriorsonthepositionofgalaxyN1,butGaus- galaxiesN1–N4.Whiteisodensitycontoursshowthetotallensing sian priors on N2–N4. This choice was intended to fix the massofthecluster.Theoutermostcontourcorrespondstoapro- position of galaxies that are further from (hence less con- jected density of 2×109M(cid:12)/kpc2, and values increase towards thecentrebyafactorof21/3=1.26.Blackisodensitycontoursiso- strainedby)lensedimages,sothatthepositionofN1could late each galaxy’s dark matter component. The outermost con- be explored in detail. There have since been two challenges to this line of reasoning. First, the position of the source tour corresponds to a projected density of 1.26×109M(cid:12)/kpc2 andvaluesincreasebyafactorof22/3.Thevisibleoffsetbetween should not influence a choice of prior. Second, giving only starsanddarkmatteringalaxiesN1andN4arebothstatistically galaxy N1 a flat prior can transfer freedom in the model significant;theasymmetryinthedistributionofN1’sdarkmatter from potential offsets in other galaxies to the position of isalsosignificant. N1, causing an overestimate in its offset (and possible off- sets in the other galaxies to be missed). Here, we impose flatpriorsonthepositionofdarkmatterassociatedwithall galaxies. (cid:13)c 2014RAS,MNRAS000,1–11 6 P. Taylor et al. N1 N2 2 2 y []00 01 y []00 01 1 1 0.5 0.5 0.4 0.4 0.3 0.3 0.2 0.2 0.1 0.1 160 160 []◦12800 []◦12800 40 40 km/s]223050000 km/s]223050000 [v150 [v150 0.15 0.15 s00..0105 s00..0105 d] 2.4 d] 2.4 [ras01..86 [ras01..86 1 0 1 2 1 0 1 2 0.10.20.30.40.5 40 80 120160 150 200 250 300 0.150.000.15 0.8 1.6 2.4 1 0 1 2 1 0 1 2 0.10.20.30.40.5 40 80 120160 150 200 250 300 0.150.000.15 0.8 1.6 2.4 x [00] y [00] [◦] v [km/s] s s [rad] x [00] y [00] [◦] v [km/s] s s [rad] N3 N4 2 2 y []00 01 y []00 01 1 1 0.5 0.5 0.4 0.4 0.3 0.3 0.2 0.2 0.1 0.1 160 160 []◦12800 []◦12800 40 40 [km/s]v122350500000 [km/s]v122350500000 0.15 0.15 s00..0105 s00..0105 [rad]s012...864 [rad]s012...864 1 0 1 2 1 0 1 2 0.10.20.30.40.5 40 80 120160 150 200 250 300 0.150.000.15 0.8 1.6 2.4 1 0 1 2 1 0 1 2 0.10.20.30.40.5 40 80 120160 150 200 250 300 0.150.000.15 0.8 1.6 2.4 x [00] y [00] [◦] v [km/s] s s [rad] x [00] y [00] [◦] v [km/s] s s [rad] Figure 4. Posterior probability distribution of the distribution of dark matter associated with galaxies N1 (top left), N2 (top right), N3(bottomleft)andN4(bottomright).Contoursshowthe68%,95%and99.7%contourlevels.Bluelinesindicatethebestfitmodel, butnotethatthiscanbenoisy,andweinsteadusethepeakofthesmoothedposteriordistribution.Positionshavebeenrecenteredsuch that(x,y)=(0,0)isthepeakofthestellarluminosity.Offsetsbetweenstarsanddarkmatteraremeasuredat>3σforgalaxiesN1and N4.Askewisdetectedat>1σ forgalaxyN1,inadirectionconsistentwiththespatialoffset. 4.2.2 Stellar mass components 4.2.3 Identification of new lensed images Adding constraints from the two new lensed images Ao7 andAo8tightensconstraintsonnearbygalaxiesN3andN4. Galaxiesdefinitelycontainstars,andthosestarshavemass. These (demagnified) images are unresolved, and any of the Explicitlyincludingastellarcomponentintogalaxies’mass featuresinthebackgroundspiralcouldbeassignedtothem. models is intended to limit bias in inferred skewness, es- We have tried relabelling one or both of the demagnified pecially for galaxies whose stellar mass appears to be off- images as either the bulge, Aon, or one of the two bright- set from their dark matter. In practice, as in M15, we find est knots of star formation, Aan or Abn. Lenstool’s out- that including the stellar mass component (or even multi- puts are statistically indistinguishable. In all cases, the en- plying/dividingitsmassbyafactor2)doesnotchangeany tire background spiral galaxy is predicted to be lensed onto other results, within their statistical errors. both the northern and the southern demagnified images. (cid:13)c 2014RAS,MNRAS000,1–11 Dark matter in galaxy cluster Abell 3827 7 2 ] y [00 01 1 1 N N1 s000...011055 ad] 2.4 [rs 1.6 1 0.8 N 2 ] x [00 01 4 1 N 2 ] y [00 01 4 1 N N4 s000...011055 ad] 2.4 [rs 1.6 4 0.8 N 10 x []00 05 C 5 10 y []00 05 C 5 1 0 1 2 1 0 1 2 0.150.00 0.15 0.8 1.6 2.4 1 0 1 2 1 0 1 2 0.150.00 0.15 0.8 1.6 2.4 5 0 5 10 5 0 5 10 N1 x [00] N1 y [00] N1 s N1 s [rad] N4 x [00] N4 y [00] N4 s N4 s [rad] C x [00] C y [00] Figure 5. Posterior probability distribution showing (minimal) correlations between the position and asymmetry of the dark matter associated with key galaxies, and between those galaxies and the cluster-scale dark matter(denoted C). Contours show the 68%, 95% and99.7%contourlevels.PositionshavebeenrecenteredasinFigure4. 5 DISCUSSION (Kahlhoeferetal.2015).Weagreewiththisinterpretation4, butnotethatthecross-sectioncanbelowerifthegalaxyhas 5.1 Galaxy N1 completedmultipleorbits;itscurrentdirectionofmotionis The previously-detected offset between galaxy N1’s stars unknown. and dark matter persists at > 3σ in our new analysis. Wealsofind(atmuchlowerstatisticalsignificance)that Adding two free parameters for the asymmetry of its dark galaxy N1’s dark matter is skewed in a direction consistent matter slightly increases uncertainty in its position. The withtheSIDMinterpretationofitsoffset.Figure6showsthe modal offset is (∆x,∆y) = (cid:0)−0.22+0.25,−0.81+0.16(cid:1) for an posteriorPDFoftheskewontothevectorpointingfromthe unskewed model, and (∆x,∆y) = (cid:0)−−0.014.29+0.25,−−0.017.71+0.30(cid:1) DMpeaktothestellarluminosityinthefiducialmodel,such −0.14 −0.16 if skewness is allowed. that a positive skew corresponds to the direction predicted The consistency between these suggests that the mea- sured position of the density peak is robust to the ex- act shape of the DM halo. If the offset is entirely due 4 WehavealsorepeatedKahlhoeferetal.(2015)’scalculationof to dark matter self-interactions, it implies a momentum- σ˜/mDM but integrating the effect of the restoring force on the transfer interaction cross-section σ˜/mDM (cid:38) 1cm2g−1, as- entire distribution rather than just the peak. The difference is suminggalaxyN1isfallingintotheclusterforthefirsttime notsignificant. (cid:13)c 2014RAS,MNRAS000,1–11 8 P. Taylor et al. for scattered (self-interacting) dark matter in colliding sys- y sit tems. More generally, it will also be useful to investigate n e claims of dynamically-induced asymmetry (Prasad & Jog d y 2016; Chemin et al. 2016), or tidal tails (which are asym- bilit metricifthesizeofabodyislargecomparedtoitsdistance a In direction of offset b from the centre of potential). o erior pr Odirrtehcotigoonn oafl otoffset tfliaotnpWorifeormhsaafvsosertainhlseogpapolarsexitsyieonnctleuodsftaaelrlnAgewablaemlxlioe3ds8’e2ld7af.orrkOtmuhraetcdtheiosrticlreeiabduos-f ost to a detected offset between a second galaxy’s stars and P its dark matter. New VLT/MUSE observations tighten the -0.3 -0.2 -0.1 0 0.1 0.2 0.3 constraints on that offset. Neither measured offset changes Projected skewness of galaxy N1 significantly if the models are allowed extra freedom to be- come skewed. Figure 6. The posterior of the skew vectors in galaxy N1 pro- We find tantalising, but low significance evidence that jectedontotheoffsetvector(black)andorthogonaltothis(grey). the galaxy closest to multiply lensed images (and therefore The Jacobian has been accounted for and the priors have been adjusted as described in section 3.4. Dashed lines indicate the thebestconstrained)hasanasymmetricdistributionofdark posterior peak and 1σ confidence intervals. The blue line indi- matter,skewedinthesamedirectionasitsoffsetfromstars. cates the (noisy) best fit value. There is a preference for a skew More work will be needed to determine the significance of that is consistent with SIDM. No such preference is shown for a thisresult:whetheritisphysical,oranartefactofsystemat- skewcomponentintheorthogonaldirection. ics in parametric lens modelling. Even in mock data where the true skew is known, skewness cannot be measured to high precision in a system as complex as Abell 3827. This by SIDM. The peak of the posterior and 1σ errors are s = is probably because the effect of skewness on lensed image 0.23+0.05. If we individually project the skewness onto the −0.22 positionsissmallerthantheeffectsofotherfreeparameters. offset direction individually in all MCMC samples, we find A promising direction for future investigation may be s=0.26+0.03.Inbothcases∼70%oftheposteriorPDFlies −0.22 provided by pairs of field galaxy in the SLACS survey, one ats>0.Astheposteriorpeakisneartheedgeoftheprior, of which has already been found to have an offset between which is chosen coservatively (see §3.2), a parametric halo dark and luminous matter (Shu et al. 2016). Whilst the model that does not break down for large skew parameters SIDMmodelpredictsalargest(mosteasilyobservable)off- could result in a stronger detection. setingalaxiesmovingthroughadenseintraclustermedium, itmaybepossibletomoretightlyconstrainanyasymmetry ofdarkmatterinthesesimplersystems.Ifthedirectionsof 5.2 Galaxies N2 and N3 theirdarkmattertailscorrelatewiththedirectionsoftheir ThedarkmatterassociatedwithgalaxiesN2andN3appears offsets,thisevidencewouldsupportthehypothesisofSIDM symmetric, and coincident with the stars. This result does with a high-mass mediator particle. not preclude offsets from existing along the line of sight. Furthermore,evenifthegalaxiesaremovingintheplaneof the sky, they could be behind or in front of the most dense regions of the cluster core, and therefore passing through a ACKNOWLEDGEMENTS less dense medium, experiencing less drag. Theauthorsaregratefulforconstructiveconversationswith Benjamin Cl´ement, Eric Jullo, Felix Kahlhoefer, Thomas Kitching, Marceau Limousin and Subir Sarkar. Figure 4 5.3 Galaxy N4 was produced using the python module corner (Foreman- ThemeasuredpositionofN4’sdarkmatterisintriguing.Ac- Mackey 2016). PT is supported by the STFC. RM is sup- countingpurelyforstatisticalerrorbars,thankstothecon- ported by the Royal Society. AR and MJ are supported by firmedpositionsofdemagnifiedimages,wefindthatgalaxy the UK Science and Technology Facilities Council (grant N4isoffsetfromthegalaxy’sstarsatthe3σlevel.However, numbers ST/H005234/1 and ST/L00075X/1). the offset position is mildly degenerate with the position of Facilities:ThispaperusesdatafromobservationsGO-12817 thecluster-scaledarkmatter(Figure5),thusaflatprioron (PI: R.Massey) with the NASA/ESA Hubble Space Tele- theclusterscalehalocouldleadtoadifferentoffsetmeasure- scope, obtained at the Space Telescope Science Institute, ment. Furthermore, the measured ellipticity of the galaxy which is operated by AURA Inc, under NASA contract light is contaminated by light from an adjacent Milky Way NAS 5-26555. This paper also uses data from observations star,anditspositionmayalsobe.Thereareinsufficientcon- made with ESO Telescopes at the La Silla Paranal Ob- straintstodeterminewithstatisticalsignificancewhetheror servatory under programmes 294.A-5014 and 295.A-5018 not the distribution of dark matter is asymmetric. (PI: R.Massey). We thank the Paranal Science Operations teamforcarryingoutthoseobservations.OurLenstoolop- timisationusedtheDiRACDataCentricsystematDurham University,operatedbytheInstituteforComputationalCos- 6 CONCLUSION mology on behalf of the STFC DiRAC HPC Facility (www. Wehavedevelopedaparametriclensmodelsforasymmetri- dirac.ac.uk).ThisequipmentwasfundedbyBISNational cally skewed mass distributions. This can be used to search E-infrastructure capital grant ST/K00042X/1, STFC cap- (cid:13)c 2014RAS,MNRAS000,1–11 Dark matter in galaxy cluster Abell 3827 9 ital grant ST/H008519/1, STFC DiRAC Operations grant APPENDIX A: ALTERNATIVE METHOD TO ST/K003267/1 and Durham University. DiRAC is part of GENERALISE LENS MASS DISTRIBUTIONS the UK National e-Infrastructure. Anotherwaytointroduceasymmetryistoapplyaweighting function w(r;{a }) to an elliptical lensing potential i ψ(r)→ψ(cid:48)(r)≡w(r;{a })ψ(r), (A1) REFERENCES i where {a } are a set of parameters. The deflection and sur- BourassaR.,KantowskiR.,1975,TheAstrophysicalJour- i face mass density are readily computed by differentiating. nal, 195, 13 We consider weighting functions of the form Carrasco E. R., et al., 2010, ApJL, 715, L160 CheminL.,Hur´eJ.-M.,SoubiranC.,ZibettiS.,CharlotS., w(r;{a })=1+sf(r,θ) (A2) i Kawata D., 2016, A&A, 588, A48 where s is a skew parameter and f(r,θ) is written in polar CloweD.,GonzalezA.,MarkevitchM.,2004,ApJ,604,596 coordinates. The second term acts as a perturbation away DespaliG.,GiocoliC.,BonamigoM.,LimousinM.,Tormen from elliptical symmetry of O(s), with s=0 corresponding G., 2016, arXiv preprint arXiv:1605.04319 to the elliptically symmetric case. We chose f(r,θ) to meet Foreman-Mackey D., 2016, The Journal of Open Source the following criteria: Software, 24 GiocoliC.,MeneghettiM.,BartelmannM.,MoscardiniL., • To ensure that the space about s=0 is explored sym- Boldrin M., 2012, Monthly Notices of the Royal Astro- metricallyinLenstool,sothatanon-zeroskewisnotartifi- nomical Society, 421, 3343 cially recovered, we require that sf(r,θ)=−sf(r,θ+π). Harvey D., Massey R., Kitching T., Taylor A., Tittley E., • To avoid difficulties near the origin, we require f(r,θ) 2015, Science, 347, 1462 tobeanincreasingfunctionofr.Thisisalsophysicallymo- Harvey D., Robertson A., Massey R., Kneib J.-P., 2017, tivated, as it is difficult to scatter particles from the centre Monthly Notices, 464, 3991 of the potential well at r=0. Hernquist L., 1990, The Astrophysical Journal, 356, 359 • To ensure that the surface mass density remains posi- Jauzac M., et al., 2014, Monthly Notices, 443, 1549 tive (or becomes negative only for large r well outside any Joseph R., Courbin F., Starck J.-L., 2016, A&A, 589, A2 region of interest), we require f(r,θ) to be bounded. JulloE.,KneibJ.-P.,LimousinM.,EliasdottirA.,Marshall P., Verdugo T., 2007, New Journal of Physics, 9, 447 KahlhoeferF.,Schmidt-HobergK.,FrandsenM.T.,Sarkar A1 The Weighting Function Pseudo-Isothermal S., 2014, Monthly Notices, 437, 2865 Skewed Mass Distribution KahlhoeferF.,Schmidt-HobergK.,KummerJ.,SarkarS., 2015, Monthly Notices: Letters, 452, L54 MeetingtheaboveconditionsweformtheweightedPseudo- Kassiola A., Kovner I., 1993, The Astrophysical Journal, Isothermal Skewed Mass Distribution (wPISMD) by apply- 417, 450 ing the weighting function Kim S. Y., Peter A. H. G., Wittman D., 2016, preprint, (cid:18) (cid:19) r (arXiv:1608.08630) w(r;s,rs,β,φs)=1+s tan−1 βr cos(θ−φs), (A3) Lage C., Farrar G., 2014, ApJ, 787, 144 s Limousin M., Kneib J., Bardeau S., Natarajan P., Czoske to the PIEMD potential, where rs is a new scale radius, β O., Smail I., Ebeling H., Smith G., 2007, Astronomy & sets the radial dependence of the skew, and φs is the skew Astrophysics, 461, 881 angle.5 Thisisnowavailableaspotential812inLenstool. Massey R., et al., 2015, Monthly Notices, 449, 3393 The resulting surface mass densities are shown in Fig- Planck Collaboration 2016, A&A, 594, A13 ure A1. The qualitative shape of the isodensity contours PrasadC.,JogC.J.,2016,preprint, (arXiv:1610.01702) changesinsideoroutsidethescaleradius(owingtothesign Randall S. W., Markevitch M., Clowe D., Gonzalez A. H., changeinsecondderivativeoftheinversetangent).Thisfea- Bradaˇc M., 2008, ApJ, 679, 1173 turecouldbeusedtoisolateabehaviourthatbestmatches Robertson A., Massey R., Eke V., 2016, Monthly Notices numericalsimulations,byfixingverylargeorverysmallrs, of the Royal Astronomical Society, 465, 569 or to capture more complex behaviour. Robertson A., Massey R., Eke V., 2017, Monthly Notices, The total mass of a wPISMD is identical to that of of 465, 569 a PIEMD. Since the weighting function is normalised by Rocha M., Peter A. H. G., Bullock J. S., Kaplinghat M., construction, the integrated mass density of a PIEMD and Garrison-Kimmel S., On˜orbe J., Moustakas L. A., 2013, wPISMD over a circular region are the same: Monthly Notices, 430, 81 c2 D (cid:90) c2 D (cid:90) Shu Y., Bolton A. S., Moustakas L. A., Stern D., Dey A., s ∇2ψ dA= s ∇2(wψ) dA, 8πGDD 8πGDD BrownsteinJ.R.,BurlesS.,SpinradH.,2016,TheAstro- l ls |r|<R l ls |r|<R physical Journal, 820, 43 SpergelD.N.,SteinhardtP.J.,2000,PhysicalReviewLet- 5 Theinversetangentformoftheradialdependenceisnotphysi- ters, 84, 3760 callymotivated,andotherfunctionalformsmayalsowork.While Vogelsberger M., Zavala J., Loeb A., 2012, Monthly No- itismathematicallyunnecessarytohavetwodegenerateparam- tices, 423, 3740 etersrs andβ,toavoidcomputationaldivisionsbyzero,thedis- WilliamsL.L.,SahaP.,2011,MonthlyNoticesoftheRoyal tributedLenstoolimplementationuseshardcodedrs=0.1(cid:48)(cid:48)and Astronomical Society, 415, 448 allowsβ tovary. (cid:13)c 2014RAS,MNRAS000,1–11 10 P. Taylor et al. Figure A2. Isodensity contours for a radially varying weighted sumoftwoPIEMDswithdifferentellipticities. valuesofs,sincetheweightingfunctionhasbeenappliedto the potential and not the density. We have found that the valueofswherethisoccursissensitivetoβ andthecutand coreradii.Forthisreason,werecommendtestingthebound- aries of the parameter space for a breakdown of the desired skewed behaviour before substantial future work. Nonethe- less, we tested the wPISMD against the null hypothesis of the unskewed example with images system distributed with Lenstool (see §3.3). Fixing β = 0.01 and starting with a flat prior s ∈ [0.3,0.3], Lenstool successfully recovers skewness s=0.002+0.002. −0.002 A2 Pseudo Isothermal Varying Ellipticity Mass Distribution Despali et al. (2016) predict that, even with standard cold dark matter, the ellipticity of a cluster scale halo should changeasafunctionofradius,becomingmoreelongatedfur- therfromthecentre.Thispredictioncanbetestedbyusing theweightingfunctionformalismtoaddanextraparameter tohalomodelsthatmimicsthisbehaviour.Toachievethis, we combine a weighted sum of two different mass densities withdifferentellipticitiesintoaPseudoIsothermalVarying Ellipticity Mass Distribution (PIVEMD) We suggest a mass density of the following form: Σ(r)=Σ (r)w (r)+Σ (r)w (r) (A5) (cid:15)1 1 (cid:15)2 2 Figure A1.IsodensitycontoursofwPISMDmassdistributions, inside (top) and outside (bottom) the scale radius rs/β. Thick, where Σ(cid:15)1(r) and Σ(cid:15)2(r) are two elliptical profiles with el- blacklinesshowanordinaryPIEMD;thebottom-leftiscircular, lipticity (cid:15) and (cid:15) . All the other parameters for these two 1 2 and the others have ellipticity (cid:15)=0.15, at angles θ(cid:15)=0◦, 45◦ and densitiesshouldbeshared.Inthiscase,wefinditmosteffec- 90◦inthesameorderasinFigure2.Thinner,greylinesshowthe tive (and possible) to apply the weighting function directly samedensityprofileswithskews=0.1,0.2,0.3,0.4totheright. to the mass density. To be computationally efficient within an MCMC loop, deflection angles must be computed once, using numerical integration, and stored in a look-up table. (A4) The weighting functions w should meet the following i where θ-dependence cancels. Taking the limit as | r |→ ∞, criteria: thelefthandsidewillconvergetothetotalmassofaPIEMD • To normalise the total mass, w (r)+w (r)≡1, ∀r. 1 2 withellipticity(cid:15),andtherighttothemassofanequivalent • So that one ellipticity dominates at small r and the wPISMD. However, the position of the density peak varies other at large r, let w (r) → 1 as r → ∞, w (r) → 1 1 slightly as a function of s. Care would need to be taken if 0 as r→0, w (r)→0 as r→∞ and w (r)→1 as r→0. 2 2 using a wPISMD to measure offsets of dark matter. As with the PISMD, this model breaks down for large Although this is quite a general set of conditions, we (cid:13)c 2014RAS,MNRAS000,1–11