© ESO, 2014

1. Introduction

Do the strongest observed gravitational lenses exceed the theoretical expectations of the standard cosmological model? This question has long been debated in the literature (see Bartelmann 2010; Meneghetti et al. 2013, for reviews), and was the central theme of this series of papers.

In our first work (Redlich et al. 2012, hereafter Paper I), we developed a new semi-analytic method for populating a fictitious observer’s past null cone with merging galaxy clusters. We studied the impact of cluster mergers on the optical depth for giant gravitational arcs and the distribution of Einstein radii, finding that mergers are particularly relevant for the strongest gravitational lenses in the Universe. Furthermore, we showed that strong-lensing cross sections for giant gravitational arcs are tightly correlated with Einstein radii. This correlation allowed us to infer optical depths from distributions of Einstein radii, and also justified why we focused on the statistics of Einstein radii in the subsequent works.

In Waizmann et al. (2012, hereafter Paper II), we developed a framework for applying the theory of extreme value statistics to the distribution of the single largest Einstein radius in a given cosmological volume. We found that the extreme value distribution is particularly sensitive to the precise choice of the halo mass function, lens triaxiality, the inner slope of the halo density profile, and the mass-concentration relation. Most importantly, and in contrast to many previous studies, we demonstrated that the largest known Einstein radius – which was observed in the galaxy cluster MACS J0717.5+3745 (Zitrin et al. 2009, 2011; Limousin et al. 2012; Medezinski et al. 2013) – is consistent with the theoretical expectations of the standard cosmological model.

Given the insight that the occurrence probability of the single largest Einstein radius is strongly influenced by several theoretical uncertainties, we argued that it would be desirable to formulate ΛCDM exclusion criteria that are based on a number of observations instead of a single event. Therefore, we extended our theoretical framework to higher order statistics in the third work (Waizmann et al. 2014, hereafter Paper III). More precisely, we studied exclusion criteria based on the individual and joint order distributions of the n largest Einstein radii. Most importantly, again, our computations revealed that the large Einstein radii of the twelve high-redshift (z > 0.5) MACS clusters analysed by Zitrin et al. (2011) do not exceed the expectations of the ΛCDM model.

However, so far, our analysis is still incomplete in the following sense: while we explicitly emphasised the importance of cluster mergers for the statistics of the strongest gravitational lenses in Paper I, we did not include mergers in our computations for Papers II and III. This decision was motivated by several factors. Firstly, particularly in Paper II, we investigated the properties of triaxial gravitational lenses (such as orientation and profile shape) that decisively affect the extreme value distribution of the largest Einstein radius. For these studies, it was important to isolate individual effects, and therefore it was neither desirable nor required to include cluster mergers. Secondly, we aimed to derive conservative exclusion criteria. We were able to show that the strongest observed gravitational lenses do not exceed the theoretical expectations of the standard cosmological model even if cluster mergers are neglected. Because mergers additionally boost the strong lensing efficiency, including these events should additionally consolidate this conclusion. Thirdly, including clusters mergers is computationally very expensive. A typical realisation of a mock universe requires Monte Carlo simulations of ~10⁶ merger trees. In addition, lensing computations with multiple haloes in the field of view are substantially more expensive than those with single haloes. To estimate the order statistics of the largest Einstein radii with an acceptable precision, ~10³ mock universes have to be sampled. This whole procedure is computationally very demanding, even though our semi-analytic method was specifically tailored to be fast.

In this sense, this Paper completes the present series. Based on previously developed techniques, we compute the impact of cluster mergers on the order statistics of the largest expected Einstein radii in the Universe. This is not only interesting in itself, but also because it differs from the analysis performed in Paper I, where we investigated the impact of cluster mergers on averaged quantities such as the optical depth for giant gravitational arcs. Moreover, we compare our theoretical results to another, independent set of observational data at redshifts z < 0.5 (Zitrin et al. 2012), finding that the observed gravitational lenses in this cosmological volume might be somewhat stronger than theoretically expected.

The plan for this Paper is as follows: in Sect. 2, we briefly summarise the most important theoretical concepts that were introduced in our previous works. Section 3 describes the merger tree algorithm used for this work, and also contains some information on speeding up the computations. The impact of cluster mergers on the order statistics of the largest Einstein radii is analysed in Sect. 4. Thereafter, in Sect. 5, we compare our theoretical predictions with current observational data. In Sect. 6, we summarise our main results and finally conclude.

Throughout this work, we adopt the best-fitting cosmological parameters (Ω_Λ0,Ω_m0,Ω_b0,h,σ₈) = (0.685, 0.273, 0.047, 0.673, 0.829), which were derived from the Planck 2013 data (Planck Collaboration XVI 2014).

2. Theoretical basis

For convenience, we briefly summarise the most important theoretical concepts this work is based upon. We refer to our previous publications for more thorough introductions to triaxial gravitational lenses (Paper I) as well as extreme value (Paper II) and order statistics (Paper III).

2.1. Triaxial gravitational lenses

As extensively discussed in Paper I, the triaxiality of galaxy clusters is very important for their strong-lensing properties. We therefore adopt the triaxial, generalised NFW profile proposed by Jing & Suto (2002) to model the dark matter density profile of galaxy clusters: $\begin{array}{ccc}& & \mathit{\rho }\mathrm{\left(}\mathit{R}\mathrm{\right)}\mathrm{=}\frac{{\mathit{\delta }}_{\mathrm{ce}}{\mathit{\rho }}_{\mathrm{crit}}\mathrm{\left(}\mathit{z}\mathrm{\right)}}{\mathrm{(}\mathit{R}\mathit{/}{\mathit{R}}_{\mathrm{0}}{\mathrm{)}}^{\mathit{\alpha }}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{R}\mathit{/}{\mathit{R}}_{\mathrm{0}}{\mathrm{)}}^{\mathrm{3}\mathrm{-}\mathit{\alpha }}}\mathit{,}\\ & & {\mathit{R}}^{\mathrm{2}}\mathrm{\equiv }\frac{{{\mathit{x}}^{\mathrm{\prime }}}^{\mathrm{2}}}{\mathrm{(}\mathit{a}\mathit{/}\mathit{c}{\mathrm{)}}^{\mathrm{2}}}\mathrm{+}\frac{{{\mathit{y}}^{\mathrm{\prime }}}^{\mathrm{2}}}{\mathrm{(}\mathit{b}\mathit{/}\mathit{c}{\mathrm{)}}^{\mathrm{2}}}\mathrm{+}{{\mathit{z}}^{\mathrm{\prime }}}^{\mathrm{2}} \mathrm{(}\mathit{a}\mathrm{\le }\mathit{b}\mathrm{\le }{\mathit{c}}^{\mathrm{)}}\mathit{.}\end{array}$Here, z denotes the halo redshift, R₀ is a scale radius, δ_ce is a characteristic density, and ρ_crit(z) denotes the critical density of the universe; see Paper I for more detail on these quantities, the sampling of random density profiles, and the computation of the lensing signal.

