Star-forming and gas-rich brightest cluster galaxies at z ∼ 0.4 in the Kilo-Degree Survey

Brightest cluster galaxies (BCGs) are typically massive ellipticals at the centers of clusters. They are believed to experience strong environmental processing, and their mass assembly and star formation history are still debated. We have selected three star-forming BCGs in the equatorial ﬁeld of the Kilo-Degree Survey (KiDS). They are KiDS 0920 ( z = 0 . 3216), KiDS 1220 ( z = 0 . 3886), and KiDS 1444 ( z = 0 . 4417). We have observed them with the IRAM 30m telescope in the ﬁrst three CO transitions. We remarkably detected all BCGs at high signal-to-noise ratio, S / N (cid:39) (3 . 8 − 10 . 2), for a total of seven detected lines out of eight, corresponding to a success rate of 88%. This allows us to double the number of distant BCGs with clear detections in at least two CO lines. We then combined our observations with available stellar, star formation, and dust properties of the BCGs and compared them with a sample of ∼ 100 distant cluster galaxies with observations in CO. Our analysis yields large molecular gas reservoirs M H 2 (cid:39) (0 . 5 − 1 . 4) × 10 11 M (cid:12) , high excitation ratios r 31 = L (cid:48) CO(3 → 2) / L (cid:48) CO(1 → 0) (cid:39) (0 . 1 − 0 . 3), long depletion times τ dep (cid:39) (2 − 4) Gyr, and high M H 2 / M dust (cid:39) (170 − 300) for the three targeted BCGs. The excitation ratio r 31 of intermediate-z BCGs, including RX1532 and M1932 from previous studies, appears to be well correlated with the star formation rate and e ﬃ ciency, which suggests that excited gas is found only in highly star-forming and cool-core BCGs. By performing color-magnitude plots and a red-sequence modeling, we ﬁnd that recent bursts of star formation are needed to explain the fact that the BCGs are measurably bluer than photometrically selected cluster members. To explain the global observed phenomenology, we suggest that a substantial amount of the molecular gas has been accreted by the KiDS BCGs but still not e ﬃ ciently converted into stars. KiDS 1220 also shows a double-horn emission in CO(3 → 2), which implies a low gas concentration. The modeling of the spectrum yields an extended molecular gas reservoir of ∼ 9 kpc, which is reminiscent of the mature extended-disk phase observed in some local BCGs.


