11institutetext: Department of Clinical Neurosciences, University of Cambridge
22institutetext: Department of Applied Maths and Theoretical Physics,
University of Cambridge
{22email<EMAIL_ADDRESS>
# D-CoRP: Differentiable Connectivity Refinement for Functional Brain Networks

Haoyu Hu 11 Hongrun Zhang 11 Chao Li corresponding author1122

###### Abstract

Brain network is an important tool for understanding the brain, offering
insights for scientific research and clinical diagnosis. Existing models for
brain networks typically primarily focus on brain regions or overlook the
complexity of brain connectivities. MRI-derived brain network data is commonly
susceptible to connectivity noise, underscoring the necessity of incorporating
connectivities into the modeling of brain networks. To address this gap, we
introduce a differentiable module for refining brain connectivity. We develop
the multivariate optimization based on information bottleneck theory to
address the complexity of the brain network and filter noisy or redundant
connections. Also, our method functions as a flexible plugin that is adaptable
to most graph neural networks. Our extensive experimental results show that
the proposed method can significantly improve the performance of various
baseline models and outperform other state-of-the-art methods, indicating the
effectiveness and generalizability of the proposed method in refining brain
network connectivity. The code will be released for public availability.

###### Keywords:

Functional Brain Network Graph Structure Learning Information Bottleneck.

##  1 Introduction

The human brain comprises multiple sophisticated brain regions, empowering it
to navigate intricate cognitive processes and handle everyday tasks [12].
Recent advances in neuroscience suggest that brain function is fulfilled by
the coordinated activities of brain regions through connectivities [16]. Brain
networks, graphical representations of the brain, hold the potential to
characterize regional and global brain activities [5, 13, 26]. Functional
brain networks, derived from function MRI (fMRI), is one of the most important
tools for characterizing the brain function. However, it remains challenging
to effectively leverage the information in functional brain networks due to
their non-Euclidean structure.

Graph neural networks (GNNs) [32] promise to offer effective toolkits to
analyze graph data [24]. For example, BrainGNN is proposed to effectively
aggregate the information from functional brain networks for learning tasks
[14, 25]. Graph pooling methods are developed for characterizing network
modularity, e.g., TopK, a score-based pooling method [8] and DiffPool, a
hierarchical pooling method [29]. However, these models only focus on node
features without considering edge connections, which may particularly limit
the model performance on MRI-derived brain networks that are prone to edge
noises or redundancies[27].

Previous studies proposed graph structure learning (GSL) to denoise and
optimize graph structure [33]. Despite successes, existing GSL models are
developed for specific type of graphs (e.g. community networks, population
graphs) using rigid optimization approaches (e.g., community property
optimization), and they are often limited in generalizability to other graphs
once trained. Scanty models have been proposed for brain networks so far.
Though a recent model [31] proposes to reconstruct brain networks through time
windows extracting features from fMRI, it cannot work directly on graphs
(e.g., correlation matrices).

Information bottleneck (IB) theory [19] employs a probabilistic optimization
between predictive performance and information compression by focusing on
model-processed features of the input graph. Unlike typical GSL methods, IB is
not restricted to a specific type of graph, which brings the advantages of
generalizability while effectively removing noisy and irrelevant information.
These advantages could equip IB as a competitive framework for denoising
graphs.

To denoise the input graph, previous IB-based methods perform subgraph
filtering that updates the graph structure, e.g., subgraph information
bottleneck (SIB) [30]. However, this approach treats the entire subgraph as
processed features without performing differentiable refinement, which may
lead to less desirable performance in removing edges. By contrast, other
methods perform edge and node filtering, preserving the overall structure,
e.g., VIB-GSL [18]. Nevertheless, these methods assume the optimized graph
follows a standard Gaussian distribution, which may over-simplify the true
distributions of edge connections that vary across different brain regions.
Finally, all these methods are implemented as a fixed framework, lacking
adaptability to existing GNN frameworks and hindering wider applicability.
Thus, it remains a challenge to develop effective models to refine edge
connections integrated into existing frameworks.

We propose an edge-denoising method for functional brain networks based on IB
theory, namely Differentiable Connection Refinement Plugin (D-CoRP), serving
to detect and remove noisy or irrelevant connections from functional brain
networks as an adaptable plugin module. Our main contributions include:

1) We design a new differentiable sampling method that renders the plugin
learnable, considering the probability of completely removing noisy
connections.

2) We propose a new optimization strategy to model each connection
independently with multivariate distribution targets in an IB-adopted
framework.

3) As far as we know, this is the first plugin for connectivity refinement of
brain networks, adding minimal computation to backbone training. We also
provide an efficient mode for flexible adoption.

We evaluated our plugin on three datasets of functional brain networks. D-CoRP
is shown to effectively remove noisy and irrelevant connections and enhance
the performance of the GNN backbones, achieving a better and more stable
performance compared with other state-of-the-art methods.

##  2 Methods

###  2.1 Overview

Define the target graph to be refined as $\bm{G}_{t}$, which can be the
processed graph output from any layer of a graph model. Suppose
$\bf{G}_{t}=(\bm{X},\bm{A})$, where $\bm{X}\in\mathbb{R}^{N\times F}$ is the
matrix stacking node features and $\bm{A}\in\mathbb{R}^{N\times N}$ is the
adjacency matrix, with $N$ being the node number, $F$ being the dimension of
node features. We aim to find a refined functional brain network adjacency
matrix $\bm{A}_{re}$ based on the original $\bm{A}$, formulated as,

$\bm{A}_{re}=\bm{m}\odot\bm{A},$ (1)

where $\bm{m}\in\mathbb{R}^{N\times N}$ is the learnable edge mask with binary
elements, indicating the essential connections within brain networks. To learn
the refinement mask $\bm{m}$, we propose to apply the IB theory [19, 20],
which makes $\bm{m}$ as a processed feature connected to the model’s input and
output.

Figure 1: Model architecture. D-CoRP refines the graph connection in two
steps: 1) Differentiable sampling for the masking matrix (contains
Distribution Estimator and Differentiable Sampling); 2) Information bottleneck
optimization for the refined results as shown in IB Optimizer part. D-CoRP can
be directly plugged into the graph model as a plugin.

###  2.2 Differentiable Connectivity Refinement

The binary refinement matrix $\bm{m}$ is naturally modeled as being Bernoulli
distributed. However, the binary values are not differentiable, which hinders
the gradient propagation in the learning process; thus, it is intractable. To
address this, we relax the binary matrix $\bm{m}$ to a pseudo binary matrix,
with each value in it being between $[0,1]$, which can be uniformly sampled
from a predefined function. We design this function such that a sampled value
will be mostly closed to 0 or 1, resembling to the binary case (see Eq.(4) for
the specific expression), which is parameterized by a learning probability
matrix $\bm{\mathcal{P}}$, indicating the probabilities of corresponding edges
being preserved.

Generation of Selection Probabilities (Distribution Estimator) The probability
matrix $\bm{\mathcal{P}}$ is obtained via an auxiliary graph network (GAT)
with the node features $\bm{X}$ and adjacency matrix $\bm{A}$, followed by the
operation Sigmoid,

$\bm{Z}=\textrm{GAT}(\bm{X},\bm{A}).$ (2)
$\bm{\mathcal{P}}=\textrm{Sigmoid}\left(\bm{Z}\bm{Z}^{T}\right).$ (3)

where $\bm{Z},\bm{\mathcal{P}}\in\mathbb{R}^{N\times N}$. We choose GAT to
distil critical graph information due to its efficient attention mechanism.

Differentiable Sampling for Masking Matrix  The sampling is performed via a
differentiable function adapted from the concrete relaxation of Bernoulli
distribution [18]:

$\begin{split}\mathcal{F}(\bm{\mathcal{P}}|\bm{A})=\mathcal{F}(\bm{\mathcal{P}},\bm{\Pi}:\mathcal{T})=\frac{1}{2}(1-\textrm{tanh}(\frac{\bm{\Pi}+\bm{\mathcal{P}}-1}{\mathcal{T}}))=\frac{1}{1+e^{-\frac{2(\bm{\Pi}+\bm{\mathcal{P}}-1)}{\mathcal{T}}}}\end{split}$
(4)

where $\bm{\Pi}\sim\mathcal{U}(0,1)$ provides uniformly random sampling.
$\mathcal{T}$ represents the temperature for the function. This function can
directly relate the probability matrix to the sampling process as Bernoulli
distribution. When $\mathcal{T}\rightarrow 0$, we can treat
$\mathcal{F}(\bm{\mathcal{P}},\bm{\Pi}:\mathcal{T})\approx\textrm{Ber}(\bm{\mathcal{P}})$,
a Bernoulli distribution, while
$\mathcal{F}(\bm{\mathcal{P}},\bm{\Pi}:\mathcal{T})$ always keeps continuous
in the $[0,1]$ range. Unlike the previous concrete relaxation, Eq. (4)
achieves the closed interval [0,1] for $\bm{\mathcal{P}}$, ensuring the full
or zero probability of preserving brain connections. Also, it provides a
smoother change of gradient without computing the logarithm, preventing the
training process from exponential explosion.

###  2.3 Optimization with Information Bottleneck

In D-CoRP, the adjacency matrix $\bm{A}$ is the input data and $\bm{A}_{re}$
is the processed feature. Combined with general variational bounds developed
by Alemi et al. [1], we design the optimization for brain network refinement
(the specific derivation is provided in the supplementary):

$\begin{split}\mathcal{L}_{IB}=\mathbb{E}_{\bm{A}_{re}\sim
p_{\theta}(\bm{A}_{re}|\bm{A})}[-\textrm{log}q(\bm{Y}|\bm{A}_{re})]+\beta\textrm{KL}[p(\bm{A}_{re}|\bm{A}),r(\bm{A}_{re})]\end{split}$
(5)

Further introducing the upper bound of Eq. (5)’s right side, we conduct two
tractable optimization targets for brain network connectivity refinement:

$\begin{split}\mathcal{L}_{AIB}=\mathbb{E}_{\bm{A}_{re}\sim
p_{\theta}(\bm{A}_{re}|\bm{A})}[-\textrm{log}q(\bm{Y}|\bm{A}_{re})]+\beta\textrm{KL}[\mathcal{N}(\bm{\mathcal{P}};\bm{\Sigma}),\mathcal{N}(\bm{\mu}_{\bm{m}};\bm{\Sigma})]\end{split}$
(6) $\begin{split}\mathcal{L}_{AIB-E}=\mathbb{E}_{\bm{A}_{re}\sim
p_{\theta}(\bm{A}_{re}|\bm{A})}[-\textrm{log}q(\bm{Y}|\bm{A}_{re})]+\beta\textrm{KL}[p(\bm{m}|\bm{A}),r(\bm{m})]\end{split}$
(7)

where $r(\bm{m})$ is the prior distribution of $\bm{m}$ approximating
$p(\bm{m})$; $\bm{\mu}_{\bm{m}}$ and $\bm{\Sigma}$ are the parameters
approximating the multivariate Gaussian distribution of $\bm{m}$.
$\mathcal{L}_{AIB}$ is the optimization target of D-CoRP and
$\mathcal{L}_{AIB-E}$ is the optimization target of D-CoRP: Efficient, a
lighter version of D-CoRP. The proof can be found in the Supplementary. Since
the left-sides of Eq. (6) and (7) are already done by the model comparing the
prediction with ground truth, we only optimize an additional
$\textrm{KL}[\mathcal{N}(\bm{\mathcal{P}};\bm{\Sigma}),\mathcal{N}(\bm{\mu}_{\bm{m}};\bm{\Sigma})]$
or $\textrm{KL}[p(\bm{m}|\bm{A}),r(\bm{m})]$. For D-CoRP: Efficient, the
plugin adds an additional $\mathcal{O}(NF^{2})$ in computational complexity to
the plugged model.

