11institutetext: University of Tyumen, ul. Volodarskogo 6, 625003 Tyumen,
Russia 11email<EMAIL_ADDRESS>11email:
<EMAIL_ADDRESS>”Organization of cognitive
associative systems” LLC, ul. Gertsena 64, 625000 Tyumen, Russia 22email:
<EMAIL_ADDRESS>
# A Comparative Study of Feature Types for Age-Based Text
Classification††thanks: Supported by the grant of the President of the Russian
Federation no. MK-637.2020.9.

Anna Glazkova and Yury Egorov 11 0000-0001-8409-6457 11 0000-0001-7670-5283
Maksim Glazkov 22 0000-0002-4290-2059

###### Abstract

The ability to automatically determine the age audience of a novel provides
many opportunities for the development of information retrieval tools.
Firstly, developers of book recommendation systems and electronic libraries
may be interested in filtering texts by the age of the most likely readers.
Further, parents may want to select literature for children. Finally, it will
be useful for writers and publishers to determine which features influence
whether the texts are suitable for children. In this article, we compare the
empirical effectiveness of various types of linguistic features for the task
of age-based classification of fiction texts. For this purpose, we collected a
text corpus of book previews labeled with one of two categories – children’s
or adult. We evaluated the following types of features: readability indices,
sentiment, lexical, grammatical and general features, and publishing
attributes. The results obtained show that the features describing the text at
the document level can significantly increase the quality of machine learning
models.

###### Keywords:

Text classification Fiction Corpus Age audience Content rating Text difficulty
RuBERT Neural network Natural language processing Machine Learning.

##  1 Introduction

Nowadays, there are quite a lot of approaches to text classification according
to document subjects, genre, author or according to other attributes. However,
modern challenges in the field of natural language processing (NLP) and
information retrieval (IR) increasingly require classification based on more
complex characteristics. For example, it may be necessary to determine whether
the text contains elements of propaganda or whether it has similar plot
characteristics to other texts. One of such urgent and complex classification
tasks is the division of literary texts into suitable for children and for
adults. Age-based classification tools could find wide practical application.
For instance, they would be useful in the personal selection of fiction or in
filtering content not intended for children.

Despite the fact that many scholars considered the issue of text difficulty
estimation, formally, text difficulty does not indicate the age of the
intended reader. The question of whether the features describing text
difficulty are suitable for age-based classification needs to be investigated.
In addition, the severity of the features depends on the text genre. It is
necessary to find out how these or those features are presented in the
literary text and whether they contain information about the age audience of
the text.

In this paper, we systematically evaluated different feature types on age-
based classification task. In addition to popular text difficulty features, we
consider special publishing attributes intrinsic fiction books, such as age
rating score and abstract features. We collected the corpus of Russian fiction
texts and applied two commonly used machine learning models, these are random
forest (RF) and linear support vector classifier (LSVC). For comparison, we
evaluated a transformer model based on RuBERT and a convolutional neural
network (CNN) trained on Word2Vec embeddings. Finally, we evaluated
feedforward neural network (FNN) trained on RuBERT text embeddings and age
rating scores.

The LSVC model using a combination of a baseline and publishing attributes
showed the best result of 95.77% (F1-score). RuBERT achives 90.16%. The FNN
model combining RuBERT embeddings and age ratings showed 94.78%. The results
show that the features describing the text at the document level gives an
advantage in case of long texts. Moreover, publishing attributes provide
valuable information for the age-based classifier. We also found that some
features used to determine text difficulty positively affect the quality of
age-based classification.

The paper is organised as follows. In Section 2 we present a brief review of
related works. Section 3 describes feature types evaluated in the paper.
Section 4 contains the description of our dataset. Section 5 presents the
structure of the models and the evaluation results. Finally, Section 6 is a
conclusion.

##  2 Related Works

In the modern world, the constant growth of information resources gives rise
to the need for filtering and ranking texts. One of the significant
characteristics of a text is its complexity. The question of determining text
difficulty naturally looks related to the task of age-based classification.

The task of estimating texts by complexity is not new. It appeared at the
beginning of the last century in the context of evaluating the readability of
educational texts. Further, during the XX century, researchers have proposed a
number of tests to determine readability based on the quantitative
characteristics of texts. Readability tests usually use quantitative text
features, such as counting syllables, words, and sentences. There are several
common readability texts for text difficulty estimation. For instance, these
are: the Flesch–Kincaid readability test, the Coleman–Liau index, the
automated readability index (ARI), the SMOG grade, the Dale-Chall formula [6,
10].

  * •

The Flesch–Kincaid readability test is based on the idea that the shorter the
sentences and words, the simpler the text. The specific mathematical formula
is:

$\textit{R}_{\textit{F}}=206.835-1.015\cdot\textit{ASL}-84.6\cdot\textit{ASW},$
(1)

ASL – average sentence length, ASW – average number of syllables per word
(i.e., the number of syllables divided by the number of words).

  * •

The Coleman–Liau index uses letters instead of syllables. The formula takes
into account the average number of letters per word and the average number of
words per sentence.

$\textit{R}_{\textit{C}}=0.0588\cdot L-0.296\cdot S-15.8,$ (2)

where L – average number of letters per 100 words, S – average number of
sentences per 100 words.

  * •

The ARI formula takes into account the number of letters. In the past, this
allowed the use of this index to measure the complexity of texts in real time
in electric typewriters.

$\textit{R}_{\textit{A}}=4.71\cdot\frac{characters}{words}+0.5\cdot\frac{words}{sentences}-21.43,$
(3)

where characters – number of letters and numbers, words – number of words,
sentences – number of sentences.

  * •

The main idea of the SMOG grade is that the complexity of the text is most
affected by complex words. Complex words are words with many syllables (more
than 3). The more syllables the more complicated the word.

$\textit{R}_{\textit{S}}=1.043\cdot\sqrt{polysyllable\cdot\frac{30}{sentences}}+3.1291,$
(4)

where polysyllable – number of polysyllable words, sentences – number of
sentences.

  * •

The Dale-Chall formula uses a count of ”hard” words. These ”hard” words are
words that do not appear on a specially designed list of common words familiar
to most 4th-grade students.

$\textit{R}_{\textit{D}}=0.1579\cdot\frac{difficult}{words}\cdot
100+0.0496\cdot\frac{words}{sentences},$ (5)

where difficult – number of difficult words, words – number of words,
sentences – number of sentences.

In addition to the above, there are many other readability tests that are also
actively used, e.g. the Fry Graph readability formula, the Spache index, the
Linsear Write formula and others. The values obtained from readability tests
are called readability indices. The Readability Index characterizes the
difficulty of perceiving a text or the expected level of education that is
required to understand it.

The readability formulas listed above are metrics for English texts. At the
same time, the quantitative characteristics of other languages can differ
significantly. For instance, Russian sentences are on average shorter than
English, and words are longer. Therefore, the readability formulas need to be
processed for use in other languages. Up to now, several studies have
suggested the adaptation of readability tests for Russian. For example, I.
Oborneva [27] proposed the coefficients for the Flesch–Kincaid formula for
Russian texts. The project [41] offers the adaptation of several readability
formulas. M. Solnyshkina et al. presented a new approach to reading difficulty
prediction in Russian texts [37, 38].

Readability is however only one aspect of age-based classification. Scholars
have proposed more complex techniques for text complexity estimation using
features of different nature. Thus, Yu. Tomina [42] considered the lexical and
syntactic features of the text complexity level. A. Laposhina et al. [22]
evaluated a wide range of different types of features, such as readability,
semantic, lexical, grammatical and others. M. Shafaei et al. [33] estimated
age suitability rating of movie dialogs using genre and sentiment features. L.
Flekova et al. [16] proposed an approach to describing the story complexity
for literary text. Y. Bertills [4] wrote about the features of literary
characters and named entities in books for children. Finally, in our previous
research, we evaluated the informativeness of some quantitative and
categorical features for age-based text classification [12]. The modern
methodology for text difficulty estimation is based in most cases on machine
learning approaches. Thus, R. Balyan et al. [3] showed that applying machine
learning methods increased accuracy by more than 10% as compared to classic
readability metrics (e.g., Flesch–Kincaid formula). To date, a number of
studies confirmed the effectiveness of various machine learning techniques for
text difficulty estimation, such as support vector machine (SVM) [36, 39],
random forest [26], and neural networks [2, 7, 35].

Another aspect of assessing the age category of text readers is the safety of
the information it contains. Currently, in many countries, publishers are
required to label books (including fiction) and other informational sources
[9, 11, 14, 18, 30] according to their age rating. For these purposes, there
are special laws that rank information in terms of the potential harm it can
bring. So, in Russia there is a Russian Age Rating System (RARS).

The RARS was introduced in 2012 when the Federal law of Russian Federation no.
436-FZ of 2010-12-23 ”On Protection of Children from Information Harmful to
Their Health and Development” was passed [15]. The law prohibits the
distribution of ”harmful” information that depicts violence, unlawful
activities, substance abuse, or self-harm. The RARS includes 5 categories,
such as for children under the age of six (0+), for children over the age of
six (6+), for children over the age of twelve (12+), for children over the age
of sixteen (16+), and prohibited for children (18+). As a rule, an age rating
is assigned to a book by editors or experts. As far as we know, there are
currently no published research of how age rating correlates with other
attributes of text, such as readability.

The reviewed studies and sources clearly indicate that age-based
classification of fiction texts includes several aspects. First, the research
topic is related to works on text difficulty evaluation. Text difficulty is
characterized by different features, these are lexical, semantic, grammatical
and other types. However, the measure of the difficulty of the text does not
guarantee that this text is targeted to a particular age audience. It is
required to evaluate the effectiveness of the existing text difficulty
features for age-based classification. In addition, it would be interesting to
evaluate the role of publishing attributes (for example, age rating labels) as
classification features. Finally, the studies presented thus far provide
evidence that machine learning approaches show the highest results in the task
of estimating texts by difficulty. Based on this, it is reasonable to evaluate
the text features for age-based classification using machine learning methods.

