# Exploring methods for creation of Boron-vacancies in hexagonal Boron Nitride
exfoliated from bulk crystal

T. Zabelotsky  The Center for Nanoscience and Nanotechnology, The Hebrew
University of Jerusalem, Jerusalem 91904, Israel  Dept. of Applied Physics,
Rachel and Selim School of Engineering, Hebrew University, Jerusalem 91904,
Israel  S. Singh  The Racah Institute of Physics, The Hebrew University of
Jerusalem, Jerusalem 91904, Israel  G. Haim  Dept. of Applied Physics, Rachel
and Selim School of Engineering, Hebrew University, Jerusalem 91904, Israel
School of Physics, The University of Melbourne, Parkville, Victoria 3010,
Australia  R. Malkinson  The Center for Nanoscience and Nanotechnology, The
Hebrew University of Jerusalem, Jerusalem 91904, Israel  Dept. of Applied
Physics, Rachel and Selim School of Engineering, Hebrew University, Jerusalem
91904, Israel  S. Kadkhodazadeh  DTU Nanolab, Technical University of Denmark,
Fysikvej, Kongens Lyngby 2800, Denmark  I. P. Radko  Department of Physics,
Technical University of Denmark, Kongens Lyngby 2800, Denmark  I. Aharonovich
School of Mathematical and Physical Sciences, University of Technology Sydney,
Ultimo, New South Wales 2007, Australia  ARC Centre of Excellence for
Transformative Meta-Optical Systems (TMOS), Faculty of Science, University of
Technology Sydney, Australia  H. Steinberg  The Center for Nanoscience and
Nanotechnology, The Hebrew University of Jerusalem, Jerusalem 91904, Israel
The Racah Institute of Physics, The Hebrew University of Jerusalem, Jerusalem
91904, Israel  Kirstine Berg-Sørensen  Department of Health Technology,
Technical University of Denmark, Kongens Lyngby 2800 Denmark  A. Huck  Center
for Macroscopic Quantum States (bigQ), Department of Physics, Technical
University of Denmark, 2800 Kongens Lyngby, Denmark  T. Taniguchi
International Center for Materials Nanoarchitectonics, National Institute for
Materials Science, 1-1 Namiki, Tsukuba 305-0044, Japan  K. Watanabe  Research
Center for Functional Materials, National Institute for Materials Science, 1-1
Namiki, Tsukuba 305-0044, Japan  N. Bar-Gill  The Center for Nanoscience and
Nanotechnology, The Hebrew University of Jerusalem, Jerusalem 91904, Israel
Dept. of Applied Physics, Rachel and Selim School of Engineering, Hebrew
University, Jerusalem 91904, Israel  The Racah Institute of Physics, The
Hebrew University of Jerusalem, Jerusalem 91904, Israel

###### Abstract

Boron vacancies (VB-) in hexagonal boron-nitride (hBN) have sparked great
interest in recent years, due to their electronic spin properties. Since hBN
can be readily integrated into devices where it interfaces a huge variety of
other 2D materials, boron vacancies may serve as a precise sensor which can be
deployed at very close proximity to many important materials systems. Boron
vacancy defects may be produced by a number of existing methods, the use of
which may depend on the final application. Any method should reproducibly
generate defects with controlled density and desired pattern. To date,
however, detailed studies of such methods are missing. In this paper we study
various techniques, focused ion beam (FIB), electron irradiation and ion
implantation, for the preparation of hBN flakes from bulk crystals, and
relevant post-processing treatments to create VB-s as a function of flake
thickness and defect concentrations. We find that flake thickness plays an
important role when optimising implantation parameters, while careful sample
cleaning proved important to achieve best results.

##  I Introduction

Hexagonal boron-nitride (hBN) is a graphite-analog layered material with a
hexagonal lattice where Boron and Nitrogen atoms alternate in each hexagon.
hBN is a wide band-gap insulator, which became an essencial ingredient in
stacked van der Waals (vdW) heterostructures[1]. hBN was first used in vdW
heterostructures as an ultra-flat, electrically stable substrate [2], and as a
means to fully encapsulate graphene [3]. It was soon realized that ultrathin
hBN flakes can be used as tunnel barriers [4], and that hBN capping can
stabilize air-sensitive materials [5]. Since hBN has strong adhesion to
graphene, it has enabled the use of selective tearing and stacking, the key
method used in the assembly of twist-angle-controlled devices [6].

Defects in hBN have been studied using electronic transport - where they are
identified as nanoscale quantum dots, sensitive to the spectra and
compressibility of nearby samples [7, 8]. They have also been extensively
studied using optical fluorescence. Color centers in hBN, specifically VB-s,
are stable spin defects in the hBN lattice with similar traits to the well
known Nitrogen Vacancies (NV-) in diamond [9, 10, 11, 12] and have comparable
applications [13].

Given the potential of such optically-active spin defects in vdW materials,
significant effort has been invested over the past few years in the controlled
fabrication and characterization of these color centers [14, 15, 16].
Nevertheless, a systematic study of VB- creation as a function of method
(electron irradiation vs. ion implantation), preparation protocol, and flake
thickness is still missing.

In this work we have tested different methods for generating VB-s in bulk
exfoliated hBN, while focusing on flake preparation and the effect of flake
thickness on the yield of VB-s. We describe in detail the preparation method
applied to the flakes, the methods applied for VB creation, the
characterization of the flakes by fluorescence spectroscopy and spin resonance
measurements. We present successful creation of VB-s using 12 keV ion
implantation through FIB, with both Nitrogen and Oxygen ion beams. The defects
were identified by their spectral signature with a peak emission at $\sim 800$
nm and their spin resonance at $\leavevmode\nobreak\ 3.47$ GHz at zero field
[17, 18]. We observe increasing fluorescence with increasing ion fluence, and
identify a thickness dependence, such that thicker flakes show higher
fluorescence intensities for the same fluence. Stopping Range of Ions in
Matter (SRIM) simulations were performed to better understand the defect
creation process. Our results indicate that in practice, vacancies are created
in the sample further than the ion stopping range, perhaps due to ion
channeling, which is not taken into account during the simulations. This
extends the interaction depth of the ion beam with the samples, resulting in a
strong thickness dependent fluorescence of the flakes, as observed. In
contrast, our attempts at generating VB- using electron beam in transmission
electron microscope (TEM) or electron lithography were not effective.

##  II Preparing hBN flakes

One of the simplest techniques for preparing hBN flakes is exfoliation from
bulk crystal (using blue tape, Nitto), similar to the processes originally
developed for graphene and molybdenum disulphide (MoS2) [1]. It must be noted,
that for our purposes the exfoliation should be handled carefully in a clean
environment, as it could introduce unwanted carbon contaminants to the flakes
[19], which could then impede the creation of VB-s, introduce other defects to
the material, and lead to decreased Signal to Noise Ratio (SNR) for VB- spin
readout.

Once exfoliated using the tape, flakes were transferred to a Silicon (Si)
substrate with a thin layer (285 nm) of Silicon dioxide (SiO2) and gold
markers. To clean the flakes from contaminants that were introduced during
exfoliation, the samples were heated in forming gas (95% Ar and 5% H) for 4
hours at 375∘C. We performed this cleaning process on flakes marked N1, N2, O1
and O2 (see Table 1 below). Additionally, another sample, flake N3, was
cleaned using an open air hot plate, for 15 minutes at 400∘C, the results of
the characterization for this flake are detailed in the supplementary material
[20].

##  III Defect creation

###  III.1 FIB treatment

Ion implantation [14], was the focus of the work presented here, and was
implemented on flakes N1, N2, N3, O1 and O2, summarized in Table 1. This
approach, which resulted in our most successful VB- yield, was based on using
a focused ion beam (FIB, Helios 5 Hydra Dualbeam plasma FIB/SEM instrument),
utilizing Nitrogen (N) and Oxygen (O) ion elements at 12 keV.

Flake | Flake | FIB | Cleaning | Color
---|---|---|---|---
name | thickness | fluence | method | in figures
| [nm] | $[{\frac{ion}{cm^{2}}]}$ |  |
N1 | 100 / 150 | $6.25\cdot 10^{14}$ | forming gas | Blue / Orange
N2 | 73 | dose test | forming gas | Yellow
N3 | N.A. | dose test | hot plate | N.A.
O1 | 114 | $1.8\cdot 10^{13}$ | forming gas | Red
O2 | 83 | dose test | forming gas | Green
Table 1: FIB testament description for different flakes discussed in this
manuscript.

