The XXL survey

In short, the Xamin pipeline tests four models to characterise the detected sources, generating likelihood estimates for point, extended, double point sources, and extended plus point source; this latter model, denoted AC, is intended to flag extended sources significantly contaminated by a central AGN. The coordinates of the X-ray source presented in Sect. 3 are based on the centre of the best-fit model. Cluster sources are further classified into C1 and C2 on the basis of pipeline parameters extent and extent_likelihood. The C1 sample corresponds to an almost pure sample of bright clusters, while the C2 sample, fainter, allows for 50% of misclassified point sources (i.e. up to 50% of the sample might be point sources; see Pacaud et al. 2006; XXL Paper XXIV). False C2 are routinely excluded by the examination of X­-ray/optical overlays for cluster candidates below z=1.

We mention that the choice of the numerical pipeline parameter values used to define the C1 and C2 criteria are still those based on detections performed with the V3 pipeline from simulated individual XMM observations. These criteria will be revised when the final V4 is fully validated and applied to mosaic simulations. However, we do not expect drastic changes in the class parameters, since they are based on likelihoods. This situation does not impact the current study, because it does not explicitly involve the cluster selection function at any stage.

The VIDEO observations consist of IR imaging undertaken with the VISTA telescope in the $\mathrm{YJHK_{s}}$ photometric bands. In the XMM-SERVS field, these observations reach 5$\sigma$ depths of at least 25.1, 24.7, 24.2, and 23.8 mag within 2 arcsecond circular apertures for $\mathrm{YJHK_{s}}$ respectively (Adams et al. 2020). The VIDEO catalogue also contains additional imaging data consisting of the Canada-France-Hawaii Telescope Legacy Survey Deep-1 field (CFHTLS-D1) and of the deep ‘layer’ of the Hyper Suprime-Camera (HSC) Subaru Strategic Program (HSC-SSP, Aihara et al. 2018a, b). The ultra-deep ‘layer’ of HSC-SSP overlaps with the XMM1 field. The photometric redshift analysis included in the catalogue employs an i-band selected source list where photometry in additional bands is obtained by applying SExtractor to astrometrically-matched pixel data at other wavelengths. In addition, we employ HSC-SSP $\mathrm{iz}$ and VIDEO $\mathrm{JK_{s}}$ photometry from this catalogue to study the properties of candidate cluster member galaxies directly. The HSC-SSP deep data have 5$\sigma$ limiting magnitudes of 25.4 and 24.6 respectively in the $i$ and $z$ bands, and the ultra-deep data have limiting magnitudes of 26.4 and 26.3 (Adams et al. 2020).

Photometric redshifts for sources in the VIDEO catalogue are computed using the LePhare photometric redshift code (Ilbert et al. 2006). The code employs the COSMO template set (Ilbert et al. 2009), including 32 templates from Polletta et al. (2007) and from Bruzual & Charlot (2003). Dust attenuation follows a Calzetti et al. (2000) law, and intergalactic medium absorption treatment is based on Madau (1995). Further details are provided by Adams et al. (2020).

We perform four steps to select candidate clusters from our X-ray sample, following a similar procedure to that of Willis et al. (2013). Each step is summarised below:

1. 1.

   Convolve the photometric redshift histogram for bright galaxies with a matched Gaussian filter chosen to match the properties of redshift peaks of spectroscopically confirmed clusters. Galaxies are selected to be brighter than the characteristic luminosity, $\mathrm{L^{\ast}}$, along the line-of-sight of each X-ray source (i.e. a one arcmin radius aperture centred on the X-ray best-fit model centre).

2. 2.

   Identify overdensities corresponding to a bright galaxy excess of $\gtrsim 4$, and display the position of potential members on colour images together with X-ray emission contours.

3. 3.

   Examine the $i-z$ and $z-J$ colour-magnitude diagrams of each candidate cluster.

4. 4.

   Employ a Gaussian model of the photometric redshift distribution of selected overdensities to estimate a refined mean cluster redshift and its standard deviation.

Figure 2 presents a visual summary of these four steps. Similar images of the other candidate clusters are presented in Appendix B.

The first step of the identification process is to select galaxies and active galactic nuclei (AGN) that might be associated with each X-ray sources. We select the galaxies and AGN sources employing their goodness-of-fit when compared to stellar, galactic, and AGN templates (see Jarvis et al. 2013). Because galaxies in clusters are more likely to host radio-loud AGN than field galaxies (Best et al. 2007), we keep AGN and discard only the star-like objects. Then, we select objects within one arcmin of the considered X-ray detection. We refer to these objects as the field-of-view galaxies.

