Cold Dark Matter Substructures in Early-Type Galaxy Halos

Accepted for publication in ApJ Preprint typeset using LATEX style emulateapj v. 01/23/15

COLD DARK MATTER SUBSTRUCTURES IN EARLY-TYPE GALAXY HALOS Davide Fiacconi1 , Piero Madau1,2,3 , Doug Potter1 , and Joachim Stadel1

arXiv:1602.03526v2 [astro-ph.GA] 6 May 2016

Accepted for publication in ApJ

ABSTRACT We present initial results from the “Ponos” zoom-in numerical simulations of dark matter substructures in massive ellipticals. Two very highly resolved dark matter halos with Mvir = 1.2 × 1013 M⊙ and Mvir = 6.5 × 1012 M⊙ and different (“violent” vs. “quiescent”) assembly histories have been simulated down to z = 0 in a ΛCDM cosmology with a total of 921,651,914 and 408,377,544 particles, respectively. Within the virial radius, the total mass fraction in self-bound Msub > 106 M⊙ subhalos at the present epoch is 15% for the violent host and 16.5% for the quiescent one. At z = 0.7, these fractions increase to 19 and 33%, respectively, as more recently accreted satellites are less prone to tidal destruction. In projection, the average fraction of surface mass density in substructure at a distance of R/Rvir = 0.02 (∼ 5–10 kpc) from the two halo centers ranges from 0.6% to & 2%, significantly higher than measured in simulations of Milky Way-sized halos. The contribution of subhalos with Msub < 109 M⊙ to the projected mass fraction is between one fifth and one third of the total, with the smallest share found in the quiescent host. We assess the impact of baryonic effects via twin, lower-resolution hydrodynamical simulations that include metallicity-dependent gas cooling, star formation, and a delayed-radiative-cooling scheme for supernova feedback. Baryonic contraction produces a super-isothermal total density profile and increases the number of massive subhalos in the inner regions of the main host. The host density profiles and projected subhalo mass fractions appear to be broadly consistent with observations of gravitational lenses. Keywords: cosmology: theory – dark matter – galaxies: halos – gravitational lensing: strong – methods: numerical 1. INTRODUCTION

Precision measurements across cosmic time have made the cold dark matter plus cosmological constant (ΛCDM) paradigm one of the pillars of our current understanding of the origin and evolution of structures in the universe. This model has proven to be remarkably successful at matching observations from scales that span the horizon length (e.g., Planck Collaboration et al. 2015) all the way down to the scales probed by the Lyman-α forest (e.g., Viel et al. 2013). A vast zoo of non-baryonic dark matter candidates has been proposed over the last three decades to reproduce the wealth of cosmological/astrophysical data (e.g., Feng 2010). The standard theory of structure formation requires dark matter to be cold, i.e. made of particles that become non-relativistic well before the matter domination era, and therefore clump on all scales. It is precisely on the smallest subgalactic scales that there have been persistent observational challenges to the cold, collisionless dark matter expectations. These “small-scale controversies”, predominantly found in the abundances and density profiles of dark matter-dominated dwarfs in the local universe, may simply stem from a poor understanding of the baryonic processes involved in galaxy formation (e.g. Pontzen & Governato 2014; Madau et al. 2014; Weinberg et al. 2015). They may alternatively indicate the need [email protected] 1 Center for Theoretical Astrophysics and Cosmology, Institute for Computational Science, University of Zurich, Winterthurerstrasse 190, CH-8057 Z¨ urich, Switzerland 2 Institute for Astronomy, ETH Zurich, CH-8093 Z¨ urich, Switzerland 3 Department of Astronomy & Astrophysics, University of California, 1156 High Street, Santa Cruz, CA 95064, USA

for more complex physics in the dark sector, and many modifications of the properties of the dark matter particle have been proposed to suppress small-scale power and alleviate some of these problems, including warm dark matter (WDM; Bode et al. 2001), self-interacting dark matter (Spergel & Steinhardt 2000), fuzzy dark matter (Hu et al. 2000), and superWIMPS (Cembranos et al. 2005). There is, however, no consensus on how serious these problems really are for ΛCDM, and detailed testing of the standard paradigm on small scales remains one of the most pressing issues in cosmology. Numerical simulations in ΛCDM have shown a rich spectrum of substructures in galaxy halos, the fossil remnants of a hierarchical merging process that is never complete (Moore et al. 1999; Klypin et al. 1999; Diemand et al. 2007, 2008; Springel et al. 2008; Stadel et al. 2009). Metcalf & Madau (2001) first showed that self-bound subhalos would have a readily detectable gravitational lensing effect if they were as numerous as predicted in ΛCDM. In particular, small fluctuations in the galaxy-scale lensing potential caused by substructures should result in measurable perturbations in the relative magnifications of quadruply-lensed quasar images. Discrepancies between the observed flux ratios and those predicted by a smooth lens model (“flux ratio anomalies”) have indeed been found to be common in quasar lenses (e.g., Mao & Schneider 1998; Chiba 2002; Metcalf & Zhao 2002), and can be explained if a large fraction (∼ 2%) of the projected mass at the Einstein radius is in the form of local substructures (Dalal & Kochanek 2002). To date, there is still no general agreement on whether semi-analytic models or N -body simulations in ΛCDM predict enough substructures to explain the frequency

2

Fiacconi et al.