Dose tests were performed with both Oxygen and Nitrogen ion beams. In both
cases the parameters used, ion beam energy, beam size, dwell time and beam
current, were kept constant, while changing the number of exposure repetitions
to achieve varying dosages. The VB- yield is not sensitive to changing these
parameters, as long as the beam fluence is kept unchanged.

The fluence was calculated in the following way:

$\large Dose=\frac{I_{beam}t_{dwell}N_{repeats}}{qA_{beam}},$ (1)

where $I_{beam}$ is the beam current, $t_{dwell}$ is the dwell time,
$N_{repeats}$ is the number of repetitions, $q$ is an electron’s charge and
$A_{beam}$ is the beam size area.

Flakes N1 and O1 underwent a uniform FIB exposure, with fluences of
$6.25{\cdot}10^{14}$ $cm^{-2}$ and $1.81{\cdot}10^{13}$ $cm^{-2}$
respectively, while the flakes N2, N3 and O2 underwent dose tests in an array
of square regions of size $3{\times}3\ {\mu}m^{2}$.

We note that increasing the fluence further (see supplementary material [20]),
to the range of $\sim 10^{16}$ $cm^{-2}$, creates visible damage on the flakes
(observed optically as a change in the flake’s color), although VB-s were
created. Increasing the dose even further, up to $10^{17}-10^{18}$ $cm^{-2}$
could sputter away layers without creating any VB-s leading to unwanted damage
to the flakes, although we note that there are other factors that should be
taken into account in this context such as flake thickness and beam energy.

###  III.2 Other methods

In addition to FIB, we attempted ion implantation using standard, commercial
(semiconductor industry) implanter instrumentation (Innovion Corp.) to
reproduce previously reported results [16]. In our case these trials were not
successful, resulting in measured fluorescence but not with characteristics
matching those of VB-s. We attribute the measured fluorescence to other
defects. Nevertheless, we note that for this process no cleaning was performed
on the flakes post transfer to the Si substrate (before implantation). Since
unclean flakes exhibited limited defect creation also using FIB, we deduce
that sample cleaning before implantation plays an important role.

Furthermore, we studied electron irradiation instead of ion implantation,
performed using two electron sources: e-beam lithography (100keV, Elionix
ELS-G100) and TEM (80keV and 300keV, Themis-Z). The hBN samples used for these
electron irradiation studies went through the same sample preparation as
detailed for the FIB. For TEM irradiation, the flakes were transferred to TEM
Gold grids (200 Mesh) with carbon film (R 1.2/1.3) using the following
procedure: first, hBN flakes were transferred to a Si substrate spin coated
with Poly(vinyl alcohol) (PVA), then the TEM grid’s carbon side was placed on
the Si, and finally a drop of water was added to dissolve the PVA, allowing
the hBN flakes to be transferred to the grid.

Electron irradiation under the above conditions did not result in successful
creation of VB-s. We attribute this to insufficient energy of the electron
beam. Theoretically, the maximum energy that 80 keV and 100 keV electrons can
transfer to B atoms are approximately 17.5 eV and 22.2 eV, respectively [21,
22]. Given the displacement threshold of 19.36 eV reported for B in pristine
hBN [23], it is unlikely to create VB-s in hBN flakes using 80 keV electrons.
Although the energy transfer from 100 keV electrons should be sufficient for
VB-s creation, the probability of such a process is relatively low and likely
requires very high electron beam doses [23].

##  IV Sample characterization

###  IV.1 Optical spectroscopy

Figure 1: Confocal setup sketch for spectral measurements, 532 nm laser is
used for excitation through an Olympus air objective M Plan Apo 100x/0.95. PL
is collected through a dichroic mirror either to a single photon counter
(Perkin-Elmer SPCM, model AQRH-14) for confocal scanning, or to a spectrometer
(Andor SR500i).

Given the well known spectral signature of VB- defects [17], optical
spectroscopy is an efficient and convenient tool for characterizing defect
creation.

Emission spectra of the defects was measured using a home made confocal setup
(Figure 1) equipped with a spectrometer (Andor SR500i).

Figure 2: Flake N1. (a) Confocal scan of a flake treated by N FIB implantation
(Green dashed line), the dots depict the areas for which the spectral
measurement was performed. (b) Microscopic image of the flake. (c) Comparison
of the spectral signature of two thicknesses of the hBN - 150 nm (Orange), 100
nm (Blue) and the silicone substrate for reference (Pink). (d) SRIM simulation
showing the N ion depth distribution for energy of 12 keV, $\mu=33.3\ nm$,
$\sigma=12.1\ nm$. (e-f) AFM scan of the flake showing it’s thickness in the
relevant areas.

A uniform fluence of Nitrogen FIB treatment was applied on flake N1 [as seen
in Figure 2(a)], which consists of regions with varying thicknesses [100-150
nm as measured by atomic force microscope (AFM), see Figures 2(e-f)]. FIB
treatment was carried out using the following parameters:

  * •

beam size 0.01${\mu}m^{2}$,

  * •

dwell time 1${\mu}s$,

  * •

current of 10$pA$,

  * •

exposure repetitions 1k.

These parameters translate to a dose of $6.25{\cdot}10^{14}$ $cm^{-2}$.
Spectral measurements confirm the creation of VB-s in flake N1 [Figure 2(c)].

The spectral measurements depicted in the figures were smoothed out using
triangle filters for visual purposes. Additionally, the spectral raw data were
fitted to gaussians in the relevant areas, to estimate the SNR from the VB-s,
further details of the SNR calculation can be found in the supplementary
materials [20]). The strongest PL signal and best SNR were observed from the
thicker area of the flake (150 nm, SNR $\sim$20), compared to a thinner area
(100 nm, SNR $\sim$9) [see Figure 2(c)]. While some area of the Si substrate
was exposed as well, and exhibited some photo-luminescence (PL), it was
clearly unrelated to the relevant hBN defects, as can be observed from the
spectral signature [Figure 2(c)].

Figure 3(a) depicts dose test experiments performed with Oxygen on flake O2
which is 83 nm thick, with the following parameters: energy 12 keV, beam size
0.04 ${\mu}m^{2}$, dwell times of 1 ${\mu}s$ and 10 ${\mu}s$, current of 30
$pA$ and varying repetitions. The PL intensity increases with dose, while
different dwell times don’t have a significant impact on PL intensity for the
same dose.

Figure 3: Flake O2. (a) Oxygen FIB dose test(see colorbar) spectral
measurement results. (b) SNR of the VBs spectral signal for different dosages,
for a given thickness of 83 nm and implantation energy (12 keV), the SNR grows
with the dose, starting to approach saturation, though not demonstrated.

###  IV.2 Thickness Dependence

SRIM simulations [Figures 4(b), 4(d)] predict the maximal ion penetration
depth to be approximately 70 nm (for our parameters). Therefore it is not
a-priori clear why for a given ion beam fluence, using 12 keV, a flake thicker
than 70 nm, would show stronger PL and better SNR. Therefore we conjecture
that this result suggests that mainly ion channeling or secondary ions create
the VB- defects, rather than first collision ions. Similar behavior have been
observed for NV-s in diamonds [24], which highlights the limitation of SRIM to
predict the depth of defect created in the material.

For the Nitrogen FIB we see a clear difference in PL intensity [Figure 4(a)]
and in calculated SNR, as a function of flake thickness [Figure 4(b)].
Similarly, Oxygen FIB implantation spectral measurements [Figure 4(c)]
indicate that the thicker flake O1 has greater SNR of about 4 while for the
thinner flake, O2, the SNR was less than 0.5, under similar FIB fluence.

Figure 4 shows that for the thinner (83 nm) flake, O2, a PL intensity of $\sim
250$ in the arbitrary, but fixed, intensity units, requires a higher dosage to
match the same PL that was measured for the thicker (114 nm) flake, O1. From
Figure 3(b) it can be seen that an increase of one order of magnitude in the
fluence for the thinner flake is needed to match the same SNR of the thicker
flake.

These results confirm the argument that a simple SRIM range simulation does
not predict the depth of VB-s in hBN, and potentially indicates the non-
negligible role of secondary ions in creating these defects. Practically, to
achieve a specific VB- concentration in a given hBN flake, for a certain beam
energy, the implantation dose should be calibrated with respect to the flake’s
thickness, since thicker flakes require a lower ion dose, compared to thinner
ones.