We select bright galaxies in the field-of-view employing luminosity function arguments: we compute the apparent magnitude of $\mathrm{M^{\ast}}$, assuming $\mathrm{M^{\ast}}=-22.26$ in the $\mathrm{K_{s}}$ band using Cirasuolo et al. (2010) with no evolution from $z\sim 3$ to $z\sim 0$ and a $k$-correction described by the bandwidth term. Galaxies brighter than the expected apparent $\mathrm{m^{\ast}}$ magnitude at their photometric redshift are retained. We apply a further photometric cut based on the 5$\sigma$ depth of VIDEO $\mathrm{K_{s}}$ band, discarding any galaxy fainter than 23.8. The latter cut is important beyond $z\sim 1.4$, where $\mathrm{m^{\ast}}$ is fainter than the VIDEO 5$\sigma$ limit.

Selected field-of-view galaxies are then binned in photometric redshift space, over the interval $0.2<z_{phot}<3.2$. Galaxies are sampled in bins of 0.04 in photometric redshift since this is the redshift iteration step employed in the photometric redshift analysis. We then use the catalogue distribution to perform a background subtraction, applying the same selection steps described above and scaling the number distribution by the relative size of our field-of-view compared to that of the full the catalogue, i.e.

where $N_{FOV}$ is the number of galaxies in a particular redshift bin in our field-of-view, $N_{cat}$ is the number of galaxies listed in the catalogue in the corresponding redshift bin after applying the same cuts, $r_{FOV}$ is the one arcminute radius we use to select galaxies associated with an X-ray detection, and $A_{cat}$ is the catalogue area.

To identify structures in the photometric redshift histogram along the line-of-sight to each extended X-ray source we employ a matched Gaussian filter with a FWHM equal to 0.12 in redshift (i.e. 3 bins). The properties of the Gaussian profile are based upon the unfiltered redshift peaks associated with spectroscopically confirmed clusters.

We perform a visual inspection of all C1, C2, and AC X-ray sources that display a signal consistent with $>4$ galaxies at a single photometric redshift. We typically employ $riH$ images from CFHTLS and VIDEO with candidate members indicated in addition to X-ray emission contours (see Fig. 2). We further generate $i-z$ and $z-J$ colour magnitude diagrams of each candidate cluster to determine if a red sequence is present.

The overdensity finding method provides a first estimate of the candidate cluster redshift based on the median redshift of the highest bins (see Fig. 2, bottom left panel). To refine this estimate, we model each candidate redshift signal as a Gaussian employing the 13 central bins of the non-filtered redshift signal. We then use the Gaussian mean as the cluster redshift and the standard deviations as a estimate of the uncertainty.

We processed a total of 284 extended X-ray detections within the XXL-N/VIDEO region. This parent sample generated a sample of 35 candidate distant galaxy clusters. Table 1 and Fig. 3 present these clusters, each of which represents a detection with a significance of approximately four galaxies or more in a photometric redshift bin satisfying $z_{phot}\geq 0.8$ (see also Fig. 1). Of these 35 candidate clusters, ten are spectroscopically confirmed while 15 are presented here for the first time. Ten additional candidates have been previously identified as distant clusters but, to our knowledge, never spectroscopically confirmed (see Olsen et al. 2007; Finoguenov et al. 2010; Durret et al. 2011; Wen & Han 2011; Licitra et al. 2016; XXL Paper XX).

Nine of the confirmed cluster detections were confirmed by prior spectroscopy (e.g. Pierre et al. 2006; Willis et al. 2013; XXL Paper XX), including one at $z=1.803$ (Andreon et al. 2014). With the spectroscopic redshifts listed in the CESAM database¹¹1http://cesam.lam.fr/xmm-lss/ (XXL Paper XX), we were able to confirm one additional cluster (candidate 3), bringing the total number of confirmed cluster to ten. All of these X-ray detections meet our candidate clusters criteria.

Four confirmed distant clusters in this area (Pierre et al. 2006; Papovich et al. 2010; XXL Paper XX), are not part of our sample because their coordinates do not correspond to a V4 C1, C2, or AC detection. This is the case for a $z=1.62$ cluster (Papovich et al. 2010). Although IRC-0218A and our detections at similar redshifts possess comparable masses (IRC-0218A $M_{200}$ is $7.7\pm 3.9\times 10^{13}\leavevmode\nobreak\ \mathrm{M_{\odot}}$), IRC-0218A X-ray emission is completely dominated by a point source (Pierre et al. 2012).

