# When BERT Meets Quantum Temporal Convolution Learning
for Text Classification in Heterogeneous Computing

###### Abstract

The rapid development of quantum computing has demonstrated many unique
characteristics of quantum advantages, such as richer feature representation
and more secured protection on model parameters. This work proposes a vertical
federated learning architecture based on variational quantum circuits to
demonstrate the competitive performance of a quantum-enhanced pre-trained BERT
model for text classification. In particular, our proposed hybrid classical-
quantum model consists of a novel random quantum temporal convolution (QTC)
learning framework replacing some layers in the BERT-based decoder. Our
experiments on intent classification show that our proposed BERT-QTC model
attains competitive experimental results in the Snips and ATIS spoken language
datasets. Particularly, the BERT-QTC boosts the performance of the existing
quantum circuit-based language model in two text classification datasets by
1.57% and 1.52% relative improvements. Furthermore, BERT-QTC can be feasibly
deployed on both existing commercial-accessible quantum computation hardware
and CPU-based interface for ensuring data isolation.

Index Terms—  Quantum machine learning, temporal convolution, text
classification, spoken language understanding, and heterogeneous computing

##  1 Introduction

Text classification (e.g., intent detection) from human spoken utterances [1,
2] is an essential element of a spoken language understanding [3] (SLU)
system. Learning intent information often involves various applications, such
as on-device voice assistants [4] and airplane travel information systems [5].
Recently, Bidirectional Encoder Representations from Transformers [6] (BERT)
has caused a stir in the lexical modeling community by providing competitive
results in intent classification as a common SLU application that we torched
upon in this work. However, recent works [7, 8] have raised new concerns about
data leakage from large-scale language models such as BERT, which can involve
sensitive information like personal identification. New regulation
requirements [9] (e.g., GDPR) on data protection, privacy [10], and security
can further motivate advanced investigations in designing distributed
algorithms on heterogeneous computing devices.

In recent two years, cloud-accessible quantum devices (e.g., IBM-Q) have shown
unique characteristics and empirical advantages in many applications [11, 12],
such as model compression, parameter isolation, and encryption. A noisy-
intermediate-scale-quantum (NISQ) device is a main hardware category that
empowers the quantum advantages by using only a few number of qubits (5 to
100). When the limitation of quantum resource is concerned, variational
quantum circuit (VQC) has been studied on the design of quantum machine
learning like quantum support vector machine, where the VQC algorithm requires
no quantum error correction as a matching working on NISQ devices. Projecting
classical data into high-dimensional quantum feature space has been proven its
quantum advantages [13, 14] (e.g., a quadratic speed-up [13] and better
presentation power [11]) on some classifications tasks. For example, Yang _et
al._ [11] demonstrates the quantum advantages on a speech command recognition
system by applying the VQC algorithm to extract acoustic features for
representation learning with a vertical federated learning architecture. The
proposed vertical federated learning could benefit from quantum
advantages111Quantum advantage means that a programmable quantum device [15]
can solve a specific problem that is intractable on classical computers., such
as parameter isolation with random circuit learning, for the acoustic modeling
task. However, designing an end-to-end neural model on SLU tasks is still an
open problem with a potential impact on heterogeneous computing devices.

Motivated by the very recent success [11] in the quantum circuit-based
acoustic modeling, we put forth a new gradient-based training to include pre-
trained BERT models and enhance its data protection by leveraging upon VQC
learning in this work. To strengthen the parameter isolation on BERT, our work
proposes a novel hybrid classical-quantum system with quantum temporal
convolution (QTC) learning, which is composed of a pre-trained BERT model (on
a classical computer) and a new type of one-dimensional random circuit [16,
11] (on a quantum device). In comparison with the benchmark pre-trained BERT
model, our proposed architecture can maintain competitive experimental
performance considering the SLU benchmark solution running in homogeneous
machines. Notably, in this work, compared with the classical deep neural
network (DNN) models, the use of QTC lowers the model complexity. The VQC
design of QTC is to use the quantum circuit as a temporal convolutional
features projector on the text embeddings that requires only a few qubits
(e.g., 4 to 9) to set up our quantum platform.

Fig. 1: BERT-QTC for vertical federated [11] text classification.

As shown in Figure 1, the design of such a BERT-QTC model is in a regime of
hybrid classical-quantum architecture, where the word embedding can be offline
attained by utilizing DNN models on classical computers before going through
quantum devices. Table 1 shows an overview of classical, quantum [17], and
heterogeneous approaches for SLU. To the best of our knowledge, this is the
first work to investigate a BERT-based hybrid classical-quantum model for text
classification (e.g., intent detection) task with competitive empirical
results on quantum circuit learning on the NISQ device.

Table 1: An overview of different natural language processing (NLP)
approaches: neural language model (NLM), logic programming (LP), and its
quantum variants. Our work belongs to hybrid CQ to access quantum advantages
(QA).

Approach Input Model Output Challenges Classical bits NLM & LP bits data
protection Quantum qubits Quantum LP [18] qubits hardware limits hybrid CQ
bits NLM + VQC bits model design

##  2 Related Work

###  2.1 Quantum Algorithms for Language Processing

Recent years have witnessed a rapid development [19] of quantum machine
learning on near-term quantum devices (5 to 200 qubits). Encoding classical
data into high-dimensional quantum features has proven rigorous quantum
advantages on classification tasks with recent theoretical justification [13,
14]. In particular, VQC learning [20] is utilized as parametric models to
build a QNN model. The research of VQC is of significance in itself because
they constitute the setting of quantum computational experiments on NISQ
devices. Although business-accessible quantum services (e.g., IBM-Q and Amazon
Braket) are still actively being developed, several works have attempted to
incorporate quantum algorithms for language processing tasks. For example,
[17] Liu _et al._ consider a novel classifier as the physical evolution
process and described the process with quantum mechanical equations. In [21,
22], the theory of quantum probability is proposed to build a general language
model for information retrieval.

Moreover, Blacoe _et al._ [23] leverage upon quantum superposition and
entanglement to investigate the potential of quantum theory as a practical
framework to capture the lexical meaning and model semantic processes such as
word similarity and association. However, those works aim at the use of
quantum concepts to build quantum language models in theory, and they are
different from the recent quantum circuit architectures featured with small-
to-none quantum error correction in the NISQ devices.

Recently, Meichanetzidis _et al._[24] first investigated computational logic
relationships for encoding graphs into circuit learning on NISQ computers.
Coecke _et al._ [18] investigate theoretical foundation and advantages of
encryption on using the VQC learning on linguistic modeling. In particular,
their quantum circuits are employed as a pipeline to map semantic and grammar
diagrams and pave the way to near-term applications of quantum computation
devices. Lorenz _et al._ [25] introduced a concept of compositional semantic
information for text classification tasks on the quantum hardware by simply
using a mapping from lexical segments to a quantum circuit. Nevertheless,
beyond grammar modeling, how to design a specific quantum model is still a
wide-open topic with potential impacts on “large-scale” data (e.g., spoken
language processing and applications).

###  2.2 BERT and Heterogeneous Computing

Devlin _et al._ [26] firstly put forth the architecture of BERT, which mainly
applies the masking and bidirectional training of Transformer equipped with an
attention model to language modeling. BERTs [27] can be used for a wide
variety of language tasks since different SLU tasks only require a few more
layers to the core model for fine-tuning. When the training time and data
protection become significant concerns for BERTs, heterogeneous computing
architectures, such as distributed training and federated learning [28],
provide a new perspective to ensure data isolation as protection and improve
training efficiency.

In our previous work [11], we propose a new vertical federated learning
architecture that combines new convolution-encoded randomized VQC and
recurrent neural network (RNN) (denoted as QCNN-RNN) for acoustic modeling
preserving quantum advantages. The QCNN-RNN system first uses a quantum
convolution embedding layer to encode features in high-dimensional Hilbert
space, then encode quantum features in a classical format for the
representation learning with RNN. Different with the quantum acoustic modeling
task, we want to advance a new VQC design working on the SLU task considering
the benchmark pre-trained language model (e.g., BERT) in this work.

##  3 Method

###  3.1 BERT for Text Embedding

Fig. 2: The BERT-QCC comprises with (a) BERT and (b) quantum temporal
convolutional (QTC) learning with random circuits.

We first use BERT as a pre-trained language model to extract hidden
representation features from text sequences as shown in Figure 2 (a). The BERT
based language models consist of a multi-layer bidirectional Transformer-based
encoder, where the input corresponds to a combination of WordPiece embeddings,
positional embeddings, and the segment embedding. Symbolically, given an input
token sequence $\boldsymbol{x}=\left(x_{1},\ldots,x_{T}\right)$ the BERT model
outputs are computed by:

$\mathbf{H}=\mathbf{BERT}\left(x_{1},\ldots,x_{T}\right),$ (1)

where $\mathbf{H}=\left(\boldsymbol{h}_{1},\ldots,\boldsymbol{h}_{T}\right)$
and $\boldsymbol{h}_{i}$ is contextual semantic representation embedding of
each token in Eq. (1).

Neural Network Intent Decoder. We select a hidden state of the first special
token for intent classification after a dropout layer to avoid over-fitting.
Conventionally, the hidden state is encoded by a neural network function
($\boldsymbol{f}(\mathbf{H};\theta_{\boldsymbol{f}})$) for extracting
sentence-level semantic representation. Finally, an unknown intent can be
predicted as:

$y^{i}=\operatorname{softmax}\left(\mathbf{W}^{i}\boldsymbol{f}(\mathbf{H};\theta_{\boldsymbol{f}})+\boldsymbol{b}^{i}\right),$
(2)

where $y^{i}$ is a predicted intent label for a test sequence. The model is
end-to-end trained by minimizing the cross-entropy loss between predicted text
labels and their corresponding ground-truth labels.

###  3.2 Quantum Temporal Convolution with Random Circuit

The expressive power of quantum convolution has been recently studied in
computer vision [16] and speech processing [11, 29] by using quantum space
encoding to project classical data into rich representations of features for
classification. Deploying a quantum circuit with randomized assigned (non-
trainable) parameters could beat both trainable and non-trainable
convolutional filters in the previous study [11]. We further consider a new
hybrid classical-quantum decoder design, which replace the perception network
($\boldsymbol{f}$) in Eq. (2) with a randomized quantum circuit decoder as
shown in Fig. 2 (b). We use a max length of 50 for pre-trained BERT models.

Variational Quantum Circuit Decoder. We adopt the concept from variational
quantum circuit learning and define a latent encoding function
$\hat{\boldsymbol{f}}$, where it contains angle encoding, parameterized
rotation, and quantum-to-classical state decoding. We use the VQC as shown in
Fig. 3, made by strong entanglement [30] circuit. First, the Bert embeddings
$[\mathbf{h}_{1},\cdots,\mathbf{h}_{T}]$ are feeding into a sliding window for
the temporal quantum convolutional filter (e.g., coordinated with time or
order of input), where we take window size $I=4$ for example as one desired
qubits size shown in Fig. 4. We then extract latent embedding
$[h_{1},h_{2},h_{3},h_{4}]$ as inputs to the VQC layer.

@C=0.8em @R=0.8em —0⟩ & R_x(h_1) R(α_1, β_1, γ_1) 1

—0⟩

R_x(h_2) R(α_2, β_2, γ_2)  1

—0⟩

R_x(h_3) R(α_3, β_3, γ_3)  1

—0⟩

R_x(h_4) R(α_4, β_4, γ_4)  -3 1347.7em–

Fig. 3: Deployed quantum circuit. The VQC component contains three major
parts: encoding, learning and quantum measurement. Here we use the angle
encoding scheme to encode the input values $h_{1}\cdots h_{4}$ (latent
embedding; taking four qubits, for example) by treating them as rotation
angles along $x$-axis. The learning part uses general unitary rotation $R$.
There are three parameters $\alpha$, $\beta$ and $\gamma$ in each $R$. The
controlled-NOT (CNOT) gates are used to entangle quantum states from each
qubit. The final quantum measurement part will output the Pauli-$Z$
expectation values of each qubit.  Fig. 4: Computing framework for quantum
temporal convolution (QTC) with random variational circuit from input text
embeddings.