|  |  | FCP | NKI-Rockland | ADHD-Normal
---|---|---|---|---|---
| Baseline | Plugin | MAE$\downarrow$ | RMSE$\downarrow$ | MAE$\downarrow$ | RMSE$\downarrow$ | MAE$\downarrow$ | RMSE$\downarrow$
Basic GNN Module | GCN [11] | None | 15.27$\pm$1.47 | 21.20$\pm$6.76 | 85.98$\pm$45.23 | 209.49$\pm$127.55 | 10.57$\pm$3.56 | 18.23$\pm$5.51
D-CoRP (Single) | 8.94$\pm$0.15 | 12.13$\pm$0.28 | 16.14$\pm$0.36 | 20.32$\pm$1.06 | 2.77$\pm$0.23 | 3.45$\pm$0.22
D-CoRP (Full) | 9.36$\pm$0.59 | 12.53$\pm$0.61 | 16.76$\pm$1.22 | 21.28$\pm$1.22 | 2.96$\pm$0.17 | 3.87$\pm$0.24
GAT [22] | None | 20.44$\pm$7.03 | 22.08$\pm$6.87 | 23.36$\pm$12.65 | 28.78$\pm$12.44 | 3.17$\pm$0.18 | 3.59$\pm$0.24
D-CoRP (Single) | 9.31$\pm$0.16 | 13.03$\pm$0.26 | 16.06$\pm$1.50 | 20.18$\pm$1.51 | 2.83$\pm$0.19 | 3.47$\pm$0.23
D-CoRP (Full) | 9.80$\pm$0.74 | 13.37$\pm$1.08 | 18.00$\pm$1.45 | 22.26$\pm$2.70 | 3.28$\pm$0.27 | 4.20$\pm$0.30
GIN [28] | None | 10.85$\pm$0.38 | 14.76$\pm$1.71 | 17.43$\pm$4.12 | 20.28$\pm$3.68 | 5.92$\pm$1.47 | 8.90$\pm$3.38
D-CoRP (Single) | 9.11$\pm$0.18 | 12.24$\pm$0.38 | 15.62$\pm$1.43 | 19.60$\pm$1.66 | 2.73$\pm$0.13 | 3.42$\pm$0.17
D-CoRP (Full) | 9.29$\pm$0.29 | 12.50$\pm$0.24 | 15.31$\pm$0.73 | 19.06$\pm$0.97 | 3.18$\pm$0.41 | 3.89$\pm$0.46
GraphSAGE [9] | None | 15.72$\pm$2.51 | 30.24$\pm$11.99 | 44.55$\pm$14.64 | 104.86$\pm$60.86 | 4.42$\pm$0.26 | 5.44$\pm$0.32
D-CoRP (Single) | 9.28$\pm$0.39 | 12.33$\pm$0.47 | 15.19$\pm$0.91 | 19.71$\pm$1.42 | 2.65$\pm$0.17 | 3.21$\pm$0.18
D-CoRP (Full) | 9.35$\pm$0.50 | 12.70$\pm$0.50 | 17.06$\pm$1.77 | 20.71$\pm$2.01 | 2.87$\pm$0.25 | 3.57$\pm$0.30
Basic Graph Pooling | DiffPool [29] | None | 12.10$\pm$0.79 | 18.01$\pm$4.55 | 26.96$\pm$9.63 | 39.04$\pm$19.44 | 4.06$\pm$0.21 | 6.56$\pm$0.86
D-CoRP (Single) | 8.55$\pm$0.17 | 11.83$\pm$0.60 | 16.32$\pm$2.44 | 20.39$\pm$2.32 | 2.99$\pm$0.20 | 3.75$\pm$0.18
D-CoRP (Full) | 8.84$\pm$0.70 | 11.96$\pm$0.73 | 24.62$\pm$10.24 | 28.56$\pm$10.85 | 3.08$\pm$0.26 | 3.89$\pm$0.36
TopK [8] | None | 11.05$\pm$0.69 | 13.23$\pm$1.11 | 27.49$\pm$9.42 | 37.61$\pm$10.82 | 6.14$\pm$3.62 | 6.70$\pm$3.42
D-CoRP (Single) | 12.22$\pm$6.40 | 15.21$\pm$5.72 | 19.88$\pm$8.72 | 24.10$\pm$9.11 | 2.73$\pm$0.20 | 3.78$\pm$0.10
D-CoRP (Full) | 13.98$\pm$4.46 | 16.34$\pm$4.48 | 30.13$\pm$8.31 | 35.16$\pm$9.36 | 2.84$\pm$0.20 | 3.58$\pm$0.28

Table 1: Performance of D-CoRP on Baselines. Underlines indicate local optimal
comparing None (D-CoRP), D-CoRP (Single) (added to a single GNN block), and
D-CoRP (Full) (added to all GNN blocks). Underlines indicate optimal results
of corresponding baselines and bold indicates global optimal results.

##  3 Experiments and Results

###  3.1 Task and Datasets

We focus on brain age prediction as the primary task. Brain age is shown as
associated with brain network features [17, 2], serving as an ideal benchmark
for model evaluation. The resting-state fMRI data from three datasets were
used: 1) 1000 Functional Connectomes Project (FCP) [3] (1,003 healthy subjects
aged 18-85 yrs (mean 28$\pm$13 yrs, 569 females). 2) Enhanced Nathan Kline
Institute - Rockland Sample (NKI-Rockland) (393 subjects, mean 35$\pm$20 yrs,
164 females) [15, 21], known for wide age range (4-85 yrs). 3) Normal subjects
(no diseases) from Attention-Deficit Hyperactivity Disorder-200 (ADHD-Normal)
(330 subjects, mean 12$\pm$3 yrs, 166 females) [6], focusing on the younger
group. All datasets are available in the UCLA multimodal connectivity database
[4].

Figure 2: Visualization of Refined Brain Network. One example of the brain
networks refined by D-CoRP, compared to other denoising methods all based on
GCN. (a): original brain network; (b): refined brain networks and their
differences from the original network with various connectivity refinement
methods.

###  3.2 Experiment Settings and Evaluation Metrics

We use two GNN blocks and a multi-layer perceptron (MLP) as the backbone (see
Supplementary for detailed structures). We used pooling methods (TopK and
DiffPool) with a pooling ratio of 0.25.

In section 3.3, to test the adaptability of D-CoRP, we plug it into 1) the
initial part of the model (D-CoRP (Single)), which refines the original brain
network; 2) every GNN block of the model (D-CoRP (Full)), which refines the
original and processed brain network. The $\bm{\Sigma}$ is set as
$\bm{I}_{N^{2}\times N^{2}}$, $\bm{\mu}_{\bm{m}}$ is set as
$\frac{\bm{1}_{N^{2}}}{20}$, where $\bm{1}_{N^{2}}$ is $N\times N$ matrix with
value 1 in each place; $\mathcal{T}=0.01$, and $\beta=1$. In section 3.4, for
D-CoRP: Efficient, $r(\bm{m})$ is set as a fixed probability, 0.05.

All training is performed on an Nvidia A-100 with pytorch 2.1 and pytorch-
geometric package [7]. We choose Adam as the optimizer with a cosine annealing
learning rate $\frac{1}{10}\textrm{cos}(\frac{epoch_{c}}{epoch_{T}})$, where
$epoch_{c}$ stands for the current epoch number counting from 0 and
$epoch_{T}$ stands for the training epoch set as 150. MSE is used as the loss
function. 10-fold cross-validation is applied to all datasets.

Two performance metrics are adopted for evaluation: Mean Absolute Error (MAE)
and Root Mean Square Error (RMSE), i.e.,
$\textrm{MAE}=\frac{\sum_{n=1}^{N}|y_{i}-\hat{y_{i}}|}{N}$,
$\textrm{RMSE}=\sqrt{\frac{\sum_{n=1}^{N}(y_{i}-\hat{y_{i}})^{2}}{N}}$, where
$y_{i}$ is the ground truth and $\hat{y_{i}}$ is the predicted value.

###  3.3 Performance of D-CoRP on Baselines

To test the applicability of D-CoRP, we tested the most widely used GNN
baselines: GCN [11], GAT [22], GIN [28], and GraphSAGE [9]. Due to the
importance of graph pooling methods in brain network research[23, 10], we also
tested TopK pooling and DiffPool[8, 29], known for their capability to distil
subgraphs from brain networks.

The experimental results (Table 1) show that D-CoRP significantly outperforms
baseline models in all comparisons, supporting the effectiveness of D-CoRP. Of
note, with D-CoRP integrated, the backbones exhibit reduced standard
deviations in the assessment metrics, suggesting improved stability. This
improvement is particularly relevant in brain network research due to data
heterogeneity arising from demographics and acquisition. Our results also
reveal that adding D-CoRP to all the GNN blocks, i.e., D-CoRP (Full), however,
does not perform better in most cases, indicating over-compressed information.

|  |  | FCP | NKI-Rockland | ADHD-Normal
---|---|---|---|---|---
Method | Baseline | MAE$\uparrow$ | RMSE$\downarrow$ | MAE$\downarrow$ | RMSE$\downarrow$ | MAE$\downarrow$ | RMSE$\downarrow$
| N/A-5% | GCN | 18.62$\pm$4.27 | 30.22$\pm$7.36 | 55.35$\pm$10.44 | 82.12$\pm$24.45 | 11.42$\pm$5.86 | 22.21$\pm$6.52
| GAT | 21.30$\pm$2.88 | 46.66$\pm$2.98 | 44.43$\pm$7.46 | 62.27$\pm$8.69 | 16.32$\pm$2.42 | 20.78$\pm$3.11
| GIN | 21.80$\pm$2.32 | 25.99$\pm$2.57 | 30.69$\pm$3.46 | 38.33$\pm$4.21 | 9.45$\pm$2.95 | 15.74$\pm$4.01
| SIB [30] | GCN | 27.24$\pm$0.15 | 30.06$\pm$0.20 | 33.96$\pm$1.72 | 39.39$\pm$1.46 | 11.56$\pm$0.22 | 12.01$\pm$0.23
| GAT | 27.41$\pm$0.22 | 30.21$\pm$0.41 | 35.62$\pm$0.92 | 41.18$\pm$1.23 | 11.80$\pm$0.43 | 12.25$\pm$0.98
| GIN | 27.28$\pm$1.25 | 29.81$\pm$1.31 | 35.40$\pm$1.66 | 40.80$\pm$2.53 | 11.23$\pm$1.45 | 11.67$\pm$1.64
| VIB-GSL [18] | GCN | 8.66$\pm$0.40 | 13.90$\pm$0.58 | 16.44$\pm$1.88 | 21.02$\pm$1.63 | 4.66$\pm$1.68 | 5.73$\pm$1.54
| GAT | 11.10$\pm$0.94 | 17.10$\pm$0.85 | 16.69$\pm$1.43 | 22.01$\pm$1.43 | 4.55$\pm$0.76 | 5.55$\pm$0.75
| GIN | 12.54$\pm$0.46 | 17.78$\pm$0.61 | 30.58$\pm$1.71 | 36.33$\pm$2.02 | 8.92$\pm$0.80 | 9.66$\pm$0.79
| D-CoRP: Efficient | GCN | 11.50$\pm$0.96 | 16.66$\pm$0.98 | 34.23$\pm$1.33 | 41.07$\pm$3.17 | 5.57$\pm$1.36 | 10.68$\pm$2.01
| GAT | 12.96$\pm$2.07 | 20.41$\pm$4.65 | 36.71$\pm$1.93 | 41.53$\pm$2.42 | 2.53$\pm$0.41 | 2.96$\pm$0.44
| GIN | 11.76$\pm$1.08 | 16.22$\pm$1.16 | 16.59$\pm$0.77 | 20.00$\pm$1.00 | 9.66$\pm$2.01 | 15.93$\pm$3.08
| D-CoRP | GCN | 8.94$\pm$0.15 | 12.13$\pm$0.28 | 16.14$\pm$0.36 | 20.32$\pm$1.06 | 2.77$\pm$0.23 | 3.45$\pm$0.22
| GAT | 9.31$\pm$0.16 | 13.03$\pm$0.26 | 16.06$\pm$1.50 | 20.18$\pm$1.51 | 2.83$\pm$0.19 | 3.47$\pm$0.23
| GIN | 9.11$\pm$0.18 | 12.24$\pm$0.38 | 15.62$\pm$1.43 | 19.60$\pm$1.66 | 2.73$\pm$0.13 | 3.42$\pm$0.17

