# All-in-one: Multi-task Learning for Rumour Verification

Elena Kochkina1,2, Maria Liakata1,2, Arkaitz Zubiaga1
1 University of Warwick, Coventry, United Kingdom
2 Alan Turing Institute, London, United Kingdom
{E.Kochkina, M.Liakata<EMAIL_ADDRESS>
###### Abstract

Automatic resolution of rumours is a challenging task that can be broken down
into smaller components that make up a pipeline, including rumour detection,
rumour tracking and stance classification, leading to the final outcome of
determining the veracity of a rumour. In previous work, these steps in the
process of rumour verification have been developed as separate components
where the output of one feeds into the next. We propose a multi-task learning
approach that allows joint training of the main and auxiliary tasks, improving
the performance of rumour verification. We examine the connection between the
dataset properties and the outcomes of the multi-task learning models used.

Title and Abstract in Russian

Три в одном: Многозадачное Обучение для Верификации Слухов

  * Автоматическая верификация слухов в Интернете - это сложная задача, которая может быть разбита на подзадачи, включающие в себя обнаружение слухов, отслеживание слухов, классификацию отношения пользователей к правдивости этих слухов и, в итоге, разрешение основной задачи - определение достоверности слуха. В предыдущих исследованиях модели для решения каждой из подзадач разрабатывались как отдельные компоненты, где вывод модели для задачи предыдущего этапа являлся вводдом для модели задачи следующего этапа. Мы предлагаем использование многозадачного обучения для совместного обучения основной и вспомогательных задач в одной модели, что позволяет улучшить качество модели верификации слухов. Также мы рассматриваем связь между свойствами набора данных и результатами использования моделей с многозадачным обучением.

##  1 Introduction

00footnotetext:  This work is licenced under a Creative Commons Attribution
4.0 International Licence. Licence details:
http://creativecommons.org/licenses/by/4.0/

Social media have gained popularity as platforms that enable users to follow
events and breaking news [Hu et al., 2012]. However, not all information that
spreads on social media during such events is accurate. The serious harm that
inaccurate information can cause to society in critical situations
[Lewandowsky et al., 2012] has led to an increased interest within the
scientific community to develop tools to verify information from social media
[Shu et al., 2017]. Likewise, social media platforms themselves, such as
Facebook,111https://www.recode.net/2017/2/16/14632726/mark-zuckerberg-
facebook-manifesto-fake-news-terrorism are investing significant effort in
mitigating the problems caused by misinformation. One of the features that
characterises social media is the rapid emergence and spread of new
information. This leads to the circulation of rumours, claims that are
unverified at the time of posting, and for which eventually evidence proving
their true or false nature may come out [DiFonzo and Bordia, 2007]. Therefore,
a rumour resolution system needs to go through a set of steps, from detecting
that a new circulating claim is a rumour, to the ultimate step of determining
its veracity value.

The rumour resolution process has been defined as a pipeline involving four
sub-tasks [Zubiaga et al., 2018a] (see Figure 1): (1) rumour detection,
determining whether a claim is worth verifying rather than the expression of
an opinion; (2) rumour tracking, collecting sources and opinions on a rumour
as it unfolds; (3) stance classification, determining the attitude of the
sources or users towards the truthfulness of the rumour, and (4) rumour
verification, as the ultimate step where the veracity value of the rumour is
predicted. These steps can be performed at different times in the life-cycle
of a rumour, making this a time-sensitive process. Ideally, rumours can be
resolved as either true or false. However, they can also remain unverified
when there is no sufficient evidence to determine their veracity [Caplow,
1947].

A recent body of work studies each of these four tasks separately [Lukasik et
al., 2016, Liu et al., 2016, Enayet and El-Beltagy, 2017]. However, the way
these subtasks interact and their integration into a complete rumour
resolution system is yet to be explored. In this work we assume that veracity
classification is the crucial component of a rumour resolution system, as the
final output that determines, given information collected about a rumour at a
certain point in time, if the rumour is true, false, or remains unverified. We
express the rumour resolution process as a multi-task problem that needs to
address a number of challenges, where the veracity classification task is the
main task and the rest of the components are auxiliary tasks that can be
leveraged to boost the performance of the veracity classifier.

We propose to achieve this setting by using a multi-task learning approach.
Multi-task learning [Caruana, 1998] refers to the joint training of multiple
tasks, which has gained popularity recently for a range of tasks in Machine
Learning and Natural Language Processing [Collobert and Weston, 2008] and has
been applied in a number of different tasks and machine learning
architectures. Its effectiveness is mainly attributed to learning shared
representations of closely related tasks, such that two complementary tasks
can give each other ‘hints’.

We propose joint learning of the tasks in the verification pipeline, that will
allow us to leverage relations between the tasks. We assess the effectiveness
of using a multi-task learning approach for the rumour resolution process in
four different scenarios: (1) performing single task learning to perform
veracity classification, (2) performing multi-task learning that combines
stance and veracity classification with the aim of boosting performance of the
latter, (3) performing multi-task learning that combines rumour detection and
verification, comparing the improvement brought by each of the two auxiliary
tasks, and (4) performing multi-task learning that combines rumour detection,
stance classification and veracity prediction to improve the performance for
veracity. We employ a deep learning architecture that views each subtask as
having a sequential nature. We compare our results with a state-of-the-art
veracity classification approach introduced by Enayet and El-Beltagy [Enayet
and El-Beltagy, 2017].

While a lot of work reports positive outcomes of the application of multi-task
learning to various NLP tasks [Collobert and Weston, 2008, Aguilar et al.,
2017, Lan et al., 2017], there are also studies showing that this is not
always the case [Alonso and Plank, 2017, Bingel and Søgaard, 2017]. Alonso and
Plank [Alonso and Plank, 2017] were the first to demonstrate that multi-task
learning brings benefits only for some combinations of main and auxiliary
tasks. They also investigate the relationship between the multi-task learning
outcome and the properties of the dataset. We perform a similar analysis to
examine the link between the properties of our rumour datasets and the results
of the multi-task learning approach used.

Our results show that a multi-task learning scenario that leverages all three
subtasks, where veracity classification is the main task and stance
classification and rumour detection are auxiliary tasks, leads to substantial
improvements over a standard, single task veracity classification system, as
well as a majority baseline and a state-of-the-art system. The combination of
all three subtasks also outperforms the multi-task learning scenario where
only two of the subtasks are combined.

##  2 Related work

###  2.1 Rumour classification system

Figure 1: Rumour resolution pipeline.

Rumour classification is a complex task that can be represented as a pipeline
of sub-tasks. While there are different possibilities for structuring this
pipeline, here we adopt the architecture of a rumour classification system
proposed in Zubiaga et al. [Zubiaga et al., 2018a] (see Figure 1). Thus, a
rumour classification system is expressed as a sequence of subtasks, namely
rumour detection, rumour tracking, rumour stance classification leading to
rumour verification. In this work we consider three of these tasks: rumour
detection, stance and verification, focusing on the evaluation of the latter.
Rumour tracking does not involve classification but rather consists of
collecting tweets following up on a rumour in the form of replies.

A recent survey [Zubiaga et al., 2018a] defines rumour detection as the task
which distinguishes rumours (unverified pieces of information) from non-
rumours (all other circulating information). Here we adhere to this
definition. Therefore, once rumour detection is performed, the information
classified as a rumour is then fed into the the stance and veracity
classification components of a system to ultimately determine the veracity of
the rumour.

While we are not aware of previous work tackling the entire rumour
classification pipeline, there has been a body of work tackling each of the
tasks individually.

#### Stance classification.

An increasing body of work has focused on the stance classification component,
i.e. determining whether different posts associated with a rumour support it,
deny it, query it, or just comment on it. The stance expressed by users
towards a particular rumour can be indicative of the veracity of a rumour,
furthermore, it has been shown that rumours attracting higher levels of
skepticism in the form of denials and questioning responses are more likely to
be proven false later [Mendoza et al., 2010, Procter et al., 2013, Derczynski
et al., 2014]. Previous work on stance classification [Lukasik et al., 2016,
Kochkina et al., 2017, Zubiaga et al., 2017, Zubiaga et al., 2018b] has
explored the use of sequential classifiers, treating the task as one that
evolves over time; it has been shown that these sequential classifiers
substantially outperform standard, non-sequential classifiers.

#### Rumour detection.

Work on rumour detection is more scarce. One of the first approaches was
introduced by Zhao et al. [Zhao et al., 2015], who built a rule-based approach
to identify skepticism (e.g. is this true?) and therefore determine that the
associated information is a rumour. The limitations of this approach consist
in having to wait for responses to arrive, as well as lack of generalisability
due to manually-defined rules. Zubiaga et al [Zubiaga et al., 2017] proposed a
sequential approach to leverage context from earlier posts during an event.
The sequential approach achieved significant improvements, especially in terms
of recall, where the rule-based approach proved limited.

#### Rumour verification.

Common approaches to rumour verification involve first the collection of
corpora of resolved rumours from rumour debunking websites such as snopes.com,
emergent.com, politifact.com. Wang [Wang, 2017] created a dataset based on
claims from politifact.com that are annotated for degrees of truthfulness. To
resolve these they propose a hybrid convolutional neural network that
integrates metadata with text. Twitter is a popular platform to study rumours,
while often seeds and annotations for rumours are still taken from rumour
debunking websites. Giasemidis et al. [Giasemidis et al., 2016] collected a
dataset of 72 rumours from Twitter. This work measured trustworthiness of a
claim at varying time windows. Boididou et al. [Boididou et al., 2014,
Boididou et al., 2017] have focused on tweets that are related to fake images,
bringing in the multimedia component into the verification process.

#### Sequence classification.

Works that consider tasks of rumour detection and verification (within their
own definition of the tasks, not necessarily aligning with the one used here)
also highlighted the importance of a sequential time-sensitive approach when
dealing with rumours [Ma et al., 2016, Kwon et al., 2017, Chen et al., 2017].
As previous work in the literature has highlighted the importance of
sequential classifiers in rumour detection we have decided to use an LSTM
based architecture for our experiments.

Figure 2: Example of a conversation with 3 branches from the RumourEval
dataset.

###  2.2 Multi-task learning

Multi-task learning refers to the joint learning of several related tasks with
a shared representation. In recent years, a multi-task learning approach has
been successfully applied in combination with neural networks for a variety of
NLP tasks [Collobert and Weston, 2008]. In this work we use the most common
approach to multi-task learning, namely hard parameter sharing, meaning that
different tasks are using the same hidden layer(s). The effectiveness of a
multi-task learning approach is attributed to: effectively increasing the size
of the training set by using additional datasets for related tasks and
regularisation, as the model has to learn a shared representation for multiple
tasks there is less risk of overfitting on one of them.

In multi-task learning auxiliary tasks can be used to direct the main task to
use/learn the features that it otherwise would have ignored or failed to
identify due to complex relations between the tasks and features. For example,
multi-task learning is particularly useful when potentially helpful features
are not used as such for the main task, but become labels in an auxiliary
task. This use case is very relevant to this work; stance classification could
be used as a feature in a system for veracity classification, as has indeed
been the case in previous studies, which have shown a relationship between the
two tasks. Zhao et al. [Zhao et al., 2015] and Enayet and El-Beltagy [Enayet
and El-Beltagy, 2017] have created successful models using this premise.
However these studies assume access to stance and veracity labels for the same
data, which does not apply in our case as we do not have stance labels for all
of the threads in the dataset (see section 3).

##  3 Data

We use two publicly available datasets of rumours: PHEME [Zubiaga et al.,
2016, Zubiaga et al., 2017] and RumourEval [Derczynski et al., 2017] that
contain different levels of annotation for the tasks of rumour detection,
rumour stance and veracity classification. Both datasets contain Twitter
conversation threads associated with different newsworthy events including the
Ferguson unrest, the shooting at Charlie Hebdo, the shooting in Ottawa, the
hostage situation in Sydney and the crash of a Germanwings plane. Figure 2
shows an example of a conversation discussing a rumour about Michael Essien
having contracted Ebola. The conversation consists of a source tweet conveying
a rumour and a tree of responses, expressing their opinion towards the claim
contained in the source tweet. The veracity label of the rumour is false,
while each of the responses could be tagged as either supporting, denying,
questioning or commenting on the rumour. Conversations can be decomposed into
branches, such that a branch is a linear sequence of tweets starting from a
leaf node of the conversation tree and going through its parent nodes up to
the source tweet. The conversation on the Figure 2 can be decomposed into
three branches.

###  3.1 RumourEval

RumourEval is a dataset that was released as part of the SemEval-2017 Task 8
competition [Derczynski et al., 2017]. It contains 325 Twitter threads
discussing rumours. All conversations in the RumourEval dataset are rumours,
therefore this dataset only covers the tasks of rumour stance and veracity
classification. It is split into training, testing and development sets. The
testing set contains a mix of rumours related to the same events as in the
training and development sets, with the addition of two rumours: about Marina
Joyce and the health condition of Hillary Clinton. Table 1 shows the number of
conversation threads, branches and tweets in each of the sets of the
RumourEval dataset, as well as the class distribution for both tasks. In the
stance classification task there is a strong class imbalance towards the
commenting class, while questioning and denying are minority classes, which
holds true for all subsets. For the verification task, the training set
contains more true instances than false or unverified, whereas the development
and testing sets are more balanced.

| Threads | Branches | Tweets | True | False | Unverified | S | D | Q | C
---|---|---|---|---|---|---|---|---|---|---
Development | 25 | 215 | 281 | 10 | 12 | 3 | 69 | 11 | 28 | 173
Testing | 28 | 772 | 1049 | 8 | 12 | 8 | 94 | 71 | 106 | 778
Training | 272 | 3030 | 4238 | 127 | 50 | 95 | 841 | 333 | 330 | 2734
Total | 325 | 4017 | 5568 | 145 | 74 | 106 | 1004 | 415 | 464 | 3685

Table 1: Number of threads, tweets and class distribution in RumourEval
dataset.

###  3.2 PHEME

We are making available the extended version of the PHEME dataset of rumours
and non-rumours related to nine events, each updated with veracity information
for each of the
rumours222https://figshare.com/articles/PHEME_dataset_for_Rumour_Detection_and_Veracity_Classification/6392078.

This dataset contains 3 levels of annotation. First, each thread is annotated
as either rumour or non-rumour; second, rumours are labeled as either true,
false or unverified. And third, a subset (threads used in RumourEval) is
annotated for stance classification at the tweet level through crowd-sourcing.

Rumours in this dataset were labeled as true, false and unverified by
professional journalists [Zubiaga et al., 2016]. While stance can be labeled
by non-expert workers as labels can be inferred directly from the text, the
rumour verification task is more challenging as it requires analysis of the
context, and further understanding of the rumours in order to determine if the
underlying story is true, false, or remains unverified. To assess the
difficulty of performing the verification task by a non-expert, an author of
this paper went through the rumours and annotated them for veracity. The
overlap between the non-expert annotator and the journalist was within the
range of $60-65\%$ on the rumour stories from five largest events.

Table 2 shows the size of each event in the PHEME dataset as well as the label
distribution for the tasks of rumour detection and verification. The
information about the stance classification task is in table 1. Events differ
in size drastically and have different class label proportions. Overall, the
PHEME dataset contains fewer rumours than non-rumours, while the majority
class for rumours is true.

Events | Threads | Tweets | Rumours | Non-rumours | True | False | Unverified
---|---|---|---|---|---|---|---
Charlie Hebdo | 2,079 | 38,268 | 458 | 1,621 | 193 | 116 | 149
Sydney siege | 1,221 | 23,996 | 522 | 699 | 382 | 86 | 54
Ferguson | 1,143 | 24,175 | 284 | 859 | 10 | 8 | 266
Ottawa shooting | 890 | 12,284 | 470 | 420 | 329 | 72 | 69
Germanwings-crash | 469 | 4,489 | 238 | 231 | 94 | 111 | 33
Putin missing | 238 | 835 | 126 | 112 | 0 | 9 | 117
Prince Toronto | 233 | 902 | 229 | 4 | 0 | 222 | 7
Gurlitt | 138 | 179 | 61 | 77 | 59 | 0 | 2
Ebola Essien | 14 | 226 | 14 | 0 | 0 | 14 | 0
Total | 6,425 | 105,354 | 2,402 | 4,023 | 1,067 | 638 | 697

Table 2: Number of threads, tweets and class distribution in the PHEME
dataset.

We perform two types of experiments using different subsets of the PHEME
dataset: (1) using the five largest events (see Table 2) and (2) using all
nine events. The five largest events create a more balanced dataset as those
are major crisis events during which true updates and false information on
different aspects of the event were shared and discussed, whereas the 4
smallest events only contain a single rumour story at the core of the event.
In both cases, cross-validation experiments are performed by relying on a
leave-one-event-out setting, i.e. using all the events except the one left for
testing in each case. This set up is more challenging than the one presented
in RumourEval, however it is more representative of the real world. We micro-
average the performance across events.

##  4 Models

###  4.1 Sequential approach

As the benefits of using a sequential approach were suggested by previous
studies of rumour stance classification and rumour detection tasks (discussed
in Section 2), we are following the branchLSTM approach described in Zubiaga
et al. [Zubiaga et al., 2018b]. We split the conversations into linear
branches and use them as training instances that become an input to a model
consisting of an LSTM layer followed by several dense ReLU layers and a
softmax layer that predicts class probabilities. As the stance classification
task is annotated at the tweet level, we use the output from each time step of
LSTM, whereas for the tasks of rumour detection and verification, annotated at
the thread level, we are only using outputs from the final time steps (this
idea is illustrated in Figure 3). To get per-thread predictions, we use
majority voting for each of the branches from the thread. The model is trained
using categorical cross entropy loss.

###  4.2 Multi-task learning approach

We leverage the relationship between the tasks from the rumour classification
pipeline in a joint multi-task learning setup. Figure 3 illustrates our
approach. At the base of it is a sequential approach, as discussed above,
represented by a shared LSTM layer (hard parameter sharing), which is followed
by a number of task-specific layers. The possible task combinations are shown
as dotted lines on Figure 3 that can be present or absent depending on the
combination. We perform experiments in three set ups: joint training of (1)
stance or (2) rumour detection together with veracity classification, and (3)
learning all three tasks together. The cost function in the multi-task models
is a sum of losses from each of the tasks. Datasets for each of the three
tasks are not equal in size, therefore when the training instance is lacking a
label for one of the tasks, its prediction does not add anything to the loss
function, as if it had been predicted correctly.

Figure 3: Multi-task learning models. Dotted lines represent that same set up
is used for learning one, two or three tasks. Input format is a branch of
tweets. Stance classification implies predictions per-tweet, whereas detection
and verification tasks only require per-branch output.

###  4.3 Baselines

We compare proposed sequential and multi-task learning approaches with several
baselines. First of all, majority vote, a strong baseline which results in
high accuracy due to the class imbalance in the veracity classification task.
Another strong benchmark is NileTMRG [Enayet and El-Beltagy, 2017], the best
veracity classification system from SemEval-2017 Task 8. The NileTMRG model is
based on a linear SVM that uses a bag-of-words representation of the tweet
concatenated with selected features: presence of URL, presence of hashtag and
proportion of supporting, denying and querying tweets in the thread. We have
made our own implementation NileTMRG* of the NileTMRG model based on their
SemEval paper. This model requires stance labels for each of the tweets in the
dataset, however these are not available for the PHEME dataset. To obtain the
proportion of different stance labels among responses in the thread, we use
our own implementation of the stance classification model from Kochkina et al.
[Kochkina et al., 2017]. On the RumourEval dataset our implementation
NileTMRG* achieves better results than the original NileTMRG score [Enayet and
El-Beltagy, 2017] (accuracy 0.570 vs reported 0.536) as it utilises a better
model for stance classification. The NileTMRG model is a baseline showing the
scenario where pipeline tasks are performed sequentially and the outcome of
the previous step (stance classification) is an input to the next one
(veracity classification). The NileTMRG model has set a precedent in
demonstrating that patterns of support are useful indicators of rumour
veracity.

###  4.4 Features

We perform the following pre-processing of the tweets in the datasets: remove
nonalphabetic characters, convert all words to lower case and tokenise
texts333For implementation of all pre-processing routines we use Python 2.7
with the NLTK package. Once tweet texts are pre-processed, we extract word2vec
word embeddings pre-trained on the Google News dataset (300d)444Embeddings
were retrieved from https://code.google.com/archive/p/word2vec/. Processing
was performed using gensim package [Řehůřek and Sojka, 2010].  [Mikolov et
al., 2013] for each word in a tweet and take the average, thus obtaining a
tweet representation. This representation was shown to work well for tweets
due to their short length [Kochkina et al., 2017].

##  5 Experiment setup

###  5.1 Hyperparameters

We determined the optimal set of hyperparameters by testing the performance of
our models on the development set for different parameter combinations. We
used the Tree of Parzen Estimators (TPE) algorithm555We used the
implementation of the TPE algorithm in the hyperopt package to search the
parameter space and minimise the loss function expressed as $(1-macroF)$ for
single task, and $(1-macroF_{a})(1-macroF_{b})$ or
$(1-macroF_{a})(1-macroF_{b})(1-macroF_{c})$ for multi-task learning of two
and three tasks respectively. This loss function gives equal weight to all
tasks. The parameter space is defined as follows: the number of dense ReLU
layers varies from one to four; the number of LSTM layers is {1, 2}; the mini-
batch size is 32; the number of units in the ReLU layer is {300, 400, 500,
600}, and in the LSTM layer is {100, 200, 300}; the strength of the L2
regularisation is {$10^{-4}$, $10^{-3}$} and the number of epochs is 50. We
performed 30 trials of different parameter combinations optimising for
accuracy on the development set in order to choose the best combination. We
also use 50% dropout before the output layer. Zero-padding and masks that
account for the varying lengths of the input branches were used in all our
models. Models were
implemented666https://github.com/kochkinaelena/Multitask4Veracity using Python
3 and the Keras package.

###  5.2 Evaluation

The RumourEval dataset was provided with a training/development/testing split.
We tune parameters on the development set and then retrain the model on the
combined training and development sets before evaluation on the testing set.
On the PHEME dataset we perform leave-one-event-out (LOEO) cross-validation,
which makes this task set up harder than for RumourEval but closer to the
realistic scenario where we want to verify unseen rumours. When choosing
parameters, the Charlie Hebdo event was used as the development set as it has
balanced labels. We evaluate models using accuracy and macro-averaged F-score
as the tasks in the PHEME dataset suffer from a class imbalance.

##  6 Results and Discussion

The main results of our experiments are presented in Table 3. It shows the
results of the multi-task learning models MTL2 with two tasks: Veracity+Stance
and Veracity+Detection; MTL3 with three tasks Stance+Veracity+Detection,
Majority, NileTMRG* and single task branchLSTM baselines on the main task,
veracity classification.

As the datasets contain a significant class imbalance, the majority baseline
achieves fairly high accuracy scores. However due to the nature of this task,
it is more important for a model to recognize all of the classes, especially
false rumours, therefore the macro-averaged F-score is more important for
performance evaluation. All models demonstrate improvement over the majority
baseline in terms of macro F-score. We observe improvements of multi-task
approaches over single task learning in both accuracy and macro F-score, and
adding the third task brings further improvement.

MTL2 Veracity+Detection and MTL3 experiments were performed only on the PHEME
dataset as the RumourEval dataset consists only of rumours (no rumour
detection). On the RumourEval dataset the single task model branchLSTM
outperforms the majority baseline, although it does not perform as well as
NileTMRG*, while MTL2 shows improvement over both NileTMRG* and branchLSTM.
Experiments on the PHEME dataset also show a pattern of increasing scores:
MTL2 outperforms single task models and MTL3 outperforms MTL2.

RumourEval
---
| Majority (True) | NileTMRG* | branchLSTM |  MTL2 Veracity+Stance |  MTL2 Veracity+Detection | MTL3
Macro F | 0.148 | 0.539 | 0.491 | 0.558 | - | -
Accuracy | 0.286 | 0.570 | 0.500 | 0.571 | - | -
PHEME 5 events
| Majority (True) | NileTMRG* | branchLSTM |  MTL2 Veracity+Stance |  MTL2 Veracity+Detection | MTL3
Macro F | 0.226 | 0.339 | 0.336 | 0.376 | 0.373 | 0.396
Accuracy | 0.511 | 0.438 | 0.454 | 0.441 | 0.410 | 0.492
PHEME 9 events
| Majority (True) | NileTMRG* | branchLSTM |  MTL2 Veracity+Stance |  MTL2 Veracity+Detection | MTL3
Macro F | 0.205 | 0.297 | 0.259 | 0.318 | 0.345 | 0.405
Accuracy | 0.444 | 0.360 | 0.314 | 0.357 | 0.397 | 0.405

Table 3: Comparison of performance of sequential single task approach, multi-
task learning approaches with two and three tasks with majority and NileTMRG
baselines on veracity classification task.

Comparing the performance of the MTL2 model using stance as an auxiliary task
with the model using rumour detection as an auxiliary task on the PHEME 5
events dataset we observe that both models bring an improvement over the
single task branchLSTM baseline.

When using 9 events of the PHEME dataset we observe worse performance than on
5 events. Even though we are adding more training data, the 4 additional
events are qualitatively different to the 5 large news-breaking events. Each
of these 5 large events contained rumours labeled with all classes as well as
non-rumours, whereas the 4 additional events are small and the event itself is
a false or unverified rumour. This highlights the difficulty of the rumour
verification task in a leave-one-event-out setup and the importance of high
quality data. NileTMRG* is a very strong baseline, and while the single task
branchLSTM model is competitive when we are using 5 largest events, NileTMRG*
is only outperformed by the multi-task learning approach.

###  6.1 Per-event and per-class results analysis

Here we analyse the performance of proposed models on each of the 5 largest
events. Figure 4 illustrates the comparison of macro-averaged F-scores of the
proposed models for each of the events. Multi-task learning models outperform
single task learning approaches for each event. The Ferguson event is the
hardest one for all of the models as it has a different class distribution to
all other events (see Table 2).

Table 4 shows per event and per class performance of the multi-task learning
model that incorporates all three tasks (MTL3). We have analysed similar
performance breakdown tables for other models but omit them here because of
space constraints. All models tend to predict the majority class (true) the
best. As the Ferguson event is strongly dominated by unverified rumours, it is
the only event with high performance on the unverified class. Single task
models (NileTMRG*, branchLSTM) are better at identifying false than unverified
rumours, whereas multitask models (MTL2, MTL3) are better at identifying
unverified than false rumours.

Figure 4: Comparison of macro F-score of sequential single task approach
(branchLSTM), multi-task learning approaches with two and three tasks on
different events from the PHEME dataset.

MTL3 5 events | MacroF | Accuracy | TRUE | FALSE | Unverified
---|---|---|---|---|---
Charlie Hebdo | 0.327 | 0.369 | 0.502 | 0.227 | 0.251
Sydney siege | 0.350 | 0.575 | 0.731 | 0.153 | 0.168
Ferguson | 0.189 | 0.338 | 0.058 | 0 | 0.508
Ottawa shooting | 0.352 | 0.645 | 0.789 | 0.168 | 0.100
Germanwings-crash | 0.429 | 0.420 | 0.538 | 0.358 | 0.364

Table 4: Per event and per-class results for multi-task learning approach with
3 tasks on PHEME 5 events.

###  6.2 Analysis of data properties

Alonso and Plank [Alonso and Plank, 2017] link the gains/losses in performance
from multi-task learning with information-theoretical metrics, properties of
the label distributions: entropy (indicating the amount of uncertainty in the
distribution) and kurtosis (indicating the skewness of the distribution). They
have shown that the multi-task learning setup works best for tasks which have
label distributions with lower kurtosis and relatively high entropy. Table 5
shows the kurtosis and entropy properties of the label distribution for each
of the events in the dataset, as well as token-type ratio (TTR). Each of the
events has very different properties; there is a strong difference in
properties between larger events (top) and smaller ones (bottom). Smaller
events tend to have more extreme values, which may explain that their addition
to the evaluation adds further complexity to the task of rumour verification.
In line with findings of Alonso and Plank [Alonso and Plank, 2017], in our
case both auxiliary tasks have on average lower kurtosis than the main task.
The stance classification dataset has on average higher entropy than the
rumour detection dataset.

Whereas Alonso and Plank [Alonso and Plank, 2017] were considering low-level
linguistically related tasks as auxiliary tasks, such as part-of-speech
tagging, we are working with higher level theme-related tasks. Therefore, it
is interesting that we still observe similar trends when looking at the same
data properties.

| Kurtosis | Entropy | TTR
---|---|---|---
Events | S | V | D | S | V | D | S | V | D
charliehebdo | -0.73 | -1.25 | -0.18 | 0.89 | 1.08 | 0.53 | 0.2 | 0.11 | 0.07
ottawashooting | -0.83 | 0.33 | -1.99 | 1.04 | 0.82 | 0.69 | 0.19 | 0.11 | 0.09
germanwings-crash | -1.11 | -0.86 | -1.99 | 0.99 | 0.99 | 0.69 | 0.26 | 0.16 | 0.13
sydneysiege | -0.79 | 0.71 | -1.91 | 1.01 | 0.76 | 0.68 | 0.19 | 0.11 | 0.07
ferguson | -0.5 | 17.44 | -0.64 | 0.99 | 0.28 | 0.56 | 0.17 | 0.09 | 0.06
ebola-essien | -0.22 | -3 | -3 | 1.01 | 0 | 0 | 0.38 | 0.27 | 0.27
putinmissing | -1.55 | 9.08 | -1.98 | 1.12 | 0.26 | 0.69 | 0.49 | 0.26 | 0.24
prince-toronto | -1.05 | 27.75 | 53.26 | 1.04 | 0.14 | 0.09 | 0.36 | 0.18 | 0.18
gurlitt | - | 25.5 | -1.95 | - | 0.14 | 0.68 | - | 0.31 | 0.25

Table 5: Properties of the datasets for each of the events and each of the
tasks, where S - Stance, V - Veracity, and D - Detection.

##  7 Conclusions and Future Work

We have proposed a rumour verification model that achieves improved
performance for veracity classification by leveraging task relatedness with
auxiliary tasks, specifically rumour detection and stance classification,
through a multi-task learning approach. We have compared single task learning
approaches with the proposed multi-task learning approaches that combine the
verification classifier with the stance and rumour detection classifiers
individually, as well as with both the rumour detection system and the stance
classifier. Our results show that the joint learning of two tasks from the
verification pipeline outperforms a single-learning approach to rumour
verification. The combination of all three tasks leads to further performance
improvements. We have also investigated the link between the properties of the
label distribution in the dataset and the outcomes of our multi-task learning
models. Our results support findings from previous research [Alonso and Plank,
2017].

In future work we plan to investigate whether further improvements on main and
auxiliary tasks are possible with multi-task learning by: adapting the
training schedule to account for different dataset sizes, such that none of
the tasks dominates the model; by incorporating the hierarchy between tasks
into the model. We would like to investigate the effect of adding extra
features (such as user features and interactions) on different tasks and at
which stage to incorporate them, in private or shared layers.

## Acknowledgements

This work was supported by The Alan Turing Institute under the EPSRC grant
EP/N510129/1. Cloud computing resources were kindly provided through a
Microsoft Azure for Research Award. Work by Elena Kochkina was partially
supported by the Leverhulme Trust through the Bridges Programme and Warwick
CDT for Urban Science & Progress under the EPSRC Grant Number EP/L016400/1.

## References

  * [Aguilar et al., 2017]  Gustavo Aguilar, Suraj Maharjan, Adrian Pastor López Monroy, and Thamar Solorio.  2017\.  A multi-task approach for named entity recognition in social media data.  In Proceedings of the 3rd Workshop on Noisy User-generated Text, pages 148–153.
  * [Alonso and Plank, 2017]  Héctor Martínez Alonso and Barbara Plank.  2017\.  When is multitask learning effective? semantic sequence prediction under varying data conditions.  In 15th Conference of the European Chapter of the Association for Computational Linguistics.
  * [Bingel and Søgaard, 2017]  Joachim Bingel and Anders Søgaard.  2017\.  Identifying beneficial task relations for multi-task learning in deep neural networks.  arXiv preprint arXiv:1702.08303.
  * [Boididou et al., 2014]  Christina Boididou, Symeon Papadopoulos, Yiannis Kompatsiaris, Steve Schifferes, and Nic Newman.  2014\.  Challenges of computational verification in social multimedia.  In Proceedings of the 23rd International Conference on World Wide Web, pages 743–748. ACM.
  * [Boididou et al., 2017]  Christina Boididou, Symeon Papadopoulos, Lazaros Apostolidis, and Yiannis Kompatsiaris.  2017\.  Learning to detect misleading content on twitter.  In Proceedings of the 2017 ACM on International Conference on Multimedia Retrieval, pages 278–286. ACM.
  * [Caplow, 1947]  Theodore Caplow.  1947\.  Rumors in war.  Social Forces, pages 298–302.
  * [Caruana, 1998]  Rich Caruana.  1998\.  Multitask learning.  In Learning to learn, pages 95–133. Springer.
  * [Chen et al., 2017]  Tong Chen, Lin Wu, Xue Li, Jun Zhang, Hongzhi Yin, and Yang Wang.  2017\.  Call attention to rumors: Deep attention based recurrent neural networks for early rumor detection.  arXiv:1704.05973.
  * [Collobert and Weston, 2008]  Ronan Collobert and Jason Weston.  2008\.  A unified architecture for natural language processing: Deep neural networks with multitask learning.  In Proceedings of the 25th international conference on Machine learning, pages 160–167. ACM.
  * [Derczynski et al., 2014]  Leon Derczynski, Kalina Bontcheva, Michal Lukasik, Thierry Declerck, Arno Scharl, Georgi Georgiev, Petya Osenova, Toms Pariente Lobo, Anna Kolliakou, Robert Stewart, et al.  2014\.  Pheme: computing veracity: the fourth challenge of big social data.  In Proceedings of ESWC EU Project Networking.
  * [Derczynski et al., 2017]  Leon Derczynski, Kalina Bontcheva, Maria Liakata, Rob Procter, Geraldine Wong Sak Hoi, and Arkaitz Zubiaga.  2017\.  Semeval-2017 task 8: Rumoureval: Determining rumour veracity and support for rumours.  In Proceedings of the 11th International Workshop on Semantic Evaluation (SemEval-2017), pages 69–76.
  * [DiFonzo and Bordia, 2007]  Nicholas DiFonzo and Prashant Bordia.  2007\.  Rumor, gossip and urban legends.  Diogenes, 54(1):19–35.
  * [Enayet and El-Beltagy, 2017]  Omar Enayet and Samhaa R El-Beltagy.  2017\.  Niletmrg at semeval-2017 task 8: Determining rumour and veracity support for rumours on twitter.  In Proceedings of the 11th International Workshop on Semantic Evaluation (SemEval-2017), pages 470–474.
  * [Giasemidis et al., 2016]  Georgios Giasemidis, Colin Singleton, Ioannis Agrafiotis, Jason RC Nurse, Alan Pilgrim, Chris Willis, and Danica Vukadinovic Greetham.  2016\.  Determining the veracity of rumours on twitter.  In International Conference on Social Informatics, pages 185–205. Springer.
  * [Hu et al., 2012]  Mengdie Hu, Shixia Liu, Furu Wei, Yingcai Wu, John Stasko, and Kwan-Liu Ma.  2012\.  Breaking news on twitter.  In Proceedings of the SIGCHI Conference on Human Factors in Computing Systems, pages 2751–2754. ACM.
  * [Kochkina et al., 2017]  Elena Kochkina, Maria Liakata, and Isabelle Augenstein.  2017\.  Turing at semeval-2017 task 8: Sequential approach to rumour stance classification with branch-lstm.  In Proceedings of SemEval.ACL.
  * [Kwon et al., 2017]  Sejeong Kwon, Meeyoung Cha, and Kyomin Jung.  2017\.  Rumor detection over varying time windows.  PloS one, 12(1):e0168344.
  * [Lan et al., 2017]  Man Lan, Jianxiang Wang, Yuanbin Wu, Zheng-Yu Niu, and Haifeng Wang.  2017\.  Multi-task attention-based neural networks for implicit discourse relationship representation and identification.  In Proceedings of the 2017 Conference on Empirical Methods in Natural Language Processing, pages 1299–1308.
  * [Lewandowsky et al., 2012]  Stephan Lewandowsky, Ullrich KH Ecker, Colleen M Seifert, Norbert Schwarz, and John Cook.  2012\.  Misinformation and its correction: Continued influence and successful debiasing.  Psychological Science in the Public Interest, 13(3):106–131.
  * [Liu et al., 2016]  Xiaomo Liu, Quanzhi Li, Armineh Nourbakhsh, Rui Fang, Merine Thomas, Kajsa Anderson, Russ Kociuba, Mark Vedder, Steven Pomerville, Ramdev Wudali, et al.  2016\.  Reuters tracer: A large scale system of detecting & verifying real-time news events from twitter.  In Proceedings of the 25th ACM International on Conference on Information and Knowledge Management, pages 207–216. ACM.
  * [Lukasik et al., 2016]  Michal Lukasik, PK Srijith, Duy Vu, Kalina Bontcheva, Arkaitz Zubiaga, and Trevor Cohn.  2016\.  Hawkes processes for continuous time sequence classification: an application to rumour stance classification in twitter.  In Proceedings of the 54th Annual Meeting of the Association for Computational Linguistics (Volume 2: Short Papers), volume 2, pages 393–398.
  * [Ma et al., 2016]  Jing Ma, Wei Gao, Prasenjit Mitra, Sejeong Kwon, Bernard J Jansen, Kam-Fai Wong, and Meeyoung Cha.  2016\.  Detecting rumors from microblogs with recurrent neural networks.  In IJCAI, pages 3818–3824.
  * [Mendoza et al., 2010]  Marcelo Mendoza, Barbara Poblete, and Carlos Castillo.  2010\.  Twitter under crisis: Can we trust what we RT?  In 1st Workshop on Social Media Analytics, SOMA’10, pages 71–79.
  * [Mikolov et al., 2013]  Tomas Mikolov, Kai Chen, Greg Corrado, and Jeffrey Dean.  2013\.  Efficient estimation of word representations in vector space.  arXiv preprint arXiv:1301.3781.
  * [Procter et al., 2013]  Rob Procter, Farida Vis, and Alex Voss.  2013\.  Reading the riots on twitter: methodological innovation for the analysis of big data.  International journal of social research methodology, 16(3):197–214.
  * [Řehůřek and Sojka, 2010]  Radim Řehůřek and Petr Sojka.  2010\.  Software Framework for Topic Modelling with Large Corpora.  In Proceedings of the LREC 2010 Workshop on New Challenges for NLP Frameworks, pages 45–50, Valletta, Malta, May. ELRA.  http://is.muni.cz/publication/884893/en.
  * [Shu et al., 2017]  Kai Shu, Amy Sliva, Suhang Wang, Jiliang Tang, and Huan Liu.  2017\.  Fake news detection on social media: A data mining perspective.  ACM SIGKDD Explorations Newsletter, 19(1):22–36.
  * [Wang, 2017]  William Yang Wang.  2017\.  Liar, liar pants on fire: A new benchmark dataset for fake news detection.  In Proceedings of the 55th Annual Meeting of the Association for Computational Linguistics (Volume 2: Short Papers), volume 2, pages 422–426.
  * [Zhao et al., 2015]  Zhe Zhao, Paul Resnick, and Qiaozhu Mei.  2015\.  Enquiring minds: Early detection of rumors in social media from enquiry posts.  In Proceedings of the 24th International Conference on World Wide Web, pages 1395–1405. International World Wide Web Conferences Steering Committee.
  * [Zubiaga et al., 2016]  Arkaitz Zubiaga, Maria Liakata, Rob Procter, Geraldine Wong Sak Hoi, and Peter Tolmie.  2016\.  Analysing how people orient to and spread rumours in social media by looking at conversational threads.  PloS one, 11(3):e0150989.
  * [Zubiaga et al., 2017]  Arkaitz Zubiaga, Maria Liakata, and Rob Procter.  2017\.  Exploiting context for rumour detection in social media.  In International Conference on Social Informatics, pages 109–123. Springer.
  * [Zubiaga et al., 2018a]  Arkaitz Zubiaga, Ahmet Aker, Kalina Bontcheva, Maria Liakata, and Rob Procter.  2018a.  Detection and resolution of rumours in social media: A survey.  ACM Comput. Surv., 51(2):32:1–32:36, February.
  * [Zubiaga et al., 2018b]  Arkaitz Zubiaga, Elena Kochkina, Maria Liakata, Rob Procter, Michal Lukasik, Kalina Bontcheva, Trevor Cohn, and Isabelle Augenstein.  2018b.  Discourse-aware rumour stance classification in social media using sequential classifiers.  Information Processing & Management, 54(2):273–290.


# Impact of Galaxy Clusters on Propagation of Ultra-High-Energy Cosmic Rays

Antonio Condorelli Université Paris-Saclay, CNRS/IN2P3, IJCLab, 91405 Orsay,
France  Jonathan Biteau Université Paris-Saclay, CNRS/IN2P3, IJCLab, 91405
Orsay, France  Institut universitaire de France (IUF)  Remi Adam  Université
Côte d’Azur, Observatoire de la Côte d’Azur, CNRS, Laboratoire Lagrange,
France

(Received June 19, 2023; Revised September 8, 2023)

###### Abstract

Galaxy clusters are the universe’s largest objects in the universe kept
together by gravity. Most of their baryonic content is made of a magnetized
diffuse plasma. We investigate the impact of such magnetized environment on
ultra-high-energy-cosmic-ray (UHECR) propagation. The intracluster medium is
described according to the self-similar assumption, in which the gas density
and pressure profiles are fully determined by the cluster mass and redshift.
The magnetic field is scaled to the thermal components of the intracluster
medium under different assumptions. We model the propagation of UHECRs in the
intracluster medium using a modified version of the Monte Carlo code SimProp,
where hadronic processes and diffusion in the turbulent magnetic field are
implemented. We provide a universal parametrization that approximates the
UHECR fluxes escaping from the environment as a function of the most relevant
quantities, such as the mass of the cluster, the position of the source with
respect to the center of the cluster and the nature of the accelerated
particles. We show that galaxy clusters are an opaque environment especially
for UHECR nuclei. The role of the most massive nearby clusters in the context
of the emerging UHECR astronomy is finally discussed.

cosmic ray astronomy (324) — galaxy clusters (584) — ultra high energy cosmic
rays (1733)

##  1 Introduction

Even though their existence has been known for more than a century, the nature
and origin of cosmic rays at the highest energies remains elusive.
Observations have allowed us to explore their spectral behavior and
composition in terms of atomic mass on Earth (Coleman et al., 2023), but the
sources of ultra-high-energy cosmic rays (UHECRs), i.e. cosmic rays above
$10^{18}\ \rm eV$, still remain unknown.

Rapid progress in computational high-energy astrophysics has dramatically
advanced the study of acceleration mechanisms in systems ranging from the jets
of stellar-sized objects such as gamma-ray bursts (GRBs, Sudilovsky et al.,
2013) to the large-scale shocks surrounding galaxy clusters (Norman et al.,
1995; Ryu et al., 2003; Kang et al., 1997). Galaxy clusters are the largest
virialized structures, having typical radii $R_{\rm cl}=1-2$ Mpc and total
masses $M\simeq 10^{14}-10^{15}M_{\odot}$, including both baryonic and dark
matter. Strong turbulent magnetic fields, with root mean square values
$B\simeq$ few $\mu$G, are present inside clusters, having typical coherence
lengths of $5–30$ kpc (Bonafede et al., 2010; Donnert et al., 2018). This
implies that cosmic rays accelerated in candidate sources inside the clusters,
e.g. by hypernovae or GRBs in star-forming galaxies, or in the accretion
shocks, jets and radiolobes of active galactic nuclei (AGNs), can be confined
for long times within clusters. They can undergo interactions with the
enhanced baryonic content of the intracluster medium, whose profile is
determined by Bremsstrahlung emission in X rays (Sarazin, 1986).

Upper limits on the flux of neutrinos and gamma rays at ultra-high energies
rule out a dominant origin of UHECRs from exotic particles (Abreu et al.,
2023, 2022a), which should then originate from extragalactic astrophysical
sources. An extragalactic origin is corroborated by the observation of a
dipolar anisotropy above 8 EeV (Aab et al., 2018a) and an evidenced
correlation of UHECRs above 40 EeV with extragalactic objects in the nearby
universe (Aab et al., 2018b; Abreu et al., 2022b). Some of these extragalactic
sources could be hosted or shadowed by clusters.

UHECR propagation in a specific cluster (e.g. the Virgo cluster) has been
already treated in different works (Dolag et al., 2009; Kotera et al., 2009;
Harari et al., 2016; Fang & Murase, 2018). Although some of these theoretical
works suggested that galaxy clusters are efficient UHECR calorimeters, some
authors recently claimed excesses of UHECRs from these structures (Ding et
al., 2021; Abbasi et al., 2021). Revisiting the propagation of UHECRs in
galaxy clusters is thus a timely topic. In the following, we evaluate whether
UHECRs can escape from such environments and how clusters should be accounted
for in UHECR astronomy. We provide in particular a single parametrization of
the escaping flux, which depends on the mass of the cluster and on the UHECR
features, such as energy and atomic mass.

The paper is organized as follows: we introduce the relevant properties of
galaxy clusters and detail the way we compute the most important macroscopic
quantities for our study in Section 2; the microphysics of UHECR propagation
in such environments is detailed in Section 3; we present our results and
discuss the impact of our assumptions in Section 4. We finally draw our
conclusions in Section 5.

##  2 Intracluster medium modeling

Clusters of galaxies and the filaments that connect them are the largest
structures in the present universe in which the gravitational force due to the
matter over-density overcomes the expansion of the universe. Massive clusters
have typical total masses of the order of $10^{15}\ M_{\odot}$, mostly in the
form of dark matter ($70-80$ % of the total mass), while baryonic matter is
harbored by galaxies (few %) and composes the hot ($T\sim 10^{8}$ K) and
tenuous ($n_{\rm gas}\simeq 10^{-1}-10^{-4}\,\rm cm^{-3}$) gas ($15-20$ %)
that forms the intracluster medium (ICM, Voit, 2005). To model UHECR
propagation in this environment, we need estimates of the gas density profile,
of the magnetic field profile and of the coherence length. While the gas
density is well understood and routinely derived from X-ray observations, this
is not the case for the magnetic field, for which only a handful of
measurements – sometimes model-dependent – are available in the literature
(Vacca et al., 2018). From theoretical arguments, however, the magnetic-field
strength is often assumed to scale with the ICM thermal density or pressure.

###  2.1 Density profile

An interesting feature of galaxy clusters is that they are self-similar
objects at first order, so that their physical properties can be fully
described given their mass and redshift (Kaiser, 1986). For instance, their
universal pressure profiles (UPP) and universal density profiles (UDP) are now
well constrained from observations (e.g., Arnaud et al., 2010; Pratt et al.,
2022). Following Arnaud et al. (2010), we use the UPP expressed as

$P(x)=\frac{P_{0}\times P_{500}(M_{500},z)\times
f(M_{500},z)}{(c_{500}x)^{\gamma_{\rm UPP}}\cdot(1+(c_{500}x)^{\alpha_{\rm
UPP}})^{\frac{\beta_{\rm UPP}-\gamma_{\rm UPP}}{\alpha_{\rm UPP}}}},$ (1)

with $P_{500}(M_{500},z)$ the self-similar normalization (Nagai et al., 2007),
$f(M_{500},z)$ a small mass-dependence correction, and where $P_{0}$,
$c_{500}$, $\alpha_{\rm UPP}$, $\beta_{\rm UPP}$, $\gamma_{\rm UPP}$ are
parameters that describe the shape of the profile as a function of the scaled
radius $x=r/R_{500}$.111The mass $M_{500}$ is defined within $R_{500}$, the
radius within which the cluster density is 500 times the critical density of
the universe at the cluster redshift.

Similarly, we use the UDP as measured by Pratt et al. (2022), which can be
expressed as

$n(x)=\frac{A(M_{500},z)\times
f_{0}}{\left(x/x_{s}\right)\left(1+(x/x_{s})^{\gamma_{\rm
UDP}}\right)^{\frac{3\beta_{\rm UDP}-\alpha_{\rm UDP}}{\gamma_{\rm UDP}}}}.$
(2)

The quantity $A(M_{500},z)$ describes the normalization as a function of mass
and redshift, and the parameters $f_{0}$, $x_{s}$, $\alpha_{\rm UDP}$,
$\beta_{\rm UDP}$, $\gamma_{\rm UDP}$ describe the shape (see also Ghirardini
et al., 2019, for another calibration of the UDP).

The gas density and pressure profile are expected to be connected. This
provides us with an alternative way to describe the ICM thermal density given
the pressure profile. Assuming a polytropic relation between gas density and
pressure, using a sample of massive nearby clusters, Ghirardini et al. (2019)
measured

$P(x)=C\times n(x)^{k},$ (3)

where $k=1.19$ and $C$ is a normalization constant.

With the gas density in hand, we can derive the electron-, proton- and helium-
density profiles by scaling through the mean molecular weights $\mu_{\rm
gas}=0.61$, $\mu_{e}=1.16$, $\mu_{\rm p}=1.39$ and $\mu_{\rm He}=14.6$ (see,
e.g., Adam et al., 2020). The proton density profile of the Coma cluster, as
obtained from the best-fit model describing the ROSAT data (Briel et al.,
1992), is shown in Figure 1. It is compared to the model derived from our
methodology, using the mass and redshift from the MCXC catalog (Piffaretti et
al., 2011). The red line gives our reference model, i.e., the one obtained
using the UPP profile combined with the polytropic relation. For further
comparison, the UDP profile, as calibrated by Pratt et al. (2022) and
Ghirardini et al. (2019), are given in green and orange, respectively. We can
observe that the main differences between the data and the models, and among
the models themselves, arise in the central part of the cluster. This reflects
the increased intrinsic scatter among the cluster population relative to the
self-similar approximation in the cluster cores, while the consistency
significantly improves at $r\in[0.2R_{500},R_{500}]$ (see, e.g., Ghirardini et
al., 2019, for details). More specifically in Figure 1, the Coma cluster is a
merging system with a very flat core, thus presenting a smaller central
density than that given in our mean model (we also refer to the appendix for
further examples). The impact of the choice of the reference density model on
our final results is discussed in Section 4.2.

Figure 1: Thermal proton density profiles for the Coma cluster: the blue line
gives the best-fit model to the ROSAT data (Briel et al., 1992), the red line
gives the density obtained from the Ghirardini et al. (2019) polytropic
relation combined with the Arnaud et al. (2010) UPP profile, the green and
orange lines give the UDP from Pratt et al. (2022) and Ghirardini et al.
(2019), respectively.

###  2.2 Magnetic-field profile

The profile of magnetic-field strength can be scaled to the thermal-gas
density under several assumptions. Assuming the magnetic energy density to be
proportional to the thermal energy, we have

$\langle B^{2}(r)\rangle=2\mu_{0}P(r)/\beta_{\rm pl},$ (4)

with $\mu_{0}$ the vacuum permeability. For the plasma, we set $\beta_{\rm
pl}=200$ following the results by Walker et al. (2017) on the Perseus cluster.
By combining the central magnetic field of the Coma cluster measured by
Bonafede et al. (2010) and the central pressure obtained by Planck
Collaboration et al. (2013), we would instead estimate $\beta_{\rm pl}=77$.
Alternatively, assuming that the magnetic field is frozen into the plasma and
amplifies under pure adiabatic compression with magnetic flux conservation, we
have

$\langle B^{2}(r)\rangle=B_{\rm ref}^{2}\left(\frac{n_{\rm gas}(r)}{n_{\rm
gas}(r_{\rm ref})}\right)^{4/3}.$ (5)

The normalization $B_{\rm ref}$, taken at the radius $r_{\rm ref}$ is defined
using the reference Coma cluster, for which detailed measurement are available
in Bonafede et al. (2010).

In Figure 2, we compare the magnetic-field profile of the Coma cluster
estimated from Faraday rotation measures (Bonafede et al., 2010) to our
models. The red line gives the profile estimated using Equation 4, with
$\beta_{\rm pl}=200$, combined with the UPP from Arnaud et al. (2010). The
orange line is based instead on $\beta_{\rm pl}=77$. The green line uses
Equation 5 with the density estimated from the UDP calibrated by Pratt et al.
(2009). We observe that despite strong assumptions involved in our modeling,
the prediction follows relatively well the measurement. This is also the case
in the inner region of the cluster, where the environment is expected to play
a major role for UHECR propagation.

Figure 2: Magnetic field strength profiles for Coma cluster: the blue line
(and shaded region) show the best-fit model (and the constrained range) as
obtained from Faraday rotation measure (Bonafede et al., 2010), the red and
orange lines give the model obtained when scaling the magnetic-field strength
to the UPP and using $\beta_{\rm pl}=200$ and $\beta_{\rm pl}=77$,
respectively, the orange line give our model when scaling the magnetic-field
strength to the gas density, with the UDP from Pratt et al. (2022).

In the following work, we use as a reference the UPP in Equation 1 to derive
the density through the polytropic relation in Equation 3. The reference
magnetic field is derived assuming constant magnetic to thermal energy density
(Equation 4) with $\beta_{\rm pl}=200$. We compare the different assumptions
in appendix for a set of clusters with different morphologies and discuss the
impact of these assumptions on UHECR propagation in the following section.

##  3 UHECR propagation in galaxy clusters

###  3.1 Interactions and diffusion in a cluster

We compute the typical timescales for photo-hadronic and hadronic interactions
of UHECRs in the cluster environment from a modified version of the Monte
Carlo code SimProp (see Aloisio et al., 2012, 2015, 2016). We account for
interactions with photons of the cosmic microwave and infrared backgrounds
(CMB, CIB), as well as for hadronic interactions within the ICM.

Under the assumption of a monochromatic photon field of number density
$n_{\gamma}$, the typical interaction rate between a relativistic atomic
nucleus ($A$) and a low energy photon is approximately
$\tau_{A\gamma}^{-1}\simeq c\sigma_{A\gamma}n_{\gamma}$, where
$\sigma_{A\gamma}$ represents the cross section of the process and $c$ is the
speed of light in vacuum. If a more realistic spectral energy distribution for
the photon field is considered and the dependence of the cross section on the
energy is taken into account, the interaction rate reads (Aloisio et al.,
2013):

$\frac{dN_{\text{int}}}{dt}=\frac{c}{2\Gamma}\int_{\epsilon^{\prime}_{\text{th}}}^{\infty}\sigma_{A\gamma}(\epsilon^{\prime})\epsilon^{\prime}\int_{\epsilon^{\prime}/2\Gamma}^{\infty}\frac{n_{\gamma}(\epsilon)}{\epsilon^{2}}\,d\epsilon\,d\epsilon^{\prime},$
(6)

where $\Gamma$ is the Lorentz factor of the interacting nucleus. Note that
primed symbols (e.g. $\epsilon^{\prime}$) refer to quantities in the nucleus
rest frame, whereas unmarked symbols refer to quantities in the laboratory
frame.

Though spallation processes between UHECRs and gas have negligible impact in
the extragalactic medium, their role can be substantial in the ICM considering
the effective time that relativistic particles spend in this environment. The
timescale for the spallation process reads:

$\tau_{\rm spal}={n_{\rm ICM}\,\sigma_{\rm sp}\,c}^{-1},$ (7)

where $n_{\rm ICM}$ is the ICM gas density and $\sigma_{\rm sp}$ is the cross
section for proton-proton or proton-nucleus interactions. This process has
been implemented in SimProp making use of the most recent hadronic model,
Sibyll 2.3d (Riehn et al., 2020), a hadronic event generator. Details on the
interface between the hadronic interaction model (HIM) and the in-source
version of SimProp can be found in Condorelli et al. (2023).

In addition to interactions, diffusion in the magnetic field has to be taken
into account. In fact, charged particles populating an astrophysical
environment can be confined for a long time before escaping. The diffusion
timescale reads: $t_{\rm D}=R^{2}/D$, where $R$ is the radius of the
environment , and where $D$ is the UHECR diffusion coefficient computed in the
context of quasi-linear theory (Lee et al., 2017). The expression of the
diffusion coefficient is: $D\simeq cr_{L}^{2-\delta}\,l_{\rm c}^{\delta-1}/3$,
where $r_{L}=E/qB$ is the particle Larmor radius, $l_{\rm c}$ is the coherence
length of the magnetic field and $\delta$ is the slope of the turbulence
power-spectrum, while $B$ is the strength of the turbulent magnetic field. We
assume $\delta=5/3$ as prescribed for a Kolmogorov turbulence cascade.
Following Subedi (2017) and Reichherzer et al. (2022), we additionally
consider the transition in the diffusion regime taking place when
$r_{L}\gtrsim l_{c}$. In this energy range, the diffusion coefficient is
estimated as $D=D_{0}(r_{L}/l_{c})^{2}$, where $D_{0}$ is the value of the
diffusion coefficient computed at the energy $E_{0}$ such that
$r_{L}(E_{0})=l_{c}$. At the highest energies, the particle propagates
ballistically so that the diffusion time tends to $R/c$.

Figure 3: Interaction and escaping lengths as a function of magnetic rigidity
at the center of a prototypical galaxy cluster: photo-hadronic interaction
times (dashed-dot lines), spallation times (dashed lines) and diffusion times
(solid lines) for protons (red) and nitrogen nuclei (green). The Hubble radius
(corresponding to the age of the universe) is shown as a long-dashed line.
Length scales have been calculated assuming the following parameters:
$R_{500}=1\ \rm Mpc$, $B=1\ \rm\mu G$, $l_{c}=10\ \rm kpc$, $n_{\rm
ICM}=1\cdot 10^{-4}\ \rm cm^{-3}$.

Figure 3 summarizes the typical length-scales for interactions and escape in
the source environment for a prototype cluster (see caption). The interplay
between length-scales governs the shape of the UHECR flux as well as the
nuclear composition at the escape from the cluster. The shortest length-scale
for protons is always dictated by diffusion; this means that some protons can
escape from the environment. For nuclei (e.g. nitrogen in Figure 3), photo-
interaction lengths are the shortest at high rigidities for the chosen
parameters of the cluster (see caption). Clusters with larger magnetic fields
also present higher target densities, which reduces the hadronic interaction
length and makes hadronic interactions predominant at lower rigidities.

###  3.2 Implementation of ICM propagation in SimProp

In order to model the UHECR transport in clusters, we have developed an
extension of SimProp. This software has been used so far in the context of the
extragalactic propagation of UHECRs (see for instance Aab, 2017; Halim et al.,
2023; Luce et al., 2022). SimProp implements different photo-disintegration
cross sections and different models for the CIB. In this work, we adopt TALYS
(Koning et al., 2005; Koning & Rochman, 2012) for the photo-disintegration
cross sections and the CIB model of Gilmore et al. (2012), which are both
representative of the state of the art. SimProp is a monodimensional
propagator. Assuming spherical symmetry, all the particles are propagated
along an axis of the cluster until they reach $3\times R_{500}$, a distance
beyond which the ICM has negligible impact with respect to the extragalactic
medium.

We also consider the impact of the magnetic field on UHECR propagation.
Charged particle moving through a uniform magnetic field undergo an angular
deflection upon traversing a distance, $l_{c}$, of
$\simeq\dfrac{l_{c}}{r_{L}}$. A particle of energy $E$ and charge $q=Ze$
traversing a distance $L$ suffers an overall angular deflection given by
$\theta(E,Z)\simeq\bigg{(}\dfrac{L}{l_{c}}\bigg{)}\cdot\theta$ (Hooper et al.,
2007), which depends on the properties of the environment ($B$, $L$ and
$l_{c}$) and of the particles ($E$, $Z$). Such deflections result in an
increase in the effective propagation length, $L_{\mathrm{eff}}$, in the ICM
given by (Armengaud et al., 2005):

$\dfrac{L_{\mathrm{eff}}}{L}\simeq 65\,\bigg{(}\dfrac{E/Z}{10^{20}\,\rm
eV/26}\bigg{)}^{-2}\bigg{(}\dfrac{L}{1\ \rm
Mpc}\bigg{)}\bigg{(}\dfrac{l_{c}}{10\ \rm kpc}\bigg{)}\bigg{(}\dfrac{B}{1\
\rm\mu G}\bigg{)}^{2}$ (8)

Knowing the properties of the cluster, it is possible to compute the effective
length and therefore the effective time that a particle spends in the
environment.

The propagation inside the cluster environment is determined according to the
following methodology: 1) The propagation axis is divided in a given number of
steps, $n_{\rm steps}\geq A$, with $A$ atomic mass of the injected nuclei,
sufficiently large to sample the interactions; 2) UHECRs are injected at a
given point in the cluster and the propagation is performed only along the
chosen axis; 3) The typical length-scales are dependent on the position,
according to the magnetic-field and gas density profiles. The probability of
interaction or escape changes as a function of the radius; 4) Particles are
moved to the following step if the interaction probability is smaller than the
escape one, otherwise they lose energy and their byproducts are accounted for
in the following steps of the propagation; 5) Once a particle has reached the
border, if the diffusion probability is larger than the interaction one, this
particle escapes from the cluster environment and is propagated through the
extragalactic medium; 6) Particles that spend a time greater than the age of
the universe in the environment are considered trapped and are not propagated
anymore. This is a conservative assumption: the dominant time is the minimum
between the age of the cluster and the age of the oldest accelerator inside
it, both smaller than the age of the universe.