We nevertheless applied our optical/IR detection criteria to those four clusters. Two clusters satisfied these criteria, including IRC-0218A, and were included in our photometric redshift accuracy assessment (see the following section). The clusters XLSS J022609.9-043120 and XLSSC 203 (see Pierre et al. 2006; XXL Paper XX) located respectively at $z=0.82$ and $z=1.077$, did not satisfy them. Because V4 works on co-­added and thus deeper images, we expect its source characterisation to be more reliable.

To estimate the reliability of our photometric redshift estimates for the candidate clusters, we compared the spectroscopic redshift of 12 confirmed clusters in the field (6 clusters listed in XXL Paper XX, four from Pierre et al. 2006; Papovich et al. 2010; Willis et al. 2013; Andreon et al. 2014, and 2 others confirmed clusters in the XXL-N/VIDEO area) to the photometric redshift generated by the cluster finding procedure. Figure 4 shows the result of this comparison. At $z\sim 1$, the differences between the photometric and spectroscopic redshifts fall well within the photometric redshift error estimates. These error bars represent the standard deviation of the best-fitting Gaussian model. For the two high-redshift clusters, the photometric redshifts seem to underestimate the spectroscopic values. Therefore, we calculate the root mean square (RMS) of the $z\sim 1$ and $z\gtrsim 1.5$ clusters separately. We obtained 0.02 and 0.14 respectively. For now on, we will use these RMS values as the uncertainties on the photometric redshifts.

Following a similar methodology than the second data release of XXL, we used scaling relations to provide mean parameters estimate for clusters for which the data quality is not enough to perform a direct spectral fit. A detailed description is provided in Sect. 4.3 of XXL Paper XX, and we provide a brief overview here. We estimated count-rates in the pn data in the [0.5-2] keV band within 300 kpc of the cluster centre using the Bayesian approach to fixed aperture photometry measurement outlined in Willis et al. (2018). We then converted this count-rate to the corresponding X-ray luminosity by adopting an initial gas temperature, a metallicity fixed to 0.3 times the solar value (as tabulated in Anders & Grevesse 1989), and the cluster spectroscopic (when available) or photometric redshift. With the same initial guess of the temperature, we estimate $r_{500,scal}$ from the mass-temperature relation constrained from a subset of 105 XXL clusters that have both measured HSC lensing masses and X-ray temperatures (see Umetsu et al. 2020). We stress that here we use a $M_{500}\leavevmode\nobreak\ |\leavevmode\nobreak\ T_{X}$ relation, obtained using the Bayesian regression scheme implemented in the LIRA package (Sereno 2016; Sereno et al. 2016), and not the $T_{X}-M_{500}$ relation reported in Umetsu et al. (2020). The luminosity is then extrapolated from 300 kpc to $r_{500,scal}$ assuming a $\beta$-model for the cluster emissivity with parameters $(r_{c},\beta)=(0.15r_{500,scal},2/3)$. Then a new temperature is evaluated using the best-fit result for the luminosity-temperature relation quoted in Table 6 of XXL Paper XX (XXL fit). The iteration on the gas temperature is stopped when the input and output values agree within 5%, and in general the process converges after 2-3 steps. The uncertainties on the derived parameters and in particular the masses, are obtained by propagation of errors on the scaling parameters, including the measured correlation among them.

There are 54 previously confirmed clusters either within or overlapping with the XXL-N/VIDEO field (XXL Paper XX). Of these, 47 are located at $z<0.8$ with the remaining seven clusters located at $z\geq 0.8$. These seven clusters correspond to a surface density of 1.6 distant clusters per square degree. This number represents a lower limit because some known clusters associated with an extended X-ray detection (e.g. JKCS 041, a z=1.803 confirmed cluster Andreon et al. 2014) are excluded of (XXL Paper XX) compilation.

Adding all our detections would bring up this number to approximately 8.2 clusters per square degree, with a flux limit of $1.7\times 10^{-15}\leavevmode\nobreak\ \mathrm{erg\leavevmode\nobreak\ cm{-2}\leavevmode\nobreak\ s^{-1}}$ in the [0.5-2] keV band (Chen et al. 2018). This represent 3.6 times and 0.51 times the surface densities reached by Willis et al. (2013) and Finoguenov et al. (2010) respectively, with depths of $1\times 10^{-14}\leavevmode\nobreak\ \mathrm{erg\leavevmode\nobreak\ cm{-2}\leavevmode\nobreak\ s^{-1}}$ and $2\times 10^{-15}\leavevmode\nobreak\ \mathrm{erg\leavevmode\nobreak\ cm{-2}\leavevmode\nobreak\ s^{-1}}$ in the [0.5-2] keV band.

