THE MUSIC OF CLASH: PREDICTIONS ON THE CONCENTRATION–MASS RELATION

The MUSIC-2 sample (Sembolini et al. 2013a, 2013b; Biffi et al. 2014) consists of a mass-limited sample of resimulated halos selected from the MultiDark cosmological simulation. This simulation is dark matter only and contains 2048³ (almost 9 billion) particles in a (1 h⁻¹ Gpc)³ cube. It was performed in 2010 using ART (Kravtsov et al. 1997) at the NASA Ames Research Center. All of the data of this simulation are accessible from the online MultiDark database.³⁵ The run was done using the best-fitting cosmological parameters to WMPA7+BAO+SNI (Ω_M = 0.27, Ω_b = 0.0469, Ω_Λ = 0.73, σ₈ = 0.82, n = 0.95, h = 0.7). This is the reference cosmological model used in the rest of the paper.

The halo sample was originally constructed by selecting all of the objects in the simulation box that are more massive than 10¹⁵ h⁻¹ M_☉ at redshift z = 0. In total, 282 objects were found above this mass limit. All of these massive clusters were resimulated both with and without radiative physics. The zooming technique described in Klypin et al. (2001) was used to produce the initial conditions for the resimulations. All particles within a sphere of 6 Mpc radius around the center of each selected object at z = 0 were found in a low-resolution version (256³ particles) of the MultiDark volume. This set of particles was then mapped back to the initial conditions to identify the Lagrangian region corresponding to a 6 h⁻¹ Mpc radius sphere centered at the cluster center of mass at z = 0. The initial conditions of the original simulations were generated in a finer mesh of size 4096³. By doing so, the mass resolution of the resimulated objects was improved by a factor of eight with respect to the original simulations. The parallel TREEPM+SPH GADGET code (Springel 2005) was used to run all of the resimulations. We stress that during the resimulation process, we make sure that all of the clusters that have been included in the MUSIC database are free from contamination by low-resolution particles that are outside the Lagrangian region of the resimulated area. If an object is formed close to the boundary of the high-resolution region, it might be very likely affected by the presence of particles with different spatial and mass resolution. In this case, we exclude this object from our analysis because it is not properly simulated. All of the MUSIC objects used in the analysis thus have their Lagrangian regions well inside the high-resolution regions defined by the 6 h⁻¹ Mpc spheres at z = 0.

The MUSIC-2 sample exists in two flavors. In a first set of resimulations, baryons were added to the dark matter distributions extracted from the parent cosmological box, and their physics was simulated via smoothed particle hydrodynamics (SPH) techniques, without including radiative processes. A second set of resimulations accounts for the effects of radiative cooling, UV photoionization, star formation, and supernova feedback, including the effects of strong winds from supernovae.

In this paper, we focus our analysis on the nonradiative version of these simulations. Our choice is based on the fact that radiative simulations without a proper description of energy feedback from active galactic nuclei (AGNs) generally produce unrealistically dense cores because of the well-known overcooling problem (see, e.g., Borgani & Kravtsov 2011). More recent simulations show that this problem is mitigated in simulations that simulate energy feedback from AGNs (Duffy et al. 2010; McCarthy et al. 2011; Planelles et al. 2014; Rasia et al. 2013a; Planelles et al. 2014). This physical ingredient is not yet included in the MUSIC-2 sample. Moreover, our intention is to correlate the profile measurements with the strong lensing efficiency of the simulated halos. Killedar et al. (2012), comparing simulations with different treatments of baryonic processes, find that the addition of gas in nonradiative simulations does not significantly change the strong lensing predictions. However, gas cooling and star formation together significantly increase the number of expected giant arcs and the Einstein radii by a nonrealistic amount, particularly for lower redshift clusters and lower source redshifts. Further inclusion of AGN feedback, however, reduces the predicted strong lensing efficiencies such that the lensing cross sections become closer to those obtained for simulations including only dark matter or nonradiative gas. The main requirements for this study are (1) a large number of highly resolved halos to accurately measure the profiles and determine the dependence of concentration on mass, and (2) the presence of gas in the simulations in order to allow their X-ray analysis (see Section 4). For these reasons, we choose to use the nonradiative version of the MUSIC-2 sample.

The mass resolution for these simulations corresponds to m_DM = 9.01 × 10⁸ h⁻¹ M_☉ and to m_SPH = 1.9 × 10⁸ h⁻¹ M_☉. The gravitational softening was set to 6 h⁻¹ kpc for the SPH and dark-matter particles in the high-resolution areas. Several low-mass clusters have been found close to the large ones and not overlapping with them. Thus, the total number of resimulated objects is considerably larger than originally identified in the parent cosmological box. In total, there are 535 clusters with M > 10¹⁴ h⁻¹ M_☉ at z = 0 and more than 2000 group-like objects with masses in the range 10¹³ h⁻¹ M_☉ < M_vir <10¹⁴ h⁻¹ M_☉. In this study, we use a subsample of these halos, as explained below.

We have stored snapshots for 15 different redshifts in the range 0 ⩽z ⩽ 9 for each resimulated object. The snapshots that overlap with the redshifts of the CLASH clusters are at z = 0.250, 0.333, 0.429, and 0.667.

The sample is complete above the mass thresholds given in Table 1. To extend our analysis toward smaller masses and to be able to constrain the concentration–mass relation over a wider mass range, we also analyze halos with masses below the completeness limits. In particular, we use all halos with mass M_vir > 2 × 10¹⁴ h⁻¹ M_☉. Therefore, we investigate a total of 1,419 halos, summing all halos at different redshifts.

2.2.1. Generalities

Navarro et al. (1996) argued that the density profiles of numerically simulated dark matter halos can be well fitted by an appropriate scaling of a "universal" function over a wide range of masses. The function suggested to fit these profiles was later dubbed the Navarro–Frenk–White density profile (NFW hereafter) and is given by

$$\begin{equation} \rho (r) = \frac{\rho _s}{(r/r_s)(1+r/r_s)^2}, \end{equation} \tag{ 1 }$$

where ρ_s and r_s are the characteristic density and the scale radius of the halo. The profile is characterized by a logarithmic slope that is shallower than isothermal for r ≪ r_s and steeper than isothermal for r ≫ r_s.

Subsequent numerical studies (see e.g., Navarro et al. 1997) confirmed that the NFW function is appropriate to describe the profiles of equilibrium halos, i.e., of systems that are close to being in virial equilibrium, and is now widely used to characterize the shape of cluster-sized halos both in observations and in simulations.

Along with the definition of the NFW density profile came that of the halo concentration, c_Δ = r_Δ/r_s, which is the ratio of the size of the halo, here defined as the radius enclosing a certain mean overdensity Δ above the critical density of the universe, ρ_crit(z). The most appropriate value to describe the size of an equilibrium halo is its virial radius, i.e., the radius within which the halo particles are gravitationally bound and settled into equilibrium orbits. In this case the virial overdensity, Δ_vir, is a function of cosmology and redshift (Bryan & Norman 1998; Nakamura & Suto 1997). To avoid this cosmological dependence, Navarro et al. (1996) adopted the round number of Δ = 200, which is commonly used in the literature independently of the assumed cosmological model. In this paper, we will also define the size of the halos as r₂₀₀, which is the radius enclosing a mean density $\overline{\rho }=200\rho _{{\rm crit}}(z)$. Diemer & Kravtsov (2014) recently showed that rescaling clusters to this radius returns a self-similar inner density profile.

Despite the fact that the profiles of equilibrium halos are well described by the NFW function, a large fraction of halos formed in a cosmological box are far from having reached virial equilibrium (Ludlow et al. 2012; Meneghetti & Rasia 2013). Balmès et al. (2014) discussed the dependence of this fraction on cosmology, finding that it is particularly sensitive to dark energy. The reason is simply understood: dark energy affects the formation and the growth of the cosmic structures. In the case of nonequilibrium halos, the NFW function gives a poorer description of the shape of the density profiles, and other functions involving a larger flexibility (i.e., additional free parameters) may yield a preferable result. One example is the generalized NFW profile (gNFW; Zhao 1996), which is given by