Table 2: Comparison of D-CoRP and Other Denoising Methods. The bold font
indicates the global optimal result for each dataset.

###  3.4 Comparison with Other Denoising Methods

We compare D-CoRP and D-CoRP: Efficient with conventional denoising method,
selecting top 5% connections according to correlation values. We also compare
two IB-based graph structure learning methods, i.e., SIB [30], and VIB-GSL
[18]. All comparisons involve three different backbones (i.e., GCN, GAT and
GIN) as the graph processing module.

From Table 2, we observe that D-CoRP outperforms other models in most
comparisons. In all cases, D-CoRP shows better stability over the comparison
models, indicated by lower variance. These results support the effectiveness
of D-CoRP for brain network research. Additionally, our results imply that
D-CoRP: Efficient may be more suitable for smaller datasets (e.g., ADHD-
Normal).

(a) D-CoRP: Efficient

(b) D-CoRP

Figure 3: Results of Variant Temperature.

Further, we visualize an example network derived from the ADHD-Normal dataset
in Fig. (2). The results show that our methods effectively remove noisy
connections compared to other methods. Interestingly, from the refined network
from D-CoRP, we observe the modularity preserved in the functional network,
one of the most important brain network properties, suggesting the D-CoRP’s
capability to preserve informative edges alongside denoising.

###  3.5 Hyperparameter Tuning on D-CoRP

The temperature $\mathcal{T}$ controls the slope of sampling function
$\mathcal{F}(\bm{\mathcal{P}},\bm{\Pi}:\mathcal{T})$, influencing how
effectively D-CoRP refines brain network connectivities. Six backbone models
in Table. 1 are utilized to test the effectiveness of our D-CoRP. 10-fold
cross-validation is performed on the FCP dataset. MAE is used to evaluate the
final performance. Fig. 3 shows $\mathcal{T}$ value of 0.01 for D-CoRP may
lead to superior and stable model performance in more than half of cases and
setting $\mathcal{T}$ as 1 may be generally suitable for all cases. We
recommend setting $\mathcal{T}$ as 0.01 if only binary masking is needed while
setting $\mathcal{T}$ as 1 if D-CoRP is also utilized to adjust the importance
of connections in addition to the denoised part.

##  4 Conclusion

In this study, we propose D-CoRP, a tool for differentiable connectivity
refinement in functional brain networks based on information bottleneck
theory. Specifically, we design a differentiable sampling function and
multivariate probability modeling for detecting valuable connections in the
brain network. Our experiments show that D-CoRP achieves state-of-the-art
performance in connectivity refinement. As a plugin, D-CoRP can be integrated
into existing GNN-based frameworks and provides an efficient mode for flexible
adoption. Future work will develop multivariate optimization with adaptive
covariance tailored to functional connectivity.

## References

  * [1] Alemi, A.A., Fischer, I., Dillon, J.V., Murphy, K.: Deep variational information bottleneck. arXiv preprint arXiv:1612.00410 (2016)
  * [2] Bashyam, V.M., Erus, G., Doshi, J., Habes, M., Nasrallah, I.M., Truelove-Hill, M., Srinivasan, D., Mamourian, L., Pomponio, R., Fan, Y., et al.: Mri signatures of brain age and disease over the lifespan based on a deep brain network and 14 468 individuals worldwide. Brain 143(7), 2312–2324 (2020)
  * [3] Biswal, B.B., Mennes, M., Zuo, X.N., Gohel, S., Kelly, C., Smith, S.M., Beckmann, C.F., Adelstein, J.S., Buckner, R.L., Colcombe, S., et al.: Toward discovery science of human brain function. Proceedings of the national academy of sciences 107(10), 4734–4739 (2010)
  * [4] Brown, J.A., Rudie, J.D., Bandrowski, A., Van Horn, J.D., Bookheimer, S.Y.: The ucla multimodal connectivity database: a web-based platform for brain connectivity matrix sharing and analysis. Frontiers in neuroinformatics 6, 28 (2012)
  * [5] Bullmore, E., Sporns, O.: Complex brain networks: graph theoretical analysis of structural and functional systems. Nature reviews neuroscience 10(3), 186–198 (2009)
  * [6] consortium, A..: The adhd-200 consortium: a model to advance the translational potential of neuroimaging in clinical neuroscience. Frontiers in systems neuroscience 6, 62 (2012)
  * [7] Fey, M., Lenssen, J.E.: Fast graph representation learning with pytorch geometric. arXiv preprint arXiv:1903.02428 (2019)
  * [8] Gao, H., Ji, S.: Graph u-nets. In: international conference on machine learning. pp. 2083–2092. PMLR (2019)
  * [9] Hamilton, W., Ying, Z., Leskovec, J.: Inductive representation learning on large graphs. Advances in neural information processing systems 30 (2017)
  * [10] Hilgetag, C.C., Goulas, A.: ‘hierarchy’in the organization of brain networks. Philosophical Transactions of the Royal Society B 375(1796), 20190319 (2020)
  * [11] Kipf, T.N., Welling, M.: Semi-supervised classification with graph convolutional networks. arXiv preprint arXiv:1609.02907 (2016)
  * [12] Levakov, G., Rosenthal, G., Shelef, I., Raviv, T.R., Avidan, G.: From a deep learning model back to the brain—identifying regional predictors and their relation to aging. Human brain mapping 41(12), 3235–3252 (2020)
  * [13] Li, C., Wei, Y., Chen, X., Schönlieb, C.B.: Brainnetgan: Data augmentation of brain connectivity using generative adversarial network for dementia classification. In: Deep Generative Models, and Data Augmentation, Labelling, and Imperfections: First Workshop, DGM4MICCAI 2021, and First Workshop, DALI 2021, Held in Conjunction with MICCAI 2021, Strasbourg, France, October 1, 2021, Proceedings 1. pp. 103–111. Springer (2021)
  * [14] Li, X., Zhou, Y., Dvornek, N., Zhang, M., Gao, S., Zhuang, J., Scheinost, D., Staib, L.H., Ventola, P., Duncan, J.S.: Braingnn: Interpretable brain graph neural network for fmri analysis. Medical Image Analysis 74, 102233 (2021)
  * [15] Nooner, K.B., Colcombe, S.J., Tobe, R.H., Mennes, M., Benedict, M.M., Moreno, A.L., Panek, L.J., Brown, S., Zavitz, S.T., Li, Q., et al.: The nki-rockland sample: a model for accelerating the pace of discovery science in psychiatry. Frontiers in neuroscience 6, 152 (2012)
  * [16] Park, H.J., Friston, K.: Structural and functional brain networks: from connections to cognition. Science 342(6158), 1238411 (2013)
  * [17] Stanley, M.L., Simpson, S.L., Dagenbach, D., Lyday, R.G., Burdette, J.H., Laurienti, P.J.: Changes in brain network efficiency and working memory performance in aging. PLoS One 10(4), e0123950 (2015)
  * [18] Sun, Q., Li, J., Peng, H., Wu, J., Fu, X., Ji, C., Philip, S.Y.: Graph structure learning with variational information bottleneck. In: Proceedings of the AAAI Conference on Artificial Intelligence. vol. 36, pp. 4165–4174 (2022)
  * [19] Tishby, N., Pereira, F.C., Bialek, W.: The information bottleneck method. arXiv preprint physics/0004057 (2000)
  * [20] Tishby, N., Zaslavsky, N.: Deep learning and the information bottleneck principle. In: 2015 ieee information theory workshop (itw). pp. 1–5. IEEE (2015)
  * [21] Tobe, R.H., MacKay-Brandt, A., Lim, R., Kramer, M., Breland, M.M., Tu, L., Tian, Y., Trautman, K.D., Hu, C., Sangoi, R., et al.: A longitudinal resource for studying connectome development and its psychiatric associations during childhood. Scientific Data 9(1), 300 (2022)
  * [22] Veličković, P., Cucurull, G., Casanova, A., Romero, A., Lio, P., Bengio, Y.: Graph attention networks. arXiv preprint arXiv:1710.10903 (2017)
  * [23] Vidaurre, D., Smith, S.M., Woolrich, M.W.: Brain network dynamics are hierarchically organized in time. Proceedings of the National Academy of Sciences 114(48), 12827–12832 (2017)
  * [24] Wang, P.Y., Sapra, S., George, V.K., Silva, G.A.: Generalizable machine learning in neuroscience using graph neural networks. Frontiers in artificial intelligence 4, 618372 (2021)
  * [25] Wei, Y., Chen, X., Zhu, L., Zhang, L., Schönlieb, C.B., Price, S., Li, C.: Multi-modal learning for predicting the genotype of glioma. IEEE Transactions on Medical Imaging (2023)
  * [26] Wei, Y., Li, C., Price, S.J.: Quantifying structural connectivity in brain tumor patients. In: Medical Image Computing and Computer Assisted Intervention–MICCAI 2021: 24th International Conference, Strasbourg, France, September 27–October 1, 2021, Proceedings, Part VII 24. pp. 519–529. Springer (2021)
  * [27] Williams, N., Arnulfo, G., Wang, S.H., Nobili, L., Palva, S., Palva, J.M.: Comparison of methods to identify modules in noisy or incomplete brain networks. Brain connectivity 9(2), 128–143 (2019)
  * [28] Wu, T., Ren, H., Li, P., Leskovec, J.: Graph information bottleneck. Advances in Neural Information Processing Systems 33, 20437–20448 (2020)
  * [29] Ying, Z., You, J., Morris, C., Ren, X., Hamilton, W., Leskovec, J.: Hierarchical graph representation learning with differentiable pooling. Advances in neural information processing systems 31 (2018)
  * [30] Yu, J., Xu, T., Rong, Y., Bian, Y., Huang, J., He, R.: Graph information bottleneck for subgraph recognition. arXiv preprint arXiv:2010.05563 (2020)
  * [31] Zhang, P., Wen, G., Cao, P., Yang, J., Zhang, J., Zhang, X., Zhu, X., Zaiane, O.R., Wang, F.: Brainusl: U nsupervised graph s tructure l earning for functional brain network analysis. In: International Conference on Medical Image Computing and Computer-Assisted Intervention. pp. 205–214. Springer (2023)
  * [32] Zhou, J., Cui, G., Hu, S., Zhang, Z., Yang, C., Liu, Z., Wang, L., Li, C., Sun, M.: Graph neural networks: A review of methods and applications. AI open 1, 57–81 (2020)
  * [33] Zhu, Y., Xu, W., Zhang, J., Liu, Q., Wu, S., Wang, L.: Deep graph structure learning for robust representations: A survey. arXiv preprint arXiv:2103.03036 14 (2021)