Figure 4: Flake N1 (Blue and Orange) and flake N2 (Yellow), flake O1 (Red) and
flake O2 (Green). (a) spectral measurement comparison for similar doses of
Nitrogen FIB ${\sim}5{\cdot}10^{14}cm^{-2}$, the thicker the flake the higher
the PL intensity. (b) SRIM simulation for Nitrogen at 12 keV alongside the
different flake thicknesses. (c) Spectral measurement of hBN flakes after
Oxygen FIB treatment with similar dose ${\sim}2{\cdot}10^{13}cm^{-2}$, the
thicker flake, O1, has significantly stronger PL intensity than the thinner
flake, O2. (d) SRIM depth range simulation for O ions on hBN, at energy of 12
keV, $\mu=29.9\ nm$, $\sigma=11.1\ nm$. (e) thickness dependence on the
spectral signal of the VB-s (divided by the ion fluence), where the thicker
flakes show stronger signal for both ion species Nitrogen (dashed Pink) and
Oxygen (dashed Black).

###  IV.3 Optically detected magnetic resonance

The VB-, like the NV-, is a spin one system with a triplet ground state. Its
resulting zero field splitting is $\sim$3.47 GHz [17, 18]. As such, and given
its optical response, an Optically Detected Magnetic Resonance (ODMR)
measurement can be performed to observe and characterize the VB- spin
properties [25].

The pulsed ODMR sequence consists of an optical pumping pulse to initialize
the spins to their $m_{s}=0$ state, a microwave (MW) pulse applied with
varying frequencies, followed by a final optical pulse for spin readout. At
the spin resonant MW frequency spin transitions are induced between the
$m_{s}=0$ and $m_{s}=\pm 1$ ground states, such that the resulting
fluorescence is suppressed.

Figure 5: Measured results for flake N1, ODMR with a bias field, showing the
Zeeman shift of $m_{s}=\pm 1$, this is consistent with the expected
$\gamma=2.8\ \frac{MHz}{G}$ of a free electron. (a) no Bias field. (b) Bias
field of 54 G and (c) 77 G.

To perform the ODMR we used an $\Omega$-shaped co-planar antenna, ensuring
that the flake is as close as possible to the conducting loop to achieve
higher effective MW power (and thus stronger driving). Zero field ODMR
measurements on the nitrogen FIB implanted flake are presented in figure 5(a).
ODMR with an external magnetic field at 54 Gauss and 77 Gauss are presented in
figure 5(b,c), demonstrating the linear Zeeman split between the $m_{s}=\pm 1$
ground spin states.

Similar results were achieved on Oxygen FIB flakes (see supplementary
material), exhibiting ODMR contrast variations with PL intensity on the
different flakes, related to different VB- concentrations.

##  V Summary

In summary, we have presented a study of the application of Nitrogen and
Oxygen FIB implantation for creating VB- defects in hBN flakes. The field of
atomic defects in 2D materials is growing rapidly, with clear prospects for
potential impact on various quantum technologies, including sensing and
photonics. Nevertheless, a systematic study of the processes and techniques
for the robust creation of desired defects is missing, and is the focus of
this manuscript. Our results show that both Oxygen and Nitrogen ions at 12 keV
can efficiently create VB-s in hBN flakes. We find that hBN flake thickness
plays a key role in the success of the process and in the achievable VB-
density, and for a given beam energy, a higher dose is required for thinner
flakes. Furthermore, we find that the cleaning process used on the flakes
prior to implantation is very important. Specifically, the inclusion of a
forming gas heating step or a simple hot plate cleaning proved crucial for
successful creation of VB- in hBN.

We also note that in certain cases flakes could be created by well-controlled,
clean methods not related to exfoliation [26, 27]. Nevertheless, exfoliation
from bulk crystal is a widely used technique, with the advantage of providing
access to high quality material. As such, the techniques presented here, and
specifically the cleaning processes, should be carefully considered prior to
implantation.

##  VI Acknowledgements

N.B. acknowledges support from the European Union’s Horizon 2020 research and
innovation program under Grant Agreements No. 101070546 (MUQUABIS) and No.
828946 (PATHOS), and has been supported in part by CIFAR, the Ministry of
Science and Technology, Israel, the innovation authority (Project No. 70033),
and the ISF (Grants No. 1380/21 and No. 3597/21). G.H. acknowledges support by
the Melbourne research scholarship. I.A. acknowledges financial support from
the Australian Research Council (CE200100010, FT220100053) and the Office of
Naval Research Global (N62909-22-1-2028). We acknowledge financial support
from the Danish Agency for Higher Education and Science, grant no.
9096-00023B.

## References

  * Geim and Grigorieva [2013] A. K. Geim and I. V. Grigorieva, Nature 499, 419 (2013).
  * Dean  _et al._ [2010] C. R. Dean, A. F. Young, I. Meric, C. Lee, L. Wang, S. Sorgenfrei, K. Watanabe, T. Taniguchi, P. Kim, K. L. Shepard, and J. Hone, Nature Nanotechnology 5, 722 (2010).
  * Mayorov  _et al._ [2011] A. S. Mayorov, R. V. Gorbachev, S. V. Morozov, L. Britnell, R. Jalil, L. A. Ponomarenko, P. Blake, K. S. Novoselov, K. Watanabe, T. Taniguchi, and A. K. Geim, Nano Letters 11, 2396 (2011), pMID: 21574627, https://doi.org/10.1021/nl200758b .
  * Amet  _et al._ [2012] F. Amet, J. R. Williams, A. G. F. Garcia, M. Yankowitz, K. Watanabe, T. Taniguchi, and D. Goldhaber-Gordon, Phys. Rev. B 85, 073405 (2012).
  * Cao  _et al._ [2015] Y. Cao, A. Mishchenko, G. L. Yu, E. Khestanova, A. P. Rooney, E. Prestat, A. V. Kretinin, P. Blake, M. B. Shalom, C. Woods, J. Chapman, G. Balakrishnan, I. V. Grigorieva, K. S. Novoselov, B. A. Piot, M. Potemski, K. Watanabe, T. Taniguchi, S. J. Haigh, A. K. Geim, and R. V. Gorbachev, Nano Letters 15, 4914 (2015), pMID: 26132110, https://doi.org/10.1021/acs.nanolett.5b00648 .
  * Cao  _et al._ [2018] Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero, Nature 556, 43 (2018).
  * Keren  _et al._ [2020] I. Keren, T. Dvir, A. Zalic, A. Iluz, D. LeBoeuf, K. Watanabe, T. Taniguchi, and H. Steinberg, Nature Communications 11, 3408 (2020).
  * Greenaway  _et al._ [2018] M. T. Greenaway, E. E. Vdovin, D. Ghazaryan, A. Misra, A. Mishchenko, Y. Cao, Z. Wang, J. R. Wallbank, M. Holwill, Y. Khanin, S. V. Morozov, K. Watanabe, T. Taniguchi, O. Makarovsky, T. M. Fromhold, A. Patanè, A. K. Geim, V. I. Fal’ko, K. S. Novoselov, and L. Eaves, Communications Physics 1, 94 (2018).
  * Genov  _et al._ [2020] G. T. Genov, Y. Ben-Shalom, F. Jelezko, A. Retzker, and N. Bar-Gill, Phys. Rev. Res. 2, 23 (2020).
  * Taylor  _et al._ [2008] J. M. Taylor, P. Cappellaro, L. Childress, L. Jiang, D. Budker, P. R. Hemmer, A. Yacoby, R. Walsworth, and M. D. Lukin, Nature Physics 4, 810 (2008).
  * Dolde  _et al._ [2011] F. Dolde, H. Fedder, M. W. Doherty, T. Nöbauer, F. Rempp, G. Balasubramanian, T. Wolf, F. Reinhard, L. C. L. Hollenberg, F. Jelezko, and J. Wrachtrup, Nature Physics 7, 459 (2011).
  * Pham  _et al._ [2011] L. M. Pham, D. L. Sage, P. L. Stanwix, T. K. Yeung, D. Glenn, A. Trifonov, P. Cappellaro, P. R. Hemmer, M. D. Lukin, H. Park, A. Yacoby, and R. L. Walsworth, New Journal of Physics 13, 045021 (2011).
  * Castelletto and Boretti [2021] S. Castelletto and A. Boretti, Advanced Photonics 3, 054001 (2021).
  * Kianinia  _et al._ [2020] M. Kianinia, S. White, J. E. Fröch, C. Bradac, and I. Aharonovich, ACS Photonics 7, 2147 (2020), https://doi.org/10.1021/acsphotonics.0c00614 .
  * Murzakhanov  _et al._ [2021] F. F. Murzakhanov, B. V. Yavkin, G. V. Mamin, S. B. Orlinskii, I. E. Mumdzhi, I. N. Gracheva, B. F. Gabbasov, A. N. Smirnov, V. Y. Davydov, and V. A. Soltamov, Nanomaterials 11, 10.3390/nano11061373 (2021).
  * Guo  _et al._ [2022] N.-J. Guo, W. Liu, Z.-P. Li, Y.-Z. Yang, S. Yu, Y. Meng, Z.-A. Wang, X.-D. Zeng, F.-F. Yan, Q. Li, J.-F. Wang, J.-S. Xu, Y.-T. Wang, J.-S. Tang, C.-F. Li, and G.-C. Guo, ACS Omega 7, 1733 (2022), https://doi.org/10.1021/acsomega.1c04564 .
  * Gottscholl  _et al._ [2020] A. Gottscholl, M. Kianinia, V. Soltamov, S. Orlinskii, G. Mamin, C. Bradac, C. Kasper, K. Krambrock, A. Sperlich, M. Toth, I. Aharonovich, and V. Dyakonov, Nature Materials 19, 540 (2020).
  * Ivády  _et al._ [2020] V. Ivády, G. Barcza, G. Thiering, S. Li, H. Hamdi, J.-P. Chou, Ö. Legeza, and A. Gali, npj Computational Materials 6, 41 (2020).
  * Huang  _et al._ [2020] Y. Huang, Y.-H. Pan, R. Yang, L.-H. Bao, L. Meng, H.-L. Luo, Y.-Q. Cai, G.-D. Liu, W.-J. Zhao, Z. Zhou, L.-M. Wu, Z.-L. Zhu, M. Huang, L.-W. Liu, L. Liu, P. Cheng, K.-H. Wu, S.-B. Tian, C.-Z. Gu, Y.-G. Shi, Y.-F. Guo, Z. G. Cheng, J.-P. Hu, L. Zhao, G.-H. Yang, E. Sutter, P. Sutter, Y.-L. Wang, W. Ji, X.-J. Zhou, and H.-J. Gao, Nature Communications 11, 2453 (2020).
  * [20] S. Material, .
  * Egerton [2011] R. Egerton,  _Electron Energy-Loss Spectroscopy in the Electron Microscope; 3rd ed_ (Springer, 2011).
  * Egerton  _et al._ [2004] R. Egerton, P. Li, and M. Malac, Micron 35, 399 (2004).
  * Kotakoski  _et al._ [2010] J. Kotakoski, C. H. Jin, O. Lehtinen, K. Suenaga, and A. V. Krasheninnikov, Phys. Rev. B 82, 113404 (2010).
  * Toyli  _et al._ [2010] D. M. Toyli, C. D. Weis, G. D. Fuchs, T. Schenkel, and D. D. Awschalom, Nano Letters 10, 3168 (2010), pMID: 20698632, https://doi.org/10.1021/nl102066q .
  * Baber  _et al._ [2022] S. Baber, R. N. E. Malein, P. Khatri, P. S. Keatley, S. Guo, F. Withers, A. J. Ramsay, and I. J. Luxmoore, Nano Letters 22, 461 (2022), pMID: 34958574, https://doi.org/10.1021/acs.nanolett.1c04366 .
  * Gao  _et al._ [2019] M. Gao, J. Meng, Y. Chen, S. Ye, Y. Wang, C. Ding, Y. Li, Z. Yin, X. Zeng, J. You, P. Jin, and X. Zhang, J. Mater. Chem. C 7, 14999 (2019).
  * Jana and Singh [2018] M. Jana and R. N. Singh, International Materials Reviews 63, 162 (2018), https://doi.org/10.1080/09506608.2017.1322833 .