of lens anomalies in currently available samples (e.g., Bradaˇc et al. 2004; Macci` o et al. 2006; Amara et al. 2006; Xu et al. 2009; Chen et al. 2011; Metcalf & Amara 2012; Xu et al. 2015). The discrepancy largely arises because of the small number of subhalos anticipated to survive near the radii where images form (typically around 5–10 kpc in projection). Somewhat surprisingly, most numerical simulation work has focused on present-day Milky Way-sized halos as the lens system. And yet, most gravitational lens galaxies are known to be early-type massive galaxies in relatively low density environments (e.g. Kochanek et al. 2000), while the average deflector redshift of the multiply-imaged quasars in the Cosmic Lens All-Sky Survey (CLASS) is hzlens i ≈ 0.6 (Browne et al. 2003). The abundance of substructures in a host is set by the competition between tidal disruption and new accretion. More massive hosts (as well as hosts at higher redshifts) are then expected to be more clumpy because their subhalos have been accreted more recently and managed to survive tidal destruction (Zentner et al. 2005; Klypin et al. 2011). In order to refine ΛCDM predictions for substructure lensing, we have initiated a program, dubbed “Ponos”, of very high resolution N -body and hydrodynamic cosmological simulations of early-type massive galaxies. In this Paper we present initial results on halo substructures from new “zoom-in” simulations of two ∼ 1013 M⊙ systems, each having between 155 and 286 million dark matter particles within the virial radius Rvir (defined here as the radius enclosing a mean matter density ∆ ρc , where ρc is the critical density and ∆ is the redshift-dependent virial overdensity, Bryan & Norman 1998). On these halo mass scales, the global fraction of E/S0 galaxies is estimated to be about 50% (Wilman & Erwin 2012). We focus on very high resolution calculations of a limited number of halos to study in detail the relative impact of light and heavy subhalos in strong lensing produced by elliptical galaxies, while we plan to extend this study to a larger halo sample. The two target halos are part of the AGORA Simulation Project4 (Kim et al. 2014). 2. NUMERICAL SIMULATIONS The Ponos collisionless simulations have been performed with the Tree/N -body code Pkdgrav3. Gravity was computed using a fast multipole method (Dehnen 2000, 2002) based on a 5th -order reduced expansion for faster and more accurate force calculation, and a multipole based Ewald summation technique for periodic boundary conditions (Stadel et al. 2009; Stadel 2013). An opening angle of θ = 0.55 was adopted for the gravity tree. Time-steps were assigned individually to√each particle based on the local dynamical time ∆t ≤ ξ/ Gρ, where ξ = 0.03 is an accuracy parameter and ρ is the density enclosed within the particle orbit (Zemp et al. 2007). This time-stepping strategy is both faster and more adaptive p and accurate than the conventional timestep ∆t ≤ ξ ǫg /|a|, where |a| is the magnitude of the local acceleration and ǫg is the gravitational softening. Initial conditions (ICs) were generated with the Music5 package (Hahn & Abel 2011), and are based upon a Wilkinson Microwave Anisotropy Probe cosmology with 4 5

sites.google.com/site/santacruzcomparisonproject/. www.phys.ethz.ch/~ hahn/MUSIC.

ΩM = 0.272, ΩΛ = 0.728, Ωb = 0.0455, σ8 = 0.807, ns = 0.961, and H0 = 70.2 km s−1 Mpc−1 (Hinshaw et al. 2013). The two target halos were selected from a low resolution pathfinder simulation of a 85.5 comoving Mpc box, and are characterized by very different assembly histories (Kim et al. 2014). One (“PonosQ”, Mvir = 6.5 × 1012 M⊙ ) has a relatively quiescent early merger history (i.e. few major mergers between z = 4 and 2, as defined by Kim et al. 2014), while the other (“PonosV”, Mvir = 1.2 × 1013 M⊙ ) is characterized by a more violent early merger history (i.e. many mergers between z = 4 and 2) 6 . A strong isolation criterion was imposed on the quiescent halo, one in which its 3Rvir radius circle does not intersect the 3Rvir radius circle of any halo with half or more of its mass at z = 0. A relaxed criterion was used for the most massive violent halo: 2Rvir circle instead of 3Rvir . Higher-resolution simulations were performed on new ICs re-centered on each of the target halos, for a total of six further nested spatial refinements (each by a factor of 2), corresponding to an effective resolution of 81923 particles with mass mp = 4.2 × 104 M⊙ . The final highest-resolution region is the convex hull Lagrangian volume that contains all the particles that fall within 3Rvir of the target halo at z = 0, sufficiently large to include all the structures that merge with the main host or have a significant impact on its evolution. The high resolution regions of PonosV and PonosQ were sampled with 890 and 388 million particles, respectively, and evolved with a force resolution ǫg = 210 pc (fixed in physical units after z = 9 and equal to 1/50 of the initial interparticle separation) from z = 100 to z = 0. Contamination within Rvir from coarse-level particles was found to be well below 0.1% in both number of particles and mass fraction down to the present epoch. To check for numerical convergence, both halos were resimulated at a factor two lower spatial resolution. Some basic properties of the simulated target halos are listed in Table 1. 2.1. Halo finding, merger-trees, and density profiles Dark matter field halos and subhalos were identified with the amiga halo finder (AHF; Gill et al. 2004; Knollmann & Knebe 2009). AHF employs a recursively refined grid to locate local overdensities in the density field. The identified density peaks are treated as centers of prospective halos, subhalos, etc. An iterative procedure is then used to collect those particles that are bound to those centers while removing those that are unbound. A bound clump of mass m at a distance d from the center of a more massive, M > m, halo is a subhalo of M if d < R+sr, where r and R are the virial (or tidal) radii of m and M , respectively, and s = 0.75 is a superposition parameter. Note that, according to this criterion, clumps whose centers lie just outside R can still be subhalos of M . We have checked the reliability of the identification of subhalos by comparing the AHF results with those produced by Rockstar (Behroozi et al. 2013). We find nice agreement between the two halo finders in terms of subhalo mass functions, positions, and distributions of 6 We note that, despite the naming of the simulations, PonosQ has a more active assembly history at low redshift than PonosV, as discussed in Section 2.1.

3

CDM Substructures in Early-Type Galaxies Table 1 Properties of the simulated elliptical halos. Name PonosV PonosQ

Merger History (2 < z < 4) violent quiescent

zf

Mvir (M⊙ )

Rvir (kpc)

Vhost (km s−1 )

rs (kpc)

c

Nvir

Nsub

1.10 0.74

1.2 × 1013 6.5 × 1012

600.5 489.6

348.2 265.3

50.9 55.0

11.8 9.06

285,558,928 154,489,668

48,681 26,521

Note. — Columns 3, 4, 5, 6, 7, 8, 9, and 10 give the formation redshift, present-day virial mass, virial radius, maximum circular velocity, scale radius and concentration parameter of the best-fit Navarro et al. (1997) (NFW) spherically-averaged density profile, the number of dark matter particles and the number of self-bound subhalos within the virial radius of the target halos, respectively.

Figure 2. Mass growth history of both PonosQ (thin dotted line) and PonosV (thin dashed line) as a function of redshift z. Red thick lines show the mean (solid) and median (dashed) mass growth of a 1013 M⊙ halo at z = 0 as determined by Fakhouri et al. (2010). The red shaded region shows the estimate of the standard deviation on the relation. The comparison reveals that both PonosQ and PonosV have a fairly representative (though different) assembly history, as they are both compatible with the mean growth of a 1013 M⊙ halo.

Figure 1. The “merger trees” of our simulated elliptical-sized halos, PonosQ (top) and PonosV (bottom). The main halo at the present day is plotted in the bottom-left corner, and all its progenitors (and also their histories) are depicted backward in time. The size of the symbols scales with halo mass, and the “main branch” is colored in green. Only mergers with mass ratios q ≡ M/m < 8 are shown, and all major mergers with mass ratios q < 4 are highlighted with the red color. The bottom portion of each panel shows the mass growth history of the main host, with major mergers marked by the vertical dashed lines.

virial (tidal) radii7 . Figure 1 shows the merger trees of the two main hosts, including only mergers with mass ratios q ≡ M/m < 8. We traced subhalos backward and forward in time, matching particles with the same ID inside each system for all snapshots between redshift 10 and 0. The main progenitor of a parent halo in a given snapshot was identified as the halo in the previous snapshot that maximizes 7 By default, Rockstar does not output the tidal radius of a halo. We extended the code by computing the tidal radius as the distance from the halo center of the farthest bound particle.

√ the ratio Ns / N n, where N and n are the number of particles in the parent and progenitor, respectively, and Ns is the number of particles they share (e.g., Fiacconi et al. 2015). As already mentioned above, PonosV and PonosQ have rather different assembly histories. While more massive, PonosV forms earlier (its formation redshift, defined as the redshift at which the halo had already assembled half of its mass, is zf = 1.1) and undergoes four major mergers with q < 4 at z . 3.5. PonosQ forms a bit later (zf = 0.74) and has only one major merger (with q = 3) at z ≃ 0.9. After z = 3, PonosV mostly grows within a main dark matter filament that routes the mass accretion toward the halo. On the other hand, PonosQ assembles around z ∼ 1 − 2 after the merger of a net of short dark matter filaments converging to the position of the halo. The latter grows quickly after z ∼ 1 and gathers most of the nearby smaller halos, cleaning out the surroundings consistently with the more stringent isolation criterion adopted for the selection. Figure 2 compares the mass growth of PonosQ and PonosV with the mean and median growth of a 1013 M⊙ halo at z = 0. These quantities have been calculated from the z-dependent mass accretion rates as inferred

4

Fiacconi et al.

from the data of both the Millennium and MillenniumII simulation (see equation (2) of Fakhouri et al. 2010). PonosV follows the average growth after z ∼ 3, while PonosQ mostly stays slightly below the curve (since it is lighter), but then it gets closer to the fit after z ∼ 0.9, as it assembles later than PonosV. Although we do not know the dispersion on the fitted relations, we can roughly estimate p the standard deviation σ from the relation |m−µ|/σ ≤ 3/5, where µ and m are the mean and median, respectively (Basu & DasGupta 1997). Figure 2 shows that both PonosQ and PonosV have mass growths compatible with the average at the present-day reference mass 1013 M⊙ ; this generally qualifies both of them as “representative” halos for their mass scale, though they behave rather differently in detail and that influences the properties of subhalos (see Section 3 and 4). Figure 3 shows the present-day spherically-averaged density profiles, their logarithmic slope d ln ρ/d ln r, and the fractional deviations from the best-fit NFW profiles for the two target halos. The NFW formula with scale radii of 51 and 55 kpc and concentration parameters of 11.8 and 9.1 for PonosV and PonosQ, respectively, produces a reasonable approximation to the density profiles down to the convergence radius of each simulation. At all radii between 1 and 500 kpc, deviations from the best-fit NFW matter densities are typically less than 20%. 3. SUBSTRUCTURE ABUNDANCE

We counted all substructures down to a minimum of 20 bound particles, corresponding to a minimum subhalo mass of Mres = 8.4 × 105 M⊙ . Our selection results in a sample of 48,681 individual subhalos at z = 0 in PonosV, and 26,521 subhalos in PonosQ. In Figure 4 we present the cumulative maximum circular velocity function (normalized to the maximum circular velocity of the host) of the subhalo population of each Ponos host. These are compared with the subhalo velocity functions of the Via Lactea II simulation (Diemand et al. 2008). Over the interval 0.02 < x < 0.08, all velocity functions can be approximated by the power-law 1/2

N (> x) = 3.6 × 10−3 x−3 Vhost ,

(1)

where x ≡ Vsub /Vhost and Vhost is given in units of km s−1 (c.f., Klypin et al. 2011). Residuals from the power-law regression line are typically of order 20%. At a given Vsub /Vhost , the number N (> x) of subhalos is 56% higher in PonosV compared to Via Lactea, and 86% higher in PonosQ. This is because more massive ellipticals are dynamically younger than Milky Way-sized galaxies, and have more subhalos that manage to survive tidal destruction. Note the relatively large deviations above x & 0.15 that are present in all three hosts compared to the fitting formula. Figure 5 shows the cumulative mass fraction in selfbound substructures of PonosQ and PonosV as a function of Msub , Z Msub 1 dN fsub (< Msub ) = dm, (2) m Mvir Mres dm at three different redshifts, z = 0, 0.5, 0.7. Here, dN/dm is the subhalo mass function. At the present epoch, we measure a total mass fraction fsub = 15% in the violent host and fsub = 16.5% in the quiescent one. At z = 0.7,

Figure 3. Top: present-day density profiles of PonosV (blue circles) and PonosQ (green squares). The red solid and dashed lines are the best-fit NFW profiles. The density of PonosQ has been rescaled by a factor of 1/10 for clarity. The colored arrows indicates the regions within which numerical convergence is not achieved because of two-body relaxation (Power et al. 2003). Middle: logarithmic slope of the matter density profiles of the two target halos. Bottom: the residuals between the density profile and the best-fit NFW profile, as a function of radius. In all panels, the gray shading indicates the region inside two gravitational softening lengths.

these fractions increase to 19% and 33%, respectively, as more recently accreted subhalos are more likely to survive tidal stresses. Compared to PonosV, PonosQ forms later, is less concentrated, and undergoes the majority of its mergers at lower redshifts. As a consequence, there is less time for the orbital decay and mass loss of infalling satellites, and its substructure mass fraction is higher (e.g., Chen et al. 2011). The definition of virial radius, while formally meaningful, is nevertheless rather arbitrary. We have checked that companion systems just outside the virial radius (in the spherical shell Rvir < r < 2Rvir ) would make a relatively modest contribution to the “subhalo” mass fraction, corresponding to less than a 30% increase in fsub . This is not true, however, for subhalos in the shell R200 < r < Rvir . To facilitate comparison with pre200 vious simulations, we have computed fsub , the subhalo mass fraction within R200 , the radius with mean enclosed overdensity equal to 200 times the critical value8 . At 200 the present epoch, we derive fsub =6.9% in PonosV and 8 At z = 0, we measure R 200 = 444.4 and 359.6 kpc for PonosV and PonosQ, respectively.

CDM Substructures in Early-Type Galaxies

Figure 4. Cumulative subhalo abundance at z = 0 as a function of maximum subhalo circular velocity, Vsub ≡ p max{ GMsub (< r)/r}, in units of the maximum circular velocity of the main host, Vhost . We show results for PonosV and PonosQ and, for comparison, for the Via Lactea II simulation. The dotted line shows the best-fit power-law relation of Equation (1).

Figure 5. The cumulative substructure mass fraction of PonosQ (dashed lines) and PonosV (solid lines), as a function of subhalo mass Msub , at z =0 (blue), 0.5 (red), and 0.7 (green). 200 fsub =11.4% in PonosQ. These values are consistent with the mean mass fraction inferred at z = 0 in the six Milky Way-sized halos simulated at resolution level 2 as part 200 of the Aquarius project, hfsub i = (7.2 ± 2.3)%, as reported by Xu et al. (2009). However, subhalos within R200 account for only 40% (PonosV) and 57% (PonosQ) of the total subhalo mass fraction fsub within Rvir , the remaining substructure material being located between R200 and Rvir . According to Xu et al. (2009), the scatter among the

5

six Aquarius halos is much larger than the differences between halos at z = 0 and z = 0.6. By contrast, the two Ponos hosts have a larger mass fraction in substructure 200 at z = 0.7 than at present: fsub =12.1% in PonosV and 200 =21.2% in PonosQ. These values are twice as large fsub as measured at z = 0. These large variations are associated with the late accretion of a few, relatively massive, Msub & 1011 M⊙ satellites. At z = 0.7 there are two such systems in the outer halo of PonosV, and four in PonosQ. By z = 0, only one of them survives tidal stripping and remains above 1011 M⊙ in PonosV, and two in PonosQ. We note, in passing, that late minor mergers with progenitors mass ratios q < 30 are fairly typical for descendants of mass Mvir ≈ 1013 M⊙ today: an analysis of the joint data set from the Millennium and Millennium-II simulations gives rates for such events that exceeds 1 Gyr−1 per halo at redshift 1 (Fakhouri et al. 2010). Over the interval 108 M⊙ < Msub < 1011 M⊙ , our measured subhalo mass function is well approximated by the power law −n Msub dN = N0 . (3) d ln Msub m0 At z = 0, we derive for PonosV a best-fit slope of n = 0.877 and an amplitude at the pivot mass, m0 = 108 M⊙ , of N0 = 847. At z = 0.7, the best-fit parameters are n = 0.915 and N0 = 868. A similar relation holds also for PonosQ, with slopes and amplitudes (n, N0 ) = (0.754, 362) at z = 0, and (n, N0 ) = (0.949, 344) at z = 0.7. In the redshift interval 0 ≤ z ≤ 0.7, the most massive subhalo reaches 8% of the mass of the host. The expected total mass fraction in self-bound substructure below our nominal resolution limit of Mres = 8.4 × 105 M⊙ is 1−n Mres N 0 m0 (4) fsub (< Mres ) = Mvir (1 − n) m0 (assuming a thermal free-streaming mass limit → 0). The above best-fit power laws yield between one and three percent of the virial mass in unresolved subhalos at z = 0, considerably smaller compared to the fractional mass in substructures that are already resolved in our simulations. The substructure mass fraction converges much more slowly at z = 0.7, however, as the slope of the subhalo mass function is closer to n = 1. Resolved subhalos do not trace the matter distribution of the host: tidal disruptions are most effective in the inner halo, leading to an antibias in the abundance profile of substructures relative to the smooth background (Diemand et al. 2007). At the present epoch, the subhalo mass fraction within the inner 10 kpc drops to 2 × 10−4 and 1.3 × 10−6 in the violent and quiescent host, respectively. We caution, however, about resolution effects (numerical “overmerging”) on this quantity. Projecting the subhalos relative to the total mass along the line-of-sight (see next section) will be less strongly influenced by the incompleteness of the inner regions of the host. 4. SUBSTRUCTURE SURFACE MASS DENSITY

Strong gravitational lensing occurs when the surface mass density along a given sightline exceeds a certain critical value, and is associated with high magnifications, multiple images, Einstein arcs and rings in the lens plane.

6

Fiacconi et al.

Figure 6. Projected dark matter density at z = 0, 0.5, 0.7 (from right to left) of PonosQ (top) and PonosV (bottom) in a slice of thickness 2Rvir centered on the target halos. The image brightness is proportional to the logarithm of the total dark matter surface density along the line-of-sight, ΣDM . The circles mark the virial radius Rvir .

Figure 6 shows the surface dark matter mass density at z = 0, 0.5, 0.7 of PonosQ and PonosV in a slice of thickness equal to 2Rvir centered on the target halos. As expected in ΛCDM, the two Ponos halos are teeming with self-bound subhalos on all resolved mass scales. To shed light on whether CDM substructures can account for the observed lensing flux ratio anomalies, we have measured in the simulations the surface mass fraction in subhalos within 20 different azimuthal annuli equally spaced in log(R/Rvir ) from −2 to 0, Fsub (R) =

Σsub (R) , ΣDM (R)

(5)

where R is the projected radius measured from the host center, and Σsub and ΣDM are the substructure and total dark matter projected mass densities, respectively. The surface density in substructures is computed by summing up all the subhalo bound mass that falls within the relevant annulus. To be specific, we adopt the following procedure to compute Fsub . First, we determine a lineof-sight by sampling uniformly the solid angle 4π. This is achieved by randomly drawing the azimuthal angle φ from a uniform distribution in the interval [0, 2π), and the co-latitudinal angle θ from the distribution P(θ) = sin(θ)/2. Then, we determine the line-of-sight through the versor n = (sin θ cos φ, sin θ sin φ, cos θ). We calculate the projected distance p of a particle in the plane perpendicular to n as R = |r|2 − (r · n)2 , where r is the

position vector of the particle from the centre of the halo, and we finally bin the positions of each particle, adding its mass to the proper radial bin. We repeat this for both the subhalos only and the host halo, and we finally calculate the ratio between the two to derive Fsub (R). Figure 7 shows the median subhalo surface mass fraction at z = 0, 0.5, 0.7 as a function of R/Rvir over 300 random projections for each of the two Ponos halos. The 68% scatter among the different projections is marked by the colored bands, and is larger at smaller radii as the area of the azimuthal annulus around R becomes progressively smaller in the inner regions. The value of Fsub is clearly sensitive to the assembly history of the host galaxy. The black dashed line depicts the results derived by Mao et al. (2004) at z = 0 from simulations of twelve halos of galactic, group, and cluster masses. Their surface mass fraction in substructures is somewhat lower than found here, perhaps because of resolution limitations – the number of particles per halo in their numerical investigations is about two orders of magnitudes smaller than here. The statistical study by Dalal & Kochanek (2002) of the anomalous flux ratios observed in a sample of seven lensed radio-loud quasars requires Fsub =0.6 to 7% (90% confidence, with a median of 2%) of the mass at the Einstein radius to be in substructures in order to reproduce the data. The black point in the figure indicates the Dalal & Kochanek (2002) constraint – for an assumed


the Aquarius suite of galactic halos and the Phoenix suite of cluster halos to those expected in massive ellipticals. They estimate a mean surface mass fraction of substructure at the Einstein radius that is three times bigger than in Milky Way-sized hosts, in better agreement with our findings. We note, however, that PonosQ has a mean subhalo surface mass fraction that is still a factor of ∼ 2 above the determinations by Xu et al. (2015). At z = 0, we find a mean substructure surface mass density around the Einstein radius of Σsub = 3.6 × 106 M⊙ kpc−2 in PonosV, and 8.9 × 106 M⊙ kpc−2 in PonosQ. At z = 0.7, we measure Σsub = 1.4 × 107 M⊙ kpc−2 and 1.0 × 107 M⊙ kpc−2 in PonosV and PonosQ, respectively. Let us denote with η the mean column number density of subhalos at projected halocentric distance R = 0.02 Rvir. Over the interval 108 M⊙ < Msub < 1011 M⊙ , the mean projected mass function can be written as −ℓ dη Msub = N0 . (6) d ln Msub m0

100 PonosV

Fsub(R)

10−1

10−2

100 PonosQ

Fsub(R)

10−1

10−2

z=0 z = 0.5 z = 0.7 0.25 R/Rvir, Mao et al.(2004) Dalal & Kochanek (2002)

10

−2

7

10 R/Rvir −1

100

Figure 7. Predicted substructure surface mass fraction at z = 0, 0.5, 0.7 in PonosV (top panel) and PonosQ (bottom panel). The panels show Fsub in azimuthal annuli as a function of the projected radius, R (in units of the virial radius), for all subhalos belonging to the corresponding host. The solid lines depict the median over 300 random projections for each of the two Ponos halos, while the colored bands mark the 68% scatter among the different projections. The black point indicates the median and 90% confidence level of the mass fraction in local substructure required to explain the flux ratio anomalies (for an assumed Einstein radius of 0.02 Rvir ∼5–10 kpc). The dashed black curve shows the results derived by Mao et al. (2004) at z = 0 from low resolution simulations of twelve halos of galactic, group, and cluster masses.

Einstein radius of 0.02 Rvir ∼ 5–10 kpc, which appears entirely consistent with the expectations for ΛCDM. At the Einstein radius, the mean surface mass fraction in substructure for the two Ponos halos ranges from 0.6% to 2.3%. In contrast with PonosV, there is hardly any redshift dependence in the value of this quantity for PonosQ. These projected mass fractions are considerably higher than the mean value of 0.2% (0.01 ≤ Fsub ≤ 0.7%) measured at the Einstein radius in the Aquarius simulations of Milky Way-sized halos (Xu et al. 2009). Recognizing that the substructure abundance in group-sized halo may be larger than in less massive Milky Way-sized hosts, Xu et al. (2015) recently rescaled the subhalo populations of

At z = 0.7, we measure in PonosV a best-fit slope of ℓ = 0.85 and a column number density at the pivot mass, m0 = 108 M⊙ , of N0 = 0.006 kpc−2 . In PonosQ we find a similar normalization but a steeper slope, ℓ = 1.09. At z = 0, we measure ℓ = 0.98 in PonosV and ℓ = 0.84 in PonosQ. Most subhalos at these image positions only appear along the line-of-sight because of projection effects. In Figure 8 we plot the mean and median surface mass fraction per decade of subhalo mass, ∆Fsub , again at projected halocentric distance R = 0.02 Rvir . As expected from a subhalo projected mass function dη/d ln Msub ∝ −ℓ Msub with ℓ close to 1, this distribution is relatively flat, with roughly equal contributions per decade of mass above the “completeness” mass scale of Msub = 107 M⊙ . The contribution of subhalos with Msub < 109 M⊙ to the projected substructure mass fraction is between one fifth and one third of the total, with the smallest share found in the quiescent host. Massive subhalos with Msub & 1010 M⊙ only survive in the outer, hri & 160 kpc, regions of their hosts because of tidal destruction. Their presence in the projected central ∼ 10 kpc is then typically associated with chance alignment. This explains the large disparity between the mean and median ∆Fsub in the largest mass bins. The insets in the figure show the broad probability density function (PDF) of the total Fsub around 0.02Rvir , for the 300 lines of sight. At z = 0.7, the median Fsub is about 1%, and the probability of observing values of Fsub > 0.05 can be as large as 12% (PonosQ). 5. BARYONIC CONTRACTION

The most severe limitation of our study is that the N -body very high resolution simulations used here include only dark matter. Standard dark matter halos are poor lenses because their central cusps (ρ ∝ r−1 ) are too shallow. In the adopted cosmology, and for a typical lens geometry with zlens = 0.7 and zsource = 2, the critical surface density for multiple imaging is Σcrit = 2.2 × 109 M⊙ kpc−2 . At the same redshift and R = 0.02 Rvir, our Ponos lens halos have Σ/Σcr = 0.2 − 0.4, i.e. are subcritical and unable to produce multiple images despite

8

Fiacconi et al.

Figure 8. Projected substructure mass fraction at R = 0.02 Rvir in each mass decade, for PonosV (left panel) and PonosQ (right panel). The individual horizontal bars depict the mean (solid lines) and median (dashed lines) over 300 random projections for each host. The three different redshifts inspected are listed in the left) panel. The insets show the broad probability density function of the total Fsub , for the 300 lines of sight, at the three redshifts. Median values are indicated by the vertical dotted lines.

having maximum circular velocities close to 240 km s−1 (PonosQ) and 340 km s−1 (PonosV). Baryon cooling and condensation is expected to strongly affect the inner density profiles of halos, although predictions for these effects are far less certain than forecasting the distribution of dark matter in the purely dissipationless regime. As the baryons cool, they drag some of the dark matter inward and may even dominate the mass within one Einstein radius, thereby converting a sub-critical dark matter halo into one capable of producing multiple images. Rather than putting in baryons “by hand”, for example by implanting an idealized model galaxy into a dark matter halo extracted from a collisionless simulation (Amara et al. 2006) or by modeling the main lens halo as a singular isothermal ellipsoid (e.g., Metcalf & Amara 2012; Xu et al. 2015), we have followed here a different approach and run a lower resolution version of PonosQ (“PonosQH”) with hydrodynamics using the TreeSPH code Gasoline (Wadsley et al. 2004). The code employs a subgrid model for the turbulent mixing of metals and energy following Shen et al. (2010). The ICs for PonosQH were initialized as for the collisionless run. The highresolution region contains 6 million dark matter particles and an equal number of gas particles, for a dark matter and (initial) gas particle mass of mDM = 2.3 × 106 M⊙ and mSPH = 4.5 × 105 M⊙ , respectively. This mass resolution is very similar to that of the “small group mass” halo in the FIRE simulation suite (Hopkins et al. 2014). The gravitational softening length was fixed to 785 pc (physical) for the dark matter, and to 501 pc for the gas. In high-density regions the gas smoothing length was allowed to shrink to 10% of the softening to ensure that hydrodynamic forces are resolved on 100 pc scales. A

non-thermal pressure floor (Agertz et al. 2009; Roˇskar et al. 2015) was applied to stabilize scales of order the gravitational softening against gravitational collapse and avoid artificial fragmentation. The simulation includes a non-equilibrium primordial chemistry network for atomic H and He, Compton cooling off the cosmic microwave background, and precomputed tabulated metal-line cooling rates from the photoionization code Cloudy (Ferland et al. 1998). A spatially-uniform, redshift-dependent cosmic UV background (Haardt & Madau 2012) modifies the ionization and excitation state of the gas, photoionizing away abundant metal ions and reducing the cooling efficiency. Star formation proceeds at a rate dρ⋆ /dt = 0.05(ρgas /tdyn ) ∝ 3/2 ρgas (i.e., locally enforcing the Schmidt law), where ρ⋆ and ρgas are the stellar and gas densities, and tdyn is the local dynamical time. Star particles form in cold gas (i.e. temperature below 104 K) that reaches a density threshold of 20 atoms cm−3 , and are created stochastically with an initial mass m⋆ = 1.35 × 105 M⊙ distributed following a Kroupa (2001) initial mass function. They inject energy, mass, and metals back into the interstellar medium (ISM) through Type Ia and Type II supernovae (SNe) and stellar winds, following the prescriptions of Stinson et al. (2006). A “delayed radiative cooling” scheme for Type II SN feedback was adopted, a simplified algorithm designed to extend the Sedov-Taylor phase of SN remnants and mimic the effect of energy deposited in the local ISM by multiple, clustered sources of mechanical luminosity (Stinson et al. 2006). While this approach has been found to be key in reproducing the properties of dwarfs (Governato et al. 2010; Madau et al. 2014; Shen et al. 2014), late-type spirals (Guedes et al. 2011), and the circumgalactic medium of z ∼ 3 galaxies (Shen et al.


9

Figure 9. The gas and stellar properties of PonosQH at different redshifts. From top to bottom: gas surface-density, gas temperature in a slice through the center of the host, mock images in the U , V , and J filter bands. The top three panels encompass a physical scale equal to the virial radius, while the bottom panels shows images of the inner 40 kpc (physical).

2013), many authors have recently stressed the importance of correctly accounting for the entire momentum budget of stellar feedback, including momentum injection by radiation pressure, stellar winds, and clustered SN explosions (see, e.g. Agertz et al. 2013; Hopkins et al. 2014; Keller et al. 2014; Kim & Ostriker 2015). It has been shown by Agertz et al. (2013) that simulations with maximal momentum injection suppress star formation to a similar degree than found in runs that, like ours, adopt adiabatic thermal feedback. We note, however, that our stellar feedback model is not strictly speaking adiabatic (since SN-heated gas particles exchange thermal energy with the ambient medium via turbulent mixing and therefore can cool “indirectly”), and that we did not include any feedback from a central active galactic nucleus (AGN).

Figure 9 shows gas surface-density and temperature maps at four different redshifts, together with mock images in the U , V , and J filter bands. The top three panels encompass a physical scale equal to the virial radius, while the bottom panels are images of the inner 40 kpc (physical). Each star particle was assigned a luminosity from tables of mass-to-light ratios based on the isochrones and synthetic stellar populations of Bressan et al. (2012). At the present epoch, PonosQH has an early-type morphology, with U −V = 1.5 and V −J = 0.9 colors that are typical of red sequence galaxies. Its present-day stellar mass within 20 kpc from the center, M⋆ = 2.3 × 1011 M⊙ , implies a star formation efficiency, M⋆ /Mvir = 0.035, which is slightly higher than the value inferred for massive halos by Kravtsov et al. (2014), but compatible within 2σ. The surface brightness profile fol-

10

Fiacconi et al.

Figure 10. The z = 0 rotation curve of PonosQH, for all components (black solid line), stars (dotted line), gas (dash-dotted line), and dark matter (dashed line). For comparison, the dark matter rotation curve (blue solid line) obtained in the purely collisionless simulation of the same object is also shown. The gray shading indicates the region inside two gravitational softening lengths of the baryonic component. The bump at ∼ 35 kpc from the center is caused by a massive satellite system.

lows a de Vaucouleurs law with a V -band effective radius Re,V = 3.6 kpc, while the half-light radius measured directly from the light profile is R1/2 = 3.1 kpc. The baryonic content of PonosQH is close to the universal baryon fraction. The mass of cold (T < 104 K) gas per unit Kband luminosity, M (HI)/LK = 0.03, places this system towards the high-end of the distribution observed in the ATLAS survey of nearby early-type galaxies (Serra et al. 2012). We now focus on the the impact of baryonic infall on the inner density profile and projected substructure mass fraction. p In Figure 10 we show the z = 0 rotation curve, Vc (r) = GM (< r)/r, for all components (gas, stars, and dark matter) separately. The dark matter circular velocity curve obtained in the collisionless simulation of the same object is plotted for comparison. Baryons dominate within the inner 20 kpc, and the ensuing rotation curve remains relatively flat from several kpc to several tens of kpc. The bump at ∼ 35 kpc is due to the presence of a massive, ∼ 1.5 × 1011 M⊙ , satellite. The total mass within 10 kpc has increased by a factor of three relative to the purely collisionless simulation. The effective power-law slope of the total (luminous plus dark; ρ ∝ r−γ ) mass distribution in the range 1.5-40 kpc is γ = 2.19 ± 0.07, i.e. marginally steeper than isothermal (γ = 2). This value appears to be consistent with the results of the Sloan Lens ACS Survey (SLACS) of 73 early-type galaxies with 0.08 < zlens < 0.5 and stellar masses above 1011 M⊙ : hγi = 2.078 ± 0.027 with an intrinsic scatter of 0.16 ± 0.02 (Auger et al. 2010). Recent modeling of the mass density profiles of early-type galaxies also yields super-isothermal central slopes, with hγi = 2.15 ± 0.04 (Chae et al. 2014) and hγi = 2.19 ± 0.03 (Cappellari et al. 2015). The steeper inner cusp is clearly seen in the density profile plot (Figure 11). In response to the slow addition of baryons to the center, dark matter may be pulled in-

Figure 11. The z = 0 density profile in PonosQH for all components (black solid line) and dark matter (dashed line). For comparison, the dark matter density profile (red dotted line) obtained in the purely collisionless simulation of the same object is also shown. The gray shading indicates the region inside two gravitational softening lengths of the baryonic component. We do not observe any significant evidence for adiabatic contraction or expansion in the dark matter at the center of the halo in the hydrodynamical run.

Figure 12. The z = 0 average density profiles of massive subhalos. Solid lines: PonosQH. Dashed lines: PonosQ. The difference between the two sets shows the impact of baryonic contraction. The red, black, and blue colors refer to the subhalo mass intervals 109 < Msub /M⊙ < 1010 , 1010 < Msub /M⊙ < 1011 , and Msub /M⊙ > 1011 , respectively. For reference, the dotted line shows an isothermal, ∝ r −2 , profile. The gray shading indicates the region inside two gravitational softening lengths of the baryonic component.

wards through a process known as adiabatic contraction (Blumenthal et al. 1986; Gnedin et al. 2004; Pillepich et al. 2014). Other processes may cause the dark matter to expand instead, such as the transfer of orbital energy via dynamical friction following dry minor mergers (El-Zant et al. 2001) and rapid gravitational poten-


11

Figure 13. Impact of baryonic contraction on the radial distribution of massive subhalos. Solid line: radial distribution of the cumulative mass, ∆Msub , in all PonosQH subhalos with Msub > 5 × 108 M⊙ . Dashed line: same for PonosQ. The distance to the center of the main halo is normalized to the virial radius. The subhalo mass is normalized to the total mass in subhalos > 5 × 108 M⊙ identified in the hydro and collisionless simulations. The left, middle, and right panels show the distribution at three different redshifts.

tial fluctuations tied to efficient SN feedback (Pontzen & Governato 2012; Madau et al. 2014). At the presentepoch, PonosQH shows little evidence for strong dark matter contraction or expansion. Using a generalized NFW model, we find an inner dark matter density slope, γDM = 1.1 ± 0.1, which is consistent with unity. The mean dark matter mass fraction projected within a cylinder of radius equal to the Einstein radius (0.02 Rvir) is FDM (< Rvir ) = 0.38 ± 0.01, in agreement with the average value, 0.38 ± 0.07, measured in the SLACS lens sample by Bolton et al. (2008). Baryonic contraction acts both for the host halo and for its subhalos. At z = 0, we identify in PonosQH more than 400 subhalos, 285 of which are above 108 M⊙ (the minimum subhalo mass corresponding to 20 bound dark matter particles is now Mres = 4.6 × 107 M⊙ ). Figure 12 shows the z = 0 density profiles of massive, Msub > 109 M⊙ subhalos averaged over different mass intervals. All subhalos above Msub > 1010 M⊙ form stars and are therefore “luminous”, while half of those in the range 109 < Msub /M⊙ < 1010 are dark. The difference between PonosQH and the collisionless PonosQ is clear. Baryonic infall increases the total mass within, say, 2 kpc from the satellite center by a factor that ranges from 2 (109 < Msub /M⊙ < 1010 ) to 14 (Msub /M⊙ > 1011 ) relative to the purely collisionless simulation. As a consequence, the mean maximum circular velocity of subhalos in the mass range 109 < Msub /M⊙ < 1010 rises from 28 km s−1 to 49 km s−1, and from 83 km s−1 to 122 km s−1 in the case of 1010 < Msub /M⊙ < 1011 subhalos. Similar effects of baryonic contraction in massive satellites of Milky Way-sized halos have also been recently reported by Zhu et al. (2015). And while the slope of the internal density profile of subhalos may have little effect on the frequency of flux anomalies (Metcalf & Amara 2012), baryon cooling and condensation within massive subhalos will make them more resilient to tidal disruptions (e.g. Macci`o et al. 2006). In Figure 13 we compare the radial distribution of the cumulative subhalo mass in PonosQH and PonosQ. Only subhalos with Msub > 5 × 108 M⊙ are included in the