The latent features ($H$) are the input of a quantum circuit layer,
$\mathbf{Q}$, that learns and encodes latent sequence:

$\hat{\boldsymbol{f}}=\mathbf{Q}(H;\mathbf{e},\theta_{\mathbf{q}},\mathbf{d})$
(3)

In Eq. (3), the approximation process of a quantum circuit block,
$\mathbf{Q}$, depends on the encoding initialization $\mathbf{e}$, the quantum
circuit parameters, $\theta_{\mathbf{q}}$, and the decoding measurement
$\mathbf{d}$. When we stack more VQC blocks, we expect that a gradient
vanishing issue would occur and degrade the representation learned from a
hybrid QNN model whose parameters are trained by joint DNN-VQC gradient
updates.

Table 2: Mathematical notation for VQC learning. Symbol | Meaning
---|---
$|0\rangle$ |  quantum state $0$
Rx |  rotation gate along $x$-axis
$\theta_{\mathbf{q}}$=$\\{\alpha,\beta,\gamma\\}$ | random parameters
CNOT | controlled-NOT gate

In this VQC, the elements $h_{1}\cdots h_{4}$ of the input vectors are used as
rotational angles for single-qubit quantum rotation gate $R_{x}$ on each qubit
to encode the classical values. After encoded into a quantum state, the state
is processed through a _learnable_ or _optimizable_ layer (in grouped box).
The parameters labeled with $\alpha$, $\beta$, and $\gamma$ are for the
optimization. The CNOT gates are applied to entangle the qubits. In the final
step, the quantum measurement procedure is carried out. In this part, we
retrieve the Pauli-$Z$ expectation values for each qubit. The obtained result
is a $4$-dimensional vector which can be further processed by other classical
or quantum routines.

Figure 3 summarize characterises of the variational quantum circuit, where
there are $4$ quantum channels and each channel is mutually entangled with
each other by applying the CNOT gate. Given the unitary matrix $\mathbb{U}$
representing the multiplication of all quantum gates and taking $B$ as
observable matrix, we attain that:

$\hat{\boldsymbol{f}}(h;\boldsymbol{\theta_{\text{random init}}})=\langle
h|U^{\dagger}(\boldsymbol{\theta}_{i})BU(\boldsymbol{\theta}_{i})|h\rangle.$
(4)

The hyper-parameters ($\boldsymbol{\theta_{\text{random init}}}$) for rotation
gates in Eq. (4) are randomly initialized and not considered to be updated
during the training phase, which is aimed to provide parameters protection in
heterogeneous computing architectures (or simulation API) against model
inversion attacks [31] and parameters leakages [32]. As shown in Fig. 2, the
temporal quantum convolution could be applied with $n$ filters to map features
before computing with a global max-pooling layer. As the first attempt to
construct quantum temporal convolution, we select a filters number from
$\\{1,2,3,4\\}$ under the 9 qubits requirement for common commercial quantum
hardware (e.g., IBMQ).

To evaluate the effectiveness of proposed QTC architectures, we select three
additional “random” encoder baselines (similar to QTC) with text embeddings
for text classification: (1) BERT with a random temporal convolutional network
(TCN); (2) word2vec [33] with random TCN; (3) word2vec with random QTC, all
followed the same filter numbers and global max-pooling after the deployed
random encoders.

We will study how random QTC learning benefits BERT-based heterogeneous
architecture in the experimental section. From a system-level perspective, the
proposed QTC learning reduces the risk of parameter leakage [32, 34, 35] from
inference-time attackers, and it tackles data isolation issues with its
architecture-wise advantages [36] on encryption [37] and the feature of
without accessing the data directly [38, 11].

##  4 Experiment

###  4.1 Quantum Computing Hardware

We use PennyLane as an open-source and reproducible environment, which uses
differential programming to calculate gradient to update gradient circuits. A
hardware back-end of PennyLane could integrate from CPU-simulation, QPU
(supported by Amazon Braket), and TPU (supported by Google Floq). Since fine-
tuning time of BERTs is often extensive, we first use a CPU-simulated VQC
environment to train the proposed BERT-QTC model. Then, we evaluate our hybrid
classical-quantum models on Amazon Braket and Floq hardware and report test
accuracy by 10-fold cross-validation. Referring to the established work [11]
on VQC learning, we consider a NISQ device with 4 to 16 qubits, which allows
us to perform up to 9 class predictions. Furthermore, we refine our dataset
according to the hardware settings, and clarify that NISQ’s size constraints
limit current applications of proposed BERT-QTC working toward word-level
slot-filling tasks (with 120 to 72 classes). We believe that the development
of commercial-accessible NISQ could gradually resolve this challenge. For
vertical learning setting, we use secure multiparty computation (MPC) protocol
[39] between local instance and quantum device during the virtualization.

###  4.2 Dataset and Setup

Snips Dataset: To evaluate the proposed framework in complex spoken language
corpora, we select Snips [4] dataset. Snips dataset is a collection of spoken
utterances from a personal voice assistant, whose domains cover speech
command, music searching, and smart home request. The training set includes
13,084 utterances, and the test set includes 700 utterances. We use another
700 utterances as the development set. There are 7 types of intent classes,
where the number of samples for each intent is nearly the same.

Fig. 5: ATIS7 Dataset for BERT-QTC experiments.

ATIS Dataset: The ATIS (Airline Travel Information System) dataset [5] is
widely used in NLU research containing spoken corpus, including intent
classification and slot filling labels for flights reservation. Since the
language scripts from ATIS are collected from human speakers, ATIS is selected
to highlight the data sensitivity under GDPR policy [9]. To compare with Snips
as a slightly unbalanced classification study, we further select a subset of
ATIS containing top-7 classes (denoted as ATIS7) of original ATIS. As shown in
Fig. 5, the deployed subset include 90.93% of the original data, where
training, development, and testing sets contain 4072, 545, 556 utterances,
respectively. As a remark, we have conducted text classification experiments
with a full ATIS dataset, but the results have a large variance ($\pm$ 7.32%)
on both random TCN and QRC encoders due to unbalanced intent labels in a long-
tail distribution [40].

BERT Model: We use an open source BERT [26] model to connect the QTC encoder.
The pre-trained models of BERT are frozen in the heterogeneous computing
process as a universal encoder to take sentences into embedding. We use the
official pretrained BERT (large) model from TensorFlow Hub with 1024-hidden
dimensions, 24 layers, and 12 attentions heads with 330M pre-trained
parameters.

###  4.3 BERT-QTC Performance

Since we are working on sentence-level text classification (intent as labels)
on the two deployed datasets, we first select two different text embeddings
methods from pre-trained word2vec and BERT embeddings to compared the proposed
random QTC encoder with a random TCN encoder. Furthermore, we aim to study
different hyper-parameters setup for TCN and QTC architecture. As a remark,
both TCN and QTC could be simulated with CPU-based environments or API, where
QTC is featured with better representation mapping with theoretical
justifications and an option running with quantum hardware to preserve full
quantum advantages.

As shown in Tab. 3 (Snips) and 4 (ATIS7), the BERT-QTC model in vertical
federated architecture performs the best average prediction accuracy
performance, which attains $\mathbf{96.62}$% in Snips and $\mathbf{96.98}$% in
ATIS7 compared with TCN based architectures and word2vec with QTC encoder.
Moreover, we investigate the convolution filter and kernel setups on BERTs-QTC
models from the second, fourth row in Tab. 3 and 4. BERT-QTC demonstrates more
significant improvement with its word2vec-QTC baseline when two filters have
been used. Interestingly, we find out that utilize QTC encoder shows a general
performance-boosting in the two deployed SLU datasets. Proposed BERT-QTC
federated models perform relative improvements of +1.57% in Snips and +1.52%
in ATIS7 compared with other heterogeneous computing ablations. The
statistical variances are $\leq 0.52\%$ in our experiments.

Table 3: Average accuracy on intent classification for Snips with a set of
different number (n) of convolutional filter and kernel size (k).

Embedding word2vec BERT (n,k) (1,4) (2,2) (2,3) (2,4) (1,4) (2,2) (2,3) (2,4)
TCN 82.02 83.37 82.90 83.15 95.48 95.23 95.12 95.27 QTC 83.32 83.94 83.61
$\mathbf{84.64}$ 96.41 96.42 $\mathbf{96.62}$ 96.44

Table 4: Average accuracy on intent classification for ATIS7 with a set of
different number (n) of convolutional filter and kernel size (k).

Embedding word2vec BERT (n,k) (1,4) (2,2) (2,3) (2,4) (1,4) (2,2) (2,3) (2,4)
TCN 80.09 80.22 80.91 82.34 95.18 95.03 94.95 95.23 QTC 81.42 82.49 83.82
$\mathbf{83.95}$ 96.69 96.92 96.32 $\mathbf{96.98}$

##  5 Conclusion

In this paper, we propose a novel hybrid classical-quantum architecture to
strengthen the BERT model with a quantum circuit decoder via quantum temporal
convolution (QTC) with random circuit learning. The proposed BERT-QTC models
show competitive results for text classification as one prompted finding for
standard SLU tasks. Moreover, our VQC encoders are capable of deploying on
both existing quantum hardware, and the simulator requires only a small amount
of qubits (4 to 8 qubits). The proposed QTC can enhance the data protection on
top of BERT models in the vertical federated learning setting.