##  VII APPENDIX A - hot plate

As mentioned in the manuscript, some samples were cleaned using simple hot
plate, set to $400^{\circ}C$ for 15 minutes.

Figure 6: Flake N3 dose test, varying fluences in the range of $10^{13}\sim
10^{14}cm^{-2}$

In Figure 6 the spectral signature of the FIB exposed areas, suggest that VB-s
were created in flake N3. The hot plate cleaning process did not eliminate all
contaminants from the flakes, thus some fluorescence is visible for
wavelengths below $700\ nm$, similarly to flake that went though annealing in
forming gas. We propose that the contaminants might have been introduced
during a SEM imaging session, though not thoroughly tested.

##  VIII APPENDIX B - data smoothing and SNR

For clear visibility of the spectral measurements, the raw data was smoothed
out using a triangular filter as described in the following manner:

$\large
S(\lambda)={\frac{\displaystyle\sum_{-\lfloor{\frac{w}{2}}\rfloor}^{\lfloor{\frac{w}{2}}\rfloor}({\lceil{\frac{w}{2}}\rceil}-{\lvert{m}\rvert})Y(\lambda+m\cdot\Delta\lambda)}{w+2\cdot\sum(w-j)}}$
(2)

Where $S(\lambda)$ is the smoothed out signal, $\Delta\lambda$ the wavelength
increments, $Y(\lambda)$ the raw signal and $w$ is the smoothing size of the
filter. Then to estimate an SNR value, a gaussian shape was extracted from the
following fit on the raw data of the spectral measurement (as can be seen in
Figure 7) using the following formula:

$\large
fit={A\cdot\exp{(-\frac{(\lambda-\mu)^{2}}{2\cdot\sigma^{2}})}+m\cdot{x}+C}$
(3)

While excluding data point below $650\ nm$, since a low pass filter of that
wavelength was used during acquisition.

Figure 7: Spectral measurement from flake N1 with a fit to the following
model: $A\cdot\exp{(-\frac{(\lambda-\mu)^{2}}{2\cdot\sigma^{2}})}+m\cdot{x}+C$

Finally for the SNR calculation, the ration between the area under the
gaussian and root mean square error ($\sigma_{RMSE}$) of the fit, as shown
here:

$\large
SNR={\frac{A\cdot\int{\exp{(-\frac{(\lambda-\mu)^{2}}{2\cdot\sigma^{2}})}}d\lambda}{\sigma_{RMSE}^{2}}}$
(4)

##  IX APPENDIX C - electron irradiation

There were few attempts to create VB-s using electron irradiation with
different energies, using two different techniques, e-beam (at 100 keV) and
TEM (at 300 keV), the examples shown in Figures 8 and 9, these attempts
weren’t successful as can be seen in the spectral measurements(Figure 8.c and
9.c) regardless of annealing treatment.

Figure 8: Two flakes that underwent e-beam dose test ranging from
$10^{17}cm^{-2}$ to $10^{19}cm^{-2}$ at 100 keV. (a) microscope image of the
flakes (b) confocal scan of the fluorescence from the flakes after the
irradiation. (c) spectral wavelength of exposed areas on the flakes.  Figure
9: Two flakes that underwent TEM dose test ranging from $10^{19}cm^{-2}$ to
$10^{22}cm^{-2}$ at 300 keV. (a) microscope image of the flakes (b) confocal
scan of the fluorescence from the flakes after the irradiation. (c) spectral
wavelength of exposed areas on the flakes.

##  X APPENDIX D - ion implantation

There was an attempt to create VB-s by ion implantation using commercial
(semiconductor industry) ion implanter instrumentation (Innovion Corp.). This
was unsuccessful as can be seen from Figure 10.b, but further note that in
this case the flakes were not properly cleaned prior to ion treatment.

Figure 10: Two flakes that underwent uniform fluence of ion beam
$10^{14}cm^{-2}$ at 40 keV. (a) confocal scan. (b) spectral measurement on
flake

##  XI APPENDIX E - ODMR

Figure 11: Flake O1, ODMR measurement without bias magnetic field


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

Haonan Liu1 Andrew C. Fabian1, Ciro Pinto1,2, Helen R. Russell1,3, Jeremy S.
Sanders4 and Brian R. McNamara5,6
1Institute of Astronomy, Madingley Road, CB3 0HA Cambridge, United Kingdom
2INAF-IASF Palermo, Via U. La Malfa 153, I-90146 Palermo, Italy
3School of Physics, Astronomy, University of Nottingham, University Park,
Nottingham NG7 2RD, UK
4Max-Planck-Institut für extraterrestrische Physik, Gießenbachstraße 1, 85748,
Garching, Germany
5Department of Physics and Astronomy, University of Waterloo, 200 University
Avenue West, Waterloo, ON N2L 3G1, Canada
6Perimeter Institute for Theoretical Physics, Waterloo, ON N2L 2Y5, Canada
E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

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.

###### keywords:

X-rays: galaxies: clusters - galaxies: clusters: general