##  4 UHECR flux escaping the ICM

Once particles escape from the magnetized environment, it is possible to
evaluate what is the impact of the ICM on the UHECR spectrum as a function of
the injection point. We inject $10^{4}$ particles logarithmically distributed
in the energy range $10^{17}-10^{21}\ \rm eV$. The results are shown in Figure
4, where the escaping fluxes are represented as a function of rigidity. The
spectra are normalized to the spectrum expected if interactions and diffusion
in the ICM were neglected. One can notice how the closer the injection is to
the nearest edge of the environment (at $\approx+1$ Mpc), the more the
escaping flux coincides with the injection spectrum: the gas and magnetic
field densities are low and the propagating particles are less affected. If
instead the UHECRs cross the center of the cluster ($y\leq 0$), the flux is
reduced at low energies due to the trapping by the magnetic field (the so-
called magnetic horizon for extragalactic propagation, e.g. Lemoine, 2005;
González et al., 2021). At the highest energies, fluctuations in the ICM
transmission are an artifact of the normalization procedure, in a regime where
interactions of UHECRs with the CMB are important.

Figure 4: Escaping proton spectra from a cluster of $M=10^{14}\ M_{\odot}$ as
a function of the injection point. Positive (negative) $y$ valued correspond
to positions closer to the nearest (furthest) edges of the cluster. The
spectra are normalized to that expected without interactions in the ICM. The
vertical line shows the ankle energy.