## References

  * [1] A. Stolcke, K. Ries, N. Coccaro, E. Shriberg, R. Bates, D. Jurafsky, P. Taylor, R. Martin, C. V. Ess-Dykema, and M. Meteer, “Dialogue act modeling for automatic tagging and recognition of conversational speech,” _Computational linguistics_ , vol. 26, no. 3, pp. 339–373, 2000.
  * [2] D. Wang and G. J. Brown, _Computational auditory scene analysis: Principles, algorithms, and applications_. Wiley-IEEE press, 2006.
  * [3] Y.-Y. Wang, L. Deng, and A. Acero, “Spoken language understanding,” _IEEE Signal Processing Magazine_ , vol. 22, no. 5, pp. 16–31, 2005.
  * [4] A. Coucke, A. Saade, A. Ball, T. Bluche, A. Caulier, D. Leroy, C. Doumouro, T. Gisselbrecht, F. Caltagirone, T. Lavril _et al._ , “Snips voice platform: an embedded spoken language understanding system for private-by-design voice interfaces,” _arXiv preprint arXiv:1805.10190_ , 2018\.
  * [5] C. T. Hemphill, J. J. Godfrey, and G. R. Doddington, “The atis spoken language systems pilot corpus,” in _Speech and Natural Language: Proceedings of a Workshop Held at Hidden Valley, Pennsylvania, June 24-27, 1990_ , 1990.
  * [6] J. Devlin, M.-W. Chang, K. Lee, and K. Toutanova, “Bert: Pre-training of deep bidirectional transformers for language understanding,” in _Proceedings of the 2019 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, Volume 1 (Long and Short Papers)_ , 2019, pp. 4171–4186.
  * [7] N. Carlini, C. Liu, Ú. Erlingsson, J. Kos, and D. Song, “The secret sharer: Evaluating and testing unintended memorization in neural networks,” in _28th $\\{$USENIX$\\}$ Security Symposium ($\\{$USENIX$\\}$ Security 19)_, 2019, pp. 267–284.
  * [8] N. Carlini, F. Tramer, E. Wallace, M. Jagielski, A. Herbert-Voss, K. Lee, A. Roberts, T. Brown, D. Song, U. Erlingsson _et al._ , “Extracting training data from large language models,” _arXiv preprint arXiv:2012.07805_ , 2020.
  * [9] P. Voigt and A. Von dem Bussche, “The eu general data protection regulation (gdpr),” _A Practical Guide, 1st Ed., Cham: Springer International Publishing_ , 2017.
  * [10] C.-H. H. Yang, S. M. Siniscalchi, and C.-H. Lee, “Pate-aae: Incorporating adversarial autoencoder into private aggregation of teacher ensembles for spoken command classification,” _Proc. Interspeech 2021_ , pp. 881–885, 2021\.
  * [11] C.-H. H. Yang, J. Qi, S. Y.-C. Chen, P.-Y. Chen, S. M. Siniscalchi, X. Ma, and C.-H. Lee, “Decentralizing feature extraction with quantum convolutional neural network for automatic speech recognition,” in _IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2021, pp. 6523–6527.
  * [12] F. Arute, K. Arya, R. Babbush, D. Bacon, J. C. Bardin, R. Barends, R. Biswas, S. Boixo, F. G. Brandao, D. A. Buell _et al._ , “Quantum supremacy using a programmable superconducting processor,” _Nature_ , vol. 574, no. 7779, pp. 505–510, 2019.
  * [13] Y. Liu, S. Arunachalam, and K. Temme, “A rigorous and robust quantum speed-up in supervised machine learning,” _Nature Physics_ , pp. 1–5, 2021.
  * [14] H.-Y. Huang, M. Broughton, M. Mohseni, R. Babbush, S. Boixo, H. Neven, and J. R. McClean, “Power of data in quantum machine learning,” _Nature communications_ , vol. 12, no. 1, pp. 1–9, 2021.
  * [15] S. Boixo, S. V. Isakov, V. N. Smelyanskiy, R. Babbush, N. Ding, Z. Jiang, M. J. Bremner, J. M. Martinis, and H. Neven, “Characterizing quantum supremacy in near-term devices,” _Nature Physics_ , vol. 14, no. 6, pp. 595–600, 2018\.
  * [16] M. Henderson, S. Shakya, S. Pradhan, and T. Cook, “Quanvolutional neural networks: powering image recognition with quantum circuits,” _Quantum Machine Intelligence_ , vol. 2, no. 1, pp. 1–9, 2020.
  * [17] D. Liu, X. Yang, and M. Jiang, “A novel classifier based on quantum computation,” in _Proceedings of the 51st Annual Meeting of the Association for Computational Linguistics (Volume 2: Short Papers)_ , 2013, pp. 484–488.
  * [18] B. Coecke, G. de Felice, K. Meichanetzidis, and A. Toumi, “Foundations for near-term quantum natural language processing,” _arXiv preprint arXiv:2012.03755_ , 2020.
  * [19] V. Havlíček, A. D. Córcoles, K. Temme, A. W. Harrow, A. Kandala, J. M. Chow, and J. M. Gambetta, “Supervised learning with quantum-enhanced feature spaces,” _Nature_ , vol. 567, no. 7747, pp. 209–212, 2019.
  * [20] K. Mitarai, M. Negoro, M. Kitagawa, and K. Fujii, “Quantum circuit learning,” _Physical Review A_ , vol. 98, no. 3, p. 032309, 2018.
  * [21] A. Sordoni, J.-Y. Nie, and Y. Bengio, “Modeling term dependencies with quantum language models for ir,” in _Proceedings of the 36th international ACM SIGIR conference on Research and development in information retrieval_ , 2013, pp. 653–662.
  * [22] I. Basile and F. Tamburini, “Towards quantum language models,” in _Proceedings of the 2017 Conference on Empirical Methods in Natural Language Processing_ , 2017, pp. 1840–1849.
  * [23] W. Blacoe, E. Kashefi, and M. Lapata, “A quantum-theoretic approach to distributional semantics,” in _Proceedings of the 2013 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies_ , 2013, pp. 847–857.
  * [24] K. Meichanetzidis, S. Gogioso, G. De Felice, N. Chiappori, A. Toumi, and B. Coecke, “Quantum natural language processing on near-term quantum computers,” _arXiv preprint arXiv:2005.04147_ , 2020.
  * [25] R. Lorenz, A. Pearson, K. Meichanetzidis, D. Kartsaklis, and B. Coecke, “Qnlp in practice: Running compositional models of meaning on a quantum computer,” _arXiv preprint arXiv:2102.12846_ , 2021.
  * [26] J. Devlin, M.-W. Chang, K. Lee, and K. Toutanova, “Bert: Pre-training of deep bidirectional transformers for language understanding,” _arXiv preprint arXiv:1810.04805_ , 2018.
  * [27] Z. Lan, M. Chen, S. Goodman, K. Gimpel, P. Sharma, and R. Soricut, “Albert: A lite bert for self-supervised learning of language representations,” in _International Conference on Learning Representations_ , 2019.
  * [28] S. Y.-C. Chen and S. Yoo, “Federated quantum machine learning,” _arXiv preprint arXiv:2103.12010_ , 2021.
  * [29] J. Qi, C.-H. H. Yang, and P.-Y. Chen, “Qtn-vqc: An end-to-end learning framework for quantum neural networks,” _arXiv preprint arXiv:2110.03861_ , 2021.
  * [30] S. Y.-C. Chen, C.-H. H. Yang, J. Qi, P.-Y. Chen, X. Ma, and H.-S. Goan, “Variational quantum circuits for deep reinforcement learning,” _IEEE Access_ , vol. 8, pp. 141 007–141 024, 2020.
  * [31] M. Fredrikson, S. Jha, and T. Ristenpart, “Model inversion attacks that exploit confidence information and basic countermeasures,” in _Proceedings of the 22nd ACM SIGSAC conference on computer and communications security_ , 2015, pp. 1322–1333.
  * [32] A. Duc, S. Dziembowski, and S. Faust, “Unifying leakage models: From probing attacks to noisy leakage.” in _Annual International Conference on the Theory and Applications of Cryptographic Techniques_. Springer, 2014, pp. 423–440.
  * [33] T. Mikolov, K. Chen, G. Corrado, and J. Dean, “Efficient estimation of word representations in vector space,” _arXiv preprint arXiv:1301.3781_ , 2013\.
  * [34] D. Leroy, A. Coucke, T. Lavril, T. Gisselbrecht, and J. Dureau, “Federated learning for keyword spotting,” in _IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2019, pp. 6341–6345.
  * [35] M. Chen, A. T. Suresh, R. Mathews, A. Wong, C. Allauzen, F. Beaufays, and M. Riley, “Federated learning of n-gram language models,” _arXiv preprint arXiv:1910.03432_ , 2019.
  * [36] C. Dwork, V. Feldman, M. Hardt, T. Pitassi, O. Reingold, and A. Roth, “The reusable holdout: Preserving validity in adaptive data analysis,” _Science_ , vol. 349, no. 6248, pp. 636–638, 2015.
  * [37] A. C.-C. Yao, “Quantum circuit complexity,” in _Proceedings of 1993 IEEE 34th Annual Foundations of Computer Science_. IEEE, 1993, pp. 352–361.
  * [38] Q. Yang, Y. Liu, T. Chen, and Y. Tong, “Federated machine learning: Concept and applications,” _ACM Transactions on Intelligent Systems and Technology (TIST)_ , vol. 10, no. 2, pp. 1–19, 2019.
  * [39] R. Cramer, I. B. Damgård _et al._ , _Secure multiparty computation_. Cambridge University Press, 2015.
  * [40] C.-H. H. Yang, L. Liu, A. Gandhe, Y. Gu, A. Raju, D. Filimonov, and I. Bulyko, “Multi-task language modeling for improving speech recognition of rare words,” in _2021 IEEE Automatic Speech Recognition and Understanding Workshop (ASRU)_. IEEE, 2021, pp. 1087–1093.


11institutetext: Dipartimento di Fisica e Astronomia, Universitá di Bologna,
Via Gobetti 93/2, 40122, Bologna, Italy 22institutetext: Hamburger Sternwarte,
University of Hamburg, Gojenbergsweg 112, 21029 Hamburg, Germany
33institutetext: Istituto di Radioastronomia, INAF, Via Gobetti 101, 40122,
Bologna, Italy 44institutetext: School of Physics and Astronomy, University of
Minnesota Twin Cities Minneapolis, MN, USA 55institutetext: INAF, Osservatorio
di Astrofisica e Scienza dello Spazio, via Pietro Gobetti 93/3, 40129 Bologna,
Italy

#  Exploring the relation between turbulent velocity and density fluctuations
in the stratified intracluster medium

M. Simonte 1122 F. Vazza 112233 F. Brighenti 11 M. Brüggen 22 T. W. Jones 44
M. Angelinelli 1155

(Accepted ???. Received ???; in original form ???)

###### Abstract

Context. The dynamics of the intracluster medium (ICM) is affected by
turbulence driven by several processes, such as mergers, accretion and
feedback from active galactic nuclei.

Aims. X-ray surface brightness fluctuations have been used to constrain
turbulence in galaxy clusters. Here, we use simulations to further investigate
the relation between gas density and turbulent velocity fluctuations, with a
focus on the effect of the stratification of the ICM.

Methods. In this work, we studied the turbulence driven by hierarchical
accretion by analysing a sample of galaxy clusters simulated with the
cosmological code ENZO. We used a fixed scale filtering approach to
disentangle laminar from turbulent flows.

Results. In dynamically perturbed galaxy clusters, we found a relation between
the root mean square of density and velocity fluctuations, albeit with a
different slope than previously reported. The Richardson number is a parameter
that represents the ratio between turbulence and buoyancy, and we found that
this variable has a strong dependence on the filtering scale. However, we
could not detect any strong relation between the Richardson number and the
logarithmic density fluctuations, in contrast to results by recent and more
idealised simulations. In particular, we find a strong effect from radial
accretion, which appears to be the main driver for the gas fluctuations. The
ubiquitous radial bias in the dynamics of the ICM suggests that homogeneity
and isotropy are not always valid assumptions, even if the turbulent spectra
follow Kolmogorov’s scaling. Finally, we find that the slope of the velocity
and density spectra are independent of cluster-centric radii.

###### Key Words.:

Galaxies: clusters: intracluster medium – Galaxies: clusters: general –
X-rays: galaxies: clusters

##  1 Introduction

In the current paradigm of hierarchical structure formation, accretion of
smaller structures and mergers play a key role in the evolution of galaxies
and galaxy clusters. Gas that falls into dark matter halos is heated, mostly
via shocks or turbulent dissipation, to the virial temperature of a cluster
($10^{7}-10^{8}K$). This results in a fully ionised X-ray emitting
intracluster medium (ICM). During the process of cluster formation, mergers
and accretion generate hydrodynamic instabilities that create turbulence.
Observations and simulations have investigated the properties of turbulent
motions in the ICM. In this work, we extend previous theoretical efforts by
examining the competing roles of stratification and large-scale infall
motions.

There is observational evidence that supports the idea of a turbulent ICM.
Examples of observed features used to measure ICM turbulence include X-ray
surface brightness fluctuations, interpreted as density fluctuations induced
by turbulence, or pressure fluctuations obtained from X-ray maps (Schuecker et
al. 2004, Sanders & Fabian 2012). Another method to infer information about
turbulence is the study of emission lines in the X-ray band. Sanders et al.
(2010) placed limits on the turbulent broadening of the emission lines in
cool-core clusters. One of the main findings so far for the internal
kinematics of the hot gas was obtained observing the cool-core Perseus cluster
with the Hitomi satellite (Hitomi Collaboration et al. 2016). Those authors
found a 1D velocity dispersion $\sim 160\rm\leavevmode\nobreak\
km\leavevmode\nobreak\ s^{-1}$, on a scale of $L\sim 50$ kpc in the core of
the Perseus cluster. Similar subsonic turbulent velocities are also suggested
by studies of a relatively cold ICM in cool-core clusters (e.g. Olivares et
al. 2019; Russell et al. 2019). Indeed, if the cold gas originates by hot ICM
cooling, the two phases share the same kinematics (Gaspari et al. 2018).
Turbulence requires a large Reynolds number, which implies a low viscosity.
This is also indicated by the evidence of sharp features connected to Kelvin-
Helmholtz instabilities and cold fronts in the ICM (Markevitch & Vikhlinin
2007), which otherwise would be suppressed due to thermal conduction and
viscosity (Roediger et al. 2013, Wang & Markevitch 2018). Finally, it is
possible to exploit the thermal Sunyaev-Zeldovich effect to infer turbulent
pressure maps (Khatri & Gaspari 2016).