In this work, we set the inner slope α of the density profile to the standard NFW value α = 1. As detailed in Sect. 5.2 of Paper II, we force the scaled axis ratios a_sc (minor axis) to lie within the 99% confidence interval of the probability density function (PDF) for a_sc derived by Jing & Suto (2002). By doing so, we avoid unrealistic density profiles with extremely low axis ratios and too low concentrations.

2.2. Definition of the Einstein radius

Axially symmetric gravitational lenses exhibit circular tangential critical curves, whose size can naturally be quantified by means of their Einstein radius. In contrast, the tangential critical curves of merging triaxial gravitational lenses are irregular and cannot be described by simple geometric figures (cf. Paper I). Nevertheless, it has been proven useful to generalise the concept of an Einstein radius to arbitrarily shaped critical curves.

In this work, we use a geometrically motivated definition. Let A denote the area enclosed by the tangential critical curve of an arbitrary gravitational lens. The effective Einstein radius θ_eff of this lens is then defined by ${\mathit{\theta }}_{\mathrm{eff}}\mathrm{\equiv }\sqrt{\frac{\mathit{A}}{\mathit{\pi }}}\mathit{,}$(3)so that a circle with radius θ_eff has the same area A.

2.3. Extreme value statistics

The theory of extreme value statistics (EVS) describes the stochastic behaviour of extremes and provides a rigorous framework for determining the likelihood of rare events (Gumbel 1958). In particular, the theory defines the mathematical recipe for quantitatively answering questions of the following kind: Given a cosmological model, how likely is it that the largest Einstein radius in a certain cosmological volume is as large as X″?