The first step consists of defining the appropriate area within which to select field-of-view galaxies for each cluster. Some studies (e.g. Wetzel et al. 2012; Pintos-Castro et al. 2019) analyse the quenched/star-forming fraction up to several virial radii. However, since several of our fields appear to contain secondary overdensities, we choose to restrict our analysis closer to the cluster centre, namely within a radius of 1.35 Mpc around the X-ray coordinates. A radius of 1.35 Mpc corresponds to approximately 1.5 times the value of $r_{500}$ for a $1\times 10^{14}\leavevmode\nobreak\ \mathrm{M_{\odot}}$ galaxy cluster (Chen et al. 2007) and to an angle of 2.76 arcminutes at $z=1$. We further select galaxies with a photometric redshift between 0.60 and 2.04, both in the cluster field and background catalogue, and sample the resulting distributions into 0.04 wide bins in $i-z$ colour.

We then correct this field-of-view distribution using the same method as applied in Eq. 1. Visual inspection reveals that the resulting colour distribution is bimodal (see Fig. 5). We model this bimodal colour distribution using two Gaussian functions – one representing the red sequence and the other the blue cloud. In some fields we note that the background correction is either too large or too small to yield to a clear bimodal distribution, resulting in poor fits or even in a non-convergent fitting algorithm (candidate 7). In such cases, we adjusted the background correction by up to 20% before determining whether the adjusted background results in an improved fit. Figure 5 presents the effect of an unadjusted background on the colour distribution on two typical field. Although both the too high and too low cases are presented, the former concerns only four clusters. One of these clusters is confirmed (candidate 3) and two others (including the case presented in Fig. 5) are among the sample of other photometric studies (see Olsen et al. 2007; Wen & Han 2011; Licitra et al. 2016). Thus, most of them are unlikely to be false detections.

Three of the eight clusters with an increased background lie in overlaps between two VIDEO footprints and have thus, according to Jarvis et al. (2013) longer exposure times. Moreover, two additional clusters are in the XMM1 field, which is part of the ultra-deep layer of HSC-SPP survey and thus possess 1 to 1.5 mag deeper photometry in the i and z band (see Adams et al. 2020).

The quenched fraction is then computed as the integral of the best-fitting red sequence Gaussian profile, divided by the sum of the integrals of the blue cloud and red sequence Gaussian profiles.

Some clusters were removed from the analysis. We removed 1, 4, 10, 21, and 24 since there are indications of more than one overdensity along the line-of-sight to each of these X-ray sources. Candidate clusters at $z\geq 1.4$ were also removed, resulting in a sample of 21 candidate clusters.

Next, we investigate whether the Gaussian modelling produces consistent average star formation histories, i.e. whether the red sequence colour can be described adequately by a simple stellar population model. We used Flexible Stellar Population Models (FSPS; see Conroy et al. 2009; Conroy & Gunn 2010; Foreman-Macke et al. 2014), modelling the red sequence stellar population with an exponentially decreasing star formation rate characterised by an $e$-folding timescale $\tau$, displaying solar metallicity and no dust. We generated a grid of models with two free parameters: we tested 100 formation redshifts between 4 and 16, and 100 $\tau$ values between 0.1 and 2.5 Gyr for a total of 10000 models. The model providing the best fit has a formation redshift of 16 with a characteristic time of 1.19 Gyr and a reduced $\chi^{2}_{\nu}=0.76$, corresponding to a rejection probability of 24%. The red sequence displayed on the CMD plot in Fig. 2 is based on this model. The associated uncertainty is the standard deviation of the difference of red sequence colours, as predicted by the model, and the Gaussian modelling results. We also used this fit to compute the evolutionary and K-correction associated with the characteristic luminosity used in method 4.

The second method to identify quenched galaxies employs the same colour binning and background subtraction procedure as employed by Gaussian method described above. We specify the colour where the blue and red Gaussian models are equal as the boundary between the red sequence and the blue cloud. We further restrict the colour interval by rejecting any galaxy redder than 2.5 times the standard deviation of the red Gaussian or bluer than 2.5 times the standard deviation of the blue Gaussian. Where those boundaries fall within a colour bin we perform a fractional calculation within the affected bins.

