# GABOR FILTER INCORPORATED CNN
FOR COMPRESSION

Akihiro Imamura  Embedded AI R&D department
DENSO CORPORATION
Shinagawa, Tokyo, Japan
<EMAIL_ADDRESS>Nana Arizumi  Embedded AI R&D department
DENSO CORPORATION
Shinagawa, Tokyo, Japan
<EMAIL_ADDRESS>
###### Abstract

Convolutional neural networks (CNNs) are remarkably successful in many
computer vision tasks. However, the high cost of inference is problematic for
embedded and real-time systems, so there are many studies on compressing the
networks. On the other hand, recent advances in self-attention models showed
that convolution filters are preferable to self-attention in the earlier
layers, which indicates that stronger inductive biases are better in the
earlier layers. As shown in convolutional filters, strong biases can train
specific filters and construct unnecessarily filters to zero. This is
analogous to classical image processing tasks, where choosing the suitable
filters makes a compact dictionary to represent features. We follow this idea
and incorporate Gabor filters in the earlier layers of CNNs for compression.
The parameters of Gabor filters are learned through backpropagation, so the
features are restricted to Gabor filters. We show that the first layer of
VGG-16 for CIFAR-10 has 192 kernels/features, but learning Gabor filters
requires an average of 29.4 kernels. Also, using Gabor filters, an average of
83% and 94% of kernels in the first and the second layer, respectively, can be
removed on the altered ResNet-20, where the first five layers are exchanged
with two layers of larger kernels for CIFAR-10.

###### Index Terms:

compression, CNN, Gabor

††publicationid: pubid: 978-1-6654-0645-1/21/$31.00 ©2021 IEEE

##  I Introduction

In recent years, convolutional neural networks (CNNs) have achieved remarkable
success in computer vision tasks, e.g. image classification [1], object
detection [2], semantic segmentation [3]. The networks have grown deeper [4,
5], and wider [6] to achieve higher accuracy in these tasks. Hence, the
networks are increasing to have more parameters with a high cost of inference
operations. It is a major problem for mobile devices and embedded systems
applications, which has a limitation in memory, power, and computation speed.
Also, it is an issue for real-time systems like image classification for web
services, which also has a time budget to process inferences.

Hence, there has been a significant amount of works on model compression, and
acceleration [7, 8]. Among which, network pruning from well-trained CNN has
been widely studied [9, 10, 11]. There are different granularities of pruning.
Fine-grained level (e.g., weight level) sparsity gives the highest flexibility
and has high compression. However, it also gives unstructured sparsity and
needs specialized hardware or/and sparse libraries to harvest the best
efficiency from the compression [11, 12, 13]. On the contrary, structured
pruning does not require them. Hence, there are increasingly more works on
channel-level pruning [14, 15, 16] and kernel-level pruning [17]. Note that
the channel-level pruning in one layer also makes kernel-level pruning in the
next layer. Layer-level pruning is also possible for deep CNNs [18, 19]. In
this work, we also focus on pruning without introducing unstructured sparsity.

The reason why we can prune the weights is that after training, the filters
become small enough to be ignored; however, we do not know which to be pruned
without training the network [20]. Hence, we can assume that unnecessary
filters become zero by training with weight decay. However, the filters might
still have redundancy. The success of compressing the network using low-rank
decomposition of the convolutional layers [21, 22] indicates that features
required to represent the layer are relatively small.

The recent advances in self-attention models showed that convolution filters
are preferable to self-attention filters in the earlier layers [23, 24, 25],
which indicates that stronger inductive biases are better in the earlier
layers. So, additional inductive biases might help choose the right filters in
the earlier layers. Since it has been known that the networks are trained to
become the Gabor filters in the earlier layers of CNNs [1], we can insert
additional biases by enforcing filters to be the Gabor filters.

We observed that the first layer of VGG-16 on CIFAR-10 requires an average of
13.2 $7\times 7\times 3$ channels of Gabor kernels without accuracy
degradation. On the contrary, the original VGG-16 has 64 $3\times 3\times 3$
channels in the first layer. Hence, in the first two layers, it makes 24%
computation and memory reduction. Even more, an average of 162.6 $7\times 7$
kernels out of 192 can be removed without accuracy degradation using kernel-
level pruning. Similar behavior is observed for ResNet-20 on CIFAR-10. We
extended this idea to the second layer and experimented with altered
ResNet-20, where the first five layers are exchanged with two layers of
$7\times 7$ and $5\times 5$ without shortcut connection. It shows that 83% and
94% of kernels can be removed in the first and second layers, respectively,
without accuracy degradation.

Using well-trained CNNs, we can fit and retrain the features by restricting to
Gabor filters in the earlier layers. This way, we can insert stronger
inductive biases in the earlier layers and help to choose more specific
features. As with classical feature detection, clear feature prediction can
significantly reduce the size of CNNs. Finding better features and reduce
redundant channels is an excellent way to compress CNNs, because it does not
create unstructured sparse pruning. Furthermore, any other compressing methods
can be applied after our channel-level pruning since it does not change the
structure of CNNs.

##  II Related Work

The Gabor filters are the classical method that can extract local features
[26], and it is particularly popular because of invariance in scale, rotation,
and translation. Also, the linear convolution of the Gabor filters is a good
representation of the simple cells in the human visual cortex [27, 28], which
is a good indication that the early stage of feature extraction of the human
visual system also uses Gabor filters. Analogously, the earlier layers of CNNs
are similar to Gabor filters [29, 1]. Hence, there are several studies on
utilizing Gabor filters with CNNs.

Gabor filters can be used as pre-processing of the input data [30, 31, 32], or
initialization [33]. [34] uses fixed Gabor filters for the first or the second
layer to reduce training complexity. [35] incorporated selected Gabor
orientation filters to each layer to improve the robustness of the networks.
[36] uses Gabor filters for the interpretability of the networks in applying
person re-identification. They train Gabor filters’ parameters but do not
learn the scale or the center of the filters since their interests are on the
features’ interpretability and not to remove the features. In our work, we
also used backpropagation to learn the parameters of Gabor filters in the
network to enforce the layer to be the suitable feature detection.

##  III methodology

In recent CNNs, kernel sizes are mostly $1\times 1$ and $3\times 3$. It is
mostly because of computational efficiency. However, the first layer can have
larger filters because the input images are usually one or three channels, and
the computational cost in channel size is much smaller. In fact, there are
several popular networks with larger first filters [1, 5]. We could expand the
kernels in the first layer for the other CNNs with a smaller kernel. Expanding
kernels itself should not change the accuracy drastically. Even though the
kernel size becomes larger, if the number of channels compressed, then the
total computation cost can be smaller than the original CNNs. Also, a larger
kernel size means it covers larger receptive fields, so if the first few
layers are only to expand receptive fields, then the layers can also be
pruned. It is more common to use larger kernels in earlier layers for larger
images.

Once the first layer’s kernel is expanded and trained, we can fit the Gabor
filters to the well-trained network. Then, train Gabor filter parameters using
backpropagation. This way, we can enforce learned features to be Gabor
filters. After retraining, prune channels using the L1 norm criterion as in
[37].

Figure 1: Schematic of pruning methods

Fig. 1 illustrates the schematic of this method in yellow. We also show the
other methods for comparison. Here, ”retrain” does not mean resuming the
learning. We want to compare each method, so we retrain the network with the
same hyperparameters, including learning rates. Even though the network
training is converged in the pretraining stage, it is known that the changing
learning rate increases accuracy [38].

###  III-A Gabor Fitting

The Gabor filter [39] is defined by a sinusoidal wave multiplied by the
Gaussian function. We use the real part

$G=a\,\exp(-\frac{\hat{x}^{2}+\gamma\hat{y}^{2}}{2\sigma^{2}})\cos(2\pi\frac{\hat{x}}{\lambda}+\psi),$
(1)

where

$\displaystyle\hat{x}$ $\displaystyle=$
$\displaystyle(x-x_{0})\cos(\theta)+(y-y_{0})\sin(\theta),$ (2)
$\displaystyle\hat{y}$ $\displaystyle=$
$\displaystyle-(x-x_{0})\sin(\theta)+(y-y_{0})\cos(\theta),$ (3)

$\lambda$ is the wavelength of the sinusoidal factor, $\theta$ is the
orientation of the normal to the parallel stripes, $\psi$ is the phase offset,
$\sigma$ is the standard deviation of the Gaussian envelope, $\gamma$ is the
spatial aspect ratio, $a$ is the amplitude, and $(x_{0},y_{0})$ is the center.
Consider the kernels are located at the $(x,y)$ grids from $(1,1)$ to $(k,k)$,
where $k$ is the size of the kernel.

Construct Gabor filters using the following parameters and choose the closest
to the original kernel using the L2 norm.

  * •

amplitude $a$ to be $-1,-0.5,0,0.5,1$

  * •

center $(x_{0},y_{0})$ to be all the combination between 1 to size of kernel

  * •

orientation $\theta$ and phase offset $\psi$ to be $0,\pi/4,\pi/2,3\pi/4$

  * •

standard deviation of the Gaussian envelope $\sigma$ and wavelength $\lambda$
to be $1,2,3,4,5$

  * •

aspect ratio $\gamma$ to be $0.2,0.4,0.6,0.8,1$

###  III-B Gabor Learning

Using backpropagation, we can train the parameters of Gabor filters,
$\lambda$, $\theta$, $\psi$, $\sigma$, $\gamma$, $a$, $x_{0}$, and $y_{0}$.
The details of backpropagation are written in [36]. Contrary to their method,
we did not set a range of values of parameters. Also, we add the amplitude and
the center of the filters for learnable parameters. We intend to prune the
filters, so amplitude should dissipate, and we can remove filters altogether.
We did not use regularization loss of the Gabor parameter since we did not
observe the difference. We used the Autograd package in PyTorch [40] to
implement. For our experiments, the same learning rate is applied to all
parameters of Gabor filters.

##  IV Experiments

In CNNs, input feature maps $X_{i}$ are convolved with channels $C_{i}$ for
$n_{i+1}$ times and produce output feature maps $X_{i+1}$ at the $i$th layer.
The channels $C_{i}$ has $n_{i}$ kernels of the size $k_{i}\times k_{i}$. We
use shorthand notation of $k_{i}\times k_{i},n_{i}$ for the channels. In our
experiments, all input images are 3 channels.

We used CIFAR-10 [41] image database, which has $32\times 32$ natural images.
Note, $7\times 7$ kernels can cover a large part of the image. There are 5000
training images and 1000 test images. Standard data augmentation and
normalization of input data are also applied.