##  3 Feature Types

According to the related works, we consider the following types of
classification features.

  1. 1.

General features. This type includes features that reflect the quantitative
characteristics of the text:

     * •

the average and median length of words (avg words len, med words len);

     * •

the average and median length of sentences (avg sent len, med sent len), e.g.
average or median number of symbols in each sentence;

     * •

the average number of syllables (avg count syl);

     * •

the percentage of long words with more than 4 syllables (many syllables);

     * •

the Type-Token Ration, TTR (ttr) [40]. The main idea of the metric is that if
the text is more complex, the author uses a more varied vocabulary so there’s
a larger number of unique words. So, the TTR’s value is calculated as the
number of unique words divided by the number of words. As a result, the higher
the TTR, the higher the variety of words;

     * •

the TTR for nouns (ttr n), adjectives (ttr a), and verbs (ttr v). The values
of TTR calculated separately for parts of speech;

     * •

the NAV metric (nav). The NAV metric is a TTR-based ratio of (TTR A + TTR
N)/TTR V proposed in [37].

  2. 2.

Readability features. We used the readability formulas with the coefficients
for the Russian language proposed by the project [41]. In this study, we
evaluated five types of readability indices using the following metrics: the
Flesch–Kincaid readability test (index fk); the Coleman–Liau index (index cl);
the ARI index (index ari); the SMOG grade (index SMOG); the Dale-Chall formula
(index dc).

  3. 3.

Lexical features. In this category, we included features constructed by the
evaluation of the text in accordance with frequency dictionaries. As frequency
dictionaries, we used the lists of Russian frequency words presented in [25,
34]:

     * •

the percentage of words included in the list of 5000 most frequent Russian
words (5000 proc);

     * •

the average frequency of the words included in the 5000 most frequent words
(5000 freq);

     * •

the average frequency of words per 1 million occurrences (ipm, words fr);

     * •

the average frequency of nouns, verbs, adjectives, adverbs and proper names
per 1 million occurrences (s fr, v fr, adj fr, adv fr, prop fr);

     * •

the average number of topic segments of the corpus111The frequency dictionary
was created on the basis of the modern subcorpus of the Main Corpus and the
Oral Corpus of the Russian National Corpus (1950-2007) [32] with a total
volume of 92 million tokens [25]. where the word was encountered (out of 100
possible, words r);

     * •

the average number of the corresponding topic segments for nouns, verbs,
adjectives, adverbs and proper names (s r, v r, adj r, adv r, prop r);

     * •

the average value of Juilland’s usage coefficients (words d). This Juilland’s
usage coefficient measures the dispersion of the word’s subfrequencies over
$n$ equally-sized subcategories of the corpus [17];

     * •

the average value of Juilland’s usage coefficients for nouns, verbs,
adjectives, adverbs and proper names (s d, v d, adj d, adv d, prop d);

     * •

the number of documents in the corpora in which a word occurs (averaged over
the text, words doc);

     * •

the average number of documents in the corpora in which a word occurs (for
nouns, verbs, adjectives, adverbs and proper names, s fr, v fr, adj fr, adv
fr, prop fr).

  4. 4.

Grammatical features. We evaluated the percentage of nouns, verbs, and
adjectives (count n, count v, count a).

  5. 5.

Sentiment features. These features obtained with Russian Sentiment Lexicon
[24]. We separately evaluated the percentage of positive and negative words
for each of the topic categories, these are opinion, feeling (private state),
or fact (sentiment connotation) (neg opinion, neg feeling, neg fact, pos
opinion, pos feeling, pos fact).

  6. 6.

Publishing features. Here we have included features based on publishing
attributes, i.e. on the book characteristics assigned by an editor or
publisher, such as age rating according to the RARS (age rating) and TF-IDF
scores for book abstracts.

##  4 Dataset

For feature evaluation, we collected a dataset of fiction books published in
Russian. Due to copyright restrictions, the full texts of the books are not
publicly available. Therefore, we used a collection of previews presented in
online libraries in the public domain. Typically, the preview is 5-10% of the
total book volume.

The corpus consists of 5592 texts of children’s and adult book previews. We
have divided the texts into two parts. The first part included 4492 texts. It
was used to train the models. The remaining 1000 texts were served as an
independent text sample. The main characteristics of the data is presented in
Table 1. Table 2 shows short text examples of adults and children’s
categories.

Table 1: Characteristics of the corpus. Characteristic | Training sample | Test sample
---|---|---
Children’s | Adult | Children’s | Adult
Number of text | 2108 | 2384 | 500 | 500
Avg number of symbols | 3134.38 | 3326.11 | 3048.69 | 3319.86
Avg number of tokens | 488.55 | 499.52 | 479.3 | 498.16
Avg number of sentences | 37.35 | 35.2 | 36.05 | 36.49
Table 2: Example short fragments. Category |  Age Rating |  Genre |  Fragment
---|---|---|---
Adults |  16+ |  Modern romance novels |  A tall young man dressed in jeans, an inconspicuous jacket and a baseball cap pulled down with a visor over his eyes, approached the entrance of a seventeen-story apartment building and stood as if waiting for someone, and when a mother with a stroller appeared at the door, he quickly jumped inside - he did not know the code. I walked up to the fifth floor, putting on thin gloves on the go, looked around, and then deftly opened the door of one of the apartments. On the threshold he froze and listened for a while, but it was quiet. The man turned the baseball cap over the visor and began a leisurely survey of the apartment, opening the doors of the cupboards and looking into the drawers. The first thing he did was to open the sliding wardrobe door in the hallway, and oversized men’s slippers, an empty box, and a bright purple scarf fell out onto the floor. The man grimaced and shoved everything back, muttering, ”I thought so.” In a large room that served as a bedroom, he lingered a little longer and grunted ironically at the sight of a luxurious couch with a carved back. 222Fragment from the book ”Men We Choose” by Evgeniya Perova (translated from Russian).
Adults |  12+ |  Histori- cal adventures |  What a wonderful autumn it was in Southern Poland that year! Almost without rain and cold winds, tenderly warm, quiet, crimson-gold. Fabulous autumn - in such an autumn it is good, having climbed into the spurs of the Beskydy, from dawn to noon to wander along the slopes of hills overgrown with beech and hazel, and to your fill, drunk to breathe in the cool and crystal clear mountain air. And then, on the cozy terrace of a small mountain tavern, eat a good portion of hot, fiery-spicy bigos with pork legs, washed down with icy ”okocim”. And in the evening, having walked up to aching knees, kindle a fire on a platform open to all the winds above a shallow ravine, and, sitting on unbound logs, look at the stars that suddenly poured out in incredible numbers overhead. And, peering to the north, in the transparent thickening blue of the air, distinguish the lights of distant Krakow or Nova Huta, or maybe Bochnia or Wieliczka, who knows? 333”Sold Poland” by Alexander Usovsky (translated from Russian).
Children’s |  6+ |  Child- ren’s adventures |  In a big, big city, where there are many, many houses, many, many cars and even more people, and the crows cannot be counted at all, there lived a ginger cat on a short street consisting of only two courtyards. His name was Ostrich. 444”Greetings from cutlets” by Evgenia Malinkina (translated from Russian).
Children’s |  12+ |  Child- ren’s fantastic tales |  The hands of the clock were approaching half past seven, but the setting sun, reluctantly sliding behind the houses, continued to burn the city with rays, and the approaching twilight did not promise the long-awaited coolness. Friday night was hot and stuffy, and the city roofs were so hot during the day that no sane cat would dare to run over them without burning their paws. August was coming to an end, and the sun knew that it was the most important thing in the city, so from the very morning it climbed everywhere, trying to melt the asphalt on the streets, drying the grass on the lawns and sneaking into the apartments to flood them with heat and stuffiness. 555”Vlad and the Secret Ghost” by Sasha Gotti (translated from Russian).

Table 3 shows the most informative quantitative features with their means and
standard deviation values. The informativeness is measured using the method of
cumulative frequencies [1, 44]. The main idea of this method consists in
dividing the range of feature values for each class into $n$ intervals. The
cumulative frequency of characteristic values is calculated for each interval.
The informativeness indicator is calculated as the maximum absolute value of
the difference in the accumulated frequencies for the corresponding intervals
in the classes.

Table 3: Top-10 of the most informative quantitative features (according to the method of cumulative frequencies). Feature | Mean (adult) | std (adult) | Mean (children’s) | std (children’s)
---|---|---|---|---
avg_sent_len | 105.65 | 54.51 | 88.69 | 30.64
med_sent_len | 97.9 | 59.06 | 79.69 | 34
index_dc | 7.85 | 2.64 | 6.36 | 2.09
adj_doc | 3484.28 | 1074.03 | 3669.48 | 1173.78
index_ari | 9.29 | 3.53 | 7.4 | 2.93
adj_fr | 135.3 | 53.56 | 146.96 | 60.39
s_doc | 3705.73 | 838.45 | 3370.52 | 766.39
index_fk | 9.21 | 3.66 | 7.36 | 2.95
v_doc | 6255.7 | 2243.76 | 6481.53 | 1937.8
adv_fr | 239.17 | 84.26 | 277.43 | 90.93

Figure 1 presents the distribution of the age rating labels (age rating is a
categorical feature) in the classes of the training data. It is interesting to
note that some of the books from the children’s class are labelled with the
18+ age category. Notable examples of this type of books are love stories for
teens.

Figure 1: The distribution of age rating categories.

##  5 Experiments

This section describes our feature evaluation experiments. We built two types
of baseline models and sequentially enriched them with different types of
features. Further, we compared the results obtained with our models with the
results of CNN and RuBERT. Our dataset and models are available at [5].

###  5.1 Baselines

