marginparsep has been altered.
topmargin has been altered.
marginparwidth has been altered.
marginparpush has been altered.
The page layout violates the ICML style. Please do not change the page layout,
or include packages like geometry, savetrees, or fullpage, which change it for
you. We’re not able to reliably undo arbitrary changes to the style. Please
remove the offending package(s), or layout-changing commands and try again.

Analyzing the Performance of Graph Neural Networks
with Pipe Parallelism

Anonymous Authors1

###### Abstract

Many interesting datasets ubiquitous in machine learning and deep learning can
be described via graphs. As the scale and complexity of graph-structured
datasets increase, such as in expansive social networks, protein folding,
chemical interaction networks, and material phase transitions, improving the
efficiency of the machine learning techniques applied to these is crucial. In
this study, we focus on Graph Neural Networks (GNN) that have found great
success in tasks such as node or edge classification and link prediction.
However, standard GNN models have scaling limits due to necessary recursive
calculations performed through dense graph relationships that lead to memory
and runtime bottlenecks. While new approaches for processing larger networks
are needed to advance graph techniques, and several have been proposed, we
study how GNNs could be parallelized using existing tools and frameworks that
are known to be successful in the deep learning community. In particular, we
investigate applying pipeline parallelism to GNN models with GPipe, introduced
by Google in 2018.

††footnotetext:  1Anonymous Institution, Anonymous City, Anonymous Region,
Anonymous Country. Correspondence to: Anonymous Author
<EMAIL_ADDRESS>
Preliminary work. Under review by the Machine Learning and Systems (MLSys)
Conference. Do not distribute.

##  1 Introduction

Traditional flat or sequential data delivery cannot fully satisfy many of
today’s demanding deep learning models, especially as more data structures of
interest can be better represented as high-dimensional graphs instead of low-
dimensional grids. Graph machine learning has demonstrated successful
applications in domains such as chemistry and drug design Duvenaud et al.
(2015); Mercado et al. (2020), natural language processing Vashishth (2019),
spatio-temporal forecasting Yu et al. (2017), security Zhou et al. (2020),
social networks Zhou et al. (2020), knowledge graphs Arora (2020), recommender
systems Ying et al. (2018), protein design discovery Strokach et al. (2020),
and material phase transitions Bapst et al. (2020). With its increased use,
especially on large datasets, performance and scaling challenges with the
Graph Neural Network (GNN) Zhou et al. (2018) are becoming prevalent when
using existing machine learning frameworks and accelerators because of memory
and data movement limitations Auten et al. (2020)

A variety of new solutions to address these issues have been proposed and are
highlighted in Section 3. However, as a practical consideration, leveraging
existing state-of-the-art tools and frameworks with demonstrated success at
improving deep neural network performance is valuable to push these
technologies forward. Therefore, to better understand how training and
inference can be more efficient for GNN models, we implement and analyze the
performance of parallelized GNN models compared to their unparallelized
counterparts when trained on a single CPU and GPU and multiple GPUs. As GNNs
are executed sequentially, either layer-by-layer or stage-by-stage, the
motivation for this study is to extend current techniques for improving
performance by introducing pipeline parallelism into the GNN model
architecture Zhang et al. (2020).

##  2 Graph Neural Networks

The success of deep learning on traditional grid- or sequence-based inputs,
such as images and sentences, cannot be overstated. Nevertheless, many
datasets in the real-world cannot be expressed within a Euclidean coordinate
system, and instead naturally take an arbitrary form of graphs or networks.
Various studies exist on how to generalize neural networks for the application
to arbitrary irregular graphs Bruna et al. (2013); Henaff et al. (2015);
Duvenaud et al. (2015); Li et al. (2015); Defferrard et al. (2016); Kipf &
Welling (2016), and we follow the exposition of Kipf and Welling Kipf &
Welling (2016) who first introduced the notion of a convolution architecture
for a graph.

We consider the general scenario of node classification in a graph. The input
is a graph $G=(V,E)$ where $V$ is the set of vertices and $E$ is the set of
edges. Each node (or vertex) $i$ in $V$ has a set of features $x_{i}$. Some
pairs of nodes are connected, and the connections are called edges (or links)
and form the set $E$. If $n$ is the number of nodes and there are $d$ features
in each node, then the set of all features form a $n\times d$ matrix. Some
nodes may have labels from a set of $C$ classes. The task is to classify each
node into $C$ classes using the ingrained feature information through the edge
connectivity between nodes across the graph.

Table 1: Comparisons of training dataset sizes used in this work and considered for future experimentation. Dataset | Nodes | Edges | Classes
---|---|---|---
Cora | 2,708 | 5,429 | 7
CiteSeer | 3,312 | 4,732 | 6
PubMed | 19,717 | 44,338 | 3
Reddit | 233,000 | 150,000,000 | 50
Amazon | 6,000,000 | 180,000,000 | 11

For our analysis, we use the Cora, CiteSeer, and PubMed datasets which are
well-established citation network datasets Sen et al. (2008); Yang et al.
(2016) often used in benchmark training. For comparison of these in Table 1,
we also list the sizes of two larger datasets, the Reddit post dataset
Hamilton et al. (2017) and the Amazon data dump McAuley et al. (2015). As the
data sets used in this study are small and do not require the techniques
explored here to provide training efficiency, they offer valuable benchmarks
for GNNs and measuring the baseline efficiency of parallelized models.

A GNN approaches the node classification problem by building a neural network
layer atop a simultaneous message passing paradigm. Suppose there are $L+1$
layers, $H_{0},\ldots,H_{L}$. Then, $H_{0}=X$, the input set of features. For
each layer $l$, the set of output features $H^{l+1}$ depends only on the
previous layer $H^{l}$ and the input graph $G$. So, for some efficiently
computable function $f$, we have $H^{l+1}=f(H^{l},G)$. Implementations of the
GNN model considers different choices of $f$. By setting the output of the
last layer to be a single neuron, the model computes the logits for each node,
which is then used to classify the nodes. Assuming $f$ is differentiable, this
approach can be optimized by standard gradient descent algorithms. In most
graph networks studied in the literature, the features $H^{l+1}(i)$ for node
$i$ depend on the original feature $H^{l}(i)$ and the features of neighboring
nodes $H^{l}(j)$, where $j$ being a neighbor of $i$ connected by edge
{$i$,$j$}. In some settings, these edges can be of different types, such as
being directed or undirected, as well as include features, in which case the
edge properties also contribute to the calculation of $f$. The message passing
paradigm is designed by how the output features of a layer are simultaneously
updated based on the input features of the layer, instead of through
sequential updates. For many learning tasks on a graph, earlier approaches
Dijkstra (1959); Cheriton & Tarjan (1976) usually introduced problem-specific
architectures or spectral graph theory to make predictions. However, these
algorithms are limited as they require prior knowledge of the graph structure,
and the GNN model provides a unified approach that allows for studying the
properties of a graph itself.

The experiments presented in this paper are based on the Graph Attention
Network (GAT), which is a novel GNN architecture that use attention layers on
top of graph convolutions Veličković et al. (2017). By leveraging this self-
attention mechanism, GAT can achieve state-of-the-art accuracy results on
several transductive and inductive graph benchmarks including the Cora,
CiteSeer, and Pubmed datasets. We chose this model because it can be applied
to graphs of different degrees by specifying arbitrary weights on the
neighbors, making it useful for a wide variety of graph datasets. GAT is also
applicable to inductive learning problems and has generalization capabilities
to unseen graph data. Through our experiments, it does not perform as
efficiently as simpler Graph Convolutional Networks (GCN), such as those
inspired by Kipf & Welling (2016), lending it to being an illustrative
technique to inform us on parallelization benchmarks.

##  3 Related Work

As described above, the core message-passing function of a GNN is the
aggregation of features from the neighborhood of each node. Computing a
gradient descent operation requires storing the entire graph as a single
training batch. With increasing graph size or edge density, the time and
memory complexity of this computation can grow exponentially and introduce an
information bottleneck Alon & Yahav (2020). This effect limits the scalability
of traditional GNNs, many of which, including GAT, do not address these
concerns as their benchmark datasets were of a reasonable size for the
available device capacities.

GraphSAGE Hamilton et al. (2017) was the first attempt to address graph
scalability by using a neighborhood sampling with mini-batch training to limit
the number of nodes included in a batch. This approach can introduce redundant
calculations as the same nodes may appear in multiple sampled batches and lead
to ”neighbor explosion” Zeng et al. (2020). Similar sampling techniques
responded to this challenge by batching sub-graphs instead of nodes, such as
in Chiang et al. (2019). However, graph clustering approaches are faced with
the challenge of defining sub-graphs that sufficiently preserve the edge
topology that guides the node feature updates during training, which is an
issue directly observed in the analysis of the present study. NeuGraph Ma et
al. (2019) introduced parallel computation to enable GNN scalability through a
new graph-parallel abstraction of Scatter-ApplyEdge-Gather-ApplyVertex (SAGA-
NN). This framework conveniently encapsulates almost all existing GNNs in the
literature Zhang et al. (2020), and serves as a foundation for studying
parallelized GNN performance optimization. The authors explored the design of
a GNN processing framework on top of a dataflow based deep learning system,
through which they optimized graph computations, scheduling, and parallelism
in a dataflow-based deep learning framework for graphs. Exploring
computational efficiencies at the device level, $G^{3}$ Liu et al. (2020) is a
GNN training framework that implements graph-structured parallel operations
that leverage the architectural features on GPUs. By directly utilizing graph-
aware components of the GPU, they demonstrated significant speedups in
training times over standard implementations in PyTorch and TensorFlow.
Finally, a recent approach that avoids graph sampling over nodes or sub-graphs
is the Scalable Inception Graph Neural Network (SIGN) Frasca et al. (2020).
Here, graph convolutional filters of different sizes precompute intermediate
node representations. This method enables its scaling to large graphs with
classic mini-batching because it retains sufficient expressiveness from the
node relationships for effective learning.

##  4 Pipeline Parallelism

Google Brain introduced the scalable pipeline parallelism library GPipe Huang
et al. (2019) to enable the efficient distributed training of large, memory-
consuming deep learning models on current accelerator architectures. According
to their published results, GPipe increased training times of a 557 million-
parameter model by 25 times using eight TPU devices and 3.5 times faster using
four devices.

GPipe configures a distribution of a sequentially-designed deep neural network
across multiple devices. To maximize these devices’ capability to calculate in
parallel, it further splits the input mini-batch from the training samples
into “micro-batches” to distribute across the devices. This micro-batching
technique reduces the load for the available memory on the accelerators,
resulting in effectively simultaneous training of the same batch across all
devices. This approach to pipeline parallelism is like a stacked combination
of model parallelism and small data parallelism. During the forward pass, when
each partition finishes processing a micro-batch, it shifts the output to the
next partition, then immediately begins work on the next micro-batch, enabling
partitions can overlap across GPUs. During the backward pass, the gradients
for each micro-batch are calculated using the same model parameters from the
forward pass, and are consistently accumulated at the end into the mini-batch
to update the model parameters. Therefore, the number of partitions separating
the data does not affect model quality.

Although designed for deep neural networks, the GPipe workflow is applicable
to GNNs, with some necessary adaptations that we explore in this study. To the
best of our knowledge, this is the first work to consider the idea of applying
pipeline parallelism using existing libraries to potentially optimize the
runtime of GNNs.

##  5 Implementation

Experiments included training a GAT-based multi-layer sequential neural
network on the task of node classification with the citation datasets
described above. The PubMed set was solely used to compare performance with
pipeline parallelism and graph data batching. The forward propagation model
structure remained consistent across all experiments designed with a drop-out
layer (p = 0.6) followed by a GAT layer with eight heads (attention drop-out =
0.6), a leaky ReLU activation function, a second drop-out layer (p = 0.6), a
second GAT layer (eight heads, attention drop-out = 0.6) where the outputs are
averaged, and, finally, a log softmax function. The neural network model was
implemented in PyTorch with the graph frameworks, PyTorch Geometric (PyG) Fey
& Lenssen (2019), and Deep Graph Library (DGL) Wang et al. (2019). Each
framework was compared for performance on each device architecture. Pipeline
and data parallelism through GPipe was only implemented through DGL. Trials
included performance measures for a single CPU, single GPU, and pipe parallel
distribution across four GPUs with and without micro-batching of the graph
data. For the single device benchmarks, we used an Intel(R) Xeon(R) CPU @
2.20GHz and NVIDIA Tesla T4 GPU, and four NVIDIA Tesla V100-SXM2 GPUs (DGX)
were used for the distributed pipeline parallel experiments.