Turbulent motions can indirectly affect other observables of galaxy clusters
because they affect non-thermal components, such as magnetic fields and cosmic
rays. Combining X-ray and radio observations, Eckert et al. (2017) and
Bonafede et al. (2018) found that the brightness fluctuations correlate with
diffuse radio emission, suggesting that turbulence might be linked to the
acceleration of non-thermal particles. Rebusco et al. (2006) investigated the
effect of turbulent diffusion on the metal abundance profile. Moreover,
turbulence might play a role in the measurement of cluster masses, as it
creates additional kinetic pressure that affects the hydrostatic equilibrium
typically assumed to apply to the ICM. In fact, observational and numerical
work has confirmed that the ratio between turbulent and thermal energy
typically varies between $\sim 10-30\%$ (e.g. Angelinelli et al. 2020; Eckert
et al. 2019; Ettori et al. 2019; Lau et al. 2009). Since precision cosmology
depends on very accurate total cluster mass, turbulent pressure profiles need
to be well constrained. Finally, turbulence has been proposed as a source of
heat to prevent development of catastrophic cooling flows (Zhuravleva et al.
2014a). The most promising scenario there involves turbulent feedback induced
by active galactic nuclei. The central gaseous medium is heated by the
interaction with buoyantly rising bubbles that are created by jets launched
from the central black hole. While such feedback flows can be energetically
sufficient, the process by which the mechanical energy is transferred to the
ICM and thermalised is still debated (Brüggen & Kaiser 2002; De Young 2010;
Gaspari et al. 2013a; Gaspari et al. 2013b; Randall et al. 2015; Fabian et al.
2017).

Fully cosmological, or more idealised simulations have been used to
investigate turbulence in galaxy clusters (Dolag et al. 2005; Kim & Ryu 2005;
Wittor et al. 2017; Vazza et al. 2009, 2011b, 2017). Zhuravleva et al. (2014b)
and Gaspari et al. (2014) studied the relation between density and velocity
perturbations, Vazza et al. (2012) investigated the turbulent diffusion
coefficient and Angelinelli et al. (2020), Nelson et al. (2014), Lau et al.
(2009) constrained the non-thermal pressure contribution coming from turbulent
motions. One of the difficulties in determining ICM turbulent quantities is
the disentanglement of bulk and turbulent motions. Several algorithms have
been developed for this purpose. For example, Lau et al. (2009) estimated ICM
turbulent velocities as the residual with respect to the velocity field
averaged over spherical shells. Alternatively, one can estimate the turbulent
motions by interpolating the original 3-D velocity field in order to map the
local mean field and to identify turbulent velocity fluctuations on scales
smaller than the interpolation scale (Vazza et al. 2009, 2011b). As another
alternative, a multi-scale iterative filtering approach has been implemented
in, for example, Vazza et al. (2012), Angelinelli et al. (2020),
2021MNRAS.504..510V. In the present work, we apply a simpler approach, that is
by using a fixed scale filtering method, based on previous results by our
group.

Describing the turbulent nature of the ICM is not trivial. Astrophysical
turbulent flows are often described in terms of Kolmogorov’s theory
(Kolmogorov 1941), which assumes homogeneity and isotropy. However, the ICM is
a stratified plasma close to hydrostatic equilibrium, where buoyancy may well
change the character of the turbulence. In particular, the relation between
density and velocity fluctuations, often used to estimate the turbulence
strength in the ICM (e.g. Zhuravleva et al. 2014b) might vary in the presence
of a varying stratification (and forcing mechanism). The dimensionless
Richardson number, given by the square of the ratio of the turbulent eddy
turn-over timescale to the buoyancy oscillation timescale (or, equivalently,
by the ratio of the buoyancy to inertial forces), is a measure of the
importance of the stratification in the turbulence dynamics. It is defined as:

$Ri=\frac{N^{2}_{\rm BV}}{(v_{l}/l)^{2}},$ (1)

where $v_{l}$ is the characteristic turbulent velocity on a scale $l$, and
$N_{\rm BV}$ is the Brunt-Väisälä frequency:

$N_{\rm
BV}=\sqrt{-\frac{g}{\gamma}\frac{d\ln}{dr}\left(\frac{P}{\rho^{\gamma}}\right)}.$
(2)

Here, $g(r)$ is the gravitational acceleration, $P(r)$ and $\rho(r)$ are the
pressure and density, respectively, and $\gamma$ is the adiabatic index of the
ICM. Thus, $Ri\ll 1$, implies short turbulent eddy turn-over times compared to
buoyancy times, suggesting homogeneous, isotropic turbulence unaffected by
density stratification. On the other hand, when $Ri>1$ the buoyancy times are
shorter than the turbulent eddy times. When the buoyancy force becomes
dynamically important, it suppresses radial motions and it leads to
preferentially azimuthal turbulence (Brethouwer et al. 2007).

There are previous works investigating turbulence in a stratified ICM-like
medium (Zhuravleva et al. 2014b; Gaspari et al. 2014; Valdarnini 2019; Shi &
Zhang 2019; Mohapatra et al. 2020). However, they focused on idealised plasma
flows, or were limited to the innermost regions of galaxy clusters ($\leq
500\rm\leavevmode\nobreak\ kpc$). Previous observational and numerical work
investigated the possibility of inferring the statistics of turbulence in the
ICM from gas density and temperature perturbations, which can be recovered
from existing X-ray observations (Arévalo et al. 2016; Churazov et al. 2012;
Zhuravleva et al. 2018). However, due to the finite sensitivity of X-ray
telescopes, this technique is still limited by the large number of photons
required to sample large enough cluster scales, with sufficient spatial
detail. In order to further compare the results with observational data, past
analyses of simulations were restricted to the virial region of the cluster.
Moreover, idealised simulations are also bound to miss substructures and
filaments, which can bias the analysis compared to the statistics derived from
more realistic, dynamically formed clusters.

In this paper we explore the character of the turbulence in realistic,
stratified ICM generated in high-resolution, cosmological simulations of
galaxy clusters. We mainly focus on the relation between density and velocity
fluctuations, their dependence on the Richardson number (see Mohapatra et al.
2020), and the anisotropy of the turbulent velocity field.

This paper is structured as follows: in Section 2 we describe the simulated
cluster sample and the numerical method used to estimate the relevant
variables of a turbulent ICM; in Section 3 we discuss the outcomes of the
analysis, looking also for a comparison with previous work; in Section 4 we
summarise our findings, highlighting some observational implications.

##  2 Methods

###  2.1 The Itasca Simulated Cluster sample

We analysed a sample of galaxy clusters from the ’Itasca Clusters’ set of
simulations. Each simulation was carried out at high spatial resolution using
ENZO, a (parabolic) cosmological numerical code for magneto-hydrodynamics
(Bryan et al. 2014). A key feature of the code is its Adaptive Mesh Refinement
(AMR) capability, which enables a large spatial and temporal dynamical range.

Here we applied a pre-determined set of nested grids in order to apply a
constant, comoving grid resolution. Starting from the root grid, which covers
the entire volume with a coarse, uniform grid, we placed finer grids as soon
as interesting regions started to evolve. In the Itasca suite, we just forced
the code to refine the $100\%$ of the innermost zoom-in region, up to the
highest allowed AMR level. If the root grid spacing was $\Delta x$, then the
spacing of a refined patch at level $l$ was $\Delta x/r^{l}$, where r is the
integer refinement factor. The simulations presented here employed a purely
hydrodynamical fluid solver based on the Godunov Piecewise Parabolic Method
(PPM) scheme. (Colella & Woodward 1984).

The cosmological Itasca simulations assume a WMAP7 $\Lambda$CDM cosmology
(Komatsu et al. 2011), with $\Omega_{B}=0.0445$, $\Omega_{\rm DM}=0.2265$,
$\Omega_{\Lambda}=0.728$, Hubble parameter $h=0.702$, $\sigma_{8}=0.8$ and a
primordial index of $n=0.961$. All runs were non-radiative and did not include
feedback from star-forming regions or active galactic nuclei. Even though the
effects of non-gravitational heating are small on the $\gg$ 100 kpc scale
compared to the heat caused by mergers (Valdarnini 2019), the combination of
cooling and feedback leads to a higher number of substructures with different
temperatures and densities could create small-scale turbulent motions. In this
sense, such structures are absent in our simulations and their presence in the
observed galaxy clusters might alter gas density and temperature estimates.
Our simulations do not include magnetic fields even though the interplay
between turbulence and dynamo amplification is well-established in galaxy
clusters (e.g. 2012PhRvE..85b6303S; review_dynamo; va18mhd; dom19). However,
the magnetic pressure is small, considering that the typical plasma $\beta$
parameter (i.e. the ratio between thermal and magnetic pressure, $\beta$) is
very large, $\beta\sim 10^{2}-10^{3}$. This means that we do not expect any
relevant modification in the density or velocity structure of the ICM by the
(neglected) presence of magnetic fields, except on very small scales ($\leq
30\rm\leavevmode\nobreak\ kpc$) where the tension of amplified magnetic field
lines is significant (va18mhd). Nevertheless, the gas dynamics on these small
scales is not our primary concern and we proceed with a purely hydrodynamical
approach.

In brief, the Itasca simulations were run according to the following
procedure. Firstly, independent cosmological boxes were simulated in order to
select the most massive objects (clusters) in the volume. Initial conditions
were generated separately for each simulation at redshift $z=30$. The spatial
resolution, on the coarse grid, was $L_{0}=110\rm\leavevmode\nobreak\
h\leavevmode\nobreak\ kpc^{-1}\approx 157$ kpc and the dark matter (DM) mass
resolution was $m_{\rm DM}=8.96\cdot 10^{7}\rm\leavevmode\nobreak\ M_{\odot}$.
A second more refined grid was created, centred on the selected cluster
formation region, which covered the innermost $\approx$ 31 Mpc. In this case,
the comoving spatial resolution reached $55h$ kpc${}^{-1}\approx 78.4$ kpc and
the DM mass resolution was $m_{\rm DM}=1.12\cdot 10^{7}\rm\leavevmode\nobreak\
M_{\odot}$. Finally, inside the central $(L_{0}/10)^{3}\approx 6.3$ Mpc3
volume of each box, a fixed further refinement was enforced, meaning that
$100\%$ of cells were refined. This volume was large enough to include the
virial radii of most of the clusters, while the spatial resolution was
increased to $13.8\rm\leavevmode\nobreak\ h\leavevmode\nobreak\
kpc^{-1}\approx 20kpc$.

The full Itasca cluster sample consists of 20 galaxy clusters and has been
used to investigate the dynamic and thermal properties of cluster ICMs. In
this work, we used 8 Itasca clusters at $z=0$ in the mass range $0.54\cdot
10^{14}\rm\leavevmode\nobreak\ M_{\odot}<M_{200}<3.32\cdot
10^{14}\rm\leavevmode\nobreak\ M_{\odot}$. For further information on the
Itasca cluster sample, we refer to Vazza et al. (2017), Wittor et al. (2017),
Angelinelli et al. (2020).

###  2.2 Filtering method and numerical analysis

In order to disentangle bulk from turbulent motions, we used a tailored
filtering method. Vazza et al. (2012), Vazza et al. (2017) have shown that
galaxy clusters might present a distribution of turbulent injection scales,
since ICM turbulence can be driven by many processes (e.g AGN feedback or
accreted clumps). However, Vazza et al. (2018) utilising the Itasca clusters
found that typical scales from turbulence injected by accretion, are in the
range $\sim$ 200 - 400 kpc. For this reason, we choose $L=300$ kpc as a
filtering scale in the following analysis. Nevertheless, we tested the
robustness of our analysis by using other filtering scales.