$$\begin{equation} \rho (r) = \frac{\rho _s}{(r/r_s)^\beta (1+r/r_s)^{3-\beta }}. \end{equation} \tag{ 2 }$$

Compared to the NFW model, this profile is characterized by an additional parameter, namely the logarithmic inner slope β,

$$\begin{equation} -\frac{{d}\ln \rho }{{d} \ln r} = \beta, \end{equation} \tag{ 3 }$$

which is radius-independent.

A strong debate exists in the literature about the inner slope of the density profile of simulated halos (see e.g., Moore et al. 1998; Newman et al. 2011). The advent of modern supercomputers allows us to push the mass and the spatial resolution of numerical simulations to unprecedented limits, and the new results indicate that there is a systematic deviation of the dark matter halo profiles from the form proposed by NFW (Merritt et al. 2006; Navarro et al. 2010). The function that best fits such profiles is the Einasto function (Einasto & Haud 1989; Retana-Montenegro et al. 2012),

$$\begin{equation} \rho (r)=\rho _{-2}\exp \left\lbrace -2n\left[\left(\frac{r}{r_{-2}}\right)^{1/n}-1\right]\right\rbrace, \end{equation} \tag{ 4 }$$

which is characterized by a running logarithmic slope,

$$\begin{equation} -\frac{{d}\ln \rho }{{d} \ln r}\propto r^{1/n}, \end{equation} \tag{ 5 }$$

parameterized in terms of the index n. The amplitude of the profile is set by the density ρ₋₂, which is the density at the radius r₋₂, i.e., at the radius where the logarithmic slope of the density profile is −2.

2.2.2. The Density Profiles of the MUSIC-2 Halos

To describe the structural properties of the MUSIC-2 halos, we perform an analysis of their three-dimensional density profiles based on the functional forms introduced in this section. This analysis is done by fitting the Equations (1), (2), and (4) to the azimuthally averaged density profiles of the simulated halos. The code used to perform this analysis is the same used in another CLASH paper by Merten et al. (2014).³⁶ As is common practice in the literature (e.g., Ludlow et al. 2013), we minimize the function

$$\begin{equation} R^{2}_{3D}=\frac{1}{N_{\rm dof}}\sum _{i}[\log _{10}{\rho _i}-\log _{10}\rho (r_i,\boldsymbol {p})]^2, \end{equation} \tag{ 6 }$$

where ρ_i is the density measured in the ith radial shell and $\boldsymbol {p}$ is the vector of parameters that are adjusted to derive the best-fitting function ρ(r). In the case of the NFW profile, $\boldsymbol {p}=[\rho _s,r_s]$, and in the cases of the gNFW or Einasto profiles, $\boldsymbol {p}=[\rho _s,r_s,\beta ]$ and $\boldsymbol {p}=[\rho _s,r_s,n]$, respectively. The variable N_dof is the number of degrees of freedom, i.e., the number of radii at which the profiles are evaluated minus the number of free parameters in the fit.

When analyzing these three-dimensional density profiles, we perform the fit over the radial range [$\tilde{r}_{{\rm min}}$, $\tilde{r}_{200}$], where $\tilde{r}_{{\rm min}}=0.02R_{{\rm vir}}$, and $\tilde{r}_{200}$ is the true r₂₀₀ of the halo. Of course, the choice of the radial range over which the fit is performed is important because substructures located within this range can affect the result of the fit (Meneghetti & Rasia 2013).

A similar analysis is performed on the two-dimensional profiles, i.e., on the azimuthally averaged surface-density profile, Σ_i, corresponding to an arbitrary line of sight to the halo. The details of this analysis are discussed in the paper by J. Vega et al. (in preparation). In this case, the fitting functions are the projections of the functions in Equations (1), (2), and (4):

$$\begin{equation} \Sigma (R)=2\int _0^{r_t}\rho (r=\sqrt{R^2+\xi ^2}){d}\xi, \end{equation} \tag{ 7 }$$

where ξ indicates the spatial coordinate along the line of sight, and R is the projected radius. In the formula above, r_t is a truncation radius, which is introduced to take into account that our halos are at the center of a cube with side length r_t = 6 h⁻¹ Mpc comoving. The figure-of-merit function to be minimized in this case is

$$\begin{equation} R^{2}_{2D}=\frac{1}{N_{\rm dof}}\sum _{i}[\log _{10}{\Sigma _i}-\log _{10}\Sigma (R_i,\boldsymbol {p})]^2\;\protect. \protect \protect \protect \end{equation} \tag{ 8 }$$

In order to be consistent with the analysis done on the CLASH clusters, we perform the two-dimensional fits over the radial range [20 h⁻¹ kpc, R_vir].

When projecting the cubes within which the halo particles are distributed, we expect that in particular for the smallest systems, there will be a two-halo contribution, which is not properly taken into account in the fitting procedure. To estimate if this may bias our conclusions, we repeat the fit using only the particles inside spheres with radius R_vir. The average concentrations do not change significantly, even at low masses, so we conclude that the two-halo contribution is a minor perturbation relative to the one-halo term for the radial scales we are probing (R_2D < R_vir).

For both the three- and the two-dimensional analyses, the best-fit parameters are used to compute the masses and the concentrations of the simulated halos. In the following, we identify the quantities estimated from these two analyses with the labels 3D and 2D, respectively. The best-fit masses are obviously obtained by integrating the best-fit density profiles,

$$\begin{equation} M=4\pi \int _{0}^{r_{200}}\rho (r,\boldsymbol {p}_{{\rm best}})r^2{d}r. \end{equation} \tag{ 9 }$$

The value of r₂₀₀ used here is derived by solving the equation

$$\begin{equation} \frac{\int _{0}^{r_{200}}\rho (r,\boldsymbol {p}_{{\rm best}})r^2{d}r}{r_{200}^3}=\frac{200}{3}\times \rho _{\rm crit}(z). \end{equation} \tag{ 10 }$$

Using its original definition (NFW), the concentration is the ratio between r₂₀₀ and the scale radius, r_s. For the NFW profile, the scale radius corresponds to the radius where

$$\begin{equation} -\frac{{d}\ln \rho }{{d} \ln r} = 2, \end{equation} \tag{ 11 }$$

that is, where the density profile has an isothermal slope. In the rest of the paper, we adopt the same definition also for the gNFW and Einasto profiles,

$$\begin{equation} c_{200} \equiv \frac{r_{200}}{r_{{-}2}}. \end{equation} \tag{ 12 }$$

Note that for the gNFW the following relation holds between r₋₂, the scale radius r_s, and the inner slope β:

$$\begin{equation} r_{{-}2} = (2-\beta)r_{s}. \end{equation} \tag{ 13 }$$

The lensing analysis of the MUSIC-2 halos is described in detail in Vega et al. (in preparation). For the purpose of this paper, we use their estimates of the Einstein radii over a large number of projections per cluster. We also use their convergence profiles, properly rescaled into surface-density profiles, and their mass and concentrations based on the fits of the surface-density profiles. The masses M_2D and the concentrations c_2D are equivalent to the values derived from a comprehensive lensing analysis of real observations. Hence, we compare M_2D and c_2D to Merten et al. (2014) and Umetsu et al. (2014).

For this work, we use our consolidated lensing simulation pipeline (see, e.g., Meneghetti et al. 2010a and references therein). Briefly, the following steps are involved.

• 1.  
  We project the particles belonging to each individual halo along the desired line of sight on the lens plane.
• 2.  
  Starting from the position of the virtual observer, we trace a bundle of light rays through a regular grid of 2048 × 2048 covering a region of 1.5 × 1.5 h⁻¹ Mpc around the halo center on the lens plane.