More massive clusters present more intense magnetic field at the center of the
cluster, which shortens the magnetic horizon of UHECRs. The impact of the
cluster magnetic field on the propagation of UHECRs is illustrated in Figure
5, where the escaping fluxes are shown as a function of rigidity, in this case
assuming only protons at the injection.

Figure 5: Escaping proton spectra from cluster of different magnetic field
values at the center, taken at 1 kpc (see legend), assuming injection at the
center of the cluster. The spectra are normalized to that expected without
interactions in the ICM. The vertical line shows the ankle energy.

###  4.1 Parametrization of the UHECR escaping flux

We provide a parametrization of the escaping fluxes as a function of the mass
of the cluster $M$, of the position of injection point $y$ and of the nature
of the accelerated particles (protons or nuclei) in order to describe the
escaping flux above the ankle. Four representative nuclear masses are studied:
1H, 4He, 14N, 28Si. The contribution from iron nuclei is neglected, as few, if
any, are expected from simple cosmological models that describe data from the
Pierre Auger Observatory (Aab, 2017; Halim et al., 2023; Luce et al., 2022).

We notice in Figure 4 that a cluster mostly affects the escaping spectrum when
UHECRs cross its center. In fact, it is the place where the magnetic field is
most intense and where the target density is the highest. For this reason,
sources placed at $y\leq 0$ would have an escaping fluxes shaped by the
propagation in the cluster environment, while the effect is weaker for host
sources placed at $y>0$ where the traversed magnetic field is milder. For this
reason, we assume clusters to be transparent for accelerators placed at $y>0$,
while we provide a single parametrization of the transparency of the clusters
for $y\leq 0$. We define the transparency $f(R)$ of a given cluster as the
escaping flux divided by the one expected without interactions in the ICM. We
approximate the transparency as a function of rigidity $R$ by a broken power-
law, with full transparency at the highest energies:

$\log\ f(R)=\begin{cases}\Gamma\ \log(R/\rho)&R\leq\rho,\\\
{0}&R\geq\rho.\end{cases}$ (9)

We notice that, in our equation above, the break rigidity, $\rho$, depends on
the mass of the cluster, $M$, following to first order:

$\log{\rho}=\log\rho_{0}+\xi\log(M/10^{15}M_{\odot}).$ (10)

We parametrize the low-rigidity slope, $\Gamma$, of the transparency function
so that it reaches a maximum value of $2$ at high cluster masses and softens
at lower masses:

$\Gamma=\dfrac{2}{1+\bigg{(}\dfrac{M}{M_{\rm free}}\bigg{)}^{-\sigma}}.$ (11)

We find that $\log(M_{\rm free}/M_{\odot})=14.4\pm 0.5$ is consistent with the
transparency functions of both nuclei and protons. The parameter $\sigma$
governs the evolution of the index with cluster mass. Also in this case, we
find a common value $\sigma=0.25\pm 0.10$ for both nuclei and protons.

The parameters are determined by fitting the model in Equation 9 to the
escaping fluxes for different position of the sources at $y\leq 0$ and for
different cluster masses, considering either protons only or nuclei only. We
find best parameter values of $\log(\rho_{0}/\rm V)=20.0\pm 0.2$ for protons
and $\log(\rho_{0}/\rm V)=24.3\pm 0.3$ for nuclei, while $\xi=0.6\pm 0.1$ for
protons and $\xi=1.7\pm 0.2$ for nuclei. The comparison to simulated data is
shown in Figure 6, for an injection at the center of the cluster.