Our filtering method estimates the turbulent fluctuation ($\delta v$,
$\delta\rho$) as follows:

$\delta v=v-\langle v\rangle_{L},$ (3)

$\delta\rho=\rho-\langle\rho\rangle_{L},$ (4)

Figure 1: Slice of the original density and velocity 3D distributions in one
cluster (left panels). The right panels show the turbulent density and
velocity fluctuations maps, after the application of the fixed-scale filtering
method.

where $v$ and $\rho$ are the local velocity and density and
$\langle\rho\rangle_{L}$ and $\langle v\rangle_{L}$ are the means estimated on
cubic boxes of (linear) size $L$, centred on the cell to filter. In the case
of the velocity, we use a density-weighted mean. Fig. 1 shows maps of density
and velocity fields, before and after we performed this algorithm in a cluster
of our sample. We then estimate the density contrast,
$\delta\rho/\langle\rho\rangle_{L}$ and the normalised velocity fluctuation,
$\delta v/c_{s}$, where $c_{s}$ is the sound speed at a given radius, obtained
from the (azimuthally averaged) temperature profile. We investigate the
relation between these two quantities to directly compare to the work by
Zhuravleva et al. (2014b).

In the following, we calculate the statistical quantities (for instance the
variance of the density contrast) in cubic boxes of linear size 600 kpc,
located at different radii and different directions. In this way we map the
whole simulation box to assess the spatial variation of the turbulence
properties. We also performed a similar analysis by dividing the main box in
spherical shells, with similar but somewhat less clear results. This can be
understood by the azimuthally asymmetric ICM dynamics, generated by accretion
events along the large scale cosmological filaments. Because of this we only
discuss the sub-boxes analysis in the rest of the paper.

####  2.2.1 Density and velocity fluctuations

We adopt as a measure for the typical density fluctuation the standard
deviation of the probability distribution function (PDF) of the density
contrast $\delta\rho/\langle\rho\rangle$, which we indicate in the following
with $\sigma_{\rho}$, computed in every box. In order to limit the
contamination from shock compression and substructures, we eliminate the cells
corresponding to the top 5% of the density contrast distribution in every box,
following earlier work concerned with X-ray modelling (e.g. Roncarelli et al.
2013). For testing purposes, we also experimented with other thresholds,
obtaining: $\sigma^{2}_{\rho,99}=0.26,\sigma^{2}_{\rho,95}=0.22$ and
$\sigma^{2}_{\rho,90}=0.20$, which refer to thresholds of 99%, 95%, 90%,
respectively.In a similar way, the mean normalised velocity fluctuation,
$\sigma_{v}$, is estimated as the standard deviation of $\delta v/c_{s}$,
defined in the previous section, using the same cells used to calculate
$\sigma_{\rho}$.

Finally, to compare our results on stratified turbulence with previous work,
we also consider the logarithmic density fluctuations,
$s=\ln(\rho/\langle\rho\rangle_{L})$ and calculate the standard deviation in
every box, $\sigma_{s}$. Fig. 2 shows, for each cluster of the sample, the
probability distribution function of the variables we use in this work: the
gas density contrast, $\delta\rho/\bar{\rho}$, the logarithmic gas density
fluctuations, $s=ln(\rho/\bar{\rho})$, and the velocity contrast $\delta
v/c_{s}$. The PDFs are obtained from the fluctuations within the virial radius
(defined as, $R_{200}$, i.e. the radius enclosing a $200$ overdensity with
respect to the cosmic critical density).

Figure 2: Probability distribution function of the gas density contrast,
logarithmic gas density fluctuations, and velocity fluctuations. The colours
correspond to different clusters, at z = 0. All statistics refer to within the
virial radius.

####  2.2.2 Richardson number

The Richardson number, $Ri$, is central to our study, but, in contrast to some
previous work in which the authors considered $Ri$ as an input parameter of
the simulation (Mohapatra et al. 2020, Shi & Zhang 2019), we had to compute it
a-posteriori following Eq. 1. We obtained the radial profile of the Richardson
number in 600 kpc wide boxes, starting from the gravitational acceleration and
gas pressure and density radial profile in the box that we used to estimate
the Brunt-Väisala frequency , $N_{\rm BV}$ (Eq. 2). We also evaluated the
turbulent velocity profile in the box, computing a root mean square after the
excision of the densest density fluctuations (see Sec. 2.2.1). For the sake of
simplicity, we fixed the spatial scale in Eq. 1 as the size $L$ of the
filtering scale (see Sec. 2.2). In order to characterise each box with its
Richardson number, we averaged the $Ri$ radial profile.

Furthermore, we obtained the radial profile of the Richardson number over the
whole volume of each galaxy cluster. The top panel of Fig. 3 shows the
Richardson number profile in our cluster sample for a filtering scale $L$ =
300 kpc (top panel), together with the Brunt-Väisala frequency profile (middle
panel). In all clusters, except one, buoyancy is always found to be dominant
($Ri\gg 1$) within the virial radius, while the turbulent cascade overcomes
gravity only in the outskirts. Moreover, we find that the Richardson number
has a strong dependence on the filtering scale $L$ (Fig. 3, lower panel). It
nearly follows the Kolmogorov scaling $Ri\propto L^{4/3}$. Consequently, a
local estimate of $Ri$ must always be referred to a turbulent scales.

As noted above, in the remainder of this work we restrict ourselves to the
fixed filtering scale of $L=300\rm\leavevmode\nobreak\ kpc$; this is motivated
by our previous works on the Itasca dataset. For example, in Vazza et al.
(2017) and Angelinelli et al. (2020) we measured, with an iterative multi-
scale filtering algorithm, that the kinetic-energy weighted distribution of
cells in the volume of our clusters have a local turbulent outer scale $\leq
400\rm\leavevmode\nobreak\ kpc$, with little dependence on the dynamical state
of the host cluster. Based on these earlier studies, we can conclude that
$L\geq 500\rm\leavevmode\nobreak\ kpc$ appears to be too large based on the
actual kinematic analysis of these systems, while $L\leq
200\rm\leavevmode\nobreak\ kpc$ scales are too small to capture the total
turbulent kinetic energy budget developed in the ICM of these simulated
systems.

Figure 3: Top panel: radial profile of the Richardson number in our clusters.
Middle panel: radial profile of the Brunt-Väisala frequency for each cluster
of the sample. Lower panel: radial profile of the Richardson number
considering different filtering scales $L$, for a single cluster of our
sample. The horizontal red line separates the turbulence-dominated ($Ri$ ¡ 1)
and buoyancy-dominated ($Ri$ ¿1) regimes.

In order to quantify the degree of anisotropy, possibly caused by the
buoyancy, we compute the parameter $\beta$,

$\beta=1-\frac{v^{2}_{t}}{2v^{2}_{r}},$ (5)

where $v_{t}$ and $v_{r}$ are the tangential and radial components of the
turbulent velocity distribution (that is, the r.m.s. values computed in 600
kpc wide boxes, with the filtering scale again set to $L=300$ kpc),
respectively. The anisotropy parameter can span the range $-\infty<\beta\leq
1$. Positive values mean that the dynamics is dominated by the radial motions,
whereas, when $\beta<0$ the dominant components of the velocity field are
those perpendicular to the radial direction. We refer the reader to Vazza et
al. (2018) for a further analysis of the anisotropy parameter in the Itasca
cluster sample.

####  2.2.3 Power spectra

In this section, we compare the density and velocity power spectra at various
cluster radii. We measure the power spectra of the ICM velocity field in large
boxes, $\approx$ 2.5 Mpc (linear size) wide, located at different radii using
a Fast Fourier Transform, with a periodic domain for simplicity (e.g. see
Vazza et al. 2011b, for a discussion):

$\lvert\textbf{v}(\textbf{k})\rvert^{2}=\sqrt{v_{x}(\textbf{k})^{2}+v_{y}(\textbf{k})^{2}+v_{z}(\textbf{k})^{2}},$
(6)

where ${v}_{i}(\textbf{k})$ is the Fourier transform of each component

$v_{i}(\textbf{k})=\frac{1}{\left(2\pi\right)^{3}}\int_{V}v_{i}(\textbf{x})e^{-2\pi
i\textbf{k}\cdot\textbf{x}}\,d\textbf{x}.$ (7)

The gas density fluctuations, $\lvert\rho(\textbf{k})\rvert^{2}$, are obtained
in a similar way, that is by computing the Fourier transform of the real space
density distribution. In this case we did not perform any filtering, since
filtering with a fixed spatial scale would just suppress the power on $k\leq
2\pi/L$ scales, but we select the inertial subrange highlighted with the
orange bars (Fig. 4): 150 kpc $\lesssim l\lesssim$ 1 Mpc. The small scale
limit is being fixed to avoid measuring spectral shapes in a regime already
affected by the numerical dissipation of the PPM hydro scheme. The
normalisation of the power spectra in Fig. 4 has been rescaled to compare the
slopes at various radii. In general, the normalisation of the density spectra
increases outside the virial radius, except for the power spectrum computed in
the outermost region of the cluster which shows a drop. This increasing trend
is likely due to the larger clumping factors in the outskirts (Angelinelli et
al. 2021).

We evaluated the slope of each $k$-bin within the inertial subrange and we
computed an average value for a more robust characterisation of the power
spectrum. For the velocity power spectra we evaluated the slope in each bin of
each velocity component and then we averaged over the all slope values.
Finally, we averaged the slopes of the power spectra performed on boxes at the
same radius, computing an inverse-variance weighted mean, to investigate the
radial trend of the slope of the power spectrum, compared to the $\propto
k^{-11/3}$ trend expected for the Kolmogorov model in three dimensions.

Figure 4: Density and velocity power spectra in a simulated cluster from the
Itasca sample, obtained at different radii. The vertical orange bars indicate
the selected inertial subrange used for the analysis. The normalisation of
density spectra has a radial trend, with larger fluctuations in the outer
regions.

As a preliminary check on the reliability of our cluster simulations for this
study, we computed the statistics of gas fluctuations. As previously
explained, some spurious features that can affect our numerical analysis might
be present in our galaxy clusters (e.g. shocks, clumps or filaments), as well
as to some degree also in real observations. In order to excise most of such
features not associated with true turbulence-driven fluctuations, for each
radial shell we computed the probability distribution function of density and
excluded the 5% densest cells. In order to have uniformly filled maps, we
replaced the original density values of these cells with the median density
within the entire radial shell. This method has been implemented in previous
works (see Vazza et al. 2011a, for an example). Considering that
bremsstrahlung is the dominant radiative mechanism at the typical temperature
of plasma in galaxy clusters ($10^{7}-10^{8}K$), we estimated the surface
brightness by projecting $S_{X}=n^{2}T^{1/2}$ (where $n$ and $T$ are the
plasma density and temperature in each cell) along multiple lines of sight.
Fig. 5 shows the projected X-ray surface brightness in one of our simulated
clusters before and after removing the clumps. We can thus estimate the
surface brightness fluctuations, dividing the X-ray images by their radial
average value within each radial bin. This method is commonly used in
observational studies, even though it does not consider any azimuth
asymmetries of the ICM. Thus, we obtain the 2D power spectra of X-ray
fluctuations on 1 Mpc wide boxes centred on the gas density peak of the
cluster. The characteristic amplitude is computed as follows Churazov et al.
(2012):

$A_{2D}^{S}(k)=\sqrt{2\pi P_{2D}(k)k^{2}},$ (8)