The third method does not include the area-corrected background subtraction used in the first two methods but instead selects cluster members based on their photometric redshift. Any galaxy within the field-of-view and with a redshift consistent with the cluster mean redshift plus or minus 1.5 times its standard deviation (calculated in Sect. 2.4) is considered as a cluster member. We then selected the red sequence and blue cloud members using the colour boundaries defined in method two.

All the methods described above employ the 5$\sigma$ limit of the VIDEO catalogue in J and $\mathrm{K_{s}}$ bands as a magnitude selection threshold. The fourth quenching method is essentially the same as method 2 yet with an evolving magnitude limit based upon the values of $\mathrm{L^{\ast}}$ in J and $\mathrm{K_{s}}$ bands computed using our best-fitting FSPS model for the red sequence, i.e. including both an evolutionary and $k$-correction (see Fig. 6). As in Sect. 2.3, we assume a characteristic absolute magnitude of -22.26 at $z=0$ (see also Cirasuolo et al. 2010).

A comparison between method 1 and the other methods is presented in Fig. 7. Each method generates similar quenched fractions compared to method 1. Within the limited variation of the approaches taken by each method, this indicates that the quenched fractions are robust from method to method. The few clusters for which the quenched fractions obtained by method 1 are more than 1.5 times the standard deviation above or below the mean are identified and highlighted. For now on, we will focus on method 4 results.

Figure 6 shows the quenched fraction results of the fourth method and indicates that there is a wide variety of quenching levels within the interval $0.8\leq z\leq 1.2$. The wide range of computed quenched fractions is nominally consistent with the expectation expressed in Sect. 1 that the X-ray selection of galaxy clusters should be less biased to the properties of their member galaxies than optical/IR overdensity methods.

We next test whether, despite the range of quenched fractions, there is a consistent variation of quenched fraction as a function of cluster-centric distance within the sample of clusters as a whole. We computed the quenched fraction for each cluster within three equal radial annuli out to 1.35 Mpc. We then computed the mean cluster quenched fractions as a function of radius into two groups based upon the total quenched fraction, i.e. we stacked those clusters quenched at the level greater than the median quenched fraction into one group and those quenched at less than the median level into another.

We then investigated if the quenched fraction depends upon galaxy luminosity. We computed the quenched fraction in three bins of $\mathrm{J/K_{s}}$ luminosities (0.5 to 1 $\mathrm{L^{\ast}}$, 1 to 2 $\mathrm{L^{\ast}}$, and 2 to 4 $\mathrm{L^{\ast}}$). Once again, the results for each cluster are stacked on the basis of their overall quenching level. This test is limited to $z<1.4$ because beyond this 0.5 $\mathrm{L^{\ast}}$ falls below the 5$\sigma$ magnitude cut.

Figure 8 displays no significant trend between the mean quenched fraction and the cluster-centric radius, especially for the highly quenched half of the sample. We note however a slight increase of the quenched fraction toward the centre of the lowly quenched cluster, suggesting a weak dependence of the quenched fraction with the distance. The bottom panels of Fig. 8 show that more luminous galaxies are more quenched, although the variation is less pronounced in highly quenched clusters.

Figure 7 presents a wide variety of quenched fractions, with no clear trend, despite that, at higher redshifts, only the most luminous (massive) galaxies are visible. The fourth method of computing the quenched fraction incorporates a luminosity cut specifically intended to mitigate this bias. However, although this method generates slightly higher quenched fractions (except for the two lowest redshift clusters in our sample; see Fig. 7), it does not otherwise affect the apparent diversity of quenched fractions. There is thus no evidence of evolution with redshift in Fig. 6.

The results presented in Sect. 4.3 demonstrate a link between quenched fraction and the $\mathrm{K_{s}}$ band luminosity, although the link seems weaker for the highly quenched half of our sample. $\mathrm{K_{s}}$ band luminosities are a proxy for galaxy stellar masses (assuming dust extinction is negligible). Such a link has been observed before (e.g. Peng et al. 2010, 2012; Muzzin et al. 2012; Balogh et al. 2016; Kawinwanichakij et al. 2017; Jian et al. 2018; Pintos-Castro et al. 2019), although whether this is due to mass-dependent galaxy evolution and/or to environment remains an open question.