• 3.  
  Using our code RayShoot (Meneghetti et al. 2010b), we compute the deflection $\boldsymbol {\alpha }(\boldsymbol {x})$ at each light-ray position $\boldsymbol {x}$, accounting for the contributions from all particles on the lens plane.
• 4.  
  The resulting deflection field is used to derive several relevant lensing quantities. In particular, we use the spatial derivatives of $\boldsymbol {\alpha }(\boldsymbol {x})$ to construct the convergence, $\kappa (\boldsymbol {x})$, and the shear, $\boldsymbol {\gamma }=(\gamma _1,\gamma _2)$, maps. These are defined as

  $$\begin{eqnarray} \kappa (\boldsymbol {x}) & = & \frac{1}{2}\left(\frac{\partial \alpha _1}{\partial x_1}+\frac{\partial \alpha _2}{\partial x_2}\right), \end{eqnarray} \tag{ 14 }$$

  $$\begin{eqnarray} \gamma _1(\boldsymbol {x}) & = & \frac{1}{2}\left(\frac{\partial \alpha _1}{\partial x_1}-\frac{\partial \alpha _2}{\partial x_2}\right), \end{eqnarray} \tag{ 15 }$$

  $$\begin{eqnarray} \gamma _2(\boldsymbol {x}) & = & \frac{\partial \alpha _1}{\partial x_2} = \frac{\partial \alpha _2}{\partial x_1}. \end{eqnarray} \tag{ 16 }$$
• 5.  
  The lens critical lines are defined as the curves along which the determinant of the lensing Jacobian is zero (e.g., Schneider et al. 1992):

  $$\begin{equation} \det A = (1-\kappa -|\gamma |)(1-\kappa +|\gamma |) = 0. \end{equation} \tag{ 17 }$$

  In particular, the tangential critical line is defined by the condition (1 − κ − |γ|) = 0, whereas the radial critical line corresponds to the line along which (1 − κ + |γ|) = 0. In the following sections, we will often use the term Einstein radius to refer to the size of the tangential critical line. As discussed in Meneghetti et al. (2013), there are several possible definitions for the Einstein radius. In this paper, we adopt the effective Einstein radius definition (see also Redlich et al. 2012),

  $$\begin{equation} \theta _E \equiv \frac{1}{d_{\rm L}}\sqrt{\frac{S}{\pi }}, \end{equation} \tag{ 18 }$$

  where S is the area enclosed by the tangential critical line and d_L is the angular diameter distance to the lens plane.

All of the lensing quantities are computed for a source redshift z_s = 2.

In order to increase the statistics and to take into account possible projection effects, J. Vega et al. (in preparation) study each halo under a large number of lines of sight. More precisely, they investigate 100 lines of sight for the halos above the mass completeness limits and 30 projections for those below the completeness limit. This implies that for each halo, we have a catalog containing at least 30 measurements of the Einstein radius, projected mass, and projected concentration.

We build a mock X-ray catalog by producing for each simulated cluster three Chandra event files corresponding to orthogonal projections aligned with the Cartesian axes of the simulation. Because of excessive computational demand, we cannot investigate all of the lines of sight considered in J. Vega et al. (in preparation). The images are created by the X-ray MAp Simulator (X-MAS; Gardini et al. 2004), in which we utilize the ancillary response function and redistribution matrix function proper of the ACIS-S3 detector (for a complementary X-ray analysis of the MUSIC-2 sample, we refer the reader to Biffi et al. 2014). The field of view (FOV) covers (16 arcmin)². For the cosmology and redshifts analyzed, the FOV size is equivalent to the following physical scales: 5.43 h⁻¹ Mpc at z = 0.250, 6.57 h⁻¹ Mpc at z = 0.333, 7.71 h⁻¹ Mpc at z = 0.429, and 9.57 h⁻¹ Mpc at z = 0.667. The spectral emission is generated by adopting the MEKAL model in which we fix the redshift to the simulation's value and the metallicity to a constant value equal to 0.3 times the solar metallicity as tabulated by Anders & Grevesse (1989). Finally, the contribution of the galactic absorption is introduced through a WABS model with N_H = 5 × 10²⁰ cm⁻² (see e.g., Lemze et al. 2009). The exposure time is set to 100 ks, allowing a fair comparison with observations.

One of the goals of this paper is to identify halos in the MUSIC-2 sample that closely resemble the X-ray properties of the clusters in the CLASH X-ray-selected sample. Because this sample was selected to have a high degree of regularity in the X-ray morphology, we try to find equivalents in the simulations that mimic these X-ray characteristics.

We use five parameters to measure the X-ray morphology in the soft-energy band ([0.5–2] keV) images of our halos. These morphological parameters are evaluated within a physical radius R_max = 500 kpc following the same procedure adopted in the X-ray analysis of the CLASH clusters. The results of this analysis will be presented in detail in a forthcoming paper by Donahue et al. (2014). The five parameters are:

• 1.  
  the centroid-shift, w, which assesses how much the centroid of the X-ray surface brightness moves when the aperture radius used to compute it decreases from R_max to smaller values. It is defined as

  $$\begin{equation} w=\frac{1}{R_{\max }} \times \sqrt{\frac{\Sigma (\Delta _i-\langle \Delta \rangle)^2}{N-1}}\protect, \protect \end{equation} \tag{ 19 }$$

  where N is the total number of apertures considered and Δ_i is the separation of the centroids computed within R_max and within the ith aperture;
• 2.  
  the ellipticity, e = 1 − b/a, where the axial ratio is equal to the ratio of the square root of the eigenvalues obtained by diagonalizing the inertia tensor of the X-ray surface brightness evaluated within R_max (Buote & Canizares 1992);
• 3.  
  the X-ray-brightness concentration, which is the ratio between the integral of the surface brightness S within two apertures with radii 100 kpc and R_max,

  $$\begin{equation} c_{X}=\frac{S(r<100 \ {\rm kpc})}{S(r<R_{{\rm max}})} \end{equation} \tag{ 20 }$$

  Cassano et al. (2010);
• 4.  
  and 5. the third- and fourth-order power ratios, P₃ and P₄. These are the third- and fourth-order multipoles of the surface-brightness distribution within an aperture of radius R_ap = R_max. The generic m-order power ratio (m > 0) is defined as P_m/P₀ with

  $$\begin{equation} P_m=\frac{1}{2m^2R^{2m}_{\rm ap}}\left(a_m^2 +b^2_m\right) \ \ \ \ {\rm and} \ \ \ \ P_0=a_0 \ln (R_{\rm ap})\protect, \protect \end{equation} \tag{ 21 }$$

  where a₀ is the total intensity within the aperture radius Buote & Tsai (1996). The generic moments a_m and b_m are expressed in polar coordinates, R' and ϕ', and given by

  $$\begin{equation} a_m(r)=\int _{R^{\prime } \le R_{\rm ap}} S(x^{\prime }) R^{\prime } \cos (m \phi ^{\prime }) d^2 x^{\prime }, \end{equation} \tag{ 22 }$$

  and

  $$\begin{equation} b_m(r)=\int _{R^{\prime } \le R_{\rm ap}} S(x^{\prime }) R^{\prime } \sin (m \phi ^{\prime }) d^2 x^{\prime }. \end{equation} \tag{ 23 }$$

For a review of X-ray morphological parameters, we refer to Rasia et al. (2013b).

The five morphological parameters introduced above are combined to define a global degree of X-ray regularity. Such a quantity is measured with respect to the mean of the simulated sample. Note that with reference to the X-ray appearance, we use the term "regular" to indicate halos with unperturbed surface brightness distributions (Rasia et al. 2012). Very often, these halos are called "relaxed. " We do not use this term to differentiate from the classification discussed in Section 5.2. Regular halos have small centroid shift, ellipticity, and power ratios. In addition, they have large surface-brightness concentrations. Thus, we define the regularity parameter