# Yet another test of Radial Acceleration Relation for galaxy clusters

S. Pradyumna1<EMAIL_ADDRESS>Sajal Gupta2
<EMAIL_ADDRESS>Sowmya Seeram3
<EMAIL_ADDRESS>Shantanu Desai4 E-mail<EMAIL_ADDRESS>1,3,4Department of Physics, Indian Institute of Technology, Hyderabad,
Telangana-502285, India  2Department of Physical Sciences, IISER-Kolkata,
Mohanpur, West Bengal-741246, India

###### Abstract

We carry out a test of the radial acceleration relation (RAR) for galaxy
clusters from two different catalogs compiled in literature, as an independent
cross-check of two recent analyses, which reached opposite conclusions. The
datasets we considered include a Chandra sample of 12 clusters and the X-COP
sample of 12 clusters. For both the samples, we find that the residual scatter
is small (0.11-0.14 dex), although the best-fit values for the Chandra sample
have large error bars. Therefore, we argue that at least one of these cluster
samples (X-COP) obeys the radial acceleration relation. However, since the
best-fit parameters are discrepant with each other as well as the previous
estimates, we argue that the RAR is not universal. For both the catalogs, the
acceleration scale, which we obtain is about an order of magnitude larger than
that obtained for galaxies, and is agreement with both the recent estimates.

##  I Introduction

Recently, two independent groups Chan and Del Popolo (2020); Tian et al.
(2020a) carried out a test of the correlation between baryonic and total
acceleration for two different cluster samples. This was motivated by the
discovery of a strong deterministic relation, between the baryonic ($a_{bar}$)
and total acceleration ($a_{tot}$) using 2693 data points from 153
rotationally supported spiral galaxies of different sizes and morphology, with
a small scatter of 0.13 dex McGaugh et al. (2016); Lelli et al. (2016). This
scatter is dominated by uncertainties in galaxy distance, disk inclination and
the stellar mass-to-light ratio. Once these are marginalized over, the
residuals around the fits have negligible scatter (of about 0.057 dex) Li et
al. (2018). Note however that this aforementioned scatter was obtained by
assuming that every galaxy obeys the mean RAR, followed by doing a separate
fit for each of the 175 galaxies, followed by marginalizing over the
aforementioned uncertainties. A similar relation has also been recently found
for elliptical galaxies (Chae et al., 2020). This relation has been dubbed as
the radial acceleration relation (RAR) and can be written as follows:

$a_{tot}=\frac{a_{bar}}{1-e^{-\sqrt{a_{bar}/a_{0}}}},$ (1)

where $a_{0}\sim 1.2\times 10^{-10}m/s^{2}$ McGaugh et al. (2016). This
equation is closely related to the Mass Discrepancy-Acceleration relation
McGaugh (2004), and is a trivial consequence of the MOND paradigm Famaey and
McGaugh (2012); Milgrom (2016), and also subsumes the baryonic Tully-Fisher
relation Lelli et al. (2017a); Wheeler et al. (2019). We note that a few
groups have have questioned the existence of a fundamental acceleration scale
in the same data Rodrigues et al. (2018a); Zobnina and Zasov (2020); Marra et
al. (2020); Zhou et al. (2020). Although, some of these criticisms have been
countered McGaugh et al. (2018); Kroupa et al. (2018), the debate is still on-
going on some of these issues Rodrigues et al. (2018b); Rodrigues et al.
(2020). There are mixed claims in the literature on whether this relation can
be easily postdicted using the standard $\Lambda$CDM paradigm Desmond (2017);
Navarro et al. (2017); Dai and Lu (2017); Keller and Wadsley (2017); Garaldi
et al. (2018); Ludlow et al. (2017); Stone and Courteau (2019), with the same
negligible scatter as observed for real data. This relation can also be
reproduced using alternative models such as Self-Interacting Dark matter Ren
et al. (2019), superfluid dark matter Berezhiani et al. (2018), Moffat’s MOG
Green and Moffat (2019), conformal gravity O’Brien et al. (2018), but is in
tension with Verlinde’s entropic gravity Lelli et al. (2017b).

Motivated by these considerations, two independent groups Chan and Del Popolo
(2020); Tian et al. (2020a) carried out a systematic observational test of the
RAR using galaxy clusters. Galaxy clusters are the most massive collapsed
objects in the universe, and have proved to be a wonderful laboratory for a
wide range of topics from cosmology, tests of modified gravity, fundamental
Physics (such as constraining neutrino and graviton mass) to galaxy evolution
Voit (2005); Vikhlinin et al. (2014); Allen et al. (2011); Desai (2018).
However, it has been known since a long time, that MOND does not work well for
both relaxed as well as merging galaxy clusters The and White (1988); Gerbal
et al. (1992); Aguirre et al. (2001); Sanders (2003); Pointecouteau and Silk
(2005); Angus et al. (2008); Natarajan and Zhao (2008), and is also at odds
with shape of the matter power spectrum Dodelson (2011). Many relativistic
theories of MOND are also in tension due to the coincident gravitational wave
and electro-magnetic observations from GW170817 (Boran et al., 2018).

Nevertheless, there was no detailed characterization of the RAR for clusters,
prior to the two works Chan and Del Popolo (2020); Tian et al. (2020a).
Despite the failures of MOND on cluster scales, a detailed characterization of
the RAR for a large ensemble of galaxy clusters would be an acid test for some
of the myriad models used to explain the same relation for galaxies. However,
there are also no clear cut predictions for whether the RAR should be
satisfied or not for galaxy clusters, except for a few works Famaey et al.
(2018); Navarro et al. (2017), which predict that this relation should not be
seen for galaxy clusters. We briefly discuss the results from these two
observational tests of RAR for clusters.

Chan and Del Popolo Chan and Del Popolo (2020) (CDP20, hereafter) considered
52 non-cool-core clusters in hydrostatic equilibrium from the HIFLUGCS X-ray
sample, based on ROSAT and ASCA observations Chen et al. (2007), with core-
radius ($r_{c}$) greater than 100 kpc. The density of the hot gas was modeled
using a single-$\beta$ profile Cavaliere and Fusco-Femiano (1978) and a
constant temperature profile was assumed. Then assuming hydrostatic
equilibrium, the total mass was calculated assuming spherical symmetry Allen
et al. (2011). From the total mass, the total acceleration was estimated
assuming Newtonian gravity Chan and Del Popolo (2020). Similarly, the baryonic
acceleration was obtained from the baryonic mass (estimated by integrating the
gas density profile). We note that the stellar mass was not included in their
baryonic mass budget. CDP20 then checked for a correlation between the
baryonic and total acceleration at four different positions ($r_{c}$,
$2r_{c}$, $3r_{c}$, $r_{500}$) for their sample. They found that the best-fit
value for $a_{0}$ was $9.5\times 10^{-10}m/s^{2}$, about ten times higher than
the same found for spiral galaxies. Furthermore, their estimated residual
scatter which they obtained when compared to the RAR is about 0.18 dex, and
much higher than that for galaxies.

Around the same time, Tian et al Tian et al. (2020a) (T20, hereafter) tested
this relation for a sample of 20 high-mass clusters from the CLASH survey
Postman et al. (2012). These clusters have a combination of strong-lensing,
weak lensing shear, and magnification data using HST and Suprime-Cam
observations Umetsu et al. (2016); and gas temperatures and density profiles
from both Chandra and XMM observations Donahue et al. (2014). The total mass
was modeled as a sum of four distinct components: the total dark matter, X-ray
emitting gas, non-BCG stellar mass and BCG stellar mass. The baryonic mass was
obtained from the sum of the last three components. The dark matter mass was
obtained from the NFW Navarro et al. (1997) density profile and mass-
concentration relations from Umetsu et al. (2016). The X-ray gas mass was
obtained from Donahue et al. (2014). The non-BCG stellar mass was obtained
from the stellar mass fraction derived using the SPT sample Chiu et al.
(2018). The BCG stellar mass was estimated by assuming an Hernquist profile
Hernquist (1990) for the stellar mass distribution and assuming the the
stellar mass estimates from Ref. Cooke et al. (2016). They looked for RAR by
transforming $a_{tot}$ and $a_{bar}$ into log-log space, thereby converting it
to a linear regression problem. They obtained an extremely tight scaling
relation between the total and baryonic acceleration, with an intrinsic
scatter of about 14%. This is comparable to the scatter obtained for the SPARC
sample McGaugh et al. (2016). However, the best-fit value for $a_{0}$ was
$2\times 10^{-9}m/s^{2}$, about 20 times higher. They also explained the
observed scatter using a semi-analytic galaxy formation model within
$\Lambda$CDM Olamaie et al. (2012). Therefore, although both of these works
get about the same value for the acceleration scale, they disagree on whether
the RAR exists for galaxy clusters with negligible scatter. We note that most
recently, this group has also found evidence for a tight Baryonic Faber-
Jackson relation for the HIFLUGCS cluster sample Tian et al. (2020b)

Given the contradictory findings between the two groups, we carry out a
similar test of the RAR for two different cluster samples. The first dataset
we analyze is a sample of 12 galaxy clusters with pointed Chandra and archival
ROSAT observations which were carefully studied by Vikhlinin et al. (2005,
2006). This dataset has been used for testing a plethora of modified gravity
theories (Rahvar and Mashhoon, 2014; Hodson and Zhao, 2017; Hodson et al.,
2017; Edmonds et al., 2018, 2017; Bernal et al., 2017). We have also
previously used this sample to constrain the graviton mass, test of constant
halo surface density, and GR corrections to hydrostatic masses Gupta and Desai
(2019, 2020); Gopika and Desai (2020).

We now provide a brief clarification of what we mean by a test of RAR for
clusters. Although the moniker “RAR” (for spiral galaxies) usually refers to
the empirical relation in Eq. 1 with the same value of $a_{0}$ as found in
McGaugh et al. (2016), one can get very similar deterministic relations
between $a_{tot}$ and $a_{bar}$ for other proposed MOND interpolating
functions Famaey and McGaugh (2012). A test of these relations for spiral
galaxies would also yield a very tight scatter, although the difference
between these different functions could be discerned using the Halo
Acceleration relation (Tian and Ko, 2019). Similarly, it is possible that for
galaxy clusters, Eq. 1 could fit the data, albeit with a higher value for
$a_{0}$, similar to the value for $a_{0}$ found in the first ever test of MOND
with clusters The and White (1988). To carry out a generic model-independent
test, we note that Eq. 1 can be re-written as a linear regression equation
between the logarithm of the total and baryonic acceleration (cf. Eq. 9). This
is the same fit done in T20. Therefore, for the rest of this manuscript by
“test of RAR”, implies a test of whether there exists a tight empirical
scaling relation between the baryonic and total accelerations in logarithmic
space, throughout the galaxy cluster. We shall also characterize the scatter
in the residuals between the two accelerations, as well as determine the
intrinsic scatter in the scaling relations. A separate question is whether
this relation is universal, or in other words the different cluster samples
yield the same best-fit values. We shall also test for the universality of
this relation for our cluster samples.