We conducted experiments, as shown in Fig. 1 on two networks, VGG [42] and
ResNet [5]. Both of which are widely popular networks and commonly used for
pruning experiments. We followed VGG-16 from [37] and ResNet-20 from [43] for
pre-processing, hyper-parameters, and learning schemes. It has been known that
the network training is sensitive to initialization and order of training
[44]. So, we ran all the experiments 5 times and took an average. We pruned
channels and kernels that do not reduce accuracy by more than 0.2%. The
channels and the kernels are chosen from the minimum according to the L1-norm.
The testing for each case is operated with the same 79 batches with 128
images.

###  IV-A Gabor Fitting of first layer

For both VGG-16 and ResNet-20, we expanded 3 $\times$ 3 kernels of the first
layer to $7\times 7$. As for ResNet-20, we remove two layers to fit the
receptive field since it only has 0.5% degradation of accuracy. We did not do
so for VGG-16, because it causes more than 2% of accuracy loss. First, we
train VGG-16 for 350 Epochs twice with the same hyperparameters and ResNet-20
for 170 Epochs once, then fit the first layer to the Gabor filters. After
retraining using conventional convolution and Gabor learning, we compare the
accuracy as shown in Tab. I. Fitting does not degrade the accuracy.

TABLE I: Accuracy of Gabor Learning first layer | $3\times 3$ | $7\times 7$
---|---|---
VGG-16
pre-train | 93.37 $\pm$ 0.09 | 93.23 $\pm$ 0.13
without fit + retrain | 93.44 $\pm$ 0.18 | 93.30 $\pm$ 0.13
fit + retrain | - | 93.39 $\pm$ 0.08
fit + Gabor learning | - | 93.11 $\pm$ 0.08
ResNet-20
pre-train | 91.25 $\pm$ 0.14 | 90.80 $\pm$ 0.07
without fit + retrain | 91.89 $\pm$ 0.09 | 91.48 $\pm$ 0.10
fit + retrain | - | 91.41 $\pm$ 0.08
fit + Gabor learning | - | 91.42 $\pm$ 0.10

The percentile of prunable kernels and channels are shown on Tab. II. For
VGG-16, the percentile of channel pruning is the same as the results with
$3\times 3$ kernels without fitting presented in [37]. Note that for
ResNet-20, two layers are also removed, but that is not shown on the table.

TABLE II: First layer pruning | VGG-16 | ResNet-20
---|---|---
| channel | kernel | channel | kernel
retrain | 49% | 55% | 14% | 21%
Gabor learning | 79% | 91% | 50% | 75%
Figure 2: Examples of kernels in the first layer channels of ResNet-20 for
both conventional convolutional learning and Gabor learning after Gabor fit.
Inside the 3 $\times$ 3 block, the first row shows after the Gabor fit, the
second row shows the Gabor learning after Gabor fit, and the third row shows
the conventional convolutional learning after Gabor fit. Each column shows
filters for RGB images.

To show how the kernels were learned, Fig. 2 shows examples of kernels in the
first layer channels of ResNet-20 for both conventional convolutional learning
and Gabor learning after Gabor fit. Each kernel is normalized, and the
colormap is shown in 6 $\sigma$. For each block with 3 $\times$ 3 kernels, the
first row shows after the Gabor fit, the second row shows the Gabor learning
after Gabor fit, and the third row shows the conventional convolutional
learning after Gabor fit. Both of the learnings are starting from the first
row. The columns are corresponding to the 3 input channels of RGB images. It
clearly shows that conventional convolutional learning does not change much
from the Gabor fit. Still, many of the kernels after Gabor learning
dissipates, and some form different structures. However, once the kernels are
close to zero in initialization, it does not change to other structures. It
seems that unnecessary filters dissipate by weight decay.

###  IV-B Gabor Fitting of first two layers

TABLE III: Accuracy of expanded ResNet-20 after second layer Gabor fitting | 7$\times$7 | 5$\times$5 | accuracy
---|---|---|---
5$\times$5 fitting | retrain | retrain | 91.92 $\pm$ 0.13
Gabor | retrain | 91.88 $\pm$ 0.07
retrain | Gabor | 91.70 $\pm$ 0.06
Gabor | Gabor | 91.60 $\pm$ 0.42
no fitting | Gabor | retrain | 91.96 $\pm$ 0.14
Figure 3: Schematic of second layer Gabor learning TABLE IV: Second layer prunings | kernel | level | pruned
---|---|---|---
First run
7 $\times$ 7 fitting | retrain \raisebox{-0.9pt}{1}⃝ |  7 $\times$ 7 | kernel | 10% $\pm$ 3%
channel | 5% $\pm$ 5%
GL \raisebox{-0.9pt}{2}⃝ |  7 $\times$ 7 | kernel | 63% $\pm$ 11%
channel | 35%$\pm$ 3%
Second run
5 $\times$ 5 fitting (7 $\times$ 7 results from GL) |  7 $\times$ 7 retrain 5 $\times$ 5 retrain \raisebox{-0.9pt}{3}⃝ | 7 $\times$ 7 | kernel | 50%$\pm$ 2%
channel | 49%$\pm$ 3%
5 $\times$ 5 | kernel | 61%$\pm$ 4%
channel | 51%$\pm$ 11%
7 $\times$ 7 GL 5 $\times$ 5 retrain \raisebox{-0.9pt}{4}⃝ | 7 $\times$ 7 | kernel | 78%$\pm$ 4%
channel | 59%$\pm$ 3%
5 $\times$ 5 | kernel | 59%$\pm$ 2%
channel | 41%$\pm$ 8%
7 $\times$ 7 retrain 5 $\times$ 5 GL \raisebox{-0.9pt}{5}⃝ | 7 $\times$ 7 | kernel | 70%$\pm$ 5%
channel | 66%$\pm$ 7%
5 $\times$ 5 | kernel | 94%$\pm$ 1%
channel | 71%$\pm$ 13%
7 $\times$ 7 GL 5 $\times$ 5 GL \raisebox{-0.95pt}{6}⃝ | 7 $\times$ 7 | kernel | 83%$\pm$ 3%
channel | 60%$\pm$ 6%
5 $\times$ 5 | kernel | 94%$\pm$ 1%
channel | 63%$\pm$ 6%
5 $\times$ 5 no fitting |  7 $\times$ 7 GL 5 $\times$ 5 retrain \raisebox{-0.9pt}{7}⃝ | 7 $\times$ 7 | kernel | 79%$\pm$ 5%
channel | 56%$\pm$ 6%
5 $\times$ 5 | kernel | 61%$\pm$ 7%
channel | 41%$\pm$ 8%

We have shown that learning the Gabor parameters in the first layer chooses
fewer features. We expand this idea to the second layer. The schematic of the
training and pruning is in Fig. 3. Each number in the figure indicates the
experiments in the following table.

The first five layers of ResNet-20 are all $3\times 3$ kernels with 16 output
channels, which we denote $3\times 3,16$. We expanded the kernels and removed
layers to make it to two layers of $7\times 7,16$, and $5\times 5,16$ without
shortcut connection. We trained the network four times, with the same
hyperparameters as the previous run, and used this as a reference. Note, the
original ResNet-20 has an accuracy of $92.38\pm 0.18$ after four runs, and it
also shows convergence. The Gabor fitting is the same as in previous cases,
except the scale $a$, which is scaled with each kernel’s maximum absolute
value.

We first prune unnecessary $7\times 7\times 3$ input channels from Gabor
learnings. To do so, we also pruned kernels in the second layer. Then, we
repeat the Gabor fitting for $5\times 5$ kernels.

Table III shows the accuracy of all the combinations of learnings with 5
$\times$ 5 Gabor fitting after the first Gabor learning. All of them have
similar results, i.e., less than 0.4% difference. Table IV shows the total
number of average pruning after the second run.Their results include the
pruning of $7\times 7\times 3$ channels in the first Gabor learnings, and it
shows that further prunings of the first layer are possible. In the second
layer, input features are no longer correlated to each other like in the case
of RGB inputs; hence the standard deviation of prunable channels is large.

Using second Gabor learning in either the first or the second layer,
\raisebox{-0.9pt}{4}⃝, \raisebox{-0.9pt}{5}⃝, and \raisebox{-0.9pt}{6}⃝, makes
significantly more prunable channels/kernels in the first layer. It might be
because the first Gabor learning makes the first layer close enough to Gabor
filters.

To check the influence on Gabor fitting in the second layer, we retrained the
network without 5 $\times$ 5 Gabor fitting and compared it with the network
with 5 $\times$ 5 Gabor fitting, i.e., comparing \raisebox{-0.9pt}{4}⃝ and
\raisebox{-0.9pt}{7}⃝. It shows that 5 $\times$ 5 Gabor fitting does not make
a clear difference in the number of prunable kernels without Gabor learning.
It implies that the structure of parameter space, not just the initial values,
is essential for learning.

Same as the first layer, 5 $\times$ 5 Gabor learning, \raisebox{-0.9pt}{5}⃝,
and \raisebox{-0.95pt}{6}⃝, can make significantly more prunable
channels/kernels. It could be possible to extend Gabor learning to further
layers.

##  V CONCLUSION AND DISCUSSION

In this paper, we propose to retrain networks by Gabor learning in the earlier
layers of CNNs to decrease the number of features/channels. Using strong
inductive biases by restricting features to Gabor filters help networks to
choose specific filters. As with classical image processing, we showed
choosing the right features for filters can make a smaller number of filters
after learning.

Our method does not compete with the other compressing methods since it does
not alter the structure of the networks. Hence, it can be used as a good
preprocessing for network compression.