where $P_{2D}(k)$ is the 2D power spectrum and $k$ is the wavenumber in
$kpc^{-1}$ units. The distribution of surface brightness fluctuations in our
cluster sample is shown in Fig. 6. Eckert et al. (2017) computed the
distribution of the amplitude in surface brightness fluctuations for 51 galaxy
clusters, and looked for the possible correlation between larger fluctuations
and the presence of non-thermal radio emission. The results can be used to
investigate the connection between turbulence and the re-acceleration of
relativistic electrons via Fermi II mechanism (e.g. Brunetti & Jones 2014, and
references therein). The red line in Fig. 6 represents the mean value of the
X-ray surface brightness fluctuations in Eckert et al. galaxy cluster sample,
which is 0.9. We observe that our fluctuations are by a factor $\sim 2$
larger.
We should mention that our statistics are not affected by observational
limitations. The surface-brightness fluctuations can be reconstructed for the
entire cluster volume with equal fidelity and ignoring such limiting factors
as the number of emitted X-ray photons. Real observations, however, are
affected by the uneven effective spatial resolution, and by the fixed
sensitivity to X-ray brightness fluctuations. Moreover, it is not
straightforward to implement exactly the same procedure for the masking of
point-like X-ray sources in observations, which can bias the observed
amplitude low by removing a part of gas density enhancement associated with
the surroundings of excised galaxies and AGN.

Figure 5: Projected X-ray surface brightness in a cluster of our sample before
(top panel) and after (bottom panel) filtering out the clumps. Figure 6:
Distribution of the X-ray surface brightness fluctuations, $A_{2D}^{S}$, in
the Itasca cluster sample.

####  2.2.4 Structure function

For the sake of completeness, we also computed the third-order structure
function of the velocity and density for the same reference galaxy cluster
used above for the power spectra analysis, as usually this is more commonly
done in study of other turbulent astrophysical environments (e.g.
2007ApJ...665..416K). Fig. 7 shows the third-order structure function of the
velocity component tangential (upper panel) and longitudinal (central panel)
to the separation, $l$, and of the particle density distribution (lower
panel). The analysis was carried out on 800 kpc wide boxes, located at
multiple distances from the centre of the cluster. We do not use any
previously applied filtering method because the influence of the laminar flow
is dominant mainly at the largest scales ($\gtrsim 1\rm\leavevmode\nobreak\
Mpc$). The maximum separation distance explored here is
$l\sim\rm\leavevmode\nobreak\ 800kpc$. We refer the reader to Vazza et al.
(2011b), (2017) for further details on how the analysis was performed. Our
analysis shows that, overall, both the gas density and velocity third-order
structure functions have a Kolmogorov-like scaling, $S_{3}(l)\propto l$,
within the separation range 100-800 kpc, similarly to the velocity power
spectra. As in previous work (Vazza et al. 2011b), the steepening at the
smallest separations is mostly driven by the increasing dissipation of the PPM
scheme, which dampens fluctuations faster than in the Kolmogorov model. The
normalisation of the velocity structure function increases from the outskirts
to the centre of the cluster, with some variance related to the presence of
substructures. The normalisation of the density structure function, on the
other hand, decreases in a more steady way towards the periphery, following
the drop of the cluster gas density profile (see Fig. 1), which contributes to
$S_{3}(l)$ with the density cube.

Recent studies investigated the structure function in galaxy clusters
(2020ApJ...889L...1L; 2021MNRAS.504..898W; 2021arXiv210901771M). Overall, the
structure functions have a Kolmogorov slope ($\propto l^{1/3}$ for first-order
functions), or steeper. In those studies the turbulence is driven by the black
hole energy output and the observed scales range between 1 and 100 kpc, which
cannot be probed by our simulations. On the other hand, our simulations are
more suitable to capture the development of turbulence injected by mass
accretion on $\sim 0.1-1\rm\leavevmode\nobreak\ Mpc$ scales, and the
complementary view of the structure function confirms that the dynamics of gas
motions are reasonably close to a Kolmogorov model.

Figure 7: Structure function of the tangential (upper panel), transverse
(middle panel) and numerical density (lower panel). The slope is consistent
with Kolmorov’s theory in the range between 100 - 800 kpc, while it has a
steeper behaviour at the smallest scales ($l<100\rm\leavevmode\nobreak\ kpc$).
The normalisation increases with cluster-centric distance, especially for the
density structure function due to the stratification of the ICM.

###  2.3 Testing our method in idealised simulations

Next, we validate our method against more idealised simulations. As previously
discussed, due to the different nature of the simulations, it is not feasible
to apply the very same algorithm that Mohapatra et al. 2020, 2021 used in
their studies to our cluster sample. For instance, in Mohapatra et al. and Shi
& Zhang, the Mach number and the Richardson number were parameters of the
simulation. However, in cosmological clusters these quantities cannot be
estimated a priori. Furthermore, the algorithm to evaluate the density
fluctuations must be revisited as well, due to the presence of self-
gravitating structures inherent in forming clusters, including the Itasca
clusters (see Sec. 2.2). Therefore, we applied our analysis to data from the
simulations by Shi & Zhang. In their work, they used the FLASH code (Fryxell
et al. 2000) to simulate 3D local boxes with
$L_{x}=L_{y}=200\rm\leavevmode\nobreak\ kpc$,
$L_{z}=250\rm\leavevmode\nobreak\ kpc$ and a resolution of $256\times
256\times 320$ cells. The aim of their study was to reproduce local ICM
conditions in various regions for both relaxed and disturbed galaxy clusters.
For this reason, they performed four simulation runs with increasing levels of
stratification corresponding to different Froude numbers ($Fr=\sqrt{1/Ri}$).
They used the Froude number to quantify the role of the buoyancy in each
simulation. The initial velocity field was generated with the purpose to mimic
a 3D Kolmogorov turbulence, with an injection scale $L=43.6$ kpc. They also
removed the compressive modes to generate merely solenoidal turbulence. Given
the initial conditions, the evolution of the velocity distribution is fully
dictated by the Mach and Froude number. In this sense, these simulations are
quite similar to those of Mohapatra et al..

For comparison, we analysed the most stratified cluster data cube at 2 Gyr and
estimated all the relevant quantities, such as the Richardson number and the
logarithmic density fluctuations. In table 1 we compare our results with the
parameters of the snapshot (S&Z method), both analysing the whole cube (full
data cube) and boxes with $80^{3}$ cells located in various spot of the cube
itself. To estimate the Froude number we used the injection scale of their
simulations, which is $L=43.6$ kpc.

| $N_{\rm BV}({\rm Gyr}^{-1})$ | Froude number
---|---|---
S&Z method | 10.6 | 0.34
Full data cube | $13.44\pm 6.86$ | $0.37\pm 0.23$
1st box | $9.08\pm 2.98$ | $0.41\pm 0.11$
2nd box | $9.21\pm 3.04$ | $0.29\pm 0.06$
3rd box | $9.13\pm 3.92$ | $0.45\pm 0.13$
4th box | $8.70\pm 3.48$ | $0.39\pm 0.09$
5th box | $8.70\pm 2.35$ | $0.36\pm 0.06$
6th box | $7.22\pm 2.97$ | $0.47\pm 0.21$
7th box | $7.08\pm 2.85$ | $0.44\pm 0.14$
8th box | $7.18\pm 3.13$ | $0.45\pm 0.17$
9th box | $7.28\pm 3.57$ | $0.48\pm 0.17$
Table 1:  Comparison between the results of our analysis, applied to the Shi &
Zhang (2019) simulation, and their output parameters.

After filtering for density fluctuations, we obtained the logarithmic
perturbations both performing a fixed-scale filtering method, with a scale
$L\approx 43.6$ kpc and computing the mean density on a slice at a given $z$
coordinate, which is the method that Mohapatra et al. 2020, 2021 used in their
studies. We plot the results in Fig. 8 for a comparison between the two
filtering algorithms. The black points are the results of the analysis
performed in a box with $80^{3}$ cells. Finally, we estimated the mean
logarithmic density fluctuation considering the full data cube, obtaining
$\sigma^{2}_{s}=0.048$. This result is in close agreement with Fig. 8 of
Mohapatra et al. (2021), which predicts $\sigma^{2}_{s}\approx 0.04-0.05$ for
a Froude number $Fr\approx 0.3-0.4$. In summary, our results match reasonably
well those in the literature, which instills confidence in our filtering
methods, both for idealised as well as cosmological simulations.

Figure 8: A comparison between the logarithmic density fluctuations estimated
with a fixed scale filtering method and calculating the density mean in slices
at a given $z$ coordinate. The red line is the bisector, which highlights the
one-to-one relation.

##  3 Results

First, we focus on the relation between the density and velocity fluctuations
in our cluster sample. Then, we investigate how the stratification can affect
this relation and test the dependence of the density fluctuation on the
Richardson number. A large part of the analysis was performed within the
virial radius of the galaxy clusters in order to compare to observations.

###  3.1 The $\sigma_{\rho}-\sigma_{v}$ relation

In observations, turbulence in galaxy clusters has been investigated using the
relation between the gas density and velocity fluctuations. This has been
followed up by numerical work in, both, cosmological and idealised simulations
(Zhuravleva et al. 2014b, Gaspari et al. 2014). The goal of those studies was
to derive the amplitude of gas turbulent motions solely from the analysis of
surface brightness fluctuations. Those studies found a one-to-one scaling
relation suggesting that such a relation is valid in, both, perturbed and
relaxed clusters.

Fig. 9 shows the $\delta\rho/\rho-\delta v/c_{s}$ relation in our sample for,
both, relaxed and disturbed clusters. The classification of the dynamical
state is based on the centroid shift parameter, $w$, which measures the offset
between the centre of mass and the position of the gas density peak of the
cluster (Rasia et al. 2013). The centroid shift of each cluster, $\langle
w_{i}\rangle$, was estimated averaging the parameter values of multiple lines
of sight. Galaxy clusters are then classified comparing $\langle w_{i}\rangle$
with the median of the centroid shift of the sample, which is $w_{c}=0.0067$
(which is half of the value used in Cassano et al. 2011), with relaxed
clusters having $w\leq w_{c}$. The relation was obtained within the virial
volume by applying a fixed scale filtering method with scale $L$ = 300 kpc and
averaging the density and velocity fluctuations over 600 kpc wide boxes (see
Sec. 2.2). The outcome is a linear relation that depends on the dynamical
state. Relaxed clusters have a slope $m=1.15\pm 0.06$, where
$\sigma_{v}\propto m\sigma_{\rho}$, consistent with Zhuravleva et al. (2014b).
Contrary to expectations, disturbed clusters show a flatter relation, with
$m=0.40\pm 0.05$.

In order to test the goodness of the fit, we used the Pearson’s correlation
coefficient

$R=\frac{\sum\limits_{i}(x_{i}-\bar{x})(y_{i}-\bar{y})}{\sqrt{\sum\limits_{i}(x_{i}-\bar{x})^{2}\sum\limits_{i}(y_{i}-\bar{y})^{2}}},$
(9)

where $\bar{x}$ and $\bar{y}$ are the means of the variables used for the fit.
$R$ is a measure of linear correlation between two sets of data, but it
ignores other types of relationship. The correlation coefficient ranges
between -1 and 1. These limits represent a perfect (anti)correlation, while a
value of 0 implies that the variables are uncorrelated. For this fit, we
obtained $R=0.43$ and $R=0.82$ for perturbed and relaxed clusters,
respectively.

Figure 9: $\sigma_{\rho}-\sigma_{v}$ relation in our sample, which shows a
dependence on the dynamical state of the galaxy cluster. The classification is
based on the centroid shift parameter. The threshold value is $w_{c}=0.0067$.
The blue line denotes the 1-to-1 relation found by Zhuravleva et al. (2014b).
Figure 10: $\sigma_{\rho}-\sigma_{v}$ relation for the total cluster
sample.The blue line denotes the 1-to-1 relation found by Zhuravleva et al.
(2014b).