$$\begin{eqnarray} M& = &\left(\frac{\log _{10}(w)-\langle \log _{10}(w)\rangle }{\sigma _{\log _{10}w}}\right)+\left(\frac{e-\langle e \rangle }{\sigma _{e}}\right) \nonumber \\ & &+ \left(\frac{\log _{10}(1/c_{X})-\langle \log _{10}(1/c_{X})\rangle }{\sigma _{\log _{10}1/c_{X}}}\right) \nonumber \\ & &+ \left(\frac{\log _{10}(P_3)-\langle \log _{10}(P_3)\rangle }{\sigma _{\log _{10}P_3}}\right) \nonumber \\ & &+ \left(\frac{\log _{10}(P_4)-\langle \log _{10}(P_4)\rangle }{\sigma _{\log _{10}P_4}}\right) \end{eqnarray} \tag{ 24 }$$

similarly to the M parameter derived in Rasia et al. (2013b). In the formula above, each morphological parameter, p_i, is compared to its mean over the simulated halos, 〈p_i〉, and rescaled by the standard deviation of its distribution, $\sigma _{p_i}$.

By plugging the parameters p_{CLASH, i} measured on the X-ray images of the CLASH clusters into Equation (24), we use the M parameter to quantify the regularity of the CLASH clusters with respect to the simulations. The M parameters of the CLASH X-ray-selected clusters are listed in Table 3. To construct a sample of CLASH-like clusters, we select the simulated halos having a regularity parameter similar to the observed clusters.

For the purpose of matching simulated halos to each individual CLASH cluster, we define the parameter C_X, which is defined as the distance, in parameter space, between each CLASH cluster and the simulated halos:

$$\begin{eqnarray} C_X & = & \sum _{i=1,5} \left(\frac{p_i-p_{{\rm CLASH},i}}{\sigma _{p_i}}\right)^2, \end{eqnarray} \tag{ 25 }$$

where p_i = [log₁₀(w), e, −log₁₀(c_X), log₁₀P₃, log₁₀P₄] are the morphological parameters discussed above and $\sigma _{p_i}$ their standard deviations. As a result, the sample constructed via the M parameter has an X-ray regularity similar to the CLASH sample. When we match halos using C_X, we identify only the simulated halos closest to each individual CLASH cluster in the morphological parameter space.

While our choice to use the nonradiative version of the MUSIC-2 halos is motivated by the need to avoid biases caused by overcooling, it is well known that hydrodynamical simulations like those employed here poorly describe several X-ray properties of real clusters (Borgani & Kravtsov 2011; Kravtsov & Borgani 2012). For this reason, we do not use gas temperatures or X-ray luminosities to match the CLASH clusters in our simulations. Our comparison is based solely on the X-ray morphology.

To evaluate how the morphological parameters used in this work are influenced by the treatment of the gas, we use the hydrodynamical simulations described in Fabjan et al. (2010) and in Bonafede et al. (2011; see also Killedar et al. 2012; Planelles et al. 2014). These simulations, performed in the framework of a cosmological setting similar to that of the MUSIC-2 simulations, exist both in nonradiative and radiative versions. Contrary to the MUSIC-2 simulations, the effects of AGN feedback are also included in the radiative case. The sample is significantly smaller, though. Seventy of these halos were recently processed with the X-MAS simulator, both in the nonradiative and radiative versions. We use this analysis to quantify the impact of radiative processes on the morphological parameters.

The distributions derived from the two simulated sets are consistent for all morphological parameters computed within 500 kpc, with the exception of the light concentration, which is lower in the radiative simulation because part of the central gas is turned into star and contributes less to the X-ray central emission. Applying the selection method based on the parameter C_X on the halos in these two data sets for a few CLASH clusters, we obtained identical matches. Therefore, we can assume that our X-ray selection method can safely be used on the nonradiative simulations.

Notice that the two samples are characterized by similar concentration distributions, as shown in Rasia et al. (2013a). In particular, that paper found the following results: (1) the c–M relations have similar slopes independent of the physics, and (2) the normalization of the c–M relation in radiative simulations with no AGN feedback is ∼20% higher than that of the nonradiative simulations; the c–M relation from simulations including AGN feedback has a normalization similar to that of the nonradiative simulations.

To illustrate how our selection based on the X-ray morphology performs, we discuss the case of A383 (Allen et al. 2008), which is the first cluster observed in the framework of the CLASH program. A383 is a galaxy cluster at redshift z = 0.189 (see, e.g., Zitrin et al. 2011). In the X-ray, it exhibits a very regular morphology, with nearly circular surface brightness contours (ellipticity ∼0.04; Postman et al. 2012). An X-ray image taken from the Archive of Chandra Cluster Entropy Profile Tables (ACCEPT) is shown in the small inset at the center of Figure 1. The image subtends ∼345.

Zoom In Zoom Out Reset image size

The four largest panels of Figure 1 show a sequence of simulated Chandra observations of MUSIC-2 halos corresponding to increasing values of C_X, which are annotated on the images. The top left panel shows the X-ray morphology of the halo that best matches A383 (C_X = 0.2). The X-ray morphology is indeed very similar to that of the observed cluster. As C_X increases, the differences between the simulated and the true X-ray morphologies become more significant. On the basis of this and other visual inspections, we verified that C_X ∼ 0.4 represents a good limit to select the halos "similar" to the true cluster.

Our selection successfully identifies simulated halos that also closely resemble more perturbed clusters. For example, this is the case for MACSJ 1149 (Ebeling et al. 2007), which is one of the CLASH clusters identified as high-magnification clusters, i.e., not included in the X-ray-selected sample. A comparison between the true X-ray morphology and that of a simulated halo with C_X = 0.18 is shown in Figure 2, where we show the true Chandra image of the cluster in the smaller inset on the right.

Zoom In Zoom Out Reset image size

Clearly, the degree of asymmetry and of elongation of the surface-brightness distribution in the simulated observation matches very closely that of MACSJ 1149.

In this section we differentiate between relaxed and unrelaxed halos on the basis of a few criteria that are commonly used in the literature. Following the most restrictive approach proposed by Neto et al. (2007), we classify as strictly relaxed (or super-relaxed, as we dub them later in the paper) those objects satisfying the following properties.

• 1.  
  Their center of mass displacement, defined as the offset between the center of mass (determined using all of the particles within the virial radius) and the minimum of the potential, in units of the virial radius, is $s=(\boldsymbol {r}_{\rm cm}-\boldsymbol {r}_{\phi })/r_{\rm vir}<0.07$.
• 2.  
  Their virial ratio is η = 2T/|U| < 1.35, where T is the kinetic energy and U is the gravitational energy, computed using the particles within the virial radius.
• 3.  
  Their substructure mass fraction, computed as the mass in resolved substructures within the virial radius, is f_sub < 0.1.

Applying these selection criteria to the MUSIC-2 halos results in a fraction of relaxed halos of about 14.9% at redshift z = 0.25. The fraction is reduced to 11.7% at redshift 0.333, and it further drops to 10.4% and 8.9% at redshifts 0.429 and 0.667, respectively.

Other authors use less restrictive or alternative criteria to identify the relaxed systems (e.g., Skibba & Macciò 2011; Skibba et al. 2011). For example, Bhattacharya et al. (2013) only use the center-of-mass displacement. In their paper, they report that the addition of the two other conditions on η and f_sub does not significantly affect the selection. On the contrary, we find that using only the center-of-mass displacement we end up with a significantly higher fraction of halos being classified as relaxed. This fraction amounts to ∼60% at z = 0.250 and decreases to ∼51% at z = 0.667. Such fractions are compatible with those quoted by Bhattacharya et al. (2013) (see also Biffi et al. 2014). Sembolini et al. (2013a) recently used the center-of-mass displacement in combination with the virial ratio to identify relaxed systems in simulations. They report that the relation between η and s becomes flat for s ≲ 0.1, thus indicating that η does not impact severely on the selection of relaxed systems. For our sample, the combination of s and η yields a fraction of relaxed halos corresponding to 47% at z = 0.250, which decreases to 29% at z = 0.667.

