Suppressed cooling and turbulent heating in the core of X-ray luminous clusters RXCJ1504.1-0248 and Abell 1664

We present the analysis of XMM-Newton observations of two X-ray luminous cool core clusters, RXCJ1504.1-0248 and Abell 1664. The Reflection Grating Spectrometer reveals a radiative cooling rate of $180\pm 40\, \rm M_{\odot}\rm\,yr^{-1}$ and $34\pm 6\, \rm M_{\odot}\rm\,yr^{-1}$ in RXCJ1504.1-0248 and Abell 1664 for gas above 0.7 keV, respectively. These cooling rates are higher than the star formation rates observed in the clusters, and support simultaneous star formation and molecular gas mass growth on a timescale of 3$\times 10^8$ yr or longer. At these rates, the energy of the X-ray cooling gas is inadequate to power the observed UV/optical line-emitting nebulae, which suggests additional strong heating. No significant residual cooling is detected below 0.7 keV in RXCJ1504.1-0248. By simultaneously fitting the first and second order spectra, we place an upper limit on turbulent velocity of 300 km$\rm s^{-1}$ at 90 per cent confidence level for the soft X-ray emitting gas in both clusters. The turbulent energy density is considered to be less than 8.9 and 27 per cent of the thermal energy density in RXCJ1504.1-0248 and Abell 1664, respectively. This means it is insufficient for AGN heating to fully propagate throughout the cool core via turbulence. We find the cool X-ray component of Abell 1664 ($\sim$0.8 keV) is blueshifted from the systemic velocity by 750$^{+800}_{-280}$ km$\rm s^{-1}$. This is consistent with one component of the molecular gas in the core and suggests a similar dynamical structure for the two phases. We find that an intrinsic absorption model allows the cooling rate to increase to $520\pm 30\, \rm M_{\odot}\rm\,yr^{-1}$ in RXCJ1504.1-0248.