We built two classifiers for model evaluation. The first one was a Random
Forest Classifier trained on bootstrap samples. The number of trees in the
forest was equal to 100 and the Gini impurity was implemented to measure the
quality of a split. The second model was a Linear Support Vector Classifier
with the ”l2” penalty and the squared hinge loss function. Both models were
implemented using Scikit-learn [29] and Python 3.6.

###  5.2 Preprocessing

To preprocess our data, we used min-max normalization. Moreover, it is obvious
that some of the features are correlated. For instance, most readability
indices show a cross-correlation greater than 0.8. Another example of
correlated feature pairs is average and median length of sentences or TTR
values for all words and for particular parts of speech (Figure 2). To reduce
the influence of feature correlation on the LSVC model, we applied linear
dimensionality reduction using Singular Value Decomposition of the data with
the minimum number of principal components such that 95% of the variance is
retained.

Figure 2: Correlation matrices for readability and general features.

###  5.3 Experiments and Results

We used models trained of TF-IDF vectors as baselines. Further, we
systematically evaluated each type of features.

To begin this process, we connected the TF-IDF vector of the text with the
corresponding vector of features of a certain type. Book previews are rather
long texts. Since the RuBERT model that participated in the comparison can
only process a sequence of limited length, we used the same fragments of 256
tokens to train both neural networks and to construct TF-IDF vectors666The
maximum sequence length for BERT is 512 tokens. However, due to the rather
large volume of the corpus, we were also limited in computational resources..
TF-IDF vectors were built over the top 2000 words ordered by term frequency
across the corpus. At the same time, the values of additional features were
calculated for the full preview texts.

To build TF-IDF vectors and CNN, the texts were pre-processed. The
preprocessing included the following steps: special character removal,
lowercase translation, lemmatization, stop-word removal. Text preprocessing
was implemented using NLTK [23] and Pymorphy2 [19].

Table 4: Age-based classification results (%). Model | Accuracy | F1-score | Precision | Recall
---|---|---|---|---
Baselines
RF (TF-IDF) | 85.8 | 86.37 | 90 | 83.03
LSVC (TF-IDF) | 83.7 | 84.01 | 85.6 | 82.47
Readability features
RF baseline + readability | 86.5 | 86.91 | 89.6 | 84.37
LSVC baseline + readability | 84.9 | 85.63 | 90 | 81.67
RF (readability) | 60.6 | 61.6 | 63.2 | 60.08
LSVC (readability) | 60.8 | 61.53 | 52.6 | 74.11
Sentiment features
RF baseline + sentiment | 85.4 | 86.04 | 90 | 82.42
LSVC baseline + sentiment | 83.8 | 84.09 | 85.6 | 82.63
RF (sentiment) | 59.4 | 62.2 | 66.8 | 58.19
LSVC (sentiment) | 68 | 67.01 | 65 | 69.15
Lexical features
RF baseline + lexical | 83.2 | 84.23 | 89.9 | 79.23
LSVC baseline + lexical | 84.2 | 84.45 | 85.8 | 83.14
RF (lexical) | 63.1 | 64.62 | 67.4 | 62.06
LSVC (lexical) | 61.8 | 62.91 | 64.8 | 61.13
Grammatical features
RF baseline + grammatical | 85.5 | 86.26 | 91 | 81.98
LSVC baseline + grammatical | 83.7 | 84.03 | 85.6 | 82.47
RF (grammatical) | 56.3 | 57.86 | 60.1 | 55.87
LSVC (grammatical) | 59.6 | 62.1 | 66.2 | 58.48
General features
RF baseline + general | 89.8 | 87.57 | 90.2 | 85.09
LSVC baseline + general | 87.9 | 88.01 | 88.8 | 87.23
RF (general) | 60.8 | 61.72 | 63.2 | 60.31
LSVC (general) | 68.8 | 70.51 | 74.6 | 66.85
Publishing attributes
RF baseline + abstracts | 87.4 | 87.77 | 90.41 | 85.28
LSVC baseline + abstracts | 85.7 | 85.89 | 87.01 | 84.8
RF baseline + age rating | 90.1 | 90.06 | 91.21 | 88.93
LSVC baseline + age rating | 89.9 | 88.94 | 88.42 | 89.46
RF baseline + publ. attr. | 90.4 | 90.94 | 96.4 | 86.07
LSVC baseline + publ. attr. | 93 | 92.9 | 91.6 | 94.24
All features
RF baseline + all features | 94.9 | 94.83 | 93.6 | 96.1
LSVC baseline + all features | 95.8 | 95.77 | 95 | 96.54
RF (all features) | 94.7 | 94.67 | 94.2 | 95.15
LSVC (all features) | 94.2 | 94.09 | 92.4 | 95.85
RF baseline+all features–publ.attr. | 86.1 | 86.41 | 88.4 | 84.51
LSVC baseline+all features–publ.attr. | 87.3 | 87.54 | 89.2 | 85.93
RuBERT | 90.5 | 90.16 | 84.28 | 93.55
FNN (RuBERT embs + age rating) | 94.8 | 94.78 | 94.4 | 95.31
CNN | 82.1 | 80.2 | 89.6 | 72.59

Table 4 shows the results obtained for each type of features (e.g. RF baseline
+ readability, LSVC baseline + readability) and the results of the models
trained only on additional features without TF-IDF vectors (e.g. RF
(readability), LSVC (readability)). For publishing attributes, we evaluated
two separate types of models. The first type used book abstracts as
supplementary information. In other words, we added the texts of abstract to
the book preview and built new TF-IDF vectors. The second type used baseline
TF-IDF vectors with an additional feature of age rating. Finally, we evaluated
three types of combined models, such as using all considered features, only
all additional features, and all features with the exception of editorial
attributes. The results obtained were compared with the results of three
neural models:

  * •

RuBERT [20], based on BERT architecture [8]. BERT showed state-of-the-art
results on a wide range of NLP tasks. RuBERT was trained on the Russian part
of Wikipedia and news data. The model was implemented using PyTorch [28] and
Transformers [43] libraries, it was trained for 3 epoches;

  * •

FNN trained on fine-tuned RuBERT [20] text embeddings obtained with PyTorch
library [28]. Text embeddings were calculated by averaging the token vectors
of the last hidden state. Age rating was presented as a one-hot numeric array
which is the most widely used coding scheme [31]. FNN consisted of three
layers including an input layer, a 1024 hidden layer with hyperbolic tangent
activation function, and an output layer with softmax activation function. We
also used Adam as an optimizer and binary cross-entropy loss. The FNN model
was implemented using Keras [13] library;

  * •

CNN trained on Word2Vec embeddings [21]. CNN consisted of four building units
including three convolutional units (CU) and a fully connected unit. Each CU
contained the following sequence of layers $C-BN-C-BN-P$ where $C$ is a
convolutional layer (CL), $BN$ is a batch normalization layer, $P$ is a
pooling layer. After every CL the LeakyRelu activation function was applied.
At the first CU we used 512 filters $7\times 7$. At the second CL we applied
1024 filters $5\times 5$. As a pooling strategy at each layer we used max
polling with a kernel $2\times 2$. The fully connected layer consisted of the
following sequence of layers $FL_{1}-BN-FN_{2}$ where $FL_{1}$ is a hidden
layer with 32 neurons, $FL_{2}$ is an output layer. We applied ReLU as an
activation function and used stochastic gradient descent with Nesterov
momentum and learning rate equal to $5\times 10^{-2}$ as optimization
parameters. The model was implemented with PyTorch [28].

The results show that additional features in most cases improve the quality of
baselines. According to F1-score, this concerns readability features, general
features and publishing attributes. We assume that the advantage of these
features is that they describe the text at the document level and allow the
model to evaluate the whole text, and not just a fragment. It also can be seen
that the using of abstracts and the age rating feature significantly improves
the quality of the classification. The best results was obtained by the LSVC
model using all considered features (95.8% of accuracy, 95.77% of F1-score,
95% of precision, and 96.54% of recall). These values are shown in bold in
Table 3. Among the models that did not use publishing attributes, the best
results were shown by RuBERT (90.5% of accuracy, 90.16% of F1-score, and
93.55% of recall) and the LSVC baseline with grammatical features (91% of
precision).

##  6 Conclusion

The purpose of the current study was to evaluate different types of features
for the task of age-based text classification. The results of this
investigation show that features used in text difficulty evaluation can
improve the quality of age-based classification. In addition, in this study,
we considered publishing attributes (such as book abstracts and age ratings)
as classification features. The results showed that the use of these
attributes in digital libraries and recommendation systems could significantly
improve the quality of machine learning approaches. Our further research will
focus on studying other types of features, such as named entity analysis or
plot and character features.