Moreover, in Fig. 11 we show the distribution of the ratio of the rms values
of density and velocity fluctuations which depends depends slightly on the
cluster-centric radius, regardless of the dynamical state of the cluster. Blue
and black lines (or points) are related to relaxed and disturbed clusters,
respectively, while the red line represents the relation found by Zhuravleva
et al. (2014b).
The Itasca simulation are non-radiative. This implies that some features of
cool-core clusters, such as the peak (drop) in the density (temperature)
radial profile in the centre of galaxy cluster, are missing. For this reason,
it is not easy to make a proper classification, dividing our clusters in
relaxed or disturbed, in this case. Due to this uncertainty, we performed the
fit of the $\sigma_{\rho}-\sigma_{v}$ relation also considering the total
cluster sample (Fig. 10). We found a flatter relation, compared to a one-to-
one relation, with a slope $m=0.63\pm 0.04$ and a Pearson’s correlation
coefficient $R=0.60$.
Finally, we noticed that some fluctuations that are not driven by turbulence
might affect the $\sigma_{\rho}-\sigma_{v}$ relation, even after filtering. As
a result, some of the $\sigma_{v}$ and $\sigma_{\rho}$ values in Fig. 9 might
be caused by such contamination. Angelinelli et al. (2021) recently studied
the radial profile of the clumping factor and the number density of clump in
the Itasca sample. Here, we performed a similar analysis. We divided the inner
region (within $3R_{500}$) into six bins. For each of them we calculated the
number of clumps as a function of the respective volume. The mass threshold
for the clumps is $10^{8}M_{\odot}$ and we set their maximum size to be 300
kpc. Fig. 12 shows the radial profile of the median clump number density. The
shadowed regions represents the 16th and 84th clump number density
distribution percentile boundaries. The number of clumps in the innermost
region is higher in the perturbed clusters, making it less straightforward to
study the relation between density and velocity fluctuations. However, the
presence of self-gravitating structures is a problem for X-ray observations,
too. In observations the identification of clumps is even more complicated,
due to their low surface brightness. A solution might be to perform the
analysis using a smaller filtering scale, in order to filter out also the
smallest clumps. Nevertheless, the driver of turbulence in this simulation is
the radial accretion. Thus, a smaller filtering scale would not be physically
motivated since the injection scales generated by such event are likely larger
than 200 kpc (Vazza et al. 2018; Angelinelli et al. 2020).

Figure 11: Top panel: ratio of the amplitudes of density and velocity
fluctuations as function of the distance from the cluster centre. Bottom
panel: the distribution of the ratio between density and velocity rms. The
black points (line) correspond to relaxed clusters and blue points (line) to
the perturbed ones. The classification into relaxed and perturbed clusters is
based on the centroid shift parameter, whose discriminating value is
$w_{c}=0.0067$.  Figure 12: Radial profile of the number density of clumps for
both relaxed (blue) and perturbed (red) galaxy clusters. The shadowed region
represent the 16th and 84th percentiles of the bin’s number density
distribution.

###  3.2 The $Ri-\sigma^{2}_{s}$ correlation

In the previous section we discussed the relation between density and velocity
perturbations (Fig. 9). However, the stratification and, consequently, the
buoyancy make the interpretation of gas fluctuations less straightforward.
Hence, we computed the mean of the Richardson number and the logarithmic
density fluctuations over 600 kpc wide boxes using a filtering method with a
fixed scale of $L=300\rm\leavevmode\nobreak\ kpc$ (see Sec. 2.2). In Mohapatra
et al. (2020), the authors found a tight relation between density fluctuations
and the Richardson number (highlighted by the orange line in Fig. 13), which
we cannot detect in our sample (black points in Fig. 13). In particular, our
simulated clusters show higher density fluctuations, up to three orders of
magnitude, for $Ri<1$, and a lower limit on
$\log_{10}(\sigma^{2}_{s})\approx-1.5$, regardless of $Ri$.

We have already presented in Sec. 4 our tests on the reliability of the
density fluctuations computed in the Itasca cluster. Even though our
distribution does not perfectly match that of Eckert et al. (2017), the
deviation in the surface brightness fluctuations is not sufficient to obtain
density fluctuations similar to those of Mohapatra et al. (2020). Moreover, we
have also shown in Sec. 2.3 that when applied to a similar (idealised)
simulation setup, our pipeline correctly recovers the expected density
fluctuations and Froude or Richardson number. We are therefore led to the
conclusion that the setup in Mohapatra et al. 2020 and 2021 is much more
idealised than more realistic ICM as produced by the Itasca sample, both
concerning the solenoidal forcing of turbulence and the absence of self-
gravitating gas substructures and bulk motions associated with accretion. In
our cluster sample, the accretion of clumps is ubiquitous and it produces both
solenoidal and compressive turbulence. This results in a distribution of flow
velocities with typical averaged (on 600 kpc wide boxes) Mach numbers of
$\approx 0.4-0.5$ and likely locally supersonic motions (see Fig. 14). Another
crucial difference is the distribution of the logarithmic density
fluctuations, which is found to be significantly different from a log-normal
distribution in our work. In summary, in our simulations we find that density
fluctuations do not depend on $Ri$. Therefore, estimating the density
perturbations as the variance of the probability distribution function, and
relating it to the physical properties of the underlying turbulent flow,
appears to be challenging in realistic situations, even without taking into
account observational effects and sampling issues.

Figure 13: $Ri$ \- $\sigma^{2}_{s}$ plot. Density perturbations are
independent of the Richardson number and much higher than those found by
Mohapatra et al. (orange line). The analysis was performed on the total
cluster sample.

###  3.3 The $Ri-\beta$ correlation

The stratification and buoyancy might also affect the velocity distribution of
the plasma, suppressing the radial motion. As a result, in the buoyancy-
dominated regime ($Ri>1$), the gas dynamics should be dominated by tangential
motions (see Brethouwer et al. 2007, Mohapatra et al. 2020, 2021). Therefore,
we look for a relation between the Richardson number and the anisotropy
parameter. Since for $Ri<1$ the hypotheses of homogeneity and isotropy are
still valid, we should expect a mean $\beta\sim 0$. On the other hand, a mean
$\beta<0$ should arise when we move to buoyancy-dominated regimes. However,
looking at Fig. 15, we see that the dynamics of the ICM is biased by radial
motions despite the stratification. As already noticed in previous work on
this sample (e.g. Vazza et al. 2018), this ubiquitous dynamical feature of the
ICM is most likely a consequence of radial accretion flows that are important
drivers of ICM turbulence. Finally, we have explored the sensitivity of our
results on the choice of the filtering scale and box size, by considering
larger and smaller filtering scales. We find that qualitatively the results
still hold (see Appendix 1 for more details).

The important conclusion is that the Richardson number is not able to
characterise the effect of ICM stratification for cosmologically evolving
clusters. Since accretion occurs largely via the infall of gas clumps whose
motion is predominantly radial (e.g. Angelinelli et al. 2021), a predominantly
radial (turbulent) velocity is found at most cluster-centric distances in all
clusters of our sample. This radial turbulence is unlike the stratified
turbulence expected for a $Ri\geq 1$ flow.

Figure 14: Distribution of the Mach numbers in our cluster sample, using
averaging over 600 kpc. The tail of the distribution at high Mach numbers is
caused by local, supersonic turbulence.

###  3.4 Radial profile of the power spectra slope

The velocity power spectrum quantifies the strength of turbulent fluctuations
on different scales, and the slope of the power spectrum is customarily
compared to Kolmogorov’s theory.

Figure 15: $Ri-\beta$ relation in our cluster sample. The radial bias is an
imprint that the dynamics is dominated by radial motions, driven by the
accretion. In this sense, the Richardson number does not have a predictive
power about the real dynamics of the ICM. The orange line is used to highlight
the turbulence (left) and buoyancy (right) dominated regime, while the blue
line distinguishes a radially ($\beta>0$) from a tangentially ($\beta<0$)
bias.

We compute the radial profile of slope of the density and velocity power
spectra in the Itasca cluster sample, which is shown in Fig. 16. The error-
bars represent the $1\sigma$ error, which were estimated as the standard error
of the weighted mean
$\sigma=\left(\sqrt{\sum\limits_{i}\frac{1}{\sigma_{i}^{2}}}\right)^{-1}$,
where $\sigma_{i}$ are the $1\sigma$ errors of each power spectrum at a given
distance from the centre of the cluster. The slope of the velocity spectrum
shows a constant radial profile and agrees with the Kolmogorov’s model,
indicated by the blue line. On the other hand, the density power spectra have
steeper slopes, albeit with a bigger scatter. The slope of the velocity
spectrum is flatter compared to Vazza et al. (2017), who used the same cluster
sample. The reason for this difference is explained by the varying
normalisation of the power spectrum, which decreases as we move to larger
radii (see Fig. 4). While in previous work we measured a unique power spectrum
across a large volume (hence mixing local spectra with a slightly decreasing
normalisation going to the cluster outskirts), here we computed spectra in
smaller boxes. On scales of $\leq 600\rm\leavevmode\nobreak\ kpc$ the
turbulence appears even more Kolmogorov-like, despite the presence of radial
biases and a non-stationary forcing. At the same time, the larger scatter in
the slope of the density power spectrum at all radii questions whether any
simple relation can be robustly derived between the two quantities.

Figure 16: Radial profile of the slope of density and velocity power spectra
in our cluster sample. The velocity slope is in agreement with Kolmogorov’s
theory and shows a constant profile as well as the slope of the density
spectra, which is steeper. The blue line traces Kolmogorov’s prediction for
the velocity slope of a 3D spectrum.

##  4 Discussion and conclusion

We have investigated whether the turbulent gas velocity fluctuations in the
ICM can be constrained using projected X-ray surface brightness fluctuations.
Previous work (Gaspari et al. 2014; Zhuravleva et al. 2014b; Schuecker et al.
2004; Sanders & Fabian 2012) have already shown that the pattern of surface
brightness fluctuations can, under certain limiting assumptions, be translated
into density fluctuations. However, in a stratified ICM, the
$\delta\rho-\delta v$ relation is affected by buoyancy effects. Mohapatra et
al. (2020), (2021) further derived an accurate relation between the two,
parameterised by the Richardson number and the Mach number of the flow, in the
very simplified case of turbulence-in-a-box simulations. So far, such analyses
have been limited to idealised simulations or to the innermost cluster
regions. Here, we used the Itasca cluster sample, which is a set of galaxy
clusters simulated with the cosmological code ENZO with a uniform comoving
spatial resolution of about $\approx 20$ kpc. In order to compare the results
to observations, we limited our analysis of the ICM to within the virial
radius.

Fig. 4 shows that while the density fluctuations increase as we move to larger
radii, the velocity power spectra remain fairly constant. The increasing
clumping factor in the outer region of galaxy clusters is the likely driver
for this trend. As a consequence, a more accurate algorithm to mask clumps is
required in observational and numerical analyses. This will be particularly
challenging for observations, where the masking small-scale and faint clumps
will remain limited. We can summarise our results as follows:

  * $\bullet$

Consistent with previous work (Zhuravleva et al. 2014b, Gaspari et al. 2014),
we find that a linear relation exists between the standard deviation of the
density contrast, $\sigma_{\rho}$ and the normalised velocity fluctuations,
$\sigma_{v}$, albeit with a significant scatter across the sample. However,
the slope of this correlation does not generally agree with the findings of
the aforementioned studies. For the whole Itasca sample (Fig. 10) it results
$\sigma_{v}\approx 0.6\sigma_{\rho}$, with smaller velocity fluctuations than
expected by applying the Gaspari et al. (2014) relation.

When we separate the sample in relaxed and perturbed clusters, we find that
relaxed objects show a steeper slope, consistent with 1, and again a large
scatter. Perturbed clusters follow instead a flatter, weaker relation,
$\sigma_{v}\approx 0.4\sigma_{\rho}$, with a low correlation coefficient.

Fig. 12 shows that this relation may be contaminated by unfiltered clumps.
This demonstrates the problem in using this relation to infer the turbulence
in realistic clusters since this could lead to incorrect estimates, for
example, of the turbulent heating rate (Zhuravleva et al. 2014a) and the
kinetic energy available to explain radio halos
($P_{H}\propto\rho\sigma_{v}^{3}$, e.g. Eckert et al. 2017).

  * $\bullet$