##  II Basic considerations

We review the Physics, common to all galaxy cluster samples, needed to
estimate the total acceleration as well as the baryonic acceleration.

In the first step, we calculate the total mass from equation of hydrostatic
equilibrium, which can be applied if the clusters are relaxed. The systematic
uncertainty in this assumption is about 15-20% Biffi et al. (2016). If the
parametric forms for the gas temperature and gas density profiles are given,
the total acceleration at given radius $r$ from the cluster center can be
derived from positing an ideal gas equation of state Allen et al. (2011):

$a_{tot}(r)=-\frac{k_{b}T}{r\mu m_{p}}\left(\frac{d\ln\rho_{gas}}{d\ln
r}+\frac{d\ln T}{d\ln r}\right),$ (2)

where $\rho_{gas}$ denotes the gas mass density, $T$ is the cluster
temperature, $m_{p}$ is the mass of the proton, and $\mu$ is the mean
molecular weight of the cluster in a.m.u. and is approximately equal to 0.6
Vikhlinin et al. (2005); Rahvar and Mashhoon (2014). Alternately, if the total
total mass ($M_{tot}$) is known, then the total acceleration is given by:

$a_{tot}=\frac{GM_{tot}}{r^{2}}$ (3)

In order to obtain the baryonic acceleration, we need to estimate the total
mass in baryons, which consists of hot diffuse gas in the intra-cluster medium
and stars. The gas mass is obtained by integrating the gas density profile:

$M_{gas}=\int 4\pi r^{2}\rho_{gas}(r)dr$ (4)

We note that the aforementioned estimates of gas mass and the total mass
assume spherical symmetry. Errors due to this assumption can be upto 5% Gopika
and Desai (2020). There are multiple methods to obtain the star mass in galaxy
clusters Lin et al. (2012); Chiu et al. (2018). We use the method in Lin et
al. (2012). This entails calculating the total mass ($M_{500}$) within a
radius with an overdensity of 500, called $r_{500}$. We then use the following
empirical relation to obtain the star mass at $r=r_{500}$ Lin et al. (2012).

$\frac{M_{star}(r=r_{500})}{10^{12}M_{\odot}}=1.8\left(\frac{M_{500}}{10^{14}M_{\odot}}\right)^{0.71}$
(5)

This equation yields the star mass at $r=r_{500}$. The star mass at any other
radius can be obtained by assuming an isothermal profile Rahvar and Mashhoon
(2014). We note that since the stellar mass contribution is negligible, errors
due to stellar mass estimate should not dominate the total budget and are
neglected.

The baryonic acceleration is then given by

$a_{bar}=\frac{G(M_{gas}+M_{stars})}{r^{2}}$ (6)

##  III Chandra X-ray cluster sample

Vikhlinin et al Vikhlinin et al. (2005, 2006) provided gas density and
temperature profiles for 13 different clusters, using measurements from the
pointed or archival observations with the Chandra and ROSAT X-ray satellites.
For our analysis we use 12 of these clusters, viz. A133, A262, A383, A478,
A907, A1413, A1795, A1991, A2029, A2390, MKW4, RXJ1159+55531, since errors for
one of them were not available. The redshifts of these clusters range
approximately upto $z=0.2$. These measurements extended up to very large radii
of about $r_{500}$ for some of the clusters. The typical exposure times ranged
from 30-130 ksecs. The temperatures span the range between 1 and 10 keV and
masses from $(0.5-10)\times 10^{14}M_{\odot}$. More details of the
observations, data reduction, and estimation of temperature and density
profiles have been discussed in Vikhlinin et al. (2006). This dataset has been
used for testing a plethora of modified gravity theories (Rahvar and Mashhoon,
2014; Hodson and Zhao, 2017; Hodson et al., 2017; Edmonds et al., 2018, 2017;
Bernal et al., 2017; Gupta and Desai, 2019). We note that although the first
test of RAR for this sample was done in Edmonds et al. (2018), no error
analysis in the observed accelerations was taken into account in that
analysis. They considered 250 radii points starting at 10 pc and extending to
2.5 kpc in increments of 10 pc (D. Edmonds, private communication). Also the
stellar mass was not included in their estimation of baryonic mass. Based on a
visual inspection of the scatter plot between $a_{tot}$ and $a_{bar}$, Edmonds
et al. (2018) asserted that there is a large scatter about the RAR. However,
since no error bars were included in that work, no quantitative estimate of
the intrinsic scatter was done.

The gas density ($\rho_{g}$) model used for this sample is the double-beta
profile Vikhlinin et al. (2006) and is equal to
$1.624m_{p}\sqrt{n_{p}(r)n_{e}(r)}$, where $m_{p}$ is the proton mass and the
product of electron and proton densities ($n_{e}n_{p}$) is given by the
double-$\beta$ profile:

$n_{e}(r)n_{p}(r)=\frac{(r/r_{c})^{-\alpha^{\prime}}}{(1+r^{2}/r_{c}^{2})^{3\beta-\alpha^{\prime}/2}}\frac{n_{0}^{2}}{(1+r^{\gamma}/r_{s}^{\gamma})^{\epsilon/\gamma}}+\frac{n_{02}^{2}}{(1+r^{2}/r_{c^{\prime}}^{2})^{3\beta^{\prime}}}.$
(7)

The explanation and values for each of the terms in Eq. 7 can be found in the
original works by Vikhlinin Vikhlinin et al. (2009). The temperature profile
is given by:

$T(r)=T_{0}\frac{(x_{0}+T_{min}/T_{0})}{x_{0}+1}\frac{(r/r_{t})^{-a^{\prime}}}{\left[1+(r/r_{t})^{b}\right]^{c^{\prime}/b}},$
(8)

where $x_{0}=\left(\frac{r}{r_{cool}}\right)^{a_{cool}}$. The physical
explanation of the eight free parameters and their values for the clusters can
be found in Vikhlinin et al. (2006); Rahvar and Mashhoon (2014). The observed
values of $T(r)$ along with their error bars at various points from the
cluster center were provided for the 12 clusters by A. Vikhlinin (private
communication).

We now plug in the values for $\rho(r)$ and $T(r)$ from Eqns. 7 and 8
respectively in Eq. 2 to calculate the total acceleration $a(r)$. The errors
were obtained by propagating the $T(r)$ errors provided in V06. We then
estimate the baryonic mass from the sum of gas mass (obtained from Eq. 4) and
star mass (obtained from Eq. 5).

Therefore, from Eqns. 2 and 6, $a_{tot}$ and $a_{bar}$ can be calculated at
any distance $r$ from the cluster center. As mentioned earlier, CDP20 have
done this test at different radii for the various clusters at integer
multiples of $r_{c}$, whereas T20 have stacked data from the CLASH cluster
sample (of about 10-18 clusters) at four different radii, with a total of 64
datapoints.

Here, we do a joint fit to the RAR using the data for $a_{tot}$ and $a_{bar}$
using data at three different radii: 200 kpc, 500 kpc, and 1000 kpc . Similar
to T20, we do a fit in log-acceleration space by doing a linear regression to
the relation:

$y=mx+b$ (9)

where $y=\ln a_{tot}$ and $x=\ln a_{bar}$. We then maximize a log-likelihood
given by :

$\displaystyle-2\ln L$ $\displaystyle=$ $\displaystyle\large{\sum_{i}}\ln
2\pi\sigma_{i}^{2}+\large{\sum_{i}}\frac{[y_{i}-(mx_{i}+b)]^{2}}{\sigma_{i}^{2}}$
(10) $\displaystyle\sigma_{i}^{2}$ $\displaystyle=$
$\displaystyle\sigma_{y_{i}}^{2}+m^{2}\sigma_{x_{i}}^{2}+\sigma_{int}^{2}$
(11)

Here, $\sigma_{y_{i}}$ and $\sigma_{x_{i}}$ denote the errors in $\ln a_{tot}$
and $\ln a_{bar}$ respectively; $\sigma_{int}$ is the intrinsic scatter for
our linear relation between the logarithms of the accelerations, which is a
free parameter. We used the emcee MCMC sampler Foreman-Mackey et al. (2013) to
maximize Eq. 10. Our best-fit values for $m$ and $b$ are shown in Table 1,
along with the corresponding fits from CDP20 and T20. The 68%, 90% and 99%
marginalized confidence intervals for $m$, $b$, and $\sigma_{int}$ can be
found in Fig. 2. Although the intrinsic scatter is small, the error on the
intercept is quite large. As we can see, the 1-D marginalized credible
interval for $\ln(\sigma_{int})$ remains flat. We also verified that alternate
samplers such as Nestle111http://kylebarbary.com/nestle/ and analytical
optimizaion routines return the same best-fit as emcee to within $1\sigma$
error bars, as well as similar credible intervals. As the contribution to the
total scatter from the best-fit value of $\sigma_{int}=2\times 10^{-4}$ is
negligible, smaller values for $\sigma_{int}$ will not change the likelihood.
That explains the flatness of the marginalized distribution of
$\ln(\sigma_{int})$. Our slope and intercept are discrepant compared to T20 at
about $2.3\sigma$ and $2.2\sigma$ respectively.

In order to compare the tightness of our RAR relation for this sample, we also
calculate the residual scatter in dex to directly compare with the results in
CDP20 and T20 (and also McGaugh et al. (2016); Li et al. (2018)), who used a
similar procedure to characterize the scatter. The residuals are given by
$\log_{10}(a_{tot})-\log_{10}(a_{exp})$, where $a_{exp}$ is the expected total
acceleration obtained from the best-fit RAR. These binned residuals (using
four bins) are then fitted to a Gaussian and the standard deviation of this
Gaussian distribution gives the residual scatter in dex. When we fit these
residuals for the Chandra sample to a Gaussian distribution, we get
$\sigma=0.13$ dex. However, we find that the residuals are asymmetric with
longer tails at negative values for the residuals. Therefore, instead of the
Gaussian $\sigma$, we use the scale parameter obtained using the normalized
median absolute deviation, $\sigma_{MAD}$, (computed using the astropy Astropy
Collaboration et al. (2018) library) to parametrize the residual scatter. This
value is given by $\sigma_{MAD}=0.14$. Therefore, the residual scatter for the
Chandra sample is equal to 0.14 dex.

We overlay our best-fit values along with the same from T20 and CDP20, as well
as the RAR found from the SPARC sample McGaugh et al. (2016) in Fig. 1.
Because of the large errors, both T20 and CDP20 fits pass through the data
points. However, despite the large errors, the RAR obtained from SPARC sample
is below most of the data points.

To determine the acceleration scale $a_{0}$, similar to T20, we fix the slope
in the regression relation in Eq. 9 to be equal to 0.5. This way the RAR
relation (Eq. 1) in the low acceleration limit becomes
$a_{dyn}=\sqrt{a_{0}a_{bar}}$. From this relation it is trivial to derive the
Baryonic Tully-Fisher or Faber-Jackson relations Famaey and McGaugh (2012).
The resulting value of $a_{0}$ is equal to $(9.26\pm 1.66)\times
10^{-10}m/s^{2}$. Therefore, similar to T20 and CDP20, the value of $a_{0}$ is
about an order of magnitude larger than obtained using galaxies McGaugh et al.
(2016). Our estimated value for $a_{0}$ agrees within our $1\sigma$ errors
with the estimate in CDP20.