## References

  * [1] A. Krizhevsky, I. Sutskever, and G. E. Hinton, “Imagenet classification with deep convolutional neural networks,” in _Advances in neural information processing systems_ , 2012, pp. 1097–1105.
  * [2] R. Girshick, J. Donahue, T. Darrell, and J. Malik, “Rich feature hierarchies for accurate object detection and semantic segmentation,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2014, pp. 580–587.
  * [3] J. Long, E. Shelhamer, and T. Darrell, “Fully convolutional networks for semantic segmentation,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2015, pp. 3431–3440.
  * [4] C. Szegedy, W. Liu, Y. Jia, P. Sermanet, S. Reed, D. Anguelov, D. Erhan, V. Vanhoucke, and A. Rabinovich, “Going deeper with convolutions,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2015, pp. 1–9.
  * [5] K. He, X. Zhang, S. Ren, and J. Sun, “Deep residual learning for image recognition,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2016, pp. 770–778.
  * [6] S. Zagoruyko and N. Komodakis, “Wide residual networks,” _arXiv preprint arXiv:1605.07146_ , 2016.
  * [7] Y. Cheng, D. Wang, P. Zhou, and T. Zhang, “A survey of model compression and acceleration for deep neural networks,” _arXiv preprint arXiv:1710.09282_ , 2017.
  * [8] J. Cheng, P.-s. Wang, G. Li, Q.-h. Hu, and H.-q. Lu, “Recent advances in efficient computation of deep convolutional neural networks,” _Frontiers of Information Technology & Electronic Engineering_, vol. 19, no. 1, pp. 64–77, 2018.
  * [9] Y. LeCun, J. S. Denker, and S. A. Solla, “Optimal brain damage,” in _Advances in neural information processing systems_ , 1990, pp. 598–605.
  * [10] B. Hassibi and D. G. Stork, “Second order derivatives for network pruning: Optimal brain surgeon,” in _Advances in neural information processing systems_ , 1993, pp. 164–171.
  * [11] S. Han, H. Mao, and W. J. Dally, “Deep compression: Compressing deep neural networks with pruning, trained quantization and huffman coding,” _arXiv preprint arXiv:1510.00149_ , 2015.
  * [12] S. Han, J. Pool, J. Tran, and W. Dally, “Learning both weights and connections for efficient neural network,” in _Advances in neural information processing systems_ , 2015, pp. 1135–1143.
  * [13] Y. Guo, A. Yao, and Y. Chen, “Dynamic network surgery for efficient dnns,” in _Advances in neural information processing systems_ , 2016, pp. 1379–1387.
  * [14] Y. He, X. Zhang, and J. Sun, “Channel pruning for accelerating very deep neural networks,” in _Proceedings of the IEEE International Conference on Computer Vision_ , 2017, pp. 1389–1397.
  * [15] Z. Liu, J. Li, Z. Shen, G. Huang, S. Yan, and C. Zhang, “Learning efficient convolutional networks through network slimming,” in _Proceedings of the IEEE International Conference on Computer Vision_ , 2017, pp. 2736–2744.
  * [16] J.-H. Luo, J. Wu, and W. Lin, “Thinet: A filter level pruning method for deep neural network compression,” in _Proceedings of the IEEE international conference on computer vision_ , 2017, pp. 5058–5066.
  * [17] S. Anwar, K. Hwang, and W. Sung, “Structured pruning of deep convolutional neural networks,” _ACM Journal on Emerging Technologies in Computing Systems (JETC)_ , vol. 13, no. 3, pp. 1–18, 2017.
  * [18] G. Huang, Y. Sun, Z. Liu, D. Sedra, and K. Q. Weinberger, “Deep networks with stochastic depth,” in _European conference on computer vision_. Springer, 2016, pp. 646–661.
  * [19] W. Wen, C. Wu, Y. Wang, Y. Chen, and H. Li, “Learning structured sparsity in deep neural networks,” in _Advances in neural information processing systems_ , 2016, pp. 2074–2082.
  * [20] J. Frankle and M. Carbin, “The lottery ticket hypothesis: Finding sparse, trainable neural networks,” _arXiv preprint arXiv:1803.03635_ , 2018.
  * [21] M. Denil, B. Shakibi, L. Dinh, M. Ranzato, and N. De Freitas, “Predicting parameters in deep learning,” in _Advances in neural information processing systems_ , 2013, pp. 2148–2156.
  * [22] X. Zhang, J. Zou, K. He, and J. Sun, “Accelerating very deep convolutional networks for classification and detection,” _IEEE transactions on pattern analysis and machine intelligence_ , vol. 38, no. 10, pp. 1943–1955, 2015\.
  * [23] S. d’Ascoli, H. Touvron, M. Leavitt, A. Morcos, G. Biroli, and L. Sagun, “Convit: Improving vision transformers with soft convolutional inductive biases,” _arXiv preprint arXiv:2103.10697_ , 2021.
  * [24] P. Ramachandran, N. Parmar, A. Vaswani, I. Bello, A. Levskaya, and J. Shlens, “Stand-alone self-attention in vision models,” _arXiv preprint arXiv:1906.05909_ , 2019.
  * [25] A. Srinivas, T.-Y. Lin, N. Parmar, J. Shlens, P. Abbeel, and A. Vaswani, “Bottleneck transformers for visual recognition,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2021, pp. 16 519–16 529.
  * [26] A. K. Jain, N. K. Ratha, and S. Lakshmanan, “Object detection using gabor filters,” _Pattern recognition_ , vol. 30, no. 2, pp. 295–309, 1997.
  * [27] J. G. Daugman, “Uncertainty relation for resolution in space, spatial frequency, and orientation optimized by two-dimensional visual cortical filters,” _JOSA A_ , vol. 2, no. 7, pp. 1160–1169, 1985.
  * [28] J. P. Jones and L. A. Palmer, “An evaluation of the two-dimensional gabor filter model of simple receptive fields in cat striate cortex,” _Journal of neurophysiology_ , vol. 58, no. 6, pp. 1233–1258, 1987.
  * [29] A. Coates, A. Ng, and H. Lee, “An analysis of single-layer networks in unsupervised feature learning,” in _Proceedings of the fourteenth international conference on artificial intelligence and statistics_ , 2011, pp. 215–223.
  * [30] H. Yao, L. Chuyi, H. Dan, and Y. Weiyu, “Gabor feature based convolutional neural network for object recognition in natural scene,” in _2016 3rd International Conference on Information Science and Control Engineering (ICISCE)_. IEEE, 2016, pp. 386–390.
  * [31] A. Calderón, S. Roa, and J. Victorino, “Handwritten digit recognition using convolutional neural networks and gabor filters,” _Proc. Int. Congr. Comput. Intell_ , 2003.
  * [32] B. Kwolek, “Face detection using convolutional neural networks and gabor filters,” in _International Conference on Artificial Neural Networks_. Springer, 2005, pp. 551–556.
  * [33] S.-Y. Chang and N. Morgan, “Robust cnn-based speech recognition with gabor filter kernels,” in _Fifteenth annual conference of the international speech communication association_ , 2014.
  * [34] S. S. Sarwar, P. Panda, and K. Roy, “Gabor filter assisted energy efficient fast learning convolutional neural networks,” in _2017 IEEE/ACM International Symposium on Low Power Electronics and Design (ISLPED)_. IEEE, 2017, pp. 1–6.
  * [35] S. Luan, C. Chen, B. Zhang, J. Han, and J. Liu, “Gabor convolutional networks,” _IEEE Transactions on Image Processing_ , vol. 27, no. 9, pp. 4357–4366, 2018.
  * [36] Y. Yuan, J. Zhang, and Q. Wang, “Deep gabor convolution network for person re-identification,” _Neurocomputing_ , vol. 378, pp. 387–398, 2020.
  * [37] H. Li, A. Kadav, I. Durdanovic, H. Samet, and H. P. Graf, “Pruning filters for efficient convnets,” _arXiv preprint arXiv:1608.08710_ , 2016.
  * [38] L. N. Smith, “Cyclical learning rates for training neural networks,” in _2017 IEEE Winter Conference on Applications of Computer Vision (WACV)_. IEEE, 2017, pp. 464–472.
  * [39] D. Gabor, “Theory of communication. part 1: The analysis of information,” _Journal of the Institution of Electrical Engineers-Part III: Radio and Communication Engineering_ , vol. 93, no. 26, pp. 429–441, 1946.
  * [40] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga _et al._ , “Pytorch: An imperative style, high-performance deep learning library,” in _Advances in neural information processing systems_ , 2019, pp. 8026–8037.
  * [41] A. Krizhevsky, G. Hinton _et al._ , “Learning multiple layers of features from tiny images,” 2009.
  * [42] K. Simonyan and A. Zisserman, “Very deep convolutional networks for large-scale image recognition,” _arXiv preprint arXiv:1409.1556_ , 2014.
  * [43] N. Zmora, G. Jacob, L. Zlotnik, B. Elharar, and G. Novik, “Neural network distiller,” Jun. 2018. [Online]. Available: https://doi.org/10.5281/zenodo.1297430
  * [44] P. Krähenbühl, C. Doersch, J. Donahue, and T. Darrell, “Data-dependent initializations of convolutional neural networks,” _arXiv preprint arXiv:1511.06856_ , 2015.


# Temporal Properties of the Compressible Magnetohydrodynamic Turbulence

Ka Ho Yuen  Theoretical Division, Los Alamos National Laboratory, Los Alamos,
NM 87545, USA  Corresponding Author<EMAIL_ADDRESS>,ORCID:
0000-0003-1683-9153 Hui Li  Theoretical Division, Los Alamos National
Laboratory, Los Alamos, NM 87545, USA  Huirong Yan  Deutsches Elektronen-
Synchrotron DESY, Germany  Institut für Physik & Astronomie, Universität
Potsdam, Germany

###### Abstract

The temporal property of the compressible magneto-hydrodynamic (MHD)
turbulence remains a fundamental unsolved question. Recent studies based on
the spatial-temporal analysis in the global frame of reference suggest that
the majority of fluctuation power in turbulence does not follow any of the MHD
wave dispersion relations but has very low temporal frequency with finite
wavenumbers. Here, we demonstrate that the Lorentzian broadening of the
dispersion relations of the three MHD modes where the nonlinear effects act
like the damping of a harmonic oscillator can explain many salient features of
frequency spectra for all MHD modes. The low frequency fluctuations are
dominated by modes with the low parallel wavenumbers that have been broadened
by the nonlinear processes. The Lorentzian broadening widths of the three MHD
modes exhibit scaling relations to the global frame wavenumbers and are
intrinsically related to energy cascade of each mode. Our results provide a
new window to investigate the temporal properties of turbulence which offers
insights for building a comprehensive understanding of the compressible MHD
turbulence.

## Main