The two parameters that influence the rigidity at which the transition
happens, i.e. $\rho_{0}$ and $\xi$, are larger for nuclei than for protons.
This is due to the fact that the nuclei interact more than protons in the ICM,
as shown in Figure 3; for this reason the transition to $f(R)=1$ happens at
higher rigidities for nuclei.

Figure 6: Transparency as a function of energy for protons, helium and
nitrogen nuclei for different cluster magnetic field (see legend), assuming an
an injection point at the center of the environment. The points show the
results obtained from the simulations with errors resulting from the number of
injected particles. The solid lines display the proposed parametrization. The
vertical line shows the ankle energy.

Our simulations show that clusters of mass $M=10^{14}M_{\odot}$ or
$M=10^{15}M_{\odot}$, with central magnetic fields of $3$ and $9\,\mu$G
respectively, are able to trap nearly all protons up to the ankle. For lower
magnetic fields, the effect of the ICM on protons is quite negligible above
the ankle. Similar conclusions can be drawn for nuclei; nonetheless it is
important to stress that they are fully disintegrated up to at least the ankle
for clusters with central magnetic fields larger than $1\,\mu$G.

The proposed parametrization describes well the impact of the galaxy cluster
on the escaping flux above the ankle. The approximation of considering the
environment as transparent for sources at $y>0$ describes well the results of
the simulation for weakly magnetized clusters. For clusters with $B\geq 3\
\mu$G, the proposed parametrization for $y>0$ overestimates the escaping
fluxes on average by 0.4 dex.