††pubyear: 2020††pagerange: Suppressed cooling and turbulent heating in the
core of X-ray luminous clusters RXCJ1504.1-0248 and Abell 1664–Suppressed
cooling and turbulent heating in the core of X-ray luminous clusters
RXCJ1504.1-0248 and Abell 1664

##  1 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\rm\,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. 2012, 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 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$\pm$10 km $\rm 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. 2010; Sanders et al.
2011; Sanders & 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, $\Delta\lambda=0.138\Delta\theta/m$ Å, where $\Delta\lambda$
is the broadening in wavelength, $\Delta\theta$ is the angular offset from the
central source in arcmin and $m$ 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 $\rm s^{-1}$ is found in A1835 and 246
km $\rm s^{-1}$ in A2204 (Bambic et al. 2018). 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$\alpha$ emission in most cool
core clusters. Many studies of the inner cluster core have shown that
H$\alpha$ 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. 2003a; Fabian et al. 2006, 2016; Crawford et al. 2005b). No strong
evidence of significant cooling below $\sim$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$\alpha$ nebula and creates fast and energetic
particles (Fabian et al. 2011). 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 $\sim$10 $\rm M_{\odot}\rm\,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$\alpha$
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 $\sim$ 50K and/or NIR H2 lines at $\sim$ 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 $\rm
M_{\odot}\rm\,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$\alpha$ 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
$z=0.2153$ with $M_{500}=1.25\times 10^{15}\rm M_{\odot}$ (Piffaretti et al.
2011) and a high X-ray bolometric luminosity of 4.1$\times 10^{45}$ $\,\rm
ergs^{-1}$ and a classical cooling rate 111In this work, we use the definition
of the classical cooling rate as the ratio of gas mass enclosed in a radius
with a radiative cooling time of 7.7 Gyr to the cooling time (see e.g.
McDonald et al. 2018). of 1500-1900 $\rm M_{\odot}\rm\,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$\alpha$ luminosity of 3.2$\times 10^{43}$ $\,\rm
ergs^{-1}$ making it one of the most luminous optical nebulae. The observed UV
flux indicates a strong star formation rate of 130 $\rm M_{\odot}\rm\,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$\pm$0.1$\times 10^{10}$ $\rm
M_{\odot}$ (Vantyghem et al. 2018). 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$\pm 15$ km $\rm
s^{-1}$ for that gas and the central molecular filament shows a velocity range
of $\sim$260 $\pm 11$ km $\rm s^{-1}$ in RXCJ1504. A dynamically young gas
structure also shows an offset velocity of $\sim$ -250 km $\rm 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 $z=0.1283$ with $M_{500}=4.06\times
10^{14}\rm M_{\odot}$ (Piffaretti et al. 2011). It has a classical cooling
rate of 100$\pm$10 $\rm M_{\odot}\rm\,yr^{-1}$(McDonald et al. 2018), and
hosts a bright H$\alpha$ nebula of 1.5$\times 10^{42}$ $\,\rm ergs^{-1}$
(Wilman et al. 2006). The star formation rate is estimated to be 14 $\rm
M_{\odot}\rm\,yr^{-1}$ in IR or 4.3 $\rm M_{\odot}\rm\,yr^{-1}$ in FUV (O’Dea
et al. 2010). The BCG has a total molecular gas mass of
1.1$\pm$0.1$\times\,10^{10}$ $\rm M_{\odot}$ (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 $\sim$590 km $\rm s^{-1}$. On a
larger scale of $\sim$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: $H_{0}=73\rm\
km^{-1}Mpc^{-1}$, $\Omega_{\rm M}=0.27$, $\Omega_{\rm\Lambda}=0.73$.

##  2 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 $\chi^{2}$ in low count statistics. We
adopt 68.3 per cent confidence level (1$\sigma$ uncertainty at $\Delta C=1$)
for measurements. For upper/lower limits, we only quote the 90 per cent
confidence level (2$\sigma$ uncertainty at $\Delta C=2.71$) uncertainty,
unless otherwise stated.

Table 1: Observational details for RXCJ1504.1-0248 and A1664.

Name | Redshift |  $D_{L}$ (Mpc) | Scale (kpc/arcsec) | Obsid | Total RGS clean time (ks) | $N_{\rm H}(10^{20}\rm\,cm^{-2})$
---|---|---|---|---|---|---
RXCJ1504.1-0248 | 0.21530 | 1030 | 3.36 | 0840580101/201 | 220 | 8.34
A1664 | 0.12832 | 579 | 2.20 | 0840580301/401 | 222 | 12.8

The redshifts are taken from the NED database (https://ned.ipac.caltech.edu/).
The total Galactic column density $N_{\rm H}$ is taken from the UK Swift
Science Data Centre (see http://www.swift.ac.uk/analysis/nhtot/; Kalberla et
al. 2005; Willingale et al. 2013).

##  3 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 $T$ and emission measure $EM=n_{\rm e}n_{\rm H}V$, where
$n_{\rm e}$ and $n_{\rm H}$ are electron and proton densities and $V$ is the
volume of the emitting gas. We use the default value of $n_{\rm 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

$\frac{dEM(T)}{dT}=\frac{5\dot{M}k}{2\mu m_{\rm H}\Lambda(T)}$ (1)

where $k$ is the Boltzmann constant, $\mu$ is the mean particle weight,
$m_{\rm H}$ is the proton mass and $\Lambda(T)$ 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 $s$ and the wavelength shift $\Delta\lambda$ are the
free parameters. The scale factor fit for the amount of line broadening and
the wavelength shift corrects for the centroid of emission.

###  3.1 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), one-stage 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 LABEL:tab:A1664cf. 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 one-stage model. The other fit parameters are also
consistent between these two models. We hence conclude a cooling rate to 0.7
keV of 180$\pm$40 $\rm M_{\odot}\rm\,yr^{-1}$ and a residual cooling rate of
from 0.7 keV less than 53$\rm M_{\odot}\rm\,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 C-stat 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, one-
stage, cooling rate of 34$\pm$6$\rm M_{\odot}\rm\,yr^{-1}$ from the current
data (Model 1).

Figure 1: Stacked first order RGS spectrum of RXCJ1504.1-0248 in rest
wavelength. The RGS spectrum is shown in black and the best fit two-stage
cooling models is seen in blue. The background is seen in grey. Strong
emission lines are labelled in red. The spectrum is overbinned by a factor of
6 for plotting purposes.  Figure 2: As Fig. 1, the stacked RGS spectrum of
A1664 with the best fit complete cooling flow model in red in rest wavelength.
Table 2: XMM/RGS fit parameters for RXCJ1504.

| Model 1 | Model 2 | Model 3 | Model 4
---|---|---|---|---
Spectral components | 1cie | 1cie+1cf | 1cie+1cf | 1cie+2cf
C-stat/DoF | 841/682 | 833/681 | 818/680 | 818/680
Fe/H | 0.68$\pm$0.06 | 0.74$\pm$0.07 | 0.76$\pm$0.08 | 0.76$\pm$0.07
O/H | 0.41$\pm$0.06 | 0.39$\pm$0.07 | 0.44$\pm$0.07 | 0.44$\pm$0.07
$T_{\rm H}$ (keV) | 5.4$\pm$0.2 | 5.6$\pm$0.2 | 6.3$\pm$0.4 | 6.3$\pm$0.4, 0.7
$T_{\rm min}$ (keV) | n/a | 0.01 | 0.7$\pm$0.1 | 0.7, 0.01
$\dot{M}_{\rm H}$ ($\rm M_{\odot}\rm\,yr^{-1}$) | n/a | 50$\pm$20 | 190$\pm$60 | 180$\pm$40
$\dot{M}_{\rm C}$ ($\rm M_{\odot}\rm\,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. $T_{\rm H}$ and
$T_{\rm min}$ are the cie temperature and the minimum temperature of the
associated cf component, respectively. $\dot{M}_{\rm H}$ is the cooling rate
between $T_{\rm H}$ and $T_{\rm min}$. $\dot{M}_{\rm C}$ is the residual
cooling rate between 0.7 and 0.01 keV in the two-stage model.

Table 3: XMM/RGS fit parameters for A1664. Models, parameters and labels are the same as Table 2. | Model 1 | Model 2 | Model 3 | Model 4
---|---|---|---|---
Spectral components | 1cie | 1cie+1cf | 1cie+1cf | 1cie+2cf
C-stat/DoF | 929/680 | 909/677 | 909/676 | 909/676
Fe/H | 0.37$\pm$0.03 | 0.49$\pm$0.05 | 0.49$\pm$0.06 | 0.50$\pm$0.06
O/H | 0.23$\pm$0.05 | 0.28$\pm$0.06 | 0.58$\pm$0.06 | 0.29$\pm$0.06
$T_{\rm H}$ (keV) | 1.95$\pm$0.07 | 2.2$\pm$0.1 | 2.2$\pm$0.1 | 2.1$\pm$0.1, 0.7
$T_{\rm min}$ (keV) | n/a | 0.01 |  $<$0.7 | 0.7, 0.01
$\dot{M}_{\rm H}$ ($\rm M_{\odot}\rm\,yr^{-1}$) | n/a | 34$\pm$6 | 34$\pm$6 | 40$\pm$20
$\dot{M}_{\rm C}$ ($\rm M_{\odot}\rm\,yr^{-1}$) | n/a | n/a | n/a | 30$\pm$10

###  3.2 Turbulence

####  3.2.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$\pm 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.

Figure 3: Left panel: The surface brightness profile of RXCJ1504 (Data) and
the components of three-Gaussian fit. Right panel: The surface brightness
profile and Gaussian components of A1664.

####  3.2.2 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
$s=1$ for the first order and $s=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 single-temperature (1cie) model and fit the micro-
turbulent velocity ($v_{\rm turb}$) of the cie components. The 1-dimensional
turbulent velocity is then $v_{\rm 1D}$=$v_{\rm mic}$/$\sqrt{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 $\rm s^{-1}$.

Table 4: Turbulent velocity limits for RXCJ1504 and A1664.

|  | 1st order spectrum |  |  | 1st and 2nd order spectra |  |
---|---|---|---|---|---|---|---
|  | Total width | Central | Central $\&$ intermediate | Total width | Central | Central $\&$ intermediate
RXCJ1504 | v1D (km $\rm s^{-1}$) | 600$\pm$100 | 500$\pm$100 |  $<$260 | 550$\pm$90 | 300$\pm$100 |  $<$310
| C-stat/DoF | 840/682 | 839/682 | 839/682 | 3111/2093 | 3082/2093 | 3087/2093
A1664 | v1D (km $\rm s^{-1}$) | 800$\pm$200 |  $<$530 |  $<$420 | 700$\pm$100 |  $<$300 |  $<$320
| C-stat/DoF | 941/680 | 931/680 | 938/680 | 2934/2089 | 2919/2089 | 2941/2089

The columns of total width represent the line width measurements without the
correction of line broadening. The next column ’Central’ is the turbulent
velocity measured by correcting the spatial broadening using the central
Gaussian only. The velocity limits in the last column ’Central $\&$
intermediate’ are corrected by using both the central and intermediate
Gaussians in the spatial profile.

##  4 Discussion

###  4.1 Soft X-ray and cooler gas

It is possible to achieve a 3$\sigma$ measurement of the cooling rate by the
two-stage model (Model 4) or better for many nearby, $z<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$\sigma$ confidence level.
Both targets are already well-studied 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; A1664: Edge 2001; Wilman et al. 2006; Kirkpatrick et al. 2009; O’Dea
et al. 2010; Russell et al. 2014; Calzadilla et al. 2019). It is then 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$\alpha$ 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<z<0.6$) with known optical nebula and will report elsewhere.

The H$\alpha$ 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$\alpha$ 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$\alpha$ 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$\alpha$
nebula (e.g. Perseus: Fabian et al. 2003a; Fabian et al. 2006; Centaurus:
Crawford et al. 2005b, Fabian et al. 2016; A1795: Fabian et al. 2001, Crawford
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$\alpha$
luminosity

$15\times\,L_{\rm H\alpha}=3/2\times\dot{M}_{\rm neb}(\frac{kT}{\mu m_{p}}),$
(2)

which simplifies to

$\dot{M}_{\rm neb}=0.99\times(\frac{L_{\rm H\alpha}}{10^{40}\,\rm erg\,\rm
s^{-1}})(\frac{kT}{\rm 1keV})^{-1}\rm M_{\odot}\,\rm yr^{-1}.$ (3)

For our targets, $\dot{M}_{\rm neb}$ is 4.6$\times 10^{3}\rm
M_{\odot}\rm\,yr^{-1}$ for RXCJ1504 and 2.1$\times 10^{2}\rm
M_{\odot}\rm\,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$\alpha$ 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$\alpha$ 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
$z=0.596$. It has a similar molecular gas mass of 2$\times 10^{10}$ $\rm
M_{\odot}$ (Russell et al. 2017) embedded in an optical line-emitting nebula
with an H$\alpha$ luminosity of 8.52$\pm 0.5\times 10^{43}$ $\,\rm ergs^{-1}$
(McDonald et al. 2014). For the Phoenix cluster, Pinto et al. (2018) reported
a cooling rate of 350${}^{+150}_{-120}$ $M_{\odot}\rm\,yr^{-1}$ below 2 keV at
68 per cent confidence level. Both the H$\alpha$ luminosity and the cooling
rate are twice of those measured in RXCJ1504. We calculate $\dot{M}_{\rm neb}$
to be 4.3$\times 10^{3}\rm M_{\odot}\rm\,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
$M_{\odot}\rm\,yr^{-1}$ that is comparable to the observed cooling rate at the
1$\sigma$ 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$\alpha$ 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$\times
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$\times
10^{45}\rm\,ergs^{-1}$ and 2.1$\times 10^{43}\rm\,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$\times
10^{43}\rm\,ergs^{-1}$ (Pulido et al. 2018). 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$\sigma$ 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 $\rm M_{\odot}\rm\,yr^{-1}$ and 0.012 $\rm
M_{\odot}\rm\,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}\rm M_{\odot}$. 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 $\rm
M_{\odot}\rm\,yr^{-1}$ in RXCJ1504 and 20-30 $\rm M_{\odot}\rm\,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.

###  4.2 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 $\rm s^{-1}$ and $1310^{+570}_{-670}$ km $\rm s^{-1}$
at 68 per cent confidence level, respectively (Sanders et al. 2011; Sanders &
Fabian 2013). Our results using the new data show a much tighter limit (of 300
km $\rm s^{-1}$ ). The turbulent velocity of our targets is comparable to the
velocity measured in many bright clusters, e.g. $<$211 km $\rm s^{-1}$ in
A1835 (Sanders & Fabian 2013; Bambic et al. 2018), $<$400 km $\rm s^{-1}$ in
2A0335+096 (Pinto et al. 2015), $\sim$164 km $\rm s^{-1}$ in Perseus (Hitomi
Collaboration et al. 2016) and $<$370 km $\rm s^{-1}$ in the Phoenix cluster
(Pinto et al. 2018).

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$\pm 17$ km $\rm 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 $c_{\rm s}=\sqrt{\gamma kT/\mu\rm
m_{\rm p}}$, where $\gamma$ is the adiabatic index, which is 5/3 for ideal
monatomic gas, $\mu=0.6$ is the mean particle mass and $\rm m_{\rm p}$ is the
proton mass. This gives a sound speed of 1300 km $\rm s^{-1}$ and 750 km $\rm
s^{-1}$ for RXCJ1504 and A1664, respectively. In this work, we calculate the
1-D Mach number for turbulence $M=V_{\rm 1D}/c_{\rm s}$. The turbulent
velocity then has a Mach number $M$ 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

$\frac{\epsilon_{\rm turb}}{\epsilon_{\rm therm}}=\gamma M^{2}.$ (4)

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 $\rm 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 ($\sim$ 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. Fujita et al. (2020) presents a
time-dependent 1D simulation of heating in cool core clusters with outbursts
reaching $10^{46}$ erg$\rm s^{-1}$ from the AGN on Gyr timescale. The central
density and temperature profiles make large excursions on this timescale with
energy advected during the outbursts. However, Panagoulia et al. (2014a),
Hogan et al. (2017) and Babyk et al. (2018) show a universal inner entropy
shape which would not be seen if the central gas properties are cycling up and
down. The issue of how the energy is replaced or flows remains open, with
sound waves remaining a possibility.

###  4.3 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$\pm 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$\pm 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 $\rm 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 $\rm 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$\pm 400$ km $\rm s^{-1}$ in the Phoenix cluster. A similar velocity of
$\sim$1000 km $\rm s^{-1}$ is found in the non cool core Coma cluster, while
gas in the Perseus cluster is more relaxed at 480$\pm 210$ km $\rm 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 $\rm s^{-1}$ (Russell et al.
2014). The blueshifted component is seen at a velocity of 571$\pm 7$ km $\rm
s^{-1}$ from CO(3-2) in our line-of-sight, with a FWHM of 190$\pm 20$ km $\rm
s^{-1}$. Given the large uncertainty in the X-ray measurements, the velocities
of the molecular and X-ray gas are consistent within 1$\sigma$. 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 $\rm s^{-1}$. The molecular gas will be shocked unless it
comoves with the soft X-ray emitting gas.

Figure 4: The energy dependent surface brightness profile of A1664.

###  4.4 Embedded multilayer cooling flow in RXCJ1504

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.

There are larger amounts of cold obscuring material in the core of galaxy
clusters, e.g. the Centaurus cluster (A3526) (Crawford et al. 2005b; Sanders &
Fabian 2006). 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\times 10^{21}\rm
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 & 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 $F_{\lambda}$. The
fraction of the emitted energy transmitted through one sheet of absorbing gas
is $f_{\lambda}$. We can write this transmission fraction as

$f_{\lambda}=e^{-\sigma(E)\Delta n_{\rm H}},$ (5)

where $-\sigma(E)$ is the absorption cross-section and $\Delta n_{\rm H}$ is
the column density of one sheet of absorbing gas. The total observed flux is
then

$F_{\rm
tot}=f_{\lambda}F_{\lambda}+f_{\lambda}^{2}F_{\lambda}+...+f_{\lambda}^{n_{\rm
sheet}}F_{\lambda}=F_{\lambda}\sum_{m=1}^{n_{\rm sheet}}f_{\lambda}^{m},$ (6)

where $n_{\rm 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

$F_{\rm tot}=F_{\lambda}\frac{1-f_{\lambda}^{n_{\rm
sheet}}}{1-f_{\lambda}^{{}_{\rm\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,}}}.$ (7)

This suggests only $n_{\rm sheet}$ and the total column density $n_{\rm
H,tot}$ are the additional free parameters. Note that $n_{\rm H,tot}=n_{\rm
sheet}\Delta n_{\rm H}$. In the large $n_{\rm sheet}$ limit, Equation 5 can be
expanded and Equation 7 rewritten as

$F_{\rm tot}=F_{\lambda}n_{\rm sheet}\frac{1-e^{-\sigma(E)\Delta n_{\rm
H}}}{\sigma(E)\Delta n_{\rm H}}.$ (8)

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
$n_{\rm H,tot}/n_{\rm 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 $n_{\rm 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 $n_{\rm sheet}$ and $n_{\rm H,tot}$ parameter space. We
apply our model to RXCJ1504 and fit for minimal C-stat for each pair of
$n_{\rm sheet}$ and $n_{\rm 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\times 10^{21}\rm cm^{-2}$. However, the difference
between the C-stats is less than 0.3 on this valley. This means $n_{\rm
sheet}$ and $n_{\rm 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 $n_{\rm sheet}\sim 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$\pm 90$ $\rm M_{\odot}\rm\,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\times 10^{22}\rm cm^{-2}$, the cooling rate is 520$\pm
30$ $\rm M_{\odot}\rm\,yr^{-1}$. The intrinsic absorption corrected cooling
rate is consistent between 1 sheet and 10 sheets model at the 1$\sigma$ level.
For an order of magnitude increase in the cooling rate above 0.7 keV in
RXCJ1504, it can contribute $\sim$40 per cent of the energy required to power
the UV/optical line-emitting nebula, but not all. In other massive clusters
with $\dot{M}_{\rm neb}>100\rm M_{\odot}\rm\,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\times 10^{21}\rm cm^{-2}$, we
only measure a cooling rate of 31$\pm 7$ $\rm M_{\odot}\rm\,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 $n_{\rm
sheet}$ and $n_{\rm 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 $\rm M_{\odot}\rm\,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).

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.

##  5 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 $z=0.2153$
and A1664 at $z=01283$. The cooling rate is measured to be 180$\pm$40$\rm
M_{\odot}\rm\,yr^{-1}$ and 34$\pm$6$\rm M_{\odot}\rm\,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 53$\rm M_{\odot}\rm\,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$\times 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 $\rm 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, 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$\pm 0.08$ keV is
blueshifted from the systemic velocity of the cluster at 750${}^{+800}_{-280}$
km $\rm 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$\pm 30$ $\rm
M_{\odot}\rm\,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.

## Acknowledgements

This work is based on observations obtained with XMM-Newton, an ESA science
mission funded by ESA Member States and USA (NASA). HL thanks the co-authors
for their contributions and the referee for the support.

## Data Availability

No new data were generated or analysed in support of this research.

## References

  * Allen & Fabian (1997) Allen S. W., Fabian A. C., 1997, MNRAS, 286, 583
  * Allen et al. (1996) Allen S. W., Fabian A. C., Edge A. C., Bautz M. W., Furuzawa A., Tawara Y., 1996, MNRAS, 283, 263
  * Ascasibar & Markevitch (2006) Ascasibar Y., Markevitch M., 2006, ApJ, 650, 102
  * Babyk et al. (2018) Babyk I. V., McNamara B. R., Nulsen P. E. J., Russell H. R., Vantyghem A. N., Hogan M. T., Pulido F. A., 2018, ApJ, 862, 39
  * Bambic et al. (2018) Bambic C. J., Pinto C., Fabian A. C., Sanders J., Reynolds C. S., 2018, MNRAS, 478, L44
  * Barret et al. (2018) Barret D., et al., 2018, in den Herder J.-W. A., Nikzad S., Nakazawa K., eds, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series Vol. 10699, Space Telescopes and Instrumentation 2018: Ultraviolet to Gamma Ray. p. 106991G (arXiv:1807.06092), doi:10.1117/12.2312409
  * Bîrzan et al. (2004) Bîrzan L., Rafferty D. A., McNamara B. R., Wise M. W., Nulsen P. E. J., 2004, ApJ, 607, 800
  * Blanton et al. (2004) Blanton E. L., Sarazin C. L., McNamara B. R., Clarke T. E., 2004, ApJ, 612, 817
  * Böhringer et al. (2005) Böhringer H., Burwitz V., Zhang Y. Y., Schuecker P., Nowak N., 2005, ApJ, 633, 148
  * Bondi (1952) Bondi H., 1952, MNRAS, 112, 195
  * Calzadilla et al. (2019) Calzadilla M. S., et al., 2019, ApJ, 875, 65
  * Cavagnolo et al. (2009) Cavagnolo K. W., Donahue M., Voit G. M., Sun M., 2009, ApJS, 182, 12
  * Churazov et al. (2002) Churazov E., Sunyaev R., Forman W., Böhringer H., 2002, MNRAS, 332, 729
  * Churazov et al. (2013) Churazov E., Ruszkowski M., Schekochihin A., 2013, MNRAS, 436, 526
  * Crawford et al. (2005a) Crawford C. S., Sanders J. S., Fabian A. C., 2005a, MNRAS, 361, 17
  * Crawford et al. (2005b) Crawford C. S., Hatch N. A., Fabian A. C., Sanders J. S., 2005b, MNRAS, 363, 216
  * Dunn & Fabian (2006) Dunn R. J. H., Fabian A. C., 2006, MNRAS, 373, 959
  * Edge (2001) Edge A. C., 2001, MNRAS, 328, 762
  * Edge et al. (2002) Edge A. C., Wilman R. J., Johnstone R. M., Crawford C. S., Fabian A. C., Allen S. W., 2002, MNRAS, 337, 49
  * Fabian (1994) Fabian A. C., 1994, ARA&A, 32, 277
  * Fabian (2012) Fabian A. C., 2012, ARA&A, 50, 455
  * Fabian et al. (2001) Fabian A. C., Sanders J. S., Ettori S., Taylor G. B., Allen S. W., Crawford C. S., Iwasawa K., Johnstone R. M., 2001, MNRAS, 321, L33
  * Fabian et al. (2003a) Fabian A. C., Sanders J. S., Crawford C. S., Conselice C. J., Gallagher J. S., Wyse R. F. G., 2003a, MNRAS, 344, L48
  * Fabian et al. (2003b) Fabian A. C., Sanders J. S., Allen S. W., Crawford C. S., Iwasawa K., Johnstone R. M., Schmidt R. W., Taylor G. B., 2003b, MNRAS, 344, L43
  * Fabian et al. (2006) Fabian A. C., Sanders J. S., Taylor G. B., Allen S. W., Crawford C. S., Johnstone R. M., Iwasawa K., 2006, MNRAS, 366, 417
  * Fabian et al. (2011) Fabian A. C., Sanders J. S., Williams R. J. R., Lazarian A., Ferland G. J., Johnstone R. M., 2011, MNRAS, 417, 172
  * Fabian et al. (2016) Fabian A. C., et al., 2016, MNRAS, 461, 922
  * Fabian et al. (2017) Fabian A. C., Walker S. A., Russell H. R., Pinto C., Sanders J. S., Reynolds C. S., 2017, MNRAS, 464, L1
  * Ferland et al. (2009) Ferland G. J., Fabian A. C., Hatch N. A., Johnstone R. M., Porter R. L., van Hoof P. A. M., Williams R. J. R., 2009, MNRAS, 392, 1475
  * Fujita et al. (2020) Fujita Y., Cen R., Zhuravleva I., 2020, MNRAS, 494, 5507
  * Gaspari et al. (2017) Gaspari M., Temi P., Brighenti F., 2017, MNRAS, 466, 677
  * Giacintucci et al. (2011) Giacintucci S., Markevitch M., Brunetti G., Cassano R., Venturi T., 2011, A&A, 525, L10
  * Hitomi Collaboration et al. (2016) Hitomi Collaboration et al., 2016, Nature, 535, 117
  * Hlavacek-Larrondo et al. (2012) Hlavacek-Larrondo J., Fabian A. C., Edge A. C., Ebeling H., Sanders J. S., Hogan M. T., Taylor G. B., 2012, MNRAS, 421, 1360
  * Hogan et al. (2017) Hogan M. T., et al., 2017, ApJ, 851, 66
  * Juett et al. (2006) Juett A. M., Schulz N. S., Chakrabarty D., Gorczyca T. W., 2006, ApJ, 648, 1066
  * Kaastra et al. (2001) Kaastra J. S., Ferrigno C., Tamura T., Paerels F. B. S., Peterson J. R., Mittaz J. P. D., 2001, A&A, 365, L99
  * Kalberla et al. (2005) Kalberla P. M. W., Burton W. B., Hartmann D., Arnal E. M., Bajaja E., Morras R., Pöppel W. G. L., 2005, A&A, 440, 775
  * Kirkpatrick et al. (2009) Kirkpatrick C. C., et al., 2009, ApJ, 697, 867
  * Liu et al. (2019) Liu H., Pinto C., Fabian A. C., Russell H. R., Sanders J. S., 2019, MNRAS, 485, 1757
  * Liu et al. (2020) Liu H., Fabian A. C., Pinto C., 2020, MNRAS, 497, 1256
  * Lodders & Palme (2009) Lodders K., Palme H., 2009, Meteoritics and Planetary Science Supplement, 72, 5154
  * McDonald et al. (2012) McDonald M., et al., 2012, Nature, 488, 349
  * McDonald et al. (2013) McDonald M., Benson B., Veilleux S., Bautz M. W., Reichardt C. L., 2013, ApJ, 765, L37
  * McDonald et al. (2014) McDonald M., et al., 2014, ApJ, 784, 18
  * McDonald et al. (2018) McDonald M., Gaspari M., McNamara B. R., Tremblay G. R., 2018, ApJ, 858, 45
  * McNamara & Nulsen (2012) McNamara B. R., Nulsen P. E. J., 2012, New Journal of Physics, 14, 055023
  * McNamara et al. (2006) McNamara B. R., et al., 2006, ApJ, 648, 164
  * O’Dea et al. (2008) O’Dea C. P., et al., 2008, ApJ, 681, 1035
  * O’Dea et al. (2010) O’Dea K. P., et al., 2010, ApJ, 719, 1619
  * Ogorzalek et al. (2017) Ogorzalek A., et al., 2017, MNRAS, 472, 1659
  * Ogrean et al. (2010) Ogrean G. A., Hatch N. A., Simionescu A., Böhringer H., Brüggen M., Fabian A. C., Werner N., 2010, MNRAS, 406, 354
  * Olivares et al. (2019) Olivares V., et al., 2019, A&A, 631, A22
  * Panagoulia et al. (2014a) Panagoulia E. K., Fabian A. C., Sanders J. S., 2014a, MNRAS, 438, 2341
  * Panagoulia et al. (2014b) Panagoulia E. K., Fabian A. C., Sanders J. S., Hlavacek-Larrondo J., 2014b, MNRAS, 444, 1236
  * Peterson et al. (2003) Peterson J. R., Kahn S. M., Paerels F. B. S., Kaastra J. S., Tamura T., Bleeker J. A. M., Ferrigno C., Jernigan J. G., 2003, ApJ, 590, 207
  * Piffaretti et al. (2011) Piffaretti R., Arnaud M., Pratt G. W., Pointecouteau E., Melin J. B., 2011, A&A, 534, A109
  * Pinto et al. (2013) Pinto C., Kaastra J. S., Costantini E., de Vries C., 2013, A&A, 551, A25
  * Pinto et al. (2015) Pinto C., et al., 2015, A&A, 575, A38
  * Pinto et al. (2016) Pinto C., et al., 2016, MNRAS, 461, 2077
  * Pinto et al. (2018) Pinto C., Bambic C. J., Sanders J. S., Fabian A. C., McDonald M., Russell H. R., Liu H., Reynolds C. S., 2018, MNRAS, 480, 4113
  * Polles et al. (2021) Polles F. L., et al., 2021, arXiv e-prints, p. arXiv:2103.09842
  * Pulido et al. (2018) Pulido F. A., et al., 2018, ApJ, 853, 177
  * Rafferty et al. (2006) Rafferty D. A., McNamara B. R., Nulsen P. E. J., Wise M. W., 2006, ApJ, 652, 216
  * Rafferty et al. (2008) Rafferty D. A., McNamara B. R., Nulsen P. E. J., 2008, ApJ, 687, 899
  * Russell et al. (2014) Russell H. R., et al., 2014, ApJ, 784, 78
  * Russell et al. (2017) Russell H. R., et al., 2017, ApJ, 836, 130
  * Russell et al. (2019) Russell H. R., et al., 2019, MNRAS, 490, 3025
  * Ruszkowski et al. (2018) Ruszkowski M., Yang H. Y. K., Reynolds C. S., 2018, ApJ, 858, 64
  * Salomé & Combes (2003) Salomé P., Combes F., 2003, A&A, 412, 657
  * Sanders & Fabian (2006) Sanders J. S., Fabian A. C., 2006, MNRAS, 370, 63
  * Sanders & Fabian (2013) Sanders J. S., Fabian A. C., 2013, MNRAS, 429, 2727
  * Sanders et al. (2008) Sanders J. S., Fabian A. C., Allen S. W., Morris R. G., Graham J., Johnstone R. M., 2008, MNRAS, 385, 1186
  * Sanders et al. (2010) Sanders J. S., Fabian A. C., Smith R. K., Peterson J. R., 2010, MNRAS, 402, L11
  * Sanders et al. (2011) Sanders J. S., Fabian A. C., Smith R. K., 2011, MNRAS, 410, 1797
  * Sanders et al. (2020) Sanders J. S., et al., 2020, A&A, 633, A42
  * Tozzi et al. (2015) Tozzi P., et al., 2015, A&A, 580, A6
  * Urdampilleta et al. (2017) Urdampilleta I., Kaastra J. S., Mehdipour M., 2017, A&A, 601, A85
  * Vantyghem et al. (2018) Vantyghem A. N., et al., 2018, ApJ, 863, 193
  * Werner et al. (2009) Werner N., Zhuravleva I., Churazov E., Simionescu A., Allen S. W., Forman W., Jones C., Kaastra J. S., 2009, MNRAS, 398, 23
  * Willingale et al. (2013) Willingale R., Starling R. L. C., Beardmore A. P., Tanvir N. R., O’Brien P. T., 2013, MNRAS, 431, 394
  * Wilman et al. (2006) Wilman R. J., Edge A. C., Swinbank A. M., 2006, MNRAS, 371, 93
  * Wilman et al. (2009) Wilman R. J., Edge A. C., Swinbank A. M., 2009, MNRAS, 395, 1355
  * XRISM Science Team (2020) XRISM Science Team 2020, arXiv e-prints, p. arXiv:2003.04962
  * Zhuravleva et al. (2018) Zhuravleva I., Allen S. W., Mantz A., Werner N., 2018, ApJ, 865, 53
  * de Plaa et al. (2012) de Plaa J., Zhuravleva I., Werner N., Kaastra J. S., Churazov E., Smith R. K., Raassen A. J. J., Grange Y. G., 2012, A&A, 539, A34