GPipe was incorporated into the GNN models for each framework with the
torchgpipe Kim et al. (2020) library, an implementation of the GPipe framework
in PyTorch. The defined model is wrapped in a method that takes as parameters
a defined distribution of the model layers across the available GPUs as the
number of micro-batches (called “chunks”) to be applied. A value of one
corresponds to the data parallelism feature being disabled.

In Listing LABEL:code:gpipe, `g` contains the complete graph, `numfeats`
represents the number of features per node, and `nclasses` is the number of
classes in the classification task. The customizable `balance` array specifies
how many layers from the sequence to distribute to each GPU. From this, GPipe
automatically manages the necessary movements of the model and data across
devices. An automated distribution algorithm is also available to optimize
this layer assignment. However, for the uniform analysis presented in this
paper, we manually set the layer distribution across four devices to ensure
consistency for all experiments. With `chunks > 1`, the complete dataset or
batches from the training are split into micro-batches by GPipe to increase
device parallelism. After a partition completes its processing of a micro-
batch, it passes the output to the next partition and begins on the next
incoming micro-batch in the pipeline. Through this approach, the multiple
devices effectively process the same mini-batch (or entire dataset)
simultaneously during a training epoch.

Listing 1: Illustrative GPipe implementation with torchgpipe.

⬇

import torch.nn as nn

from torchgpipe import GPipe

# Define a sequential model

model = nn.Sequential(

nn.Dropout(0.6),

GAT(g, numfeats, 8),

nn.ELU(),

nn.Dropout(0.6),

GAT(g, 8 * 8, nclasses, take_mean = True),

nn.LogSoftmax(1))

# Wrap the model for pipeline parallelism management

model = GPipe(model, balance = [1, 2, 1, 2], chunks=4)

A key challenge with this implementation for a GNN is that a sequential module
is required for the network layers. A cascade of additional restrictions
results, beginning with only a single input of features may be passed through
the layers. However, the graph convolution layer expects as input the graph
data object and its corresponding features. For our experiments that did not
incorporate model parallelism across multiple GPUs, this condition did not
pose an issue because we could simply include the full graph data object, `g`,
into the GNN model definition and pass the single tensor of features. However,
when model parallelism is activated, GPipe applies micro-batching to this
feature tensor, and the corresponding subset of graph nodes must instead be
presented to the graph convolution layer, instead of the full graph data
object.

As a workaround for enabling micro-batching, we exploited the option that the
sequential module can pass a single tuple comprised of multiple tensors. Then,
we pass the node indices of the graph as the first tensor along with the
corresponding features in a second tensor. GPipe applies its micro-batching to
each tensor in the tuple, and a subset of graph nodes with the corresponding
features are passed along the sequence of layers, as needed. When the graph
convolution layer receives the passed tuple, our adapted code extracts the
node tensor comprised of the sub-graph as determined by the micro-batch from
GPipe. Both DGL and PyG graph frameworks include a method to re-build a graph
structure from a subset of graph nodes, which requires the full graph data
object, `g`, for the re-build. The output is then a sub-graph structure
expected by the graph convolution layer. The second tensor of the passed tuple
that includes the features is subsequently extracted in the graph convolution
layer and used in the forward calculation. Upon completion, the two-tensor
tuple is reformed with the original nodes of the sub-graph and the updated
features to be passed along through the remaining layers of the sequence.

##  6 Results

Table 2: Benchmark results on multiple compute architectures and graph frameworks for the standard, small graph datasets. Compute | _Ave. epoch (ms) $|$ Test accuracy_
---|---
Package | Cora | CiteSeer | PubMed
CPU – PyG | 104.4 $|$ 0.717 | 182.9 $|$ 0.696 | 798.5 $|$ 0.718
CPU – DGL | 71.3 $|$ 0.785 | 153.4 $|$ 0.710 | 338.6 $|$ 0.723
GPU – PyG | 7.7 $|$ 0.796 | 8.4 $|$ 0.720 | 10.9 $|$ 0.718
GPU – DGL | 13.3 $|$ 0.721 | 12.4 $|$ 0.641 | 12.5 $|$ 0.682

All experimental training runs were performed for 300 epochs on the same GNN
model structure. This model was not optimized for best training performance
but remained consistent for all scenarios so that direct comparisons focusing
on the graph frameworks, hardware, and parallelism approach could be observed
independent of the model structure.

###  6.1 Benchmarks

Table 3: Benchmarks on different compute architectures and graph frameworks for the sequential GAT model with the PubMed dataset. *The full graph was defined in the GNN model instead of being passed through as a subset of nodes to be re-built as a sub-graph. Framework | Compute | Epoch 1 (s) | Epochs 2–300 (s) | Ave. Epoch (s) | Train Loss | Train Acc. | Val Acc.
---|---|---|---|---|---|---|---
DGL | Single CPU | 0.3555 | 101.2 | 0.3386 | 0.2000 | 0.9833 | 0.7520
DGL | Single GPU | 0.2254 | 3.736 | 0.0125 | 0.2030 | 1.000 | 0.7520
PyG | Single CPU | 0.7946 | 238.7 | 0.7985 | 0.1567 | 0.9833 | 0.7910
PyG | Single GPU | 0.2509 | 3.260 | 0.0109 | 0.2131 | 1.000 | 0.7920
DGL | DGX with GPipe Chunk = 1* | 6.985 | 3.755 | 0.0126 | 0.1984 | 1.000 | 0.7660
PyG | DGX with GPipe Chunk = 1* | 7.312 | 3.407 | 0.0114 | 0.2097 | 1.000 | 0.7840
DGL | DGX with GPipe Chunk = 1 | 7.294 | 15.62 | 0.0522 | 0.1879 | 0.9500 | 0.7780
DGL | DGX with GPipe Chunk = 2 | 7.192 | 12.30 | 0.0411 | 0.4283 | 0.8333 | 0.6000
DGL | DGX with GPipe Chunk = 3 | 7.281 | 15.29 | 0.0511 | 0.5204 | 0.7667 | 0.4920
DGL | DGX with GPipe Chunk = 4 | 7.712 | 18.06 | 0.0604 | 0.6016 | 0.7500 | 0.4580

As a first comparison benchmark, we trained the GNN model on single devices
with the citation datasets of Cora, CiteSeer, and PubMed. The training time
and test accuracy results are summarized in Table 2. As expected, training
times, as measured by the average time per training epoch, on the GPU are
faster for both graph frameworks across all datasets. Interestingly, DGL
trained on average 35% faster than PyG on a CPU, while PyG trained on average
29% faster than DGL on a GPU. This outcome suggests, at least for our applied
GNN model training on a single device, PyG may be better optimized for a GPU
and DGL for a CPU. Training accuracy remained within a range of 15.5%, with
PyG averaging 2.4% better than DGL over all datasets.

Next, we compare the average training time per epoch on three compute
architectures, including the single CPU and GPU, as previously measured, with
the DGX system comprised of four GPUs leveraging GPipe pipeline parallelism
without micro-batching. In each case presented in Figure 1, the complete graph
data object was included in the graph convolution layer during each training
epoch. The comprehensive benchmark report in Table 3 implements the GAT model
on the PubMed dataset across combinations of frameworks and compute clusters
for an expanded analysis and comparison of the runtimes and accuracy. We also
include the duration of the first epoch in the reported training times to
provide a complete comparison of the graph frameworks and hardware, which
varied slightly across graph frameworks and architectures. The remaining
training epochs ran on the order of 80–100 times faster on the single GPU
compared to the single CPU.

Figure 1: Benchmark training times for DGL and PyG on the PubMed dataset
comparing the single devices to multiple devices with pipeline parallelism.
Here, data parallelism is disabled.

Surprisingly, no significant performance improvement in training time is
observed in the four GPU system using GPipe with a “chunk size” = 1 (i.e., no
micro-batching) compared to a single GPU. The PubMed dataset used in these
experiments is considered small compared to those that are intended to benefit
from pipeline parallelism. Therefore, the added cost of shifting data across
the four GPUs may overtake the minimal speedup provided by GPipe. This may
also suggest that the additional feature of data parallelism (via the data
“chunks”) provided by GPipe is crucial to realizing meaningful performance
improvements. We also measured the training accuracy resulting from both graph
machine learning frameworks applied with GPipe across four GPUs, but without
micro-batching, exactly as in the timing measurements. As plotted in Figure 2,
each framework converged similarly in accuracy over 300 training epochs in
this configuration.

Figure 2: Training accuracy with the DGL and PyG frameworks with pipe
parallelism across four GPUs with no graph data batching.

###  6.2 Increased training time

To investigate the impact of data parallelism within GPipe, we activated
micro-batching and ran the training with the DGL graph framework to compare
total training times between a single GPU and multiple distributed GPUs. As
seen in Figure 3, the training times dramatically increase with micro-batching
enabled at two, three, and four batches, as generated by GPipe.

Figure 3: Increased training time with GPipe applied pipeline parallelism with
an increasing number of graph micro-batches.

Our approach adapts the forward training to pass the graph information along
with the features through a tuple of tensors into the sequential model. The
GPipe micro-batching splits each tensor within the tuple so that only a subset
of nodes is passed through, along with its corresponding set of features, as
expected. The first convolution layer receives this subset of nodes indices
but still must have a complete graph structure as its input with the feature
tensor. So, a re-build of a graph is first performed with a DGL framework-
delivered method. This sub-graph creation from the provided subset of nodes
requires the full graph data object as a reference. However, DGL necessitates
that the full graph, `g`, must remain on the CPU. To generate the sub-graph
within the convolution layer, a copy of the subset node tensor must first be
moved from the GPU onto the CPU, then the sub-graph is built and moved back
onto the GPU. This data flow across devices was performed twice because our
model includes two convolution layers. So, significant overhead was added to
the total time just to enable the basic training calculations. As the chunk
size increased, more micro-batches were generated, resulting in even more sub-
graph build steps. Fortunately, the feature tensor extracted from the passed
tuple could remain on the GPU. However, the updated values were still re-
packaged into a tuple with the original sub-graph nodes to be returned into
the forward pass of the model sequence.

###  6.3 Degraded accuracy

We also observed that the training accuracy suffered severely with an
increasing number of micro-batches. Although GPipe micro-batching can be
disabled, as configured for our benchmark tests (Figures 1 and 2), the
expected benefit of pipeline parallelism requires micro-batching. We next ran
the same DGL framework-based model with GPipe across four GPUs, sequentially
distributed as before, to observe the effects of micro-batching in our adapted
implementation.

The intended design of the GPipe micro-batching through the torchgpipe library
implementation is to separate the features tensor into uniform batches. This
challenges our adaptation that passes a tuple containing both a node tensor
and feature tensor, as we observed the micro-batching being applied to each
tensor by sequentially selecting the tensor indices into a number of batches
equal to the set chunk size parameter. This sequential separation preserves
the nodes of the resulting sub-graph with their corresponding features.
However, the edge relationships between the nodes are lost. Although edges are
re-established during the sub-graph re-build in the convolution layers, the
original graph structure is not expected to store its edges sequentially.
Therefore, separating the graph this way during the micro-batching likely
eliminates crucial node relationships that need aggregated during the graph
convolution layer calculations. As expected from such a potential for
significant information loss during the GPipe micro-batching, as the number of
batches generated increases, the training accuracy drops, as seen in Figure 4.

Figure 4: Accuracy drop-off with GPipe and graph micro-batching with
comparisons to the previous training accuracy results without batching.

##  7 Conclusion and Future Work

By analyzing the performance of GNNs using pipeline parallelism via GPipe, our
first results suggest that although GPipe has demonstrated great success in
optimizing deep neural networks, its ability to deliver efficiency for a graph
neural network remains limited with the adapted implementation presented here.

An immediate scope for future work is to determine how to customize the GPipe
data parallelism to utilize intelligent graph batching instead of a sequential
separation by index. Such an improvement is expected to increase accuracy to
match the benchmark levels while benefiting from the parallelism in runtime
efficiency. The SIGN technique described in Section 3 may be the best batching
approach to consider for parallelizing GNNs with our implementation because it
avoids the experienced pitfalls of node and graph sampling and instead
provides precomputed node representations that may be straightforwardly mini-
batched by GPipe.