Overall, we can conclude that only a few percent of nuclei can escape from
clusters with $B\geq 1\,\mu$G up to energies of $10^{19}$ eV for He and
$10^{19.5}$ eV for N. Protons are strongly suppressed as well in the most
massive clusters: only 40% escape at $10^{19}\ \rm eV$ for a central magnetic
field of $3\,\mu$G while practically none escape at this energy for $B\geq
9\,\mu$G. Galaxy clusters are thus hostile environments for UHECRs. The
filtering is more intense for nuclei, which are fully disintegrated in the
most massive clusters even in the outer regions of the environment.

###  4.2 Impact of our assumptions

In this investigation, many assumptions have been made. This section aims to
discuss their impact.

The most impacting assumption is the parametrization of the diffusion time,
which is based in this work on the scaling laws expected from diffusion theory
(see Sec. 3.1). Using instead the equation 20 of Harari et al. (2014), which
is based on Monte-Carlo simulations, results in an even more opaque
environment, with a transparency reduced by a factor $\simeq 2$ above the
ankle for protons injected in the center of the cluster.

The coherence length of the magnetic field in the ICM also influences our
results. In this work, all the clusters are assumed to have the same coherence
length, $7$ kpc, based on observations of Coma. More detailed constraints on
this quantity would be instrumental in determining the UHECR transparency of
clusters on a case-by-case basis.

Another important topic to be discussed is the assumption of the magnetic-
field and gas density profiles. Instead of the reference model detailed in
Section 2, we performed the same analysis as in Section 4.1 using the best-fit
models shown as blue lines in Figure 1 and 2 (see also appendix) for three
different clusters: Virgo (Planck Collaboration et al., 2016), Coma (Adam et
al., 2021) and Perseus (Churazov et al., 2004). In the three examined
clusters, the differences in transparency are at maximum of the order of
$1\%$, irrespective of their morphology, therefore it is possible to affirm
that the assumption of a UPP does not influence the main results of this work.
The self-similar framework is largely driven by the cluster mass, which can be
a difficult quantity to measure. The accuracy of cluster-mass estimates is
thus expected to be a primary source of uncertainty. One should nonetheless
note that, under our reference approach, the magnetic-field strength scales to
first order as $B\propto\sqrt{P_{500}}\propto M_{500}^{1/3}$, so that only an
order-of-magnitude uncertainty has a strong impact on the UHECR transparency
illustrated in Figure 6.

We only considered interactions with CMB and CIB photons, neglecting the
contribution of stellar and dust-grain emission in the cluster. Harari et al.
(2016) estimated such galactic contributions to be comparable to the CIB.
Neglecting the galactic emission does not affect our results. This can be
understood looking at Figure 3, where the change of slope at around
$10^{19.6}\ \rm eV$ in the curves labeled “photo-interaction” corresponds to
the transition from lower-energy interactions with CIB to higher-energy
interactions with CMB. For protons, the CMB is the only relevant field within
the Hubble time. For nuclei, we investigated the impact of doubling the CIB
density to model the galactic emission. The transparency changes only by few
percents. For this reason, as in Harari et al. (2016); Hussain et al. (2021),
the galactic emission is neglected.