comparison. The distance to the center of the main halo is normalized to the virial radius (slightly larger in PonosQH), while the subhalo mass is normalized to the total mass in all Msub > 5 × 108 M⊙ subhalos identified in the hydro and collisionless simulations. The left, middle, and right panels show subhalos at three different redshifts. The substructure mass in the inner regions is consistently larger in PonosQH. This notwithstanding the fact that the increased mass concentration at the center of the main host induces stronger tidal forces that can potentially destroy subhalos. At z = 0.7, for example, there are 4 subhalos more massive than 109 M⊙ within 30 kpc of the center of PonosQH, compared to none in PonosQ. At z = 0.5, the innermost Msub > 1010 M⊙ satellite in PonosQH is at 17 kpc from the center, versus 46 kpc for PonosQ. The effect of baryon cooling and condensation on the predicted projected mass fraction is depicted in Figure 14. We have modified the numerator and denominator of Equation (5) to account for the impact of baryons on the main halo and subhalos as Σsub (R) Fsub (R) = . (7) ΣDM (R) + Σbar (R) Here, Σsub is the substructure projected mass (baryons + dark matter) density, while ΣDM and Σbar are the total dark matter and baryon surface densities of the main host, all measured in PonosQH. PonosQH does not have enough resolution to generate a realistic substructure population at masses Msub . 5 × 108 M⊙ , and so this figure should be interpreted as the surface mass fraction in substructure at intermediate and large mass scales only. As in Figure 7, we have plotted the median substructure mass fraction rather than the mean, for better comparison with the Dalal & Kochanek (2002) constraints. Note how, even without small-scale power, the predicted median at the Einstein radius and z & 0.5 is consistent with the data. 6. DISCUSSION AND CONCLUSIONS