There has been a long history in studying the temporal properties of MHD
turbulence [1, 2, 3], particularly the origin of low frequency temporal
fluctuations[4]. These low-frequency fluctuations have implications for
several problems in both space physics and astrophysics [5], including the
heating of the solar corona [6], the low-frequency “$1/f$ noise" in the solar
wind [7, 8], the formation and evolution of stars and molecular clouds in the
interstellar medium [9], as well as the propagation and acceleration of cosmic
rays [10, 11, 12]. Some of the earliest theoretical models came from the “2D
plus slab" model in the nearly incompressible magnetohydrodynamics [13, 14],
suggesting that a perpendicular cascade (i.e., the 2D fluctuations with global
frame parallel wavenumber $k_{\parallel}=0$111 The notion of $k_{\parallel}=0$
is defined in simulations that have periodic boundary condition parallel to
the global mean field direction. The fluctuation of $k_{\parallel}=0$ refers
to the infinite integration of variations along the global mean field.
However, in realistic astrophysical settings, all astrophysical objects are
bounded by a certain characteristic length $L$ [15] which defines the minimal
wavenumber $k_{min}=2\pi/L$, therefore the notion of $k_{\parallel}=0$ has
questionable physical meaning. However, the role of $k_{\parallel}=0$
fluctuations in the energy transfer of both in weak [16, 17] and strong [18]
turbulence should not be neglected. ) could generate nearly zero-frequency
fluctuations. Alternatively, the non-resonant three-wave interaction in strong
Alfvénic turbulence could generate MHD fluctuations at $k_{\parallel}\approx
0$ and $\omega\approx 0$ [18, 19]. Different views have fueled the debate
whether to treat these low frequency fluctuations as waves or nonlinear
structures[20]. Several physical pictures were put forth to interpret the low-
frequency fluctuations such as those from damped harmonic oscillators [21,
22], sweeping modes[20, 23], magnetosonic modes [24] and also transition from
weak to strong turbulence [25].

One interesting approach is to perform spatio-temporal analysis on simulated
turbulence fluctuations in order to extract their properties [7, 8, 26, 21,
25, 23, 24, 27, 28, 29] and gain unique insights on how turbulence evolves in
space and time. In particular, recent numerical studies [30, 31] have
quantified the spatio-temporal distribution of velocity, magnetic field and
density variations, showing that they not only deviate from the simple
dispersion relations for compressible MHD modes but also hold a dominant
fraction in low temporal frequencies, which is also observed from satellite
observations[32, 33, 34, 35].

In this paper, we explain quantitatively how the ubiquitous low frequency
$\omega$ fluctuations are physically generated in magnetized turbulence,
including all the compressible modes via the mode analysis [36] in the global
frame of reference[19]. We propose a new broadened Lorentzian profile that can
fit the simulation results. This profile allows us to quantify the
contributions by different mode groups, enabling a better understanding of the
frequency properties.We also discuss the implications of this new model for
the temporal behavior of MHD turbulence.

## A broadened Lorentzian model

In the linear wave theory [37, 36], small amplitude turbulence fluctuations
can be viewed as simple harmonic oscillators with natural frequencies
corresponding to the response frequencies of the three MHD waves (see Eq. 16).
The nonlinear terms (${\bf v}\cdot\nabla{\bf v}$, ${\bf\delta
B}\cdot\nabla{\bf\delta B}$) can be modeled as damping terms in the equation
of motion [8]. In this scenario, the turbulence system can be seen as a
collection of damped harmonic oscillators (resembling an oscillating Langevin
antenna in the case of incompressible Alfvénic turbulence [22]). By modeling
the nonlinear term in the equation of motion as $\omega_{nl}{\bf v}$, with the
exact functional form of $\omega_{nl}({\bf k})$ to be discussed in detail, the
spatial-temporal energy distribution function for a selected global frame
wavevector ${\bf k}=(k_{\parallel},k_{\perp})$ and a MHD mode is a broadened
Lorentzian distribution (see Supplementary material for deviation):

$E({\bf k},\omega)=\int dte^{i\omega t}E({\bf
k},t)\propto\frac{\omega_{nl}^{2}\omega_{\text{wave}}}{(\omega^{2}-\omega_{\text{wave}}^{2})^{2}+\omega_{nl}^{2}(\omega+v_{A}|{\bf
k}|\mu)^{2}}$ (1)

where $\omega_{\text{wave}}$ is the wave frequency of the individual mode (see
Eq. 16), $v_{A}$ is the Alfvén speed, and $\mu=\hat{\bf k}\cdot\hat{\bf B}$.
Consequently, the dispersion relations of the three MHD modes no longer follow
the linear form (Eq. 16), but are modulated by the nonlinear term (see
Supplementary Material).

Figure 1:  Temporal power $E(\omega)$ vs $\omega$ of velocity fluctuations
showing how MHD turbulence with different $\beta$ (upper row: $\beta\ll 1$,
simulation A1; middle row $\beta\ll 1$, A2; lower row, $\beta\gg 1$,
simulation A0; Blue: Alfvén, Red: Slow and Black: Fast mode, see Table 1)
produces significant fraction of low frequency fluctuations from Lorentzian
broadening (the scattered points in each panel, c.f. Eq.13) for three
different regimes of $k$. From the left: $k_{\parallel}=0$,
$k_{\parallel}<k_{\perp}$, $k_{\parallel}\geq k_{\perp}$. Each curve is
normalized by its own fluctuation power for a particular choice of wavenumber
$(k_{\parallel},k_{\perp})$ and fitted with Eq.1 (dash lines in each panel) .
The x-axis of each panel is normalized with respect to the Alfvénic frequency
($\tau_{A}^{-1}$) of the corresponding simulation.

Fig. 1 shows the $E(\omega)-\omega$ diagram, each curve normalized by its own
spectral power for a given wavevector (See Eq.17). Three MHD modes at selected
values of ${\bf k}=(k_{\parallel},k_{\perp})$ and different plasma $\beta$
(ratio of thermal to magnetic pressure, see Tab.2 for the definition of
symbols) are plotted, where the fluctuations are separated according to the
mode decomposition algorithm[36]. The mode decomposition algorithm assumes
negligible contributions from the nonlinear terms, which is not true
particularly when $k_{\parallel}$ is small, decreasing the accuracy of the
mode classifications (See Supplementary Material for discussions of the caveat
of the mode decomposition method in the strongly nonlinear systems). In Fig.1,
we intentionally include two cases for the same $\beta$ to illustrate that
higher turbulence levels (indicated by higher sonic Mach number $M_{s}$ and
Alfvénic Mach number $M_{A}$) produce more broadened Lorentzians (First two
row of Fig.1). All $E(\omega)$ curves are fitted by Eq.1 where we have made
$\omega_{nl}({\bf k})$ a fitting variable for different wave modes.
Furthermore, the peaks of the Lorentzian profiles correspond to the wave
eigenfrequency of different wave modes calculated from plasma $\beta$ and
wavevector ${\bf k}$ (Eq.10). The broadening behavior is generally consistent
with previous literature[8, 38, 22] and numerical simulations [39, 30, 31].
The nonlinear broadening by $\omega_{nl}({\bf k})$ for each peak has a strong
effect, extending the frequency distribution to both very low and high
frequency limits. Notice that the analysis is performed in the global frame,
how the local frame fluctuations are mapped into the global frame fluctuations
measured both numerically [7, 8, 28, 30] and observationally [34, 35] will be
addressed in the later section.

Figure 2:  $E(\omega)-\omega$ curves for the frequency power spectra of the
run A1 (See Tab. 1) where only Alfvénic modes are retained for different
choices of ${\bf k}$. The vertical dashed lines denote the wave eigen-
frequencies.

## Non-zero low frequency fluctuations are produced by nonlinear interactions

Figure 3: (a) The definition of low and high frequency fluctuations with
respect to the wave peak given a particular half-width. In this example the
half-width is $\omega_{nl}$. (b) Relative energy fraction of Alfvén modes
classified by frequencies in simulation A1. (c) Relative energy fraction for
integrated temporal power, where modes with global $k_{\parallel}=0$ are also
considered. (d) A comparison between the measured $\frac{E({\bf
k},\omega\rightarrow 0)}{E({\bf k},\omega=v_{A}k_{\parallel})}$ and the
theoretically predicted value. The red line denotes the equality in Eq.2.

Fig.1 also highlights another important aspect, namely the origin of the low
$\omega$ fluctuations. For instance, it can be seen that the $k_{\parallel}=0$
modes in all cases have an increasing power towards low $\omega$. To quantity
their contributions more clearly, Fig. 2 shows the frequency power spectra of
the run A1 (See Tab.1) where we extract only the Alfvén mode powers. For
frequency significantly less than $\tau_{A}^{-1}=v_{A}/L_{box}$ , there are
two types of contributions: those with $k_{\parallel}>0$ and those with
$k_{\parallel}=0$. Note that we have excluded the contributions by the ${\bf
k}$ modes within the injection region. For the $k_{\parallel}>0$ modes, their
contribution at low frequency, e.g., $\tau_{A}\omega=0.1$, is very small, and
they are mainly from the low frequency wing of the Lorentzian broadened
fluctuations. Modes with $k_{\parallel}=0$ (e.g., $k_{\perp}=3,4$) dominate
the low $\omega$ power. This implies that these finite frequency fluctuations
above $\omega=0$ are a result of the nonlinear interaction since their Alfvén
wave frequencies are zero when $k_{\parallel}=0$. The broadening of
$E(\omega)$ at $k_{\parallel}=0$ is always significant as long as
$\omega_{nl}$ is non-zero. The non-stationary nature (i.e. $E(\omega)\neq 0$
when $\omega\neq 0$) of the $k_{\parallel}=0$ mode is consistent with the
earlier literature that suggests the nonlinear terms (${v}\cdot\nabla{v}$,
${\delta B}\cdot\nabla{\delta B}$) require the $k_{\parallel}=0$ mode for
efficient energy transfer and the purely perpendicular 2D cascade [40, 18, 16,
8].

For finite $k_{\parallel}$ Alfvénic fluctuations, one can further quantify the
ratio of nonlinear component versus the linear (wave) component. Using the
incompressible Alfvénic fluctuations as an example, the analytical model (Eq.
1) fits the simulations very well. For the ease of discussion, we will define
the wave-like component as the integrated power around the linear wave
frequency $\omega_{\rm wave}({\bf k})$ within $\pm\Delta\omega({\bf k})$,
where $\Delta\omega({\bf k})$ is the half-width. We further denote the
fluctuations below $\omega_{\rm wave}-\Delta\omega$ and above
$\omega>\omega_{\rm wave}+\Delta\omega$ as low and high frequency
fluctuations, respectively (Panel (a) of Fig.3). For quantitative analysis, we
consider two choices of $\Delta\omega({\bf k})=0.5$ and $1$ $\omega_{nl}({\bf
k})$. Using simulation A1, we present the relative fraction of the wave-like,
low and high frequency fluctuations for different combinations of non-zero
$(k_{\parallel},k_{\perp})$ in Panel (b) of Fig.3. It can be seen that
significant fraction resides in the low and high frequency ranges, and the
wave-like fraction increases when using the larger $\Delta\omega$. Integrating
over all ${\bf k}$ modes (including $k_{\parallel}=0$ but excluding the
injection scale), we plot the relative contributions from the low, wave-like,
high, and $k_{\parallel}=0$ components, respectively, in Panel (c), again for
the two choices of $\Delta\omega$. Note that the relative fraction of the
$k_{\parallel}=0$ modes is a strong function of the minimum frequency used in
the analysis (which is chosen to be $0.1\omega_{A}$ in this plot). The
fraction contributed by the $k_{\parallel}=0$ modes is expected to increase if
a smaller minimum $\omega$, i.e. a longer time series, is employed in the
analysis.