We observe no significant evidence that the quenching fraction varies with the cluster-centric distance, although the quenched fraction is slightly higher toward the centre of the lowly quenched clusters in our sample – an observation in agreement with previous studies (e.g. Muzzin et al. 2012; Pintos-Castro et al. 2019). The more quenched half of our sample possesses levels of quenching above 70% in all bins, which might be because only galaxies above $\mathrm{L^{\ast}}$ are included in this profile and that bright galaxies tend to be slightly more quenched at all radii.

The $J-K_{s}$ colours of the candidate BCGs in the distant cluster sample are not consistent with the properties of a single, passively-evolving stellar population. There is tentative evidence for a bimodal distribution of BCG colours. Separating the sample employing the fiducial stellar population model of Lidman et al. (2012), we find that the colours of the redder BCGs in our sample are consistent with an older, instantaneous burst of metal-rich stars whereas the bluer BCG colours are more diverse and possibly related to a younger component in the stellar population.

Despite this variety of colours, only one of the BCGs is bluer than the red sequence model computed in Sect. 4.1 (see also Fig. 9). This might indicate that BCG stellar populations have different properties (such as metallicity) than the average quenched galaxies in cluster, as suggested by some low-redshift studies (e.g. Von Der Linden et al. 2007; Loubser et al. 2009).

Figure 11 shows the degeneracy between metallicity and formation redshift as applied to the stellar population models describing the two BCG samples. We observe that the 98% confidence regions of the two BCG samples do not overlap, supporting the idea that the stellar populations of the blue and red BCGs are indeed different. The confidence intervals for the blue sample are centred on lower metallicities than the red sample regions while formation redshift also tends to be lower – although the confidence regions of both models extend over a wide range of redshifts.

Based on the stellar population analysis, one might infer that blue BCGs have formed more recently than the red ones. However, Li & Han (2007) report that when more than one stellar population is present, the age and metallicity obtained by colours alone might be biased toward the younger and more metal-rich stars. Therefore, the blue BCGs might be bluer because they have experienced either an extended star formation episode or more than one bursts in the past. In such scenarios, differences in the duration of the star formation or in the relative importances and epochs of the secondary bursts would account for the spread in colour observed in the bluer BCGs.

One factor that may trigger a star formation episode is a cooling flow. Several studies have reported the existence of blue, moderately star-forming BCGs in cool core clusters at low to moderate redshifts (e.g. Egami et al. 2006; Bildfell et al. 2008; Stott et al. 2008; Loubser et al. 2009, 2016; Pipino et al. 2009; Rawle et al. 2012; Green et al. 2016). Since X-ray selected clusters are biased toward relaxed clusters (Rossetti et al. 2017), this might partially explain why the blue BCGs represent a significant part of our sample.

Alternatively, statistical studies based on infrared (e.g. Webb et al. 2015; Bonaventura et al. 2017) or SZ-selected clusters (e.g. McDonald et al. 2016) have shown that in-situ star formation is an important mechanism at $z\gtrsim 1$, possibly triggered by galaxy interactions (McDonald et al. 2016). This suggests that a range of BCG colours might be part of every high-redshift sample, regardless of the sample selection.

One limitation of the current analysis is that the fitted stellar population models are unable to accommodate the slope of the high-redshift end of the red sample J-$\mathrm{K_{s}}$ colour variation with redshift (see Fig. 9) and that no BCGs are red beyond $z=1.05$. This suggest that the stars making up the red BCGs formed later or evolve faster than predicted by our model. The former is supported by the degeneracy between formation redshift and metallicity. The best-fitting formation redshift, $z_{form}\sim 12$, might be considered high compared with the predictions of the most recent simulations (e.g. Ragone-Figueroa et al. 2018; Rennehan et al. 2020), but – alternatively – consistent with recent observations (Hashimoto et al. 2018; Willis et al. 2020).

Dust extinction might be evoked to explain why these objects are so red. Indeed, a fit including dust provides an equivalent statistical description of the red BCGs (see Table 2), although we note that the colour of a model stellar population is degenerate between star formation history, metallicity, and dust obscuration. For this reason, we prefer a simple, dust free model for the analysis of the red BCG population.

Ultimately, the current data available for the BCGs sample are unable to provide a definitive explanation of the blue/red dichotomy. The acquisition of rest-frame optical spectroscopy of several BCGs of both groups is needed in order to get a more accurate and thorough picture of the star formation history of these objects (e.g. Lonoce et al. 2015, 2020; Belli et al. 2017; Saracco et al. 2019).