Pipeline parallelism is intended to benefit neural network training on very
large datasets, much greater in size than the PubMed set used here to
establish our adapted implementation without overburdening memory resources.
We anticipate that runtime performance will increase for training on extremely
large graphs after memory bottlenecks or computational complexity overtake the
capability offered by a single GPU. Extending the current implementation to
massive datasets, on the scale of millions of nodes and a billion edges, such
as the Reddit post dataset Hamilton et al. (2017), the Amazon data dump
McAuley et al. (2015), and others available with data loaders for DGL and PyG
through the Open Graph Benchmark (OGB) Hu et al. (2020), will better
illustrate the impact of GPipe parallelism on GNNs, and provide a deeper
understanding for potential enhancements to our current implementation for
improving training performance with these existing tools and frameworks.

## Acknowledgments

Thank you to Prof. Ian Foster, Prof. Rick Stevens, and Peng Ding for the
guidance and feedback on this paper. We also thank the DGX team, Daniel
Murphy-Olson and Ryan Aydelott, and the Computing, Environment, and Life
Sciences directorate at Argonne National Laboratory.

## References

  * Alon & Yahav (2020) Alon, U. and Yahav, E.  On the bottleneck of graph neural networks and its practical implications.  _arXiv preprint arXiv:2006.05205_ , 2020.
  * Arora (2020) Arora, S.  A survey on graph neural networks for knowledge graph completion, 2020\.
  * Auten et al. (2020) Auten, A., Tomei, M., and Kumar, R.  Hardware acceleration of graph neural networks.  In _2020 57th ACM/IEEE Design Automation Conference (DAC)_ , pp. 1–6. IEEE, 2020.
  * Bapst et al. (2020) Bapst, V., Keck, T., Grabska-Barwińska, A., Donner, C., Cubuk, E., Schoenholz, S., Obika, A., Nelson, A., Back, T., Hassabis, D., and Kohli, P.  Unveiling the predictive power of static structure in glassy systems.  _Nature Physics_ , 16(4):448–454, 2020.
  * Bruna et al. (2013) Bruna, J., Zaremba, W., Szlam, A., and LeCun, Y.  Spectral networks and locally connected networks on graphs.  _arXiv preprint arXiv:1312.6203_ , 2013.
  * Cheriton & Tarjan (1976) Cheriton, D. and Tarjan, R.  Finding minimum spanning trees.  _SIAM Journal on Computing_ , 5(4):724–742, 1976\.
  * Chiang et al. (2019) Chiang, W.-L., Liu, X., Si, S., Li, Y., Bengio, S., and Hsieh, C.-J.  Cluster-gcn: An efficient algorithm for training deep and large graph convolutional networks.  In _Proceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ , pp. 257–266, 2019.
  * Defferrard et al. (2016) Defferrard, M., Bresson, X., and Vandergheynst, P.  Convolutional neural networks on graphs with fast localized spectral filtering.  In _Advances in neural information processing systems_ , pp. 3844–3852, 2016.
  * Dijkstra (1959) Dijkstra, E.  A note on two problems in connexion with graphs.  _Numerische Mathematik_ , 1:269–271, 1959.
  * Duvenaud et al. (2015) Duvenaud, D. K., Maclaurin, D., Iparraguirre, J., Bombarell, R., Hirzel, T., Aspuru-Guzik, A., and Adams, R. P.  Convolutional networks on graphs for learning molecular fingerprints.  In _Advances in neural information processing systems_ , pp. 2224–2232, 2015.
  * Fey & Lenssen (2019) Fey, M. and Lenssen, J. E.  Fast graph representation learning with pytorch geometric.  _arXiv preprint arXiv:1903.02428_ , 2019.
  * Frasca et al. (2020) Frasca, F., Rossi, E., Eynard, D., Chamberlain, B., Bronstein, M., and Monti, F.  Sign: Scalable inception graph neural network.  _arXiv preprint arXiv:2004.11198_ , 2020.
  * Hamilton et al. (2017) Hamilton, W., Ying, Z., and Leskovec, J.  Inductive representation learning on large graphs.  In _Advances in neural information processing systems_ , pp. 1024–1034, 2017.
  * Henaff et al. (2015) Henaff, M., Bruna, J., and LeCun, Y.  Deep convolutional networks on graph-structured data.  _arXiv preprint arXiv:1506.05163_ , 2015.
  * Hu et al. (2020) Hu, W., Fey, M., Zitnik, M., Dong, Y., Ren, H., Liu, B., Catasta, M., and Leskovec, J.  Open graph benchmark: Datasets for machine learning on graphs.  _arXiv preprint arXiv:2005.00687_ , 2020.
  * Huang et al. (2019) Huang, Y., Cheng, Y., Bapna, A., Firat, O., Chen, D., Chen, M., Lee, H., Ngiam, J., Le, Q. V., Wu, Y., et al.  Gpipe: Efficient training of giant neural networks using pipeline parallelism.  In _Advances in neural information processing systems_ , pp. 103–112, 2019.
  * Kim et al. (2020) Kim, C., Lee, H., Jeong, M., Baek, W., Yoon, B., Kim, I., Lim, S., and Kim, S.  torchgpipe: On-the-fly pipeline parallelism for training giant models.  2020\.
  * Kipf & Welling (2016) Kipf, T. N. and Welling, M.  Semi-supervised classification with graph convolutional networks.  _arXiv preprint arXiv:1609.02907_ , 2016.
  * Li et al. (2015) Li, Y., Tarlow, D., Brockschmidt, M., and Zemel, R.  Gated graph sequence neural networks.  _arXiv preprint arXiv:1511.05493_ , 2015.
  * Liu et al. (2020) Liu, H., Lu, S., Chen, X., and He, B.  G3: when graph neural networks meet parallel graph processing systems on gpus.  _Proceedings of the VLDB Endowment_ , 13(12):2813–2816, 2020.
  * Ma et al. (2019) Ma, L., Yang, Z., Miao, Y., Xue, J., Wu, M., Zhou, L., and Dai, Y.  Neugraph: parallel deep neural network computation on large graphs.  In _2019 $\\{$USENIX$\\}$ Annual Technical Conference ($\\{$USENIX$\\}$$\\{$ATC$\\}$ 19)_, pp. 443–458, 2019.
  * McAuley et al. (2015) McAuley, J., Targett, C., Shi, Q., and Van Den Hengel, A.  Image-based recommendations on styles and substitutes.  In _Proceedings of the 38th international ACM SIGIR conference on research and development in information retrieval_ , pp. 43–52, 2015.
  * Mercado et al. (2020) Mercado, R., Rastemo, T., Lindelöf, E., Klambauer, G., Engkvist, O., Chen, H., and Bjerrum, E. J.  Graph Networks for Molecular Design.  8 2020.  doi: 10.26434/chemrxiv.12843137.v1.  URL https://chemrxiv.org/articles/preprint/Graph_Networks_for_Molecular_Design/12843137.
  * Sen et al. (2008) Sen, P., Namata, G., Bilgic, M., Getoor, L., Galligher, B., and Eliassi-Rad, T.  Collective classification in network data.  _AI magazine_ , 29(3):93–93, 2008.
  * Strokach et al. (2020) Strokach, A., Becerra, D., Corbi-Verge, C., Perez-Riba, A., and Kim, P.  Fast and flexible protein design using deep graph neural network.  _Cell Systems_ , 11(4):402–411.e4, 2020.
  * Vashishth (2019) Vashishth, S.  Neural graph embedding methods for natural language processing.  _CoRR_ , abs/1911.03042, 2019.  URL http://arxiv.org/abs/1911.03042.
  * Veličković et al. (2017) Veličković, P., Cucurull, G., Casanova, A., Romero, A., Lio, P., and Bengio, Y.  Graph attention networks.  _arXiv preprint arXiv:1710.10903_ , 2017.
  * Wang et al. (2019) Wang, M., Zheng, D., Ye, Z., Gan, Q., Li, M., Song, X., Zhou, J., Ma, C., Yu, L., Gai, Y., Xiao, T., He, T., Karypis, G., Li, J., and Zhang, Z.  Deep graph library: A graph-centric, highly-performant package for graph neural networks.  _arXiv preprint arXiv:1909.01315_ , 2019.
  * Yang et al. (2016) Yang, Z., Cohen, W., and Salakhudinov, R.  Revisiting semi-supervised learning with graph embeddings.  In _International conference on machine learning_ , pp. 40–48. PMLR, 2016.
  * Ying et al. (2018) Ying, R., He, R., Chen, K., Eksombatchai, P., Hamilton, W. L., and Leskovec, J.  Graph convolutional neural networks for web-scale recommender systems.  _CoRR_ , abs/1806.01973, 2018.  URL http://arxiv.org/abs/1806.01973.
  * Yu et al. (2017) Yu, B., Yin, H., and Zhu, Z.  Spatio-temporal graph convolutional neural network: A deep learning framework for traffic forecasting.  _CoRR_ , abs/1709.04875, 2017.  URL http://arxiv.org/abs/1709.04875.
  * Zeng et al. (2020) Zeng, H., Zhou, H., Srivastava, A., Kannan, R., and Prasanna, V.  Graphsaint: Graph sampling based inductive learning method.  _arXiv preprint arXiv:1907.04931_ , 2020.
  * Zhang et al. (2020) Zhang, Z., Leng, J., Ma, L., Miao, Y., Li, C., and Guo, M.  Architectural implications of graph neural networks.  _IEEE Computer Architecture Letters_ , 19(1):59–62, 2020.
  * Zhou et al. (2018) Zhou, J., Cui, G., Zhang, Z., Yang, C., Liu, Z., Wang, L., Li, C., and Sun, M.  Graph neural networks: A review of methods and applications.  _arXiv preprint arXiv:1812.08434_ , 2018.
  * Zhou et al. (2020) Zhou, J., Xu, Z., Rush, A. M., and Yu, M.  Automating botnet detection with graph neural networks.  _arXiv preprint arXiv:2003.06344_ , 2020.


# Measurement of the scaling slope of compressible magnetohydrodynamic
turbulence by synchrotron radiation statistics

Xue-Wen Zhang,1 Jian-Fu Zhang,1,2 Ru-Yue Wang 1 and Fu-Yuan Xiang1,2
1Department of Physics, Xiangtan University, Xiangtan 411105, China,
2Key Laboratory of Stars and Interstellar Medium, Xiangtan University,
Xiangtan 411105, China E-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Based on magnetohydrodynamic turbulence simulations, we generate synthetic
synchrotron observations to explore the scaling slope of the underlying MHD
turbulence. We propose the new $Q$-$U$ cross intensity $X$ and cross-
correlation intensity $Y$ to measure the spectral properties of magnetic
turbulence, together with statistics of the traditional synchrotron $I$ and
polarization $PI$ intensities. By exploring the statistical behavior of these
diagnostics, we find that the new statistics $X$ and $Y$ can extend the
inertial range of turbulence to improve measurement reliability. When focusing
on different Alfvénic and sonic turbulence regimes, our results show that the
diagnostics proposed in this paper not only reveal the spectral properties of
the magnetic turbulence but also gain insight into the individual plasma modes
of compressible MHD turbulence. The synergy of multiple statistical methods
can extract more reliable turbulence information from the huge amount of
observation data from the Low-Frequency Array for Radio astronomy and the
Square Kilometer Array.

###### keywords:

ISM: general – ISM: turbulence—magnetohydrodynamics (MHD) — methods: numerical
— polarization

††pubyear: 2023††pagerange: Measurement of the scaling slope of compressible
magnetohydrodynamic turbulence by synchrotron radiation
statistics–LABEL:lastpage

##  1 Introduction

The magnetized turbulent fluids in astrophysical environments can be usually
described by magnetohydrodynamic (MHD) turbulence theory, which plays a
critical role in many astrophysical processes such as star formation (Mac Low
& Klessen 2004), heat conduction (Narayan & Medvedev 2001), magnetic
reconnection (Beresnyak 2017), and acceleration of cosmic rays (Yan & Lazarian
2008; Zhang & Xiang 2021; Zhang et al. 2023). Therefore, studying the
properties of MHD turbulence helps to advance the theory of MHD turbulence and
to understand astrophysical processes associated with MHD turbulence.