The strength of the Lorentzian broadening is determined by the ratio of
$\omega_{nl}\omega_{\cal{A}}$ and $\omega^{2}_{wave}$ (see Eq. 1). When
$k_{\parallel}L_{inj}\gg 1$, the wave propagation frequency dominates over the
nonlinear feature, resulting in a sharper peak in the $E(\omega)$ distribution
(e.g. middle row of Fig. 1). Meanwhile, $E(\omega\rightarrow 0)$ is a non-zero
constant when $k_{\parallel}>0$, giving

$\frac{E({\bf k},\omega\rightarrow 0)}{E({\bf
k},\omega=v_{A}k_{\parallel})}=\frac{4\omega_{nl}^{2}/\omega_{\cal{A}}^{2}}{1+\omega_{nl}^{2}/\omega_{\cal{A}}^{2}}=\frac{4\chi^{2}}{1+\chi^{2}}\leavevmode\nobreak\
\leavevmode\nobreak\ ,$ (2)

where $\chi=\omega_{nl}({\bf k})/\omega_{\cal{A}}({\bf k})$. Notice that Eq.2
can be larger than 1 when $\omega_{nl}\gg\omega_{wave}$, indicating that the
low frequency fluctuations could have a higher amplitude even compared to that
at wave eigenfrequencies. To verify Eq.2, we compute the numerical value of
$E({\bf k},\omega\rightarrow 0)/E({\bf k},\omega=v_{A}k_{\parallel})$ in
simulation A1 and compare them to the predicted values using Eq.(2) in Panel
(d) of Fig.3, where we have used the $E(\omega)$ curves from simulation A1
with $2<|{\bf k}|<10$. We observe a reasonable agreement between the numerical
data to the theoretical prediction (Eq.2), though some data points are $\sim
2$ larger than the theoretical values. This deviation is mainly due to the
fact that the mode decomposition [36] is increasingly inaccurate as $\chi$
increases.

## Lorentzian broadening of the compressible modes and their low-frequency
contributions

The Lorentzian broadening exists in all three MHD modes but the broadening
strength and behavior are different (Fig .1). For the case of
$k_{\parallel}=0$ modes, the low frequency fluctuations are mostly contributed
by both Alfvén and slow modes, where the exact ratio is determined by the
relative energy fraction of the two modes (left column of Fig. 1). As
discussed previously, the broadening width at $k_{\parallel}=0$ only depends
on the value of $\omega_{nl}$. The similar slopes in $\omega$ for both Alfvén
and slow modes’ temporal power spectrum $E(\omega)$ suggest that their
$\omega_{nl}$ is similar in magnitude (See also Fig.4).

For modes with $k_{\parallel}>0$, slow mode contributes more low-frequency
fluctuations than the other two modes due to its lower $\omega_{wave}$ (middle
and right columns of Fig. 1). This behavior is more amplified in low $\beta$
where the slow wave speed is significantly smaller than the Alfvén speed for
the same ${\bf k}$. Different from the case of $k_{\parallel}=0$, the relative
fraction of low frequency fluctuations for $k_{\parallel}>0$ is governed by
$\chi$ parameter (Eq.2), which is the largest for slow modes. Therefore the
Lorentzian broadening is effectively stronger for slow modes, albeit the mode
fractions between slow and other two modes have to be taken into account[29].
In contrast, fast mode plays a negligible role in low-frequency fluctuations
because its wave speed is significantly higher than those of Alfvén and slow
modes, and is always non-zero unless $|{\bf k}|=0$.

Figure 4:  $\omega_{nl}$ (in units of $\tau_{A}^{-1}$) vs. $k_{\perp}$ (with
$k_{\parallel}=3$) or $|{\bf k}|$ for three different modes (Blue: Alfvén,
Red: Slow and Black: Fast mode). The red dash line in each panel corresponds
to $\omega_{nl}\propto k_{\perp}^{2/3}$ while the black line corresponds to
$\omega_{nl}\propto k^{1}$ for the left two panels, and $k^{3/4}$ for the
right panel.). The calculations are based on runs A0, A3 and A1, respectively.

## Scalings of nonlinear frequency $\omega_{nl}$ for different modes

To quantify the nonlinear broadening for each MHD mode and their dependence on
turbulence properties, we can use Eq. 1 to extract $\omega_{nl}$ from the
Lorentzian profiles. Fig. 4 shows the nonlinear frequencies of three subsonic,
sub-Alfvénic simulations with various $\beta$ as a function of $k_{\perp}$ or
$k=|{\bf k}|$. The effect of the local and global reference frame [41] is less
significant when applying the mode decomposition technique [36] since $M_{A}$
is small for our simulation (Tab.1). In the case of high $\beta$ (Left panel
of Fig. 4), both Alfvén and slow modes have their $\omega_{nl}$ scale in the
same way, $\propto k_{\perp}^{2/3}$. Fast modes cascade isotropically [36],
and therefore we plot the $\omega_{nl}$ for fast modes along $k$. We find a
scaling of $\omega_{nl}\propto k^{1}$ best fit our data for high and
intermediate $\beta$, suggesting a cascade of $E_{k}\propto k^{-2}$ for these
two cases. For low $\beta$ (Right panel of Fig.4), the data points are too
scattered to be conclusive, despite an apparent trend-line of $k^{3/4}$ is
observed for a limited range of $k$.

How do we understand the scaling trend in Fig.4? It is commonly assumed in the
case of strong Alfvénic turbulence that $\omega_{nl}\sim k_{\perp}\delta
v_{k}$ [18, 42]. However, the nonlinear time $\tau_{nl}\sim\omega_{nl}^{-1}$
does not necessarily correspond to the cascade time. The dependence of the
cascade time is actually a function of $\chi({\bf k})$, as noted in earlier
literature [43, 3, 42, 36, 44]. A rigorous closure calculation by Tripathi et
al. (in prep.) shows that the cascade time ($\omega_{tr}$) has the following
dependence on $\chi$:

$\omega_{tr}\approx\frac{\omega_{nl}^{2}}{\omega_{\text{wave}}(1+\chi)}$ (3)

In particular, for extreme cases of $\chi$:

$\displaystyle\omega_{tr}$
$\displaystyle\sim\frac{\omega_{nl}^{2}}{\omega_{\text{wave}}}$
$\displaystyle(\chi\ll 1)$ (4) $\displaystyle\omega_{tr}$
$\displaystyle\sim\omega_{nl}$ $\displaystyle(\chi\gg 1)$

where the first expression is the Iroshnikov–Kraichnan cascade rate[43, 3] and
the latter is commonly adopted for strong turbulence cascade [18]. The
constancy of the energy cascade rate $\delta v_{k}^{2}\omega_{tr}$ allows for
a quick quantification of the relation between $\omega_{nl}$ and spectral
power $E_{k}$ as functions of $k$. Writing:

$\text{const}\approx\begin{cases}\frac{\delta
v_{k}^{2}\omega_{nl}^{2}}{\omega_{wave}}\quad(\chi\ll 1)\\\ \delta
v_{k}^{2}\omega_{nl}\quad(\chi\gg 1)\end{cases}$ (5)

using $E_{k}\sim\delta v_{k}^{2}/k$ gives:

$E_{k}\propto\begin{cases}\omega_{nl}^{-2}\quad(\chi\ll 1)\\\
k^{-1}\omega_{nl}^{-1}\quad(\chi\gg 1)\end{cases}$ (6)

In the incompressible Alfvén and high-$\beta$ slow mode limit, where
$E_{k}\propto k_{\perp}^{-5/3}$, which implies that $\omega_{nl}\propto
k_{\perp}^{2/3}$, which are observed in Fig.4 for both Alfvén and slow modes
in all choices of $\beta$. For fast modes, the Iroshnikov–Kraichnan spectrum
[43, 3, 44] ($E_{k}\propto k^{-3/2}$) suggests $\omega_{nl}\propto k^{3/4}$,
which is only observed in the case of low $\beta$ case. For the measured
scaling $\omega_{nl}\propto k^{1}$ in the left and middle panels of Fig. 4,
Eq. 6 gives $E_{k}\propto k^{-2}$, which is commonly proposed as the
alternative scaling of fast modes [29].

## Discussion

The Lorentzian broadening effect is intrinsically caused by the nonlinear
effects in turbulence. However, the results are obtained in the global frame
with respect to the mean background magnetic field direction. In the case of
incompressible Alfvénic turbulence, many MHD turbulence studies [18, 42, 45]
have emphasized that the nonlinear fluctuations are generated perpendicular to
the local magnetic fields, characterized by the local wavevector
$(k_{\parallel,L},k_{\perp,L})$. These fluctuations, when viewed in a global
frame, could lead to a broadening effect in $\omega$. This can be understood
using the following picture: the spectral power of a particular global
wavevector $(k_{\parallel},k_{\perp})$ measured in the global frame could
collect many eddies of different sizes in the local frame [19]. These eddies
will contribute different spectral weights, with particularly higher weights
for the local eddies that have dimensions
$(k_{\parallel,L}^{-1},k_{\perp,L}^{-1})$ close to
$(k_{\parallel}^{-1},k_{\perp}^{-1})$. Each of these local wavevectors is
projected onto the global $k_{\parallel}$ axis, which gives a measured
$\omega$ as $\sim k_{\parallel,L}v_{A}$, along with its contribution to the
$E(\omega)$ spectrum. The variation of $k_{\parallel}$ due to the projections
of the local wavevectors is then translated to the $\omega$ space, resulting
in a broadening of the $E(\omega)$ spectrum.