In this Paper, we have used new very high resolution N -body and hydrodynamical cosmological simulations of

12

Fiacconi et al. 100 PonosQH

Fsub(R)

10−1

10−2 z=0 z = 0.5 z = 0.7 0.25 R/Rvir, Mao et al.(2004) Dalal & Kochanek (2002)

10

−2

10 R/Rvir −1

100

Figure 14. Impact of baryonic contraction on the substructure surface mass fraction. Here, we have measured the total (baryons + dark matter) substructure surface density in the hydro simulation PonosQH. Legend is the same as Figure 7. Because PonosQH only resolve substructures at intermediate and large mass scales, the sightline-to-sightline scatter has increased compared to Figure 7.

early-type massive galaxies to refine ΛCDM predictions for substructure gravitational lensing. It has long been known that, if dark matter substructure can survive to constitute a percent or more of the surface mass density at small impact parameters, it will cause detectable anomalies in the magnification ratios of multiply-imaged QSOs (Metcalf & Madau 2001). In our dissipationless simulations, we have found that the total mass fraction in self-bound subhalos increases from about 15% at the present epoch to 20-30% at redshift 0.7. In projection, the average fraction of surface mass density in substructure around the Einstein radius at z = 0.7 exceeds 2%. More massive hosts (as well as hosts at higher redshifts) are predicted to be dynamically younger and therefore more clumpy, as their subhalos are accreted more recently and tend to survive tidal destruction. Indeed, our two Ponos Mvir ∼ 1013 M⊙ halos have significantly higher, by as much as a factor of 10, mean projected substructure mass fractions at z = 0.7 than measured at the present-epoch in the Aquarius and Via Lactea simulations of Milky Way-sized systems. As a result, the frequency of flux ratio anomalies predicted by the richer substructure population of early-type galaxy halos increases noticeably compared to estimates based on Milky Way-sized halos (e.g., Amara et al. 2006; Macci`o et al. 2006; Xu et al. 2009). We have incorporated the effects of baryonic contraction on the host halo using twin, lower-resolution hydrodynamical simulations that include metallicitydependent gas cooling, a star formation recipe based on a high gas density threshold, and a delayed-radiativecooling scheme for feedback by SNe, but no AGN feedback. The inclusion of the baryonic component produces at the present epoch a red sequence galaxy with a superisothermal central slope in the total (luminous plus dark) mass distribution, a rotation curve that is relatively flat from several kpc to several tens of kpc, a projected dark

matter mass fraction inside the Einstein radius of 40%, and little evidence for strong dark matter contraction or expansion. Such properties of the matter density profile appear broadly consistent with observations of early-type galaxies (e.g., Bolton et al. 2008; Auger et al. 2010; Chae et al. 2014; Cappellari et al. 2015). Baryonic contraction increases the number of massive subhalos in the inner regions of the main host. Our ΛCDM simulations, both in the purely collisionless case or after accounting for the impact of baryonic contraction, appear at first sight to predict enough substructure to explain the frequency of lens anomalies in currently available samples (Dalal & Kochanek 2002). While this is indeed promising, numerical calculations of the lensing potential, deflection angles, and magnifications are needed to generate theoretical flux-ratio probability distributions for comparison with the observations. Two such studies have been performed recently. Metcalf & Amara (2012) created a large number of simulated lenses with finite source sizes, compared the predicted flux-ratio probability distributions in the presence of substructure to an observational sample of seven lenses, and found approximate consistency with ΛCDM N -body simulations. Xu et al. (2015) constructed samples of lens potentials by adding (rescaled to those expected in massive ellipticals) subhalo populations from the galaxyscale Aquarius and the cluster-scale Phoenix simulation suites, and matched the resulting flux ratio distributions to the best available sample of radio lenses. They reached the conclusion that CDM substructures cannot account for all the observed anomalies. Such detailed investigations are beyond the scope of this work, and we defer them to a future paper. Here, we only note that, according to Figure 9 in Metcalf & Amara (2012), a subhalo surface mass density of Σsub = 107 M⊙ kpc−2 would imply a 15% chance of a clear outlier in the general distribution of ∆θ and Rcusp values predicted in the absence of substructures, consistent with the frequency (one out of seven) actually observed9. A substructure surface mass density of 107 M⊙ kpc−2 is comparable to that measured in the Ponos halos. Contrary to the model of Metcalf & Amara (2012), however, most of this projected mass density is associated with subhalos more massive than Msub = 109 M⊙ . Which brings up the next point. Magnifications, which depend on the second spatial derivative of the lensing potential, are affected about equally by all mass scales provided the Einstein radius of the deflector is larger than the size of the background source. The constraints provided by the frequency of flux anomalies cannot therefore discriminate between different substructure mass scales. On the other hand, distinguishing CDM from keV-mass WDM – in which the free-streaming cutoff occurs on dwarf galaxy scales – requires measuring the subhalo 9 Here, ∆θ and R cusp are two parameters that characterize 4-image lenses. Sources near a cusp in the caustic produce “cusp” configurations with three of the images (“triplet”) lying close together on one side of the lens galaxy. The parameter ∆θ is the angular separation between the close triplet, and is small when the source is near one of the cusps in the caustic. In any smooth lensing potential, the three close images satisfy an asymptotic magnification relation (the “cusp-caustic relation”) Rcusp = |µA + µB + µC |/(|µA | + |µB | + |µC |) → 0, with the total absolute magnification |µA | + |µB | + |µC | → ∞ (e.g. Keeton et al. 2003).

CDM Substructures in Early-Type Galaxies mass function well below a mass of ∼ 109 M⊙ (e.g., Li et al. 2015). We find that, in the Ponos hosts, the contribution of subhalos with Msub < 109 M⊙ to the total projected mass fraction is sub-dominant, between one fifth and one third of the total. This fact highlights the potential importance of dwarf-sized systems in the anomalous flux ratio problem. Small dark substructures may actually play a lesser role in causing flux anomalies than the more massive subhalos, and some of these massive perturbers may actually contain visible dwarf galaxies. Indeed, three out of the seven gravitational lens systems used by Dalal & Kochanek (2002) to statistically detect dark substructure around early-type galaxies show evidence for additional mass structure in the form of luminous dwarf satellites (see, e.g., McKean et al. 2007, and references therein). In the case of the lens systems B2045+265 and MG 2016+112, the dwarf perturber constitutes about 1% of the total projected lens mass, which seems consistent with the predictions of our Ponos simulations. Moreover, other components of the lensing galaxy, like an edge-on disk, may also account for some anomalies, as suggested by recent observations of CLASS B1555+375 (Hsueh et al. 2016). It remains therefore unclear how effective a probe of the matter power spectrum on sub-galactic scales the current, small sample of lensed radio-loud quasars may actually be. Another complicating factor makes the comparison between ΛCDM predictions with the observations still rather preliminary. In addition to substructures within the halo of the lensing galaxy, dwarf-sized perturbers along the line-of-sight to the lensed quasar may also affect the lens potential and give rise to flux ratio anomalies. And while lensing by dark matter clumps near the lens galaxy may be more effective, the cumulative effect of all intergalactic halos along the line-of-sight could be significant (Chen et al. 2003; Wambsganss et al. 2005; Metcalf 2005; Inoue & Takahashi 2012; Inoue 2016). With the next generation of wide-field optical surveys capable, e.g., of increasing the samples of multiplyimaged quasars by two orders of magnitude (Oguri & Marshall 2010), however, strong gravitational lensing appears poised to dramatically improve our understanding of dark matter substructure in galaxy halos. Flux ratio anomalies are only one means of detecting substructure. The most massive subhalos may also visibly perturb the deflection angle of lensed images, causing “astrometric anomalies” (Metcalf & Madau 2001; Chiba 2002; Chen et al. 2007). High-precision measurements of time delay perturbations between images in strong gravitational lens systems may complement flux ratio and astrometric anomalies as they depend on a different moment 2 (Msub ) of the subhalo mass function (Keeton & Moustakas 2009). An alternative technique, the “direct gravitational imaging” of individual clumps based on perturbations to the surface brightness of highly magnified Einstein rings, has been developed by Koopmans (2005). Detailed studies of individual Einstein ring systems have led Vegetti et al. (2012) to detect a (1.9 ± 0.1) × 108 M⊙ dark subhalo in an extended optical galaxy-galaxy lens system at z = 0.88. From a search of mass clumps in a sample of 11 lens galaxies from the SLACS, Vegetti et al. (2014) infer a mean projected substructure mass fraction Fsub = 0.0076+0.0208 −0.0052 (68% confidence level) around the Einstein radius of massive early-type host galaxies at

13

Figure 15. The differential projected mass function of subhalos at the Einstein radius. The errorbars indicate the 95% confidence limits on the projected number density of subhalos around the dusty galaxy SDP.81 (Hezaveh et al. 2016). For comparison, the shaded grey rectangles show the 90% confidence region in a decade of mass at z = 0.5 from the combined PonosQ and PonosV simulation sets. The red circle with the red error band shows the mean and 95% confidence limits in PonosQH in the mass bin 0.5 − 2 × 109 M⊙ .

hzlens i = 0.2. This fraction is again consistent with the expectations from ΛCDM simulations. Recently, Hezaveh et al. (2016, see also Inoue et al. 2016) have used ALMA observations of the strongly lensed dusty galaxy SDP.81 to find evidence for a Msub = 108.96±0.12 M⊙ subhalo near one of the images, and produce constraints on the projected abundance of substructure near the Einstein radius. Figure 15 shows the resulting constraints on the differential projected subhalo mass function, compared with predictions from the Ponos simulation sets. Again, the observations appear consistent with theoretical expectations. More importantly, the study by Hezaveh et al. (2016) shows that future ALMA data have the potential of constraining the abundance of dark matter subhalos down to Msub ∼ 2 × 107 M⊙ . Independently, Woldesenbet & Williams (2015) have reached similar conclusions (i.e. the necessity of & 108 M⊙ substructures projected near the Einstein radius) by using a model-free analysis of image positions in the relative angle space. Our hydro simulations do not have enough resolution to generate a realistic substructure population at masses Msub . 5 × 108 M⊙ , and higher resolution runs with dissipation are in the making in order to extend our analysis to smaller-mass subhalos. It is also necessary to perform more N -body simulations of early-type galaxy halos in order to obtain unbiased samples and reliable statistics of the subhalo population. We plan to perform a statistical exploration of flux anomalies induced by substructures near the Einstein radius on a larger halo sample. Since our results emphasize the relevance of dwarf-sized systems, such exploration may rely on a large suite of less challenging lower resolution simulations, as long as they can robustly resolve subhalos with masses & 108 M⊙ . This is instrumental to broadly asses the impact of substructures on the occurrence of flux anomalies in lensed

14

Fiacconi et al.

systems. We thank the anonymous Referee for useful comments that improved the quality of this work. We thank Lucio Mayer for useful discussions. We thank Oliver Hahn and the AGORA collaboration for help with the initial conditions of the simulations, and Yashar D. Hezaveh for providing the data plotted in Figure 15. Support for this work was provided to P.M. by the NSF through grant AST-1229745, by NASA through grant NNX12AF87G, and by the Pauli Center for Theoretical Studies. D.F. is supported by the Swiss National Science Foundation under grant 200021 140645. REFERENCES Agertz, O., Kravtsov, A. V., Leitner, S. N., & Gnedin, N. Y. 2013, ApJ, 770, 25 Agertz, O., Teyssier, R., & Moore, B. 2009, MNRAS, 397, L64 Amara, A., Metcalf, R. B., Cox, T. J., & Ostriker, J. P. 2006, MNRAS, 367, 1367 Auger, M. W., Treu, T., Bolton, A. S., et al. 2010, ApJ, 724, 511 Basu, S., & DasGupta, A. 1997, Theory of Probability & Its Applications, 41, 210 Behroozi, P. S., Wechsler, R. H., & Wu, H.-Y. 2013, ApJ, 762, 109 Blumenthal, G. R., Faber, S. M., Flores, R., & Primack, J. R. 1986, ApJ, 301, 27 Bode, P., Ostriker, J. P., & Turok, N. 2001, ApJ, 556, 93 Bolton, A. S., Treu, T., Koopmans, L. V. E., et al. 2008, ApJ, 684, 248 Bradaˇ c, M., Schneider, P., Lombardi, M., et al. 2004, A&A, 423, 797 Bressan, A., Marigo, P., Girardi, L., et al. 2012, MNRAS, 427, 127 Browne, I. W. A., Wilkinson, P. N., Jackson, N. J. F., et al. 2003, MNRAS, 341, 13 Bryan, G. L., & Norman, M. L. 1998, ApJ, 495, 80 Cappellari, M., Romanowsky, A. J., Brodie, J. P., et al. 2015, ApJL, 804, L21 Cembranos, J. A. R., Feng, J. L., Rajaraman, A., & Takayama, F. 2005, Physical Review Letters, 95, 181301 Chae, K.-H., Bernardi, M., & Kravtsov, A. V. 2014, MNRAS, 437, 3670 Chen, J., Koushiappas, S. M., & Zentner, A. R. 2011, ApJ, 741, 117 Chen, J., Kravtsov, A. V., & Keeton, C. R. 2003, ApJ, 592, 24 Chen, J., Rozo, E., Dalal, N., & Taylor, J. E. 2007, ApJ, 659, 52 Chiba, M. 2002, ApJ, 565, 17 Dalal, N., & Kochanek, C. S. 2002, ApJ, 572, 25 Dehnen, W. 2000, ApJL, 536, L39 —. 2002, Journal of Computational Physics, 179, 27 Diemand, J., Kuhlen, M., & Madau, P. 2007, ApJ, 657, 262 Diemand, J., Kuhlen, M., Madau, P., et al. 2008, Nature, 454, 735 El-Zant, A., Shlosman, I., & Hoffman, Y. 2001, ApJ, 560, 636 Fakhouri, O., Ma, C.-P., & Boylan-Kolchin, M. 2010, MNRAS, 406, 2267 Feng, J. L. 2010, ARA&A, 48, 495 Ferland, G. J., Korista, K. T., Verner, D. A., et al. 1998, PASP, 110, 761 Fiacconi, D., Feldmann, R., & Mayer, L. 2015, MNRAS, 446, 1957 Gill, S. P. D., Knebe, A., & Gibson, B. K. 2004, MNRAS, 351, 399 Gnedin, O. Y., Kravtsov, A. V., Klypin, A. A., & Nagai, D. 2004, ApJ, 616, 16 Governato, F., Brook, C., Mayer, L., et al. 2010, Nature, 463, 203 Guedes, J., Callegari, S., Madau, P., & Mayer, L. 2011, ApJ, 742, 76 Haardt, F., & Madau, P. 2012, ApJ, 746, 125 Hahn, O., & Abel, T. 2011, MNRAS, 415, 2101 Hezaveh, Y. D., Dalal, N., Marrone, D. P., et al. 2016, ArXiv e-prints, arXiv:1601.01388 Hinshaw, G., Larson, D., Komatsu, E., et al. 2013, ApJS, 208, 19 Hopkins, P. F., Kereˇs, D., O˜ norbe, J., et al. 2014, MNRAS, 445, 581 Hsueh, J.-W., Fassnacht, C. D., Vegetti, S., et al. 2016, ArXiv e-prints, arXiv:1601.01671 Hu, W., Barkana, R., & Gruzinov, A. 2000, Physical Review Letters, 85, 1158 Inoue, K. T. 2016, ArXiv e-prints, arXiv:1601.04414 Inoue, K. T., Minezaki, T., Matsushita, S., & Chiba, M. 2016, MNRAS, 457, 2936

Inoue, K. T., & Takahashi, R. 2012, MNRAS, 426, 2978 Keeton, C. R., Gaudi, B. S., & Petters, A. O. 2003, ApJ, 598, 138 Keeton, C. R., & Moustakas, L. A. 2009, ApJ, 699, 1720 Keller, B. W., Wadsley, J., Benincasa, S. M., & Couchman, H. M. P. 2014, MNRAS, 442, 3013 Kim, C.-G., & Ostriker, E. C. 2015, ApJ, 802, 99 Kim, J.-h., Abel, T., Agertz, O., et al. 2014, ApJS, 210, 14 Klypin, A., Kravtsov, A. V., Valenzuela, O., & Prada, F. 1999, ApJ, 522, 82 Klypin, A. A., Trujillo-Gomez, S., & Primack, J. 2011, ApJ, 740, 102 Knollmann, S. R., & Knebe, A. 2009, ApJS, 182, 608 Kochanek, C. S., Falco, E. E., Impey, C. D., et al. 2000, ApJ, 543, 131 Koopmans, L. V. E. 2005, MNRAS, 363, 1136 Kravtsov, A., Vikhlinin, A., & Meshscheryakov, A. 2014, ArXiv e-prints, arXiv:1401.7329 Kroupa, P. 2001, MNRAS, 322, 231 Li, R., Frenk, C. S., Cole, S., et al. 2015, ArXiv e-prints, arXiv:1512.06507 Macci` o, A. V., Moore, B., Stadel, J., & Diemand, J. 2006, MNRAS, 366, 1529 Madau, P., Shen, S., & Governato, F. 2014, ApJL, 789, L17 Mao, S., Jing, Y., Ostriker, J. P., & Weller, J. 2004, ApJL, 604, L5 Mao, S., & Schneider, P. 1998, MNRAS, 295, 587 McKean, J. P., Koopmans, L. V. E., Flack, C. E., et al. 2007, MNRAS, 378, 109 Metcalf, R. B. 2005, ApJ, 629, 673 Metcalf, R. B., & Amara, A. 2012, MNRAS, 419, 3414 Metcalf, R. B., & Madau, P. 2001, ApJ, 563, 9 Metcalf, R. B., & Zhao, H. 2002, ApJL, 567, L5 Moore, B., Ghigna, S., Governato, F., et al. 1999, ApJL, 524, L19 Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493 Oguri, M., & Marshall, P. J. 2010, MNRAS, 405, 2579 Pillepich, A., Kuhlen, M., Guedes, J., & Madau, P. 2014, ApJ, 784, 161 Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2015, ArXiv e-prints, arXiv:1502.01589 Pontzen, A., & Governato, F. 2012, MNRAS, 421, 3464 —. 2014, Nature, 506, 171 Power, C., Navarro, J. F., Jenkins, A., et al. 2003, MNRAS, 338, 14 Roˇskar, R., Fiacconi, D., Mayer, L., et al. 2015, MNRAS, 449, 494 Serra, P., Oosterloo, T., Morganti, R., et al. 2012, MNRAS, 422, 1835 Shen, S., Madau, P., Conroy, C., Governato, F., & Mayer, L. 2014, ApJ, 792, 99 Shen, S., Madau, P., Guedes, J., et al. 2013, ApJ, 765, 89 Shen, S., Wadsley, J., & Stinson, G. 2010, MNRAS, 407, 1581 Spergel, D. N., & Steinhardt, P. J. 2000, Physical Review Letters, 84, 3760 Springel, V., Wang, J., Vogelsberger, M., et al. 2008, MNRAS, 391, 1685 Stadel, J. 2013, PkdGRAV2: Parallel fast-multipole cosmological code, Astrophysics Source Code Library, ascl:1305.005 Stadel, J., Potter, D., Moore, B., et al. 2009, MNRAS, 398, L21 Stinson, G., Seth, A., Katz, N., et al. 2006, MNRAS, 373, 1074 Vegetti, S., Koopmans, L. V. E., Auger, M. W., Treu, T., & Bolton, A. S. 2014, MNRAS, 442, 2017 Vegetti, S., Lagattuta, D. J., McKean, J. P., et al. 2012, Nature, 481, 341 Viel, M., Becker, G. D., Bolton, J. S., & Haehnelt, M. G. 2013, Phys. Rev. D, 88, 043502 Wadsley, J. W., Stadel, J., & Quinn, T. 2004, Nature, 9, 137 Wambsganss, J., Bode, P., & Ostriker, J. P. 2005, ApJL, 635, L1 Weinberg, D. H., Bullock, J. S., Governato, F., Kuzio de Naray, R., & Peter, A. H. G. 2015, Proceedings of the National Academy of Science, 112, 12249 Wilman, D. J., & Erwin, P. 2012, ApJ, 746, 160 Woldesenbet, A. G., & Williams, L. L. R. 2015, MNRAS, 454, 862 Xu, D., Sluse, D., Gao, L., et al. 2015, MNRAS, 447, 3189 Xu, D. D., Mao, S., Wang, J., et al. 2009, MNRAS, 398, 1235 Zemp, M., Stadel, J., Moore, B., & Carollo, C. M. 2007, MNRAS, 376, 273 Zentner, A. R., Berlind, A. A., Bullock, J. S., Kravtsov, A. V., & Wechsler, R. H. 2005, ApJ, 624, 505 Zhu, Q., Marinacci, F., Maji, M., et al. 2015, ArXiv e-prints, arXiv:1506.05537