Here, we briefly describe three significant advances made in earlier research
on turbulence. The first is about incompressible non-magnetized turbulence. By
using a self-similarity assumption of the turbulence cascade, Kolmogorov
(1941, henceforth K41) derived a power-law relation of $E(k)\sim k^{-5/3}$ in
the inertial range, which is called a classic Kolmogorov spectrum. The second
is about incompressible magnetized turbulence. Iroshnikov & Kraichnan (1963;
1965, henceforth IK65) obtained the power-law scaling of $E(k)\sim k^{-3/2}$
in the inertial range by introducing nonlinear energy cascade. Although IK65
advanced the K41 theory by considering the effect of magnetic fields, it
ignored a critical issue that the turbulence should be anisotropic in the
magnetized fluids (Montgomery & Turner 1981). The third is still about
incompressible magnetized turbulence but focuses on the anisotropy of MHD
turbulence due to turbulent magnetic fields. Focused on the nonlinear energy
cascade of incompressible strong MHD turbulence, Goldreich & Sridhar (1995,
hereafter GS95) provided the theoretical predictions on the power-law scaling
and anisotropic relationship, as described in more detail in Section 2.1.

At present, many numerical simulations have significantly increased our
knowledge on the scaling, anisotropy and compressibility of MHD turbulence
(e.g., Cho & Lazarian 2002; see textbook by Beresnyak & Lazarian 2019; and a
recent review by Beresnyak 2019 ). The properties obtained by the simulation
of MHD turbulence can understand the acceleration and propagation of cosmic
rays (Yan & Lazarian 2002). In particular, the turbulent reconnection model
proposed in Lazarian & Vishniac (1999, hereafter LV99), which provides a new
interpretation for GS95 theory from the perspective of eddies, has been
applied to various astrophysical environments such as gamma-ray bursts
(Lazarian et al. 2003; Zhang & Yan 2011), microquasars (de Gouveia dal Pino &
Lazarian 2005), active galactic nuclei (Kadowaki et al. 2015) and radio
galaxies (Brunetti & Lazarian 2016).

Due to the large scale of astrophysical system with a high Reynolds number
$R_{\rm e}>10^{10}$, it is challenging to simulate a realistic astrophysical
environment by direct numerical simulation. The currently available 3D MHD
simulations can achieve the case of the Reynolds number $R_{\rm e}\simeq
10^{5}$ (e.g., Beresnyak & Lazarian 2019). A distinctive feature is that the
realistic inertial range in astrophysical turbulence is much greater than that
revealed by numerical simulations. It is more effective to move away from
direct numerical simulations and develop statistical techniques using
observational data, to explore the properties of MHD turbulence.

When relativistic electrons motion in turbulent magnetic fields, they produce
synchrotron radiation fluctuations providing information on the magnetic
fields (Schnitzeler et al. 2007; Lazarian & Pogosyan 2012, 2016 henceforth
LP12 and LP16, Iacobelli et al. 2013; Van Eck et al. 2017; West et al. 2020;
Sun et al. 2022 ). Based on the modern understanding of MHD turbulence and
synchrotron radiation theory, LP12 explored the properties of MHD turbulence
by statistics of synchrotron total intensity fluctuations. They predicted that
synchrotron intensity fluctuations are anisotropic with a long extension along
the direction of magnetic fields. Using the ratio between quadrupole and
monopole components can determine the anisotropy of MHD turbulence, which is
sensitive to the compressibility of underlying turbulence. These theoretical
predictions on anisotropy have been confirmed successfully using numerical
simulations (Herron et al. 2016; Lee et al. 2019; Wang et al. 2022). These
studies are opening new avenues for exploring MHD turbulence using
observational data.

Moreover, LP16 proposed to recover the properties of MHD turbulence using
synchrotron polarization intensity fluctuations (including Faraday rotation
fluctuations). They developed two main techniques from the perspective of
analytical theory, i.e., polarization frequency analysis and polarization
spatial analysis. The former, using a variance of synchrotron polarization
intensity (or its derivative) as a function of the square of the wavelength,
was tested by Zhang et al. (2016). The latter, making use of spatial
correlations of synchrotron polarization intensity (or its derivative) at the
fixed wavelength as a function of the spatial separation $R$, was tested by
Lee et al. (2016) and Zhang et al. (2018). As confirmed, these two methods
obtain the scaling index of the underlying turbulence cascade in the inertial
range. Compared with synchrotron radiation, polarized radiation can reveal not
only information about magnetic fields in the plane of the sky but also that
parallel to the line of sight (LOS).

From the perspective of synthetic observations, numerical dissipation
inevitably limits the extension of the power-law range, and the greater the
inertial range, the higher the reliability of the measurement. However, from
an observational point of view, the scaling index measurement of MHD
turbulence is also limited by the telescope’s resolution and data noise. It is
necessary to synergize multiple techniques to reveal the properties of MHD
turbulence and enhance the reliability of the measurement results.

This paper aims to advance the study of the power-law scaling properties of
MHD turbulence. With the power spectrum (PS) and structure function (SF)
methods for synchrotron diagnostics, we propose two new statistical quantities
to explore the scaling slope properties of compressible MHD turbulence. The
paper is structured as follows. Section 2 describes theory aspects involving
the basic theory of MHD turbulence, synchrotron radiative process and
statistical methods. Section 3 introduces the setups of the numerical
simulation of MHD turbulence. Sections 4 and 5 present the numerical results.
In Sections 6 and 7, we provide our discussion and summary, respectively.

##  2 Theoretical description

###  2.1 MHD turbulence theory

GS95 theory is generally considered the basis for MHD turbulence. Note that
GS95 theory focused on incompressible strong MHD turbulence with Alfvénic Mach
number $M{\rm{}_{A}}=V_{\rm L}/V_{\rm A}\simeq 1$, where $V_{\rm L}$ is the
injection velocity at the injection scale $L_{\rm inj}$ and $V_{\rm A}$ is the
Alfvénic velocity. This theory combined the motions of the eddies
perpendicular to the magnetic field with those parallel to the magnetic field
by the critical balance condition $l_{\perp}/v_{\perp}=l_{\parallel}/V_{\rm
A}$, where $v_{\perp}$ is the velocity at the scale $l$, and the scales
$l_{\parallel}$, $l_{\perp}$ represent the parallel and perpendicular scales
of eddies, respectively. They found that the motions of eddies perpendicular
to the magnetic field have similar properties to Kolmogorov turbulence with
the spectrum of $E(k_{\perp})\propto{\epsilon}^{2/3}k_{\perp}^{-5/3}$, and the
velocity-scale relation of $v_{\perp}\propto(\epsilon l_{\perp})^{1/3}$, where
$k_{\perp}$ is the wave-vector component perpendicular to the magnetic field
and $\epsilon$ is the rate of energy cascade. According to the velocity-scale
relation and critical balance condition, they predicted an anisotropic
relationship of

$l_{\parallel}\sim V_{\rm A}{\epsilon^{-1/3}}l_{\perp}^{2/3},$ (1)

which delineates the dependencies between the perpendicular and parallel
scales of the eddies.

Later, the GS95 theory was generalized from the trans-Alfvénic turbulence to
sub-Alfvénic and super-Alfvénic ones, respectively (LV99; Lazarian 2006). For
the former, $M_{\rm A}<1$, that is, the turbulence drives with the injection
velocity $V_{\rm L}$ less than the Alfvénic velocity $V_{\rm A}$, LV99 found
that the turbulence cascade corresponds to two regimes. The first regime is a
weak turbulence cascade ranging from the injection scale $L_{\rm inj}$ to the
transition scale $L_{\rm tr}=L_{\rm inj}M_{\rm A}^{2}$. The second one is
strong turbulence from the transition scale $L_{\rm tr}$ to the dissipation
scale $L_{\rm diss}$, where the energy cascade perpendicular to the magnetic
field is analogous to the hydrodynamic Kolmogorov cascade. In this strong
turbulence regime, they derived the turbulence velocity as

$v_{\perp}\approx V_{\rm L}L_{\rm{inj}}^{-1/3}M_{\rm A}^{1/3}l_{\perp}^{1/3},$
(2)

and the anisotropic relation as

$l_{\rm\parallel}\approx L_{\rm inj}^{1/3}M_{\rm A}^{-4/3}l_{\perp}^{2/3}.$
(3)

When taking $M_{\rm A}=1$, the above equations will return to the relevant
expressions of GS95 theory.

As for the latter, $M_{\rm A}>1$, the MHD turbulence starting from the
injection scale $L_{\rm inj}$ is almost no constraint of the magnetic field
and has properties similar to those of hydrodynamic turbulence. With the
cascade of turbulence, it experiences a transition from hydrodynamic-like
turbulence to MHD one at the scale $L_{\rm A}=L_{\rm inj}M_{\rm A}^{-3}$.
However, from the scale $L_{\rm A}$ to $L_{\rm diss}$, the turbulence again
follows the characteristics of GS95 theory, having the velocity-scale relation
of

$v_{\rm\perp}\approx V_{\rm L}L_{\rm inj}^{-1/3}l_{\perp}^{1/3},$ (4)

and the anisotropy of

$l_{\parallel}\approx L_{\rm inj}^{1/3}M_{\rm A}^{-1}l_{\perp}^{2/3}.$ (5)

At present, the properties of compressible MHD turbulence have become an
important part of the modern understanding of MHD turbulence theory.
Compressible MHD turbulence can be decomposed into three modes, namely Alfvén,
slow and fast modes, as confirmed by numerical simulations (Cho & Lazarian
2002, 2003; Kowal & Lazarian 2010). Specifically, they found that Alfvén and
slow modes follow the GS95-type scaling law, namely $E(k_{\perp})\propto
k_{\perp}^{-5/3}$ and the scale-dependent anisotropy, while fast mode presents
the scaling law of $E(k_{\perp})\propto k_{\perp}^{-3/2}$ and the isotropy. In
addition, for compressible MHD turbulence, the Alfvén mode is incompressible,
while the slow and fast modes, called magnetosonic modes, are compressible.
111 When focusing on the compressible MHD turbulence as done in this work,
one, for the sake of simplicity, can call the slow mode as a compressible
mode. However, for the incompressible MHD turbulence with the plasma parameter
$\beta\gg 1$, the slow mode is a pseudo-Alfvén mode with a purely solenoidal
feature.

Despite the progress made in the development of MHD turbulence theory, there
are still a lot of controversial issues. For example, Maron & Goldreich (2001)
numerically studied the incompressible MHD turbulence and found a shallow
energy spectral index of $k^{-3/2}$ different from $k^{-5/3}$ given by GS95.
Subsequently, to explain this shallow index, Boldyrev (2006) proposed the
dynamic alignment model to modify the GS95 scaling index from $-5/3$ to
$-3/2$. Later, Beresnyak & Lazarian (2010) and Beresnyak (2014) thought that
the spectral index $-5/3$ cannot extend to the entire inertial range, but
deviate near the part of the injection scale (see also Beresnyak & Lazarian
2019 for the recent review). This can explain why the low-resolution numerical
simulations generate a shallower spectral index, while the results of higher-
resolution numerical simulations are consistent with GS95. However, some
recent studies, e.g., Chandran et al. (2015), agreed with the dynamical
alignment theory.

By analyzing the power spectrum of super-sonic turbulence from the solenoidal
and compressive driving ways, Federrath (2013) found the velocity spectral
indices satisfy with $k^{-2}$. In the case of solenoidal driving, the spectrum
of the density-weighted velocity $\rho^{1/3}v$ satisfies with $k^{-1.74}$,
while in the case of compressive driving, the slope is significantly steeper
and close to $k^{-2.1}$. This result is consistent with the compressible
turbulence theory (Galtier & Banerjee 2011), which predicts the scaling of
density-weighted velocity $k^{-19/9}$. Recently, Mallet & Schekochihin (2017)
proposed the intermittency model to modify MHD turbulence theory at scales
close to the dissipation scales. However, because of the limitations of
numerical simulations, it is difficult to confirm.

Until now, many attempts have not significantly changed the framework of the
GS95 theory. Although our study below is based on the GS95 theory, the change,
in theory, does not affect our results based on synthetic synchrotron
observations.

###  2.2 Synchrotron emission fluctuations