The broadening width can be approximated as follows: The wavevector difference
between the local and global frames in the case of Alfvénic turbulence is
strongly correlated to the Alfvénic Mach number $M_{A}$ [41], i.e., $\delta
k_{\parallel}/k_{\parallel}\sim M_{A}\sim\delta v/v_{A}$. As a result, the
average dispersion of $\omega$ is $\delta\omega\sim\delta
k_{\parallel}v_{A}\sim k_{\parallel}\delta v$, centered at the global frame
$\omega\sim k_{\parallel}v_{A}$. In other words, the frame transformation
naturally generates a broadening in $\omega$ as long as $M_{A}$ is non-zero.
In comparison, the dispersion relation of fast modes has only weak dependence
on the wave vector direction. The projection effect for fast modes is
therefore marginal, accounting partially for their weaker nonlinear behavior.

As a remark, the concept of critical balance[18] is also closely connected to
the nature of non-zero low-frequency fluctuations. It has been argued that the
nonlinear timescale of incompressible Alfvén mode is approximately equal to
its wave propagation timescale, i.e., $\chi=1$, forming the basis of modern
MHD turbulence theory [18, 36], despite dissent views remain [25, 46, 47].
Refinement of the critical balance theory has been proposed in the community
in the incompressible limit [48]. Understanding the low-frequency fluctuations
provides valuable insights into the nonlinear timescales, as there is a
designated anisotropy scaling ($k_{\parallel}\propto k_{\perp}^{2/3}$ in the
case of Alfvén modes) related to the critical balance. However, whether such
balance exists for the compressible modes is still an unresolved question
[44].

## Data Availability

The data used in this work are listed in Tab.1. Data and the input file will
be available upon request.

Table 1: Table of simulations used in the current work  Model Name | $M_{s}$ | $M_{A}$ | $\beta$ | Energy Injection Rate
---|---|---|---|---
A0 | 0.20 | 0.50 | 12.5 | 0.01
A1 | 0.92 | 0.22 | 0.11 | 0.01
A2 | 0.17 | 0.04 | 0.13 | 0.0001
A3 | 0.45 | 0.46 | 2.06 | 0.001
The default parameters are:
$c_{s}=1$, $L_{box}=1$, $L_{inj}\geq 1/2$, $\langle\rho\rangle=1$,
Injection frequency $\tau_{s}\omega_{inj}=100$, resolution = $512^{3}$.
Table 2: Table of symbols used in this paper Symbol | Meaning
---|---
Shorthand Operators
$\langle x\rangle$ | Spatial average of x.
$\hat{\bf x}$ |  Unit vector of ${\bf x}$.
$\tilde{x}$ | Spatial-temporal Fourier transform of x.
Plasma Parameters
$\rho$ | Fluid density.
${\bf v}$ |  Fluid velocity vector, R.M.S. value of it =$v_{turb}$.
${\bf B}_{0}$ |  Mean Magnetic Field vector, its unit vector is sometimes denoted as $\hat{\lambda}$.
${\bf\delta B}$ |  $={\bf B}-{\bf B}_{0}$, the mean-subtracted magnetic field.
$v_{inj}$ | Injection velocity amplitude.
$v_{A}$ |  $=\frac{|{\bf B_{0}}|}{\sqrt{4\pi\langle\rho\rangle}}$, Alfvénic Speed. The vectorized form: ${\bf v}_{A}=\frac{\bf{\bf B}_{0}}{\sqrt{4\pi\langle\rho\rangle}}$
$c_{S}$ | Sonic Speed.
$M_{s}$ |  $=v_{inj}/c_{s}$, Sonic Mach number.
$M_{A}$ |  $=v_{inj}/v_{A}$, Alfvénic Mach number.
$\beta$ |  $=2M_{A}^{2}/M_{s}^{2}$, plasma compressibility.
$L_{box}$ |  Simulation domain size, always $=1$ in our set-up.
$L_{inj}$ |  Injection scale, $L_{box}/2$.
Fourier space & Mode parameters
$E(\omega)$ |  = $\tilde{E}(t)$, temporal energy spectrum.
$\omega$ | The temporal frequency.
${\bf k}$ | Wavevector (3D vector).
$\delta v_{k}$ |  Spatial Fourier amplitude of $v$ at ${\bf k}$.
$k_{\parallel},k_{\perp}$ | The parallel and perpendicular wavenumber in the global frame of reference.
$\mu=\cos\theta$ |  $=\hat{\bf k}\cdot\hat{\bf B_{0}}$, the cosine of the angle $\theta$ between wavevector and mean magnetic field.
${\cal A},{\cal S},{\cal F}$ | Alfvén, Slow, Fast Mode subscripts.
$\alpha$ |  $=\beta/2$.
$\cal{D}$ |  $=(1+\alpha)^{2}-4\alpha\cos^{2}\theta$.
$\zeta_{{\cal A},{\cal S},{\cal F}}$ |  Linearized Alfvénic, Slow and Fast wave vectors (c.f. Eq.10).
Timescales and nonlinearity parameters
$\tau_{s},\tau_{A}$ |  = $L_{box}/c_{s},L_{box}/v_{A}$, the sound crossing time and Alfvénic time, respectively.
$\omega_{nl}$ | Nonlinear frequency.
$\omega_{tr}$ |  Cascade frequency, c.f. Eq.3.
$\omega_{wave}$ |  Wave frequency, c.f .Eq.16.
$\omega_{inj}$ | Injection Frequency.
$\chi$ |  Natural nonlinear parameter, defined as $\omega_{nl}/\omega_{wave}$.

## Code Availability

### Numerical simulations

The numerical simulations are performed with Athena++ [49].We summarize the
simulations in Tab. 1. Our data are time series of three-dimensional, triply
periodic, isothermal MHD simulations with continuous force driving via direct
spectral injection 222The reason why we did not employ the Ornstein–Uhlenbeck
forcing is because we would like to have more control on the injected values
of ${\bf k}$ and the injection frequency $\omega_{inj}$. unless specified. We
run our simulations for at least 10 sound crossing time
($\tau_{s}=L_{box}/c_{s}$) and take snapshots at $\Delta\tau=\tau_{s}/100$ to
ensure the time-axis sampling satisfies the condition that:

$\Delta\tau_{required}<\frac{L_{box}}{v_{fastest}}$ (7)

where $L_{box}$ is the size of the simulation domain, and $v_{fastest}$ is the
fastest speed in the numerical simulations. The typical parameters of our
simulations are listed in Tab.1. The injection is performed so that we only
have eddies at scales $L_{inj}/L_{box}\geq 1/2$, which corresponds to $|{\bf
k}|\leq 2$. All simulations are driven solendoially. All numerical simulations
are truncated in Fourier space to $128^{3}$ regardless of its original size to
save computational resources. That will not change the statistics of spatio-
temporal spectrum for $|{\bf k}|<128$.

### Analysis

Analysis are performed by Julia with the packages available upon request.