Figure 1: A test of the radial acceleration relation for 12 clusters from the
Chandra sample (black dots) Vikhlinin et al. (2006) using data stacked at 200
kpc, 500 kpc, and 1000 kpc. We also overlay our best-fit along with that from
T20, CDP20 as well as the best-fit RAR deduced from the SPARC sample McGaugh
et al. (2016). We see that the SPARC-based RAR does not agree with the Chandra
data. Figure 2: Plot showing 68%, 90% and 99% marginalized credible intervals
for the slope, intercept and (natural log) intrinsic scatter for the linear
relation between $\ln a_{tot}$ and $\ln a_{bar}$ for the Chandra cluster
sample.

##  IV X-COP sample

X-COP (XMM-Newton Cluster Outskirts Project) (PI: D.Eckert) is a very large
program on XMM-Newton dedicated to studying X-ray emission in cluster
outskirts. The X-COP project targeted a sample of massive 12 galaxy clusters
($M_{500}>3\times 10^{14}M_{\odot}$) in the redshift range of $0.04-0.1$
selected from the Planck all-sky Sunyaev-Zel’dovich (SZ) catalog Planck
Collaboration et al. (2014). These selected clusters had a signal-to-noise
ratio cutoff of 12. The detailed cluster selection criterion is outlined in
Ghirardini et al. (2019). The final clusters included in this project are -
A85, A644, A1644, A1795, A2029, A2142, A2255, A2319, A3158, A3266, RXC1825,
and ZW1215. Eckert et al. (2017). The main goal of this project is to
reconstruct the ICM properties upto $R_{200}$. More details on the XCOP sample
and associated science results using this sample can be found in Ettori et al.
(2019); Ghirardini et al. (2019, 2018); Ettori et al. (2020); Ghizzardi et al.
(2020a). This cluster sample has also been used to test modified theories of
gravity, which dispense with dark matter such as MOND and Verlinde’s emergent
gravity Ettori et al. (2019).

For all the clusters in the X-COP sample, the gas density has been obtained
from the deprojected XMM surface brightness profile. The temperature is
obtained from the deprojected XMM spectral estimates. The SZ (Sunyaev-
Zel’dovich) Sunyaev and Zeldovich (1970) effect provides a direct measurement
of the thermal electron pressure integrated along the line of sight. The
pressure profile is recovered from the SZ signal measured in the all-sky
survey by the Planck SZ survey. Planck Collaboration et al. (2016); Ghirardini
et al. (2018). The gas pressure has a smooth spatial distribution along the
azimuth except for the clusters which are undergoing a merger. From the ideal
gas relation, one can deduce the pressure of the gas, when only spectral
measurements are available. Similarly the temperature can be obtained using
pressure obtained from deconvolution of Planck measurements, when spectral
measurements are not available Ghirardini et al. (2019). Parametric
temperature and pressure profiles for the XCOP samples can be found in
Ghirardini et al. (2019).

Mass profiles have been calculated using multiple methods for the X-COP
sample, as documented in Ghirardini et al. (2018). The total mass and gas mass
profiles, obtained using the backward method (with NFW model as reference),
can be downloaded from the web, and are readily available in .fits files for
all the clusters222X-COP Data: https://dominiqueeckert.wixsite.com/xcop/data.
The relative errors in mass measurements were less than 8% Ettori et al.
(2019). For our work, these were used to obtain the total mass. To obtain the
gas mass, the data in the files on the aforementioned website titled Gas mass
and gas fraction profile were used. The gas density can be obtained from
$n_{e}$ using $\rho_{gas}={\mu}{m_{p}}(n_{e}+n_{p})$, with $n_{e}=1.17n_{p}$
Ettori et al. (2019); Eckert et al. (2019). The electron number density
($n_{e}$), as mentioned earlier has been estimated from the deprojected XMM
surface brightness profile.

The stellar mass profiles are available for seven of the galaxy clusters
Ghizzardi et al. (2020b). Accordingly, we used them for these seven clusters,
whereas for the remaining five clusters, Eq. 5 has been used to obtain the
stellar mass. In order to obtain the (baryonic as well as total) mass
estimates at any arbitrary radius, we interpolate the different quantities
using quadratic spline interpolations, and then calculated their values at
100, 200, 400, and 1000 kpc. Then, $a_{bar}$ is obtained using Eq. 6 and
$a_{tot}$ is obtained using Eq. 3. The XCOP website also provides $1\sigma$
errors for all the measurements, which we include for our analysis. We note
that these errors do not account for non-thermal contributions to gas pressure
and departures from spherical symmetry and hydrostatic equilibrium Ettori et
al. (2019).

To find the best-fit values for the linear regression relation between
$\ln(a_{tot})$ and $\ln(a_{bar})$, we maximize the same likelihood as in Eq
10. Our best-fit values are shown in Table 1. The 68%, 90%, and 99% confidence
intervals for the slope, intercept, and natural log of intrinsic scatter are
shown in Fig. 4. Both the slope and intercept differ from T20 by about
$6.7\sigma$. We find an intrinsic scatter of about 21%, which is of the same
order of magnitude as found in T20. We also computed the residual scatter
using the same method as in Sect. III, by fitting a Gaussian distribution to
difference in $\log_{10}$ of the total and baryonic accelerations and
computing its standard deviation. For computation of the residual scatter for
XCOP clusters, we binned the residuals into four bins. Then, residual scatter
was calculated as the standard deviation of the best fit gaussian to the
binned residuals. We find that residual scatter computed using this method is
equal to 0.11 dex, which is also comparable to T20. Therefore, we find that
the X-COP sample obeys the RAR.

The best-fit RAR for our cluster sample along with those found in T20 and C20
are shown in Fig. 3. We also show that the RAR found using the cluster sample.
If we fix the slope to 0.5, the best-fit value of $a_{0}$ is equal to
$(1.12\pm 0.11)\times 10^{-9}m/s^{2}$. This value is mid-way between the
estimates in CDP20 and T20, and is also about an order of magnitude larger
than that for galaxies.

Figure 3: The best-fit straight line (green) for $a_{tot}$ vs $a_{bar}$ of the XCOP set of 12 clusters. The dots indicate the data points. The $a_{bar}$ and $a_{tot}$ have been plotted at 100, 200, 400, 1000 kpc for each of the clusters. The RAR fits from McGaugh et al McGaugh et al. (2016), T20, CDP20 have been overlaid for reference. The RAR deduced from the SPARC sample (McGaugh et al., 2016) also does not agree with the X-COP data. Figure 4: Plot showing 68%, 90% and 99% marginalized credible intervals for the slope, intercept and (ln) intrinsic scatter for the linear relation between $\ln a_{tot}$ and $\ln a_{bar}$ for the X-COP cluster sample. Cluster Sample | Slope | Intercept | Intrinsic Scatter | Residual Scatter (dex) | $a_{0}(m/s^{2})$
---|---|---|---|---|---
Chandra (this work) | $0.77\pm 0.1$ | $-3.5^{+2.6}_{-2.7}$ | ($0.0002\pm 0.018$)% | 0.14 | $(9.26\pm 1.66)\times 10^{-10}$
X-COP (this work) | $1.09_{-0.07}^{+0.07}$ | $4.50_{-1.78}^{+1.79}$ | $0.21_{-0.02}^{+0.03}$ | 0.11 | $(1.12\pm 0.11)\times 10^{-9}$
CLASH Tian et al. (2020a) | $0.51^{0.04}_{-0.05}$ | $-9.8^{+1.07}_{-1.08}$ | $0.147^{+0.029}_{-0.028}$ | 0.11 | $(2.02\pm 0.11)\times 10^{-9}$
HIFLUGCS | $0.72_{-0.06}^{+0.06}$ | $-5.49_{-1.39}^{+1.39}$ | $0.27_{-0.02}^{+0.02}$ | 0.13 | $9.5\times 10^{-10}$
Table 1: Summary of our results for linear regression between $\ln(a_{tot})$
as a function of $\ln(a_{bar})$ (Eq. 9), for Chandra and X-COP samples
analyzed, as well as our re-analysis of the HIFLUGCS sample. Also shown for
comparison are the results from T20 Tian et al. (2020a). CDP20 obtain a higher
residual scatter of 0.18 dex for this sample compared to our estimate of 0.13
dex. The value of $a_{0}$ for the first three samples was obtained by assuming
that the slope in the linear regression relation is equal to 0.5. For the
HIFLUGCS sample, only $a_{0}$ shown here is from CDP20.

##  V Comparison of results

In Sects. III and IV, we carried out a test of RAR in the same way as in T20
for the Chandra and X-COP cluster samples, respectively by fitting to Eq. 9.
We have already compared our results for the slope and intercept of these
samples to those for the CLASH sample obtained in T20. However, a direct
comparison of our results with the analysis in CDP20 cannot be trivially done,
since CDP20 used the original RAR relation for spiral galaxies (Eq. 1),
instead of Eq. 9.

Therefore, in order to bring the results from all the cluster samples to a
common footing, we reanalyze the same dataset as in CDP20, using the same
procedure as for the other samples. We select a sample of 52 non-cool core
clusters from the dataset analyzed in Chen et al. (2007), by applying the
condition that the cooling time, $t_{cool}>14$ Gyr (Chen et al., 2007) and
$r_{c}\geq 100r_{c}$, similar to CDP20. For the baryonic mass, we also
included the stellar mass using Eq. 5. Then, we evaluated $a_{tot}$ and
$a_{bar}$ at 100, 200, 400, and 1000 kpc. The best-fit slope, intercept, and
intrinsic scatter were then determined using Eq. 9. Our best-fit values for
the slope, intercept, and intrinsic scatter for this sample are shown in Table
1. The residual scatter, which we get for this sample is 0.13 dex. Note
however that $a_{0}$ in this table for this sample is the same as the estimate
in CDP20. We get an intrinsic scatter of about 27%. Our slope and intercept
for the Chandra sample is consistent with the HIFLUGCS sample to within
$1\sigma$. However, the X-COP slope and intercept are discrepant with respect
to our re-estimate of the HIFLUGCS sample at about 4$\sigma$ and 4.3$\sigma$
respectively.

Therefore, we find that the X-COP sample and the CLASH sample analyzed in T20
show comparable scatter, indicating that RAR is obeyed for the X-COP sample.
However, the HIFLUGS sample (analyzed in CDP20) and the Chandra samples
(analyzed in this work) either show a high intrinsic scatter in the RAR, or
have large error bars for the slope and the intercept. However the slopes and
intercepts for the X-COP and HIFLUGS sample are discrepant with respect to
each other, which implies that RAR is not universal for clusters.

We now discuss, if the discrepancy between the CDP20 result and the other
analyses could be attributed to different assumptions. CDP20 did not include
the stellar mass in the baryonic mass budget. However, our re-analysis of
CDP20, which includes the stellar mass contributions shows a residual scatter
comparable to the CLASH and X-COP sample. Another difference between the T20
and the other three analyses is that T20 included the BCG stellar mass
contribution. However, the contribution of BCG stellar mass is dominant only
at distances close to the cluster center, whereas the RAR analysis for the
other three samples is done at distances far from the center, where the BCG
stellar mass contribution is less than 1% (Tian et al., 2020a; Umetsu et al.,
2016).