For the sake of simplicity, this work assumes that relativistic electrons
interacting with the turbulent magnetic field satisfy a homogeneous and power-
law energy distribution of $N(E)=N_{0}E^{-p}$, where $p$ and $E$ represent the
spectral index and energy of relativistic electrons, respectively. Here,
$N_{0}$ is the normalization constant of electrons. According to the classic
textbooks (Rybicki & Lightman 1979; Longair 2011), observable Stokes
parameters under the condition of no Faraday rotation effect can be expressed
as follows (see also, e.g.,Waelkens et al. 2009; LP16):

${I}({\bm{X}})=\int_{0}^{L}dz(B_{\rm x}^{2}({\bm{x}})+B_{\rm
y}^{2}({\bm{x}}))^{\frac{p-3}{4}}(B_{\rm x}^{2}({\bm{x}})+B_{\rm
y}^{2}({\bm{x}})),$ (6) $Q_{0}({\bm{X}})=\int_{0}^{L}dz(B_{\rm
x}^{2}({\bm{x}})+B_{\rm y}^{2}({\bm{x}}))^{\frac{p-3}{4}}(B_{\rm
x}^{2}({\bm{x}})-B_{\rm y}^{2}({\bm{x}})),$ (7)
$U_{0}({\bm{X}})=\int_{0}^{L}dz(B_{\rm x}^{2}({\bm{x}})+B_{\rm
y}^{2}({\bm{x}}))^{\frac{p-3}{4}}(2B_{\rm x}({\bm{x}})B_{\rm y}({\bm{x}})),$
(8)

where $L$ is the integral depth along the LOS, $B_{\rm x}$ and $B_{\rm y}$ the
components of the magnetic field perpendicular to the LOS, and
${\bm{X}}=({x},{y})$ the spatial coordinate in the plane of the sky.

Focusing on linear polarization synchrotron radiation, we have a complex
vector

${\bm{P}}({\bm{X}},\lambda^{2})=Q+iU=\int_{0}^{L}d{z}P_{\rm
in}({\bm{X}},z)e^{2i{\rm\phi}({\bm{X}},z)},$ (9)

describing polarization states in the plane of the sky. In this equation, the
exponential factor involves Faraday rotation effect. The observed polarization
angle $\phi$ is expressed by

$\phi=\phi_{0}+{\lambda^{2}}\rm RM,$ (10)

where the angle $\phi_{0}$ is the intrinsic angle. The Faraday rotation
measure $\rm RM$ is written as

${\rm RM}(\bm{X},z)=0.81\int_{0}^{z}n_{\rm
e}({\bm{X}},z^{\prime})B_{\parallel}({\bm{X}},z^{\prime})dz^{\prime}~{}{\rm
rad}~{}{\rm m^{-2}},$ (11)

where $n_{\rm e}$ is the density of thermal electrons, $B_{\parallel}$ the
component of the magnetic field along the LOS. The integral length $L$ is
along the LOS from the position of the source at $z$ to the observer.
Moreover, the part $P_{\rm in}$ of integrated function in Equation (9)
represents the intrinsic polarization intensity density and can be expressed
by $P_{\rm in}\equiv(Q_{0},U_{0})$. After including Faraday rotation effect,
the new Stokes parameters $Q$ and $U$ can be rewritten as

$Q({\bm{X}},{\rm\lambda^{2}})=Q_{\rm 0}{\rm\cos 2\phi}+U_{\rm 0}{\rm\sin
2\phi},$ (12) $U({\bm{X}},{\rm\lambda^{2}})=U_{\rm 0}{\rm\cos 2\phi}-Q_{\rm
0}{\rm\sin 2\phi},$ (13)

from which we obtain the synchrotron polarization intensity of

$PI=\sqrt{Q^{2}+U^{2}}.$ (14)

A complete description of synchrotron radiation can be encoded by a
polarization matrix, e.g., as done in Equation (E1) of LP12. This paper
focuses on the correlation statistics between $Q$ and $U$. The first is the
$Q$-$U$ cross intensity defined by 222 This definition was used to explore the
anisotropy of MHD turbulence by the structure function (LP12) and to trace
magnetic field directions by gradient techniques (Lazarian & Yuen 2018).

$X^{2}=QU,$ (15)

which is related to the relative importance of $Q$ and $U$. In general, both
the different turbulence properties and the level of Faraday rotation
depolarization will lead to different $Q$ and $U$ values, and different ratios
of $Q$ and $U$. As done in LP16, when considering the correlation function of
the polarization complex vector ${\bm{P}}$, we have

$\langle P(\bm{X_{1}})P^{*}(\bm{X_{2}})\rangle=\langle
Q(\bm{X_{1}})Q(\bm{X_{2}})+U(\bm{X_{1}})U(\bm{X_{2}})\rangle\\\ +i\langle
U(\bm{X_{1}})Q(\bm{X_{2}})-Q(\bm{X_{1}})U(\bm{X_{2}})\rangle,$ (16)

which is split into real and imaginary parts that are separately invariant
with respect to frame rotation. The symmetric real part carries the most
straightforward information about the magnetized turbulent medium and has been
numerically studied in Zhang et al. (2016). LP16 predicted that the
antisymmetric imaginary part reflects helical correlations of the magnetic
field, which still needs numerical testing. In addition, analytical studies
demonstrated that the anisotropy of the MHD turbulence can generate the
observable antisymmetric correlations. Based on the antisymmetric imaginary
part in Equation (16), we rewrite the cross-correlation intensity as

$Y^{2}(\bm{X}^{\prime})=\int{d^{2}{\bm{X}}[U({\bm{X}})Q({\bm{X}}+\bm{X}^{\prime})-Q({\bm{X}})U({\bm{X}}+\bm{X}^{\prime})]}.$
(17)

Adopting Equations (15) and (17), we will explore the scaling property of MHD
turbulence by comparing the traditional $PI$ and $I$ statistics. Note that the
cross-correlation intensity $Y$ is covariant variable during rotation and
translation transformation in the Stokes frame, while the cross-intensity $X$
are unchanged only when the Stokes frame is translated.

###  2.3 Statistical methods

Although turbulence is a complex and chaotic process, it allows us to use
statistical methods to reveal its underlying properties. In this paper, we
focus on the SF and PS methods. We first consider the simplest and often used
correlation function for an arbitrary 2D physical quantity $\zeta$. According
to the textbook by Monin & Yaglom (1975), correlation and structure functions
are written as

${\rm CF}({\bm{R}})=\langle\zeta{({\bm{X}}+{\bm{R}})}\zeta({\bm{X}})\rangle,$
(18) ${\rm
SF}({\bm{R}})=\langle(\zeta{({\bm{X}}+{\bm{R}})}-\zeta({\bm{X}}))^{2}\rangle,$
(19)

and they satisfy the following relation

${\rm SF}({\bm{R}})=2[{\rm CF(0)-CF}(\bm{R})],$ (20)

where ${\bm{R}}$ is a separation vector, and $\langle...\rangle$ represents
the average through the whole volume space.

The PS, a common statistical tool in the study of turbulence, can provide
information on the energy cascade of MHD turbulence, such as the spectral
shape and index, the source and the sink. The PS of a two-dimensional physical
quantity is expressed by

$P_{\rm
2D}({\bm{K}})=\frac{1}{(2\pi)^{2}}\mathrm{\int}\langle\zeta({\bm{X}})\zeta({\bm{X}}+{\bm{R}})\rangle
e^{{-i{\bm{K}}}\cdot{\bm{R}}}d{\bm{R}}$ (21)

by the Fourier transform of the correlation function. The ring-integrated 1D
spectrum for a 2D variable follows

$E_{\rm{2D}}({K})=\int_{K-0.5}^{K+0.5}P_{\rm{2D}}(K)d{K}.$ (22)

Note that there is a direct connection between PS and SF by the scaling slope:
$E_{\rm 2D}(K)\propto K^{-m}$ and ${\rm SF}(R)\propto R^{m-1}$ (LP12; see also
numerical confirmation in Lee et al. 2016; Zhang et al. 2018), where $m$ is
equal to 8/3 for Kolmogorov power spectrum in two dimensions.

Table 1: The information of data cubes arising from the simulation of compressible MHD turbulence. Relevant parameters used to characterize data cubes are given as —- $B_{0}$: mean magnetic field along the $x$ coordinate; $\beta$: plasma parameter; $L_{\rm tr}$: transition scale of strong turbulence in sub-Alfvénic regime; $L_{\rm A}$: transition scale of strong turbulence in the super-Alfvénic regime. Run | $B_{0}$ | $M_{\rm A}$ | $M_{\rm s}$ | $\beta=2M^{2}_{\rm A}/M_{\rm s}^{2}$ |  $L_{\rm inj}$[2<$k$<3] | $L_{\rm inj}[k=2.5]$ |  $L_{\rm tr}$($L_{\rm A}$)[2<$k$<3] |  $L_{\rm tr}$($L_{\rm A}$)[$k$ = 2.5]
---|---|---|---|---|---|---|---|---
1 | 1.00 | 0.70 | 0.87 | 1.30 | [170.6, 256.0] | 204.8 | [83.59, 125.44] | 100.35
2 | 1.00 | 0.55 | 4.46 | 0.03 | [170.6, 256.0] | 204.8 | [51.61, 77.44] | 61.95
3 | 1.00 | 0.65 | 0.48 | 3.67 | [170.6, 256.0] | 204.8 | [72.09, 108.16] | 86.53
4 | 0.10 | 1.69 | 3.11 | 0.60 | [170.6, 256.0] | 204.8 | [35.34, 53.04] | 42.43
5 | 0.10 | 1.72 | 0.45 | 29.30 | [170.6, 256.0] | 204.8 | [33.53, 50.31] | 40.25

##  3 MHD turbulence simulations

To generate synchrotron observations, we use a third-order accurate hybrid,
essentially non-oscillatory (ENO) scheme (Cho & Lazarian 2002) to solve ideal
isothermal MHD equations in a periodic box of size $2\pi$:

$\frac{\partial\rho}{\partial t}+\nabla\cdot(\rho{\bm{v}})=0,$ (23)
$\rho[\frac{\partial{\bm{v}}}{\partial
t}+({\bm{v}}\cdot\nabla){\bm{v}}]+\nabla
p-{\bm{J}}\times\frac{\bm{B}}{4\pi}={\bm{f}},$ (24)
$\frac{\partial{\bm{B}}}{\partial t}-\nabla\times({\bm{v}}\times{\bm{B}})=0,$
(25) $\nabla\cdot{\bm{B}}=0,$ (26)

where $\rho$ is the density, $p=c^{2}_{s}\rho$ the thermal gas pressure,
$\bm{J}$ the current density, and $\bm{f}$ an external driving force. The
turbulence is driven by a solenoidal driving force at the wave number of
$2<k<3$ (corresponding to the mean wavenumber of $k\approx 2.5$) (Cho &
Lazarian 2003). By setting the initial mean magnetic field (along the $x$
axis) and the gas pressure, we run several simulations with the resolution of
$512^{3}$ covering different MHD turbulence regimes. When the running reaches
a statistically steady state, from output data cubes of the magnetic field,
velocity, and density, we calculate the Alfvénic and sonic Mach numbers to
characterize the properties of each simulation. Specifically, Alfvénic Mach
number is obtained by $M_{\rm A}=\langle\bm{V}_{\rm L}/\bm{V}_{\rm A}\rangle$
and sonic Mach number by $M_{\rm s}=\langle\bm{V}_{\rm L}/c_{\rm s}\rangle$,
where $c_{\rm s}=\sqrt{p/\rho}$ is sound speed, and $V_{\rm
A}\approx\frac{\bm{B}}{\sqrt{\rho}}$ is Alfvénic velocity. The resulting
numerical values are listed in Table 1.

The compressible MHD turbulence can be decomposed into Alfvén, slow and fast
modes by using the following theoretical procedures (Cho & Lazarian 2002; Cho
& Lazarian 2003)

$\hat{\xi}_{\rm
s}\propto(1+\frac{\beta}{2}-\sqrt{D})(k_{\perp}{\hat{\bm{k}}}_{\perp})+(-1+\frac{\beta}{2}-\sqrt{D})(k_{\parallel}{\hat{\bm{k}}}_{\parallel}),$
(27) $\hat{\xi}_{\rm
f}\propto(1+\frac{\beta}{2}+\sqrt{D})(k_{\perp}{\hat{\bm{k}}}_{\perp})+(-1+\frac{\beta}{2}+\sqrt{D})(k_{\parallel}{\hat{\bm{k}}}_{\parallel}),$
(28) $\hat{\xi}_{\rm
A}\propto-{\hat{\bm{k}}}_{\perp}\times{\hat{\bm{k}}}_{\parallel},$ (29)