## References

  * [1] Aivazyan, S.A., Bukhshtaber, V.M., Enyukov, I.S. et al.: Applied Statistics: Classification and Dimension Reduction: A Handbook. Finansy i statistika, Moscow (1989).
  * [2] Azpiazu, I.M., Pera, M.S.: Multiattentive recurrent neural network architecture for multilingual readability assessment. Trans. of the Association for Comp. Ling. 7, 421–436. /10.1162/tacl a 00278
  * [3] Balyan, R., McCarthy. K.S., McNamara, D.S.: Applying Natural Language Processing and Hierarchical Machine Learning Approaches to Text Difficulty Classification. Int. J. of Art. Intelligence in Education, 1–34 (2020). 10.1007/s40593-020-00201-7
  * [4] Bertills, Y.: Beyond identification : proper names in children’s literature. Abo Akademi University Press, Turku (2003).
  * [5] Corpus and Baselines for Age-Based Text Classification, https://github.com/oldaandozerskaya/age_based_classification. Last accessed 24 Sep 2020.
  * [6] Crossley, S.A., Skalicky, S., Dascalu, M. et al.: Predicting text comprehension, processing, and familiarity in adult readers: New approaches to readability formulas. Discourse Processes 54, 5–6 (2017). 10.1080/0163853x.2017.1296264
  * [7] Cuzzocrea, A., Bosco, G. L., Pilato, G. et al.: Multi-class text complexity evaluation via deep neural networks. LNCS 11872, 313–322 (2019). 0.1007/978-3-030-33617-2 32
  * [8] Devlin, J., Chang, M.W., Lee, K. et al.: Bert: Pre-training of deep bidirectional transformers for language understanding. arXiv preprint arXiv:1810.04805 (2018).
  * [9] Dogruel, L., Joeckel, S.: Video game rating systems in the US and Europe: Comparing their outcomes. International Communication Gazette 757, 672–692 (2013).
  * [10] Didegah, F., Thelwall, M.: Which factors help authors produce the highest impact research? Collaboration, journal and document properties. J. of informetrics 7(4), 861–873 (2013). 10.1016/j.joi.2013.08.006
  * [11] Grealy, L., Driscoll, C., Cather, K.: A history of age-based film classification in Japan. Japan Forum (2020). 10.1080/09555803.2020.1778058
  * [12] Glazkova, A.: An Approach to Text Classification based on Age Groups of Addressees. SPIIRAS Proceedings 52(3), 51–69 (2017). 10.15622/sp.52.3
  * [13] Gulli, A., Pal, S. Deep learning with Keras. Packt Publishing Ltd (2017).
  * [14] Kim, S.W. et al. A Global Comparative Study on the Game Rating System. Journal of Digital Convergence 1712, 91-108 (2019).
  * [15] Federal Law of Dec. 29, 2010 N 436-FZ ”On the Protection of Children from Information Harmful to Their Health and Development”, http://www.consultant.ru/document/cons_doc_LAW_108808/. Last accessed 23 Jul 2020.
  * [16] Flekova, L., Stoffel, F., Gurevych, I. et al.: Content-based Analysis and Visualization of Story Complexity. In: Visualisierung sprachlicher Daten, pp. 185–223. Heidelberg: Heid. Univ. Publishing (2018).
  * [17] Juilland, A.G., Brodin, D.R., Davidovitch, C.: Frequency dictionary of French words. Hague, Paris (1971).
  * [18] Hamid, R.S., Shiratuddin, N.: Age Classification of the Existing Digital Game Content Rating System Across the World: A Comparative Analysis. In Proceedings of KMICe, pp. 218-222 (2018).
  * [19] Korobov, M.: Morphological analyzer and generator for Russian and Ukrainian languages. In: International Conference on Analysis of Images, Social Networks and Texts, pp. 320-332. Springer, Cham (2015). 10.1007/978-3-319-26123-2 31
  * [20] Kuratov, Y., Arkhipov, M.: Adaptation of deep bidirectional multilingual transformers for Russian language, arXiv preprint arXiv:1905.07213. 2019.
  * [21] Kutuzov, A., Kuzmenko, E.: WebVectors: A Toolkit for Building Web Interfaces for Vector Semantic Models. CCIS 661, 155–161 (2017). 10.1007/978-3-319-52920-2 15
  * [22] Laposhina, A.N., Veselovskaya, T.S., Lebedeva, M.U. et al.: Automated text readability assessment for Russian second language learners. In: Komp. Lingv. i Intel. Tehn., pp. 396–406 (2018).
  * [23] Loper, E., Bird, S.: NLTK: the natural language toolkit, arXiv preprint cs/0205028 (2002).
  * [24] Loukachevitch, N., Levchik, A.: Creating a General Russian Sentiment Lexicon. In: Proc. of LREC-2016, pp. 1171-1176 (2016).
  * [25] Lyashevskaya, O. N., Sharov S. A.: Frequency Dictionary of the Modern Russian Language (based on the materials of the National Corps of the Russian Language). Azbukovnik, Moscow (2009).
  * [26] Mukherjee, P., Leroy, G., Kauchak D.: Using Lexical Chains to Identify Text Difficulty: A Corpus Statistics and Classification Study. J. of Biomed. and Health Informatics 23(5), 2164–2173 (2019). 10.1109/jbhi.2018.2885465
  * [27] Oborneva, I. V.: Automated estimation of complexity of educational texts on the basis of statistical parameters. Pedagogy Cand. Diss. Moscow (2006).
  * [28] Paszke, A., Gross, S., Massa, F. et al.: Pytorch: An imperative style, high-performance deep learning library. In: Adv. in neural information processing systems, pp. 8026-8037 (2019).
  * [29] Pedregosa, F., Varoquaux, G., Gramfort, A. et al.: Scikit-learn: Machine learning in Python. The J. of Machine Learning research 12, 2825–2830 (2011).
  * [30] Piasecki, S., Malekpour, S.: Morality and religion as factors in age rating computer and video games: ESRA, the Iranian games age rating system. Online-Heidelberg Journal of Religions on the Internet, 11.
  * [31] Potdar, K., Pardawala, T. S., Pai, C. D.: A comparative study of categorical variable encoding techniques for neural network classifiers. International journal of computer applications 1754, 7–9 (2017). 10.5120/ijca2017915495
  * [32] Russian National Corpus, https://ruscorpora.ru/new/en/index.html. Last accessed 23 Jul 2020.
  * [33] Shafaei, M., Samghabadi, N.S., Kar, S., Solorio, T.: Age Suitability Rating: Predicting the MPAA Rating Based on Movie Dialogues. In: Proceedings of The 12th Language Resources and Evaluation Conference, pp. 1327-1335 (2020).
  * [34] Sharoff, S.: Meaning as use: exploitation of aligned corpora for the contrastive study of lexical semantics. In: Proc. of LREC02, pp. 447–452. Las Palmas, Spain (2002).
  * [35] Schicchi, D., Pilato, G., Bosco, G.L.: Deep Neural Attention-Based Model for the Evaluation of Italian Sentences Complexity. In: 2020 IEEE 14th ICSC, pp. 253-256. 10.1109/icsc.2020.00053
  * [36] Schwarm, S. E., Ostendorf, M. Reading level assessment using support vector machines and statistical language models. In: Proc. of ACL’05, pp. 523–530 (2005). 10.3115/1219840.1219905
  * [37] Solnyshkina, M., Ivanov, V., Solovyev, V.: Readability Formula for Russian Texts: A Modified Version. In: Proc. of MICAI 2018, part II, pp. 142–145. Springer, Cham (2018). 10.1007/978-3-030-04497-8 11
  * [38] Solovyev, V., Solnyshkina, M., Ivanov, V. et al.: Prediction of reading difficulty in Russian academic texts. J. of Int. & Fuzzy Systems 36(5), 4553–4563 (2019). 10.3233/jifs-179007
  * [39] Sung, Y.T., Chen, J.L., Cha, J.H. et al.: Constructing and validating readability models: the method of integrating multilevel linguistic features with machine learning. Behavior research methods 47(2), 340–354 (2015). 10.3758/s13428-014-0459-x
  * [40] Templin, M. C.: Certain language skills in children; their development and interrelationships. Univ. of Minnesota Press, Minneapolis (1957).
  * [41] Text readability rating, http://readability.io/. Last accessed 23 Jul 2020.
  * [42] Tomina, Yu. A.: Objective assessment of the language difficulty of texts (description, narration, reasoning, proof). Pedagogy Cand. Diss. Moscow (1985).
  * [43] Wolf, T., Debut, L., Sanh, V. et al.: HuggingFace’s Transformers: State-of-the-art Natural Language Processing. ArXiv, arXiv-1910 (2019).
  * [44] Zagoruiko, N.G.: Applied methods of data and knowledge analysis. Izd-vo IM SO RAN, Novosibirsk (1999).


# Are there any extragalactic high speed dark matter particles in the Solar
neighborhood?

Isabel Santos-Santos  Nassim Bozorgnia  Azadeh Fattahi  and Julio F. Navarro

###### Abstract

We use the APOSTLE suite of cosmological hydrodynamical simulations of the
Local Group to examine the high speed tail of the local dark matter velocity
distribution in simulated Milky Way analogues. The velocity distribution in
the Solar neighborhood is well approximated by a generalized Maxwellian
distribution sharply truncated at a well-defined maximum “escape” speed. The
truncated generalized Maxwellian distribution accurately models the local dark
matter velocity distribution of all our Milky Way analogues, with no evidence
for any separate extragalactic high-speed components. The local maximum speed
is well approximated by the terminal velocity expected for particles able to
reach the Solar neighborhood in a Hubble time from the farthest confines of
the Local Group. This timing constraint means that the local dark matter
velocity distribution is unlikely to contain any high-speed particles
contributed by the Virgo Supercluster “envelope”, as argued in recent works.
Particles in the Solar neighborhood with speeds close to the local maximum
speed can reach well outside the virial radius of the Galaxy, and, in that
sense, belong to the Local Group envelope posited in earlier work. The local
manifestation of such envelope is thus not a distinct high-speed component,
but rather simply the high-speed tail of the truncated Maxwellian
distribution.

##  1 Introduction

Evidence for the gravitational interaction of dark matter (DM) with ordinary
matter exists on various astronomical scales, from galaxies to the largest
scale structures in the Universe (see e.g. refs. [1, 2, 3, 4]). If DM consists
of new elementary particles, we can search for its non-gravitational
interactions with ordinary matter in the laboratory. Direct detection
experiments search for DM particles with mass in the range of $\sim$ [keV –
TeV] in underground detectors, by measuring the small recoil energy of a
target nucleus or electron after interacting with a DM particle passing
through the detector. Axion searches, on the other hand, directly search for
the conversion of light DM ($\sim[10^{-6}-10^{-3}]$ eV mass) into photons in
the detector. A proper understanding of the DM density and kinematics in the
Solar neighborhood is needed in order to interpret the results of these
searches.