## References

  * [1] Kolmogorov, A.  The Local Structure of Turbulence in Incompressible Viscous Fluid for Very Large Reynolds’ Numbers.  _Akademiia Nauk SSSR Doklady_ 30, 301–305 (1941).
  * [2] Iroshnikov, P. S.  Turbulence of a Conducting Fluid in a Strong Magnetic Field.  _Soviet Ast._ 7, 566 (1964).
  * [3] Kraichnan, R. H.  Inertial-Range Spectrum of Hydromagnetic Turbulence.  _Physics of Fluids_ 8, 1385–1387, DOI: 10.1063/1.1761412 (1965).
  * [4] Matthaeus, W. H.  Turbulence in space plasmas: Who needs it?  _Physics of Plasmas_ 28, 032306, DOI: 10.1063/5.0041540 (2021).
  * [5] Zank, G. P. _et al._ The Transport of Low-frequency Turbulence in Astrophysical Flows. I. Governing Equations.  _ApJ_ 745, 35, DOI: 10.1088/0004-637X/745/1/35 (2012).
  * [6] Matthaeus, W. H., Zank, G. P., Oughton, S., Mullan, D. J. & Dmitruk, P.  Coronal Heating by Magnetohydrodynamic Turbulence Driven by Reflected Low-Frequency Waves.  _ApJ_ 523, L93–L96, DOI: 10.1086/312259 (1999).
  * [7] Dmitruk, P. & Matthaeus, W. H.  Low-frequency 1/f fluctuations in hydrodynamic and magnetohydrodynamic turbulence.  _Phys. Rev. E_ 76, 036305, DOI: 10.1103/PhysRevE.76.036305 (2007).
  * [8] Dmitruk, P. & Matthaeus, W. H.  Waves and turbulence in magnetohydrodynamic direct numerical simulations.  _Physics of Plasmas_ 16, 062304, DOI: 10.1063/1.3148335 (2009).
  * [9] McKee, C. F. & Ostriker, E. C.  Theory of Star Formation.  _ARA &A_ 45, 565–687, DOI: 10.1146/annurev.astro.45.051806.110602 (2007).  0707.3514.
  * [10] Zank, G. P., Matthaeus, W. H., Bieber, J. W. & Moraal, H.  The radial and latitudinal dependence of the cosmic ray diffusion tensor in the heliosphere.  _J. Geophys. Res._ 103, 2085–2098, DOI: 10.1029/97JA03013 (1998).
  * [11] Yan, H. & Lazarian, A.  Scattering of Cosmic Rays by Magnetohydrodynamic Interstellar Turbulence.  _Phys. Rev. Lett._ 89, 281102, DOI: 10.1103/PhysRevLett.89.281102 (2002).  astro-ph/0205285.
  * [12] Shalchi, A., Bieber, J. W., Matthaeus, W. H. & Schlickeiser, R.  Parallel and Perpendicular Transport of Heliospheric Cosmic Rays in an Improved Dynamical Turbulence Model.  _ApJ_ 642, 230–243, DOI: 10.1086/500728 (2006).
  * [13] Matthaeus, W. H., Goldstein, M. L. & Roberts, D. A.  Evidence for the presence of quasi-two-dimensional nearly incompressible fluctuations in the solar wind.  _J. Geophys. Res._ 95, 20673–20683, DOI: 10.1029/JA095iA12p20673 (1990).
  * [14] Zank, G. P. & Matthaeus, W. H.  Waves and turbulence in the solar wind.  _J. Geophys. Res._ 97, 17189–17194, DOI: 10.1029/92JA01734 (1992).
  * [15] Chepurnov, A. & Lazarian, A.  Extending the Big Power Law in the Sky with Turbulence Spectra from Wisconsin H$\alpha$ Mapper Data.  _ApJ_ 710, 853–858, DOI: 10.1088/0004-637X/710/1/853 (2010).  0905.4413.
  * [16] Ng, C. S. & Bhattacharjee, A.  Interaction of Shear-Alfven Wave Packets: Implication for Weak Magnetohydrodynamic Turbulence in Astrophysical Plasmas.  _ApJ_ 465, 845, DOI: 10.1086/177468 (1996).
  * [17] Nazarenko, S.  2D enslaving of MHD turbulence.  _New Journal of Physics_ 9, 307, DOI: 10.1088/1367-2630/9/8/307 (2007).
  * [18] Goldreich, P. & Sridhar, S.  Toward a Theory of Interstellar Turbulence. II. Strong Alfvenic Turbulence.  _ApJ_ 438, 763, DOI: 10.1086/175121 (1995).
  * [19] Cho, J. & Vishniac, E. T.  The Anisotropy of Magnetohydrodynamic Alfvénic Turbulence.  _ApJ_ 539, 273–282, DOI: 10.1086/309213 (2000).  astro-ph/0003403.
  * [20] Zhou, Y., Matthaeus, W. H. & Dmitruk, P.  Colloquium: Magnetohydrodynamic turbulence and time scales in astrophysical and space plasmas.  _Reviews of Modern Physics_ 76, 1015–1035, DOI: 10.1103/RevModPhys.76.1015 (2004).
  * [21] Howes, G. G. & Nielson, K. D.  Alfvén wave collisions, the fundamental building block of plasma turbulence. I. Asymptotic solution.  _Physics of Plasmas_ 20, 072302, DOI: 10.1063/1.4812805 (2013).  1306.1455.
  * [22] TenBarge, J. M., Howes, G. G., Dorland, W. & Hammett, G. W.  An oscillating Langevin antenna for driving plasma turbulence simulations.  _Computer Physics Communications_ 185, 578–589, DOI: 10.1016/j.cpc.2013.10.022 (2014).  1305.2212.
  * [23] Lugones, R., Dmitruk, P., Mininni, P. D., Wan, M. & Matthaeus, W. H.  On the spatio-temporal behavior of magnetohydrodynamic turbulence in a magnetized plasma.  _Physics of Plasmas_ 23, 112304, DOI: 10.1063/1.4968236 (2016).  1608.08244.
  * [24] Andrés, N. _et al._ Interplay between Alfvén and magnetosonic waves in compressible magnetohydrodynamics turbulence.  _Physics of Plasmas_ 24, 102314, DOI: 10.1063/1.4997990 (2017).  1707.09190.
  * [25] Meyrand, R., Galtier, S. & Kiyani, K. H.  Direct Evidence of the Transition from Weak to Strong Magnetohydrodynamic Turbulence.  _Phys. Rev. Lett._ 116, 105002, DOI: 10.1103/PhysRevLett.116.105002 (2016).  1509.06601.
  * [26] Svidzinski, V. A., Li, H., Rose, H. A., Albright, B. J. & Bowers, K. J.  Particle in cell simulations of fast magnetosonic wave turbulence in the ion cyclotron frequency range.  _Physics of Plasmas_ 16, 122310, DOI: 10.1063/1.3274559 (2009).
  * [27] Lugones, R., Dmitruk, P., Mininni, P. D., Pouquet, A. & Matthaeus, W. H.  Spatio-temporal behavior of magnetohydrodynamic fluctuations with cross-helicity and background magnetic field.  _Physics of Plasmas_ 26, 122301, DOI: 10.1063/1.5129655 (2019).  1911.03400.
  * [28] Yang, L. P. _et al._ Energy occupation of waves and structures in 3D compressive MHD turbulence.  _MNRAS_ 488, 859–867, DOI: 10.1093/mnras/stz1747 (2019).
  * [29] Makwana, K. D. & Yan, H.  Properties of Magnetohydrodynamic Modes in Compressively Driven Plasma Turbulence.  _Physical Review X_ 10, 031021, DOI: 10.1103/PhysRevX.10.031021 (2020).  1907.01853.
  * [30] Gan, Z., Li, H., Fu, X. & Du, S.  On the Existence of Fast Modes in Compressible Magnetohydrodynamic Turbulence.  _ApJ_ 926, 222, DOI: 10.3847/1538-4357/ac4d9d (2022).  2201.07965.
  * [31] Fu, X., Li, H., Gan, Z., Du, S. & Steinberg, J.  Nature and Scalings of Density Fluctuations of Compressible Magnetohydrodynamic Turbulence with Applications to the Solar Wind.  _ApJ_ 936, 127, DOI: 10.3847/1538-4357/ac8802 (2022).  2207.09490.
  * [32] Narita, Y., Sahraoui, F., Goldstein, M. L. & Glassmeier, K. H.  Magnetic energy distribution in the four-dimensional frequency and wave vector domain in the solar wind.  _Journal of Geophysical Research (Space Physics)_ 115, A04101, DOI: 10.1029/2009JA014742 (2010).
  * [33] Narita, Y., Gary, S. P., Saito, S., Glassmeier, K. H. & Motschmann, U.  Dispersion relation analysis of solar wind turbulence.  _Geophys. Res. Lett._ 38, L05101, DOI: 10.1029/2010GL046588 (2011).
  * [34] Zhao, L. L., Zank, G. P., Nakanotani, M. & Adhikari, L.  Observations of Waves and Structures by Frequency-Wavenumber Spectrum in Solar Wind Turbulence.  _ApJ_ 944, 98, DOI: 10.3847/1538-4357/acb33b (2023).
  * [35] Zhao, S., Yan, H., Liu, T. Z., Yuen, K. H. & Wang, H.  Satellite Observations of the Alfvénic Transition from Weak to Strong Magnetohydrodynamic Turbulence.  _arXiv e-prints_ arXiv:2301.06709, DOI: 10.48550/arXiv.2301.06709 (2023).  2301.06709.
  * [36] Cho, J. & Lazarian, A.  Compressible magnetohydrodynamic turbulence: mode coupling, scaling relations, anisotropy, viscosity-damped regime and astrophysical implications.  _MNRAS_ 345, 325–339, DOI: 10.1046/j.1365-8711.2003.06941.x (2003).  astro-ph/0301062.
  * [37] Cho, J. & Lazarian, A.  Compressible Sub-Alfvénic MHD Turbulence in Low- $\beta$ Plasmas.  _Phys. Rev. Lett._ 88, 245001, DOI: 10.1103/PhysRevLett.88.245001 (2002).  astro-ph/0205282.
  * [38] TenBarge, J. M. & Howes, G. G.  Evidence of critical balance in kinetic Alfvén wave turbulence simulations a).  _Physics of Plasmas_ 19, 055901, DOI: 10.1063/1.3693974 (2012).  1201.0056.
  * [39] Brodiano, M., Andrés, N. & Dmitruk, P.  Spatiotemporal Analysis of Waves in Compressively Driven Magnetohydrodynamics Turbulence.  _ApJ_ 922, 240, DOI: 10.3847/1538-4357/ac2834 (2021).  2107.09212.
  * [40] Montgomery, D. & Matthaeus, W. H.  Anisotropic Modal Energy Transfer in Interstellar Turbulence.  _ApJ_ 447, 706, DOI: 10.1086/175910 (1995).
  * [41] Yuen, K. H., Yan, H. & Lazarian, A.  Anomalous compressible mode generation by global frame projections of pure Alfven mode.  _MNRAS_ 521, 530–545, DOI: 10.1093/mnras/stad287 (2023).  2301.13344.
  * [42] Lazarian, A. & Vishniac, E. T.  Reconnection in a Weakly Stochastic Field.  _ApJ_ 517, 700–718, DOI: 10.1086/307233 (1999).  astro-ph/9811037.
  * [43] Iroshnikov, P. S.  Turbulence of a Conducting Fluid in a Strong Magnetic Field.  _AZh_ 40, 742 (1963).
  * [44] Galtier, S.  Fast magneto-acoustic wave turbulence and the Iroshnikov-Kraichnan spectrum.  _Journal of Plasma Physics_ 89, 905890205, DOI: 10.1017/S0022377823000259 (2023).  2303.00643.
  * [45] Chandran, B. D. G.  Scattering of Energetic Particles by Anisotropic Magnetohydrodynamic Turbulence with a Goldreich-Sridhar Power Spectrum.  _Phys. Rev. Lett._ 85, 4656–4659, DOI: 10.1103/PhysRevLett.85.4656 (2000).  astro-ph/0008498.
  * [46] Oughton, S. & Matthaeus, W. H.  Critical Balance and the Physics of Magnetohydrodynamic Turbulence.  _ApJ_ 897, 37, DOI: 10.3847/1538-4357/ab8f2a (2020).  2006.04677.
  * [47] Chhiber, R., Matthaeus, W. H., Oughton, S. & Parashar, T. N.  A detailed examination of anisotropy and timescales in three-dimensional incompressible magnetohydrodynamic turbulence.  _Physics of Plasmas_ 27, 062308, DOI: 10.1063/5.0005109 (2020).  2005.08815.
  * [48] Mallet, A., Schekochihin, A. A. & Chandran, B. D. G.  Refined critical balance in strong Alfvenic turbulence.  _MNRAS_ 449, L77–L81, DOI: 10.1093/mnrasl/slv021 (2015).  1406.5658.
  * [49] Stone, J. M., Tomida, K., White, C. J. & Felker, K. G.  The Athena++ Adaptive Mesh Refinement Framework: Design and Magnetohydrodynamic Solvers.  _ApJS_ 249, 4, DOI: 10.3847/1538-4365/ab929b (2020).  2005.06651.

## Competing interests

The authors declare no competing interests.

## Supplementary information

## Acknowledgements

Research presented in this article was supported by the Laboratory Directed
Research and Development program of Los Alamos National Laboratory under
project numbers 20220700PRD1 and 20220107DR, and a U.S. Department of Energy
Fusion Energy Science project. This research used resources provided by the
Los Alamos National Laboratory Institutional Computing Program, which is
supported by the U.S. Department of Energy National Nuclear Security
Administration under Contract No. 89233218CNA000001. This research also used
resources of the National Energy Research Scientific Computing Center (NERSC),
a U.S. Department of Energy Office of Science User Facility located at
Lawrence Berkeley National Laboratory, operated under Contract No. DE-
AC02-05CH11231 using NERSC award HEP-ERCAP-m2407,m3122 (PI: Fan Guo, LANL) and
FES-ERCAP-m4239 (PI: KHY, LANL). Inspiring discussions with Jungyeon Cho
(Chungnam National University), Patrick Diamond (UCSD), Fan Guo (LANL), Alex
Lazarian (UW Madison), Willian Matthaeus (U. Delaware), Bindesh Tripathi (UW
Madison), Ka Wai Ho (UW Madison/LANL) and Ethan Vishniac (JHU) are
acknowledged.