in the Fourier space, where $D=(1+\frac{\beta}{2})^{2}-2\beta\cos^{2}\theta$,
and $\cos\theta={\hat{\bm{k}}}_{\parallel}\cdot{\hat{\bm{B}}}_{0}$. By
projecting the magnetic field into the displacement vectors $\hat{\xi}_{\rm
f}$, $\hat{\xi}_{\rm A}$ and $\hat{\xi}_{\rm s}$, we get the individual
components of three modes for magnetic field and then convert them to volume
space by the Fourier inverse transform. Later, Kowal & Lazarian (2010)
proposed an optimized decomposition method by introducing a discrete wavelet
transform, and decomposing each component of the magnetic fields into
orthogonal wavelets using a discrete wavelet transform. This method, which
depends on the local magnetic field rather than the mean magnetic field, is
universal for decomposition in the super-Alfvénic turbulence with weak
magnetic fields. Since only sub-Alfénic turbulence is involved in our work
when decomposing plasma modes, we use the Fourier decomposition method.

When exploring below the influence of the angle between the mean magnetic
field and LOS on the power spectrum, we adopt the Euler rotation algorithm to
rotate data cubes (Lazarian & Yuen 2018; Carmo et al. 2020; Wang et al. 2022;
Malik et al. 2023). The components of the rotation matrix
$\bm{\hat{F}}=\hat{\bm{F}}_{x}\hat{\bm{F}}_{y}\hat{\bm{F}}_{z}$ are expressed
as follows:

${\bm{\hat{F}_{x}}}=\left[\begin{array}[]{cccc}1&0&0&\\\ 0&\rm
cos(\varphi_{x})&-\rm sin(\varphi_{x})&\\\ 0&\rm sin(\varphi_{x})&\rm
cos(\varphi_{x})&\end{array}\right],$ (30)
${\hat{\bm{F}}_{y}}=\left[\begin{array}[]{cccc}\rm cos(\varphi_{y})&0&\rm
sin(\varphi_{y})&\\\ 0&1&0&\\\ -\rm sin(\varphi_{y})&0&\rm
cos(\varphi_{y})&\\\ \end{array}\right],$ (31)
${\hat{\bm{F}}_{z}}=\left[\begin{array}[]{cccc}\rm cos(\varphi_{z})&-\rm
sin(\varphi_{z})&0&\\\ \rm sin(\varphi_{z})&\rm cos(\varphi_{z})&0&\\\
0&0&1&\\\ \end{array}\right],$ (32)

where $\varphi_{m=x,y,z}$ is the rotation angle along the $x$, $y$, $z$ axis,
respectively. For data cube $\Re(\bm{r})$ of components of magnetic filed, the
rotated data cube is obtained by
${\hat{\bm{F}}\Re({\bm{\hat{}}{\bm{F}}}^{-1}}{\bm{r}})$ transformation. Since
the rotation of the cube is equivalent to the rotation of the observation
frame in the opposite direction, we perform an inverse transformation of the
position vector ${\bm{r}}$.

Figure 1: Images of synchrotron radiation diagnostics: $Q$-$U$ cross ($X$),
cross-correlation ($Y$), linear polarization $(PI)$ and total ($I$)
intensities. The calculations without the Faraday rotation effect are based on
run3 listed in Table 1. The mean magnetic field is along the horizontal
direction.

##  4 Statistical Results: Compressible MHD Turbulence

When generating synthetic synchrotron observations, we adopt dimensionless
units. In the case of involving the Faraday rotation effect, we quantify
relevant physical quantities: magnetic fields as $B=1.23\rm\ \mu G$, thermal
electron density as $n_{\rm e}=0.01\rm\ cm^{-3}$, and the box size as
$L=100\rm\ pc$, which correspond to the Galactic ISM environment. In addition,
we set the spectral index of relativistic electrons to be $p=2.5$.333 Electron
spectral indices are associated with a specific acceleration mechanism. For
instance, de Gouveia dal Pino & Lazarian (2005) predicted that turbulent
reconnection acceleration provided a steeper spectral index $p=2.5$, confirmed
numerically by Zhang et al. (2023). The spectral index determined from
observations will change in different astrophysical environments. Our previous
studies demonstrated that the change in spectral indices could not impede the
turbulence measurement using the synchrotron polarization technique (see LP12
for theoretical predictions; Lee et al. 2016 and Zhang et al. 2018 for
numerical confirmations).

###  4.1 Slope measurement of MHD turbulence without Faraday rotation

Before performing statistical analysis, let us illustrate the map structures
of the statistical diagnostics considered in this paper. Based on run3 listed
in Table 1, we obtain the intensities of different synchrotron radiation
diagnostics, namely $Q$-$U$ cross $X$, cross-correlation $Y$, linear
polarization $PI$ and total $I$ intensities. Here, we do not consider the
Faraday rotation effect. The imaging of these diagnostics is plotted in Figure
1 (using the real part of $X$, $Y$ to exhibit their maps), from which we can
see that map structures of the $Q$-$U$ cross $X$, cross-correlation $Y$ are
elongated along the direction perpendicular to the mean magnetic field while
maps of linear polarization $PI$ and total $I$ intensities have nearly similar
structures extending along the horizontal direction, i.e., the direction of
the mean magnetic field. The perpendicular distribution from $X$ and $Y$
should be dominated by the Stokes parameter $U$, while the horizontal
structure from $PI$ and $I$ is dominated by the Stokes parameter $Q$. In
general, map structures in the Stokes parameter $Q$ are aligned with the
direction of the mean magnetic field, while those of $U$ are perpendicular to
the mean magnetic field. Moreover, the intensities of $PI$ and $I$ diagnostics
have larger amplitudes than those of $X$ and $Y$. This is caused by amplitude
changes in $Q$ and $U$, the intensity of which is associated with the magnetic
field strength in the plane of the sky (see Equations (7) and (8)).

The SFs of $X$, $Y$, $PI$, and $I$ are plotted in Figure 2 for four different
turbulence regimes: sub-Alfvénic and supersonic (left upper panel), sub-
Alfvénic and subsonic (right upper), super-Alfvénic and supersonic (left
lower) and super-Alfvénic and subsonic (right lower). As shown, SFs cannot
recover the scaling of MHD turbulence in the regime ranging from the injection
scale $L_{\rm inj}$ to the transition scale $L_{\rm tr}$ for $M_{\rm A}<1$ (or
$L_{\rm A}$ for $M_{\rm A}>1$). These numerical results are in agreement with
theoretical predictions of MHD turbulence cascade due to weak turbulent
interaction (LV99; Lazarian 2006). At the scale less than the transition scale
$L_{\rm tr}$ (or $L_{\rm A}$), i.e., in the strong turbulence regime, these
four diagnostics present the power-law distributions predicted by LV99 and
Lazarian (2006). From the figures, We can see that: (1) in the case of sub-
Alfvénic turbulence (upper panels), the measurements from $X$ and $Y$ are
closer to the slope index 5/3 than those from $PI$ and $I$; and (2) in the
case of super-Alfvénic turbulence scenario (lower panels), the $Y$ statistics
can better determine the slope index 5/3 compared with the other three
statistics $X$, $PI$ and $I$. Comparing sub-Alfvénic and super-Alfvénic
turbulence, we find that super-Alfvénic turbulence has a shorter inertial
range for $X$, $PI$, and $I$. Interestingly, we find that statistics $Y$ can
well reflect the scaling of 5/3 with a wide inertial range and does not depend
on specific turbulence properties.

Figure 2:  Structure function of synchrotron radiation diagnostics: $Q$-$U$
cross ($X$), cross-correlation ($Y$), linear polarization $(PI)$ and total
($I$) intensities in different turbulence regimes. The yellow and green
vertical dashed lines represent the injection and transition scales,
respectively. Figure 3: Power spectra of the synchrotron radiation
diagnostics: $Q$-$U$ cross ($X$), cross-correlation ($Y$), linear polarization
$(PI)$ and total ($I$) intensities in different turbulence regimes. The yellow
and green vertical dashed lines represent the injection and transition scales,
respectively.

Based on data cubes used in Figure 2, we plot the PS of $X$, $Y$, $PI$ and $I$
in Figure 3. As seen, the scaling indices of PS of $X$ and $Y$ are consistent
with $-8/3$ in four turbulence regimes, with an extended inertial range for
MHD turbulence. This is a valuable finding in this paper. Due to the presence
of numerical dissipation at large wavenumber, the measurements of $PI$ and $I$
can only provide a narrow power-law range, which limits their flexibility to
determine the scaling slope of the underlying MHD turbulence. Therefore, the
new statistics $X$ and $Y$ have advantages over traditional statistics $PI$
and $I$ in measuring the spectral index and inertial range of MHD turbulence.

In addition, the amplitudes of PS of $X$, $Y$, $PI$, and $I$ in the sub-
Alfvénic turbulence regime are greater than those in the super-Alfvénic
turbulence regime. The reason is that large mean magnetic fields for sub-
Alfvénic turbulence produce more synchrotron radiative information, enhancing
Stokes parameters $Q$ and $U$ intensities. In the case of sub-Alfvénic
turbulence (upper panels), the amplitudes of $X$ are greater than those of
$PI$ and $I$ in the inertial range, while it is the opposite in the case of
super-Alfvénic turbulence (lower panels). Our studies on PS demonstrate that
the four statistics explored can measure the scaling slope of MHD turbulence.
We emphasize their synergistic measurement abilities to enhance reliability.
At the same time, by comparing the amplitudes of different quantities, we can
understand their magnetic field strength. It provides a new way to measure
magnetization strength, $M_{\rm A}$; further research is necessary in the
future.

###  4.2 Slope measurement of MHD turbulence with Faraday rotation

####  4.2.1 Effect of radiative frequency

In this section, we explore how the radiation frequency and the angle between
the mean magnetic field and the LOS influence the PS of $X$, $Y$, and $PI$
($I$ independent of the frequency) in the presence of the Faraday rotation
effect. To explore the influence of radiation frequency, we first set the mean
magnetic field parallel and perpendicular to the LOS, respectively. Based on
the run3 listed in Table 1, we show the numerical results in Figure 4 for the
mean magnetic field parallel (left column) and perpendicular (right column) to
the LOS. As is shown in the left column, the PS of $X$, $Y$, and $PI$ follow
the scaling law of $-8/3$ in the scales of $<L_{\rm tr}$ for simulations at
high frequency (about $\nu\geq 0.1\ \rm GHz$), while in the case of low
frequencies, they downward (upward) deviate from $-8/3$ at the small (large)
wavenumber regions. This is because, in the high-frequency range, the effect
of the Faraday rotation depolarization on the PS is small. With decreasing the
frequency, the strong Faraday rotation depolarization leads to a weaker
correlation of the radiation signal. As the frequency decreases, the
appearance of noise gradually fills the entire synthetical map of Stokes
parameters $Q$ and $U$, severely downward distorting the PS statistics at
large scales (small wavenumbers). The increase of noise at small scales leads
to the upward deviation of the spectral distribution.

The right column of Figure 4 shows the results of the PS at different
frequencies, for which we stress that the mean magnetic field is in the
direction perpendicular to the LOS. As seen, the PS of three synchrotron
diagnostics $X$, $Y$, and $PI$ satisfy the power law of $-8/3$ at the higher
frequencies (about $\nu\geq 0.1\ \rm GHz$) compared with the left panels. At
the low-frequency regimes, they show a distribution similar to those of the
left column. In addition, we see that the amplitudes of PS of $X$, $Y$, and
$PI$ have significant changes at the lower frequencies. In addition to the
case of $M_{\rm A}<1$ studied above, we also consider other possible scenarios
with $M_{\rm A}>1$ (the relevant results not shown in this paper). When the
frequency is higher than 100$~{}\rm MHz$, the PS of $X$, $Y$, and $PI$ can
reveal the spectral index of the underlying MHD turbulence.

As a result, in the case of moderate depolarization, $X$ and $Y$ have more
advantages than $PI$ for measuring the scaling slope of MHD turbulence. When
the mean magnetic field is parallel to the LOS, it is more helpful to use
these statistics to reveal the magnetic turbulence information in the case of
a lower frequency (down to about $\nu\simeq 0.1\ \rm GHz$ for our parameter
selections).