Introduction
Brightest cluster galaxies (BCGs) are among the most massive galaxies in the Universe, are typically massive ellipticals of cD type, and are often radio galaxies (Zirbel 1996). There exists a strong, albeit still debated, coevolution between BCGs, their cluster environments (Lauer et al. 2014), and the accretion processes onto the supermassive black holes residing at galaxy centers (Fabian 2012).
Because of their exceptional masses and luminosities, BCGs are believed to evolve via phenomena such as dynamical friction (White 1976), galactic cannibalism (Hausman & Ostriker 1978), and interactions with the intracluster medium (ICM; Stott et al. 2012). Active galactic nucleus (AGN) feedback-regulated cool-Spectra in Fig. 3 are available in electronic form at the CDS via anonymous ftp to cdsarc.u-strasbg.fr (130.79.128.5) or via http://cdsweb.u-strasbg.fr/cgi-bin/qcat?J/A+A/ e-mail: gianluca.castignani@unibo.it ing in the ICM is another mechanism commonly invoked to explain the evolution of BCGs and in particular the fact that large reservoirs of molecular gas are rarely observed in BCGs (e.g., Edge 2001;Salomé & Combes 2003). While condensation of the ICM occurs in the central and low-entropy regions of the cluster, the energy injected into the ICM via radio jets propagating from the BCGs themselves may prevent the formation of the strong (∼ (100 − 1000) M yr −1 ) cooling flows that were theoretically predicted (e.g., Fabian 1994;McNamara et al. 2000;Peterson et al. 2003;Peterson & Fabian 2006;McNamara & Nulsen 2012).
Most of the BCG stellar mass is assembled at high redshifts (z ∼ 3 − 5) from smaller sources that are later swallowed by the BCG (De Lucia & Blaizot 2007). Indeed, BCGs have doubled their stellar mass since z ∼ 1 (Lidman et al. 2012), which is consistent with evolution via the dry accretion of satellites (Collins et al. 2009;Lavoie et al. 2016). However, recent studies have proposed that high star formation is present in distant BCGs (McDonald et al. 2016;Webb et al. 2015a,b;Fig. 1: Color composite 40 × 40 KiDS images (red = i band; green = r band; blue = g band) with OmegaCAM at VST centered around the BCGs KiDS 0920 (left), KiDS 1220 (center), and KiDS 1444 (right). North is up, and east is to the left. Bonaventura et al. 2017), a scenario where star formation is fed by rapid gas deposition, possibly regulated by cooling flows (Salomé & Combes 2003;Olivares et al. 2019). Such results imply that at least some distant BCGs should experience strong environment-driven gas processing and quenching mechanisms, such as strangulation, ram pressure stripping, and galaxy harassment (e.g., Larson et al. 1980;Moore et al. 1999), which are still, however, substantially unexplored for distant BCGs.
In support of the latter scenario, some studies report growing evidence for large reservoirs of molecular gas ( 10 10 M ) feeding the star formation in some intermediate-to high-redshift BCGs (McDonald et al. 2014;Russell et al. 2017;Castignani et al. 2020a;D'Amato et al. 2020;Dunne et al. 2021;O'Sullivan et al. 2021). However, a number of other distant star-forming BCGs have only robust upper limits on the order of ∼ 10 10 M in their H 2 gas masses (Castignani et al. 2020a,b,c). In Castignani et al. (2020b), in particular, we found gas-rich companions around a distant BCG at z = 1.7, which is instead gas-poor.
These results suggest that gas-rich BCGs in the distant Universe (z ∼ 0.3 − 2) belong to a rare population of star-forming cluster ellipticals that experience gas deposition and a subsequent conversion of gas into star formation. As suggested in Castignani et al. (2020a) for distant BCGs of the Cluster Lensing And Supernova survey with Hubble (CLASH) sample (Postman et al. 2012), this is an environment-driven process. The detection of large molecular gas reservoirs is likely due to the strong cooling of the ICM at the cluster cores, which ultimately favors condensation into molecular gas, similar to what has been observed in some local BCGs (Salomé et al. 2006;Tremblay et al. 2016).
To understand the effect of the cluster environments on the mass assembly and star formation of distant BCGs, it is essential to start building a sample of star-forming BCGs detected at different molecular gas transitions, which probe different gas densities, and possibly different degrees of environmental processing. However, this is challenging as it requires a large sample of well-defined clusters and a robust characterization of their BCGs at multiple wavelengths to enable the stellar and star formation properties to be accurately inferred.
To this aim, in this work we report the study of three starforming BCGs that we selected from the multiwavelength widefield Kilo Degree Survey (KiDS; de Jong et al. 2017;Kuijken et al. 2019) and that we observed in the first three CO transitions with the Institut de Radioastronomie Millimétrique (IRAM) 30m telescope. The present study is part of a larger campaign targeting galaxies in and around distant clusters in CO, mainly BCGs (Castignani et al. 2018(Castignani et al. , 2019(Castignani et al. , 2020a(Castignani et al. ,b,c,d, 2022a, with the final goal of evaluating the impact of dense megaparsec-scale environments in processing the galaxies' gas reservoirs. The paper is structured as follows. In Sect. 2 we describe the BCGs that are the subject of this work, in Sect. 3 we describe the molecular gas observations and data analysis, in Sect. 4 we describe our comparison samples, in Sect. 5 we present the results, and in Sect. 6 we summarize them and draw our conclusions. Stellar mass and star formation rate (SFR) estimates reported in this work for the three targeted BCGs rely on the Chabrier (2003) initial mass function (IMF). Magnitudes are reported in the AB system. Throughout this work we adopt a flat Λ-cold dark matter (ΛCDM) cosmology with matter density Ω m = 0.30, dark energy density Ω Λ = 0.70, and Hubble constant h = H 0 /(100 km s −1 Mpc −1 ) = 0.70.

The selection of the BCGs in KiDS
KiDS is a wide field imaging survey that in its final release will cover around 1,350 deg 2 of the sky (de Jong et al. 2017;Kuijken et al. 2019), equally divided between an equatorial field (KiDS-N) and a southern field (KiDS-S). VLT Survey Telescope (VST) ugri optical observations and Visible and Infrared Survey Telescope for Astronomy (VISTA) telescope ZYJHK infrared imaging enable the detection of ∼ 5 × 10 7 galaxies down to a limiting AB magnitude r ∼ 25.0 (5σ in 2 aperture). Optical-infrared photometry enabled the determination of accurate photometric redshifts up to z ∼ 1 (Hildebrandt et al. 2021;Benítez 2000), which are complemented by spectroscopic information as KiDS is also covered by Galaxy And Mass Assembly (GAMA; Baldry et al. 2010) and Sloan Digital Sky Survey (SDSS; York et al. 2000) spectroscopic surveys.
This rich data set allows the detection of distant clusters in KiDS and the characterization of their galaxy populations. Distant clusters up to z ∼ 1 have been indeed detected at signal-tonoise ratio S/N> 3 within KiDS (Maturi et   Objects (AMICO) cluster finder (Bellagamba et al. 2018), which looks for overdensities of galaxies using photometric redshifts and a matched filtering. For this work we consider the parent sample of 5,251 BCGs (of which 1,484 have spectroscopic redshifts) in the equatorial KIDS-N area, which can be targeted with observational facilities in the northern hemisphere. These BCGs were selected by Radovich et al. (2020) within the Data Release (DR) 3 of KiDS (de Jong et al. 2017) as the brightest galaxies in the r band within the cluster cores, on the basis of AMICO probabilistic cluster memberships.
We aim at targeting actively star-forming BCGs, to investigate their molecular gas reservoirs feeding the star formation, and explore their properties in connection with those of the host galaxies and their environments. We thus limit ourselves to the 684 BCGs in DR3 KiDS-N with spectroscopic redshifts z > 0.3. This redshift cut allows us to focus on distant BCGs, and on their evolution, while BCGs at lower redshifts, including those of Abell clusters, have been extensively studied in previous work (Edge 2001;Salomé & Combes 2003;Olivares et al. 2019).
Similarly to the target selection we did in Castignani et al. (2019) we further limit ourselves to those sources with clear counterparts in the all-sky data release of the Wide-field Infrared Survey Explorer (WISE; Wright et al. 2010) at 22 µm in the observer frame (i.e., the W4 channel). These are thus good starforming BCG candidates in the distant Universe. We found only nine such BCGs with WISE W4 counterparts at S/N> 3, of which six had S/N> 3.3. Among the latter in this work we consider the three with highest values in the range S/N∼ 4.1 − 6.5, while the remaining three have WISE W4 S/N∼ 3.1 − 3.3 and will be the subject of a forthcoming paper.
The three intermediate-z BCGs of this work are KiDS 0920 (i.e., G198598 at z = 0.3216), KiDS 1220 (G172934 at z = 0.3886), and KiDS 1444 (G363753 at z = 0.4417). In parentheses, we report the ID numbers of the sources from GAMA DR3 (Baldry et al. 2018) and the corresponding spectroscopic redshifts. 1 The three sources belong to moderately rich clusters, which have been detected with S/N = 3.7, 4.7, and 3.7 (Maturi et al 2019; Bellagamba et al. 2019), while they have intrinsic richness λ = 20 ± 3, 33 ± 6, and 19 ± 3, respectively. The corresponding richness-based cluster masses range between log(M 200 /M ) = 13.8−14.2 (see Table 1). Here λ is the intrinsic cluster richness as defined in Maturi et al (2019). Based on these properties we infer, on average, a small probability of ∼ 0.3% of the three clusters being false positives (see Fig. 12 of Maturi et al 2019). We additionally searched for our clusters in published catalogs of clusters. KiDS 1220 matches with WHL J121959.1-022611, which is an optically identified cluster in SDSS, found by both the redMaPPer (Rykoff et al. 2014;Rozo et al. 2015) and the Wen et al. (2012) cluster finding algorithms. On the other hand, our search did not produce positive results when looking in catalogs of X-rays and Sunyaev-Zel'dovich clusters. This is not surprising as our considered clusters have relatively low masses and are at intermediate redshifts. Figure 1 displays the color composite images for the three KiDS BCGs. The images show that our targets are massive elliptical galaxies, bulge-dominated, but they also show complex and disturbed morphology. Indeed, KiDS 0920 shows an elongated reddish emission, possibly associated with the stellar disk, along the northeast to southwest direction. Similarly, KiDS 1444 shows disturbed morphology, with evidence of a faint elongated tail to the south. Both KiDS 1220 and KiDS 1444 have nearby companions at a projected separation of ∼ 3 arcsec (i.e., ∼17 kpc at the redshifts of the sources).

Stellar and star formation properties
We now investigate the multiwavelength properties of the three BCGs of this work.  Driver et al. (2018). Photometric data include GALEX (Martin et al. 2005;Morrissey et al. 2007) in the ultraviolet, SDSS (York et al. 2000) in optical, as well as VISTA Kilo-degree Infrared Galaxy Survey (VIKING; Edge et al. 2013), WISE (Wright et al. 2010), and Herschel-ATLAS (Eales et al. 2010;Valiante et al. 2016) in the near-to far-infrared. Table 1 summarizes the BCG properties including their stellar, star formation, and dust properties from the SED fits. All three BCGs have stellar masses exceeding 10 11 M , which confirms that they are very massive galaxies. Furthermore, with the only exception of KiDS 1444, the remaining two BCGs have farinfrared emission from Herschel, well modeled by a dust component. This supports the scenario that the observed infrared emission is mainly due to star formation, with minimum AGN contamination. Indeed, for the three BCGs, dust masses and luminosities are in the range log(M dust /M ) 8.5 − 8.7 and log(L dust /L ) 11.5 − 11.8, respectively, typical of luminous infrared galaxies (LIRGs). The SFR estimates based on the ultraviolet to far-infrared SED modeling are also high, and in the range ∼ (20 − 30) M /yr. These are consistent with, but higher, than the SFRs of main sequence (MS) galaxies of similar mass and redshifts (Speagle et al. 2014); details are provided in Table 1.

Spectral energy distributions
We remind the reader that we selected the three distant starforming BCGs from a parent sample of 684 BCGs at z > 0.3 (Sect. 2.1). They represent 0.4% of the parent sample and are among the most star-forming BCGs in KiDS. This is a rare population of BCGs that experience stellar mass assembly and are thus caught in a special phase of their evolution (see, e.g., Castignani et al. 2020a,b,

Line diagnostics
We further investigate the star formation properties of the BCGs using the available spectroscopy in GAMA DR3. According to Kennicutt (1998), the SFR is directly proportional to both Hα line and [O II] forbidden-line doublet (3726 Å, 3729 Å) luminosities L Hα and L [O II] , respectively. We thus use these relations to estimate the SFR, calibrated to a Chabrier (2003) IMF. We also correct for dust attenuation using the empirical relation by Garn & Best (2010), which relates dust attenuation to stellar mass as follows: where M = log(M /M ) − 10. With these conventions, the Hαand [O II]-based SFRs can be expressed as follows (see, e.g., Gilbank et al. 2010;Zeimann et al. 2013;Old et al. 2020): where r lines is the ratio of extinguished [O II] to Hα fluxes. We found GAMA DR3 spectra for all three BCGs, taken with the AAOmega/2dF multi-fiber spectrograph on the 3.9 m Anglo-Australian Telescope. Gaussian fits to several lines are also available (Gordon et al. 2017). Among the three BCGs, Hα emission is detected for KiDS 0920 only, with a flux of F Hα = (208.6 ± 18.5) × 10 −17 erg s −1 cm −2 . At a spectral resolution of R 1300 the Hα line is not blended with [N II]. For the other two BCGs Hα is redshifted to wavelengths longer than 9100 Å, thus outside the wavelength range of 3750 Å λ 8850 Å covered by the spectrograph (Hopkins et al. 2013 .3), (72.6 ± 15.8), and (88.0 ± 20.1), respectively, in units of 10 −17 erg s −1 cm −2 . The line ratio r lines is equal to r lines = 0.9 for KiDS 0920, while we fix it to r lines = 0.5 for the other two BCGs for which we do not have Hα measurements, as in Gilbank et al. (2010).
Using these values, Eq. 2 yields SFR Hα = (20.9 ± 1.9) M /yr for KiDS 0920. We obtain a similar SFR Hα = (27.6 ± 2.6) M /yr if we compute the dust attenuation using the observed Hα/Hβ Balmer decrement, following Domínguez et al. (2013). Eq. 3 yields instead SFR [O II] = (24.1 ± 5.3) and (37.7 ± 8.6) M /yr for KiDS 1220 and 1444, respectively. 2 We find a good agreement, within the reported uncertainties, when comparing these line-based SFRs with those in Table 1 that are used throughout this work and come from the SED fits. This is remarkable, considering that discrepancies, even up to a factor of ∼ 3, are not uncommon when comparing different estimates of the SFR (e.g., Calzetti 2013).
Finally, we use the line ratios to identify the dominant mechanism of gas ionization using ratios of strong optical lines. We do it by means of the Baldwin, Phillips & Terlevich (BPT) diagram (Baldwin et al. 1981;Veilleux & Osterbrock 1987 Sarzi et al. 2010;Singh et al. 2013). These results are not surprising, as BCGs are often associated with AGN and radio galaxies in particular (Zirbel 1996). In this case the BCG is likely a Type 2 AGN, as the spectrum does not show any clear broad emission line (e.g., Hα and Hβ).
For the other two BCGs, KiDS 1220 and KiDS 1444, since we do not have Hα we cannot locate them precisely in the BPT diagram. However, the sources have log([O III] λ5007/Hβ) = −0.35 ± 0.19 and −0.15 ± 0.18, respectively, which are lower than the corresponding line ratio of KiDS 0920. This suggests that the two sources are more likely located in the star-forming region of the BPT diagram.
Last, we inspected the position of our sources in the infrared color-color WISE diagram (Jarrett et al. 2017). Interestingly, KiDS 0920 is classified as an intermediate disk, while KiDS 1220 and KiDS 1444 are starbursts. These findings strengthen the selection of our three targets as star-forming galaxies and suggest that any possible contamination from a circum-nuclear dusty torus at infrared wavelengths is likely negligible, which is also supported by the good agreement between the SED-and line-based SFRs.

IRAM 30m observations and data reduction
We observed the three KiDS BCGs using the IRAM 30m telescope at Pico Veleta in Spain. The observations were carried out in June, August, and October 2021 (ID: 079-21; PI: G. Castignani). We used the Eight Mixer Receiver (EMIR) to observe CO(J→J-1) emission lines from the BCGs, where J is a positive integer denoting the total angular momentum. For each galaxy, the specific CO(J→J-1) transitions were chosen to maximize the likelihood of the detection, in terms of the ratio of the predicted signal to the expected rms noise.
The E090, E150, and E230 receivers, operating between ∼1-3 mm, offer 4×4 GHz instantaneous bandwidth covered by the correlators. For KiDS 1220 and KiDS 1444 BCGs we used both E090 and E150 receivers, simultaneously, to observe CO(1→0) and CO(2→1), redshifted to 3 mm and 2 mm in the observer frame, respectively. For all three targets we also used both E090 and E230 receivers to target simultaneously CO(1→0) and CO(3→2) lines, redshifted to 3 mm and 1 mm, respectively. We note that we could not observe KiDS 0920 in CO(2→1) as the line is redshifted to 174.440 GHz in the observer frame, too close to the 183 GHz water vapor line, and thus prohibited because of low atmospheric transmission.
The IRAM 30m half power beam width (HPBW) is ∼ 16 arcsec λ obs 2 mm (Kramer et al. 2013), where λ obs is the observer frame wavelength. All targets are thus unresolved by our observations. We refer to Sect. 5.1 for further discussion. The wobblerswitching mode was used for all the observations to minimize the impact of atmospheric variability. The Wideband Line Multiple Autocorrelator was used to cover the 4×4 GHz bandwidth, in each linear polarization. We also simultaneously recorded the data with the fast Fourier transform spectrometers, as a backup, at 200 kHz resolution.
During our June and August observations, we had variable weather conditions, while in October we had good weather during our observations. Overall we had average system temperatures typical of summer season and equal to T sys 120 − 130 K, 190-300 K, and 290-440 K at 3, 2, and 1 mm, respectively. For pointing and focus we used the following strong sources: Mars, Mercury, Venus, 3C 084, 3C 273, 3C 279, and OJ +287.
We summarize our observations in Table 2. The data reduction and analysis were performed using the CLASS software of the GILDAS package 3 . For each source and frequency tuning, we merged all available observations. We flagged only a few bad scans. At a velocity resolution of 60 km s −1 the resulting rms values of the antenna temperature (Ta * ) are similar among our three targets, at their corresponding frequency. They range between rms = 0.20 − 0.31 mK at 1 mm and between rms = 0.57 − 0.59 mK at 2 mm. At 3 mm we have instead rms = 0.59, 0.38, and 0.96 mK for KiDS 0920, 1220, and 1444, respectively. Our results are presented in Sect. 5, where we further describe the analysis of our IRAM 30m spectra.

The comparison samples
In this section we describe samples of galaxies with available observations of their molecular gas that we use as a comparison for our analysis.

Distant cluster galaxies observed in CO
We want to put our galaxies in a broader context, in comparison with existing samples of distant cluster galaxies. To this aim, we consider the compilation built by Castignani et al. (2020b), comprising 120 (proto-)cluster galaxies at 0.2 z 5.0 with M > 10 9 M and CO observations from the literature. In particular, we refer to Tables A.1, A.2 of Castignani et al. (2020b), where stellar, star formation, and gas properties are listed. Since in this work we consider massive intermediate-redshift BCGs, we limit ourselves to sources in the compilation with stellar masses M > 10 10 M and redshifts z < 2 (that is, we exclude the most distant clusters and still assembling protoclusters). As we have detected all three of our KiDS BCGs in CO, concerning the comparison sample we similarly considered sources with CO detections, that is, we excluded sources with only upper limits to the H 2 masses. The only exception is represented by the SpARCS 104922.6+564032.5 BCG at z = 1.7, further discussed below, which is one of the few BCGs at z > 1 with observations in CO, as well as M and SFR estimates. This selection yields 58 galaxies over 16 clusters at moderate to high redshifts: z ∼ 0.2−0.5 (Geach et al. 2011;Jablonka et al. 2013;Cybulski et al. 2016); z ∼ 1 − 1.2 (Wagg et al. 2012;Castignani et al. 2018); z ∼ 1.5 − 2.0 (Aravena et al. 2012;Rudnick et al. 2017;Noble et al. 2017Noble et al. , 2019Hayashi et al. 2018;Kneissl et al. 2019;Coogan et al. 2018;Castignani et al. 2020c).
For our comparison, as further discussed below, we include additional 24 BCGs drawn from the CLASH and SpARCS surveys. They span the redshift range z ∼ 0.2 − 1.7 and have been observed in CO by recent studies (Fogarty et al. 2019;Castignani et al. 2020a,c;Dunne et al. 2021). With the inclusion of  Table 1 (bottom x axis) and the observer frame frequency (top x axis). these sources and the three BCGs presented in this work the full list comprises 84 cluster galaxies over 42 clusters at z 0.2 − 2 and stellar masses in the range log(M /M ) 10 − 12. We use this sample in Sect. 5.2 to perform a comparison with the three BCGs of this work in terms of stellar mass, SFR, molecular gas content, and depletion time.

Distant BCGs observed in CO
Among the distant cluster galaxies with CO observations and stellar mass estimates in the comparison sample, considered above, there are 24 BCGs, as outlined in the following.
The SpARCS 104922.6+564032.5 BCG, located at z = 1.7. Recent CO(1→0) Jansky Very Large Array (JVLA) observations by Barfety et al. (2022) yield a large molecular gas reservoir of M H 2 ∼ 10 11 M , in agreement with the value previously found with single-dish observations by Webb et al. (2017). However, the CO(1→0) peak is offset by ∼ 15 kpc to the southeast of the BCG, which suggests that the molecular gas reservoir may be associated with cluster core BCG companions. NOrthern Extended Millimeter Array (NOEMA) observations (Castignani et al. 2020c) supported this scenario, with the detection in CO(4→3) of two gas-rich BCG companions within 20 kpc from the BCG. The observations allowed us to set a robust upper limit to the molecular gas reservoir M H 2 < 1.1 × 10 10 M for the BCG that we use in this work.
There are 17 BCGs from the recent SpARCS survey, detected in CO(2→1) and S/N 3 by Dunne et al. (2021). They span the redshift range z ∼ 0.2 − 1.2 and have molecular gas masses in the range M H 2 ∼ (10 10 − 10 11 ) M .

Galaxies with CO(1→0) and CO(3→2) observations
In addition to cluster galaxies with observations in CO, we considered a second compilation of sources with molecular gas observations in the first and third CO transitions. We stress that the number of such sources is limited in the literature as detecting galaxies in both these transitions is expensive in terms of observational resources. We use these galaxies in Sect. 5.3 to investigate the relation between the excitation ratios and both star formation and gas properties, in comparison with the KiDS BCGs. We are particularly interested in CO(1→0) and CO(3→2) as we observed all three BCGs of this work at these transitions. The galaxies considered for the comparison, outlined below, are massive, log(M /M ) ≥ 10, similarly to our BCGs. They also cover broad ranges in redshift and SFR, which allow us to effectively search for possible relations and trends.
We additionally included five 2.5 < z < 2.7 (ultra)LIRGs with infrared luminosities in the range ∼ (11.6 − 12.    Tacconi et al. (2018). The horizontal dashed lines correspond to y-axis values equal to unity, while the dotted lines denote the fiducial uncertainties associated with the MS. The uncertainty is chosen to be equal to ±0.48 dex for both the left and right panels since the MS is commonly identified by 1/3 < SFR/SFR MS < 3. For the central panel, the plotted uncertainties are estimated at redshift z = 0.5. The color-coding for the data points is the same as in Fig. 4.

Molecular gas properties
With our observations we clearly detect all three targeted BCGs in multiple transitions, with a signal-to-noise ratio (S/N) in the range S/N (3.8 − 10.2), as outlined in the following. Our observations thus increase the still limited sample of rare distant starforming BCGs with detections of their molecular gas. Two out of the three targeted BCGs, namely KiDS 1220 and KiDS 1444 were detected in CO(1→0), CO(2→1), and CO(3→2). The third one, KiDS 0920, was detected in CO(1→0), while we set an upper limit to the CO(3→2) emission. The corresponding spectra are shown in Fig. 3.
These results imply that we double the still limited sample of distant (z > 0.2) BCGs with clear detections in two distinct CO lines. To the best of our knowledge there are in fact only other two distant BCGs detected in multiple CO transitions. These are the intermediate-redshift (z ∼ 0.4) BCGs RX1532 and M1932 from CLASH (Postman et al. 2012). RX1532 was detected both in CO(1→0) and CO(2→1) with our recent CLASH IRAM 30m campaign (Castignani et al. 2020a). M1932 was instead detected in CO(1→0), CO(3→2), and CO(4→3) with ALMA by Fogarty et al. (2019).
Based on the optical morphologies illustrated in Fig. 1 the three targeted KiDS BCGs are unresolved by our observations, assuming a CO-to-optical size ratio ∼ 0.5 (Young et al. 1995) and given that the IRAM 30m beam is ∼ 16 arcsec at 2 mm (Kramer et al. 2013, see our Sect. 3). We detect all targeted lines except CO(3→2) of KiDS 0920, this yields an exceptional success rate of 88%, that is, seven detections out of eight targeted lines, higher than that of recent surveys of distant BCGs (Castignani et al. 2020a,b;Dunne et al. 2021). We then fit each detected CO line with a Gaussian curve, after removing the baseline with a linear fit. We used instead a combination of two Gaussian curves to fit the CO(3→2) line of KiDS 1220, as the spectrum shows a clear double-horn profile. Hints of a double-peaked line are also present in the CO(1→0) and CO(2→1) spectra.
We believe that the observed double-peaked CO emission is originated from the KiDS 1220 BCG only and we further discuss the underlying structural parameters in Sect. 5.5. A possible alternative is that the nearby companion of KiDS 1220 BCG contributes to the observed CO(3→2) emission. However, this is less likely as the source is at a projected separation of 2.7 arcsec, so that with a HPBW of 9.6 arcsec (Kramer et al. 2013) the efficiency is suppressed down to 80% the maximum. With a similar argument, we believe that the CO emission observed for KiDS 1444 comes from the BCG only, with null or negligible contamination from the nearby companion within the beam. However, it is possible that the KiDS 1220 and KiDS 1444 BCGs are interacting with their companions, which may favor the replenishment of the BCGs' gas reservoirs, thus ultimately contributing to sustain the high level of star formation activity observed in the targeted BCGs.
We detect KiDS 0920 in CO(1→0), and it is likely that the detected molecular gas comes from the elongated disk-like structure seen in Fig. 1. However, we do not find any evidence of CO(3→2) emission, so that we first removed the baseline in the CO(3→2) spectrum with a linear fit. We then estimated the rms noise level for the antenna temperature (Ta * ). Using standard conversions, outlined below, we converted the rms into a 3σ upper limit to the integrated CO flux, at a resolution of 300 km/s, which corresponds to the typical full width at half maximum for massive galaxies and BCGs in particular (Edge 2001;Dunne et al. 2021).
In Table 2, we report the results of our observations, where we apply standard efficiency corrections to convert i) Ta * into the main beam temperature T mb and then ii) T mb into the corresponding CO line flux, with the conversion factor of 5 Jy/K. In particular, we adopt the following efficiency corrections T mb /T a * = 1.2, 1.4 and 2.0 for ∼3, 2, and 1 mm observations, respectively. 5 We then derive the CO(J→J-1) luminosity L CO(J→J−1) , in units K km s −1 pc 2 , from the velocity integrated CO(J→J-1) flux 5 https://www.iram.es/IRAMES/mainWiki/Iram30mEfficiencies S CO(J→J−1) ∆ , in units Jy km s −1 . To this aim we used Eq. (3) from Solomon & Vanden Bout (2005): where ν obs is the observer frame frequency, in GHz, of the CO(J→J-1) transition, D L is the luminosity distance in Mpc, and z is the redshift of the BCG. All three BCGs are detected in CO(1→0): this enables us to derive the total H 2 gas mass as M H 2 = α CO L CO(1→0) . The use of the CO(1→0) transition is advantageous with respect to higher-J ones as it does not require any assumption on the excitation ratio r J1 = L CO(J→J−1) /L CO(1→0) .
The three KiDS BCGs have similar specific star-formation rates, sSFR 1.6 × sSFR MS (see Table 1). They are all less than the value of 3 × sSFR MS , below which the galaxy sSFR is within the typical scatter of the MS. To estimate H 2 gas masses we therefore assume a Galactic CO-to-H 2 conversion factor of α CO = 4.36 M (K km s −1 pc 2 ) −1 , commonly adopted for starforming galaxies at the MS. The use of a single α CO factor also allows us to do a homogeneous comparison in terms of gas content.
As all three BCGs are detected in CO(1→0), we use both LO and LI sidebands to search for the CN(N=1→0) doublet, which is redshifted to ∼(79-86) GHz in the observer frame, close to the CO(1→0) line. However, none of the three BCGs shows evidence for any of the CN(N=1→0) J=1/2→1/2 and CN(N=1→0) J=3/2→1/2 lines of the doublet. This is not surprising, since they are usually much fainter than CO(1→0). In our recent study we performed a similar search for the CN(N=1→0) doublet in the BCG RX1532 at z ∼ 0.4, which also led to a non-detection (Castignani et al. 2020a). However, thanks to the long integration times of several hours at 3 mm (Table 2) we reach low rms values of 1.7, 1.3, and 2.3 mJy for the KiDS 0920, 1220, and 1444 BCGs, respectively, at a resolution of 60 km/s. We do not attempt to set any upper limit to the continuum emission of the BCGs by using the available total 15 GHz (LI, LO, UI, UO) bandwidths, for each polarization. In fact, the BCGs are faint and there is significant intrinsic atmospheric instability at millimeter wavelengths that prevents us from robustly determining the continuum levels or estimating upper limits.
We also estimate the depletion timescale τ dep = M H 2 /SFR, which is the characteristic time to deplete the molecular gas reservoirs. To compute τ dep we use the SED-based SFR estimates reported in Table 1. Similarly, we also estimate the molecular-gas-to-stellar-mass ratio, M H 2 /M , and the ratio between the molecular gas and the dust masses. For a comparison we compute the depletion time τ dep,MS and the molecular gas to stellar mass ratio M H 2 M MS for MS galaxies in the field with redshift and stellar mass equal to those of our BCGs, by using empirical prescriptions by Tacconi et al. (2018), calibrated to the CO-to-H 2 conversion factor of α CO = 4.36 M (K km s −1 pc 2 ) −1 used in this work.
In Table 3 we summarize the molecular gas properties for the three KiDS BCGs. As both KiDS 1220 and KiDS 1444 are detected in multiple CO transitions, in the Table we also list their r 21 and r 31 excitation ratios, while for KiDS 0920 we report the upper limit to r 31 .

Gas fraction, star formation, and depletion time
We now exploit the compilation of 84 cluster galaxies with available H 2 gas masses, SFR, and M estimates presented in Sect. 4.1. Figure 4 displays the ratio of molecular gas to stellar mass M H 2 /M as a function of both redshift and specific SFR for the compilation of cluster galaxies with observations in CO and stellar mass estimates.
In Fig. 5 we show instead the SFR, M H 2 /M , and the depletion time (τ dep ), all normalized to their MS values, as a function of the stellar mass. In both figures, to allow a homogeneous comparison, we have used α CO = 4.36 M (K km s −1 pc 2 ) −1 for all data points and the overplotted Tacconi et al. (2018) relations. We also highlight the KiDS BCGs (red dots) and BCGs from the literature, while the remaining cluster galaxies in our comparison sample are shown in background (gray dots). By looking at Fig. 5 it is striking that our BCGs have high stellar masses, log(M /M ) ∼ 11.4, that are systematically higher than the SpARCS BCGs in Dunne et al. (2021). This result strengthens the BCG association for our targets, further discussed in Sect. 5.4.
Interestingly, the figures also show that all three KiDS BCGs are gas-rich, with H 2 gas masses largely exceeding 10 10 M and H 2 -to-stellar mass ratios in the range ∼ 0.2 − 0.6, which are values similar to those found by previous studies of star-forming BCGs at intermediate redshifts (SpARCS and CLASH survey;Fogarty et al. 2019;Castignani et al. 2020a;Dunne et al. 2021), highlighted in the Figures.
Furthermore, while our KiDS BCGs have M H 2 /M ratios higher than those of MS galaxies, they also have MS levels of star formation activity (Fig. 5 left). These results imply that our targeted BCGs convert the molecular gas into stars quite inefficiently, having depletion timescales τ dep = M H 2 /SFR in the range ∼ (2 − 4) Gyr. As illustrated in Fig. 5 (right) these values are still consistent with the depletion timescales of MS galaxies, given the uncertainties, but higher than those found for starforming CLASH and SpARCS BCGs. This result is intriguing and may suggest that the gas fueling in the three KiDS BCGs differs somehow from that of previously studied star-forming BCGs at intermediate redshifts. Similarly, as outlined in Table 3, our three KiDS BCGs have H 2 -to-dust mass ratios in the range ∼ 170 − 300, thus exceeding the ratio of ∼ 100, typical of dis-tant star-forming galaxies (Berta et al. 2016;Scoville et al. 2014Scoville et al. , 2016. To explain both the high τ dep and the high M H 2 /M dust , a possibility is that a substantial amount of the H 2 gas has been recently accreted, but still not efficiently converted into stars. This accretion may have occurred via filaments (as it may be the case of KiDS 0920, which shows an elongated morphology), or via interaction with the companions (for KiDS 1220 and 1444). Alternatively, adopting a lower α CO conversion factor than the Galactic one, as usually done for starburts, would lead to lower M H 2 and τ dep . Motivated by these results, in the following Section we further investigate the differences of the selected BCGs with respect those studied in previous works by means of the excitation ratio.

Excitation ratios
To the best of our knowledge there exist only a few distant BCGs with CO detections in different transitions, as outlined in Sect. 5.1. These are RX1532 (z = 0.361, Castignani et al. 2020a) and M1932 (z = 0.353, Fogarty et al. 2019), in addition to the KiDS BCGs of this work. As the sources are all detected in CO(1→0) as well as in CO(2→1), CO(3→2), or CO(4→3) we consider the r 21 , r 32 , and r 43 excitation ratios, whenever available.
As pointed out above, all these intermediate-z BCGs are starforming and span a quite broad range in the SFR (20 − 200) M /yr. On the other hand the BCGs are gas-rich, with similar gas content in terms of M H 2 (0.5 − 1.4) × 10 11 M and M H 2 /M (0.1 − 0.6), which are less scattered than the SFR. Consistently with this result, we do not find any trend when plotting the excitation ratios r J1 against the M H 2 or M H 2 /M . Figure 6 displays the excitation ratios r J1 of the BCGs against the sSFR (left) and the depletion time (right). For the latter we used α CO = 4.36 M (K km s −1 pc 2 ) −1 for all sources in the plot, to allow a homogeneous comparison. The plots show a clear dichotomy between the CLASH BCGs (M1932 and RX1532) that show the highest SFR 100 M /yr and the KiDS BCGs consid- Fig. 7: Color magnitude plots for the sources in the field of the BCGs KiDS 0920 (top), KiDS 1220 (center), and KiDS 1444 (bottom). Field galaxies (gray points), photometrically (blue points) and spectroscopically (red points) selected cluster members, and the BCGs (green stars) are highlighted. Red-sequence models are shown as horizontal lines. See the legend and text for further details. ered in this work, with lower SFR (20 − 30) M /yr. The latter not only have higher depletion times, as discussed in Sect. 5.2, but also lower excitation ratios.
As all three BCGs have observations in both CO(1→0) and CO(3→2) we include in Fig. 6 the ASPECS, BASS, and xCOLD GASS sources presented in Sect. 4.3, for a comparison. When considering the r 31 values of the KiDS BCGs and the comparison galaxies altogether we find that the excitation ratio turns out to be correlated with sSFR and anticorrelated with the depletion time, at a level of 3.1σ (p − value = 1.81 × 10 −3 ) and 4.7σ (p − value = 2.73 × 10 −6 ), as found with the Spearman test, respectively. We also note that, at fixed BCG, it holds r 41 < r 32 < r 21 as indeed higher-J CO transitions are associated with denser gas, more difficult to excite.
These findings suggest that highly excited gas is preferentially found in highly star-forming and cool-core BCGs such as RX1532 and M1932, as well as in star-forming galaxies of the comparison samples (e.g., ASPECS). As further discussed in Castignani et al. (2020a) the cool-core cluster environments likely favor the condensation and the inflow of gas onto the BCGs themselves. The present work further supports this scenario, while suggesting meanwhile that SFRs lower than ∼ 100 M /yr appear to be not sufficient to excite the gas to the highest levels. The latter is true for the KiDS BCGs, as well as for lower mass xCOLD GASS and BASS galaxies of the comparison sample. Interestingly, trends of the excitation ratios as a function of star formation and gas properties were also reported in previous studies of local galaxies (Lamperti et al. 2020;Yao et al. 2003;Leech et al. 2010), while in this work we find that they are also valid for star-forming and distant BCGs.
Although the excitation ratios of the BCGs nicely follow the expected trends as a function of the sSFR and the depletion time, we note that the r J1 values of the KiDS BCGs are lower than those of distant star-forming galaxies, for which r 31 0.6, on average (Carilli & Walter 2013;Bothwell et al. 2013;Genzel et al. 2015;Boogaard et al. 2020). We remind that the IRAM 30m HPBW is ∼24, 16, and 10 arcsec for the CO(1→0), CO(2→1), and CO(3→2) transitions, respectively. As all KiDS BCGs have nearby companions within 10 arcsec in projection, the latter ones might contribute to the observed CO(1→0) and CO(2→1) emission, while the contamination is likely negligible for CO(3→2), as already discussed in Sect. 5.1. The relatively low r J1 values of the KiDS BCGs might therefore be explained if the BCG companions are rich in gas. High resolution interferometric observations are however needed to have a definitive answer, as in the case of SpARCS 104922.6+564032.5 BCG at z = 1.7 (Castignani et al. 2020c, see our Sect. 4.2).

Color-magnitude plots
We now study the cluster galaxy population of the three KiDS clusters of this work by means of color-magnitude (CM) plots, with a particular focus on the BCGs. To this aim we select sources in the KiDS DR3 within a separation of 2 arcmin from each of the three targeted BCGs.
In Fig. 7 we report the corresponding g-r versus r plots, where g and r are AB apparent magnitudes, obtained with the OmegaCAM camera on the VST (see Sect. 2.1). At the redshifts of our targeted BCGs the rest frame 4000 Å break is redshifted to ∼ (5300 − 5800) Å, between the g and r filters, whose effective wavelengths are 4800 Å and 6250 Å, respectively. The chosen filters are thus optimal to detect the red sequence in the considered clusters.
In each CM plot we highlight the BCGs (green stars). We also distinguish field galaxies (gray dots) from fiducial cluster members, which have been selected as follows. Photometrically selected cluster members (filled blue dots) are those with AM-ICO cluster membership probabilities greater than 50%, which implies that the galaxies are formally more likely to belong to the clusters than the field. The selected threshold potentially allows us to reach a completeness level of 90% and at the same time limit the contamination from the field (e.g., George et al. 2011;Castignani & Benoist 2016). Similarly to previous studies (e.g., Rozo et al. 2015), we also selected spectroscopic cluster members (open red symbols) as those with a line-of-sight velocity separation δ < 2000 km s −1 with respect to the BCG, where δ = c |z spec −z BCG | 1+z BCG , c is the speed of light, z BCG is the BCG redshift, and z spec is the spectroscopic redshift of the fiducial cluster member.
The CM plots show that our targeted BCGs are indeed the brightest (and the most massive) among the selected cluster members. There is only one exception, which is a photometrically selected member with r=18.0 and a projected separation of 63 arcsec (∼0.3 Mpc) in projection from the KiDS 1220 BCG. The galaxy is thus much (∼ 1 mag) brighter than BCG and we suspect it is a foreground source, despite its photometric redshift of 0.46 ± 0.07.
In all CM plots we also appreciate a clustering of red sources at g-r∼ 1.5. However, the observed clustering is mitigated by the photometric redshift selection of cluster members. Interestingly, the targeted BCGs are systematically bluer. These findings suggest that the considered cluster cores are sufficiently mature with their galaxy populations well evolved in terms of color properties. On the other hand, the bluer colors of the BCGs are consistent with the fact that we selected them to be star-forming so that they are caught in a rare phase of their evolution. These results motivated us to perform a red-sequence modeling of the photometric data, similarly to previous studies (Kotyla et al. 2016;Castignani et al. 2019), as described in the following.
To model the red sequence, we used the Galaxy Evolutionary Synthesis Models (GALEV) tool 6 (Kotulla et al. 2009). We adopt chemically consistent models of elliptical galaxies, according to the GALEV prescriptions. We impose a galaxy formation at high redshift (z = 8), consistent with formation scenarios of massive ellipticals (Partridge & Peebles 1967;Larson 1975;De Lucia et al. 2006) We verified that changing the formation redshift between z = 6.5 − 20 does not significantly change the red-sequence model in the CM plot.
We also assume that the total initial mass of the galaxy is between 1 × 10 10 M and 5 × 10 11 M . For our red-sequence models we further rely on a Kroupa (2001) IMF between 0.1 and 100 M , which agrees well with the Chabrier (2003) IMF. We note that the SED modeling of the BCGs (Sect. 2.2.1) relies on the Chabrier (2003) IMF, which is, however, not included in the GALEV library.
With the above assumptions we performed different redsequence models, which we report as horizontal lines in Fig. 7, as further outlined in the following. For all three CM plots we show a model assuming no burst, which yields the reddest red sequence possible at the cluster redshift. We also report different models assuming a single strong burst of star formation, exponentially declining with time. For these models we assume a burst duration of 2 Gyr (e-folding time) and a burst strength of 100%, which is the fraction of the available gas at the onset of burst that is transformed into stars during the burst. We re-port different bursty red-sequence models, corresponding to single bursts of star formation occurring between z = 0.6 − 2.5, as shown in the legend of Fig 7. Interestingly, models with a star formation burst at z = 2.5 (dashed lines in Fig 7), that is, approximately at the peak of the cosmic SFR density (Madau & Dickinson 2014), yield colors that are quite similar to those inferred from models with no burst of star formation (solid lines in Fig 7). Furthermore, while these two models reproduce fairly well the colors of red cluster members, those of the star-forming BCGs are bluer. This analysis motivates us to perform additional evolutionary models. We indeed find that recent bursts at z ∼ 0.6 − 0.7, z ∼ 0.8 − 0.9, and z ∼ 1.2 − 1.3 are needed to better reproduce the observed g-r colors, and also the observed r-band magnitudes of the three KiDS 0920, 1220, and 1444 BCGs, respectively. For KiDS 1444 BCG at the highest redshift, a burst at z = 2.5 still reproduces fairly well the BCG g-r color. However, a later burst at z ∼ 1.2 − 1.3 appears to be favored as it implies both a bluer color and a brighter magnitude, more similar to those of the BCG. However, we note that neither of the adopted red-sequence models is able to reach the low r-band magnitude of KiDS 1444 BCG, despite the high 5 × 10 11 M initial mass adopted in the modeling. This possibly suggests that a faster mass growth (e.g., via mergers) than that predicted by the close-box GALEV evolutionary models would be a viable alternative for the BCG.
Overall, our analysis confirms that the three considered clusters have some evidence of a red sequence, well modeled assuming either no burst or an old burst of star formation at cosmic noon (z = 2.5). The associated BCGs show instead bluer colors, which are better reproduced assuming a recent burst of star formation. These results agree well with the fact that the BCGs are star-forming and gas-rich. With our observations and the modeling we are not able to explain precisely the physical origin for the observed high star formation, the blue colors, and the large gas reservoirs. However, recently accreted pristine gas may explain the observed properties, a scenario that we further discussed in Sect. 5.2.

Double-horn modeling for K1220. Structural parameters
We now specifically investigate the structural properties of KiDS 1220 BCG. As discussed in Sect. 5.1 this BCG shows evidence for a double-horn profile in the CO(3→2) emission. This is indicative of a low concentration of the molecular gas toward the center of the galaxy, which results in the observed flux dimming at small velocities in Fig. 3. Asymmetric spectra as the one observed for the KiDS 1220 BCG are commonly observed in HI in many galaxies (Richter & Sancisi 1994;Haynes et al. 1998;Matthews et al. 1998;Yu et al. 2020), as the atomic gas reservoirs are more extended and fragile than the molecular gas. Indeed the latter is usually more centrally peaked, with an exponentially decaying profile (Nishiyama et al. 2001;Regan et al. 2001;Leroy et al. 2008).
However, asymmetric CO spectra similar to the one found in KiDS 1220 are also reported in previous studies, mostly for lowz galaxies (Wiklind et al. 1997;Combes et al. 2007;Castignani et al. 2022b), including BCGs (Salomé et al. 2015). Following Wiklind et al. (1997) and Salomé et al. (2015) we attempt to link the observed double-horn spectrum to the structural parameters of the gas reservoir, as follows.
The observed velocity integrated S CO(3→2) ∆ spectrum can be related to the H 2 gas density profile n(r) and the rotational According to the model, both the scale length and the central density can be related to the DM mass (Salucci et al. 2007;Memola et al. 2011), and consequently V rot,DM (r) 2 . In order to reproduce a flattening of the rotational velocity profile at a large distance r ∼ 100 kpc, as it is the case of massive galaxies such as KiDS 1220 BCG, we assume a DM mass of M DM = 5×10 12 M , which implies ρ 0 = 2.40 × 10 −25 g cm −3 , r 0 = 44.2 kpc, a stellar to DM halo mass ratio of M /M DM = 6% and a flattening of the rotation velocity to a value of ∼ 300 km s −1 .
To determine the structural parameters of the gaseous reservoir we then vary the disk scales d, d 1 , and d 2 from 50 pc up to 20 kpc. From the integration in Eq. 5 over all radii and over the velocity support of the CO(3→2) line we can relate the observed spectrum to the one inferred from our modeling. The overall normalization is fixed by the height of the line. Increasing the inclination angle i has the net effect of increasing the velocity separation between the two peaks in the spectrum (Wiklind et al. 1997). As M H 2 is fixed, decreasing the disk scales d, d 1 has the net effect of increasing the central density n 0 , and thus the concentration, which ultimately reduces the observed depth in the spectrum at small velocities, with respect to the peaks. Given these relations, we find that both a ring and a disk with an inclination of i = 40 deg reproduce the observed spectrum well. The inferred sizes are d = d 1 = 9 kpc and d 1 − d 2 = 1 kpc. With these parameters we get peak H 2 -gas-density values of 800 and 300 M /pc 2 in the case of a disk and a ring, respectively, which are typical values found for LIRGs, like our BCGs, along the Kennicutt-Schmidt relation. Figure 8 (top) displays the gas surface density profiles. In the central panel we report the rotational velocity profile in the case of a d = 9 kpc disk, while similar curves are obtained in the case of a ring with the above parameters. In the bottom panel of Fig. 8 we show both the observed spectrum and our modeling in both cases of a gaseous disk or a ring. As the central density scales as the square of the inverse of disk scale, we stress that by adopting a smaller radius for the molecular gas reservoir we significantly increase the central density up to high values 1000 M /pc 2 , which are less likely. At the same time we also increase the concentration and we are not able anymore to reproduce well the flux dimming at small velocities in the spectrum, even assuming an extreme edge-on orientation (i = 90 deg).
Overall, our results thus yield an extended molecular gas reservoir that is reminiscent of the extended disk phase described in Olivares et al. (2022) for local BCGs. During this phase the star formation and so the molecular gas is more likely to be detected. Similar to what was suggested by these authors, the angular momentum may prevent the cold gas from being rapidly accreted, which would explain the low concentration and the large extent of molecular gas in KiDS 1220. Following the evolutionary scenario proposed by Olivares et al. (2022), the BCG may then experience a rapid compaction phase, similarly to that previously suggested for distant star-forming galaxies (Barro et al. 2013(Barro et al. , 2017, including BCGs (Castignani et al. 2020b).

Summary and conclusions
In this work we have investigated the molecular gas and star formation properties of three star-forming BCGs at intermediate redshifts, drawn from the Radovich et al. (2020) sample of BCGs within the third data release of KiDS (de Jong et al. 2017;Kuijken et al. 2019). We have selected three targets from a large sample of 684 spectroscopically confirmed z > 0.3 BCGs in the equatorial KiDS-N field, namely KiDS 0920 (z = 0.3216), KiDS 1220 (z = 0.3886), and KiDS 1444 (z = 0.4417), that show significant infrared emission (Sect. 2.1). Our sources thus represent a small fraction (0.4%) of the parent sample and are among the most star-forming BCGs in KiDS. They have line-based (Hα, [O II]) and SED-based SFR (20 − 30) M , as well as dust masses and luminosities typical of LIRGs (Sects. 2.2.1 and 2.2.2). By performing BPT line diagnostics, we find KiDS 0920 is classified as a dusty AGN, while the other two BCGs likely fall in the star-forming region of the BPT diagram (Sect. 2.2.2). The selected BCGs thus constitute a rare subpopulation of starforming BCGs that are experiencing stellar mass assembly and are thus caught in a special phase of their evolution, similar to other distant star-forming BCGs (e.g., Castignani et al. 2020a,b).
In order to investigate the mass growth and the fueling of gas feeding star formation, we observed the three BCGs with the IRAM 30m telescope, targeting the first three CO transitions (Sect. 3). We clearly detect at S/N = (3.8−10.2) all three sources in multiple transitions. Our campaign thus yields a high success rate of 88%, much higher than that of previous studies that targeted distant BCGs in CO (Sect. 5.1). With the present work we also double the number of distant BCGs that are clearly detected in at least two CO transitions.
We then compared the molecular gas properties of the KiDS BCGs with those of a sample of 84 distant cluster galaxies (Sect. 4) with molecular gas observations from the literature as well as with scaling relations of MS field galaxies. Our BCGs appear to be rich in gas, with M H 2 (0.5 − 1.4) × 10 11 M and M H 2 /M (0.2 − 0.6). These results place them among the most gas-rich BCGs at intermediate redshifts.
The KiDS BCGs of this work have a high content of molecular gas, well above the MS. Conversely, they also have star formation activities typical of MS galaxies (Sect. 5.2). This implies that our targeted BCGs accumulate large molecular gas reservoirs that are not efficiently converted into stars, having depletion timescales τ dep (2 − 4) Gyr, higher than those previously found in distant star-forming BCGs from the CLASH and SpARCS surveys (Fogarty et al. 2019;Castignani et al. 2020a;Dunne et al. 2021).
All these results are intriguing and suggest that the gas fueling the three KiDS BCGs differs somehow from that of previously studied star-forming BCGs at intermediate redshifts. Consistent with this interpretation, we find a correlation between the excitation ratio and the star formation activity of the BCGs (Sect. 5.3), in the sense that the highest excitation ratios are found only in highly star-forming SFR(> 100 M ) BCGs, which also have short depletion times of ∼ (0.5 − 0.7) Gyr (Fogarty et al. 2019;Castignani et al. 2020a). Our KiDS BCGs instead have lower excitation ratios, lower SFRs ∼ (20 − 30) M , and longer τ dep .
We further investigated the star formation properties of the KiDS BCGs by performing a CM analysis (Sect. 5.4). We find that the BCGs are bluer and brighter than photometrically selected cluster members. Red-sequence modeling suggests that recent bursts of star formation are needed to better reproduce both the observed g-r colors and the r-band magnitudes of the three BCGs.
One possible explanation for the observed phenomenology is that the H 2 gas, or a substantial amount of it, was recently accreted by the KiDS BCGs but still not efficiently converted into stars. The infall of pristine gas may also explain the M H 2 /M dust 170 − 300 ratios that we find for our KiDS BCGs. This accretion may have occurred via filaments (as may be the case of KiDS 0920, which shows an elongated morphology) or via interaction with companions, which are in fact found for KiDS 1220 and 1444.
Interestingly, for KiDS 1220 we have additionally found a double-horn emission in CO(3→2), which implies low gas concentration toward the galaxy center. The CO(3→2) emission is well modeled with an extended molecular gas reservoir of ∼ 9 kpc (Sect. 5.5), which is reminiscent of a mature extended disk phase observed in local BCGs by Olivares et al. (2022). This phase may precede a compaction phase (Barro et al. 2013(Barro et al. , 2017Castignani et al. 2020b), and later the quenching.