In the following sections, we will study the properties of the MUSIC-2 halos, dividing them into three subsamples. First, we will consider all halos, regardless of their relaxation state. Second, we will set the limit defined above on the center-of-mass displacement to construct the subsample of relaxed halos. Third, we will further downsize the sample by using all three criteria described above to identify the super-relaxed halos.

As explained in Section 2.2, we fit the density profiles of the MUSIC-2 halos using the functions in Equations (1), (2), and (4). In Figure 3, we show the results of the fitting procedure. We quantify the goodness of fit by means of the residuals given in Equations (6) and (8).

Zoom In Zoom Out Reset image size

The upper left panel shows the distributions of the fit residuals for the entire MUSIC-2 sample. When all halos are considered, regardless of their relaxation state, the NFW profile is the worst-fitting model, i.e., the one with the largest residuals (see also Meneghetti & Rasia 2013). This is not surprising given that the NFW model has one free parameter less than the gNFW or the Einasto profiles. However, this result highlights the difficulty of fitting all profiles with a universal law. Because the gNFW and the Einasto functions generally provide better fits to the profiles, we may use the statistical distributions of their residuals to identify the halos deviating significantly from the NFW form. As can be seen from Figure 3, such distributions are nearly log-normal, which suggests that halos having too-large NFW residuals compared to the Einasto and gNFW models may be identified via their deviation δ = ln R_{3D, NFW} − 〈ln R_{3D, x}〉, where 〈ln R_{3D, x}〉 is the mean value of ln R_3D for either the Einasto or the gNFW model. Using this criterion, we find that about 40% of the halos have NFW fits resulting in too-large residuals compared to what is typically found by fitting with more flexible profiles.

This fraction drops to ∼19% and ∼6% if only relaxed and super-relaxed halos are considered. The distributions of the fit residuals for the super-relaxed subsample are shown in the upper right panel of Figure 3. For these halos, the NFW model is only a slightly worse fit compared to the gNFW and Einasto models.

In Figure 4 we see that the profiles that most deviate from the NFW form have inner slopes β (resulting from the gNFW fits) that significantly differ from unity: their profiles are steeper or shallower than the NFW model. There is a slight indication for preferring a steep over a shallow slope (see also Figure 5). Indeed, the mean value of the inner slope β, measured for the whole sample, is 〈β〉 = 1.03 ± 0.31, where the error is the rms in the sample. We also find that the goodness of the gNFW fit is not correlated with the inner slope β, i.e., shallow or steep inner slopes are not systematically the result of a bad gNFW fit.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

When fitting the surface-density profiles, we find again that the NFW model is generally the worst-fitting function among the three models employed in this work. This is shown in the bottom left panel of Figure 3. Again, we find that restricting the analysis to the relaxed halos reduces the differences between the residual distributions of the NFW and gNFW or Einasto fits. However, from the results shown in the bottom right panel of Figure 3, it appears that a fraction of halos that are well fitted by NFW models in 3D are not NFW-like in projection. This result must be caused by the halo triaxiality and by the effects of substructures and additional matter along the line of sight. The work of J. Vega et al. (in preparation), from which the 2D analysis shown here is taken, investigates the effects of triaxiality on the shape of the surface-density profiles of the CLASH clusters. We refer the reader to that paper for more details. We note that the halo surface-density profiles were derived by using all of the particles in a cylinder centered on the halo and with depth 6 h⁻¹ Mpc.

The distributions of the inner slopes obtained from the gNFW fits of the surface-density profiles are shown in Figure 5 (thick histograms). We find that a large number of halos have rather flat profiles in 2D. The mean value of β is 〈β〉 = 0.89 ± 0.47. About 33% (15%) of the halo projections are fitted with β ⩽ 0.8 (⩽0.5). The red histograms show the distributions of the Einasto indexes 1/n. The indexes obtained from the fit of the density profiles are slightly smaller than what are obtained from the fit of the surface-density profiles. The smaller the 1/n, the steeper the inner profile. The mean values are 〈1/n〉 = 0.21 ± 0.07 and 〈1/n〉 = 0.24 ± 0.09 for the 3D and 2D distributions, respectively. Such Einasto slopes appear to be in excellent agreement with the recent results of Dutton & Macciò (2014).

To summarize, the halos in the MUSIC-2 sample span a wide range of structural parameters. As expected, the density profiles can differ significantly from the NFW form, and their shape can be better described with more flexible functions, such as the Einasto or gNFW models. When projecting the mass distributions, the scatter in the profile parameters and the deviation from the NFW model become even larger.

Having determined the level of diversity among density and surface-density profiles of the MUSIC-2 halos, we consider now how precisely the halo masses are derived from the profile fits. We consider both the cases of 2D and 3D masses, the former being the masses derived by deprojecting the best-fit models of the surface-density profiles under the assumption of spherical symmetry, and the latter those derived from the fits of the density profiles. Note that when measuring the 2D masses, we are not simulating any lensing analysis at this stage. In particular, we are not considering additional sources of systematics that may depend on the particular method to derive the mass from the weak and the strong lensing signals. Other works have shown that different methods of analysis may introduce systematic errors that are due, for example, to the presence of substructures inside and outside the clusters (Meneghetti et al. 2010b; Becker & Kravtsov 2011; Rasia et al. 2012) and to the Bright Central Galaxy (Giocoli et al. 2013). Nevertheless, this exercise gives us important information on the intrinsic limits of the mass measurements based on the analyses of azimuthally averaged density or surface-density profiles.

We begin with the 3D masses. The distributions of the ratios between such masses and the true halo masses are shown in the left panel of Figure 6. The results are shown for the three fitting functions employed in this work (black, red, and orange histograms). We find that the masses recovered from the azimuthal fits of the density profiles are generally in good agreement with the true masses. The best agreement is obtained with the Einasto and gNFW profiles, with a slight preference for the first. These fits provide ratios around unity with rms = 0.06 and 0.05, respectively. The masses estimated through the NFW fits are also in good agreement with the true masses. In this case, the median (mean) ratio is 0.98 (0.97) and the distribution is twice as broad as in the two previous cases. The purple histogram is constructed by choosing, for each cluster, the mass estimate derived from the fitting function leading to the smallest residuals. In other words, we choose the most reliable mass estimate among those obtained with the three fitting functions. In most cases, the best model is the Einasto profile. Thus, the purple and the orange histograms are nearly coincident.

Zoom In Zoom Out Reset image size

The histograms shown here refer to the whole halo sample, regardless of the relaxation state. As shown in the previous section, the density profiles of the relaxed halos are generally equally well fitted by NFW, gNFW, or Einasto models. Indeed, restricting the analysis to these halos, we find smaller rms for all three kinds of fit (≲ 0.03), with mean and median ratios very close to unity. Despite the fact that the fraction of relaxed halos varies with redshift, we find that the mean mass ratios and their scatter remain constant as a function of redshift.

Even when fitting the surface-density profiles, the mass estimates (M_2D) deviate only slightly from the true masses. The 2D masses appear to be underestimated by ∼5% on average, with the NFW and gNFW fits being slightly more biased than the Einasto fits. However, the scatter is much larger (∼13%–14%) than for the 3D masses. The larger scatter is expected, given that the masses are derived under the assumption of spherical symmetry. Halos are generally triaxial, and projection effects can easily cause the mass to be over- or underestimated by a significant amount, depending on the halo orientation (see, e.g., Meneghetti et al. 2010b). As reported by Giocoli et al. (2012a), the halo prolateness may also cause a systematic underestimation of the mass derived from the 2D analysis. Assuming the triaxial model of Jing & Suto (2002), they estimate this bias to be of order ∼10%.