####  4.2.2 Effect of noise and view angle

We here explore how the angle between the mean magnetic field and the LOS
affects the distribution of PS from relevant diagnostics in Figure 5, plotted
at the frequency of $0.14~{}\rm GHz$. As is seen in this figure, the inertial
range and the amplitude are decreased with decreasing the angle for $X$ and
$PI$, resulting in the measured spectral index deviating from -8/3 at a large
scale. The reason is that with decreasing angle, the Faraday rotation measure
makes synchrotron-polarized signal depolarization. For $Y$, only its amplitude
changes rather than its inertial range.

Figure 4: Power spectra of synchrotron radiation diagnostics: $Q$-$U$ cross
($X$), cross-correlation ($Y$) and linear polarization $(PI)$ intensities at
different frequencies. Left column: the mean magnetic field is along the LOS.
Right column: the mean magnetic field is perpendicular to the LOS. The yellow
and green vertical dashed lines are plotted to represent the injection and
transition scales, respectively. Our calculations are based on the run3 listed
in Table 1.  Figure 5: Power spectra of synchrotron radiation diagnostics:
$Q$-$U$ cross ($X$), cross-correlation ($Y$) and linear polarization $(PI)$
intensities at different angles between the mean magnetic field and the LOS at
the frequency $\nu$ = 0.14$~{}\rm GHz$. The yellow and green vertical dashed
lines are the injection and transition scales, respectively. Our calculations
are based on the run3 listed in Table 1.  Figure 6: The influence of the noise
on the power spectra of synchrotron radiation diagnostics: $Q$-$U$ cross
($X$), cross-correlation ($Y$) and linear polarization $(PI)$ calculated at
the frequency $\nu$ = 0.14$\rm~{}GHz$. The symbol $\sigma$ indicates a
standard deviation of Gaussian noise and accounts for a fraction of the mean
synchrotron intensities. The yellow and green vertical dashed lines denote the
injection and transition scales, respectively.

Figure 6 explores the influence of the noise on the scaling index of PS of
$X$, $Y$, and $PI$ at the frequency of $0.14~{}\rm GHz$. We generate a
Gaussian noise map with the resolution of $512^{2}$ and add it to the original
image to study the noise effect. Here, the standard deviation of Gaussian
noise accounts for the fraction of the mean synchrotron intensity. The figure
clearly shows that the PS of $X$, $Y$, and $PI$ adding the Gaussian noise
deviate upward from those without noise in the large-$K$ regime. This is
because adding noise increases the random fluctuation of the original image,
resulting in an increase in the PS in the small-scale region. In addition, we
can see that under the same noise level, the inertial range measured by $X$ is
wider than those of $PI$ and $Y$. And the higher the level of Gaussian noise,
the more obvious the deviation. Therefore, increasing the level of Gaussian
noise makes the power-law inertial range narrower.

###  4.3 The influence of numerical resolution on the results

To test the influence of numerical resolution on the measurement of the
scaling slope of turbulence, we simulate the data cube with the low numerical
resolution of $256^{3}$ in the same way as the run1 listed in Table 1. The
difference of resolution results in slightly different Mach numbers: $M_{\rm
A}=0.70$ and $M_{\rm s}=0.74$ for $256^{3}$, as well as $M_{\rm A}=0.70$ and
$M_{\rm s}=0.87$ for $512^{3}$.

Figure 7: Structure functions of synchrotron radiation diagnostics: $Q$-$U$
cross ($X$), cross-correlation ($Y$), linear polarization $(PI)$ and total
($I$) intensities. The results are presented in two numerical resolutions
$256^{3}$ (blue) and $512^{3}$ (red), with the run1 listed in Table 1. The
vertical lines represent the injection and transition scales, respectively.

Still focusing on the SFs of $X$, $Y$, $PI$ and $I$, we provide numerical
results in Figure 7. From this figure, we find that: (1) the measurements from
high resolution $512^{3}$ show better results, namely, closer to $R^{5/3}$,
than those from low resolution $256^{3}$, as expected; (2) For the data cubes
with the same numerical resolution, the structure functions of $X$ and $Y$ are
more advantageous than those of $PI$ and $I$ in measuring the scaling index
and inertial range of MHD turbulence.

In addition, we here explore the influence of numerical resolution on the PS
of synchrotron radiation diagnostics using two data cubes with the resolutions
$256^{3}$ and $512^{3}$, the parameters of which correspond to those of data
cubes used in Figure 7. The numerical results are presented in Figure 8, from
which we see that the resolution $512^{3}$ can better determine the spectral
index of $-8/3$ expected in the inertial range. Importantly, we find that the
measurements of $Q$-$U$ cross $X$, and cross-correlation $Y$ extend the width
of the inertial range, compared with the traditional statistics of linear
polarization $PI$ and total $I$ intensities.

Comparing the upper and lower panels of Figure 8, we find that the
measurements of $PI$ and $I$ have a dissipation at large wavenumbers, i.e.,
small scales, while $X$ and $Y$ show a weak dissipation to extend the inertial
range. Therefore, when recovering the scaling slope of MHD turbulence from
observational data, we recommend $X$ and $Y$ statistics. However, we should be
particularly cautious that in the largest wavenumber, there is an effect of
numerical noise. In this regard, interested readers are advised to refer to
Zhang et al. (2016) and Zhang et al. (2018) who found that the noise would
cause the spectrum to reverse upwards.

Figure 8: Power spectra of synchrotron radiation diagnostics: $Q$-$U$ cross
($X$), cross-correlation ($Y$), linear polarization $(PI)$ and total ($I$)
intensities. The vertical lines denote the injection and transition scales of
the underlying MHD turbulence. Numerical results are plotted in resolutions of
$256^{3}$ and $512^{3}$, using run1 listed in Table 1.

##  5 Statistical Results: Decomposition of Compressible MHD turbulence

With the procedures described in Section 3, we decompose data cubes of
compressible MHD turbulence, namely, run3 listed in Table 1, and then explore
the PS of $X$, $Y$, $PI$, and $I$ for three modes based on the post
decomposition data cubes. As is shown in Figure 9, the PS of $X$, $Y$, $PI$,
and $I$ follows the power law index of $-8/3$ for Alfvén and slow modes, and
$-5/2$ for fast mode. Meanwhile, it can be seen that the PS of new statistics
$X$ and $Y$ has a more extended inertial range than that of traditional $PI$
and $I$ statistics. The properties of PS obtained for these plasma modes are
consistent with the theoretical prediction of compressible MHD turbulence
described in Section 2.

Figure 9: Power spectra of synchrotron radiation diagnostics: $Q$-$U$ cross
($X$), cross-correlation ($Y$), linear polarization $(PI)$ and total ($I$)
intensities for Alfvén, slow and fast modes. The yellow and green vertical
dashed lines represent the injection and transition scales, respectively. The
decomposition of data cubes is based on the run3 listed in Table 1.

We further explore how the frequency influences the PS of different
diagnostics arising from three modes, as shown in Figure 10. Each row
corresponds to the PS of the same diagnostic for Alfvén (left), slow (middle),
and fast (right) modes, respectively, while in each column, the PS of
different diagnostics $X$ (upper), $Y$ (middle) and $PI$ (lower) for the same
mode. The PS of three diagnostics ($X$, $Y$ and $PI$) follows the scaling
slope of $K^{-8/3}$ for Alfvén (left column) and slow modes (middle column),
and $K^{-5/2}$ for fast mode (right column), at the frequency $\nu>0.05$ GHz,
while they deviate from the expected values ($-8/3$ or $-5/2$) at lower
frequencies, particularly, for small wavenumbers. It can be seen that the
amplitudes of PS of three diagnostics go up with decreasing frequency in the
large-$K$ part. For the slow and fast modes, the physical interpretation of
the above phenomena is similar to that of Section 4.2.1. But for the Alfvén
mode, the Stokes parameters of pure Alfvén mode at $\varphi=90^{\circ}$
projects quicker than random walk (Lazarian et al. 2022) in the case of low
$M_{\rm A}$. Physically it means the Alfvén mode without Faraday rotation
self-projects to zero if the integration length is large enough. FR destroys
this phenomenon and creates fluctuations that are not canceling itself.

Significantly, we find that the PS of various diagnostics arising from three
modes depends on the frequency. Note that the PS of different statistics for
the slow mode has a slightly weaker dependence on the frequency than the other
two modes. The reason is that the slow mode has a small Faraday rotation
measure value.

Our research by simulation data confirms that scaling slopes of compressible
plasma modes can be obtained from observational data. In practice, using
observational data to extract the properties of three plasma modes is a
challenging subject. Since the scaling index of $-5/2$ for the fast mode is
different from the $-8/3$ for the Alfvén and slow modes, one could first
obtain the properties of the fast mode (Zhang et al. 2020a). However, the more
challenging is how to effectively distinguish Alfvén and slow modes from the
observational data. This deserves more exploratory efforts.

Figure 10: Power spectra of synchrotron radiation diagnostics at different
frequencies. In the order of the rows: $Q$-$U$ cross $X$ (first row), cross-
correlation $Y$ (second), and linear polarization $PI$ (third). In the order
of the columns: Alfvén (first column), slow (second), and fast (third) modes.
The yellow and green vertical dashed lines represent the injection and
transition scales, respectively. The decomposition of data cubes is based on
the run3 listed in Table 1.

##  6 Discussion

Synchrotron radiation is an important source of magnetic field information in
the interstellar medium environment. Statistics of the total intensity $I$ can
provide the information of magnetic field in the plane of the sky. Compared
with the intensity $I$, the linear polarization intensity $PI$ can reflect
more information about the magnetic field, i.e., not only in the plane of the
sky but also in the LOS. The disadvantage is that $P$ can only provide the
total polarization information but not the relative importance of $Q$ and $U$,
i.e., the relative changes in $Q$ and $U$ values. In this paper, we propose
two new diagnostics $X$ and $Y$ together with the traditionally used $PI$ and
$I$ to explore the scaling properties of MHD turbulence by PS and SFs and
found that PS of $X$ and $Y$ have a larger measurable inertial range compared
with $PI$ and $I$. In fact, each technique has its advantages and limitations
when measuring turbulence properties. A synergy of various techniques can
obtain more comprehensive turbulence information, enhancing the reliability of
the turbulence measurement.

Although the PS of synchrotron diagnostics can not provide more information on
the spatial structure of MHD turbulence, it is an advantageous statistical
method for studying the source, sink and scaling slope of turbulence. As
studied in Section 4.2.1, we found that the PS has different amplitudes when
the LOS is perpendicular and parallel to the mean magnetic field, so we expect
that it can also be an alternative tool for measuring magnetization.

For compressible MHD turbulence, one can also explore density fluctuation
information within MHD turbulence by introducing Faraday rotation. However,
the biggest difficulty in measuring magnetic fields through involvement in
Faraday rotation studies is that there is not currently a good way to decouple
the coupling between vector and scalar. With a proper understanding of the
magnetic field through synergistically related techniques, one can gain
insight into the density information. Various MHD turbulence modes have
important effects on many astrophysical processes. Therefore, the study of
plasma modes is helpful to understand the contribution of different modes to
these physical processes, such as the acceleration and diffusion process of
cosmic rays (Zhang & Xiang 2021; Sampson et al. 2023).

This paper focused on the scaling properties of magnetic fields by PS and SF
statistics. Notice that the latter can also be used to recover properties of
magnetic field structure and eddy, as done in Wang et al. (2020) and Zhang et
al. (2020b). To understand other aspects of MHD turbulence, many other
synergistic techniques have been developed based on synchrotron radiation.
These techniques include the kurtosis and skewness exploring the anisotropy of
MHD turbulence (Herron et al. 2016) and constraining the sonic Mach number
(Burkhart et al. 2012), the quadrupole moment revealing the anisotropy (Herron
et al. 2016; Lee et al. 2019; Wang et al. 2022), as well as the gradient
statistics measuring magnetic field directions (Lazarian et al. 2017; Lazarian
& Yuen 2018; Zhang et al. 2019a; Zhang et al. 2019b, 2020b; Wang et al. 2021;
Liu et al. 2023) and magnetization (Carmo et al. 2020; Lazarian et al. 2022).
In addition, the PS of the tension force can diagnose the spatial structure of
the magnetic structures (Schekochihin et al. 2004; Waelkens et al. 2009; Sun
et al. 2014).