Nevertheless, T20 have found that the total stellar mass is important for the
CLASH cluster sample within 200 kpc. If only the gas mass is used for their
analysis, then the slope of the RAR is different from the estimated baryonic
acceleration and the scatter will also be larger (Yong Tian, private
communication). To investigate the effect of stellar mass, we redid our
analysis for the X-COP sample (as it had smaller errors), without including
the contribution of stellar mass to the total baryonic mass budget. We find
that the new slope is given by $m=1.29\pm 0.01$ and $b=9.33\pm 0.01$, which is
about $3\sigma$ (slope) and $2\sigma$ (intercept) different, compared to the
inclusion of stellar mass. Without stellar mass, the intrinsic scatter which
we obtain is 39% and the residual scatter is 0.17 dex. Therefore, both the
scatters are also increased if the stellar mass contribution is not included.
We also computed the residual scatter in our re-analysis of HIFLUGS data, by
excluding the stellar mass contribution. Upon doing that, we get a residual
scatter of 0.15 dex, based on the Gaussian standard deviation, which is larger
than the 0.13 dex residual scatter which we obtain after inclusion of stellar
mass. Therefore, we can conclude that excluding the stellar mass does increase
the scatter and could be one of the reasons for the large scatter in CDP20.

##  VI Conclusions

Recently, two independent groups (CDP20 and T20) carried out a systematic test
of a possible deterministic relation between the baryonic and total
acceleration for galaxy clusters, known in the astrophysical literature as
RAR. This was motivated by the very tight scatter found from a similar test
done for spiral galaxies using the SPARC sample McGaugh et al. (2016). T20 did
this test with the CLASH cluster sample using a combination of strong lensing,
weak lensing and X-ray measurements, whereas CDP20 used the hydrostatic masses
of the non-cool core clusters from the ROSAT dataset. T20 found a tight
empirical correlation between the two, with a slope consistent with the low-
acceleration limit found for galaxies. However, they find an acceleration
scale about ten times higher than for galaxies. However, CDP20 found a large
residual scatter about the RAR of about 0.18 dex, and a higher acceleration
scale about 4-10 times larger. Therefore, the CDP20 analysis shows that galaxy
clusters do not obey the RAR.

Given the somewhat contradictory nature of the two results, we carry out an
independent test of the same relation using two different cluster samples. The
first sample is a dataset of 12 clusters imaged using Chandra telescope by
Vikhlinin et al. The second sample which we use is the X-COP sample, in which
a subset of 12 Planck-SZ clusters with a high SNR were imaged using the XMM-
Newton telescope. For both these samples, we carry out a linear regression
relation between the natural log of baryonic and the same for total
acceleration. We also re-analyzed the HIFLUGCS cluster sample used in CDP20,
using the same regression relation used for the Chandra and X-COP sample, and
also by including the stellar mass contribution in the baryonic mass budget.
We characterize the tightness of the RAR by computing both the intrinsic
scatter and the residual scatter (in dex).

For the Chandra sample, we find that the slope and intercept disagree with the
estimates in T20 at about 2.2$\sigma$. However, they agree with the sample
analyzed in CDP20 to within 1$\sigma$. The residual scatter which we get is
about 0.14 dex. The X-COP sample (which has much smaller mass errors compared
to the Chandra sample) disagrees with the estimates for both the slope and
intercept, when compared to the T20 and CDP20 estimates at about $6.7\sigma$
and $4\sigma$, respectively. The residual scatter which we get for this sample
is about the same as T20 (0.11 dex). The intrinsic scatter which we obtain for
the X-COP sample is about 20%, which is about the same as that obtained in T20

Therefore, we conclude that the X-COP sample obeys the RAR in agreement with
the conclusions of T20. The best-fit errors for the Chandra sample are
somewhat large, and hence a decisive statement cannot be made for this sample.

However, given the differences in the best fits for the slopes and intercepts
between the different cluster samples, we find that this relation is not
universal for clusters.

We also find that exclusion of stellar mass changes the intercept and slope by
about $2-3\sigma$, and also increases both the residual and intrinsic scatter.
This may explain the high residual scatter seen in CDP20.

If we assume that the total acceleration scales as the square root of the
baryonic acceleration, then the characteristic acceleration scale we obtain is
about the same order of magnitude as that obtained in T20 and CDP20, which is
about an order of magnitude larger than what found using the SPARC sample.

Looking to the future, a more definitive test should soon be possible with the
recent launch of the eROSITA X-ray satellite Predehl et al. (2020), which
should provide X-ray measurements for about 100,000 clusters, along with
lensing and other mass estimates from follow-up programs, which would allow
high precision measurements of both the baryonic as well as the total mass.

###### Acknowledgements.

We are grateful to Stefano Ettori and Dominik Eckert for providing the data
and answering all our questions regarding the X-COP sample; Alexey Vikhlinin
for the same regarding the Chandra sample; Doug Edmonds, Yong Tian, Man-Ho
Chan, and Antonio Del Popolo for helpful correspondence. The work of Sajal
Gupta was supported by a DST-INSPIRE fellowship. We also thank the anonymous
referee for several constructive comments on our manuscript.