As in the left panel, the purple histogram in the right panel of Figure 6 shows the distribution of the ratios between the best 2D mass estimate and the true mass. Again, the distribution is close to that obtained by fitting with the Einasto profile.

On the basis of this result, we conclude that we should expect a modest negative bias of ∼5% on the mass estimates obtained by fitting the surface-density (or the convergence) profiles of galaxy clusters. This is due to the prolate shape of the halos, which are more frequently elongated on the plane of the sky than along the line of sight. The choice of the NFW or gNFW models to fit the halos tend to slightly increase the bias, and the opposite occurs with the Einasto profile.

If we repeat this analysis on the samples of relaxed and super-relaxed halos, we find that the mass bias tends to become smaller. In fact, the 2D masses deviate from the true masses by only ∼1%–2% in these cases. If the bias originates from halo triaxiality, this suggests that the most relaxed systems must generally be more spherical. In Figure 7, we show the distribution of the axis ratios b/a and c/a of all of the MUSIC-2 halos (color map). Here a, b, and c are the semi-axes of the inertial ellipsoid fitting the mass distribution of the halos with a > b > c. This fit is done using all particles within the virial radius. It is clear from this plot that the relaxed (yellow dashed contours) and super-relaxed systems (white contours) generally have higher values of both b/a and c/a. Thus, their shape is closer to spherical than that of nonrelaxed halos, in agreement with Lemze et al. (2012).

Zoom In Zoom Out Reset image size

The concentration–mass (c − M − z) relation is derived by means of nonlinear least-squares fitting using a Levenberg–Marquardt algorithm. The fitting function we employ is

$$\begin{equation} c(M,z) = A\left(\frac{1.34}{1+z}\right)^B\left(\frac{M}{8\times 10^{14}\,h^{-1}\ M_{\odot }}\right)^C, \end{equation} \tag{ 26 }$$

which was also used by Duffy et al. (2008) and De Boni et al. (2013), although using a different pivot mass and redshift. We perform this analysis for the three fitting models considered and report the corresponding best-fit parameters and errors in Table 2. The results are reported for the full sample as well as for the subsamples of relaxed and super-relaxed halos. We use Equation (26) to fit the c–M–z relations derived from the analyses of the density profiles.

5.4.1. Comparison between Fitting Models

In the following, we consider the concentrations obtained from the NFW fit of the density profiles as a reference when making comparisons with the concentrations derived from the gNFW and Einasto fits. The yellow and green histograms in the upper panel of Figure 8 show the the distributions of the ratios c_{3D, gNFW}/c_{3D, NFW} and c_{3D, Einasto}/c_{3D, NFW} obtained from our analysis. In both cases, we find that the distributions peak at values around ∼0.9–0.95, with the Einasto concentrations being generally smaller than the NFW ones. This result is in agreement with the recent findings of Dutton & Macciò (2014), who also find that the Einasto concentrations are ∼10–15% smaller than the NFW concentrations on the mass scale of the MUSIC-2 halos. The halos with the smallest concentrations are of course the unrelaxed systems, for which we already pointed out that the NFW model is generally a bad fit. An example of such profiles is shown in the bottom panel of Figure 8. In this case, the best-fit NFW concentration is c_{3D, NFW} = 2.5, and the gNFW and Einasto concentrations are c_{3D, gNFW} = 10⁻² and c_{3D, Einasto} = 0.1, respectively. Considering only the relaxed or the super-relaxed halos, the ratios between fitted and true concentrations are much closer to unity. For example, the mean ratios of c_{3D, gNFW}/c_{3D, NFW} and c_{3D, Einasto}/c_{3D, NFW} for the super-relaxed systems are 1.0 and 0.99, respectively. We want to stress that the concentration of the Einasto profile being smaller than the NFW does not imply necessarily that the halos are less concentrated. For the Einasto profile, the mass inside the scale radius also depends on the 1/n parameter. A halo with the same mass ratio between two radii as given by the NFW model can be fitted with a smaller concentration and a larger n.

Zoom In Zoom Out Reset image size

In Figure 9, we show the c–M relations obtained from fitting the density profiles of the MUSIC-2 with the NFW, gNFW, and Einasto models (upper, middle, and bottom panels, respectively). The results are displayed for the halos at the lowest redshift investigated in this work (z = 0.250). Each circle corresponds to a halo, and the solid, dashed, and dotted lines indicate the best-fit c–M–z relations for the full, relaxed, and super-relaxed samples. At fixed mass, the distribution of NFW halo concentrations is reasonably well fitted by a log-normal distribution and have a standard deviation σ_c ∼ 0.25, compatible with the findings of several previous works (see e.g., Dolag et al. 2004). The concentrations derived from the gNFW and Einasto fits are characterized by a larger scatter. In all cases we find that the dependence of the concentration on mass is very shallow. For the NFW profile, c∝M^{−0.057 ± 0.017} for the full sample. Instead, for the gNFW and the Einasto profiles, the logarithmic slope of the c–M relation is slightly positive. For the relaxed and super-relaxed halos, all of the c–M relations have logarithmic slopes that are negative or consistent with zero. As expected, we find that the more relaxed the halos are, the higher their concentrations (Zhao et al. 2009; Giocoli et al. 2012b). This result holds regardless of the fitting function. At the lowest masses, the relative change in typical concentrations between the full and the relaxed (or super-relaxed) samples is larger for the gNFW and Einasto fits. In fact, we find that a larger fraction of small-mass unrelaxed halos are fitted with lower concentrations using these two fitting models than with the NFW profile. These halos are responsible for the positive logarithmic slope of the c–M relation when fitting with the gNFW or Einasto profiles.

Zoom In Zoom Out Reset image size

As can be seen from the B parameters listed in Table 2, the normalization of the 3D c–M relation has an almost negligible redshift dependence for the full sample. For example, in the case of the NFW profile, c∝(1 + z)^{−0.29 ± 0.08}. For the gNFW and Einasto profiles, the redshift evolution is even shallower. We notice, however, that the dependence of the concentration on redshift appears to be stronger for the most relaxed systems. In particular, for the super-relaxed halos, we find B ∼ 0.52, regardless of the fitting function.

5.4.2. The NFW Concentration–Mass Relation

There are several parameterizations of the c–M relation in the literature, mostly derived from fitting simulated halos using NFW profiles. In the upper panel of Figure 10, we show the NFW c–M–z relation derived from the 3D analysis for the whole sample of MUSIC-2 halos (solid lines). We use different colors to show how the relation evolves with redshift. We find a rather shallow dependence of the concentrations on mass and redshift. Over the mass range [4–12 × 10¹⁴ h⁻¹ M_☉], the concentrations vary by less than 10%, decreasing as a function of mass as M^{−0.058 ± 0.017}. The amplitude of the c−M relation scales with redshifts as (1 + z)^{−0.29 ± 0.08}. Other authors find that the c–M relation of massive halos is rather flat. For example, Zhao et al. (2009), studying an ensemble of numerical simulations in the context of various cosmological models, find that the concentration is strongly correlated with the age of the universe when the halo progenitor on the mass accretion history first reaches 4% of its current mass. According to this correlation, they find that the concentration is nearly constant for halos with mass M ≳ 10¹⁴ h⁻¹ M_☉. They also predict a very shallow redshift evolution of the c–M relation. In a recent work, De Boni et al. (2013) also find concentrations that scale with mass and redshift, similar to our results. Their concentrations scale with mass and redshift as M^−0.07 and (1 + z)^−0.26, respectively.

Zoom In Zoom Out Reset image size