For completeness, our work explored how to get the scaling index of the
turbulence by a synchrotron signal in the cases of both subsonic and
supersonic turbulence. Indeed, the hot/warm ionized diffuse media with the low
$M_{\rm s}$ can be probed by radio synchrotron emission (such as the Galactic
ISM with $M_{\rm s}\leq 2$, see Gaensler et al. 2011), while some environments
still have a large $M_{\rm s}$, such as the regions of active galactic nuclei
and supernova remnants interacting with the surrounding cold molecular cloud.
Therefore, the $M_{\rm s}$ we explored in this paper were not much greater
than 1. For turbulence regimes much larger than 1, one can use alternative
approaches, such as velocity channel analysis and velocity correlation
spectrum (e.g., Lazarian & Pogosyan 2004; Yuen & Lazarian 2017; Yang et al.
2021). In this work, we did not involve the effect of self-absorption. This
process will become important when the magnetic field interacts with
relativistic electrons at low-frequency regimes. In the presence of self-
absorption, the PS of these statistics may vary not only in the scaling index
but also in the inertial range, which provides us with a new research
perspective to recover the 3D magnetic field structure.

##  7 Summary

In this paper, we proposed two new synchrotron diagnostics: the cross
intensity $X$ and cross-correlation intensity $Y$ to reveal the MHD turbulence
properties. Using their PS and SF together with traditional diagnostics $PI$
and $I$, we have well understood the spectral properties of the underlying
compressible MHD turbulence. We focused on exploring how Mach numbers, noise,
Faraday depolarization, and numerical resolution affect the spectral
measurement of magnetic turbulence. The main results are summarized as
follows.

  * •

The SF of statistics $X$, $Y$, $PI$, and $I$ can determine the scaling slope
of MHD turbulence in sub-Alfvénic regimes. Interestingly, new statistics $Y$
could better measure the scaling slope compared with other statistics $X$,
$PI$, and $I$ in the different Alfvénic regimes.

  * •

The noise does not impede the recovery of the scaling index of MHD turbulence,
and the inertial range of PS measured by $X$ is wider than that by $PI$ and
$Y$ at the same noise level.

  * •

In the case of moderate Faraday depolarization, they still improve the scaling
slope measurements since the statistics $X$ and $Y$ extend the inertial range.
The influence of numerical resolution does not change our conclusions.

  * •

The change of angle between the mean magnetic field and the LOS does not
affect the measurement of the scaling index, but the inertial range and
amplitude.

  * •

Using the synchrotron radiation diagnostics ($X$, $Y$ and $PI$) can measure
the spectral properties of Alfvén, slow and fast modes.

## ACKNOWLEDGMENTS

We thank the anonymous referee for valuable comments that significantly
improved the quality of the paper. J.F.Z. thanks to the support from the
National Natural Science Foundation of China (grant Nos. 11973035), the Hunan
Province Innovation Platform and Talent Plan-HuXiang Youth Talent Project (No.
2020RC3045), and the Hunan Natural Science Foundation for Distinguished Young
Scholars (No. 2023JJ10039). F.Y.X. acknowledges the support from the Joint
Research Funds in Astronomy U2031114 under a cooperative agreement between the
National Natural Science Foundation of China and the Chinese Academy of
Sciences.

## DATA AVAILABILITY

The data underlying this paper can be shared on reasonable request to the
corresponding author.

## References

  * Beresnyak (2014) Beresnyak A., 2014, ApJ, 784, L20
  * Beresnyak (2017) Beresnyak A., 2017, ApJ, 834, 47
  * Beresnyak (2019) Beresnyak A., 2019, Living Reviews in Computational Astrophysics, 5, 2
  * Beresnyak & Lazarian (2010) Beresnyak A., Lazarian A., 2010, ApJ, 722, L110
  * Beresnyak & Lazarian (2019) Beresnyak A., Lazarian A., 2019, Turbulence in Magnetohydrodynamics, Studies in Mathematical Physics, 12, De Gruyter
  * Boldyrev (2006) Boldyrev S., 2006, Phys. Rev. Lett., 96, 115002
  * Brunetti & Lazarian (2016) Brunetti G., Lazarian A., 2016, MNRAS, 458, 2584
  * Burkhart et al. (2012) Burkhart B., Lazarian A., Gaensler B. M., 2012, ApJ, 749, 145
  * Carmo et al. (2020) Carmo L., et al., 2020, ApJ, 905, 130
  * Chandran et al. (2015) Chandran B. D. G., Schekochihin A. A., Mallet A., 2015, ApJ, 807, 39
  * Cho & Lazarian (2002) Cho J., Lazarian A., 2002, Phys. Rev. Lett., 88, 245001
  * Cho & Lazarian (2003) Cho J., Lazarian A., 2003, MNRAS, 345, 325
  * Federrath (2013) Federrath C., 2013, MNRAS, 436, 1245
  * Gaensler et al. (2011) Gaensler B. M., et al., 2011, Nature, 478, 214
  * Galtier & Banerjee (2011) Galtier S., Banerjee S., 2011, Phys. Rev. Lett., 107, 134501
  * Goldreich & Sridhar (1995) Goldreich P., Sridhar S., 1995, ApJ, 438, 763
  * Herron et al. (2016) Herron C. A., Burkhart B., Lazarian A., Gaensler B. M., McClure-Griffiths N. M., 2016, ApJ, 822, 13
  * Iacobelli et al. (2013) Iacobelli M., et al., 2013, A&A, 558, A72
  * Iroshnikov (1963) Iroshnikov P. S., 1963, Azh, 40, 742
  * Kadowaki et al. (2015) Kadowaki L. H. S., de Gouveia Dal Pino E. M., Singh C. B., 2015, ApJ, 802, 113
  * Kolmogorov (1941) Kolmogorov A., 1941, Akademiia Nauk SSSR Doklady, 30, 301
  * Kowal & Lazarian (2010) Kowal G., Lazarian A., 2010, ApJ, 720, 742
  * Kraichnan (1965) Kraichnan R. H., 1965, Physics of Fluids, 8, 1385
  * Lazarian (2006) Lazarian A., 2006, ApJ, 645, L25
  * Lazarian & Pogosyan (2004) Lazarian A., Pogosyan D., 2004, ApJ, 616, 943
  * Lazarian & Pogosyan (2012) Lazarian A., Pogosyan D., 2012, ApJ, 747, 5
  * Lazarian & Pogosyan (2016) Lazarian A., Pogosyan D., 2016, ApJ, 818, 178
  * Lazarian & Vishniac (1999) Lazarian A., Vishniac E. T., 1999, ApJ, 517, 700
  * Lazarian & Yuen (2018) Lazarian A., Yuen K. H., 2018, ApJ, 865, 59
  * Lazarian et al. (2003) Lazarian A., Petrosian V., Yan H., Cho J., 2003, arXiv e-prints, pp astro–ph/0301181
  * Lazarian et al. (2017) Lazarian A., Yuen K. H., Lee H., Cho J., 2017, ApJ, 842, 30
  * Lazarian et al. (2022) Lazarian A., Yuen K. H., Pogosyan D., 2022, ApJ, 935, 77
  * Lee et al. (2016) Lee H., Lazarian A., Cho J., 2016, ApJ, 831, 77
  * Lee et al. (2019) Lee H., Cho J., Lazarian A., 2019, ApJ, 877, 108
  * Liu et al. (2023) Liu M., Hu Y., Lazarian A., Xu S., Soida M., 2023, MNRAS, 519, 1068
  * Longair (2011) Longair M. S., 2011, High Energy Astrophysics. Cambridge University Press, Cambridge
  * Mac Low & Klessen (2004) Mac Low M.-M., Klessen R. S., 2004, Reviews of Modern Physics, 76, 125
  * Malik et al. (2023) Malik S., Yuen K. H., Yan H., 2023, arXiv e-prints, p. arXiv:2303.17282
  * Mallet & Schekochihin (2017) Mallet A., Schekochihin A. A., 2017, MNRAS, 466, 3918
  * Maron & Goldreich (2001) Maron J., Goldreich P., 2001, ApJ, 554, 1175
  * Monin & Yaglom (1975) Monin A. S., Yaglom A. M., 1975, Statistical Fluid Mechanics: Mechanics of Turbulence, Vol. 2. The MIT Press, Cambridge, Massachusetts
  * Montgomery & Turner (1981) Montgomery D., Turner L., 1981, Physics of Fluids, 24, 825
  * Narayan & Medvedev (2001) Narayan R., Medvedev M. V., 2001, ApJ, 562, L129
  * Rybicki & Lightman (1979) Rybicki G. B., Lightman A. P., 1979, Radiative processes in astrophysics. Wiley, New York
  * Sampson et al. (2023) Sampson M. L., Beattie J. R., Krumholz M. R., Crocker R. M., Federrath C., Seta A., 2023, MNRAS, 519, 1503
  * Schekochihin et al. (2004) Schekochihin A. A., Cowley S. C., Taylor S. F., Maron J. L., McWilliams J. C., 2004, ApJ, 612, 276
  * Schnitzeler et al. (2007) Schnitzeler D. H. F. M., Katgert P., Haverkorn M., de Bruyn A. G., 2007, A&A, 461, 963
  * Sun et al. (2014) Sun X. H., Gaensler B. M., Carretti E., Purcell C. R., Staveley-Smith L., Bernardi G., Haverkorn M., 2014, MNRAS, 437, 2936
  * Sun et al. (2022) Sun X.-H., Gao X.-Y., Reich W., Jiang P., Li D., Yan H., Li X.-H., 2022, Research in Astronomy and Astrophysics, 22, 125011
  * Van Eck et al. (2017) Van Eck C. L., et al., 2017, A&A, 597, A98
  * Waelkens et al. (2009) Waelkens A. H., Schekochihin A. A., Enßlin T. A., 2009, MNRAS, 398, 1970
  * Wang et al. (2020) Wang R.-Y., Zhang J.-F., Xiang F.-Y., 2020, ApJ, 890, 70
  * Wang et al. (2021) Wang R.-Y., Zhang J.-F., Lazarian A., Xiao H.-P., Xiang F.-Y., 2021, MNRAS, 505, 6206
  * Wang et al. (2022) Wang R.-Y., Zhang J.-F., Lazarian A., Xiao H.-P., Xiang F.-Y., 2022, ApJ, 940, 158
  * West et al. (2020) West J. L., Henriksen R. N., Ferrière K., Woodfinden A., Jaffe T., Gaensler B. M., Irwin J. A., 2020, MNRAS, 499, 3673
  * Yan & Lazarian (2002) Yan H., Lazarian A., 2002, Phys. Rev. Lett., 89, 281102
  * Yan & Lazarian (2008) Yan H., Lazarian A., 2008, ApJ, 673, 942
  * Yang et al. (2021) Yang B., Zhang J.-F., Lazarian A., de Medeiros J. R., 2021, MNRAS, 503, 768
  * Yuen & Lazarian (2017) Yuen K. H., Lazarian A., 2017, ApJ, 837, L24
  * Zhang & Xiang (2021) Zhang J.-F., Xiang F.-Y., 2021, ApJ, 922, 209
  * Zhang & Yan (2011) Zhang B., Yan H., 2011, ApJ, 726, 90
  * Zhang et al. (2016) Zhang J.-F., Lazarian A., Lee H., Cho J., 2016, ApJ, 825, 154
  * Zhang et al. (2018) Zhang J.-F., Lazarian A., Xiang F.-Y., 2018, ApJ, 863, 197
  * Zhang et al. (2019a) Zhang J.-F., Lazarian A., Ho K. W., Yuen K. H., Yang B., Hu Y., 2019a, MNRAS, 486, 4813
  * Zhang et al. (2019b) Zhang J.-F., Liu Q., Lazarian A., 2019b, ApJ, 886, 63
  * Zhang et al. (2020a) Zhang H., Chepurnov A., Yan H., Makwana K., Santos-Lima R., Appleby S., 2020a, Nature Astronomy, 4, 1001
  * Zhang et al. (2020b) Zhang J.-F., Hu K., Cho J., Lazarian A., 2020b, ApJ, 895, 20
  * Zhang et al. (2023) Zhang J.-F., Xu S., Lazarian A., Kowal G., 2023, Journal of High Energy Astrophysics, submitted
  * de Gouveia dal Pino & Lazarian (2005) de Gouveia dal Pino E. M., Lazarian A., 2005, A&A, 441, 845