INTRODUCTION
The inner core region of relaxed clusters of galaxies shows a complex structure of different gas phases. Most of the gas mass is collisionally ionised and cooling via thermal bremsstrahlung and line emission in X-rays onto the central Brightest Cluster Galaxy (BCG). The low temperature and high density of the cool core are indicative of a short radiative cooling time of < 1 Gyr (e.g. Fabian 1994;Fabian et al. 2006;Blanton et al. 2004;McNamara et al. 2006;Panagoulia et al. 2014b;Liu et al. 2019). It predicts a massive cooling flow in the most massive clusters, such as A1835 and the Phoenix cluster (e.g. Allen et al. 1996;McDonald et al. 2012McDonald et al. , 2014. However, spectral evidence from the Reflection Grating Spectrometers (RGS) onboard XMM-Newton only supports a mild cooling rate in cool cores, typically less than 10 per cent of the predicted rate in the absence of heating (e.g. Kaastra et al. 2001;Peterson et al. 2003;Liu et al. 2019). It requires a heating process that needs an energy source and an efficient mechanism to transport the energy throughout the core. The central ★ E-mail: hl479@cam.ac.uk AGN interacts with its host environment. For cool core clusters at a low Eddington fraction, AGN feedback operates in the kinetic mode, where gas accretion generates powerful relativistic jets which inflate bubbles (Fabian 2012;McNamara & Nulsen 2012). Bubbles rise buoyantly with a mechanical power similar to the cooling rate in the absence of heating (e.g. Churazov et al. 2002;Bîrzan et al. 2004;Dunn & Fabian 2006;Rafferty et al. 2006;Hlavacek-Larrondo et al. 2012). The temperature map of clusters is roughly isotropic, which suggests heating occurs away from the jet axis. The mode of such energy transfer is still under debate. While the energy can be propagated azimuthally by gravity waves, it can not transport the energy radially. An alternative mode of powerful sound waves can provide the required velocity for radial energy transport (Fabian et al. 2003b(Fabian et al. , 2017, but a suitable energy dissipation mechanism needs to be developed. Turbulent heating of the gas has also been of interest for this purpose. The Hitomi Soft X-ray Spectrometer (SXS) made an accurate measurement of the level of turbulence at 164±10 km s −1 in the Perseus cluster (Hitomi Collaboration et al. 2016). The energy density of isotropic turbulence is only 4 per cent of the thermal energy density there which is too low to reach the full cooling core. Turbulence alone is insufficient to offset radiative cooling.
It is possible to measure an upper limit to the level of turbulence using the RGS in many X-ray peaked clusters (e.g. Sanders et al. 2010Sanders et al. , 2011Sanders & Fabian 2013). Since the RGS is a slitless spectrometer, the spectral lines are broadened by the spatial extent of the source in addition to other broadening processes. The spatial broadening follows the RGS dispersion law, Δ = 0.138Δ / Å, where Δ is the broadening in wavelength, Δ is the angular offset from the central source in arcmin and is the spectral order. It contributes significantly to the total line width in nearby sources (e.g. Pinto et al. 2015). To obtain a conservative limit, this artificial broadening can be corrected for by using the surface brightness profile of the European Photon Imaging Camera (EPIC) image (e.g. Pinto et al. 2018;Bambic et al. 2018). A tight 90 per cent upper limit of 244 km s −1 is found in A1835 and 246 km s −1 in A2204 . These limits are similar to the level of turbulence in the Perseus cluster found by the Hitomi Collaboration et al. (2016).
Another interesting feature is the presence of H emission in most cool core clusters. Many studies of the inner cluster core have shown that H emission in the form of filaments is spatially aligned with the soft X-ray emitting gas and the two gas phases are likely mixing (e.g. Fabian et al. 2003aFabian et al. , 2006Fabian et al. , 2016Crawford et al. 2005b). No strong evidence of significant cooling below ∼1 keV suggests that the soft X-ray component is likely not cooling radiatively, but is mixing and powering the observed optical/IR emission due to the atomic and partially-ionised gas (Fabian et al. 2003a). This situation can occur if the hot X-ray component interpenetrates the cold H nebula and creates fast and energetic particles . The fast particles can then heat and excite the cold gas, powering the observed nebulosity (Ferland et al. 2009). In a previous work, we have shown that the thermal energy of the radiative cooling gas is sufficient as the power source for the optical/UV nebula in clusters with a cooling rate below ∼10 M yr −1 , but the most luminous clusters are likely powered by hotter gas or otherwise (Liu et al. 2019). Churazov et al. (2013) argued that buoyant bubbles stretch fluid elements to form gaseous filaments with amplified magnetic field. The release of magnetic energy allows dissipation into filaments. Alternatively, H filaments can also powered by Cosmic Rays (Ruszkowski et al. 2018).
The origin and fate of the molecular gas is another mystery. A massive cold molecular gas reservoir is often present in the core and seen by CO lines from a component at ∼ 50K and/or NIR H 2 lines at ∼ 2000K (Edge 2001;Edge et al. 2002;Salomé & Combes 2003;Wilman et al. 2009;Olivares et al. 2019;Russell et al. 2019). Star formation of up to hundreds of M yr −1 in the most massive clusters is a major consumer of the molecular gas deposit. At the higher rates, the observed molecular gas reservoir will be depleted by star formation in 10 8 -10 9 yr if not replenished (e.g. Pulido et al. 2018). On the other hand, this timescale is comparable to the central radiative cooling time, which suggests the molecular gas cools from the hot X-ray atmosphere (e.g. Russell et al. 2019). The molecular gas filaments have a smaller spatial extent and are often embedded in the H nebula and hence the soft X-ray gas (e.g. Olivares et al. 2019;Russell et al. 2019). Surprisingly, the RGS spectra have revealed that the molecular gas mass is comparable to the X-ray gas mass emitting below 1 keV in a sample of nearby luminous clusters, e.g. 2A0335+096, A2052, A3581 (Liu et al. 2020). These two gas phases are likely intermingled and the structural integrity is held by magnetic fields.
Both of our targets, RXCJ1504 and A1664, are remarkably luminous in both the X-ray and optical bands, and possess a massive molecular gas reservoir. RXCJ1504 is one of the most massive low redshift cool core clusters at = 0.2153 with 500 = 1.25×10 15 M (Piffaretti et al. 2011) and a high X-ray bolometric luminosity of 4.1×10 45 ergs −1 and a classical cooling rate 1 of 1500-1900 M yr −1 (Böhringer et al. 2005). Giacintucci et al. (2011) discovered a minihalo of 140 kpc in radius at the centre of the cluster confined to the cool core. This suggests a tight connection between the X-ray emitting cool core and the relativistic plasma. It also has an H luminosity of 3.2×10 43 ergs −1 making it one of the most luminous optical nebulae. The observed UV flux indicates a strong star formation rate of 130 M yr −1 (Ogrean et al. 2010). The inner 5 kpc of the cool core contains most gas from the massive molecular gas reservoir of 1.9±0.1×10 10 M . The kinematics of the molecular gas is complex as revealed by ALMA CO observations (e.g. Vantyghem et al. 2018). Vantyghem et al. (2018) infer a turbulent velocity of 335±15 km s −1 for that gas and the central molecular filament shows a velocity range of ∼260 ±11 km s −1 in RXCJ1504. A dynamically young gas structure also shows an offset velocity of ∼ -250 km s −1 from the rest of the gas in the BCG.
A1664 is one of the first cool core clusters in which CO emission was observed (e.g. Edge 2001). It has a redshift of = 0.1283 with 500 = 4.06 × 10 14 M (Piffaretti et al. 2011 (Wilman et al. 2006). The star formation rate is estimated to be 14 M yr −1 in IR or 4.3 M yr −1 in FUV (O'Dea et al. 2010). The BCG has a total molecular gas mass of 1.1±0.1× 10 10 M (Russell et al. 2014). The molecular gas is also seen disturbed within 10 kpc of the core (Russell et al. 2014). The CO(1-0) and CO(3-2) lines are well resolved into two Gaussian components with a velocity difference of ∼590 km s −1 . On a larger scale of ∼50 kpc, cold fronts are observed in the X-ray atmosphere produced by sloshing (Calzadilla et al. 2019), and it is possible that core sloshing can affect lower temperature gas. If the X-ray and molecular gas are related, they are likely sharing a similar velocity structure (for theoretical modelling, see e.g. Gaspari et al. 2017).
At the present time, the XMM-Newton RGS can place the most accurate constraint on the velocity of the soft X-ray emitting gas. The dispersive nature of the RGS means the spectral resolution is improved with lower photon energies, and surpasses the Hitomi/XRISM resolution below 1 keV (12.4 Å) for point-like and extended sources below 1 arcmin of spatial extent.
In this work, we present recent deep XMM-Newton/RGS observations of these two X-ray luminous clusters: RXCJ1504.1-0248 (RXCJ1504) and Abell 1664 (A1664). We measure radiative cooling rates and place constraints on turbulent velocity at the 90 per cent confidence level, which is an important proxy for heat propagation. The structure of the paper is as follows. Section 2 provides observational details of the clusters and the data reduction procedures. Section 3 introduces the spectral models used to measure the cooling rate and place the upper limit on turbulent velocity. Section 4 discusses the implications of our results and we try to correct for intrinsic absorption of the clusters. We assume the following cosmological parameters: 0 = 73 km −1 Mpc −1 , Ω M = 0.27, Ω Λ = 0.73.

OBSERVATIONS AND DATA REDUCTION
The XMM-Newton observatory observed each of the clusters RXCJ1504 and A1664 for two orbits (PI Fabian). The observational details are listed in Table 1. RXCJ1504 was observed between 15-Aug-2019 and 17-Aug-2019 and between 09-Feb-2020 and 10-Fen-2020. The offset of the roll-angle of the pointings between observations is 171.65 degrees. A1664 was observed between 28-Jan-2020 and 29-Jan-2020 and between 30-Jan-2020 and 31-Jan-2020. The offset of the roll-angle of the pointings is 0.65 degrees.
Here we used data from the RGS and the EPIC onboard XMM-Newton. We follow the standard data reduction procedure with the latest XMM-Newton Science Analysis System v 17.0.0. We extract the first and second order RGS spectra by the SAS task rgsproc. The second order spectra possess twice the spectral resolution and hence are used for turbulent velocity measurements. We set the xpsfincl mask to include 90% of the point spread function. This is equivalent to a narrow 0.9 arcmin region. We use template background files based on count rates in CCD 9 to produce background-subtracted spectra. To achieve the highest S/N ratio, the RGS 1 and 2 spectra of both observations are stacked using the task rgscombine and then processed by the task trafo to be analysed by SPEX. We check that the pointing of both observations is consistent to avoid spurious broadening of emission lines.
The spatial broadening of the RGS spectra is corrected by the surface brightness profile of the MOS image. The MOS cameras are aligned with the associated RGS detectors and have slightly better spatial resolution than the pn detectors. We only use MOS1 images from the earlier observation for each object (0840580101 for RXCJ1504 and 0840580301 for A1664). The images are produced by the SAS task emproc. We extract the surface brightness profiles in the 0.5-1.8 keV energy band using the task rgsvprof.
We used SPEX version 3.05.00 for spectral analysis with its default proto-Solar abundances of Lodders & Palme (2009) and ionisation balance (Urdampilleta et al. 2017). The spectral fitting uses C-statistic (C-stat) minimisation which is equivalent to 2 in low count statistics. We adopt 68.3 per cent confidence level (1 uncertainty at Δ = 1) for measurements. For upper/lower limits, we only quote the 90 per cent confidence level (2 uncertainty at Δ = 2.71) uncertainty, unless otherwise stated.

RESULTS
The stacked RGS spectra are binned by a factor of 3 to be consistent with the spectral resolution and preserve most spectral information. We fit the first order spectra over the 7-28 Å band where the background is lower than the continuum. We include the 7-20 Å band for the second order spectra to use the most spectral information.
We use the collisional ionisation equilibrium component (cie) and the cooling flow component (cf ) available in SPEX to construct our cooling flow models as described in Liu et al. (2019). The cie component represents a plasma with a free temperature and emission measure = e H , where e and H are electron and proton densities and is the volume of the emitting gas. We use the default value of e in SPEX. It is typically used in single-temperature and two-temperature models to describe a cluster. The cf component consists of a set of cie components and calculates the differential emission measure to match that of the required cooling rate where is the Boltzmann constant, is the mean particle weight, H is the proton mass and Λ( ) is the cooling function. The maximum temperature of the cf component is coupled to cie component and we assume both components have the same abundances.
To reduce the number of free parameters, we fix the Ne/Fe and Mg/Fe ratios for both clusters. We set Ne/Fe=0.8 and Mg/Fe=0.75 in RXCJ1504 and Ne/Fe=Mg/Fe=0.6 in A1664. These ratios are measured from a 1cie+1cf model (Model 2) and do not change with additional cf components. The abundances of the other elements are coupled to Fe. The cie and cf components are modified by redshift, cold Galactic absorption with solar abundances and spatial broadening (lpro; Pinto et al. 2015). The lpro component uses the surface brightness profile as the input. The scale factor and the wavelength shift Δ are the free parameters. The scale factor fit for the amount of line broadening and the wavelength shift corrects for the centroid of emission.

Cooling flow analysis
To construct the cooling flow models, we first model the hot plasma (> 2 keV) in the multi-phased intracluster medium (ICM) with a cie component (Model 1). Three cooling flow models are then considered combining cie and cf component: complete (one-stage), onestage with a free minimum temperature and two-stage models. We define the 'complete' cooling rate as the rate measured from the cie temperature down to the minimum temperature of the cf component of 0.01 keV (Model 2). This minimum temperature of 0.01 keV is the lowest possible value allowed in SPEX. This one-stage cooling flow model is often sufficient for the spectra of clusters and groups with low statistics and a low cie temperature (e.g. Liu et al. 2019), but not necessarily for RXCJ1504.1-0248 and A1664. We then free the minimum temperature of the cf component to include the possibility that the ICM stops cooling radiatively in X-rays at a higher temperature (Model 3). This also leads to a 'two-stage' cooling flow model that has two cf components, where the cooling rates are measured between the cie temperature and 0.7 keV and between 0.7 keV and 0.01 keV, respectively. We refer to the cooling rate between 0.7 keV and 0.01 keV as the residual cooling rate in this work (Model 4). In a previous work, Liu et al. (2019) discussed the effect of the transition temperature between two cooling flow components on the cooling rates in the two-stage model. For a high transition temperature up to 0.9 keV, the cooling rate above the transition temperature is likely increased by 20 per cent. For a low transition temperature, the cooling rate is likely decreased by 10 per cent, while the residual cooling rate is over-predicted due to a narrow temperature range. We found that the transition temperature of 0.7 keV is suitable for fitting the Fe XVII lines and its forbidden-to-resonance line ratio. It is also consistent with the one-stage cooling flow model with a free minimum temperature.
The cooling rates, cie temperatures and O and Fe abundances of three cooling flow models are detailed in Table 2 and 3. We show the stacked RGS spectra in Fig. 1 and 2 with the best fit cooling flow models.
In RXCJ1504, we find that both the one-stage cooling flow model with a free minimum temperature (Model 3) and the two-stage model (Model 4) yield the minimum C-stat for the same number of degrees of freedom (DoF). The transition temperature of the two-stage model is consistent with the free minimum temperature of the onestage model. The other fit parameters are also consistent between these two models. We hence conclude a cooling rate to 0.7 keV of 180±40 M yr −1 and a residual cooling rate of from 0.7 keV less than 53M yr −1 at 90 per cent confidence level. The Fe XVII resonance line is seen in the spectrum and mixed with a broad feature at 15 Å in rest wavelength (18.3 Å in observed wavelength). This indicates a cooling flow is present at around 0.7 keV. There are several possibilities for the nature of the broad feature. First, the Fe XVII resonance line is suppressed in the line of sight and re-emitted from the outer region. The spatial extent of the gas emitting the Fe XVII resonance line is broadened which results in a broader line at 15 Å due to the fact that the RGS detectors are slitless. Second, the gaseous neutral iron in the interstellar medium has 2 deep and broad absorption edges at 17.2 Å and 17.5 Å . However, most iron is in dust, which has a different edge shape and position. A spectral modelling of the iron edge which does not account for dust might introduce some systematic effects including a spurious bump around 18 Å (see, e.g., Juett et al. 2006;Pinto et al. 2013).
The O VII triplet is not observed. Due to the high continuum emission of the hot gas, the mass of cold gas at 0.2 keV is difficult to detect. The distinction between the complete and two-stage models is statistically significant. While some fit parameters are consistent such as the cie temperature and metallicities, the two-stage model gives 3.6 times higher cooling rate above 0.7 keV than the complete cooling rate. By comparison in Liu et al. (2019), we find such a ratio of 2.2, 1.7, 3.8 in A262, Centaurus and M87, respectively, all of which show a significant statistical improvement in the two-stage model. Given the 20-40 per cent uncertainty on the measurements, this ratio is broadly consistent with other nearby cool core clusters (Liu et al. 2019).
In A1664, we find that the cooling flow models are improved by using a second line broadening component for the cf components. The second lpro component uses the same surface brightness profile and we fit the scale factor and wavelength shift as the lpro component for the cie. The extra line broadening component improves the Cstat by 5 in comparison with using just one broadening component in the complete cooling model. We find that all the cooling flow models provide a similar fit to the spectrum with the same C-stat and consistent cooling rate. As a result, we only report a complete, onestage, cooling rate of 34±6M yr −1 from the current data (Model 1).

Spatial broadening
The observed line broadening in RGS spectra is the sum of thermal broadening, turbulent motion and spatial broadening. Thermal broadening is already calculated in the thermal components such as cie or cf. To place constraints on the turbulent velocity, we need to estimate the level of spatial broadening. Since the scale of the hot ICM is much larger than that of the cool core, using the full spatial profile over-predicts the contribution to the spatial broadening and hence underestimate the turbulence. We follow the method used in Bambic et al. (2018) and Pinto et al. (2018) for a more accurate estimate of spatial broadening due to the cool gas.
The SPEX task rgsvprof gives a cumulative flux of the surface brightness profile of the MOS images as a function of wavelength. This can be inverted into a Gaussian shaped profile as expected from the image. Such a profile can be modelled by the sum of three Gaussians. The central narrowest Gaussian represents the coolest gas in the core. The bremsstrahlung continuum from the hot ICM is seen in the broadest outer Gaussian. The remaining intermediate Gaussian provides the transition between the ICM and the cool core. As we try to measure the turbulence in the cool core, the central and intermediate Gaussians are the relevant components in the estimation of spatial broadening.
The surface brightness profiles of RXCJ1504.1-0248 and A1664 are seen in Fig. 3. We find that the profile of A1664 is skewed and the asymmetry is seen in the DETY direction of the MOS 1 image. The separation between the centre of the central and intermediate Gaussian is 0.045±0.001 Å. From the RGS dispersion law, such a wavelength separation corresponds to a physical separation of 70 kpc. This means the intermediate Gaussian component is indeed at the rim of the cool core.
We reconstruct two profiles of cumulative flux that can be used in SPEX, which include either using the central Gaussian alone or using both the central and intermediate Gaussians. C (M yr −1 ) n/a n/a n/a <53 Model 1 is the single-temperature (1cie) model, model 2 is the complete cooling (1cie+1cf ) model, model 3 is the one-stage model with a free minimum temperature and model 4 is the two-stage (1cie+2cf ) model. H and min are the cie temperature and the minimum temperature of the associated cf component, respectively. H is the cooling rate between H and min . C is the residual cooling rate between 0.7 and 0.01 keV in the two-stage model. C (M yr −1 ) n/a n/a n/a 30±10 Figure 2. As Fig. 1, the stacked RGS spectrum of A1664 with the best fit complete cooling flow model in red in rest wavelength.

Turbulent velocity measurements
We simultaneously fit the first and the second order spectra to measure the turbulent velocity. The observed spatial profile is replaced by the profiles reconstructed from the Gaussian approximations in the lpro component. To conserve the RGS dispersion law, we then set the scale factor of the lpro to = 1 for the first order and = 0.5 for the second order spectra. We also fit the wavelength shift parameter in the lpro component to adjust for redshift. We use the singletemperature (1cie) model and fit the micro-turbulent velocity ( turb ) of the cie components. The 1-dimensional turbulent velocity is then 1D = mic / √ 2. The fit parameters between the two sectors (first and second order spectra) are then coupled.
We summarise the velocity limits in Table 4. The total line width due to turbulence and spatial broadening is calculated by using the full spatial profile and setting the scale factor to 0. The most accurate velocity limit is measured by simultaneously fitting the first and second order spectra. We find that correcting the spatial broadening both using the central Gaussian alone and using both the central and intermediate Gaussians give consistent velocity limits. We report that the best 90 per cent upper limit for RXCJ1504 and A1664 are both 300 km s −1 .

Soft X-ray and cooler gas
It is possible to achieve a 3 measurement of the cooling rate by the two-stage model (Model 4) or better for many nearby, < 0.01, clusters (Liu et al. 2019), but only a few at the redshift similar to our targets or higher (e.g. Tozzi et al. 2015;Pinto et al. 2018). From the analysis of the deep observations of the luminous clusters RXCJ1504 and A1664, we can provide reliable measurements of the cooling rate at the 4-5 confidence level. Both targets are already wellstudied in other energy bands as well as spatially resolved analysis in Chandra (e.g. RXCJ1504: Böhringer et al. 2005;Cavagnolo et al. 2009;Ogrean et al. 2010;Sanders et al. 2011;Vantyghem et al. 2018   of great interest to understand the role of such X-ray cooling rate in cluster evolution at intermediate redshifts. To be more precise, in this section, we discuss the connection between the soft X-ray emitting gas and the cooler materials including the H nebula, molecular gas reservoir, gas consumed by star formation activities and AGN accretion. We are currently analysing the archival spectra of luminous clusters at intermediate redshifts (0.1 < < 0.6) with known optical nebula and will report elsewhere.
The H nebulae of RXCJ1504 and A1664 are exceedingly luminous for intermediate redshift clusters. To power these partially ionised nebulae, at least a factor of 15 times the observed H luminosity is required to include other UV/IR emission due to the same gas (Fabian et al. 2003a;Ferland et al. 2009). O'Dea et al. (2010) and Ogrean et al. (2010) reported that the energy of stellar photoionisation is comparable to the H luminosity in our targets. This means that additional sources of energy are required to power the remaining emission. Fabian et al. (2003a) suggested that the soft X-ray gas can provide sufficient energy for the nebulae. This is supported by the spatial coincidence between the soft X-ray components and the H nebula (e.g. Perseus: Fabian et al. 2003aFabian et al. , 2006 et al. 2005a). Since most soft X-ray gas stops cooling radiatively below 0.7 keV as seen in the spectra, it can release a significant amount of energy if it continues to cool non-radiatively. For nearby clusters, we found that the energy of the 0.7 keV gas is sufficient only for the less luminous nebulae, while the most luminous nebulae require a much warmer gas (Liu et al. 2019). O'Dea et al. (2008) found the same conclusion for the 1 keV gas to power 5 times the IR luminosity. To calculate the energy required by the nebulae in our clusters into a mass inflow rate, we assume 15 times the H luminosity For our targets, neb is 4.6×10 3 M yr −1 for RXCJ1504 and 2.1×10 2 M yr −1 for A1664. Both of these values are much larger than the observed cooling rate between the ICM temperature and 0.7 keV. They are also 2-3 times of the classical cooling rate predicted in the absence of heating. The cooling flow at 0.7 keV is therefore insufficient to power the observed H nebulae. Alternatively, if the nebulae are powered by the warmer gas of the same rate as radiatively cooling, equation 3 suggests the temperature of the gas is 25 keV and 6.3 keV for RXCJ1504 and A1664, respectively. The temperature is much hotter than the temperatures of the hot ICM in RXCJ1504 and is not expected in the cool core. Therefore, other sources of significant energy are required to power the H nebulae of at least RXCJ1504 in addition to stellar photoionisation and the soft X-ray cooling flow (see Section 4.4 for additional energy in an alternative cooling flow model).
The gas properties of RXCJ1504 compare well to those of the Phoenix cluster at = 0.596. It has a similar molecular gas mass of 2×10 10 M   yr −1 below 2 keV at 68 per cent confidence level. Both the H luminosity and the cooling rate are twice of those measured in RXCJ1504. We calculate neb to be 4.3×10 3 M yr −1 for the 2 keV gas. This means the soft X-ray gas and stellar photoionisation are also insufficient as the power source. However, the Phoenix cluster has a star burst of 500-800 yr −1 that is comparable to the observed cooling rate at the 1 confidence level (McDonald et al. 2013(McDonald et al. , 2014. This suggests the molecular gas reservoir is likely growing slowly at the young age of the cluster. Although the energy produced by the cooling rates does not match the energy required by the H nebulae, the fate of the mass of the cooling gas still needs to be accounted for. The condensation of X-ray cooling gas is strongly linked to both the massive molecular gas reservoir and the star formation in the BCG, which are only present when the radiative cooling time falls below a Gyr (Rafferty et al. 2008;Pulido et al. 2018). Russell et al. (2019) and Liu et al. (2020) found that the mass of the soft X-ray gas is consistent with the molecular gas mass in the inner 10 kpc. If the X-ray cooling flow is indeed a major source of gas for the molecular gas reservoir and then star formation, we can calculate the timescale for forming the reservoir. Without any star formation activity and AGN gas accretion, the molecular gas requires 10 8 yr to accumulate in RXCJ1504 and 3.2×10 8 yr in A1664. However, both of our targets are extremely star forming clusters and may have a strong AGN activity. Our results show that RXCJ1504 is cooling at 10 per cent and A1664 is cooling at 34 per cent of the classical cooling rate predicted rate in the absence of heating. This means most radiative cooling is suppressed by AGN feedback. The amount of heating required can be deduced from the luminosity of the cooling flow component above 0.7 keV, which is available in SPEX. This indicates at least 2.25×10 45 ergs −1 and 2.1×10 43 ergs −1 for RXCJ1504 and A1664, respectively. Such energy is about a third of the mechanical power in A1664 and hence AGN feedback must supply the larger power of 6.8×10 43 ergs −1 . However, the required energy is 10 times larger than the mechanical power from the AGN in RXCJ1504, which suggests the AGN has been much more active than now observed. Ogrean et al. (2010) found the same conclusion in RXCJ1504 using the 3 upper limit of the cooling rate measured from archival EPIC/RGS spectra. Assuming accretion efficiency of 0.1, the required energy is equivalent to a black hole growth rate of 0.39 M yr −1 and 0.012 M yr −1 for RXCJ1504 and A1664, respectively. If the AGN is powered by Bondi accretion from the X-ray emitting gas, it requires a cool component of about 0.5 keV in RXCJ1504 (Bondi 1952;Ogrean et al. 2010). Although our cooling flow models find most gas is above 0.7 keV, the detection of the Fe XVII resonance line shows that it is likely to have some cool gas at around 0.5 keV. Liu et al. (2020) measured the mass of 0.7 keV in nearby cool core clusters of 10 8 − 10 9 M . RXCJ1504 is likely to have a higher gas mass at this temperature, since the luminosity of the 0.7 keV gas in the two-temperature model is 9 times larger than that of 2A0335+096, which has the largest gas mass below 1 keV. Such a cool gas can fuel the AGN on the timescale of a few 10 9 yr. The ratio of the black hole growth rate to the star formation rate is 0.003 and 0.00086 for RXCJ1504 and A1664, respectively. The relation of black hole growth and star formation is in good agreement with other clusters (Rafferty et al. 2006). Finally, the ratio of the radiative cooling rate to the star formation rate is between 1.5 and 2.5, which is smaller than most moderate star forming clusters but consistent with more luminous clusters such as A1835 (Rafferty et al. 2006;O'Dea et al. 2008;Liu et al. 2019). We find a net mass deposition rate of 50 M yr −1 in RXCJ1504 and 20-30 M yr −1 in A1664. These increase the molecular gas formation timescale by 2-3 times. Nevertheless, these timescales are consistent with the typical radiative cooling time of cool core clusters (e.g. A1664, Kirkpatrick et al. 2009), which suggests a strong link between the X-ray cooling gas and the molecular gas.

Turbulence versus heat propagation
The archival XMM-Newton observations of A1664 (ObsID: 0302030201/0302030201) did not point at the centre of the cluster and no turbulent velocity measurement was made by previous works. The previous spectroscopic analyses and Monte Carlo simulation of turbulence found a velocity of 670 +600 −360 km s −1 and 1310 +570 −670 km s −1 at 68 per cent confidence level, respectively Sanders & Fabian 2013). Our results using the new data show a much tighter limit (of 300 km s −1 ). The turbulent velocity of our targets is comparable to the velocity measured in many bright clusters, e.g. <211 km s −1 in A1835 (Sanders & Fabian 2013;Bambic et al. 2018), <400 km s −1 in 2A0335+096 (Pinto et al. 2015), ∼164 km s −1 in Perseus (Hitomi Collaboration et al. 2016) and <370 km s −1 in the Phoenix cluster . It is worth noting that resonant scattering can place constraint on the level of turbulence in elliptical galaxies in galaxy groups (de Plaa et al. 2012;Ogorzalek et al. 2017). For elliptical galaxies with a temperature below 1 keV, strong Fe XVII lines are usually seen in the RGS spectra. Ogorzalek et al. (2017) measured a mean turbulent velocity of 107±17 km s −1 in 13 elliptical galaxies. This is lower than the upper limit in our targets, but individual galaxies can have a higher turbulence (e.g. NGC 5044, de Plaa et al. 2012). The temperature of BCG in clusters is typically above 1.5-2 keV and Fe XVII lines are not detected in all clusters. It is also difficult to measure the Fe XVII resonance-to-forbidden ratio due to the high continuum. In the case of RXCJ1504, the Fe XVII forbidden line is redshifted to the RGS chip gap and therefore not detected by the RGS. It is ideal to place constraint on turbulence with resonant scattering in clusters with Fe L lines, which is beyond the scope of this work.
We can calculate the adiabatic sound speed s = √︁ / m p , where is the adiabatic index, which is 5/3 for ideal monatomic gas, = 0.6 is the mean particle mass and m p is the proton mass. This gives a sound speed of 1300 km s −1 and 750 km s −1 for RXCJ1504 and A1664, respectively. In this work, we calculate the 1-D Mach number for turbulence = 1D / s . The turbulent velocity then has a Mach number of 0.23 in RXCJ1504 and 0.4 in A1664. To calculate the ratio of the energy density in turbulence to the thermal energy of the plasma, we follow equation 11 in Werner et al. (2009) and obtain turb therm = 2 .
We find that the energy density ratio is less than 8.9 per cent in RXCJ1504, which is comparable to the ratio of 4 per cent in Perseus (Hitomi Collaboration et al. 2016) and 13 per cent in A1835 (Sanders et al. 2010). In A1664, the turbulence energy is less than 27 per cent of thermal energy. Bambic et al. (2018) and Pinto et al. (2018) calculated the minimum propagation velocity required to balance radiative cooling as a function of radius in 4 cool core clusters. The gas properties of RXCJ1504 are similar to those of the Phoenix cluster and A1835, while the core of 1664 is similar to that of A2204. The upper limits of turbulent velocity of 300 km s −1 are lower than required in both clusters at more than 15 kpc from the core. It is clear that the energetics of the turbulent motion of hot gas can not fully balance radiative cooling throughout the cool core. We now discuss the problem of energy transport and dissipation. Zhuravleva et al. (2018) argued that the large (∼ 10 per cent) surface brightness fluctuations in the X-ray images are isobaric and/or isothermal on spatial scales of 10-60 kpc and are likely associated with slow gas motions and bubbles of relativistic plasma (X-ray cavities). Bubbles tend to propagate along an axis but heating is also needed in directions away from that axis. This requires a faster propagation than turbulence alone. Internal waves or g-modes (buoyancy waves) are invoked on energetic grounds, but these waves do not propagate fast enough.
An alternative is to invoke time variability.

Blueshifted component in A1664
The best fit cooling flow models of A1664 show that the spectrum is well fitted by two lpro components (see section 3.1). This is also seen in some other clusters, e.g. Centaurus, where the Fe XVII lines are narrower than emission lines from hot-gas (Pinto et al. 2016;Liu et al. 2019). We find that although the scale factor is consistent in the lpro components in A1664, the wavelength shift is different. To understand the nature of this shift, we adopt a two-temperature model (2 cie) and each cie is associated with a separate lpro component. We find that the cooler component has a temperature of 0.80±0.08 keV and blueshifted by 0.046 +0.049 −0.017 Å from the hot gas. Such a difference can be achieved by either a blueshifted gas component or the different centroids of the hot and cool gas phase. We extract the surface brightness profile of the MOS1 image in 0.5-1 keV and 1-3 keV energy bands. These bands cover most emission seen in the core RGS spectrum. The profiles are shown in Fig. 4. We find that the centroids of different gas phases are separated by 0.0017±0.0006 Å. This only accounts for 4 per cent of the observed wavelength shift. Therefore, the blueshift is due to the motion of the cool gas. Assuming the blueshift is driven by the Fe XVII resonance line, we estimate the blueshifted velocity of 750 +800 −280 km s −1 . Independently, we also decouple the redshift of the two cie components and convolve both with the same lpro component. We find a consistent blueshifted velocity of 1000 +500 −300 km s −1 . Note that the difference of the roll-angle of the pointings between the two observations is small. By simultaneously fitting the spectra of individual observations, we find that the fit parameters are consistent with the stacked spectrum of both observations. A different line-of-sight velocity of different gas phases is seen in other clusters. By decoupling the redshift, Pinto et al. (2018) found a velocity of 1000±400 km s −1 in the Phoenix cluster. A similar velocity of ∼1000 km s −1 is found in the non cool core Coma cluster, while gas in the Perseus cluster is more relaxed at 480±210 km s −1 (e.g. Sanders et al. 2020). It is possible to drive cool gas by sloshing due to minor mergers (e.g. Ascasibar & Markevitch 2006). The shift in velocity is then seen in spatial coincidence with cold fronts (Sanders et al. 2020).
In A1664, the molecular gas system in the centre is divided in 2 roughly equal clumps with a velocity separation of 600 km s −1 (Russell et al. 2014). The blueshifted component is seen at a velocity of 571±7 km s −1 from CO(3-2) in our line-of-sight, with a FWHM of 190±20 km s −1 . Given the large uncertainty in the X-ray measurements, the velocities of the molecular and X-ray gas are consistent within 1 . Although it is unclear whether the blueshifted molecular gas lies in front or behind the BCG along the line of sight, the system is only a few kpc from the core in the transverse direction. This molecular gas is likely embedded in the soft X-ray gas cloud so these two gas phases may be related. The 0.8 keV gas has a sound speed of 460 km s −1 . The molecular gas will be shocked unless it comoves with the soft X-ray emitting gas.

Embedded multilayer cooling flow in RXCJ1504
There are larger amounts of cold obscuring material in the core of galaxy clusters, e.g. the Centaurus cluster (A3526) (Crawford et al. 2005b;. This suggests intrinsic absorption of the target galaxy is likely important and can reduce the amount of observed emission from a radiative cooling flow. Sanders et al. (2008) reported a factor of 3 larger cooling rate in the Centaurus cluster if there is a 4 × 10 21 cm −2 column density intrinsic to the cluster. This represents a significant amount of extra intrinsic luminosity available for powering emission due to cold gas.
In this section, we reintroduce a simple multilayer, intrinsically absorbed, cooling flow model, which was first proposed by Allen Figure 5. The schematic diagram of embedded cooling flow with three absorbing sheets and three cooling flow sheets. In the cluster, the blue columns represent sheets of absorbing gas and red columns represent sheets of radiative cooling flow. The black lines are the boundary between the columns. Blue arrows originated from the cluster represent emission from the associated sheet of cooling flow. & Fabian (1997). The schematic diagram of the model is shown in Fig. 5. For simplicity, we assume the cool core consists of several parallel sheets of material. Identical sheets of radiatively cooling gas in X-rays are placed in-between identical sheets of absorbing gas. The absorbing gas is assumed to be cold and neutral. An X-ray cooling sheet is absorbed by all absorbing sheets along the line-of-sight. This means the cooling sheet closest to the observer is absorbed once, and the furthest cooling sheet is absorbed three times in Fig. 5. The physical depth of these sheets is irrelevant in this work.
We assume each cooling gas sheet is emitting a flux . The fraction of the emitted energy transmitted through one sheet of absorbing gas is . We can write this transmission fraction as where − ( ) is the absorption cross-section and Δ H is the column density of one sheet of absorbing gas. The total observed flux is then where sheet is the number of sheets of absorbing gas components. Since absorption is a multiplicative process, it is possible to use the geometric series and equation 6 becomes This suggests only sheet and the total column density H,tot are the additional free parameters. Note that H,tot = sheet Δ H . In the large sheet limit, Equation 5 can be expanded and Equation 7 rewritten as Unfortunately, the geometric series implementation of absorption components is not yet available in SPEX. Nevertheless, it is possible to implement a brute-force model combining the existing spectral component of absorption and cooling flow. First, we select the number of absorbing sheets and choose a total column density. Each absorption component has the same column density of H,tot / sheet and is fixed in the spectral fitting. We assume the temperature of the absorbing gas is 0.5 eV, and the abundances are coupled to the X-ray cooling gas. The X-ray cooling gas is modelled by sheet cf components. Each of these cf components will be modified by a different number of absorption components before Galactic absorption. We couple the fit parameters of all cf components to one cf component. Then we measure the cooling rate from the cie temperature down to 0.01 keV as the complete cooling flow model described in section 3.1. The intrinsic absorption corrected cooling rate is then the total cooling rate of all cf components. This reconstructs the multilayer cooling flow model, which is equivalent to the complete cooling flow model with intrinsic absorption. We use a 15x15 grid in sheet and H,tot parameter space. We apply our model to RXCJ1504 and fit for minimal C-stat for each pair of sheet and H,tot . The improvement of C-stat from the complete cooling model is seen in Fig. 6, where we have included the contour at 68 per cent, 90 per cent and 95 per cent confidence levels. We search for the total column density that gives the minimum C-stat at any given number of sheets of absorbing gas. We find a valley of minimal C-stats in the parameter space. The absolute minimum of C-stat occurs for 1 absorption component with a column density of 6 × 10 21 cm −2 . However, the difference between the C-stats is less than 0.3 on this valley. This means sheet and H,tot are highly degenerate. Allen & Fabian (1997) found that the multilayer cooling flow model typically overpredicts the intrinsic column density by a factor of 1.5-3 in a sample of low redshift clusters. In conjunction with the valley of minimal C-stats in Fig. 6, this suggests that the true value of the intrinsic column density corresponds to a complex multilayer model with sheet ∼ 10.
We can compare the intrinsic absorption corrected cooling rate between the simplest 1 absorbing sheet model and the 10 sheets model. For the 1 sheet model, we measure a cooling rate of 430±90 M yr −1 . This is 8 times higher than the complete cooling rate without intrinsic absorption (see Model 2 in Table 2). For the 10 sheet model with a total column density of 1.5 × 10 22 cm −2 , the cooling rate is 520±30 M yr −1 . The intrinsic absorption corrected cooling rate is consistent between 1 sheet and 10 sheets model at the 1 level.
For an order of magnitude increase in the cooling rate above 0.7 keV in RXCJ1504, it can contribute ∼40 per cent of the energy required to power the UV/optical line-emitting nebula, but not all. In other massive clusters with neb > 100M yr −1 (such as 2A0335+096, A1835 and A2597, see Liu et al. 2019), 10 times the cooling rate means the soft X-ray emitting gas can power the UV/optical nebula alone.
We also apply the embedded cooling flow model to A1664. For 10 sheets of absorbing gas with a total column density of 3 × 10 21 cm −2 , we only measure a cooling rate of 31±7 M yr −1 . This is consistent with the cooling rate without intrinsic absorption. No significant change of the cooling rate is detected for other combinations of sheet and H,tot . Note that the effect of embedded absorption on the optical line-emitting gas is explored by Polles et al. 2021.
We also need to reexamine the role of absorption corrected cooling rate in the AGN feedback. It is 24 per cent of the predicted rate in the absence of heating. This only slightly reduces the amount of heating from the AGN and the black hole growth rate. On the other hand, the cooling rate is 3.3 times higher than the star formation rate. The difference between the cooling and star formation rates is 300 M yr −1 . This suggests the molecular gas reservoir is growing 6 times faster than the unabsorbed cooling model with a formation timescale less than 10 8 yr.
In future work, it will be interesting to reduce the degeneracy between intrinsic column density and the multilayer structure of embedded cooling flow model. Consideration will be needed for scattering of resonance lines (see studies of cool X-ray emitting gas in groups and elliptical galaxies by Pinto et al. 2016, Ogorzalek et al. 2017. Such lines can be absorbed by the cold gas as they scatter around in the plasma. The observational situation will be improved with the future high-spectral-resolution mission XRISM (XRISM Science Team 2020), with its non-dispersive calorimeter, and later by the X-IFU of Athena (Barret et al. 2018).

CONCLUSIONS
We have performed a multiphase cooling flow analysis on deep XMM-Newton RGS observations of two X-ray luminous cool core clusters RXCJ1504 at = 0.2153 and A1664 at = 01283. The cooling rate is measured to be 180±40M yr −1 and 34±6M yr −1 for RXCJ1504.1-0248 and A1664, respectively. It is higher than the observed star formation rate in both clusters. We detect an upper limit of residual cooling rate below 0.7 keV of 53M yr −1 at 90 per cent confidence level in RXCJ1504.1-0248. The energy of the cooling gas is insufficient to power the UV/optical line-emitting nebula in both clusters and additional sources of energy are required. If the molecular gas reservoir is accumulating mass from the condensation of the radiatively cooling gas, the formation timescale is 1-3×10 8 yr from the observed cooling rate but is likely longer due to the high star formation activities.
We also place a tight constraint on turbulence in the core. An upper limit of 300 km s −1 of 1-D turbulent velocity at 90 per cent confidence level is measured in both clusters. These velocities correspond to a Mach number of 0.23 and 0.4 for RXCJ1504.1-0248 and A1664, Figure 6. The improvement of C-stat over the complete cooling model without intrinsic absorption. The yellow curve represents the maximum improvement of C-stat at each number of sheets of absorbing gas.
respectively. The energy density of turbulence is equivalent to 8.9 per cent and 27 per cent of the thermal energy density, which is inadequate to fully transfer AGN heating throughout the cooling core. We find the cool component of 0.80±0.08 keV is blueshifted from the systemic velocity of the cluster at 750 +800 −280 km s −1 in A1664. This is consistent with the velocity of the blueshifted component in the molecular gas, but we cannot rule out an origin within a sloshing cold front for the blueshifted X-ray gas.
We reintroduce a multilayer, intrinsically-absorbed, cooling flow model. In RXCJ1504.1-0248, we find that the cooling rate increases to 520±30 M yr −1 using the 10 absorbing sheet model. This is an order of magnitude higher than the cooling rate measured without intrinsic absorption. The intrinsically absorbed cooling rate of A1664 is unaffected and consistent with the current measurement.
In the future, XRISM and Athena will help to unveil the connection between molecular and X-ray emitting gas phases and determine the influence of intrinsic absorption on cooling flows.