The Standard Halo Model (SHM) [5] is the most commonly adopted model for the
DM distribution in the Milky Way (MW). It assumes that DM is distributed in an
isotropic and isothermal sphere around the center of the Galaxy, and that the
“local” DM velocity distribution (i.e., in the Solar neighbourhood) may be
modelled as a Maxwell-Boltzmann distribution with a peak speed close to the
local circular speed. Recent high resolution simulations of galaxy formation,
which include both DM and baryons, find that a Maxwellian distribution fits
well the local DM velocity distribution of simulated MW-like halos [6, 7, 8,
9, 10, 11, 12]. However, simulations show large halo-to-halo scatter, which
leads to a large range of best fit peak speeds for the local Maxwellian
velocity distribution.

In addition, simulations studying the effect of large satellite galaxies such
as the Large Magellanic Cloud (LMC) on the MW show that the LMC significantly
impacts the high velocity tail of the local DM velocity distribution, shifting
it to higher speeds [13, 14, 15]. This has strong consequences for the
interpretation of results from direct detection experiments, leading to
considerable shifts of the direct detection exclusion limits towards smaller
DM-nucleon or DM-electron cross sections and DM masses. In general, the impact
of the LMC on direct detection limits is more significant for low mass DM,
where the experiments probe the high speed tail of the local DM velocity
distribution.

Prior work has suggested a further effect that may impact the highest
velocities expected for DM particles in the Solar neighborhood, namely the
possibility that “extragalactic” DM particles from the “envelopes” of the
Local Group (LG) or the Virgo Supercluster may add a population of particles
in a narrow range of velocities (about the MW local escape speed for the LG
envelope, roughly $600$ km/s, or exceeding $1000$ km/s for the Virgo
Supercluster) [16, 17, 18, 19]. Such high speed DM particles in the Solar
neighborhood can significantly impact the sensitivity of DM direct detection
experiments [18, 19].

If such component exists, it should be present in cosmological simulations of
the formation of the LG, which capture the highly dynamic nature of the
assembly of the halos of the MW and M31 galaxies and their time-evolving
configuration in the proper cosmological context. In this work, we use the
APOSTLE hydrodynamical simulations [20, 21] of LG volumes to investigate
whether a high speed extragalactic DM particle component is found in the Solar
neighborhood of simulated MW-like halos. We consider MW analogues in both
dark-matter-only (DMO) and hydrodynamical (Hydro) runs of the APOSTLE
simulations at three resolution levels, and study in detail their local DM
velocity distributions, with particular emphasis on their high-speed tail. In
particular, we focus on what sets the maximum speed of DM particles in the
Solar neighborhood, and on whether it may be affected by the presence of
extraneous matter of extragalactic origin.

This paper is structured as follows. In section 2 we provide details on the
APOSTLE simulations. In section 3 we present the local DM velocity
distribution extracted from the simulated halos, model its high speed tail,
and present a physical interpretation for the origin of the maximum speed
expected in each halo. We discuss and summarize our main results in section 4.

##  2 Simulations

The APOSTLE project [20, 21] consists of a suite of cosmological
N-body/hydrodynamical simulations of LG-like volumes selected from a
$100^{3}\,$Mpc3 periodic box and resimulated at different resolutions using
the zoom-in technique [22]. The LG-like volumes were identified based on the
kinematic properties of the MW-Andromeda pair and of the surrounding Hubble
flow. More precisely, each volume includes a pair of $\sim 10^{12}\,{\rm
M}_{\odot}$ halos with separation of $[600-1000]$ kpc, and relative radial and
tangential velocities of $[-250,0]$ km/s and $[0-100]$ km/s, respectively. The
selection also attempts to match, as much as possible, the observed Hubble
flow out to $3$ Mpc. Like the observed LG, the halo pairs are relatively
isolated, with no other halos of comparable or greater mass within 2.5 Mpc
from their barycentre. Details of the LG volume selection process can be found
in ref. [21].

The simulations were performed using a modified version of the Tree-PM SPH
code, P-Gadget3 [23], with the EAGLE galaxy formation model [24]111using the
“Reference” model in the language of ref. [24]. which has been shown to
reproduce a wide range of galaxy properties such as stellar mass function,
mass-size relation and observed rotation curves [24, 25, 26]. The simulation
volumes have also been run in DMO mode, where all the simulation particles are
treated as gravitating collisionless fluid. We use both the Hydro and DMO runs
in this study.

Halos and galaxies in the simulations were identified using the Friends-of-
Friends (FoF) algorithm with the linking length of 0.2 times the mean
interparticle separation, followed by the SUBFIND algorithm to find self-bound
halos and subhalos [23]. Cosmological parameters of the simulations are
adopted from WMAP-7 [27] and are those of a flat $\Lambda$CDM universe with
$h=0.704$, $\Omega_{0}=0.272$, $\Omega_{\Lambda}=0.728$, $\Omega_{\rm
bar}=0.0455$, $\sigma_{8}=0.81$.

The APOSTLE suite includes 12 LG-like volumes, but for this study we only use
4 volumes that have been simulated at all three resolution levels, both as DMO
and Hydro. This enables us to check the sensitivity of our results to
numerical limitations. DM and gas particles have masses of $\sim 5\times
10^{4}\,{\rm M}_{\odot}$ and $\sim 10^{4}\,{\rm M}_{\odot}$, respectively, for
the highest resolution (HR) runs. The particle mass increases by $\sim
10\times$ and $\sim 100\times$ in the case of medium resolution (MR) and low
resolution (LR), respectively, relative to HR222These resolution levels are
labelled as L1/L2/L3 for LR/MR/HR, respectively, in ref. [21].. The
gravitational softening length is $134$, $307$, and $711$ kpc for HR, MR, and
LR runs, respectively.

Amongst the 8 primary halos in the four APOSTLE volumes studied here (AP-1,
AP-4, AP-6, and AP-11), 4 of them have a satellite as massive as the LMC,
although they may not match the current position and velocity of the LMC [28,
13]. There has been no attempt in APOSTLE to match other observational
constraints on larger scales, like the presence of nearby groups such as M81
and Cen A; the Virgo cluster; or the Supergalactic Plane. As we discuss below,
we do not believe that this omission should have a substantial impact on our
conclusions.

##  3 Analysis and Results

###  3.1 Local dark matter velocity distribution in APOSTLE halos

We select at random one of the eight primary halos in the four APOSTLE volumes
to illustrate our analysis of the local DM velocity distribution. This halo is
the most massive of the primary pair in APOSTLE AP-6. The $z=0$ circular
velocity profiles, $v_{c}(r)=\sqrt{GM(<r)/r}$, of this halo are shown in
figure 1 for the DMO runs at three resolution levels (HR in thick solid, MR in
dashed, and LR in dotted line types). Assuming spherical symmetry, circular
velocities only depend on the total enclosed mass, $M(<r)$, within radius $r$,
measured from the halo centre, which is identified with the position of the
gravitational potential minimum.

Figure 1: Circular velocity profile of a DMO APOSTLE halo shown for the three
different resolution levels. HR, MR, and LR are shown in thick solid, dashed,
and dotted line types, respectively. Note the excellent numerical convergence
between the three different runs. The LR line becomes thinner at radii where
the difference with HR is greater than 10%. The gray vertical band indicates
the range of “Solar neighborhood” radii considered when extracting the “local”
DM velocity distributions. The circular velocity at that radius shows
excellent convergence for all three resolutions. The vertical dotted lines
indicate the virial radius, $r_{200}$. Downward pointing arrows mark $2\times$
and $8\times r_{200}$.

Note the excellent agreement between the LR, MR, and HR runs, except for the
innermost regions of the lowest resolution run, which shows a small deficit of
mass due to the limited number of particles and the use of a gravitational
softening. Outside $0.6$ kpc (for MR) and $3.3$ kpc (for LR) (marked by the
transition from thin to thick linetype in figure 1) all profiles agree very
well. The gray vertical shaded area indicates the position of the “Solar
neighborhood” in the simulations, defined as a spherical shell with an inner
and outer radius of 9 kpc and 11 kpc, respectively. The vertical dotted line
indicates the virial333Virial quantities are computed at a radius encompassing
a mean density contrast of $200$ times the critical density of the Universe,
and are identified by a “200” subscript. radius, $r_{200}$, of the halo.
Downward pointing arrows indicate radii equal to $2\times$ or $8\times$ the
virial radius, for future reference. All three simulations show good
convergence for the circular velocity in the Solar neighborhood. The number of
particles in the shell vary from $2,214$ for the LR halo to $24,946$ for MR,
and $37,560,165$ for HR.

Figure 2: Speed distribution of DM particles in the Solar neighborhood of the
same DMO halo shown in figure 1. The Solar neighborhood is defined as a
spherical shell at $r=10\pm 1$ kpc. The full range of DM speeds is shown on
the left panel; the right panel zooms in on the high speed tail of the
distribution. The red, cyan, and green histograms show the distributions for
the HR, MR, and LR levels, respectively, with $1\sigma$ Poisson error bars.
The maximum DM particle speed, $v_{\rm lmx}$, reached in each case is quoted
in the legend. The solid and dashed black lines show the generalized
Maxwellian (eq. (3.1)) and the _truncated_ generalized Maxwellian (eq. (3.2))
models that best fit the peak and shape of the local DM speed distribution for
the HR level, with parameters as listed on the right panel legend.

Figure 2 shows the speed distribution of DM particles in the Solar
neighborhood for the same halo shown in the previous figure. The left panel
shows the full range of DM particle speeds, whereas the right panel zooms in
on the high speed tail of the distribution. Velocities are measured in the
halo-centric rest frame (computed by SUBFIND as the center-of-mass velocities
of the halo DM particles), and the velocity modulus distribution (i.e. speed
distribution), $f(|\bf v|)$, is normalized to one, such that $\int f(|{\bf
v}|)\,dv=1$. Histograms of different colors correspond to different resolution
levels in the simulations. The error bars represent the $1\sigma$ Poisson
uncertainties on the data points.

There is clearly excellent agreement between the three different resolution
runs; the best-fit velocity distribution peaks at $190$ km/s and decays
rapidly at high speeds, vanishing beyond a well-defined maximum speed of
roughly $\sim 600$ km/s. More precisely, the local maximum DM particle speed
is $v_{\rm lmx}=554$, $577$, and $594$ km/s for the LR, MR, and HR resolution
levels, respectively.