The method developed in Paper II is based on the so-called Gnedenko approach, which models the distribution of block maxima (Fisher & Tippett 1928; Gnedenko 1943). In simple words, it can be summarised as follows: suppose that we intend to calculate the extreme value distribution of the largest Einstein radius in a certain cosmological volume. To this end, we populate the volume with a mock catalogue of triaxial gravitational lenses and compute their Einstein radii. We note the largest Einstein radius ${\mathit{\theta }}_{\mathrm{eff}}^{\mathrm{1}}$ of the current realisation, and repeat this random process n times. After this procedure, we have a set of n largest Einstein radii, $\mathit{M}\mathrm{=}{}^{\mathrm{(}}\mathit{\theta }_{\mathrm{eff}}^{\mathrm{1}}\mathit{,}\mathit{...}\mathit{,}{{\mathit{\theta }}_{\mathrm{eff}}^{\mathit{n}}}^{\mathrm{)}}$. As was shown in Paper II, for large n, the cumulative distribution function (CDF) of these maxima converges to the general extreme value (GEV) distribution ${\mathit{G}}_{\mathit{\gamma ,}\hspace{0.17em}\mathit{\beta ,}\hspace{0.17em}\mathit{\alpha }}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{=}\{\begin{array}{c}\\ \exp \left\{\mathrm{-}{\left[\mathrm{1}\mathrm{+}\mathit{\gamma }\left(\frac{\mathit{x}\hspace{0.17em}\mathrm{-}\hspace{0.17em}\mathit{\alpha }}{\mathit{\beta }}\right)\right]}^{\mathrm{-}\mathrm{1}\mathit{/}\mathit{\gamma }}\right\}\mathit{,}& \mathrm{for} \mathit{\gamma }\ne \mathrm{0}\mathit{,}\\ \exp \left\{{\mathrm{e}}^{\mathrm{-}\left(\frac{\mathit{x}\hspace{0.17em}\mathrm{-}\hspace{0.17em}\mathit{\alpha }}{\mathit{\beta }}\right)}\right\}\mathit{,}& \mathrm{for} \mathit{\gamma }\mathrm{=}\mathrm{0}\mathit{,}\end{array}$(4)with the location, scale and shape parameters α, β, and γ.

Fitting the GEV to the empiric CDF has several advantages. Amongst others, the fit allows us to derive analytic relations for probabilities and exclusion criteria, and also enables us to smoothly extrapolate to regions that are only sparsely sampled by the empiric data.

2.4. Order statistics

In Paper III, we argued that the GEV of the single largest Einstein radius is sensitive to several theoretical uncertainties and that it might therefore be advisable to formulate more robust exclusion criteria based on n observations instead of a single one. To do this, we extended our framework to the order statistics of the n largest Einstein radii in a given cosmological volume.

For a random sample of n Einstein radii ${\mathit{\theta }}_{\mathrm{eff}}^{\mathrm{1}}\mathit{,}{\mathit{\theta }}_{\mathrm{eff}}^{\mathrm{2}}\mathit{,}\mathit{...}\mathit{,}{\mathit{\theta }}_{\mathrm{eff}}^{\mathit{n}}$, the order statistic is given by the ordered set ${\mathit{\theta }}_{\mathrm{eff}}^{\mathrm{\left(}\mathrm{1}\mathrm{\right)}}\mathrm{\le }{\mathit{\theta }}_{\mathrm{eff}}^{\mathrm{\left(}\mathrm{2}\mathrm{\right)}}\mathrm{\le }\mathit{...}\mathrm{\le }{\mathit{\theta }}_{\mathrm{eff}}^{\mathrm{\left(}\mathit{n}\mathrm{\right)}}$, so that ${\mathit{\theta }}_{\mathrm{eff}}^{\mathrm{\left(}\mathrm{1}\mathrm{\right)}}$ and ${\mathit{\theta }}_{\mathrm{eff}}^{\mathrm{\left(}\mathit{n}\mathrm{\right)}}$ denote, respectively, the smallest (minimum) and the largest (maximum) Einstein radius in the sample. The CDF of the ith order is ${\mathit{F}}_{\mathrm{\left(}\mathit{i}\mathrm{\right)}}\mathrm{\left(}{\mathit{\theta }}_{\mathrm{eff}}\mathrm{\right)}\mathrm{=}\sum _{\mathit{k}\mathrm{=}\mathit{i}}^{\mathit{n}}\left(\begin{array}{c}\mathit{n}\\ \mathit{k}\end{array}\right){\left[\mathit{F}\mathrm{\left(}{\mathit{\theta }}_{\mathrm{eff}}\mathrm{\right)}\right]}^{\mathit{k}}{\left[\mathrm{1}\mathrm{-}\mathit{F}\mathrm{\left(}{\mathit{\theta }}_{\mathrm{eff}}\mathrm{\right)}\right]}^{\mathit{n}\hspace{0.17em}\mathrm{-}\hspace{0.17em}\mathit{k}}\mathit{,}$(5)where $\mathit{F}{}^{\mathrm{(}}\mathit{\theta }_{\mathrm{eff}}{}^{\mathrm{)}}$ denotes the CDF the sample of Einstein radii was drawn from. As discussed above (and, more extensively, in Paper III), the distribution function of the maxima ${\mathit{F}}_{\mathrm{\left(}\mathit{n}\mathrm{\right)}}{}^{\mathrm{(}}\mathit{\theta }_{\mathrm{eff}}{}^{\mathrm{)}}$ converges to the GEV for large sample sizes ${\mathit{F}}_{\mathrm{\left(}\mathit{n}\mathrm{\right)}}\mathrm{\left(}{\mathit{\theta }}_{\mathrm{eff}}\mathrm{\right)}\mathrm{=}{\left[\mathit{F}\mathrm{\left(}{\mathit{\theta }}_{\mathrm{eff}}\mathrm{\right)}\right]}^{\mathit{n}}\mathrm{\to }{\mathit{G}}_{\mathit{\gamma ,}\hspace{0.17em}\mathit{\beta ,}\hspace{0.17em}\mathit{\alpha }}\mathrm{\left(}{\mathit{\theta }}_{\mathrm{eff}}\mathrm{\right)}\mathit{,}$(6)which, in a certain domain, allows us to relate the underlying distribution function $\mathit{F}{}^{\mathrm{(}}\mathit{\theta }_{\mathrm{eff}}{}^{\mathrm{)}}$ to the GEV and, using Eq. (5), derive analytic expressions for the CDFs of the second largest, third largest, etc. Einstein radius.

3. Algorithm for including cluster mergers

3.1. Extended Press-Schechter merger tree algorithm

In Paper I, we developed a semi-analytic method that allows us to project merging galaxy clusters onto a fictitious observer’s past null cone and calculate their strong-lensing properties. Our method can be briefly summarised as follows: First, we specify the redshift range [z_min,z_max] and the sky fraction f_sky that we are about to analyse, and compute the size of the comoving volume that corresponds to these parameters. Next, we envisage this cosmological volume at redshift z = 0 and populate it with a representative halo population by sampling from a halo mass function. We assign each halo a (properly weighted) random observation redshift z_obs${}^{\mathrm{(}}\mathit{z}_{\min }\mathit{<}{\mathit{z}}_{\mathrm{obs}}\mathit{<}{\mathit{z}}_{\max }{\mathrm{)}}^{}$ and use a merger tree algorithm to evolve the halo backwards in time up to z_obs. A merger tree simulation is a Monte-Carlo process that starts with the initial halo mass M₀ at z = 0, takes a discrete time step Δz backwards in time, splits the halo into smaller progenitors M₀ → M₁ + M₂ (where the masses M₁ and M₂ are randomly drawn from a probabilistic merger rate), and iteratively applies the same procedure to all arising progenitors until z_obs is reached. The major splitting events that occur during these simulations represent cluster mergers and are highly important for the strong-lensing statistics. Using a simplistic model for the kinematics during cluster mergers, we determine the trajectories of all progenitors that arise during the merger tree simulation. In particular, we compute the angular positions of the progenitors on the observer’s sky at redshift z_obs. Applying this procedure to all initial haloes, we effectively populate the past null cone with spatially correlated haloes, some of which are actively merging.

In Paper I, we emphasised that we employed two completely independent algorithms to compute the strong-lensing statistics: (1) a standard Monte Carlo approach to sample isolated haloes from a given mass function; and (2) our new method that incorporates cluster mergers. We did so in preparation of the following works, but also to independently cross-check the new merger algorithm. For self-consistency, the two methods had to be based on the original mass function of Press & Schechter (1974).

On the other hand, in Paper II we found that the precise choice of the halo mass function has a significant impact on the GEV distribution of the largest Einstein radius. Because it is based on the overly idealistic theory of spherical collapse, the mass function of Press & Schechter (1974) is known to underestimate the abundance of high-mass haloes, which in turn leads to lower expectation values for the largest Einstein radius in a given cosmological volume. Conversely, the mass function of Sheth & Tormen (1999) over-predicts the number of high-mass haloes, which shifts the GEV distribution to too high values. The Tinker et al. (2008) (Tinker) mass function is a compromise between these two extremes and is broadly accepted as a more accurate representation of the mass function determined from N-body simulations. Therefore, all calculations of Paper III were based on the Tinker mass function.

In this work, we aim to calculate the order statistics of the largest Einstein radii as realistically as possible, so that our results can be compared with real observations. It would therefore be desirable to replace the original merger algorithm by an alternative algorithm that generates halo catalogues consistent with the Tinker mass function.

Jiang & van den Bosch (2014) recently compared several different merger tree codes, finding that the algorithm proposed by Parkinson et al. (2008) (PCH) reproduces the merger rates measured in N-body simulations most accurately. The basic idea of the PCH algorithm is quite simple: although merger algorithms based on the original extended Press-Schechter formalism are slightly inaccurate, they are certainly not completely unacceptable. Instead, they exhibit statistical properties (such as trends with mass and redshift) that agree well with those of merger trees constructed from N-body simulations. Parkinson and co-workers thus proposed to take the Press-Schechter merger rate as a starting point, and slightly perturb this quantity using an empirically motivated function. The probability for drawing a progenitor of mass M₁ from a parent halo of mass M₀ at redshift z is then modified according to $\frac{\mathrm{d}\mathit{N}}{\mathrm{d}{\mathit{M}}_{\mathrm{1}}}\mathrm{\to }\mathit{G}\mathrm{\left(}{\mathit{M}}_{\mathrm{0}}\mathit{,}{\mathit{M}}_{\mathrm{1}}\mathrm{\right)}\hspace{0.17em}\frac{\mathrm{d}\mathit{N}}{\mathrm{d}{\mathit{M}}_{\mathrm{1}}}\mathit{,}$(7)with the empirical function $\mathit{G}\mathrm{\left(}{\mathit{M}}_{\mathrm{0}}\mathit{,}{\mathit{M}}_{\mathrm{1}}\mathrm{\right)}\mathrm{=}{\mathit{G}}_{\mathrm{0}}\hspace{0.17em}{\left[\frac{\mathit{\sigma }\mathrm{\left(}{\mathit{M}}_{\mathrm{1}}\mathrm{\right)}}{\mathit{\sigma }\mathrm{\left(}{\mathit{M}}_{\mathrm{0}}\mathrm{\right)}}\right]}^{{\mathit{\gamma }}_{\mathrm{1}}}\hspace{0.17em}{\left[\frac{\mathit{\delta }\mathrm{(}{\mathit{z}}^{\mathrm{)}}}{\mathit{\sigma }\mathrm{\left(}{\mathit{M}}_{\mathrm{0}}\mathrm{\right)}}\right]}^{{\mathit{\gamma }}_{\mathrm{1}}}\mathit{,}$(8)where δ(z) is the linear density threshold for collapse and $\mathit{\sigma }{}^{\mathrm{\right(}}\mathit{M}^{\mathrm{\left)}}$ denotes the rms linear density fluctuation extrapolated to redshift z = 0 in spheres containing mass M (Lacey & Cole 1993). G₀, γ₁ and γ₂ are free empirical parameters of this ansatz, which PCH constrained by fitting their algorithm to merger trees extracted from the Millennium simulation (Springel et al. 2005).

While testing our implementation of the PCH algorithm, we observed that the algorithm calibrated with the best-fitting parameters given by PCH still notably under-predicts the abundance of high-mass haloes (~10¹⁵ M_⊙/h). More accurate results can be achieved with the parameters given by Benson (2008), who employed a slightly different fitting procedure. In this work, we thus used the values derived by Benson (2008): G₀ = 0.605, γ₁ = 0.375 and γ₂ = −0.115.

The performance of our algorithm for projecting merging galaxy clusters onto the past null cone is visualised in Fig. 1. We show a random halo catalogue in the redshift range 0.5 <z < 1.0 (full sky) that was generated with the PCH merger algorithm. As a reference, we also plot the theoretically expected values of the mass functions of Press-Schechter, Tinker, and Sheth-Tormen. The initial haloes at redshift z = 0 were drawn from the Tinker mass function. Figure 1 confirms that the sampled catalogue agrees well with the expectations of the Tinker mass function at the considered redshift range. Yet, even with the parameters from Benson (2008), the generated mass function falls slightly short of high-mass haloes. The deviation is only moderate, however, and to fix this problem, one would have to run large N-body simulations containing sufficiently many clusters with masses M > 10¹⁵ M_⊙/h and re-fit the PCH algorithm. Clearly, this would go beyond the scope of the present work. For the following sections, it is only important to remember that the generated halo catalogues lack a few high-mass objects, which renders our estimates of the order statistics more conservative. Moreover, if one is only interested in haloes with redshifts z ≫ 0, the problem of too few high-mass haloes can be mitigated by starting the merger tree simulations at higher initial redshifts, z_ini > 0.

We cannot use an independent approach here to sample single haloes (without mergers) because the PCH algorithm does not exactly reproduce the Tinker mass function. Whenever we compute the statistics of Einstein radii neglecting cluster mergers, we simply take the catalogues generated by the merger algorithm instead and ignore the spatial correlations of the haloes, that is, we simply treat all haloes as isolated objects.

Fig. 1 Comparison between a Monte-Carlo sampled halo catalogue in the redshift range 0.5 <z < 1.0 (full sky) and the corresponding theoretically expected values of different halo mass functions. The random halo catalogue was generated with the merger tree algorithm of Parkinson et al. (2008). The red dashed, solid black, and green dotted curves indicate the mass functions proposed by Press & Schechter (1974; PS), Tinker et al. (2008; Tinker), and Sheth & Tormen (1999; ST).

3.2. Suggestions to reduce the computing time

Increasing computational complexity is one main difficulty in estimating the order statistics of the largest Einstein radii including clusters mergers. As was shown in Paper II, about 10³ mock universes have to be sampled for accurate fits of the GEV distribution. For each realisation, typically ~10⁶ merger trees need to be simulated, and subsequently, their strong-lensing signals have to be evaluated. This whole procedure is computationally quite demanding.

However, there are several steps that greatly reduce the required wall-clock time. Since future studies will probably face similar problems, it might be useful to list some of these steps here:

• Mass thresholds: as shown previously (and again verified for this work), only massive haloes contribute to the order statistics of the largest Einstein radii. We therefore only sample haloes with masses M > 10¹⁴ M_⊙/h at the initial redshift. While simulating their merger trees, we discard all progenitors whose mass falls below 5 × 10¹² M_⊙/h, because subhaloes at least need to exceed 5% of the main halo mass to notably perturb the Einstein radius (cf. Paper I). As described above, for each initial halo, our merger routine projects a list of progenitors onto the past null cone. If the total mass of all progenitors falls below 10¹⁴ M_⊙/h, we discard the system. Otherwise, we sort the progenitors in descending mass order, determine the Einstein radius of the most massive halo, and flag all progenitors that have been enclosed by the tangential critical curve. Next, we check if the total mass of the remaining unenclosed progenitors is larger than 10¹⁴ M_⊙/h. If that is the case, we again sort them in descending mass order and compute the next Einstein radius. We repeat this procedure until no relevant unenclosed progenitors remain.

  Fig. 2 Mass and redshift distributions of the gravitational lenses that produce the largest and the twelfth-largest Einstein radius. The distributions were extracted from 2000 mock realisations of Einstein radii in the redshift interval 0.5 ≤ z ≤ 1.0 on the full sky. The calculations incorporated the impact of cluster mergers. The plotted mass denotes the total mass of all haloes enclosed by the corresponding tangential critical curves.

• Maximum separation: when computing the tangential critical curve that surrounds a certain progenitor, we only take the neighbouring haloes into account (superposition of deflection angles) whose distance is smaller than the sum of both virial radii, d ≤ (r_vir,1 + r_vir,2). Haloes that are farther away can safely be neglected because they do not notably perturb the lensing signal.

• Einstein radius threshold: by experience, we know that the nth largest Einstein radius of a mock realisation will certainly be larger than a lower threshold ${\mathit{\theta }}_{\mathrm{eff}}^{\min }$. Then, before computing the detailed shape of the tangential critical curve of a lens system, we first quickly estimate the size of the resulting Einstein radius using the following steps: We loop over all relevant haloes in the field of view and treat them as isolated gravitational lenses. Because of their triaxial density profile, the projected surface mass density is ellipsoidal. We determine the major axis of the isodensity ellipses and compute the radius ${\mathit{\theta }}_{\mathrm{E}}^{\mathrm{est}}$ at which the tangential eigenvalue λ_t vanishes. A circle with radius ${\mathit{\theta }}_{\mathrm{E}}^{\mathrm{est}}$ encloses the full tangential critical curve of the halo and accordingly certainly over-estimates the real effective Einstein radius, ${\mathit{\theta }}_{\mathrm{eff}}^{\mathrm{est}}\mathit{>}{\mathit{\theta }}_{\mathrm{eff}}$. If the sum of all estimated Einstein radii (plus some tolerance) is smaller than the threshold ${\mathit{\theta }}_{\mathrm{E}\mathit{,}\mathrm{eff}}^{\min }$, we can safely discard the system and skip the expensive detailed lensing computations.

• Shoelace formula: the effective Einstein radius is derived from the area enclosed by the tangential critical curve. This area A can efficiently be computed using Stoke’s theorem in two dimensions $\mathit{A}\mathrm{=}{\mathrm{\int }}_{\mathit{A}}\mathrm{d}\mathit{A}\mathrm{=}{\mathrm{\int }}_{\mathit{A}}\mathrm{rot}\hspace{0.17em}{v}\mathrm{d}\mathit{A}\mathrm{=}{\mathrm{\int }}_{\mathit{\partial A}}{v}\hspace{0.17em}\mathrm{d}\mathit{s,}$(9)where v must be a vector field with rot v = − ∂_yv_x + ∂_xv_y = 1, such as v = 1 / 2 (− y,x)^T. The calculation of the area thus reduces to a simple line integral, which is far cheaper than standard connected-component labelling algorithms. Practically, the line integral is implemented by means of a loop over all boundary points, summing up the contributions using the so-called shoelace formula.

• Parallelisation: the calculations can optimally (and trivially) be parallelised with OpenMP/MPI. The wall-clock time decreases inversely proportional to the number of CPUs used.

4. Extreme value and order statistics with cluster mergers

We now study the impact of cluster mergers on the order statistics of the largest Einstein radii in a certain cosmological volume. For comparison with our previous work (Paper III), we exemplarily consider full-sky realisations of the cluster population in the redshift range 0.5 ≤ z ≤ 1.0, which contains the twelve high-redshift (z > 0.5) MACS clusters analysed by Zitrin et al. (2011). In agreement with these authors, we assume a constant source redshift of z_s = 2.0 throughout this work. We sampled 2000 mock realisations of the cluster population and collected the largest Einstein radius of each run. These data allowed us to accurately estimate the order statistics of the twelve largest Einstein radii.

Figure 2 shows the mass and redshift distributions of the clusters that produce the largest and the twelfth-largest Einstein radius. Because our analysis also includes merging clusters, it is important to stress that the mass referred to in Fig. 2 is the total mass of all haloes enclosed by the corresponding tangential critical curves. Generally, we note that cluster mergers do not significantly alter the mass and redshift distributions of the strongest gravitational lenses. The plots shown here are almost identical to those presented in Paper III, which did not include the impact of cluster mergers. This is why we only show the plots of the first and the twelfth order here. Most importantly, however, Fig. 2 confirms that the lower mass threshold of M_min = 10¹⁴ M_⊙/h discussed in Sect. 3.2 is well justified.

Fig. 3 Cumulative distribution functions (CDF) of the twelve largest Einstein radii, including the impact of cluster mergers. The CDFs were extracted from 2000 mock realisations of Einstein radii in the redshift range 0.5 ≤ z ≤ 1.0 on the full sky. The black dashed curve indicates the fit of the general extreme value (GEV) distribution to the CDF of the largest Einstein radius.

Figure 3 shows the CDFs of the twelve largest Einstein radii extracted from the 2000 mock realisations sampled. While cluster mergers notably shift the CDFs to larger Einstein radii, the general characteristics discussed in Paper III are conserved. The CDFs steepen with increasing order, implying that higher orders are, in principle, more constraining. The black dashed line indicates the excellent fit of the GEV to the CDF of the 2000 maxima sampled.

Fig. 4 Impact of cluster mergers on the cumulative distribution function (CDF) and the probability density function (PDF) of the largest Einstein radius in the redshift range 0.5 ≤ z ≤ 1.0 on the full sky. The general extreme value (GEV) distribution was fitted to the empiric CDF of maxima collected from 2000 mock realisations of the cosmological volume considered.

Next, we quantify the impact of cluster mergers on the GEV distribution of the largest Einstein radius. With mergers, the best-fitting parameters of the GEV distribution are given by (α,β,γ) = (50.85 ± 0.02,6.61 ± 0.03,0.100 ± 0.005). Excluding mergers, the best-fitting parameters are (α,β,γ) = (45.88 ± 0.02,6.52 ± 0.02,0.105 ± 0.004). The CDFs and the PDFs of these distributions are compared in Fig. 4. In both cases, the shape parameters γ are positive, indicating that the GEV distributions are bounded from below. Most importantly, mergers increase the location parameter α by ~5″ (~11%), which means that the GEV distribution is significantly shifted to higher values. The mode of the GEV distribution – which is the most likely value – is given by ${\mathit{x}}_{\mathrm{0}}\mathrm{=}\mathit{\alpha }\mathrm{+}\frac{\mathit{\beta }}{\mathit{\gamma }}\left[{\left(\mathrm{1}\mathrm{+}\mathit{\gamma }\right)}^{\mathrm{-}\mathit{\gamma }}\mathrm{-}\mathrm{1}\right]\mathit{,}$(10)and increases from 45.2″ without mergers to 50.2″ with mergers. The expectation value value of the GEV distribution reads ${\mathit{E}}_{\mathrm{GEV}}\mathrm{=}\mathit{\alpha }\mathrm{-}\frac{\mathit{\beta }}{\mathit{\gamma }}\mathrm{+}\frac{\mathit{\beta }}{\mathit{\gamma }}\mathrm{\Gamma }\left(\mathrm{1}\mathrm{-}\mathit{\gamma }\right)\mathit{,}$(11)where Γ denotes the Gamma function. Again, the expectation value also increases from 50.4″ to 55.4″.

Fig. 5 Box-and-whisker diagram comparing the order statistics of the twelve largest Einstein radii including and excluding cluster mergers. For each order, in the upper plot the red lines indicate the median (Q50), the blue bordered grey boxes give the interquartile range (IQR), and the black whiskers mark the range between the 2 and 98 percentile (Q2, Q98) of the order statistics including the impact of cluster mergers. The green errors bars indicate the 68% confidence intervals of the order statistics excluding mergers. The lower plot shows the percentage increase of the medians caused by mergers.

The impact of cluster mergers on the order statistics of the twelve largest Einstein radii is summarised by the box-and-whisker diagram presented in Fig. 3. We chose this representation to simplify the comparison with the results of Paper III, but also because these diagrams compactly visualise important properties of statistical distributions. All medians of the twelve (six) highest orders are higher than 36″ (43″), indicating that the ΛCDM model predicts a dozen Einstein radii as large as ~35″−55″ in the considered cosmological volume. If the 98 percentile were defined as exclusion criterion, one would need to observe approximately ten Einstein radii with θ_eff ≳ 41″, five with θ_eff ≳ 48″ or one large system with θ_eff ≳ 80″, to claim tension with the expectations of the ΛCDM model. Current observational data certainly do not exceed these expectations (Zitrin et al. 2011). Finally, the lower panel of Fig. 5 reveals that not only the expected value of the largest Einstein radius, but also those of the higher orders shown in the diagram, increases by 10−12% because of the impact of cluster mergers.

In summary, the results of this section agree well with the findings of Paper I: cluster mergers are significant for the statistics of the strongest gravitational lenses. As a rule of thumb, mergers increase the expected values of the largest Einstein radii by ~10%.

5. Comparison with observational data

5.1. Galaxy cluster MACS J0717.5+3745

As discussed in Paper II, the X-ray luminous galaxy cluster MACS J0717.5+3745 is a remarkable system: it is extremely massive, actively merging, and exhibits the largest known Einstein radius (Zitrin et al. 2009, 2011; Limousin et al. 2012). The cluster was re-observed as part of the Cluster Lensing And Supernova survey with Hubble (CLASH; Postman et al. 2012), a 524-orbit Hubble Space Telescope (HST) multi-cycle treasury program. Thereupon, Zitrin and collaborators revised their initial mass model from 2009 as published in Medezinski et al. (2013). The latest strong-lensing analysis of MACS J0717.5+3745 was conducted as part of the Hubble Frontier Fields program¹, yielding an Einstein radius of θ_eff = (55 ± 6)″ for a source redshift z_s = 2.0 (Zitrin et al., in prep.; priv. comm.).

Although our case study in Paper II already revealed that the Einstein radius of MACS J0717.5+3745 is not in tension with the expectations of the ΛCDM model, for completeness, we briefly update our results with the latest cosmological parameters from Planck (Planck Collaboration XVI 2014), and include the impact of cluster mergers. The cosmological volume analysed in the previous section was purposely chosen to contain the twelve MACS clusters discussed before. We can therefore use the best-fitting GEV parameters derived in Sect. 4.

Figure 6 confirms that cluster mergers further confirm our main conclusion from Paper II: the large Einstein radius of MACS J0717.5+3745 clearly does not exceed the theoretical expectations of the ΛCDM model. In contrast, the probability for observing a maximum Einstein radius even larger than 55″ amounts to ~42%. Neglecting cluster mergers, this probability decreases to ~24% (as can easily be verified by evaluating the GEV distribution with the best-fitting parameters given in Sect. 4). The 98 percentile of the GEV distribution is located at θ_eff = 82″, indicating that extraordinarily strong lenses would have to be observed to claim disagreement with the ΛCDM model. Interestingly, when cluster mergers are included, the Einstein radius of MACS J0717.5+3745 (accidentally) coincides with the expectation value θ_eff = 55″ of the theoretical GEV distribution (cf. Eq. (11)).

Fig. 6 Comparison between the theoretically expected general extreme value (GEV) distribution of the largest Einstein radius and the Einstein radius observed in the galaxy cluster MACS J0717.5+3745. The GEV distribution was fitted to the cumulative distribution of the largest Einstein radii extracted from 2000 mock realisations of the cluster population in the redshift range 0.5 ≤ z ≤ 1.0 (full sky).

There is no need to revise the order statistics of the twelve MACS clusters discussed in Paper III, because we already showed that most Einstein radii of the higher orders lie below the theoretically expected values, even if cluster mergers are excluded. As indicated by Fig. 5, the inclusion of mergers would not add any new insight.

5.2. SDSS clusters

Hao et al. (2010) applied a cluster-finding algorithm to the Sloan Digital Sky Survey (SDSS) Data Release 7 (DR7) data (Abazajian et al. 2009) and assembled a large optical galaxy cluster catalogue consisting of over 55 000 rich clusters in the redshift range 0.1 ≤ z ≤ 0.55. Zitrin et al. (2012) proposed a new method for performing a simplistic strong-lensing analysis of these clusters in an automated way, and applied their technique to a subsample of 10 000 clusters to estimate their Einstein radii.

The method proposed by Zitrin et al. (2012) is based on the simple assumption that the light distribution observed in galaxy clusters generally traces their mass distribution well. The starting points are the red cluster member galaxies, which are assigned simple parametric mass profiles. The superposition of these individual profiles represents the galaxy component of the mass distribution. The dark matter component is constructed by smoothing the galaxies’ distribution with a two-dimensional cubic spline interpolation. The sum of the two components serves as an indicator for the total projected matter density map of the cluster. The crucial point of the method is the calibration of the mass-to-light (M/L) ratio, which sets the normalisation of the projected surface mass density. Zitrin et al. (2012) calibrated this M/L ratio using a subsample of ten well-studied galaxy clusters that were covered by both high-quality HST images and the SDSS, and assumed that this calibration is (approximately) valid for the entire SDSS sample. Using this procedure, Zitrin and co-workers automatically (and blindly) processed the 10 000 SDSS clusters, derived simple lens models from the photometry of the cluster member galaxies and estimated the corresponding Einstein radii for sources at redshift z_s = 2.0.

To compare their results with the theoretically expected order statistics of the largest Einstein radii, we sampled 2000 mock realisations of the redshift range 0.1 ≤ z ≤ 0.55, assuming full sky coverage. Note that this choice was purposely too optimistic, because the SDSS DR7 only covered approximately one fourth of the full sky and, additionally, Zitrin et al. (2012) analysed only ~20% of the discovered clusters to reduce the required computing time. Our theoretical estimates should therefore be expected to exceed the distribution of Einstein radii determined for the subsample of SDSS clusters.

However, Fig. 7 and Table 1 reveal exactly the opposite trend. The mean estimated Einstein radii extracted from the SDSS sample significantly exceed the expected values, and the difference increases with higher orders. Starting with the second order, the mean estimated Einstein radii even exceed the 98 percentiles of the theoretical distributions. On the other hand, the errors of the SDSS Einstein radii are still large. Zitrin et al. (2012) suggested to assume an uncertainty of at least ~17% for the 1σ-boundary. From this point of view, all observed Einstein radii agree with the theoretically expected values within the 3σ confidence interval. However, we recall again that our theoretical estimates were computed for a substantially larger cosmological volume.

Fig. 7 Box-and-whisker diagram comparing the theoretically expected order statistics of the twelve largest Einstein radii in the redshift range 0.1 ≤ z ≤ 0.55 (full sky) and the largest Einstein radii in the SDSS sample analysed by Zitrin et al. (2012). For each order, the red lines indicate the median (Q50), the blue bordered grey boxes give the interquartile range (IQR), and the black whiskers mark the range between the 2 and 98 percentile (Q2, Q98) of the order statistics. The green error bars indicate the estimated 68% confidence intervals of the Einstein radii in the SDSS sample.

These results may imply that the observed gravitational lenses at low redshifts are somewhat stronger than theoretically expected, which would agree with the findings of Horesh et al. (2011), for example. Given the current state of the observational data, however, this conclusion would certainly be premature. The uncertainties of the method used by Zitrin et al. (2012) are still significant. Most importantly, the reliability of the assumed functional form for the M/L ratio needs to be verified and re-calibrated with a larger sample of well-studied gravitational lenses. Furthermore, Zitrin and co-workers mainly focused on the universal distribution of Einstein radii, instead of analysing the strongest lenses of their sample in detail (A. Zitrin, priv. comm.).

More meaningfully, we therefore conclude that the results of this section might indicate a discrepancy between theory and observations at low redshifts, and that it should certainly be interesting to carefully re-analyse the SDSS sample. Given the current state of the data and also the theoretical uncertainties influencing the modeling of the order statistics, there is, however, no reliable statistical evidence for claiming that these observations seriously challenge the predictions of the standard cosmological model.

6. Conclusions

This series of papers focused on the question whether or not the strongest observed gravitational lenses exceed the theoretical predictions of the standard cosmological model. In this last work, we combined all previously developed techniques and calculated the impact of cluster mergers on the order statistics of the largest Einstein radii. Moreover, we compared the theoretical results with observational data at different redshifts.

In the first part, we demonstrated that cluster mergers shift the GEV distribution of the largest Einstein radius to significantly higher values, typically increasing the location parameter and expected values by ~ 10%. Furthermore, we showed that the order statistics of the n largest Einstein radii are also shifted by a similar amount. This confirmed our findings of Paper I, where we argued that cluster mergers are particularly important for the statistics of the strongest gravitational lenses.

In the second part, we compared the order statistics of the largest Einstein radii, including the impact of cluster mergers, with recent observational data. As already shown in Papers II and III, at redshifts z > 0.5, we see no evidence for a tension between the strength of observed gravitational lenses and the theoretical predictions of the ΛCDM model. On the contrary, we actually find that the largest known Einstein radius at redshifts z > 0.5, which was observed in the galaxy cluster MACS J0717.5+3745 (Zitrin et al. 2009, 2011; Limousin et al. 2012; Medezinski et al. 2013), agrees excellently well with the theoretical expectation value.

At redshifts z < 0.5, the situation is a little more controversial. We compared the largest Einstein radii of the 10 000 SDSS clusters analysed by Zitrin et al. (2012) to the theoretically expected order statistics, finding that the observed gravitational lenses in this redshift range appear to be stronger than expected. However, in this context, it is important to stress that the errors of the Einstein radii estimated for the SDSS sample are still large. Thus, given the premature state of current data, there is still no reliable statistical evidence for claiming disagreement with the ΛCDM model. Nevertheless, the SDSS sample might contain extraordinarily strong gravitational lenses, and it should certainly be interesting to carefully re-analyse these systems in future studies.

In this last Paper of our series, it is probably worth to highlight some important conclusions that should be drawn from our studies.

Firstly, we explained that the theory of extreme value statistics – and more generally, the order statistics of the n largest Einstein radii – provides a rigorous mathematical framework for discussing the likelihood of the largest observed Einstein radii, and also illustrated how exclusion criteria can be formulated. We developed fast, semi-analytic methods for computing theoretical forecasts, particularly emphasising that including cluster mergers is significant for accurate estimates of the order statistics. For this reason, future semi-analytic studies inevitably need to take the impact of cluster mergers into account before questioning the validity of the ΛCDM model.

Secondly, we extensively discussed that the order statistics of the n largest Einstein radii is strongly influenced by many model parameters. This property can both be interpreted as a strength and a weakness. It can be seen as a strength because it renders the order statistics a sensitive cosmological probe that reacts to changes of the cosmological parameters (e.g., Ω_m and σ₈), details of structure formation (e.g., merger rates and halo mass function) and properties of triaxial dark matter haloes (e.g., concentration and axis ratios). On the other hand, it can be seen as a weakness because due to this sensitivity, the order statistics is currently still subject to many theoretical uncertainties. Given the current state of the theory, it is certainly difficult to formulate reliable exclusion criteria that seriously challenge the underlying cosmological model.

Finally – and this is probably the most important result of this series of papers – we can conclude that there is no statistical evidence at present for claiming that the largest observed Einstein radii exceed the theoretical expectations of the ΛCDM model.

Acknowledgments

We thank Adi Zitrin for valuable comments and for kindly providing the data set of 10 000 Einstein radii. We also thank Shaun Cole and Andrew J. Benson for helpful discussions. This work was supported in part by the German Deutsche Forschungsgemeinschaft, DFG (DFG), and in part by contract research Internationale Spitzenforschung II-1 of the Baden-Württemberg Stiftung. Most simulations required for this work were performed on the bwGRiD cluster (http://www.bw-grid.de), member of the German D-Grid initiative, funded by the Ministry for Education and Research (Bundesministerium für Bildung und Forschung) and the Ministry for Science, Research and Arts Baden-Wuerttemberg (Ministerium für Wissenschaft, Forschung und Kunst Baden-Württemberg).

References

All Tables

All Figures

Fig. 1 Comparison between a Monte-Carlo sampled halo catalogue in the redshift range 0.5 <z < 1.0 (full sky) and the corresponding theoretically expected values of different halo mass functions. The random halo catalogue was generated with the merger tree algorithm of Parkinson et al. (2008). The red dashed, solid black, and green dotted curves indicate the mass functions proposed by Press & Schechter (1974; PS), Tinker et al. (2008; Tinker), and Sheth & Tormen (1999; ST).

In the text

Fig. 2 Mass and redshift distributions of the gravitational lenses that produce the largest and the twelfth-largest Einstein radius. The distributions were extracted from 2000 mock realisations of Einstein radii in the redshift interval 0.5 ≤ z ≤ 1.0 on the full sky. The calculations incorporated the impact of cluster mergers. The plotted mass denotes the total mass of all haloes enclosed by the corresponding tangential critical curves.

In the text

	Fig. 3 Cumulative distribution functions (CDF) of the twelve largest Einstein radii, including the impact of cluster mergers. The CDFs were extracted from 2000 mock realisations of Einstein radii in the redshift range 0.5 ≤ z ≤ 1.0 on the full sky. The black dashed curve indicates the fit of the general extreme value (GEV) distribution to the CDF of the largest Einstein radius.
In the text

	Fig. 4 Impact of cluster mergers on the cumulative distribution function (CDF) and the probability density function (PDF) of the largest Einstein radius in the redshift range 0.5 ≤ z ≤ 1.0 on the full sky. The general extreme value (GEV) distribution was fitted to the empiric CDF of maxima collected from 2000 mock realisations of the cosmological volume considered.
In the text

Fig. 5 Box-and-whisker diagram comparing the order statistics of the twelve largest Einstein radii including and excluding cluster mergers. For each order, in the upper plot the red lines indicate the median (Q50), the blue bordered grey boxes give the interquartile range (IQR), and the black whiskers mark the range between the 2 and 98 percentile (Q2, Q98) of the order statistics including the impact of cluster mergers. The green errors bars indicate the 68% confidence intervals of the order statistics excluding mergers. The lower plot shows the percentage increase of the medians caused by mergers.

In the text

Fig. 6 Comparison between the theoretically expected general extreme value (GEV) distribution of the largest Einstein radius and the Einstein radius observed in the galaxy cluster MACS J0717.5+3745. The GEV distribution was fitted to the cumulative distribution of the largest Einstein radii extracted from 2000 mock realisations of the cluster population in the redshift range 0.5 ≤ z ≤ 1.0 (full sky).

In the text

Fig. 7 Box-and-whisker diagram comparing the theoretically expected order statistics of the twelve largest Einstein radii in the redshift range 0.1 ≤ z ≤ 0.55 (full sky) and the largest Einstein radii in the SDSS sample analysed by Zitrin et al. (2012). For each order, the red lines indicate the median (Q50), the blue bordered grey boxes give the interquartile range (IQR), and the black whiskers mark the range between the 2 and 98 percentile (Q2, Q98) of the order statistics. The green error bars indicate the estimated 68% confidence intervals of the Einstein radii in the SDSS sample.

In the text