The normalization of our c–M–z relation is higher than that found by some other authors, like De Boni et al. (2013; dot-dashed lines in the upper panel of Figure 10) or Duffy et al. (2008). In these cases, the differences can be explained in terms of different cosmological settings. For example, De Boni et al. (2013) analyze halos evolved in the framework of a WMAP3 cosmological model and adopt a rather small normalization of the matter power spectrum, σ₈ = 0.72. If we consider other analyses in the literature in the context of WMAP7 normalized cosmologies, the agreement is much better. For example, the c–M relation that best fits our data at low redshift is in rather good agreement with the results of Bhattacharya et al. (2013) for nonrelaxed halos. For comparison, their c−M relation is overplotted in the upper panel of Figure 10 (dashed lines). At z = 0.250, the concentrations we measure at a given mass are only ≲ 6% higher than found by Bhattacharya et al. (2013). However, their c–M relation has a stronger redshift evolution. Between z = 0.250 and z = 0.667, their concentrations at a fixed mass decrease by ∼17%, while ours vary only by ∼10%.

Potentially important differences between this work and Bhattacharya et al. (2013) are (1) our simulations include baryons, while the halos studied by Bhattacharya et al. (2013) are made only of dark matter; (2) our analysis focuses on a limited mass range, and the volume we sample is smaller than in the simulations employed by Bhattacharya et al. (2013); (3) the mass resolution of our simulations is roughly two orders of magnitude better; (4) Bhattacharya et al. (2013) fit their halos over a different radial range, [0.1–1r_vir] versus [0.02–1r₂₀₀]; and, finally, (5) Bhattacharya et al. (2013) fit the mass profiles instead of the density profiles, as we do. Given that our simulations are nonradiative, it is unlikely that the differences between the c–M relations arise from baryonic effects. De Boni et al. (2013) show that concentrations are higher by 5%–15% in radiative simulations compared to dark matter-only simulations. This result, however, was obtained using hydrodynamical simulations that are known to suffer from the overcooling problem. It has been shown by other authors that halos in adiabatic simulations develop density profiles pretty similar to those of pure dark matter halos (Killedar et al. 2012). The different mass range, volume, and resolution of the simulations may have a larger impact on the results. Because our halos are sampled with a larger number of particles, the profiles are better resolved. Thus, the measurements of the individual concentrations should be more robust and allow us to resolve smaller radial scales. On the other hand, because their volume is bigger, Bhattacharya et al. (2013) have a larger number of massive halos to constrain the c–M relation at the cluster scales. In contrast, Bhattacharya et al. (2013) fit halos over three orders of magnitude in mass. It may be possible that the strong redshift evolution of their c–M relation is driven by the smallest halos. Overall, it is likely that the higher normalization of our c–M is largely due to the better resolution of the MUSIC-2 sample compared to the simulation sets used in Bhattacharya et al. (2013).

The bottom panel in Figure 10 shows another comparison between our best-fit NFW c–M–z relation and the results of Bhattacharya et al. (2013). The solid and the dotted lines show our relations for relaxed and the super-relaxed samples, respectively. The most striking difference from Bhattacharya et al. (2013) (dashed lines) is that we find a much stronger dependence of concentration on the halo dynamical state. While the normalization of our c–M–z relation increases by ∼10% between the full and relaxed samples, Bhattacharya et al. (2013) find that concentrations of relaxed halos increase only by ∼3%.

5.4.3. The c–M Relation and Temperature Selection

5.4.4. The Concentration–Mass Relation in 2D

The 2D concentration–mass relation of the MUSIC-2 halos will be discussed in detail in an upcoming paper (J. Vega et al., in preparation). We briefly summarize some properties of this c–M relation that are relevant for the following discussion. Projection effects do affect concentrations, which are generally found to be smaller than in 3D. This effect of triaxiality, also discussed in Giocoli et al. (2012a), is illustrated in Figure 11, where we show the distribution of the MUSIC-2 halos in the (c_2D − c_3D)/c_3D versus (M_2D − M_3D)/M_3D plane. The 2D histograms show that regardless of the fitting model, the masses and concentrations derived from fitting the surface-density profiles tend to be smaller than measured from fitting the density profiles. The trend is in qualitative agreement with the findings of Giocoli et al. (2012a), although the amplitude of both the concentration and mass biases found here is smaller. The white contours overlaid on the 2D histograms show the intensity levels corresponding to 10%, 50%, and 90% of the peaks of the distributions. The red contours indicate the same intensity levels for the subsample of super-relaxed halos. As explained in Section 5.3, the bias is reduced for the relaxed halos because these systems are typically more spherical than the unrelaxed halos.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The best-fit parameters of the 2D c–M–z relation are listed in Table 2. The relations for halos at z = 0.250 are given by the solid lines in Figure 13. Interestingly, the c–M relation is very flat and characterized by an inverted slope compared to the c–M relation in 3D. This suggests that the 2D concentrations underestimate the 3D ones more significantly at the lowest than at the highest masses. One possible explanation is that the halo triaxiality is somehow biased below the completeness limits listed in Table 1. However, we checked that the c–M relation obtained only from halos above the completeness limits does not differ significantly from what we obtain using the extended sample. However, the constraints on its slope are obviously weakened. In addition, we notice that Giocoli et al. (2012a) also find indications for a 2D concentration bias that decreases as a function of mass. This will be discussed in J. Vega et al. (in preparation).

Zoom In Zoom Out Reset image size

As explained above, the CLASH cluster sample is composed of 25 galaxy clusters, of which only five were selected on the basis of their SL strength. The remaining 20 clusters are not SL selected, and they were chosen on the basis of their X-ray morphology. We will discuss this selection method in the next section. Nevertheless, SL features (multiple images and arcs) have been securely detected in all CLASH clusters except RXJ 1532.8+3021. The analyses of these SL features have allowed the creation of detailed lens models and the measurement of their Einstein radii (Zitrin et al. 2014). The Einstein radii for sources at redshift z_s = 2 are within the range 5–55 arcsec.

We construct the c–M–z relation of SL galaxy clusters by selecting those projections where we measure an Einstein radius compatible with those measured in the CLASH sample (θ_E > 5''). As explained, the Einstein radius is defined as in Equation (18).

In Figure 13, we show the concentration–mass relations at z = 0.250 derived from SL halos in the MUSIC-2 sample. The relations are displayed for the NFW, gNFW, and Einasto models (upper, middle, and bottom panels). The corresponding parameters are listed in Table 2. The dotted, dashed, and dot–dashed lines indicate the relations obtained for the full, the relaxed, and the super-relaxed samples, respectively. For comparison, we show also the c–M relation for the full sample, including also the nonstrong lenses, as discussed in Section 5.4.4. By requiring that the halos are strong lenses in their projections, we remove a large fraction of halos with low concentrations, obtaining relations characterized by a larger normalization. In particular, an increasingly larger number of halos of small mass are unable to produce an appreciable SL signal. By removing them from the initial catalog, we restore the negative logarithmic slope of the c–M relation. Because of this selection, the concentration scales with mass as c∝M^{−0.214 ± 0.018}. This result is in very good agreement with the theoretical predictions of Giocoli et al. (2013) and Oguri et al. (2012), who estimated the lensing bias of the c–M via semianalytic calculations employing triaxial halos. For the gNFW and Einasto models, the concentration–mass relations are slightly flatter.

Even for the SL halos, the normalization of the c–M–z relation depends on the relaxation state. The most relaxed systems have the largest concentrations. The differences between the c–M–z relations of relaxed and unrelaxed halos are smaller than found earlier for the whole sample including nonstrong lenses, though. This is because in the SL-selected sample the fraction of relaxed and super-relaxed halos is pretty high. At z = 0.250, about 75% of the SL projections belong to relaxed halos. The fraction of super-relaxed halos in this sample is ∼27%. At z = 0.667 the fractions of relaxed and super-relaxed halos are ∼60% and ∼13%, respectively.

Finally, we find that the redshift evolution of the c−M relation of SL halos is stronger than for non-SL halos. The values for the B parameters listed in Table 2 are in the range [0.48–0.64] for the three fitting models.