###  3.2 The truncated generalized Maxwellian distribution

As shown in ref. [6], a generalized Maxwellian distribution,

$f(|{\bf v}|)\propto|{\bf v}|^{2}\exp\left[-(|{\bf
v}|/v_{0})^{2\alpha}\right],$ (3.1)

with free parameters $v_{0}$ and $\alpha$, reproduces the local DM speed
distribution of simulated MW analogues quite well. Notice that $\alpha=1$
recovers the standard Maxwellian form. The solid black curve in the left panel
of figure 2 shows that a generalized Maxwellian distribution with $v_{0}=124$
km/s and $\alpha=0.6$ fits the bulk of the local DM speed distributions from
the simulations quite well, except for the high-speed tail, where the
simulated distributions show a much sharper truncation than expected from the
generalized Maxwellian model.

One important point to note is that there is no hint of a separate
“extragalactic” component contributing an excess of high-speed DM particles;
if anything, it is the opposite: there are a lot fewer high-speed particles in
the simulated distribution than expected from the tail of the generalized
Maxwellian. Although we have shown this in detail for only one APOSTLE halo in
figure 2, this is a general result that applies to all the $8$ APOSTLE halos
we have analyzed, run either in DMO or Hydro mode.

The differences between the simulated local DM speed distributions and the
generalized Maxwellian fit at high-speeds can be better appreciated in the
right panel of figure 2: there is a clearly defined maximum speed of about
$600$ km/s for DM particles in this region, regardless of resolution. Accurate
fits to the local DM speed distribution, therefore, require that the
generalized Maxwellian model be modified to include a maximum speed. The
dashed black line in figure 2 corresponds to eq. (3.1) multiplied by a
_correction factor_ ,

$g(\lvert\textrm{\bf{v}}\rvert)=\textrm{tanh}\left(\eta\left[1-\frac{|v|}{v_{\rm
term}}\right]^{2}\right).$ (3.2)

Here $v_{\rm term}$ corresponds to the maximum (“terminal”) velocity beyond
which the distribution vanishes, and $\eta$ characterizes the range of speeds
over which the correction to eq. (3.1) needs to be applied. For the case shown
in figure 2, $v_{\rm term}=604$ km/s. The simulated distributions start
deviating from the generalized Maxwellian at $\sim 400$ km/s (see the right
panel of figure 2), which implies $\eta=30$. Below $400$ km/s the correction
factor (eq. (3.2)) becomes negligible.

Figure 3: Top: the median speed of DM particles in the tail of the local
velocity distribution, shown as a function of the fraction of mass in the
tail, for the same halo as in figure 2. As in that figure, HR is shown in red,
MR in cyan, and LR in green. The total number of DM particles in the Solar
neighborhood shell ($r=10\pm 1$ kpc) is quoted in the legend for each
resolution level. The solid and dashed black lines show the expected median
speeds for the generalized Maxwellian and for the truncated generalized
Maxwellian fits shown in figure 2, respectively. Bottom left: Best fit
generalized Maxwellian (solid) and truncated generalized Maxwellian (dashed)
distributions shown in figure 2. Downward pointing arrows indicate the
location of the tails of the distribution for various mass fractions in the
tail, as labelled. Bottom right: a zoomed-in version of the bottom left panel
highlighting the high-speed tail.

###  3.3 Local maximum speed and numerical resolution

Is this maximum speed just the result of limited numerical resolution, or is
it a physically meaningful local characteristic of the halo? In support of the
latter interpretation, we see from the legends in figure 2 that maximum speed
estimates ($v_{\rm lmx}$) increase only slightly with resolution, by roughly
$10\%$ as the number of particles is increased by a factor of $100$. This is
about what is expected for a velocity distribution which is actually truncated
at around $\sim 600$ km/s and sampled with different particle numbers.

We show this in the top panel of figure 3, where we plot the median speed of
particles in the high-speed tail of the local DM velocity distribution as a
function of the fraction of the mass in the tail, decreasing from left to
right. For the generalized Maxwellian fit to the data in figure 2, this median
speed is expected to increase without limit as the mass fraction in the tail
decreases. From the solid black curves in the top and bottom panels of figure
3 we see that the median speed is expected to increase from $444$ to $668$ to
$863$ km/s as the high-speed mass fraction goes from $10\%$ to $1\%$ to
$0.1\%$ of the total.

This is clearly not what is seen in the simulations; the median speed of high-
speed particles grow only from $397$ to $507$ to $555$ km/s for the same mass
fractions as above, following the increase expected for the truncated
Maxwellian distribution shown by the dashed lines in figures 2 and 3. Notice
that we use median speeds here rather than the maximum speed to reduce
numerical noise in the estimate. The good agreement between the three colored
curves obtained from the simulation data at the three different resolutions
and the truncated Maxwellian distribution, thus, strongly suggests that a
maximum speed of $\sim 600$ km/s is a true physical property of this APOSTLE
halo in the Solar neighborhood.

###  3.4 The origin of the local maximum speed

The next question that we would like to address is what sets this maximum
speed. For an equilibrium isolated system of finite mass, the maximum speed
could be safely identified with the local “escape speed”; i.e., the minimum
speed needed for a particle to reach infinity: higher-speed particles would
simply leave the system to never return. Unfortunately, DM halos in a
cosmological context are neither isolated nor fully in equilibrium, but rather
constantly evolving through accretion and mergers. In addition, total masses
in cosmology are not well defined, as enclosed masses increase without bound
with increasing radius, confounding the concept of escape speed.

Maximum speeds in non-linear objects formed by accretion (also referred to as
“secondary infall” in the literature) result from the depth of the potential
well as well as from the finite age of the Universe, which restricts the
maximum radius from where a particle could have been accreted by the present
time [29]. Particles coming from even further away could in principle reach
higher speeds, but they have not yet had time to reach a particular location.
The finite age of the Universe then imprints a maximum turnaround radius (and
thus a maximum speed) for particles at a given radius.

We can use this idea to provide a simple estimate of the maximum speeds
expected at any radius of a halo by using its present-day mass distribution
and the age of the Universe. For simplicity, the estimate assumes that the
mass distribution does not evolve in time and is spherically symmetric; this
is clearly a rather crude approximation which, however, works quite well at
reproducing the simulation results.

To illustrate the model, figure 4 shows the orbit of a test particle on a
radial orbit that reaches the centre of the halo at the present day after
evolving for $13.8$ Gyr. The acceleration on the particle is computed from the
sphericalized mass distribution of the halo and its surroundings at $z=0$;
i.e., for the same $M(<r)$ profile as in the simulation. Two cases are shown;
one (dashed curves) where the particle was initially at $r=0$ and moving
radially outwards with enough speed to reach a turnaround radius of $\sim 850$
kpc before returning to the centre after $13.8$ Gyr. The second case (solid
green curves) contemplates the case of a particle that starts at rest at $t=0$
from the maximum radius ($r\approx 1.6$ Mpc) from which such particle may
reach the centre by the present time. As may be seen from the inset in the
bottom panel of figure 4, the terminal velocity of the particle in case 2 is
only slightly higher than in case 1, indicating that much of the acceleration
occurs well inside $800$ kpc, the smaller of the two turnaround radii.

Figure 4: Radial distance (top panel) and radial velocity (bottom panel) of a
test particle as a function of time, orbiting within the gravitational
potential of the same DMO halo as in figure 1. Two radial orbits are shown,
both arriving at the halo centre for the first time after $13.8$ Gyr (the age
of the Universe). The orbits assume an unevolving spherical Newtonian
potential with the same $M(<r)$ profile as the halo and its surroundings at
$z=0$. The dashed curve corresponds to an orbit that is initially at $r=0$,
the solid line corresponds to an orbit which starts at rest at $t=0$. The
circles show the results at redshift $z=0$ in each case. The inset panel zooms
in on the “terminal velocity” of both orbits at $z=0$ (solid circles).

We shall hereafter adopt case 2 to define “terminal velocities” ($v_{\rm
term}$) at any radius of the halo, although we emphasize that adopting case 1
would give similar results. The same model, applied to the Solar neighborhood
(i.e., $r=10$ kpc, rather than the halo centre, as in the example shown in
figure 4), yields terminal velocities of order $600$ km/s, in excellent
agreement with the local maximum speeds obtained in the simulations and with
the asymptotic limit shown by the dashed curve in the bottom panel of figure
3.

Actually, the model reproduces quite well the maximum speeds obtained at all
radii, as may be seen in figure 5. Here the red curves show the maximum speeds
as a function of radius in the halo for the LR (dotted), MR (dashed) and HR
(solid) runs. The thick green curve shows the result of the terminal velocity
model presented above, which seems to capture well the radial trend of the
maximum speed at all radii. Only near the centre does the terminal velocities
slightly overestimate the simulation results, but even there the difference
between the predicted maximum speed and that obtained in the HR simulation is
less than $3\%$.

Figure 5: Maximum speeds as a function of radius for the same halo selected
for analysis in the previous figures. The grey vertical band identifies the
location of the Solar neighborhood, as adopted in this study. Red curves
indicate the maximum speeds for the HR (solid), MR (dashed), and LR (dotted)
runs, computed in concentric spherical shells. Maximum speeds increase
slightly with increasing resolution, and approach the terminal velocity curve
(shown in solid green) computed using the model illustrated in figure 4. The
terminal velocities agree quite well with the Newtonian “escape speeds”
computed assuming that the $z=0$ mass distribution of the halo is truncated at
$2\,r_{200}$. Such escape speeds are ill defined in cosmology, and depend on
the truncation radius assumed: blue curves show, for example, escape speeds
for a truncation radius of $8\,r_{200}$.