In this work, we adopt Sibyll2.3d (Riehn et al., 2020) as HIM. The systematics
related to the use of this specific HIM cannot be explored currently, because
no other HIM are currently implemented in SimProp. This investigation should
be discussed in future works.

##  5 Discussion

In this work, we develop a detailed model to explore the extent to which
galaxy clusters impact UHECR propagation. In particular, the modelling of the
cluster environment and the use of a HIM for propagation in this environment
represent novelties for UHECR propagation studies.

We work under the assumption of self-similarity. From this assumption, it is
possible to derive the important quantities for UHECR propagation, namely the
magnetic-field and the gas density profiles given the mass and the redshift of
the clusters. We find that the cluster environment acts as a high-pass filter,
allowing a fraction of UHE protons to escape while the UHE nuclei interact
with the gas and photons present in the ICM.

This work presents some advances with respect to the previous literature. The
use of a software dedicated to the treatment of the cluster environment is new
in UHECR physics; the conclusions of this work are in line with other works
which predicted that galaxy cluster are hostile environment for UHE nuclei
(Kotera et al., 2009; Harari et al., 2016; Fang & Murase, 2018), while they
have a weaker although non-negligible effect on the propagation of UHE
protons. For example, Harari et al. (2016) suggest in their Figure 6 that, for
a cluster with a central magnetic field of $B=1\ \mu G$ and coherence length
of $l_{\rm c}=10\ \rm kpc$, the environment is completely transparent for
protons above the ankle energy ($>5\cdot 10^{18}\ \rm eV$), while it affects
slightly the escaping flux of intermediate-mass nuclei (70% transparency for
carbon nuclei at 10 EeV) and heavy nuclei (9% for iron nuclei at 10 EeV). In
our case, neglecting the differences in the profiles and parametrization of
the diffusion coefficient, we can compare these results with those obtained
for a structure of mass $M=10^{13}\ M_{\odot}$, which corresponds to a central
magnetic field of ${\simeq}\,1\,\mu G$. The two results are comparable for
protons, while our model predicts a significantly larger depletion of
intermediate nuclei: 10% transparency in our work instead of 70% in Harari et
al. (2016). We confirmed through simulations that this difference arises from
the different treatment of the hadronic interactions; the use of a HIM,
instead of a simple analytical model, can strongly influence the propagation
of the cascade, affecting both the fraction of energy lost and the
fragmentation of heavy nuclei. No direct comparison can be performed with the
work of Fang & Murase (2018), which shows the cumulative spectra of UHECRs at
Earth that escaped from a population of sources, nor with that of Kotera et
al. (2009), where the UHECR flux escaping from a single cluster is arbitrarily
normalized.

The present work leads to important conclusions for the emerging field of
UHECRs astronomy. Two different trends can be observed in the mass composition
of UHECRs measured with the Pierre Auger Observatory (Yushkov, 2020): a
transition from heavy to light mass composition is observed up to
$10^{18.3}\,\rm eV$, while data at higher energies suggest a transition to
intermediate-heavy masses. Based on our simulations, we should not observe UHE
nuclei coming from the inner regions of massive galaxy clusters above the
ankle energy. This includes in particular the Virgo cluster, the closest
galaxy cluster to us ($d\simeq 16\ \rm Mpc$, $M\simeq 1.2\cdot 10^{14}\
M_{\odot}$, from Planck Collaboration et al., 2016). The Pierre Auger
Collaboration indeed does not see any indication of excess in this direction,
which could have an important implication as pointed in Biteau et al. (2021).
Assuming that the UHECR production rate follows the star-formation rate or
stellar mass of nearly half million of galaxies, Biteau (2021) found that the
computed skymaps should show some excess in the direction of the Virgo
cluster, not present in the observed skymaps (Abreu et al., 2022b). Our work
confirms that this tension is lifted by magnetic trapping of UHECRs in Virgo,
as was already hypothesized in Biteau et al. (2021) through a more naive
argument (confinement time greater than the ballistic one). The result of our
work reduces the discrepancies between the arrival direction model and the
data, justifying the lack of UHE nuclei in the directions of the galaxy
clusters and thus suggesting interesting pathways to investigate composition
anisotropies.

Another application of our work is related to the dipole observed by the
Pierre Auger Observatory above 8 EeV, whose direction is qualitatively
explained from the distribution of local extragalactic matter and UHECR
deflections in the Galactic magnetic field (Aab et al., 2017). The strong
contribution to the dipole from the Virgo cluster inferred e.g. by Ding et al.
(2021) assuming that the UHECR production rate follows the distribution of
matter should be significantly lowered when accounting for magnetic trapping
and shadowing in the ICM. This statement is true also for the Perseus cluster
($d\simeq 74\ \rm Mpc$, $M\simeq 5.8\cdot 10^{14}\ M_{\odot}$, from Urban et
al., 2014), in the direction of which the Telescope Array collaboration claims
an indication of excess above $5.7\cdot 10^{19}\ \rm eV$ (Telescope Array
Collaboration et al., 2021). From the analysis, it cannot be excluded that the
UHECRs come from the vicinity or outer shocked region of the cluster. This
work tends to exclude the possibility that Telescope Array and the Pierre
Auger Observatory see UHE nuclei accelerated by a host source close to the
center of the Perseus or Virgo cluster; either they have to come from the
cluster outskirts or they have to be UHE protons, both primaries or
secondaries due to the fragmentation of heavy nuclei in the environment
surrounding the accelerator.

An interesting result of this work concerns the role of filaments of the
cosmic web (Kotera & Lemoine, 2008). In these regions (Carretti et al., 2022),
the turbulent component (inferred to be at the $\simeq\rm nG$ level) is weaker
than the regular magnetic field ($\simeq 30\,\rm nG$), both much weaker than
central magnetic fields of clusters. This means that while UHECRs are trapped
in the central regions of galaxy clusters, they can escape from filaments as
stated in Kim et al. (2019). If, as suggested by the authors, UHECRs are
correlated with filaments connected to the Virgo cluster, they should escape
from galaxies in the filaments.

A possible critical aspect beyond the scope of this work concerns the
secondary production. In fact, interactions of UHECRs lead to an excess of
secondaries, namely secondary cosmic rays, neutrinos and photons, which can
escape from the environment and could be in tension with the current
measurements. It should be noted that secondary protons produced by the
fragmentation of heavy nuclei would remain trapped in the environment, so that
they would not show up at lower energies. A natural step forward in this
analysis would concern the multi-messenger connection, by taking into account
the emission and propagation of photons and neutrinos in the environment. In
this way, it would be possible to compare the escaping gamma rays with the
possible excess observed by Fermi-LAT in the direction of the Coma cluster
($d\simeq 100$ Mpc, $M\simeq 7\cdot 10^{14}$ M⊙, from Planck Collaboration et
al., 2013), as well as to determine the expected sensitivity of upcoming
gamma-ray and neutrino facilities at higher energies.

## Acknowledgments

The authors would like to thank the reviewer for constructive suggestions,
which helped to improve the quality of this work. AC and JB gratefully
acknowledge funding from ANR via the grant MultI-messenger probe of Cosmic Ray
Origins (MICRO), ANR-20-CE92-0052. AC acknowledges co-developers of SimProp
for useful feedback in the development of the software.

In this appendix, we show how the models of the gas-density and magnetic-field
profiles discussed in Section 2 match the available measurements. We use
reference clusters for which magnetic field profile estimates are available in
the literature (Vacca et al., 2018), using the mass and redshift taken from
the MCXC catalog (Piffaretti et al., 2011). For a sufficiently large sampling
of mass and dynamical state, we select the Perseus, Coma, A 194, Hydra A, A
2634 clusters. For each of them, a best-fit model (or a model that was matched
to the data, in the case of the magnetic field) is reported: Churazov et al.
(2004); Taylor et al. (2006); Walker et al. (2017) for Perseus, Briel et al.
(1992) for Coma, Govoni et al. (2017) for A 194, Vogt & Enßlin (2003); Enßlin
& Vogt (2006); Laing et al. (2008); Kuchar & Ensslin (2011) for Hydra A and
Schindler & Prieto (1997) for A 2634. These are compared to experimental data
from the Plank/ROSAT project (Eckert et al., 2012), the XCOP project (Eckert
et al., 2017), and the ACCEPT data base (Cavagnolo et al., 2009) whenever
available and to the models outlined in Section 2. The red lines correspond to
the reference model of this work: UPP and polytropic assumption for
determining the thermal gas density, constant magnetic to thermal energy
density for magnetic field. The other possible choices are detailed in Section
2.

In Figure 7.1, the best-fit models of the thermal-gas density profiles are
matched to X-ray data. These models are beta models (or the sum of beta
models), which present a flat core by construction and thus may be over-
simplistic. Magnetic-field estimates derive from the matching of cluster
simulations to radio data and/or Faraday-rotation measure, both generally
assuming a scaling between the magnetic-field strength and the thermal-gas
density. The uncertainties in these estimations are often not well defined.
Some of the magnetic-field profiles from the literature may not be fit to the
data. We conclude that our reference model is in acceptable agreement with
measurements from the literature in the explored mass and redshift range,
given the above-mentioned caveats in the measurement uncertainties and
simplifying modeling assumptions.

Fig. Set7. Thermal-proton and magnetic-field densities of reference clusters.

Figure 7.1: Perseus cluster ($M\simeq 5.8\cdot 10^{14}\ M_{\odot}$) : typical
relaxed cluster with a dense core. The complete figure set (11 images) is
available in the online journal.