## References

  * Chan and Del Popolo (2020) M. H. Chan and A. Del Popolo, Mon. Not. R. Astron. Soc. p. 218 (2020), eprint 2001.06141.
  * Tian et al. (2020a) Y. Tian, K. Umetsu, C.-M. Ko, M. Donahue, and I. N. Chiu, Astrophys. J. 896, 70 (2020a), eprint 2001.08340.
  * McGaugh et al. (2016) S. S. McGaugh, F. Lelli, and J. M. Schombert, Phys. Rev. Lett. 117, 201101 (2016), eprint 1609.05917.
  * Lelli et al. (2016) F. Lelli, S. S. McGaugh, and J. M. Schombert, Astron. J. 152, 157 (2016), eprint 1606.09251.
  * Li et al. (2018) P. Li, F. Lelli, S. McGaugh, and J. Schombert, Astron. & Astrophys. 615, A3 (2018), eprint 1803.00022.
  * Chae et al. (2020) K.-H. Chae, M. Bernardi, H. Dominguez Sanchez, and R. K. Sheth, arXiv e-prints arXiv:2010.10779 (2020), eprint 2010.10779.
  * McGaugh (2004) S. S. McGaugh, Astrophys. J. 609, 652 (2004), eprint astro-ph/0403610.
  * Famaey and McGaugh (2012) B. Famaey and S. S. McGaugh, Living Reviews in Relativity 15, 10 (2012), eprint 1112.3960.
  * Milgrom (2016) M. Milgrom, arXiv e-prints arXiv:1609.06642 (2016), eprint 1609.06642.
  * Lelli et al. (2017a) F. Lelli, S. S. McGaugh, J. M. Schombert, and M. S. Pawlowski, Astrophys. J. 836, 152 (2017a), eprint 1610.08981.
  * Wheeler et al. (2019) C. Wheeler, P. F. Hopkins, and O. Doré, Astrophys. J. 882, 46 (2019), eprint 1803.01849.
  * Rodrigues et al. (2018a) D. C. Rodrigues, V. Marra, A. del Popolo, and Z. Davari, Nature Astronomy 2, 668 (2018a), eprint 1806.06803.
  * Zobnina and Zasov (2020) D. I. Zobnina and A. V. Zasov, Astronomy Reports 64, 295 (2020), eprint 2003.08845.
  * Marra et al. (2020) V. Marra, D. C. Rodrigues, and Á. O. F. de Almeida, Mon. Not. R. Astron. Soc. 494, 2875 (2020), eprint 2002.03946.
  * Zhou et al. (2020) Y. Zhou, A. Del Popolo, and Z. Chang, Physics of the Dark Universe 28, 100468 (2020), eprint 2008.04065.
  * McGaugh et al. (2018) S. S. McGaugh, P. Li, F. Lelli, and J. M. Schombert, Nature Astronomy 2, 924 (2018).
  * Kroupa et al. (2018) P. Kroupa, I. Banik, H. Haghi, A. H. Zonoozi, J. Dabringhausen, B. Javanmardi, O. Müller, X. Wu, and H. Zhao, Nature Astronomy 2, 925 (2018), eprint 1811.11754.
  * Rodrigues et al. (2018b) D. C. Rodrigues, V. Marra, A. Del Popolo, and Z. Davari, Nature Astronomy 2, 927 (2018b), eprint 1811.05882.
  * Rodrigues et al. (2020) D. C. Rodrigues, V. Marra, A. Del Popolo, and Z. Davari, Nature Astronomy 4, 134 (2020), eprint 2002.01970.
  * Desmond (2017) H. Desmond, Mon. Not. R. Astron. Soc. 464, 4160 (2017), eprint 1607.01800.
  * Navarro et al. (2017) J. F. Navarro, A. Benítez-Llambay, A. Fattahi, C. S. Frenk, A. D. Ludlow, K. A. Oman, M. Schaller, and T. Theuns, Mon. Not. R. Astron. Soc. 471, 1841 (2017), eprint 1612.06329.
  * Dai and Lu (2017) D.-C. Dai and C. Lu, Phys. Rev. D 96, 124016 (2017), eprint 1712.01654.
  * Keller and Wadsley (2017) B. W. Keller and J. W. Wadsley, Astrophys. J. Lett. 835, L17 (2017), eprint 1610.06183.
  * Garaldi et al. (2018) E. Garaldi, E. Romano-Díaz, C. Porciani, and M. S. Pawlowski, Phys. Rev. Lett. 120, 261301 (2018), eprint 1712.04448.
  * Ludlow et al. (2017) A. D. Ludlow, A. Benítez-Llambay, M. Schaller, T. Theuns, C. S. Frenk, R. Bower, J. Schaye, R. A. Crain, J. F. Navarro, A. Fattahi, et al., Phys. Rev. Lett. 118, 161103 (2017), eprint 1610.07663.
  * Stone and Courteau (2019) C. Stone and S. Courteau, Astrophys. J. 882, 6 (2019), eprint 1908.06105.
  * Ren et al. (2019) T. Ren, A. Kwa, M. Kaplinghat, and H.-B. Yu, Physical Review X 9, 031020 (2019), eprint 1808.05695.
  * Berezhiani et al. (2018) L. Berezhiani, B. Famaey, and J. Khoury, JCAP 2018, 021 (2018), eprint 1711.05748.
  * Green and Moffat (2019) M. A. Green and J. W. Moffat, Physics of the Dark Universe 25, 100323 (2019), eprint 1905.09476.
  * O’Brien et al. (2018) J. G. O’Brien, T. L. Chiarelli, and P. D. Mannheim, Physics Letters B 782, 433 (2018), eprint 1704.03921.
  * Lelli et al. (2017b) F. Lelli, S. S. McGaugh, and J. M. Schombert, Mon. Not. R. Astron. Soc. 468, L68 (2017b), eprint 1702.04355.
  * Voit (2005) G. M. Voit, Reviews of Modern Physics 77, 207 (2005), eprint astro-ph/0410173.
  * Vikhlinin et al. (2014) A. A. Vikhlinin, A. V. Kravtsov, M. L. Markevich, R. A. Sunyaev, and E. M. Churazov, Physics Uspekhi 57, 317-341 (2014).
  * Allen et al. (2011) S. W. Allen, A. E. Evrard, and A. B. Mantz, Ann. Rev. Astron. Astrophys. 49, 409 (2011), eprint 1103.4829.
  * Desai (2018) S. Desai, Physics Letters B 778, 325 (2018), eprint 1708.06502.
  * The and White (1988) L. S. The and S. D. M. White, Astron. J. 95, 1642 (1988).
  * Gerbal et al. (1992) D. Gerbal, F. Durret, M. Lachieze-Rey, and G. Lima-Neto, Astron. & Astrophys. 262, 395 (1992).
  * Aguirre et al. (2001) A. Aguirre, J. Schaye, and E. Quataert, Astrophys. J. 561, 550 (2001), eprint astro-ph/0105184.
  * Sanders (2003) R. H. Sanders, Mon. Not. R. Astron. Soc. 342, 901 (2003), eprint astro-ph/0212293.
  * Pointecouteau and Silk (2005) E. Pointecouteau and J. Silk, Mon. Not. R. Astron. Soc. 364, 654 (2005), eprint astro-ph/0505017.
  * Angus et al. (2008) G. W. Angus, B. Famaey, and D. A. Buote, Mon. Not. R. Astron. Soc. 387, 1470 (2008), eprint 0709.0108.
  * Natarajan and Zhao (2008) P. Natarajan and H. Zhao, Mon. Not. R. Astron. Soc. 389, 250 (2008), eprint 0806.3080.
  * Dodelson (2011) S. Dodelson, International Journal of Modern Physics D 20, 2749 (2011), eprint 1112.1320.
  * Boran et al. (2018) S. Boran, S. Desai, E. O. Kahya, and R. P. Woodard, Phys. Rev. D 97, 041501 (2018), eprint 1710.06168.
  * Famaey et al. (2018) B. Famaey, J. Khoury, and R. Penco, JCAP 2018, 038 (2018), eprint 1712.01316.
  * Chen et al. (2007) Y. Chen, T. H. Reiprich, H. Böhringer, Y. Ikebe, and Y. Y. Zhang, Astron. & Astrophys. 466, 805 (2007), eprint astro-ph/0702482.
  * Cavaliere and Fusco-Femiano (1978) A. Cavaliere and R. Fusco-Femiano, Astron. & Astrophys. 70, 677 (1978).
  * Postman et al. (2012) M. Postman, D. Coe, N. Benítez, L. Bradley, T. Broadhurst, M. Donahue, H. Ford, O. Graur, G. Graves, S. Jouvel, et al., Astrophys. J. Suppl. Ser. 199, 25 (2012), eprint 1106.3328.
  * Umetsu et al. (2016) K. Umetsu, A. Zitrin, D. Gruen, J. Merten, M. Donahue, and M. Postman, Astrophys. J. 821, 116 (2016), eprint 1507.04385.
  * Donahue et al. (2014) M. Donahue, G. M. Voit, A. Mahdavi, K. Umetsu, S. Ettori, J. Merten, M. Postman, A. Hoffer, A. Baldi, D. Coe, et al., Astrophys. J. 794, 136 (2014), eprint 1405.7876.
  * Navarro et al. (1997) J. F. Navarro, C. S. Frenk, and S. D. M. White, Astrophys. J. 490, 493 (1997), eprint astro-ph/9611107.
  * Chiu et al. (2018) I. Chiu, J. J. Mohr, M. McDonald, S. Bocquet, S. Desai, M. Klein, H. Israel, M. L. N. Ashby, A. Stanford, B. A. Benson, et al., Mon. Not. R. Astron. Soc. 478, 3072 (2018), eprint 1711.00917.
  * Hernquist (1990) L. Hernquist, Astrophys. J. 356, 359 (1990).
  * Cooke et al. (2016) K. C. Cooke, C. P. O’Dea, S. A. Baum, G. R. Tremblay, I. G. Cox, and M. Gladders, Astrophys. J. 833, 224 (2016), eprint 1610.05310.
  * Olamaie et al. (2012) M. Olamaie, M. P. Hobson, and K. J. B. Grainge, Mon. Not. R. Astron. Soc. 423, 1534 (2012), eprint 1109.2834.
  * Tian et al. (2020b) Y. Tian, P.-C. Yu, P. Li, S. S. McGaugh, and C.-M. Ko, arXiv e-prints arXiv:2010.00992 (2020b), eprint 2010.00992.
  * Vikhlinin et al. (2005) A. Vikhlinin, M. Markevitch, S. S. Murray, C. Jones, W. Forman, and L. Van Speybroeck, Astrophys. J. 628, 655 (2005), eprint astro-ph/0412306.
  * Vikhlinin et al. (2006) A. Vikhlinin, A. Kravtsov, W. Forman, C. Jones, M. Markevitch, S. S. Murray, and L. Van Speybroeck, Astrophys. J. 640, 691 (2006), eprint astro-ph/0507092.
  * Rahvar and Mashhoon (2014) S. Rahvar and B. Mashhoon, Phys. Rev. D 89, 104011 (2014), eprint 1401.4819.
  * Hodson and Zhao (2017) A. O. Hodson and H. Zhao, Astron. & Astrophys. 598, A127 (2017), eprint 1701.03369.
  * Hodson et al. (2017) A. O. Hodson, H. Zhao, J. Khoury, and B. Famaey, Astron. & Astrophys. 607, A108 (2017), eprint 1611.05876.
  * Edmonds et al. (2018) D. Edmonds, D. Farrah, D. Minic, Y. J. Ng, and T. Takeuchi, International Journal of Modern Physics D 27, 1830001-296 (2018), eprint 1709.04388.
  * Edmonds et al. (2017) D. Edmonds, D. Farrah, C. M. Ho, D. Minic, Y. J. Ng, and T. Takeuchi, International Journal of Modern Physics A 32, 1750108 (2017), eprint 1601.00662.
  * Bernal et al. (2017) T. Bernal, V. H. Robles, and T. Matos, Mon. Not. R. Astron. Soc. 468, 3135 (2017), eprint 1609.08644.
  * Gupta and Desai (2019) S. Gupta and S. Desai, Classical and Quantum Gravity 36, 105001 (2019), eprint 1811.09378.
  * Gupta and Desai (2020) S. Gupta and S. Desai, Physics of the Dark Universe 28, 100499 (2020), eprint 1909.07408.
  * Gopika and Desai (2020) K. Gopika and S. Desai, Physics of the Dark Universe 30, 100707 (2020), eprint 2006.12320.
  * Tian and Ko (2019) Y. Tian and C.-M. Ko, Mon. Not. R. Astron. Soc. 488, L41 (2019), eprint 1907.04501.
  * Biffi et al. (2016) V. Biffi, S. Borgani, G. Murante, E. Rasia, S. Planelles, G. L. Granato, C. Ragone-Figueroa, A. M. Beck, M. Gaspari, and K. Dolag, Astrophys. J. 827, 112 (2016), eprint 1606.02293.
  * Lin et al. (2012) Y.-T. Lin, S. A. Stanford, P. R. M. Eisenhardt, A. Vikhlinin, B. J. Maughan, and A. Kravtsov, Astrophys. J. Lett. 745, L3 (2012), eprint 1112.1705.
  * Vikhlinin et al. (2009) A. Vikhlinin, R. A. Burenin, H. Ebeling, W. R. Forman, A. Hornstrup, C. Jones, A. V. Kravtsov, S. S. Murray, D. Nagai, H. Quintana, et al., Astrophys. J. 692, 1033 (2009), eprint 0805.2207.
  * Foreman-Mackey et al. (2013) D. Foreman-Mackey, D. W. Hogg, D. Lang, and J. Goodman, Publ. Astron. Soc. Pac. 125, 306 (2013), eprint 1202.3665.
  * Astropy Collaboration et al. (2018) Astropy Collaboration, A. M. Price-Whelan, B. M. Sipőcz, H. M. Günther, P. L. Lim, S. M. Crawford, S. Conseil, D. L. Shupe, M. W. Craig, N. Dencheva, et al., Astron. J. 156, 123 (2018), eprint 1801.02634.
  * Planck Collaboration et al. (2014) Planck Collaboration, P. A. R. Ade, N. Aghanim, C. Armitage-Caplan, M. Arnaud, M. Ashdown, F. Atrio-Barand ela, J. Aumont, H. Aussel, C. Baccigalupi, et al., Astron. & Astrophys. 571, A29 (2014), eprint 1303.5089.
  * Ghirardini et al. (2019) V. Ghirardini, D. Eckert, S. Ettori, E. Pointecouteau, S. Molendi, M. Gaspari, M. Rossetti, S. De Grandi, M. Roncarelli, H. Bourdin, et al., Astron. & Astrophys. 621, A41 (2019), eprint 1805.00042.
  * Eckert et al. (2017) D. Eckert, S. Ettori, E. Pointecouteau, S. Molendi, S. Paltani, and C. Tchernin, Astronomische Nachrichten 338, 293 (2017), eprint 1611.05051.
  * Ettori et al. (2019) S. Ettori, V. Ghirardini, D. Eckert, E. Pointecouteau, F. Gastaldello, M. Sereno, M. Gaspari, S. Ghizzardi, M. Roncarelli, and M. Rossetti, Astron. & Astrophys. 621, A39 (2019), eprint 1805.00035.
  * Ghirardini et al. (2018) V. Ghirardini, S. Ettori, D. Eckert, S. Molendi, F. Gastaldello, E. Pointecouteau, G. Hurier, and H. Bourdin, Astron. & Astrophys. 614, A7 (2018), eprint 1708.02954.
  * Ettori et al. (2020) S. Ettori, V. Ghirardini, and D. Eckert, Astronomische Nachrichten 341, 210 (2020), eprint 1911.05039.
  * Ghizzardi et al. (2020a) S. Ghizzardi, S. Molendi, R. van der Burg, S. De Grandi, I. Bartalucci, F. Gastaldello, M. Rossetti, V. Biffi, S. Borgani, D. Eckert, et al., arXiv e-prints arXiv:2007.01084 (2020a), eprint 2007.01084.
  * Sunyaev and Zeldovich (1970) R. A. Sunyaev and Y. B. Zeldovich, Astrophysics & Space Sciences 7, 3 (1970).
  * Planck Collaboration et al. (2016) Planck Collaboration, P. A. R. Ade, N. Aghanim, M. Arnaud, M. Ashdown, J. Aumont, C. Baccigalupi, A. J. Banday, R. B. Barreiro, R. Barrena, et al., Astron. & Astrophys. 594, A27 (2016), eprint 1502.01598.
  * Eckert et al. (2019) D. Eckert, V. Ghirardini, S. Ettori, E. Rasia, V. Biffi, E. Pointecouteau, M. Rossetti, S. Molendi, F. Vazza, F. Gastaldello, et al., Astron. & Astrophys. 621, A40 (2019), eprint 1805.00034.
  * Ghizzardi et al. (2020b) S. Ghizzardi, S. Molendi, R. van der Burg, S. De Grandi, I. Bartalucci, F. Gastaldello, M. Rossetti, V. Biffi, S. Borgani, D. Eckert, et al., arXiv e-prints arXiv:2007.01084 (2020b), eprint 2007.01084.
  * Predehl et al. (2020) P. Predehl, R. Andritschke, V. Arefiev, V. Babyshkin, O. Batanov, W. Becker, H. Böhringer, A. Bogomolov, T. Boller, K. Borm, et al., arXiv e-prints arXiv:2010.03477 (2020), eprint 2010.03477.