For comparison, we also show in figure 5 the radial dependence of the escape
speed from the halo computed by assuming that the gravitational potential of
the halo (which we assume is Newtonian) is given by the halo’s mass
distribution at $z=0$. This distribution needs to be truncated at some radius
for the escape speed to converge: we show two cases, truncated at either two
(cyan) or eight (blue) virial radii. Interestingly, the terminal velocities
are nearly identical to the escape speeds computed using a truncation radius
of $2\,r_{200}$. Note that, unlike the escape speed, the terminal velocity
model does not require the choice of an arbitrary truncation radius: it only
depends on the mass distribution in and around the halo, and on the age of the
Universe. In practice, as seen in figure 4, only the matter well inside $\sim
800$ kpc, or inside two to three virial radii seem to really matter.

###  3.5 Halo-to-halo scatter

Finally, we address how well the terminal velocity model proposed above
reproduces the results for other halos in different APOSTLE volumes. Each
volume has primary halos of different masses, surrounded by large scale mass
distributions that may also vary significantly from one LG realization to
another, so this is a rather demanding test of the model.

We compare in figure 6 the maximum speeds of DM particles in the Solar
neighborhood (i.e., at $r=10$ kpc) with the results of the terminal velocity
model introduced above. Results are shown for all eight APOSTLE primary halos
in 4 different volumes, run at LR, MR, and HR resolutions, both for DMO and
Hydro runs. As may be seen in this figure, the terminal velocity model does a
remarkable job at accounting for the maximum speeds seen in the HR
simulations, be them DMO or Hydro. As discussed in figure 3, LR and MR runs
give systematically (but slightly) lower estimates of the maximum speed
because their lower resolutions do not allow them to sample the highest
velocities of the distribution.

Figure 6: Maximum speed of DM particles in the Solar neighborhood (i.e., at
$r=10\pm 1$ kpc) versus the terminal velocity, $v_{\rm term}(r=10~{}{\rm
kpc})$, for the eight APOSTLE halos in the DMO (circles) and Hydro (diamonds)
simulations at HR (red), MR (cyan), and LR (green) resolution levels. Maximum
speeds increase slightly but systematically with increasing resolution, and
converge to the $1$:$1$ line (in dotted black) at the highest resolution
(shown in red). Note the excellent agreement between the terminal velocity
model and the local maximum speed for all HR halos at the Solar neighborhood.

Figure 6 also shows that maximum DM speeds in the Solar neighborhood may vary
greatly, even for halos which, like those in APOSTLE, were chosen to be
reasonably realistic realizations of MW analogues. The range spans
$\sim[350-600]$ km/s, which is much greater than the spread in circular
velocities at the same radius, which vary only between $115$ and $159$ km/s
for DMO runs, and between $129$ and $177$ km/s in Hydro runs. This may be seen
in figure 7, where we show the circular velocity profiles of all eight HR
APOSTLE halos in the DMO (left panel) and Hydro (right panel) runs. The curves
are colored by their corresponding local DM maximum speed (i.e. the values on
the $y$-axis in figure 6). Maximum circular velocities are marked with small
crosses in the figure.

The correlation between the local maximum DM speed and circular velocity is
thus rather poor, as seen by the black symbols in figure 8. This implies that
it could be quite difficult to estimate reliably the maximum DM speed expected
for the Solar neighborhood using only the circular velocity at the Solar
circle as a guide.

The main reason for this is that terminal velocities depend on the total mass
of the system out to at least twice the virial radius (and beyond), whereas
circular velocities depend only on enclosed masses. Cold DM halos follow
Navarro-Frenk-White (NFW) density profiles [30, 31], and, for NFW models, the
mass at fixed radii in the inner regions does not scale linearly with virial
mass [32]. For example, a $10^{12}\,{\rm M}_{\odot}$ NFW halo of average
concentration encloses roughly $\sim 4\times 10^{10}\,{\rm M}_{\odot}$ within
the inner $10$ kpc. Varying the virial mass by a factor of two above and below
that value leads to variations of less than $\sim 25\%$ in the DM mass
enclosed within that radius, so even small variations in circular velocity may
be accompanied by large changes in the terminal velocity at the same radius.

On the other hand, the correlation between local maximum speeds, $v_{\rm
lmx}$, and the maximum circular velocity, $V_{\rm max}$, of the system is
excellent, as may be seen from the cyan symbols in figure 8. This is mainly
because the maximum circular velocities probe better the outskirts of the
halos, which are responsible for accelerating infalling particles to the
highest speeds. A simple relation, where the maximum speed in the Solar
neighborhood is given by $v_{\rm lmx}\sim 3.75\,V_{\rm max}^{0.94}$ fits well
the simulation results. The dashed cyan line in figure 8 shows this relation;
for reference, the central escape speed of an NFW halo is $3.04\,V_{\rm max}$.
Tight constraints on the MW’s $V_{\rm max}$ are therefore needed in order to
make accurate predictions about the maximum speed expected for DM particles in
the Solar neighborhood.

Figure 7: Circular velocity profiles, $v_{c}(r)$, for each of the eight HR
APOSTLE halos in the DMO (left) and Hydro (right) simulations. Profiles are
colored according to the maximum speed of DM particles in the Solar
neighborhood, $v_{\rm lmx}$ (see color bar). Small crosses indicate $V_{\rm
max}$, the maximum circular velocity of each system. Figure 8: Local maximum
DM speeds for the DMO (circles) and Hydro (diamonds) halos in the HR run
compared with the local circular velocity (i.e., $v_{c}$ at $10$ kpc; black
symbols) and with the maximum circular velocity of each system (cyan symbols).
Local maximum DM speeds are poorly correlated with the local circular
velocity, but are tightly related to the maximum circular velocity, as shown
by the simple proportionality curve, $v_{\rm lmx}=3.75\times V_{\rm
max}^{0.94}$, shown in dashed blue.

##  4 Summary and Conclusions

We have used the APOSTLE suite of cosmological hydrodynamical simulations [20,
21] to study the high-speed tail of the DM velocity distribution in the Solar
neighborhood of MW-like halos. Our main goals are to characterize the tail of
the distribution, to determine the origin of the maximum speeds achieved by DM
particles, and to search for the presence of distinct components of high-speed
DM particles of extragalactic origin hypothesized in earlier work [16, 17, 18,
19].

In particular, we studied eight APOSTLE MW analogues at three resolution
levels (LR, MR, and HR), run in dark-matter-only mode, and in Hydro mode;
i.e., including gas and galaxy formation. We find that at the location of the
Sun (identified with $r=10$ kpc in the simulations) the distribution of DM
velocities is well approximated by a generalized Maxwellian distribution (eq.
(3.1)), sharply truncated at a clearly defined maximum speed (eq. (3.2)).

There is no sign of a separate “extragalactic” component, or an excess of
high-speed particles over the Maxwellian distribution. If anything, the
opposite is true, for all simulations show that the local DM velocity
distributions do not extend to arbitrarily large velocities but are truncated
at a clearly defined, radially dependent maximum speed. At the Solar
neighborhood, this maximum speed varies from halo to halo, but is insensitive
to the numerical resolution of the simulations.

Maximum speeds may be reliably estimated using the $z=0$ mass distribution of
the halo (assumed non-evolving) and its surroundings to compute the “terminal
velocity” of a particle in a radial orbit able to reach from rest, in a Hubble
time, the Solar circle for the first time. This simple model reproduces
remarkably well the results of all APOSTLE DMO and Hydro halos not only at
$r=10$ kpc, but at all radii.

The terminal velocities agree fairly well with the Newtonian “escape speed”
estimated when truncating the $z=0$ halo mass distribution at twice the virial
radius. This means that particles near the tail of the distribution have
apocentres well beyond the confines of the halo. In that sense, the local
high-speed tail does originate in the LG “envelope” hypothesized in earlier
work [17]. However, these envelope particles do not constitute a discernible
separate component of high-speed particles, but are rather just the local
representatives of the far outskirts of the halo, which spills beyond its
nominal virial boundary.

We see also no sign of any “extra-high-speed” DM particles contributed by
nearby groups or clusters, a result not unexpected given the success of the
terminal velocity model presented above: particles from outside the LG have
not yet had time to reach the Solar neighborhood in the finite age of the
Universe. We conclude that it is highly unlikely that the Solar neighborhood
contains any $\sim 1000$ km/s extragalactic particles contributed by the Virgo
Supercluster, as speculated in previous work [17].

The halo-to-halo scatter in the maximum DM speed at the Solar neighborhood is
much larger than expected from the scatter in the circular velocity at the
same radius. This is because the terminal velocities depend on the mass
distribution extending well outside the virial radius, whereas circular
velocities depend only on the enclosed mass within the Solar circle. This
implies that it is not straightforward to estimate local maximum speeds based
solely on the local circular velocity.

The simulations suggest, on the other hand, a strong correlation between the
maximum circular velocity of the system and the local maximum DM speed. This
is because the maximum halo circular velocity is a better probe of the mass
distribution on larger scales, which is needed to make an accurate prediction
of the local maximum DM speed. If the maximum circular velocity of the MW is
$220$ km/s [33], the simulation results suggest a local maximum DM speed of
$3.75\times 220^{0.94}=597$ km/s.

This estimate of the local maximum DM speed depends critically on the maximum
circular velocity assumed, as well as on the assumption that the Solar
neighborhood is close to equilibrium. Strong transients, such as the recent
arrival of the Magellanic Clouds into the Milky Way halo, could result in
large fluctuations in the halo, causing significant shifts in the high-speed
tail of the local DM velocity distribution that have been recently studied in
both idealized and cosmological simulations [13, 14, 15].

## Acknowledgments