## References

  * Aab et al. (2017) Aab, A., et al. 2017, Science, 357, 1266, doi: 10.1126/science.aan4338
  * Aab et al. (2018a) Aab, A., Abreu, P., Aglietta, M., et al. 2018a, The Astrophysical Journal Letters, 853, L29, doi: 10.3847/2041-8213/aaa66d
  * Aab et al. (2018b) —. 2018b, The Astrophysical Journal, 868, 4, doi: 10.3847/1538-4357/aae689
  * Aab (2017) Aab, A. e. a. 2017, Journal of Cosmology and Astroparticle Physics, 2017, 038–038, doi: 10.1088/1475-7516/2017/04/038
  * Abbasi et al. (2021) Abbasi, R. U., et al. 2021, arXiv e-prints.  https://arxiv.org/abs/2110.14827
  * Abreu et al. (2022a) Abreu, P., Aglietta, M., Albury, J. M., et al. 2022a, The Astrophysical Journal, 933, 125, doi: 10.3847/1538-4357/ac7393
  * Abreu et al. (2022b) —. 2022b, The Astrophysical Journal, 935, 170, doi: 10.3847/1538-4357/ac7d4e
  * Abreu et al. (2023) —. 2023, Phys. Rev. Lett., 130, 061001, doi: 10.1103/PhysRevLett.130.061001
  * Adam et al. (2021) Adam, R., Goksu, H., Brown, S., Rudnick, L., & Ferrari, C. 2021, A&A, 648, A60, doi: 10.1051/0004-6361/202039660
  * Adam et al. (2020) Adam, R., Goksu, H., Leingartner-Goth, A., et al. 2020, Astron. Astrophys., 644, A70, doi: 10.1051/0004-6361/202039091
  * Aloisio et al. (2013) Aloisio, R., Berezinsky, V., & Grigorieva, S. 2013, Astroparticle Physics, 41, 73, doi: 10.1016/j.astropartphys.2012.07.010
  * Aloisio et al. (2016) Aloisio, R., Boncioli, D., di Matteo, A., et al. 2016, arXiv e-prints.  https://arxiv.org/abs/1602.01239
  * Aloisio et al. (2015) —. 2015, arXiv e-prints.  https://arxiv.org/abs/1505.01347
  * Aloisio et al. (2012) Aloisio, R., Boncioli, D., Grillo, A., Petrera, S., & Salamida, F. 2012, Journal of Cosmology and Astroparticle Physics, 2012, 007–007, doi: 10.1088/1475-7516/2012/10/007
  * Armengaud et al. (2005) Armengaud, E., Sigl, G., & Miniati, F. 2005, Phys. Rev. D, 72, 043009, doi: 10.1103/PhysRevD.72.043009
  * Arnaud et al. (2010) Arnaud, M., Pratt, G. W., Piffaretti, R., et al. 2010, A&A, 517, A92, doi: 10.1051/0004-6361/200913416
  * Biteau (2021) Biteau, J. 2021, The Astrophysical Journal Supplement Series, 256, 15, doi: 10.3847/1538-4365/ac09f5
  * Biteau et al. (2021) Biteau, J., Marafico, S., Kerfis, Y., & Deligny, O. 2021, Cosmographic model of the astroparticle skies, arXiv, doi: 10.48550/ARXIV.2108.10775
  * Bonafede et al. (2010) Bonafede, A., Feretti, L., Murgia, M., et al. 2010, Astronomy and Astrophysics, 513, A30, doi: 10.1051/0004-6361/200913696
  * Briel et al. (1992) Briel, U. G., Henry, J. P., & Boehringer, H. 1992, Astronomy & Astrophysics, 259, L31
  * Carretti et al. (2022) Carretti, E., Vacca, V., O’Sullivan, S. P., et al. 2022, MNRAS, 512, 945, doi: 10.1093/mnras/stac384
  * Cavagnolo et al. (2009) Cavagnolo, K. W., Donahue, M., Voit, G. M., & Sun, M. 2009, ApJS, 182, 12, doi: 10.1088/0067-0049/182/1/12
  * Churazov et al. (2004) Churazov, E., Forman, W., Jones, C., Sunyaev, R., & Böhringer, H. 2004, Monthly Notices of the Royal Astronomical Society, 347, 29, doi: 10.1111/j.1365-2966.2004.07201.x
  * Coleman et al. (2023) Coleman, A., et al. 2023, Astroparticle Physics, 147, 102794, doi: https://doi.org/10.1016/j.astropartphys.2022.102794
  * Condorelli et al. (2023) Condorelli, A., Boncioli, D., Peretti, E., & Petrera, S. 2023, Phys. Rev. D, 107, 083009, doi: 10.1103/PhysRevD.107.083009
  * Ding et al. (2021) Ding, C., Globus, N., & Farrar, G. R. 2021, The Astrophysical Journal Letters, 913, L13, doi: 10.3847/2041-8213/abf11e
  * Dolag et al. (2009) Dolag, K., Kachelrieß, M., & Semikoz, D. V. 2009, JCAP, 01, 033, doi: 10.1088/1475-7516/2009/01/033
  * Donnert et al. (2018) Donnert, J., Vazza, F., Brüggen, M., & ZuHone, J. 2018, Space Science Reviews, 214, doi: 10.1007/s11214-018-0556-8
  * Eckert et al. (2017) Eckert, D., Ettori, S., Pointecouteau, E., et al. 2017, Astronomische Nachrichten, 338, 293, doi: 10.1002/asna.201713345
  * Eckert et al. (2012) Eckert, D., Vazza, F., Ettori, S., et al. 2012, A&A, 541, A57, doi: 10.1051/0004-6361/201118281
  * Enßlin & Vogt (2006) Enßlin, T. A., & Vogt, C. 2006, Astronomy & Astrophysics, 453, 447, doi: 10.1051/0004-6361:20053518
  * Fang & Murase (2018) Fang, K., & Murase, K. 2018, Nature Physics, 14, 396, doi: 10.1038/s41567-017-0025-4
  * Ghirardini et al. (2019) Ghirardini, V., Ettori, S., Eckert, D., & Molendi, S. 2019, A&A, 627, A19, doi: 10.1051/0004-6361/201834875
  * Ghirardini et al. (2019) Ghirardini, V., Eckert, D., Ettori, S., et al. 2019, Astronomy & Astrophysics, 621, A41, doi: 10.1051/0004-6361/201833325
  * Gilmore et al. (2012) Gilmore, R. C., Somerville, R. S., Primack, J. R., & Domínguez, A. 2012, MNRAS, 422, 3189, doi: 10.1111/j.1365-2966.2012.20841.x
  * González et al. (2021) González, J. M., Mollerach, S., & Roulet, E. 2021, Physical Review D, 104, doi: 10.1103/physrevd.104.063005
  * Govoni et al. (2017) Govoni, F., et al. 2017, Astron. Astrophys., 603, A122, doi: 10.1051/0004-6361/201630349
  * Halim et al. (2023) Halim, A. A., Abreu, P., Aglietta, M., et al. 2023, Journal of Cosmology and Astroparticle Physics, 2023, 024, doi: 10.1088/1475-7516/2023/05/024
  * Harari et al. (2014) Harari, D., Mollerach, S., & Roulet, E. 2014, Phys. Rev. D, 89, 123001, doi: 10.1103/PhysRevD.89.123001
  * Harari et al. (2016) —. 2016, Journal of Cosmology and Astroparticle Physics, 2016, 010, doi: 10.1088/1475-7516/2016/08/010
  * Hooper et al. (2007) Hooper, D., Sarkar, S., & Taylor, A. M. 2007, Astroparticle Physics, 27, 199, doi: 10.1016/j.astropartphys.2006.10.008
  * Hussain et al. (2021) Hussain, S., Alves Batista, R., de Gouveia Dal Pino, E. M., & Dolag, K. 2021, Monthly Notices of the Royal Astronomical Society, 507, 1762, doi: 10.1093/mnras/stab1804
  * Kaiser (1986) Kaiser, N. 1986, MNRAS, 222, 323, doi: 10.1093/mnras/222.2.323
  * Kang et al. (1997) Kang, H., Rachen, J. P., & Biermann, P. L. 1997, Monthly Notices of the Royal Astronomical Society, 286, 257, doi: 10.1093/mnras/286.2.257
  * Kim et al. (2019) Kim, J., Ryu, D., Kang, H., Kim, S., & Rey, S.-C. 2019, Science Advances, 5, eaau8227, doi: 10.1126/sciadv.aau8227
  * Koning & Rochman (2012) Koning, A., & Rochman, D. 2012, 113, doi: 10.1016/J.NDS.2012.11.002
  * Koning et al. (2005) Koning, A. J., Hilaire, S., & Duijvestijn, M. C. 2005, American Institute of Physics Conference Series, Vol. 769, TALYS: Comprehensive Nuclear Reaction Modeling, 1154–1159, doi: 10.1063/1.1945212
  * Kotera et al. (2009) Kotera, K., Allard, D., Murase, K., et al. 2009, The Astrophysical Journal, 707, 370, doi: 10.1088/0004-637x/707/1/370
  * Kotera & Lemoine (2008) Kotera, K., & Lemoine, M. 2008, Phys. Rev. D, 77, 123003, doi: 10.1103/PhysRevD.77.123003
  * Kuchar & Ensslin (2011) Kuchar, P., & Ensslin, T. A. 2011, Astron. Astrophys., 529, A13, doi: 10.1051/0004-6361/200913918
  * Laing et al. (2008) Laing, R. A., Bridle, A. H., Parma, P., & Murgia, M. 2008, Monthly Notices of the Royal Astronomical Society, 391, 521, doi: 10.1111/j.1365-2966.2008.13895.x
  * Lee et al. (2017) Lee, J., Wright, J., Bertelli, N., et al. 2017, Physics of Plasmas, 24, 052502, doi: 10.1063/1.4982060
  * Lemoine (2005) Lemoine, M. 2005, Phys. Rev. D, 71, 083007, doi: 10.1103/PhysRevD.71.083007
  * Luce et al. (2022) Luce, Q., Marafico, S., Biteau, J., Condorelli, A., & Deligny, O. 2022, The Astrophysical Journal, 936, 62, doi: 10.3847/1538-4357/ac81cc
  * Nagai et al. (2007) Nagai, D., Kravtsov, A. V., & Vikhlinin, A. 2007, ApJ, 668, 1, doi: 10.1086/521328
  * Norman et al. (1995) Norman, C. A., Melrose, D. B., & Achterberg, A. 1995, ApJ, 454, 60, doi: 10.1086/176465
  * Piffaretti et al. (2011) Piffaretti, R., Arnaud, M., Pratt, G. W., Pointecouteau, E., & Melin, J. B. 2011, A&A, 534, A109, doi: 10.1051/0004-6361/201015377
  * Planck Collaboration et al. (2013) Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2013, A&A, 554, A140, doi: 10.1051/0004-6361/201220247
  * Planck Collaboration et al. (2016) —. 2016, A&A, 596, A101, doi: 10.1051/0004-6361/201527743
  * Pratt et al. (2022) Pratt, G. W., Arnaud, M., Maughan, B. J., & Melin, J. B. 2022, Astron. Astrophys., 665, A24, doi: 10.1051/0004-6361/202243074
  * Pratt et al. (2009) Pratt, G. W., Croston, J. H., Arnaud, M., & Böhringer, H. 2009, A&A, 498, 361, doi: 10.1051/0004-6361/200810994
  * Reichherzer et al. (2022) Reichherzer, P., Becker-Tjus, J., Zweibel, E. G., Merten, L., & Pueschel, M. J. 2022, Monthly Notices of the Royal Astronomical Society, 514, 2658, doi: 10.1093/mnras/stac1408
  * Riehn et al. (2020) Riehn, F., Engel, R., Fedynitch, A., Gaisser, T. K., & Stanev, T. 2020, Physical Review D, 102, doi: 10.1103/physrevd.102.063002
  * Ryu et al. (2003) Ryu, D., Kang, H., Hallman, E., & Jones, T. W. 2003, The Astrophysical Journal, 593, 599, doi: 10.1086/376723
  * Sarazin (1986) Sarazin, C. L. 1986, Rev. Mod. Phys., 58, 1, doi: 10.1103/RevModPhys.58.1
  * Schindler & Prieto (1997) Schindler, S., & Prieto, M. A. 1997, A&A, 327, 37, doi: 10.48550/arXiv.astro-ph/9706239
  * Subedi (2017) Subedi, P. 2017, The Astrophysical Journal, 837, 140, doi: https://iopscience.iop.org/article/10.3847/1538-4357/aa603a
  * Sudilovsky et al. (2013) Sudilovsky, V., Greiner, J., Rau, A., et al. 2013, Astronomy & Astrophysics, 552, A143, doi: 10.1051/0004-6361/201321247
  * Taylor et al. (2006) Taylor, G. B., Gugliucci, N. E., Fabian, A. C., et al. 2006, MNRAS, 368, 1500, doi: 10.1111/j.1365-2966.2006.10244.x
  * Telescope Array Collaboration et al. (2021) Telescope Array Collaboration, Abbasi, R. U., Abu-Zayyad, T., et al. 2021, arXiv e-prints, arXiv:2110.14827, doi: 10.48550/arXiv.2110.14827
  * Urban et al. (2014) Urban, O., et al. 2014, Mon. Not. Roy. Astron. Soc., 437, 3939, doi: 10.1093/mnras/stt2209
  * Vacca et al. (2018) Vacca, V., Murgia, M., Govoni, F., et al. 2018, Galaxies, 6, 142, doi: 10.3390/galaxies6040142
  * Vogt & Enßlin (2003) Vogt, C., & Enßlin, T. A. 2003, Astronomy & Astrophysics, 412, 373, doi: 10.1051/0004-6361:20031434
  * Voit (2005) Voit, G. M. 2005, Reviews of Modern Physics, 77, 207, doi: 10.1103/RevModPhys.77.207
  * Walker et al. (2017) Walker, S. A., Hlavacek-Larrondo, J., Gendron-Marsolais, M., et al. 2017, MNRAS, 468, 2506, doi: 10.1093/mnras/stx640
  * Yushkov (2020) Yushkov, A. 2020, PoS, ICRC2019, 482, doi: 10.22323/1.358.0482