We discuss now the impact of the X-ray selection on the concentration–mass relation. In particular, we discuss the expectations for halos selected so as to resemble the X-ray morphologies of the clusters in the CLASH X-ray-selected sample. The results shown here are based on the analysis of three projections per halo, and the halos considered are those with 3D mass above the completeness limits given in Table 1. The restriction of the analysis to this smaller sample of simulated halos was dictated by the large computational time required to produce the X-ray-simulated observations. The mass range covered by these simulations is however representative of the mass range of the CLASH clusters (see both Merten et al. 2014 and Umetsu et al. 2014).

As explained in Section 4.1, the X-ray morphology is measured by means of five morphological parameters. They can be combined to quantify the degree of regularity of the halos, as shown in Equation (24). The regularity parameters of the CLASH clusters, as measured in their X-ray images, are listed in Table 3. In Figure 14, we show their distribution (red histogram), and we compare it to the distribution of the regularity parameters in the MUSIC-2 sample (black histogram). The histograms have been normalized to have the same peak value. As it emerges from these distributions, the CLASH clusters have quite typical regularity parameters in the simulations. With the exception of MACSJ 1206.2-0847, the clusters in the CLASH X-ray-selected sample have negative M parameters, indicating that they are more regular than the mean of the simulations. Even in the case of MACSJ 1206.2-0847, other works based on different analyses find that this system is not likely to be perturbed by significant substructures (Lemze et al. 2013; Biviano et al. 2013). This is expected given that these clusters were selected on the basis of their X-ray regularity. On the other hand, the comparison shows that their regularity is not extreme, in the sense that there are several simulated halos with regularity parameters exceeding the values for the CLASH clusters. In fact, the simulated sample has a tail of low M values extending well beyond those of the CLASH clusters.

Zoom In Zoom Out Reset image size

In the upper panel of Figure 15, we show that the concentration inferred from the analysis of the 2D mass distributions is correlated with the regularity parameter M. The red, green, and blue circles refer to unrelaxed, relaxed, and super-relaxed halos. The correlation was evaluated by measuring the linear Pearson correlation coefficient P between the log₁₀c_{200, 2D} and the M values. It is stronger for the super-relaxed halos, for which we measure P = −0.67. For the relaxed and the full samples, we obtain P = −0.46 and P = −0.39, respectively. The best linear fit between the two parameters is

$$\begin{eqnarray} &&\log _{10}c_{200,2D}=(0.598\pm 0.009)-(0.019\pm 0.002)\times M.\nonumber \\ \end{eqnarray} \tag{ 27 }$$

Zoom In Zoom Out Reset image size

If we refer to the average of all halos in the simulations (M = 0 by construction), for negative values of M we expect a positive concentration bias. Because the median value of the M parameters of the CLASH clusters is M_CLASH = −3.44, on the basis of Equation (27), we can give an estimate of the expected concentration bias for the CLASH X-ray-selected sample, which is

$$\begin{equation} \frac{c_{200,{\rm CLASH}}}{c_{200,2D}(M=0)}=1.11\pm 0.03. \end{equation} \tag{ 28 }$$

An interesting question is whether this concentration excess compared to the full sample arises from the selection of purely relaxed halos. The answer is already contained in the upper panel of Figure 15: a selection based on the M regularity parameter does not lead to the construction of a purely relaxed sample. Indeed, the left side of the diagram contains several red circles, indicating that unrelaxed halos can have M < 0. The composition of samples selected by means of the M parameter is shown in the bottom panel of Figure 15. The curves show the fractions of relaxed (R), nonrelaxed (NR), and super-relaxed (SR) halos in the samples with regularity parameter smaller than M. As indicated by the dotted and the solid black lines (R and NR halos), the fraction of relaxed and unrelaxed halos is nearly constant as a function of M. Thus, there is no strong correlation between X-ray regularity and halo relaxation. In particular, we find that only ∼70% of the halos among those with M < 0 are relaxed.³⁷ The remainder ∼30% of the halos are unrelaxed. As stated, this composition is very similar to that of the full sample.

In contrast, as indicated by the dashed line, the fraction of super-relaxed halos decreases as a function of M. Thus, in samples of clusters selected to have regular X-ray morphologies, we expect to have a larger fraction of super-relaxed halos. Because these typically have larger concentrations, we expect that the average concentration in an M-selected sample is higher than in the full sample.

The red solid line shows that the fraction of strong lensing (SL) halos in M-selected samples also decreases as a function of M. This trend reflects the correlation between concentration and regularity parameter. Because the halos are more concentrated, they more easily act as strong lenses. However, we notice that a correlation exists also between the concentrations and the M parameters of the unrelaxed halos, although this is weaker than for the relaxed and super-relaxed halos. For unrelaxed halos, the linear Pearson coefficient is P = −0.22, indicating also that among these halos, those with a small M parameter tend to have larger concentrations. In part, the classification of unrelaxed halos as regular is due to the different radial scales over which the relaxation and the regularity are evaluated. Whereas the former is measured using all particles inside the virial radius, the second is meant to quantify the morphology of the cluster cores, within 500 kpc. A fraction of halos with regular X-ray morphologies have significant substructures outside 500 kpc, which implies that they are classified as unrelaxed. These substructures explain the low concentrations of those unrelaxed halos that have small M. However, for ∼42% of the unrelaxed halos with M < 0, we do not find evidence for substructures outside the region where we carry over the X-ray morphological classification. These halos generally have 2D concentrations higher than the average of the sample, indicating that the selection based on X-ray morphology may lead to the inclusion of unrelaxed objects that are elongated along the line of sight. Such a sample would then be affected by a small orientation bias.

We use the M parameters to create a sample of X-ray-selected halos. These halos are drawn from the full MUSIC-2 sample so as to reproduce the distribution of the M parameters found for the CLASH clusters. In doing so, we take into account the masses and redshifts of the CLASH clusters. The masses are taken from Merten et al. (2014). A halo is selected if it has a suitable M parameter and if the mass inferred from the 2D analysis is within 3σ of the mass measured by CLASH. To account for the redshift distribution, given that the halos available for this analysis are simulated only at four redshifts, we create a match between each CLASH cluster and the nearest simulated redshift. The matches are listed in Table 3.

As explained earlier in the paper, the X-ray analysis is limited to three orthogonal lines of sight per halo. Given that many more projections are available in the 2D analysis of the MUSIC-2 halos, we can improve our statistical power by increasing the number of projections used. To do so, we identify the projections whose lines of sight are within $20\deg$ of those selected in the X-ray analysis.

Using the concentrations and masses inferred from the 2D analysis of the X-ray-selected projections, we fit the c–M–z relation for our X-ray-selected sample. The relation is shown in Figure 16 for all of the fitting models employed in this study. The best-fit parameters are listed in Table 2. Overall, the c–M–z relation for X-ray-selected halos is in good agreement with the SL c–M–z relation for a sample composed of both relaxed and unrelaxed halos. This is not surprising given that X-ray-selected halos are frequently efficient strong lenses, with only ∼8% of them not having an extended critical line for sources at z = 2. About 70% of the selected projections belong to relaxed halos. About 18% of them correspond to halos classified as super-relaxed. For the NFW model, we find that the concentrations scale with mass as c∝M^{−0.16 ± 0.11}, resulting in average concentrations that are intermediate between those predicted in 3D for relaxed and super-relaxed halos in the mass range 2 × 10¹⁴ ≲ M₂₀₀ ≲ 10¹⁵ h⁻¹ M_☉.

Zoom In Zoom Out Reset image size

Some differences between the fitting models are found with regards to the redshift evolution of the concentration–mass relation. For the NFW model, the c–M–z relation is evolving strongly. The redshift dependence is shallower in the case of the gNFW model, and it is consistent with zero evolution for the Einasto profile.