ISS acknowledges support from the European Research Council (ERC) through
Advanced Investigator grant to C.S. Frenk, DMIDAS (GA 786910). NB acknowledges
the support of the Canada Research Chairs Program, the Natural Sciences and
Engineering Research Council of Canada (NSERC), funding reference number
RGPIN-2020-07138, and the NSERC Discovery Launch Supplement, DGECR-2020-00231.
AF is supported by a UKRI Future Leaders Fellowship (grant no MR/T042362/1).
This work used the DiRAC Memory Intensive system at Durham University,
operated by ICC on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk).
This equipment was funded by BIS National E-infrastructure capital grant
ST/K00042X/1, STFC capital grant ST/H008519/1, and STFC DiRAC Operations grant
ST/K003267/1 and Durham University. DiRAC is part of the National
E-Infrastructure.

## References

  * [1] G. Bertone (ed.), Particle Dark Matter: Observations, Models and Searches, vol. 2012.  Cambridge Univ. Press, Cambridge, 2010.
  * [2] G. Jungman, M. Kamionkowski, and K. Griest, Supersymmetric dark matter, Phys.Rept. 267 (1996) 195–373, [hep-ph/9506380].
  * [3] L. Bergström, Non-baryonic dark matter: observational evidence and detection methods, Rep. Prog. Phys. 63 (apr, 2000) 793–841, [hep-ph/0002126].
  * [4] G. Bertone, D. Hooper, and J. Silk, Particle dark matter: Evidence, candidates and constraints, Phys. Rept. 405 (2005) 279–390, [hep-ph/0404175].
  * [5] A. K. Drukier, K. Freese, and D. N. Spergel, Detecting Cold Dark Matter Candidates, Phys. Rev. D33 (1986) 3495–3508.
  * [6] N. Bozorgnia, F. Calore, M. Schaller, M. Lovell, G. Bertone, C. S. Frenk, R. A. Crain, J. F. Navarro, J. Schaye, and T. Theuns, Simulated Milky Way analogues: implications for dark matter direct searches, JCAP 1605 (2016), no. 05 024, [arXiv:1601.04707].
  * [7] C. Kelso, C. Savage, M. Valluri, K. Freese, G. S. Stinson, and J. Bailin, The impact of baryons on the direct detection of dark matter, JCAP 1608 (2016) 071, [arXiv:1601.04725].
  * [8] J. D. Sloane, M. R. Buckley, A. M. Brooks, and F. Governato, Assessing Astrophysical Uncertainties in Direct Detection with Galaxy Simulations, Astrophys. J. 831 (2016) 93, [arXiv:1601.05402].
  * [9] N. Bozorgnia and G. Bertone, Implications of hydrodynamical simulations for the interpretation of direct dark matter searches, Int. J. Mod. Phys. A32 (2017), no. 21 1730016, [arXiv:1705.05853].
  * [10] N. Bozorgnia, A. Fattahi, C. S. Frenk, A. Cheek, D. G. Cerdeno, F. A. Gómez, R. J. Grand, and F. Marinacci, The dark matter component of the Gaia radially anisotropic substructure, JCAP 07 (2020) 036, [arXiv:1910.07536].
  * [11] R. Poole-McKenzie, A. S. Font, B. Boxer, I. G. McCarthy, S. Burdin, S. G. Stafford, and S. T. Brown, Informing dark matter direct detection limits with the ARTEMIS simulations, JCAP 11 (2020) 016, [arXiv:2006.15159].
  * [12] E. Rahimi, E. Vienneau, N. Bozorgnia, and A. Robertson, The local dark matter distribution in self-interacting dark matter halos, JCAP 02 (2023) 040, [arXiv:2210.06498].
  * [13] A. Smith-Orlik, N. Ronaghi, N. Bozorgnia, M. Cautun, A. Fattahi, G. Besla, C. S. Frenk, N. Garavito-Camargo, F. A. Gómez, R. J. J. Grand, F. Marinacci, and A. H. G. Peter, The impact of the Large Magellanic Cloud on dark matter direct detection signals, arXiv e-prints (Feb., 2023) arXiv:2302.04281, [arXiv:2302.04281].
  * [14] G. Besla, A. Peter, and N. Garavito-Camargo, The highest-speed local dark matter particles come from the Large Magellanic Cloud, JCAP 1911 (2019), no. 11 013, [arXiv:1909.04140].
  * [15] K. Donaldson, M. S. Petersen, and J. Peñarrubia, Effects on the local dark matter distribution due to the Large Magellanic Cloud, Mon. Not. Roy. Astron. Soc. L513 (2022), no. 1 46, [arXiv:2111.15440].
  * [16] K. Freese, P. Gondolo, and L. Stodolsky, On the direct detection of extragalactic WIMPs, Phys. Rev. D 64 (2001) 123502, [astro-ph/0106480].
  * [17] A. N. Baushev, Extragalactic dark matter and direct detection experiments, Astrophys. J. 771 (2013) 117, [arXiv:1208.0392].
  * [18] G. Herrera and A. Ibarra, Direct detection of non-galactic light dark matter, Phys. Lett. B 820 (2021) 136551, [arXiv:2104.04445].
  * [19] G. Herrera, A. Ibarra, and S. Shirai, Enhanced prospects for direct detection of inelastic dark matter from a non-galactic diffuse component, JCAP 04 (2023) 026, [arXiv:2301.00870].
  * [20] T. Sawala, C. S. Frenk, A. Fattahi, J. F. Navarro, R. G. Bower, R. A. Crain, C. Dalla Vecchia, M. Furlong, A. Jenkins, I. G. McCarthy, Y. Qu, M. Schaller, J. Schaye, and T. Theuns, Bent by baryons: the low-mass galaxy-halo relation, MNRAS 448 (Apr., 2015) 2941–2947, [arXiv:1404.3724].
  * [21] A. Fattahi, J. F. Navarro, T. Sawala, C. S. Frenk, K. A. Oman, R. A. Crain, M. Furlong, M. Schaller, J. Schaye, T. Theuns, and A. Jenkins, The APOSTLE project: Local Group kinematic mass constraints and simulation candidate selection, Monthly Notices of the Royal Astronomical Society 457 (Mar., 2016) 844–856, [arXiv:1507.03643].
  * [22] A. Jenkins, A new way of setting the phases for cosmological multiscale Gaussian initial conditions, MNRAS 434 (Sept., 2013) 2094–2120, [arXiv:1306.5968].
  * [23] V. Springel, N. Yoshida, and S. D. M. White, GADGET: a code for collisionless and gasdynamical cosmological simulations, New Astronomy 6 (Apr., 2001) 79–117, [astro-ph/0003162].
  * [24] J. Schaye, R. A. Crain, R. G. Bower, M. Furlong, M. Schaller, T. Theuns, C. Dalla Vecchia, C. S. Frenk, I. G. McCarthy, J. C. Helly, A. Jenkins, Y. M. Rosas-Guevara, S. D. M. White, M. Baes, , and . more authors, The EAGLE project: simulating the evolution and assembly of galaxies and their environments, MNRAS 446 (2015) 521–554, [arXiv:1407.7040].
  * [25] R. A. Crain, J. Schaye, R. G. Bower, M. Furlong, M. Schaller, T. Theuns, C. Dalla Vecchia, C. S. Frenk, I. G. McCarthy, J. C. Helly, A. Jenkins, Y. M. Rosas-Guevara, S. D. M. White, and J. W. Trayford, The EAGLE simulations of galaxy formation: calibration of subgrid physics and model variations, MNRAS 450 (June, 2015) 1937–1961, [arXiv:1501.01311].
  * [26] M. Schaller, C. S. Frenk, R. G. Bower, T. Theuns, A. Jenkins, J. Schaye, R. A. Crain, M. Furlong, C. D. Vecchia, and I. McCarthy, Baryon effects on the internal structure of $\Lambda$CDM haloes in the EAGLE simulations, Mon. Not. Roy. Astron. Soc. 451 (2015), no. 2 1247–1267, [arXiv:1409.8617].
  * [27] E. Komatsu, K. M. Smith, J. Dunkley, C. L. Bennett, B. Gold, G. Hinshaw, N. Jarosik, D. Larson, M. R. Nolta, L. Page, D. N. Spergel, M. Halpern, R. S. Hill, A. Kogut, M. Limon, S. S. Meyer, N. Odegard, G. S. Tucker, J. L. Weiland, E. Wollack, and E. L. Wright, Seven-year Wilkinson Microwave Anisotropy Probe (WMAP) Observations: Cosmological Interpretation, The Astrophysical Journal Supplement Series 192 (Feb., 2011) 18, [arXiv:1001.4538].
  * [28] I. M. E. Santos-Santos, A. Fattahi, L. V. Sales, and J. F. Navarro, Magellanic satellites in $\Lambda$CDM cosmological hydrodynamical simulations of the Local Group, Mon. Not. Roy. Astron. Soc. 504 (2021), no. 3 4551–4567, [arXiv:2011.13500].
  * [29] E. Bertschinger, Cosmological self-similar shock waves and galaxy formation, The Astrophysical Journal 268 (May, 1983) 17–29.
  * [30] J. F. Navarro, C. S. Frenk, and S. D. M. White, The Structure of Cold Dark Matter Halos, The Astrophysical Journal 462 (May, 1996) 563, [astro-ph/9508025].
  * [31] J. F. Navarro, C. S. Frenk, and S. D. M. White, A Universal Density Profile from Hierarchical Clustering, The Astrophysical Journal 490 (Dec., 1997) 493–508, [astro-ph/9611107].
  * [32] I. Ferrero, J. F. Navarro, M. G. Abadi, J. A. Benavides, and D. Mast, A unified scenario for the origin of spiral and elliptical galaxy structural scaling laws, Astronomy and Astrophysics 648 (Apr., 2021) A124, [arXiv:2009.03916].
  * [33] M. Cautun, A. Benitez-Llambay, A. J. Deason, C. S. Frenk, A. Fattahi, F. A. Gómez, R. J. Grand, K. A. Oman, J. F. Navarro, and C. M. Simpson, The Milky Way total mass profile as inferred from Gaia DR2, Mon. Not. Roy. Astron. Soc. 494 (2020), no. 3 4291–4313, [arXiv:1911.04557].