## Author contributions

Yuen and Li initiated and designed the project. Yuen performed the simulations
and analyzed the main result of the manuscript. Yuen and Yan derived the
nonlinear timescales. All authors contributed in writing, editing and
approving the manuscript.

## Methodology

The simulations shown in Tab.1 are in the form of time series, say ${\bf
v}({\bf r},t)$. We perform a full 4D Fourier transform:

${\bf v}({\bf r},t)=\sum_{{\bf k},\omega}\tilde{\bf v}({\bf
k},\omega)e^{i({\bf k}\cdot{\bf r}-\omega t)}$ (8)

Taking the cylindrical coordinate ${\bf
k}\rightarrow(k_{\parallel},k_{\perp})$, the spatio-temporal power spectrum
$E(k_{\parallel},k_{\perp},\omega)$ is given by:

$E(k_{\parallel},k_{\perp},\omega)=\oint k_{\perp}d\Omega|\tilde{\bf v}({\bf
k},\omega)|^{2}=2\pi k_{\perp}|\tilde{\bf v}({\bf
k}=(k_{\parallel},k_{\perp}),\omega)|^{2}$ (9)

which cautions have to be taken that the coordinate transform adds an
additional factor of $k_{\perp}$.

### Mode Decomposition when $\chi\ll 1$

The MHD modes are obtained by performing the dot product $\tilde{\bf
v}_{X}=(\tilde{\bf v}\cdot\hat{\zeta}_{X})\hat{\zeta}_{X}$ to the turbulence
variables, where
$X=\cal{A}\text{(Alfv\'{e}n)},\cal{S}\text{(Slow)},\cal{F}\text{(Fast)}$ and
$\zeta_{X}$ are[36]:

$\displaystyle\zeta_{\cal A}$ $\displaystyle\propto\hat{\bf
k}\times\hat{\bf\lambda}$ (10) $\displaystyle\zeta_{\cal S}$
$\displaystyle\propto(-1+\alpha-\sqrt{\cal{D}})({\bf
k}\cdot\hat{\bf\lambda})\hat{\bf\lambda}+(1+\alpha-\sqrt{\cal{D}})(\hat{\bf\lambda}\times({\bf
k}\times\hat{\bf\lambda}))$ $\displaystyle\zeta_{\cal F}$
$\displaystyle\propto(-1+\alpha+\sqrt{\cal{D}})({\bf
k}\cdot\hat{\bf\lambda})\hat{\bf\lambda}+(1+\alpha+\sqrt{\cal{D}})(\hat{\bf\lambda}\times({\bf
k}\times\hat{\bf\lambda}))$

where $\lambda=\langle\hat{\bf B}\rangle$, $\alpha=\beta/2$,
$\beta=2M_{A}^{2}/M_{s}^{2}$,
${\cal{D}}=(1+\alpha)^{2}-4\alpha\cos^{2}\theta$, $\cos\theta=\mu=\hat{\bf
k}\cdot\hat{\bf\lambda}$. We have to emphasize that the decomposition was
proposed assuming negligible non-linear terms [36], as they are treated as
second order terms and should be zeroed out after linear perturbations.
Furthermore, the decomposition is performed in the global frame of
reference[19]. For the case where $\chi$ is not small, we expect the
decomposition to be less accurate. Example signatures include the lower right
panel of Fig.1, where fast wave eigen-peak appears in slow mode $E(\omega)$.
However, a simple eigenvalue analysis suggests that the appearance of the non-
linear term does not change the mean direction of the eigenvectors to the
second order.

The 4D energy spectrum of a particular mode is given by (c.f. Eq.9):

$E_{X=({\cal A},{\cal S},{\cal F})}(k_{\parallel},k_{\perp},\omega)=\oint
k_{\perp}d\Omega|\tilde{\bf v}({\bf k},\omega)|^{2}=2\pi k_{\perp}|\tilde{\bf
v}({\bf k}=(k_{\parallel},k_{\perp}),\omega)\cdot\hat{\zeta}_{X}|^{2}$ (11)

###  The quantitative model of MHD turbulence

#### Deviation of the Lorentzian function (Eq.1) from MHD equations

To derive the Lorentzian function used in Eq.1, we need to make the following
approximations to the ideal MHD equations with isothermal equation of state.

  1. 1.

We assume first order perturbations on $\rho,v$ & $B$, ${\bf B}={\bf
B}_{0}+\delta{\bf B},z=ln(\rho/\langle\rho\rangle),$ and write the Alfvénic
vector ${\bf v}_{A}=\frac{{\bf B}_{0}}{\sqrt{4\pi\langle\rho\rangle}}$.

  2. 2.

We preserve all nonlinear terms ($({\bf v}\cdot\nabla){\bf v},({\bf
B}\cdot\nabla){\bf B}$) and approximate them as $\omega_{nl}{\bf v}/2$ and
$\omega_{nl}{\bf B}/2$, respectively. The approximation is exact for pure
Alfvénic turbulence [22].

  3. 3.

We perform Fourier transform, writing $\partial_{t}\rightarrow-i\omega$ and
$\nabla\rightarrow i{\bf k}$.

Then the equation of motion in the Fourier space is given by
($c_{n}^{2}=c_{s}^{2}+\frac{v_{A}^{2}}{2}$):

$(-i\omega+\frac{\omega_{nl}}{2})^{2}\tilde{\bf v}+{\bf k}c_{n}^{2}{\bf
k}\cdot{\bf\tilde{v}}-(i{\bf v}_{A}\cdot{\bf
k}+\frac{\omega_{nl}}{2})^{2}\tilde{\bf v}=0$ (12)

The Lagrangian power spectrum is defined by the squared amplitude of the
Fourier transformed velocity $\tilde{\bf v}$. Solving Eq.12 against the plane
wave solution (the Fourier transform of $\delta$-function), the Lagrangian
power spectrum is given by (using
$\omega_{{\cal{A},\cal{S},\cal{F}}}^{2}=(v_{A}|{\bf k}|\mu)^{2}+c_{n}^{2}({\bf
k}\cdot{\bf\tilde{\hat{v}}})^{2}$, c.f. Eq.16):

$E(\omega,{\bf k})=|\tilde{\bf
v}|^{2}=\frac{1}{(\omega^{2}-\omega_{\cal{A},\cal{S},\cal{F}}^{2})^{2}+\omega_{nl}^{2}(\omega+v_{A}k\mu)^{2}}$
(13)

Special cases for Incompressible regime: In the limit of incompressible
turbulence $\omega_{\cal{A}}=\omega_{\cal{S}}$ and Eq.13 reduces to the
following form:

$E(\omega,{\bf k})\propto\frac{1}{(\omega+\omega_{\cal
A})^{2}}\frac{1}{(\omega-\omega_{\cal A})^{2}+\omega_{nl}^{2}}$ (14)

Special cases for $\mu=0$: In the limit of $\mu=0$, Eq.13 reduces to the
following form:

$E(\omega,{\bf
k})\propto\frac{1}{\omega^{2}}\frac{1}{(\omega-\omega_{\cal{A},\cal{S},\cal{F}}(\mu=0))^{2}+\omega_{nl}^{2}}$
(15)

###  Linear Wave Frequencies and Dispersion Relations of MHD Modes

The dispersion relations of the three MHD waves are given by[36]
($\cal{A}=\text{Alfv\'{e}n},\cal{S}=\text{Slow},\cal{F}=\text{Fast}$):

$\displaystyle\omega_{\cal{A}}$ $\displaystyle=v_{A}\mu|{\bf k}|$ (16)
$\displaystyle\omega_{\cal{S}}$
$\displaystyle=\Big{(}\frac{1}{2}\left[v_{A}^{2}+c_{s}^{2}-\sqrt{(v_{A}^{2}+c_{s}^{2})^{2}-4v_{A}^{2}c_{s}^{2}\mu^{2}}\right]\Big{)}^{\frac{1}{2}}|{\bf
k}|$ $\displaystyle\omega_{\cal{F}}$
$\displaystyle=\Big{(}\frac{1}{2}\left[v_{A}^{2}+c_{s}^{2}+\sqrt{(v_{A}^{2}+c_{s}^{2})^{2}-4v_{A}^{2}c_{s}^{2}\mu^{2}}\right]\Big{)}^{\frac{1}{2}}|{\bf
k}|$

where $v_{A}$ is the Alfvén speed, $c_{s}$ is the sound speed and
$\mu=\cos\theta=\hat{\bf k}\cdot\langle\hat{\bf B}\rangle$.

### Lorentzian Profile

We normalize the energy spectrum by each of its integrated value:

$\bar{E}_{X}(\omega;k_{\parallel},k_{\perp})=\frac{E_{X}(\omega;k_{\parallel},k_{\perp})}{\int
d\omega E_{X}(\omega;k_{\parallel},k_{\perp})}$ (17)

The Lorentzian profiles (c.f. Figures 1,2,4) are fitted by least square method
to the functional form

$y(\omega;[A,B,C,D])={\frac{D}{(\omega^{2}-A^{2})^{2}+(\omega+B)^{2}C^{2}}}$
(18)

normalized to unity. Here $A\rightarrow\omega_{wave}$, $B\rightarrow{\bf
v}_{A}\cdot{\bf k}=v_{A}|{\bf k}|\mu$, $C=\omega_{nl}$ (c.f. Eq.1). The
fitting was done by a bounded 4-parameter regression on $\log(\omega)$ and
$\log(1/y)$ with $A,C,D\geq 0$ and $B\in\mathcal{R}$.

### Scaling Relation of $\omega_{nl}$ for incompressible Alfvénic turbulence

In the strong B-field limit, the nonlinear frequency for the Alfvén mode can
be approximated by the Kolmogorov scaling perpendicular to the mean field
(i.e., $v/v_{\text{inj}}\propto(l/L_{\text{inj}})^{1/3}$) [42]:

$\omega_{nl}=v_{\text{turb}}l_{\perp}^{-1}\sim
v_{\text{inj}}L_{\text{inj}}^{-1/3}k_{\perp}^{2/3}$ (19)

where $v_{\text{inj}}$ and $L_{\text{inj}}$ are the velocity and length scale
at injection. Notice that the nonlinear time for the three different modes is
actually slightly different (c.f. Fig. 4), contrary to the common belief [44],
as the nonlinear interactions act differently on the three modes. However, the
conclusion that the fast modes are in weak turbulence [44] when Alfvén waves
are in strong turbulence [18] in small scales is phenomenologically correct .