Although buoyancy is the dominant process in the largest part of the cluster
volume (Fig. 3), we did not find a strong relation between the logarithmic
density fluctuations and the Richardson number. The results in Sec. 3.2 show
that our fluctuations are much higher compared to those of Mohapatra et al.
(2020), (2021), who performed an analogous analysis in more idealised
simulations, and are independent of the Richardson number. We attribute the
difference to the unavoidable presence of compressive turbulence, consistent
with previous work (e.g. Miniati 2015; Vazza et al. 2017; 2021MNRAS.504..510V)
. Furthermore, the turbulent anisotropy parameter shows the effects of
substantial radial accretion (see Fig. 15), which appears to be the dominant
driver for the gas fluctuations in our simulations (e.g. Angelinelli et al.
2021).

  * $\bullet$

Turbulence in the ICM is usually assumed to be homogeneous and isotropic. We
show that this may not be valid, even if the spectra are found to follow
Kolmogorov’s scaling. This is primarily caused by radial accretion that
dominates the dynamics of the ICM (see Fig. 15). In cluster cores, jets from
central radio galaxies may produce a similar effect (Brüggen & Vazza 2015).

Finally, we point out that our cosmological simulations are optimised for a
uniform resolution across the cluster volume and cannot attain the high
spatial resolution of observations in cluster cores. Therefore, our
simulations do not extend to scales below $\leq 10\rm\leavevmode\nobreak\ kpc$
probed by some X-ray observations of nearby clusters of galaxies.

## ACKNOWLEDGEMENTS

The authors thank Xun Shi for providing data for subsection 2.3. The
cosmological simulations described in this work were performed using the ENZO
code (http://enzo-project.org), which is the product of a collaborative effort
of scientists at many universities and national laboratories. F.V.
acknowledges financial support from the European Union’s Horizon 2020 program
under the ERC Starting Grant ’MAGCOW’, no. 714196. Work by TWJ was supported
by the US National Science Foundation grant AST1714205. MB acknowledges
support by the Deutsche Forschungsgemeinschaft (DFG, German Research
Foundation) under Germany’s Excellence Strategy – EXC 2121 “Quantum Universe”
– 390833306, Germany’s Excellence Strategy – EXC-2094 “Origins” – 390783311.

## References

  * Angelinelli et al. (2020) Angelinelli M., Vazza F., Giocoli C., Ettori S., Jones T. W., Brunetti G., Brüggen M., Eckert D., 2020, MNRAS, 495, 864
  * Angelinelli et al. (2021) Angelinelli M., Ettori S., Vazza F., Jones T. W., 2021, arXiv e-prints, p. arXiv:2102.01096
  * Arévalo et al. (2016) Arévalo P., Churazov E., Zhuravleva I., Forman W. R., Jones C., 2016, ApJ, 818, 14
  * Bonafede et al. (2018) Bonafede A., et al., 2018, MNRAS, 478, 2927
  * Brethouwer et al. (2007) Brethouwer G., Billant P., Lindborg E., Chomaz J. M., 2007, Journal of Fluid Mechanics, 585, 343
  * Brüggen & Kaiser (2002) Brüggen M., Kaiser C. R., 2002, Nature, 418, 301
  * Brüggen & Vazza (2015) Brüggen M., Vazza F., 2015, Turbulence in the Intracluster Medium. p. 599, doi:10.1007/978-3-662-44625-6˙21
  * Brunetti & Jones (2014) Brunetti G., Jones T. W., 2014, International Journal of Modern Physics D, 23, 1430007
  * Bryan et al. (2014) Bryan G. L., et al., 2014, ApJS, 211, 19
  * Cassano et al. (2011) Cassano R., Brunetti G., Venturi T., 2011, Journal of Astrophysics and Astronomy, 32, 519
  * Churazov et al. (2012) Churazov E., et al., 2012, MNRAS, 421, 1123
  * Colella & Woodward (1984) Colella P., Woodward P. R., 1984, Journal of Computational Physics, 54, 174
  * De Young (2010) De Young D. S., 2010, ApJ, 710, 743
  * Dolag et al. (2005) Dolag K., Vazza F., Brunetti G., Tormen G., 2005, MNRAS, 364, 753
  * Eckert et al. (2017) Eckert D., Gaspari M., Vazza F., Gastaldello F., Tramacere A., Zimmer S., Ettori S., Paltani S., 2017, ApJ, 843, L29
  * Eckert et al. (2019) Eckert D., et al., 2019, A&A, 621, A40
  * Ettori et al. (2019) Ettori S., et al., 2019, A&A, 621, A39
  * Fabian et al. (2017) Fabian A. C., Walker S. A., Russell H. R., Pinto C., Sanders J. S., Reynolds C. S., 2017, MNRAS, 464, L1
  * Fryxell et al. (2000) Fryxell B., et al., 2000, ApJS, 131, 273
  * Gaspari et al. (2013a) Gaspari M., Ruszkowski M., Oh S. P., Churazov E., Brighenti F., Ettori S., Sharma P., Temi P., 2013a, in Markevitch M., ed., SnowCLUSTER 2013, Physics of Galaxy Clusters. p. 88
  * Gaspari et al. (2013b) Gaspari M., Brighenti F., Ruszkowski M., 2013b, Astronomische Nachrichten, 334, 394
  * Gaspari et al. (2014) Gaspari M., Churazov E., Nagai D., Lau E. T., Zhuravleva I., 2014, A&A, 569, A67
  * Gaspari et al. (2018) Gaspari M., et al., 2018, ApJ, 854, 167
  * Hitomi Collaboration et al. (2016) Hitomi Collaboration et al., 2016, Nature, 535, 117
  * Khatri & Gaspari (2016) Khatri R., Gaspari M., 2016, MNRAS, 463, 655
  * Kim & Ryu (2005) Kim J., Ryu D., 2005, ApJ, 630, L45
  * Kolmogorov (1941) Kolmogorov A., 1941, Akademiia Nauk SSSR Doklady, 30, 301
  * Komatsu et al. (2011) Komatsu E., et al., 2011, ApJS, 192, 18
  * Lau et al. (2009) Lau E. T., Kravtsov A. V., Nagai D., 2009, ApJ, 705, 1129
  * Markevitch & Vikhlinin (2007) Markevitch M., Vikhlinin A., 2007, Phys. Rep, 443, 1
  * Miniati (2015) Miniati F., 2015, ApJ, 800, 60
  * Mohapatra et al. (2020) Mohapatra R., Federrath C., Sharma P., 2020, MNRAS, 493, 5838
  * Mohapatra et al. (2021) Mohapatra R., Federrath C., Sharma P., 2021, MNRAS, 500, 5072
  * Nelson et al. (2014) Nelson K., Lau E. T., Nagai D., 2014, ApJ, 792, 25
  * Olivares et al. (2019) Olivares V., et al., 2019, A&A, 631, A22
  * Randall et al. (2015) Randall S. W., et al., 2015, ApJ, 805, 112
  * Rasia et al. (2013) Rasia E., Meneghetti M., Ettori S., 2013, The Astronomical Review, 8, 40
  * Rebusco et al. (2006) Rebusco P., Churazov E., Böhringer H., Forman W., 2006, MNRAS, 372, 1840
  * Roediger et al. (2013) Roediger E., Kraft R. P., Nulsen P., Churazov E., Forman W., Brüggen M., Kokotanekova R., 2013, MNRAS, 436, 1721
  * Roncarelli et al. (2013) Roncarelli M., Ettori S., Borgani S., Dolag K., Fabjan D., Moscardini L., 2013, MNRAS, 432, 3030
  * Russell et al. (2019) Russell H. R., et al., 2019, MNRAS, 490, 3025
  * Sanders & Fabian (2012) Sanders J. S., Fabian A. C., 2012, MNRAS, 421, 726
  * Sanders et al. (2010) Sanders J. S., Fabian A. C., Smith R. K., Peterson J. R., 2010, MNRAS, 402, L11
  * Schuecker et al. (2004) Schuecker P., Finoguenov A., Miniati F., Böhringer H., Briel U. G., 2004, A&A, 426, 387
  * Shi & Zhang (2019) Shi X., Zhang C., 2019, MNRAS, 487, 1072
  * Valdarnini (2019) Valdarnini R., 2019, ApJ, 874, 42
  * Vallés-Pérez et al. (2021) Vallés-Pérez D., Planelles S., Quilis V., 2021, MNRAS,
  * Vazza et al. (2009) Vazza F., Brunetti G., Kritsuk A., Wagner R., Gheller C., Norman M., 2009, A&A, 504, 33
  * Vazza et al. (2011a) Vazza F., Roncarelli M., Ettori S., Dolag K., 2011a, MNRAS, 413, 2305
  * Vazza et al. (2011b) Vazza F., Brunetti G., Gheller C., Brunino R., Brüggen M., 2011b, A&A, 529, A17
  * Vazza et al. (2012) Vazza F., Roediger E., Brüggen M., 2012, A&A, 544, A103
  * Vazza et al. (2017) Vazza F., Jones T. W., Brüggen M., Brunetti G., Gheller C., Porter D., Ryu D., 2017, MNRAS, 464, 210
  * Vazza et al. (2018) Vazza F., Angelinelli M., Jones T. W., Eckert D., Brüggen M., Brunetti G., Gheller C., 2018, MNRAS, 481, L120
  * Wang & Markevitch (2018) Wang Q. H. S., Markevitch M., 2018, ApJ, 868, 45
  * Wittor et al. (2017) Wittor D., Jones T., Vazza F., Brüggen M., 2017, MNRAS, 471, 3212
  * Zhuravleva et al. (2014a) Zhuravleva I., et al., 2014a, Nature, 515, 85
  * Zhuravleva et al. (2014b) Zhuravleva I., et al., 2014b, ApJ, 788, L13
  * Zhuravleva et al. (2018) Zhuravleva I., Allen S. W., Mantz A., Werner N., 2018, ApJ, 865, 53

##  5 Appendix 1: On the $Ri-\sigma^{2}_{s}$ relation at different filtering
scales

In Sec. 2.2.2 we showed that the Richardson number is very dependent on the
filtering scale. In order to validate the robustness of our results, we
performed tests with a filtering scale of 100 and 800 kpc. Due to the
resolution of our simulation, an analysis with a much smaller scale cannot be
accomplished.

Fig. A1 and A2 show the $Ri-\sigma^{2}_{s}$ and $Ri-\beta$ relation for
$L=100$ kpc and $L=800$ kpc, respectively. In essence, the results are almost
identical to those of Sec. 3.2, as we did not detect any strong relation
between the variables. It looks there is a $Ri-\sigma^{2}_{s}$ relation for
$L=800$ kpc. However, there is a lot a scatter and the fit would definitely be
different from that of Mohapatra et al. since the density fluctuations are
much higher. Furthermore, a filtering algorithm with such a large scale, it
might not be a properly choice, since we are essentially removing the bulk
motions on scale of typical size of the cluster virial radius, leaving out the
laminar flows on smaller scales. Besides, the bottom panel of Fig. A2 shows
that a radial bias is still present, confirming that the Richardson number is
not the ideal proxy to quantify the effect of stratification on turbulence in
galaxy clusters.

Fig. A1: $Ri-\sigma^{2}_{s}$ and $Ri-\beta$ plot for a filtering scale $L=100$
kpc. The orange line is used to highlight the turbulence (left) and buoyancy
(right) dominated regime, while the blue line distinguishes a radially
($\beta>0$) from a tangentially ($\beta<0$) bias.

Fig. A2: $Ri-\sigma^{2}_{s}$ and $Ri-\beta$ plot for a filtering scale $L=800$
kpc. The orange line is used to highlight the turbulence (left) and buoyancy
(right) dominated regime, while the blue line distinguishes a radially
($\beta>0$) from a tangentially ($\beta<0$) bias.