* [33] L. Silberman and A. Venkatesh, On quantum unique ergodicity for locally symmetric spaces, Geom. Funct. Anal. 17 (2007), 960-998.
  * [34] J. Sjöstrand and M. Zworski, Asymptotic distribution of resonances for convex obstacles, Acta Mathematica 183, no. 2, 191-253 (1999).
  * [35] K. Soundararajan, Quantum unique ergodicity for $SL_{2}(\mathbb{Z})/\mathbb{H}$, Ann. Math. 172, 1529-1538. (2010)
  * [36] S. Zelditch, Uniform distribution of eigenfunctions on compact hyperbolic surfaces, Duke Math. J. 55, 919-941 (1987).
  * [37] S. Zelditch, Quantum ergodicity on the sphere, Comm. Math. Phys. 164, 61-71 (1992).
  * [38] S. Zelditch and M. Zworski, Ergodicity for eigenfunctions of ergodic billiards, Comm. Math. Phys. 175, 673-682 (1996).
  * [39] M. Zworski, Semiclassical Analysis, Graduate Studies in Mathematics 138, Amer. Math. Soc. (2012).


# Converting Anyone’s Voice:
End-to-End Expressive Voice Conversion with a Conditional Diffusion Model

###### Abstract

Expressive voice conversion (VC) conducts speaker identity conversion for
emotional speakers by jointly converting speaker identity and emotional style.
Emotional style modeling for arbitrary speakers in expressive VC has not been
extensively explored. Previous approaches have relied on vocoders for speech
reconstruction, which makes speech quality heavily dependent on the
performance of vocoders. A major challenge of expressive VC lies in emotion
prosody modeling. To address these challenges, this paper proposes a fully
end-to-end expressive VC framework based on a conditional denoising diffusion
probabilistic model (DDPM). We utilize speech units derived from self-
supervised speech models as content conditioning, along with deep features
extracted from speech emotion recognition and speaker verification systems to
model emotional style and speaker identity. Objective and subjective
evaluations show the effectiveness of our framework. Codes and samples are
publicly available.††Codes & Speech Samples: https://a2023aa.github.io/DEVC/
Junchen Lu and Kun Zhou contributed to this work during their internship at UT
Dallas.

##  1 Introduction

Emotions play a vital role in natural speech, as they convey one’s feelings,
moods, and personality [1]. Speech expressiveness encompasses a diverse range
of emotions [2]. Expressive voice conversion aims to jointly perform speaker
identity and style transfer for emotional speakers, which poses huge potential
in movie dubbing, voice acting, and human-computer interaction [3, 4, 5].

Voice conversion (VC) is the task of converting speaker identity while
preserving linguistic content [6]. Early VC studies focus on learning a
statistical mapping between speech recordings of source and target speakers
[7, 8]. Recent advancements in deep learning-based approaches significantly
improve the performance of VC by learning the disentanglement across various
speech characteristics. For instance, separate encoders guided by appropriate
constraints have been explored to learn disentangled content and speaker
representations [9, 10, 11]. Several studies introduce information bottleneck
[10, 12, 13, 14], instance normalization [15, 16] and mutual information [11,
17] to achieve better disentanglement between content and speaker components.
We note that most VC studies primarily focus on neutral speech, overlooking
the prosodic variations manifested in different emotions. Expressive VC aims
to fill this gap.

Figure 1: A comparison between the conversion of different speech components
across VC[6], expressive VC[3], and emotional VC[18].

(a) Speaker representations with emotion cues

(b) Speaker-independent emotion representations

Figure 2: Visualization of speaker representations and speaker-independent
emotion representations of 50 randomly selected utterances from 4 speakers
across 4 different emotional states. Each point represents one expressive
utterance and the legend indicates speaker identity and emotional state
information.

It remains a formidable undertaking to model emotional styles for expressive
VC due to the hierarchical structure, inherent complexity, subjective nature,
and variability of human emotional speech [3]. Moreover, emotional style
contains both speaker-independent and speaker-dependent features [19, 20, 21].
Speaker-independent emotional features often exhibit consistent associations
with the same emotion across different speakers [22], while speaker-dependent
ones are manifested in speaker characteristics and are unique to each speaker
[19]. Despite their significant roles in conveying emotional style, speaker-
dependent emotional features have been overlooked in previous VC studies [5,
4]. In this paper, we aim to enhance the expressiveness of converted speech
for VC systems by exploring both speaker-dependent and speaker-independent
emotional features.

Another limitation of previous approaches [5, 4, 3] of expressive voice
conversion stems from their reliance on vocoders to reconstruct speech
waveform from acoustic features. Consequently, the quality of the synthesized
speech is intricately influenced by vocoder’s performance [23], which
typically requires a substantial amount of high-quality speech data for
training. To address these issues, we propose a fully end-to-end diffusion-
based framework for expressive VC for the first time. We also demonstrate the
flexibility of our model on any-to-any conversion.

In this paper, we introduce an expressive voice conversion framework, namely
DEVC, which enables effective speaker identity conversion for emotionally
expressive speakers. To achieve this, we employ three encoders dedicated to
handling content representation, speaker representation incorporating speaker-
dependent emotional cues, and speaker-independent emotion representation,
respectively. Content conditioning is facilitated through speech units, while
deep features derived from speaker verification (SV) and speech emotion
recognition (SER) tasks are utilized to capture speaker-dependent and speaker-
independent emotional information. Furthermore, we employ a conditional
denoising diffusion probabilistic model (DDPM) to iteratively reconstruct the
waveform from Gaussian noise, conditioned on the content, speaker, and emotion
representations. Our main contributions include:

  * •

We propose a fully end-to-end expressive voice conversion framework based on a
conditional diffusion model without the need for large-scale training data and
manual annotations;

  * •

Our findings reveal that speaker embeddings derived from an SV model pre-
trained on neutral data effectively capture speaker-dependent emotional cues,
thereby demonstrating their utility in enhancing expressive voice conversion;
and

  * •

Our proposed framework shows flexibility in identity conversion for both seen
and unseen emotional speakers, achieving any-to-any expressive voice
conversion.

The rest of this paper is organized as follows: Section 2 provides an
introduction to the related work. In Section 3, we undertake a comprehensive
analysis of the emotion and speaker representations employed in our framework.
Section 4 presents the details of our proposed framework. The experimental
results are reported in Section 2. In Section 7, we conclude the study.

##  2 Related Work

###  2.1 Expressive Voice Conversion

Expressive speech is a complex composition of various components, including
speaker identity, emotional style, and linguistic content. Previous models in
expressive voice conversion [5, 4] primarily focus on eliminating redundant
information between emotional style and other speech components, such as
content and speaker identity. This objective has been achieved through
techniques like mutual information loss [5] and prosody filters [4]. In
another work [3], a StarGAN-based model incorporates deep features from an SER
model as emotion representation and a simple one-hot vector as the speaker
representation. However, this simplistic speaker representation may not
adequately capture the nuances of emotional speakers, thus limiting its
ability to achieve any-to-any voice conversion for arbitrary speakers.

Another related task is emotional voice conversion, which focuses on
converting the emotional state of the speaker while preserving the speaker’s
identity [24, 18]. In contrast, expressive voice conversion aims to
simultaneously convert both speaker identity and speech style for speakers
displaying emotional expression [3], as shown in Figure 1. This paper will
concentrate on expressive voice conversion.

###  2.2 Denoising Diffusion Probabilistic Models in Speech Synthesis

Denoising diffusion probabilistic models [25] are a class of diffusion models
that uses a Markov chain to gradually convert a simple distribution into a
complicated data distribution. It contains two processes: a diffusion process
where clean input data is converted to an isotropic Gaussian distribution by
adding noise step by step; and a reverse process where clean input can be
recovered from Gaussian noise by predicting and removing the noise introduced
in each step of the diffusion process. Denoising diffusion probabilistic
models have achieved great success in various speech synthesis tasks, such as
vocoding [26], text-to-speech[27], and VC[28, 29, 30]. We further develop the
idea and introduce an end-to-end framework for expressive voice conversion,
which will be elaborated in Section 4.

###  2.3 Representation Learning with Self-Supervised Models

Collecting labeled speech data is difficult, expensive, and time-consuming
[31, 32, 33]. Self-supervised speech models pre-trained on a large amount of
unlabeled data can learn high-level speech representations without relying on
labeled data [34, 35, 36, 37, 38]. Self-supervised representations have been
proven to be highly effective in various speech synthesis tasks due to their
ability to capture essential speech content information [39, 40]. The majority
of existing VC methods leveraging such representations separate content from
speaker-related information through discretization [41, 42]. To retain more
linguistic content in VC, van Niekerk et al. [43] propose to learn soft speech
units by predicting a distribution over discretized representations, further
improving the intelligibility of converted speech. Inspired by this, we
incorporate soft speech units as content conditioning in our proposed method.

Table 1: Mean Euclidean distance between speaker representations from
different emotion pairs among 4 speakers. The first row and the first column
respectively represent emotion states of two groups, with each cell indicating
a unique emotion pair. Note that speakers here are the same set of speakers
from Figure 2.

| Angry | Happy | Neutral | Sad
---|---|---|---|---
Angry | 0.670 | 0.724 | 0.761 | 0.739
Happy | 0.722 | 0.719 | 0.773 | 0.754
Neutral | 0.753 | 0.765 | 0.676 | 0.703
Sad | 0.729 | 0.751 | 0.705 | 0.667

(a) Speaker: 0013 (M1)

| Angry | Happy | Neutral | Sad
---|---|---|---|---
Angry | 0.705 | 0.774 | 0.811 | 0.867
Happy | 0.784 | 0.669 | 0.761 | 0.783
Neutral | 0.832 | 0.754 | 0.693 | 0.735
Sad | 0.873 | 0.776 | 0.733 | 0.615

(b) Speaker: 0020 (M2)

| Angry | Happy | Neutral | Sad
---|---|---|---|---
Angry | 0.699 | 0.803 | 0.782 | 0.955
Happy | 0.811 | 0.751 | 0.783 | 0.904
Neutral | 0.787 | 0.774 | 0.726 | 0.896
Sad | 0.952 | 0.896 | 0.890 | 0.682

(c) Speaker: 0016 (F1)

| Angry | Happy | Neutral | Sad
---|---|---|---|---
Angry | 0.672 | 0.745 | 0.725 | 0.808
Happy | 0.756 | 0.700 | 0.732 | 0.758
Neutral | 0.724 | 0.733 | 0.668 | 0.780
Sad | 0.814 | 0.752 | 0.772 | 0.695

(d) Speaker: 0018 (F2)

##  3 Emotion and Speaker Representations: A Novel Analysis for Expressive VC

Expressive speech introduces variations in the acoustic features, resulting in
increased complexity for speaker identity among emotional speakers [44]. In
expressive voice conversion, it is expected that speaker identity contains
both speaker-dependent emotional cues and other speaker-related information.
Speaker-dependent emotional cues refer to the emotional information that is
specific to an individual speaker. These cues include variations in
intonation, rhythm, and voice quality that reflect the speaker’s emotional
state.

To incorporate such speaker-dependent emotional cues into our expressive voice
conversion model, we extract speaker representations from a pre-trained SV
model [45]. These speaker representations serve as conditioning information
for our expressive voice conversion model, allowing it to generate converted
speech that preserves both target speaker identity and their emotional
nuances.

By employing the t-SNE algorithm [46], we visualize the speaker
representations of two female speakers and two male speakers, as shown in Fig.
2(a). We observe that the speaker representations form separate clusters for
each speaker, indicating successful differentiation of speaker identities. We
also notice the formation of smaller clusters within the same speaker and
emotional state, indicating the presence of speaker-dependent emotional
information.

To further evaluate the similarity of speaker representations, we perform a
comprehensive analysis. We randomly selected 240 emotional speech utterances
from the ESD dataset [18] and divided them equally into two groups. We
calculate the Euclidean distance between utterance pairs from these groups and
obtain the mean distance of emotion pairs. As presented in Table 3(d), the
results consistently demonstrate that the distances within the same emotional
state are lower compared to those between different emotions. This consistent
pattern strongly supports our hypothesis that the speaker representations not
only capture speaker-related information but also encompass speaker-dependent
emotional cues.

In our proposed method, we adopt deep features obtained from a pre-trained
speaker-independent SER model [47] as speaker-independent emotion
representations, following the approach in [48]. To visually analyze the
speaker-independent emotion representations, we employ the t-SNE algorithm
[46] and present the results in Figure 2(b). The visualization demonstrates
that the emotion representations derived from different emotions form well-
separated clusters. Within each cluster, representations from the same emotion
but different speakers tend to be mixed up. This observation confirms that our
emotion representations encompass speaker-independent emotional style
information, allowing for generalization across different individuals.

Our analysis provides valuable insights into the nature of our emotion and
speaker representations. By incorporating these representations, our framework
enables a single DDPM to synthesize expressive speech for any given speaker by
effectively modeling both common characteristics of each emotional state
across different speakers and emotional nuances associated with each
individual. Further details on this approach will be provided in Section 4.

Figure 3: An illustration of the training phase of the proposed DEVC, where
the green boxes represent the modules that are involved in the training while
the others are not.

##  4 Proposed Method - DEVC

In this section, we introduce our fully end-to-end expressive VC model,
Diffusion-based Expressive Voice Conversion (DEVC). DEVC consists of a content
encoder, a speaker encoder, a speaker-independent emotion encoder, and a
diffusion-based decoder, as illustrated in Figure 3. These three encoders
extract content representations, speaker representation with emotional cues,
and speaker-independent emotion representations from input expressive speech,
respectively. Utilizing these auxiliary representations as conditions, the
diffusion-based decoder iteratively generates the converted expressive speech,
starting from Gaussian noise.

###  4.1 Training Stage

Given an expressive speech utterance with $S$ audio segments, the content
encoder extracts a 256-dimensional content representation
$R_{c}\in\mathbbm{R}^{S\times 256}$ from the waveform. The speaker encoder
extracts a $256$-dimension utterance-level speaker representation ${R_{s}}$
with emotional cues, as reported in Table 3(d) and Figure 2(a). The emotion
encoder extracts an utterance-level speaker-independent emotional style
representation $R_{e}$ with 128 dimensions. To align with $R_{c}$, we upsample
the two utterance-level $R_{s}$ and $R_{e}$ to
$R^{\prime}_{s}\in\mathbbm{R}^{S\times 256}$ and
$R^{\prime}_{e}\in\mathbbm{R}^{S\times 128}$ by repeating them, respectively.
We formulate the conditioning $c$ by segment-wise concatenating $R_{c}$ with
$R^{\prime}_{s}$ and $R^{\prime}_{e}$:

$c={\rm concat}(R_{c},R^{\prime}_{s},R^{\prime}_{e})$ (1)

The conditional DDPM-based decoder has two subprocesses: the diffusion process
and the reverse process. Given a segment of waveform $x_{0}$, the diffusion
process is the process of $x_{0}$ being gradually corrupted to Gaussian noise
$x_{T}$ within finite $T$ steps. Assume $\\{x_{1},...,x_{T-1}\\}$ is a
sequence of latent variables where $x_{t}$ is transformed from $x_{t-1}$ by
adding Gaussian noises at each timestep $t\in[1,T]$:

$q(x_{t}|x_{t-1})=\mathcal{N}(x_{t};\sqrt{1-\beta_{t}}x_{t-1},\beta_{t}I)\\\ $
(2)

where $\\{\beta_{1},...\beta_{T}\\}$ is a fixed variance schedule. Given clean
data $x_{0}$, sampling of $x_{t}$ can be written in a closed form:

$\begin{gathered}q(x_{t}|x_{0})=\mathcal{N}(x_{t};\sqrt{\bar{\alpha}_{t}}x_{0},(1-\bar{\alpha_{t}})I)\\\
x_{t}=\sqrt{\bar{\alpha}_{t}}x_{0}+\sqrt{1-\bar{\alpha}_{t}}\epsilon\end{gathered}$
(3)

where $\alpha_{t}=1-\beta_{t}$ and
$\bar{\alpha_{t}}=\prod_{s=1}^{t}\alpha_{s}$. Noise
$\epsilon\sim\mathcal{N}(0,I)$ has the same dimensionality as data $x_{0}$ and
latent variables $x_{1},...x_{T}$.

The reverse process generates a reverse sequence by sampling the posteriors
$q(x_{t-1}|x_{t})$, starting from a Gaussian noise sample $x_{T}$. However,
since $q(x_{t-1}|x_{t})$ is intractable, the decoder learn parameterized
Gaussian transitions $p_{\theta}(x_{t-1}|x_{t})$ with a learned mean
$\mu_{\theta}(x_{t},t,c)$ and a fixed variance $\sigma^{2}_{t}I$[25]:

$p_{\theta}(x_{t-1}|x_{t})=\mathcal{N}(\mu_{\theta}(x_{t},t,c),\sigma^{2}_{t}I)$
(4)

where $\mu_{\theta}(x_{t},t,c)$ is the function of a noise approximator
$\epsilon_{\theta}(x_{t},t,c)$.

Based on speech-conditioning pairs, we then learn the conditional diffusion-
based decoder via:

$L_{\theta}=\|\epsilon-\epsilon_{\theta}(x_{t},t,c)\|^{2}_{2}$ (5)

where $\epsilon\sim\mathcal{N}(0,I)$ is the noise and $\epsilon_{\theta}$
denotes the decoder with learnable parameters $\theta$.

###  4.2 Run-time Conversion

At run-time, DEVC takes a source utterance from one speaker and a reference
utterance from a target speaker, which is either seen or unseen during
training, as shown in Figure 4. The encoders extract content representation
and speaker-independent emotional style representation from the source
utterance and speaker representation with emotional cues from the reference
utterance and the diffusion-based decoder generates converted expressive
speech from Gaussian noise and timestep conditioned on those representations.
We note that emotion representations are derived from a speaker-independent
emotion recognizer, thereby inherently assumed to be speaker-independent (as
shown in Section 3). Through our experiments, we illustrate that the emotion
encoder can take input from either the source utterance or the reference.

Figure 4: An illustration of the run-time phase of the proposed DEVC,
highlighting that waveform generation occurs without the use of a vocoder.

##  5 Experiments

Table 2: Objective evaluations results for the conversion between speakers
seen during training (S2S), the conversion for seen to unseen speakers (S2U),
and the conversion between unseen speakers (U2U). MCD is measured in dB.

| Neutral | Angry | Happy | Sad
---|---|---|---|---
|  MCD$\downarrow$ |  SV$\uparrow$ |  FFE$\downarrow$ |  VDE$\downarrow$ |  MCD$\downarrow$ |  SV$\uparrow$ |  FFE$\downarrow$ |  VDE$\downarrow$ |  MCD$\downarrow$ |  SV$\uparrow$ |  FFE$\downarrow$ |  VDE$\downarrow$ |  MCD$\downarrow$ |  SV$\uparrow$ |  FFE$\downarrow$ |  VDE$\downarrow$
|  Baseline[3] | 8.86 | 0.72 | 0.39 | 0.34 | 8.79 | 0.75 | 0.39 | 0.34 | 8.84 | 0.72 | 0.47 | 0.37 | 9.08 | 0.73 | 0.41 | 0.34
S2S | DEVC | 8.18 | 0.84 | 0.32 | 0.29 | 7.99 | 0.87 | 0.33 | 0.30 | 8.42 | 0.82 | 0.41 | 0.34 | 8.01 | 0.80 | 0.34 | 0.30
| Baseline-U | 8.81 | 0.68 | 0.42 | 0.37 | 8.74 | 0.53 | 0.46 | 0.39 | 8.54 | 0.71 | 0.47 | 0.39 | 8.51 | 0.67 | 0.43 | 0.35
S2U | DEVC | 8.80 | 0.80 | 0.34 | 0.30 | 8.61 | 0.71 | 0.39 | 0.34 | 8.43 | 0.73 | 0.42 | 0.36 | 8.12 | 0.76 | 0.39 | 0.33
| Baseline-U | 9.02 | 0.65 | 0.42 | 0.37 | 9.26 | 0.69 | 0.47 | 0.40 | 8.72 | 0.70 | 0.48 | 0.39 | 8.62 | 0.67 | 0.45 | 0.38
U2U | DEVC | 8.67 | 0.79 | 0.35 | 0.31 | 8.62 | 0.84 | 0.38 | 0.33 | 8.39 | 0.79 | 0.41 | 0.35 | 8.19 | 0.77 | 0.40 | 0.35

###  5.1 Dataset

We evaluate the proposed DEVC on the ESD dataset[18]. For each emotion and
each speaker, there are 300 training utterances, 20 reference utterances, and
30 test utterances. In our experiments, we select four emotions (Neutral,
Angry, Happy, and Sad) and four speakers (0013, 0016, 0018, and 0020) as seen
speakers for training, while randomly selecting four other speakers as unseen
speakers. We evaluate our model in three scenarios: the conversion between
seen speakers (S2S), the conversion from seen speakers to unseen speakers
(S2U), and the conversion between unseen speakers (U2U), commonly referred to
as the any-to-any task. All speech samples are sampled at 16 kHz and saved in
16-bit format.

###  5.2 Experimental Setup

DEVC consists of a content encoder, a speaker encoder, an emotion encoder, and
a diffusion-based decoder. The content encoder is based on the HuBERT-Base
backbone network, augmented with a linear layer [43] and pre-trained on the
LibriSpeech-960 dataset [49]. The emotion encoder consists of a three-
dimensional CNN layer, a BLSTM layer, an attention layer, and a fully-
connected (FC) layer. It is initially pre-trained on the IEMOCAP dataset [50]
and subsequently fine-tuned using the ESD dataset [48]. The speaker encoder
follows the architecture in [45], utilizing a 3-layer LSTM with projection to
generate 256-dimensional representations.

For the diffusion decoder, the architecture is similar to that described in
[26]. The step encoder comprises a 128-dimensional sinusoidal position
encoding and two FC layers with Switch activations. The PreNet consists of a
Conv1$\times$1 layer followed by ReLU activation. The output of the step
encoder and latent variables are then fed into a stack of 64 residual blocks
with 128 residual channels. The skip connections from all residual layers are
summed before entering the PostNet, which includes two Conv1$\times$1 layers.
DEVC is trained on a single Nvidia 3080Ti GPU for 1.2 million steps using Adam
optimizer with a batch size of 16 and a learning rate of 0.0002.

###  5.3 Baselines

In a comparative study, we adopt the following two expressive VC models as our
baseline:

  * •

Baseline: JES-StarGAN [3], a many-to-many expressive voice conversion
framework for S2S scenarios; and

  * •

Baseline-U: We replace the one-hot speaker label in JES-StarGAN [3] with
speaker representations for S2U and U2U settings.

Table 3: The MOS results with 95% confidence interval to evaluate the speech
quality of converted speech.

Framework | S2S | S2U | U2U
---|---|---|---
Baseline[3] |  2.12$\pm$0.15 | NA | NA
Baseline-U | NA |  2.45$\pm$0.16 |  2.46$\pm$0.18
DEVC |  3.28$\pm$0.08 |  3.16$\pm$0.12 |  3.27$\pm$0.13

Table 4: Ablation studies results in S2S setting.

Model |  MCD [dB] $\downarrow$ |  SV $\uparrow$ |  F0-RMSE [Hz] $\downarrow$
---|---|---|---
$R_{c}$, $R_{s}$ | 10.81 | 0.71 | 82.38
$R_{c}$, $R_{e}$ | 8.34 | 0.79 | 55.62
$\bm{R_{c}}$, $\bm{R_{e}}$, $\bm{R_{s}}$ | 8.15 | 0.85 | 53.98

###  5.4 Objective Evaluations

Table 2 presents the objective evaluation results of our proposed framework
and the baselines. Mel-cepstral distortion (MCD)[51] is utilized to measure
spectral distortion between synthesized samples and target samples. Lower MCD
values indicate better quality and closer similarity between the synthesized
samples and the target samples [52]. Our proposed method outperforms the
baselines in all settings (S2S, S2U, and U2U) by achieving consistently lower
MCD values, which indicates a better conversion performance.

To evaluate speaker similarity, we utilize a pre-trained SV
model111https://github.com/resemble-ai/Resemblyzer and report the SV accuracy
in Table 2. Higher SV accuracy indicates better conversion of speaker identity
[31]. Our proposed framework achieves higher speaker similarity compared to
the baselines, highlighting its capability to accurately convert the
characteristics of emotional speakers.

In addition, we employ two pitch similarity metrics: Voicing Decision Error
(VDE) [53] and F0 Frame Error (FFE) [54]. VDE measures the proportion of
voiced/unvoiced decision difference between the converted and target
utterances, while FFE captures both pitch similarity and voiced/unvoiced
decision differences [42]. Our proposed method outperforms the baselines,
indicating its effectiveness in capturing and reproducing the desired
characteristics in the converted speech.

###  5.5 Subjective Evaluation

We conduct listening experiments with 13 subjects to assess the speech
quality, speaker similarity, and emotional style similarity. Each of them
listens to 192 converted utterances in total. In the Mean Opinion Score (MOS)
test, participants evaluated the converted utterances using a rating scale
ranging from 1 (poor) to 5 (excellent) to assess speech quality. The results
for S2S, S2U, and U2U are presented in Table 3. As shown in Table 3, DEVC
outperforms the baselines in terms of speech quality.

In addition, we perform ABX preference tests to evaluate speaker similarity in
the S2S, S2U, and U2U settings. The listeners compare the converted utterances
and select the utterance that sounds closer to the reference in terms of
speaker identity. The results, illustrated in Figure 5(a), demonstrate that
our proposed framework significantly outperforms the baseline, indicating that
our method excels in accurately converting speaker identities for both seen
and unseen speakers.

We further conduct ABX preference tests to evaluate the emotional style
similarity in S2S, S2U, and U2U settings. These tests require listeners to
compare the emotional styles of the reference and converted utterances. Figure
5(b) presents the results, clearly demonstrating the superior performance of
our proposed framework over the baseline. These findings highlight the
capability of our method to successfully convert emotional styles for
arbitrary speakers, further validating the effectiveness of our emotional
style modeling approach.

It is worth mentioning that the conversion between unseen speakers to unseen
speakers is particularly challenging due to the lack of training data.
However, even in this challenging scenario, our proposed method achieves
excellent results in subjective evaluation, further affirming its ability to
conduct high-quality conversion in any-to-any setting.

(a) Speaker Similarity

(b) Emotional Style Similarity

Figure 5: ABX preference results for S2S, S2U, and U2S settings to evaluate:
(a) speaker similarity and (b) emotional style similarity. We used Baseline
(JES-StarGAN) for S2S, and Baseline-U for S2U and U2S settings.

###  5.6 Ablation studies

We conduct ablation studies to analyze the impact of our speaker
representations with emotion cues and speaker-independent emotion
representations in the context of expressive voice conversion. Specifically,
we investigate three different models, all sharing the same diffusion-based
decoder architecture: the first model includes content representations and
speaker representations, denoted as $R_{c},R_{s}$. The second model
incorporates content representations, speaker-independent emotional style
representations, and one-hot speaker labels, denoted as $R_{c},R_{e}$. The
third model represents our proposed method, incorporating content
representations, speaker representations with emotion cues, and speaker-
independent emotion representations, denoted as $R_{c},R_{e},R_{s}$.

In the S2S setting, we evaluate the performance of these models by calculating
metrics MCD, SV accuracy, and F0 root mean squared errors (F0-RMSE)[55], as
presented in Table 4. From the results, our proposed method, which utilizes
speaker representations instead of one-hot speaker labels, exhibits superior
performance compared to the model incorporating one-hot speaker vectors. This
highlights the advantage of our speaker representation with emotion cues,
indicating that it captures and represents speaker-dependent emotion cues more
effectively. In addition, we observe that our proposed method outperforms the
other models. This finding suggests the effectiveness of our emotional style
modeling method in the expressive voice conversion task.

##  6 Discussion

We note that our proposed method has the potential to be extended to convert
emotion states for different speakers with various reference utterance
settings. However, in this paper, our focus is solely on expressive VC. We
intend to broaden the scope of this work to include the conversion of emotion
states for the same speaker (emotional VC) or different speakers in our future
research.

##  7 Conclusion

In this paper, we present a novel framework for expressive voice conversion
using a diffusion-based approach that supports any-to-any conversion. Our
framework leverages speech units as content representations and incorporates
deep features extracted from speech emotion recognition and speaker
verification tasks to capture emotion and speaker characteristics. A key
finding of our research is that speaker embeddings obtained from a pre-trained
SV model using neutral data inherently contain speaker-dependent emotional
features, thereby proving advantageous for the expressive voice conversion
task. Our proposed framework exhibits remarkable flexibility in converting
both seen and unseen speakers without the need for vocoders, marking it as the
first end-to-end diffusion model-based expressive voice conversion framework
known to us.

## References

  * [1] Dagmar M Schuller and Björn W Schuller,  “A review on five recent and near-future developments in computational processing of emotion in the human voice,”  Emotion Review, p. 1754073919898526, 2020.
  * [2] Kun Zhou, Berrak Sisman, Rajib Rana, Björn W. Schuller, and Haizhou Li,  “Speech synthesis with mixed emotions,”  IEEE Transactions on Affective Computing, pp. 1–16, 2022.
  * [3] Zongyang Du, Berrak Sisman, Kun Zhou, and Haizhou Li,  “Expressive voice conversion: A joint framework for speaker identity and emotional style transfer,”  in 2021 IEEE Automatic Speech Recognition and Understanding Workshop (ASRU). IEEE, 2021, pp. 594–601.
  * [4] Wendong Gan, Bolong Wen, Ying Yan, Haitao Chen, Zhichao Wang, Hongqiang Du, Lei Xie, Kaixuan Guo, and Hai Li,  “Iqdubbing: Prosody modeling based on discrete self-supervised speech representation for expressive voice conversion,”  arXiv preprint arXiv:2201.00269, 2022.
  * [5] Zongyang Du, Berrak Sisman, Kun Zhou, and Haizhou Li,  “Disentanglement of Emotional Style and Speaker Identity for Expressive Voice Conversion,”  in Proc. Interspeech 2022, 2022, pp. 2603–2607.
  * [6] Berrak Sisman, Junichi Yamagishi, Simon King, and Haizhou Li,  “An overview of voice conversion and its challenges: From statistical modeling to deep learning,”  IEEE/ACM Transactions on Audio, Speech, and Language Processing, vol. 29, pp. 132–157, 2021.
  * [7] Tomoki Toda, Alan W Black, and Keiichi Tokuda,  “Voice conversion based on maximum-likelihood estimation of spectral parameter trajectory,”  IEEE Transactions on Audio, Speech, and Language Processing, vol. 15, no. 8, pp. 2222–2235, 2007.
  * [8] Ryo Aihara, Ryoichi Takashima, Tetsuya Takiguchi, and Yasuo Ariki,  “Gmm-based emotional voice conversion using spectrum and prosody features,”  American Journal of Signal Processing, vol. 2, no. 5, pp. 134–138, 2012.
  * [9] Songxiang Liu, Jinghua Zhong, Lifa Sun, Xixin Wu, Xunying Liu, and Helen Meng,  “Voice conversion across arbitrary speakers based on a single target-speaker utterance.,”  in Interspeech, 2018, pp. 496–500.
  * [10] Kaizhi Qian, Yang Zhang, Shiyu Chang, Xuesong Yang, and Mark Hasegawa-Johnson,  “Autovc: Zero-shot voice style transfer with only autoencoder loss,”  in International Conference on Machine Learning. PMLR, 2019, pp. 5210–5219.
  * [11] Disong Wang, Liqun Deng, Yu Ting Yeung, Xiao Chen, Xunying Liu, and Helen Meng,  “VQMIVC: Vector Quantization and Mutual Information-Based Unsupervised Speech Representation Disentanglement for One-Shot Voice Conversion,”  in Proc. Interspeech 2021, 2021, pp. 1344–1348.
  * [12] Andros Tjandra, Berrak Sisman, Mingyang Zhang, Sakriani Sakti, Haizhou Li, and Satoshi Nakamura,  “Vqvae unsupervised unit discovery and multi-scale code2spec inverter for zerospeech challenge 2019,”  Proc. Interspeech 2019, pp. 1118–1122, 2019.
  * [13] Kazuhiro Kobayashi, Wen-Chin Huang, Yi-Chiao Wu, Patrick Lumban Tobing, Tomoki Hayashi, and Tomoki Toda,  “crank: An open-source software for nonparallel voice conversion based on vector-quantized variational autoencoder,”  in ICASSP 2021-2021 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). IEEE, 2021, pp. 5934–5938.
  * [14] Ziyi Chen and Pengyuan Zhang,  “Tvqvc: Transformer based vector quantized variational autoencoder with ctc loss for voice conversion.,”  in Interspeech, 2021, pp. 826–830.
  * [15] Yen-Hao Chen, Da-Yi Wu, Tsung-Han Wu, and Hung-yi Lee,  “Again-vc: A one-shot voice conversion using activation guidance and adaptive instance normalization,”  in ICASSP 2021-2021 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). IEEE, 2021, pp. 5954–5958.
  * [16] Ju chieh Chou and Hung-Yi Lee,  “One-Shot Voice Conversion by Separating Speaker and Content Representations with Instance Normalization,”  in Proc. Interspeech 2019, 2019, pp. 664–668.
  * [17] Liangjie Huang, Tian Yuan, Yunming Liang, Zeyu Chen, Can Wen, Yanlu Xie, Jinsong Zhang, and Dengfeng Ke,  “Limi-vc: A light weight voice conversion model with mutual information disentanglement,”  in ICASSP 2023 - 2023 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), 2023, pp. 1–5.
  * [18] Kun Zhou, Berrak Sisman, Rui Liu, and Haizhou Li,  “Emotional voice conversion: Theory, databases and esd,”  Speech Communication, vol. 137, pp. 1–18, 2022.
  * [19] Kusha Sridhar and Carlos Busso,  “Unsupervised personalization of an emotion recognition system: The unique properties of the externalization of valence in speech,”  IEEE Transactions on Affective Computing, vol. 13, no. 4, pp. 1959–1972, 2022.
  * [20] Kun Zhou, Berrak Sisman, and Haizhou Li,  “Vaw-gan for disentanglement and recomposition of emotional elements in speech,”  in 2021 IEEE Spoken Language Technology Workshop (SLT). IEEE, 2021, pp. 415–422.
  * [21] Kun Zhou, Berrak Sisman, Rajib Rana, Björn W. Schuller, and Haizhou Li,  “Emotion intensity and its control for emotional voice conversion,”  IEEE Transactions on Affective Computing, vol. 14, no. 1, pp. 31–48, 2023.
  * [22] Kun Zhou, Berrak Sisman, Mingyang Zhang, and Haizhou Li,  “Converting anyone’s emotion: Towards speaker-independent emotional voice conversion,”  Proc. Interspeech 2020, pp. 3416–3420, 2020.
  * [23] Bac Nguyen and Fabien Cardinaux,  “Nvc-net: End-to-end adversarial voice conversion,”  in ICASSP 2022 - 2022 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), 2022, pp. 7012–7016.
  * [24] Nirmesh Shah, Mayank Singh, Naoya Takahashi, and Naoyuki Onoe,  “Nonparallel emotional voice conversion for unseen speaker-emotion pairs using dual domain adversarial network & virtual domain pairing,”  in ICASSP 2023 - 2023 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), 2023, pp. 1–5.
  * [25] Jonathan Ho, Ajay Jain, and Pieter Abbeel,  “Denoising diffusion probabilistic models,”  Advances in Neural Information Processing Systems, vol. 33, pp. 6840–6851, 2020.
  * [26] Zhifeng Kong, Wei Ping, Jiaji Huang, Kexin Zhao, and Bryan Catanzaro,  “Diffwave: A versatile diffusion model for audio synthesis,”  in International Conference on Learning Representations, 2021.
  * [27] Rongjie Huang, Max W. Y. Lam, Jun Wang, Dan Su, Dong Yu, Yi Ren, and Zhou Zhao,  “Fastdiff: A fast conditional diffusion model for high-quality speech synthesis,”  in Proceedings of the Thirty-First International Joint Conference on Artificial Intelligence, IJCAI-22, Lud De Raedt, Ed. 7 2022, pp. 4157–4163, International Joint Conferences on Artificial Intelligence Organization,  Main Track.
  * [28] Songxiang Liu, Yuewen Cao, Dan Su, and Helen Meng,  “Diffsvc: A diffusion probabilistic model for singing voice conversion,”  in 2021 IEEE Automatic Speech Recognition and Understanding Workshop (ASRU). IEEE, 2021, pp. 741–748.
  * [29] Vadim Popov, Ivan Vovk, Vladimir Gogoryan, Tasnima Sadekova, Mikhail Sergeevich Kudinov, and Jiansheng Wei,  “Diffusion-based voice conversion with fast maximum likelihood sampling scheme,”  in International Conference on Learning Representations, 2021.
  * [30] Tasnima Sadekova, Vladimir Gogoryan, Ivan Vovk, Vadim Popov, Mikhail A Kudinov, and Jiansheng Wei,  “A unified system for voice cloning and voice conversion through diffusion probabilistic modeling.,”  in INTERSPEECH, 2022, pp. 3003–3007.
  * [31] Jheng-hao Lin, Yist Y Lin, Chung-Ming Chien, and Hung-yi Lee,  “S2vc: A framework for any-to-any voice conversion with self-supervised pretrained representations,”  arXiv preprint arXiv:2104.02901, 2021.
  * [32] Abdelrahman Mohamed, Hung-yi Lee, Lasse Borgholt, Jakob D Havtorn, Joakim Edin, Christian Igel, Katrin Kirchhoff, Shang-Wen Li, Karen Livescu, Lars Maaløe, et al.,  “Self-supervised speech representation learning: A review,”  IEEE Journal of Selected Topics in Signal Processing, 2022.
  * [33] Szu-Jui Chen, Jiamin Xie, and John H.L. Hansen,  “FeaRLESS: Feature Refinement Loss for Ensembling Self-Supervised Learning Features in Robust End-to-end Speech Recognition,”  in Proc. Interspeech 2022, 2022, pp. 3058–3062.
  * [34] Alexei Baevski, Steffen Schneider, and Michael Auli,  “vq-wav2vec: Self-supervised learning of discrete speech representations,”  arXiv preprint arXiv:1910.05453, 2019.
  * [35] Steffen Schneider, Alexei Baevski, Ronan Collobert, and Michael Auli,  “wav2vec: Unsupervised pre-training for speech recognition,”  in Proc. Interspeech 2019, pp. 3465–3469.
  * [36] Alexei Baevski, Yuhao Zhou, Abdelrahman Mohamed, and Michael Auli,  “wav2vec 2.0: A framework for self-supervised learning of speech representations,”  in Proc. NeurIPS, 2020, vol. 33, pp. 12449–12460.
  * [37] Wei-Ning Hsu, Benjamin Bolte, Yao-Hung Hubert Tsai, Kushal Lakhotia, Ruslan Salakhutdinov, and Abdelrahman Mohamed,  “Hubert: Self-supervised speech representation learning by masked prediction of hidden units,”  IEEE/ACM Transactions on Audio, Speech, and Language Processing, vol. 29, pp. 3451–3460, 2021.
  * [38] Mu Yang, Ram C. M. C. Shekar, Okim Kang, and John H. L. Hansen,  “What Can an Accent Identifier Learn? Probing Phonetic and Prosodic Information in a Wav2vec2-based Accent Identification Model,”  in Proc. INTERSPEECH 2023, 2023, pp. 1923–1927.
  * [39] Felix Kreuk, Adam Polyak, Jade Copet, Eugene Kharitonov, Tu Anh Nguyen, Morgan Rivière, Wei-Ning Hsu, Abdelrahman Mohamed, Emmanuel Dupoux, and Yossi Adi,  “Textless speech emotion conversion using discrete & decomposed representations,”  in Proceedings of the 2022 Conference on Empirical Methods in Natural Language Processing. 2022, pp. 11200–11214, Association for Computational Linguistics.
  * [40] Junchen Lu, Berrak Sisman, Mingyang Zhang, and Haizhou Li,  “High-Quality Automatic Voice Over with Accurate Alignment: Supervision through Self-Supervised Discrete Speech Units,”  in Proc. INTERSPEECH 2023, 2023, pp. 5536–5540.
  * [41] Wen-Chin Huang, Yi-Chiao Wu, and Tomoki Hayashi,  “Any-to-one sequence-to-sequence voice conversion using self-supervised discrete speech representations,”  in ICASSP 2021-2021 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). IEEE, 2021, pp. 5944–5948.
  * [42] Adam Polyak, Yossi Adi, Jade Copet, Eugene Kharitonov, Kushal Lakhotia, Wei-Ning Hsu, Abdelrahman Mohamed, and Emmanuel Dupoux,  “Speech Resynthesis from Discrete Disentangled Self-Supervised Representations,”  in Proc. Interspeech 2021, 2021.
  * [43] Benjamin van Niekerk, Marc-André Carbonneau, Julian Zaïdi, Matthew Baas, Hugo Seuté, and Herman Kamper,  “A comparison of discrete and soft speech units for improved voice conversion,”  in ICASSP 2022-2022 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). IEEE, 2022, pp. 6562–6566.
  * [44] Srinivas Parthasarathy, Chunlei Zhang, John HL Hansen, and Carlos Busso,  “A study of speaker verification performance with expressive speech,”  in 2017 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). IEEE, 2017, pp. 5540–5544.
  * [45] Li Wan, Quan Wang, Alan Papir, and Ignacio Lopez Moreno,  “Generalized end-to-end loss for speaker verification,”  in 2018 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). IEEE, 2018, pp. 4879–4883.
  * [46] Laurens Van der Maaten and Geoffrey Hinton,  “Visualizing data using t-sne.,”  Journal of machine learning research, vol. 9, no. 11, 2008.
  * [47] Mingyi Chen, Xuanji He, Jing Yang, and Han Zhang,  “3-d convolutional recurrent neural networks with attention model for speech emotion recognition,”  IEEE Signal Processing Letters, vol. 25, no. 10, pp. 1440–1444, 2018.
  * [48] Kun Zhou, Berrak Sisman, Rui Liu, and Haizhou Li,  “Seen and unseen emotional style transfer for voice conversion with a new emotional speech dataset,”  in ICASSP 2021-2021 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). IEEE, 2021, pp. 920–924.
  * [49] Vassil Panayotov, Guoguo Chen, Daniel Povey, and Sanjeev Khudanpur,  “Librispeech: An asr corpus based on public domain audio books,”  in 2015 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), 2015, pp. 5206–5210.
  * [50] Carlos Busso, Murtaza Bulut, Chi-Chun Lee, Abe Kazemzadeh, Emily Mower, Samuel Kim, Jeannette N Chang, Sungbok Lee, and Shrikanth S Narayanan,  “Iemocap: Interactive emotional dyadic motion capture database,”  Language resources and evaluation, vol. 42, pp. 335–359, 2008.
  * [51] Robert Kubichek,  “Mel-cepstral distance measure for objective speech quality assessment,”  in Proceedings of IEEE Pacific Rim Conference on Communications Computers and Signal Processing. IEEE, 1993, vol. 1, pp. 125–128.
  * [52] Zongyang Du, Kun Zhou, Berrak Sisman, and Haizhou Li,  “Spectrum and prosody conversion for cross-lingual voice conversion with cyclegan,”  in 2020 Asia-Pacific Signal and Information Processing Association Annual Summit and Conference (APSIPA ASC), 2020, pp. 507–513.
  * [53] Tomohiro Nakatani, Shigeaki Amano, Toshio Irino, Kentaro Ishizuka, and Tadahisa Kondo,  “A method for fundamental frequency estimation and voicing decision: Application to infant utterances recorded in real acoustical environments,”  Speech Communication, vol. 50, no. 3, pp. 203–214, 2008.
  * [54] Wei Chu and Abeer Alwan,  “Reducing f0 frame error of f0 tracking algorithms under noisy conditions with an unvoiced/voiced classification frontend,”  in 2009 IEEE International Conference on Acoustics, Speech and Signal Processing. IEEE, 2009, pp. 3969–3972.
  * [55] Zhaojie Luo, Jinhui Chen, Tetsuya Takiguchi, and Yasuo Ariki,  “Neutral-to-emotional voice conversion with cross-wavelet transform f0 using generative adversarial networks,”  APSIPA Transactions on Signal and Information Processing, vol. 8, pp. e10, 2019.


# Statistics of Pressure Fluctuations in Turbulent Kinetic Plasmas

S. Adhikari,1,2 W. H. Matthaeus,2,3 T. N. Parashar,2,4 M.A. Shay,2,3 and P. A.
Cassak1,2,3
1 Department of Physics and Astronomy, West Virginia University, Morgantown,
West Virginia 26506, USA
2 Department of Physics and Astronomy, University of Delaware, Newark,
Delaware 19716, USA
3 Bartol Research Institute, Department of Physics and Astronomy,University of
Delaware, Newark, Delaware 19716, USA
4 School of Chemical and Physical Sciences, Victoria University of Wellington,
Wellington 6140, NZ E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

In this study we explore the statistics of pressure fluctuations in kinetic
collisionless turbulence. A $2.5$D kinetic particle-in-cell (PIC) simulation
of decaying turbulence is used to investigate pressure balance via the
evolution of thermal and magnetic pressure in a plasma with $\beta$ of order
unity. We also discuss the behavior of thermal, magnetic and total pressure
structure functions and their corresponding wavenumber spectra. The total
pressure spectrum exhibits a slope of $-7/3$ extending for about a decade in
the ion-inertial range. In contrast, shallower $-5/3$ spectra are
characteristic of the magnetic pressure and thermal pressure. The steeper
total pressure spectrum is a consequence of cancellation caused by density-
magnetic field magnitude anticorrelation. Further, we evaluate higher order
total pressure structure functions in an effort to discuss intermittency and
compare the power exponents with higher order structure functions of velocity
and magnetic fluctuations. Finally, applications to astrophysical systems are
also discussed.

###### keywords:

Pressure spectrum, Structure functions, Intermittency, Pressure balance

††pubyear: 2023††pagerange: Statistics of Pressure Fluctuations in Turbulent
Kinetic Plasmas–Statistics of Pressure Fluctuations in Turbulent Kinetic
Plasmas

##  1 Introduction

###  1.1 Motivation

Most naturally occurring as well as laboratory plasma flows are turbulent.
Astrophysical systems of interest include the intracluster medium (ICM)
(Schuecker et al., 2004; Churazov et al., 2012; Zhuravleva et al., 2014;
Zhuravleva et al., 2019), black hole accretion disks (Balbus & Hawley, 1998;
Pessah, 2010), solar wind (Coleman Jr, 1968; Matthaeus & Goldstein, 1982;
Bruno & Carbone, 2013), and planetary magnetospheres (Sahraoui et al., 2020).
Turbulent stresses and pressure play an important role at scales ranging from
the largest structures in the universe (Simionescu et al., 2019) to laboratory
plasmas. A prominent feature of galactic plasma observations is the “great
powerlaw” spectrum of electron density fluctuations that extends over more
than ten orders of magnitude in scale (Armstrong et al., 1981) in the
interstellar medium. Under mild assumptions (Montgomery et al., 1987) this
spectrum is found to be consistent with Kolmogorov spectrum of magnetic field
fluctuations. This calculation depends on an assumption of near-
incompressibility and also implies that density and magnetic field, and
therefore mechanical and magnetic pressures, are anticorrelated and partially
cancel one another. This kind of statistical pressure balance is an element of
turbulent dynamics, in contrast to the superficially similar static pressure
balances that are special solutions of the ideal MHD equations.

Pressure balance plays an important role in astrophysics. The mass of galaxy
clusters (Kravtsov & Borgani, 2012; Planelles et al., 2015), the largest
structures in the universe, provides stringent constraints on cosmological
models (Voigt & Fabian, 2006; Carlberg et al., 1997; Pratt et al., 2019).
Hydrostatic pressure balance is a common assumption used to estimate the mass
of relaxed clusters (Lau et al., 2013). However, this assumption can easily
break in clusters far from equilibrium and can also be broken by turbulent
stress contributions (Markevitch et al., 2004; Angelinelli et al., 2020).
X-ray (Schuecker et al., 2004; Churazov et al., 2012; Zhuravleva et al., 2014)
and radio observations via the Sunayev-Zel’dovich (Khatri & Gaspari, 2016)
effect provide detailed measurements of the thermal properties of galaxy
clusters. An accurate understanding of the pressure dynamics in the
intracluster medium can help disentangle contributions from thermal pressure
and turbulent stresses.

Similarly, properties of the pressure have a major impact in the heliosphere.
Near-Earth, the solar wind plays a vital role in shaping Earth’s magnetosphere
and driving space weather. The solar wind is measured to be turbulent (Coleman
Jr, 1968; Matthaeus & Goldstein, 1982; Bruno & Carbone, 2013) and is observed
to be hotter than expected from adiabatic cooling (Wang & Richardson, 2001;
Hellinger et al., 2013). Heating implies an increase of the trace of the
pressure tensors, and this increase in the internal energy can occur through
interaction between the pressure tensors and the velocity gradient tensors
(Yang et al., 2017a). The solar wind routinely displays intervals of pressure
balance (Vellante & Lazarus, 1987; Goldstein et al., 1995; Reisenfeld et al.,
1999; Ruffolo et al., 2021) which potentially emerge as a result of the nearly
incompressible magnetohydrodynamics (NI-MHD) of the solar wind (Matthaeus et
al., 1991).

In this study, we develop a framework to discuss pressure balance in the
context of the equation governing pressure in incompressible MHD with and
without shear. We also explore various aspects of pressure dynamics in kinetic
plasmas using a fully kinetic particle-in-cell (PIC) simulation. We find the
pressure spectrum follows the hydrodynamic scaling in the incompressible MHD
regime. Although the kinetic plasmas are not fully incompressible, the
pressure spectrum shows an excellent agreement with the power scaling
following the MHD predictions, indicating that the simulation is likely in a
state of NI-MHD. The remaining paper is organized as follows: Section 2 first
reviews some additional theoretical background and then constructs the
Poisson’s equation for pressure in an incompressible MHD framework. In section
3 we report the results from PIC simulation, and finally we discuss the
implications and future directions in section 4.

##  2 Theory

###  2.1 Background

Following the seminal work on hydrodynamic pressure fluctuations by Obukhoff
and Batchelor (Obukhoff & Yaglom, 1951; Batchelor, 1951),the second-order
pressure structure function $D_{P}^{2}$ and pressure spectrum have been
studied both experimentally (Uberoi, 1953; George et al., 1984; Tsuji &
Ishihara, 2003) and via direct numerical simulations (DNS) (Schumann &
Patterson, 1978; Kim & Antonia, 1993; Pullin & Rogallo, 1994; Nelkin & Chen,
1998; Squire et al., 2023), mostly for homogeneous and isotropic turbulence.
While these studies rely on the assumptions that the velocity distribution is
joint gaussian or the small-scale statistics follow Kolmogorov’s similarity
hypothesis, $D_{P}^{2}$ as a function of lag $r$ is found to be proportional
to $r^{4/3}$. Therefore, the pressure spectrum $E_{P}(k)$ varies as $k^{-7/3}$
within the inertial range, where $k$ is the wavenumber. Attempts have also
been made to extend the pressure spectrum to the higher wavenumber
(dissipation range) (Zhao et al., 2016). However, the kinetic/dissipation
range scaling exponents are less well understood, in general and for the
pressure spectrum.

The study of nearly incompressible MHD (Klainerman & Majda, 1981) is a bit
more subtle than the hydro case. Beginning with the Montgomery et al. (1987)
explanation of the $k^{-5/3}$ spectrum of the interstellar density structure
function, several additional layers of theory have been presented (Higdon,
1984; Matthaeus & Brown, 1988a) based on the idea that the core solution of
the plasma dynamics in a perturbation theory controlled by the turbulent Mach
number is incompressible MHD, and the leading order density fluctations are
those discussed by Montgomery et al. (1987). Additional theories have been
presented that elaborate on the role of thermal fluctuations (Zank &
Matthaeus, 1990; Bayly et al., 1992), and variations in plasma beta (Zank &
Matthaeus, 1992, 1993). In light of the latter references, the original
perturbation theory of Montgomery et al. (1987) and Matthaeus & Brown (1988a)
can be seen to apply most readily to plasma beta greater than unity.

To place the current work better in the context of these prior studies, it is
important to emphasize that most of the literature focuses on mechanical
pressure statistics in both hydrodynamic turbulence and MHD. Then mechanical
pressure is linked to density and thermal fluctuations through an equation of
state. In contrast, the present study is mainly concerned with statistics of
the total pressure – magnetic plus mechanical – and contrasting this with
behavior of the magnetic and mechanical pressure separately. Therefore we are
concerned with pressure balances. We will also examine higher order statistics
related to total pressure. Notably we will study in detail these properties
using a particle-in-cell model of collisonless plasma, thus departing
substantially from the starting point of the MHD model.

###  2.2 Pressure Poisson’s equation in MHD

The incompressible MHD equations (using Einstein’s convention) are given by

$\displaystyle\partial_{t}u_{i}+u_{j}\partial_{j}u_{i}-b_{j}\partial_{j}b_{i}=-\partial_{i}P+\nu\partial_{j}\partial_{j}u_{i},$
(1)
$\displaystyle\partial_{t}b_{i}+u_{j}\partial_{j}b_{i}-b_{j}\partial_{j}u_{i}=\eta\partial_{j}\partial_{j}b_{i},$
(2)

where u is the flow velocity, b is the magnetic field in Alfvén units, $P$ is
the total pressure (magnetic plus thermal), $\nu$ is the kinematic viscosity,
and $\eta$ is the electric diffusivity. Here, the density is assumed to be
unity and the velocity is solenoidal $\nabla\cdot{\bf u}=0$. Now, in terms of
Elsässer variables $\textbf{z}^{\pm}$ defined as
$\textbf{z}^{\pm}=\textbf{u}\pm\textbf{b}$, the above equations can be
combined as

$\partial_{t}z_{i}^{\pm}+z_{j}^{\mp}\partial_{j}z_{i}^{\pm}=-\partial_{i}P+\nu\partial_{j}\partial_{j}z_{i}^{\pm},$
(3)

where we assume $\nu=\eta$. Taking the divergence of Eqn. 3 and employing
incompressibility ($\partial_{i}z_{i}=0$) we get

$-\partial_{ii}P=\partial_{i}(z_{j}^{\mp}\partial_{j}z_{i}^{\pm})=(\partial_{i}z_{j}^{\mp})(\partial_{j}z_{i}^{\pm}).$
(4)

As written, the pressure is now seen to be as a constraint force that acts to
maintain the solenoidal nature of the velocity field. The so-called Nearly
Incompressible MHD (NI-MHD) theory proceeds to compute the density as a linear
response to the pressure fluctuations (Batchelor, 1951; Montgomery et al.,
1987; Matthaeus & Brown, 1988b) that emerge from the solution of Eqn. 4.

Incompressibility does not imply constant pressure in the system, as seen in
Eqn. 4; instead, pressure is a non-linear function of the Elsässer variables
suggesting any fluctuation in the Elsässer variable creates a pressure field
in the system and vice-versa. The right hand side of Eqn. 4 can be split into
symmetric and anti-symmetric (skew-symmetric) parts as

$-\partial_{ii}P=(\partial_{i}z_{j}^{\mp})(\partial_{j}z_{i}^{\pm})=(\tau_{ij}^{\mp}+R_{ij}^{\mp})(\tau_{ij}^{\pm}-R_{ij}^{\pm}),$
(5)

where
$\tau_{ij}^{\pm}=\frac{1}{2}(\partial_{i}z_{j}^{\pm}+\partial_{j}z_{i}^{\pm})$
is the symmetric part of the decomposition, and is called the rate of strain
tensor related to the Elsässer’s variable, while
$R_{ij}^{\pm}=\frac{1}{2}(\partial_{i}z_{j}^{\pm}-\partial_{j}z_{i}^{\pm})$ is
the anti-symmetric component, called the rate of rotation tensor associated
with the curl of the Elsässer variable. One can simplify Eqn. 5 further to
obtain

$-\partial_{ii}P=\tau_{ij}^{-}\tau_{ij}^{+}-R_{ij}^{-}R_{ij}^{+}.$ (6)

The first term on the right hand side is related to the strain tensor. If one
decomposes $\tau_{ij}^{\pm}$ such that
$S_{ij}=\frac{1}{2}(\partial_{i}u_{j}+\partial_{j}u_{i})$ and
$M_{ij}=\frac{1}{2}(\partial_{i}b_{j}+\partial_{j}b_{i})$ are the velocity and
magnetic rates of strain respectively, then
$\tau_{ij}^{-}\tau_{ij}^{+}=S_{ij}S_{ij}-M_{ij}M_{ij}=S^{2}-M^{2}$, where
$S^{2}$ and $M^{2}$ are the respective squares of the velocity and magnetic
rate of strain tensors. These are the associated second tensor invariants,
equal to the sum of squares of the respective eigenvalues.

The second term on the right hand side simplifies further to yield
$R_{ij}^{-}R_{ij}^{+}=\frac{1}{4}(\nabla\times\mathbf{z}^{-})_{k}(\nabla\times\mathbf{z}^{+})_{k}=\frac{1}{4}(\omega_{k}-j_{k})(\omega_{k}+j_{k})=\frac{1}{4}(\omega_{k}\omega_{k}-j_{k}j_{k})$,
where $\mathbf{\omega}=\nabla\times\textbf{u}$ is the vorticity and
$\mathbf{j}=\nabla\times\textbf{b}$ is the current density. Therefore, Eqn. 6
can be re-written as

$-\partial_{ii}P=(S^{2}-M^{2})-\frac{1}{4}(\omega_{i}\omega_{i}-j_{i}j_{i}).$
(7)

Eqn. 7 is the Poisson equation for pressure, which shows that the pressure
fluctuations are related to the rate of velocity and magnetic field strains
and their rotations. In the absence of a magnetic field, Eqn. 7 reduces to
$-\partial_{ii}P=S^{2}-\frac{1}{4}\omega_{i}\omega_{i}$.

###  2.3 Poisson’s equation with shear

Using Reynolds’ decomposition, we can write a physical quantity (say $P$) as
the sum of mean ($\overline{P}$) and fluctuating component ($\tilde{p}$).
Substituting in Eqn. 4, and using $\partial_{i}z_{i}^{\pm}=0$, we get

$\partial_{ii}\overline{P}+\partial_{ii}\tilde{p}=-\partial_{i}\partial_{j}(\overline{Z}_{i}^{\mp}\overline{Z}_{j}^{\pm}+\overline{Z}_{i}^{\mp}\tilde{z}_{j}^{\pm}+\tilde{z}_{i}^{\mp}\overline{Z}_{j}^{\pm}+\tilde{z}_{i}^{\mp}\tilde{z}_{j}^{\pm}).$
(8)

Subtracting the ensemble average from itself we get

$\partial_{ii}\tilde{p}=-\partial_{i}\partial_{j}(\overline{Z}_{i}^{\mp}\tilde{z}_{j}^{\pm}+\tilde{z}_{i}^{\mp}\overline{Z}_{j}^{\pm})-\partial_{i}\partial_{j}(\tilde{z}_{i}^{\mp}\tilde{z}_{j}^{\pm}-\overline{\tilde{z}_{i}^{\mp}\tilde{z}_{j}^{\pm}}),$
(9)

which is the shear form of the Poisson equation for fluctuations in pressure.

The form of Eqn. 9 is the same as its hydrodynamic counterpart except that $u$
is replaced by $z$. Clearly, the Laplacian of pressure fluctuation is composed
of two different terms. The first term represents the distortion produced by
the shear as a result of interaction of shear with turbulence, while the
second term is the interaction among the turbulent fluctuations. Solving Eqns.
7 and 9 for pressure requires specific boundary conditions which is out of
scope of this paper. Instead, we focus ourselves on the statistics of
pressure, as follows. The characteristics of pressure fluctuation can be
studied using the pressure correlation functions $R_{P}(\mathbf{r})=\langle
P(\mathbf{x}+\mathbf{r})P(\mathbf{x})\rangle$ (Batchelor, 1951) or the second-
order pressure structure functions $D^{2}_{P}(\textbf{r})=\lvert\delta
P(\textbf{r})\rvert^{2}$ (Obukhoff & Yaglom, 1951) both of which are related
to each other (Monin & Yaglom, 1975). Here, x represents the positions space,
$\mathbf{r}$ represents the spatial lag, $\langle...\rangle$ represents the
ensemble average, and the pressure increment is defined as $\delta
P(\textbf{r})=P(\textbf{x}+\textbf{r})-P(\textbf{x})$.

The scaling of the pressure spectrum in the inertial range can be readily
identified by applying Kolmogorov similarity hypothesis to pressure increments
(Kolmogorov, 1991). The second-order structure function $D^{2}_{P}(r)$ depends
only on spatial lag $r$ and rate of dissipation of energy $\epsilon$, and can
be written as

$D^{2}_{P}(r)=\langle|P(\mathbf{x}+\mathbf{r})-P(\mathbf{x})|^{2}\rangle=C_{p}\epsilon^{4/3}r^{4/3},$
(10)

where $C_{p}$ is a constant. The form in Eqn. 10 is valid in the inertial
range. A Fourier transfrom of Eqn. 10 results in the pressure spectrum
$E_{P}(k)$ with the mathematical form

$E_{P}(k)\propto\epsilon^{4/3}k^{-7/3},$ (11)

where $k$ is the wavenumber. The second-order pressure structure function has
also been linked to the fourth-order velocity structure function (Hill &
Wilczak, 1995; Hill & Boratav, 1997; Nelkin & Chen, 1998) and pressure-
gradient velocity-velocity structure function (Hill & Boratav, 2001). However,
exploring such relationship is beyond the scope of this paper, where our
emphasis is on kinetic plasmas.

##  3 Simulation and Results

To study the statistics of pressure fluctuations in kinetic plasmas, we
analyze a $2.5$D fully kinetic particle-in-cell (PIC) simulation of turbulence
performed using the P3D code (Zeiler et al., 2002). The simulation follows the
normalization where length is normalized to the ion-inertial length
$d_{i}=c/\omega_{pi}$, where $c$ is the speed of light and $\omega_{pi}$ is
the plasma frequency for ions; time is normalized to the inverse of the
cyclotron frequency $\omega_{ci}$ and velocity is normalized to Alfvén speed
$v_{A0}=B_{0}/\sqrt{4\pi m_{i}n_{0}}$, where $m_{i}$ is the mass of ions,
$B_{0}$ and $n_{0}$ are the normalizing parameters for magnetic field and
number density, respectively. Similarly, temperature is normalized to
$m_{i}v_{A0}^{2}$ and pressure is normalized to $n_{0}m_{i}v_{A0}^{2}$.

The simulation is a periodic square domain of length $L=149.6d_{i}$ with grid
points of $4096^{2}$, and $3200$ particles initially in each grid. The mass
ratio of ions to electrons is set to $25$ with an initial background density
of $1$, uniform initial temperature of $T_{e}=T_{i}=0.3$ with a total plasma
beta $\beta=1.2$, and an initial out-of-plane magnetic field of strength $1$.
The system is initially populated with Fourier modes for magnetic field and
velocity (both ions and electrons) within $2k_{0}\leq|k|\leq 4k_{0}$, where
$k_{0}=2\pi/L$ and allowed to evolve for about $370\,t\omega_{ci}$ without
external forcing (see Parashar et al. (2018), Adhikari et al. (2021) for
details).

###  3.1 Pressure balances

Fig. 1 shows the time evolution of the ensemble averaged (i.e., volume
integrated) change in magnetic $\Delta P_{mag}$ and total thermal $\Delta
P_{th}$ pressure. The change in the pressure is calculated with respect to the
initial value at $t=0$. Initially, the magnetic pressure increases rapidly for
a small period of time and then falls, while the thermal pressure continues to
increase gradually throughout the simulation. Note, the change in the total
pressure stays roughly constant once it reaches the maximum value; as the
thermal pressure compensates for the fall in the magnetic pressure. The change
in the magnetic and thermal pressure intersects at $t\omega_{ci}=116.5$ when
the mean square current is maximum in the system (not shown). This time, when
MHD dissipation typically reaches a maximum, is often regarded as the
turbulence having reached a fully developed state. One finds that in this
state, the volume-integrated magnetic pressure and the volume-integrated
thermal pressure are anticorrelated, producing a near-constant volume
integrated total pressure. Since we are interested in the turbulence
properties of pressure fluctuations, we focus our analysis around this time.

Figure 1:  Time evolution of the change in thermal $P_{th}$, magnetic
$P_{mag}$ and total pressure $P_{tot}$ calculated with respect to their
initial values. $P_{mag}$ and $P_{th}$ have opposite behavior once the
magnetic pressure reaches the maximum value. The red star represents the time
of analysis. Figure 2: Two dimensional image at $t\omega_{ci}=116.5$ of the
magnetic $P_{mag}$, thermal $P_{th}$, and total $P_{tot}$ pressure
fluctuations normalized to their respective mean values. Figure 3: Joint
probability distribution of the thermal $P_{th}$ and magnetic $P_{mag}$
pressure at $t\omega_{ci}=116.5$. A dashed line of slope $-1$ is drawn for
reference and $r_{P}$ represents the Pearson correlation coefficient.

To understand how the global pressure balance impacts the local pressures, in
Fig. 2 we show the two-dimensional view of the thermal, magnetic, and total
pressure fluctuations at $t\omega_{ci}=116.5$, each normalized to their
respective mean value. It is observed that most of the regions near the
current sheets (see Fig.10) have low magnetic pressure and high thermal
pressure indicating a negative correlations between the two. As a result the
total pressure in these regions stay relatively uniform as seen on the right
panel of Fig. 2. On the other hand, the magnetic islands have similar
magnitude of $P_{mag}$ and $P_{th}$. This is certainly due to the magnetic
compression and heating within the islands. The total pressure balance is
further quantified by calculating the variance of the pressure fluctuations
plotted in Fig. 2, and is shown in Table. 1. The variance of the magnetic and
thermal pressure fluctuations have similar magnitude, while the variance of
the total pressure fluctuation is much smaller ($\sim 40\%$ of $P_{th}$ or
$P_{mag}$).

Next, we plot the joint probability distribution function (PDF) of the thermal
and magnetic pressure in Fig. 3, where a line of slope $-1$ is drawn for
reference. The joint PDF clearly follows the reference line suggesting a
negative correlation between the two pressures. In addition, we compute the
Pearson correlation coefficient $r_{P}$ (Cohen et al., 2009) between $P_{th}$
and $P_{mag}$ to be $-0.257$ supporting the trend seen in Fig. 3. This
anticorrelation quantifies the tendency towards statistical (in contrast to
pointwise) pressure balance.

Table 1: Variance about the spatial mean of the magnetic, thermal and total pressure as plotted in Fig. 2. Variance | $P_{th}$ | $P_{mag}$ | $P_{tot}$
---|---|---|---
$\sigma^{2}$ | $0.0107$ | $0.0101$ | $0.0039$
Figure 4: (Top) Pressure spectra for magnetic pressure (dashed line), thermal
pressure (dash-dot) and total pressure (solid green line) as a function of
wavenumber at the time of analysis. Reference lines with slopes $-5/3$, $-7/3$
and $-3$ are provided. (Bottom) Magnetic energy spectra at the same time.
Lines of slope $-5/3$, $-8/3$, and $-11/3$ are drawn for reference. Both total
spectrum and component spectra are shown (See legend).

We now discuss the pressure spectra. In the top panel of Fig. 4, we plot the
magnetic, thermal, and total pressure spectra. The dotted vertical lines
represent the wavenumber corresponding to the ion inertial length $d_{i}^{-1}$
, electron inertial length $d_{e}^{-1}$ and the Debye length
$\lambda_{D}^{-1}$. Additional lines of slope $-5/3$, $-7/3$, and $-3$ are
drawn for reference. The thermal and magnetic pressure spectra exhibit very
similar behavior except at the higher wavenumbers. In the inertial range
$kd_{i}\lesssim 1$, both these spectra exhibit similar spectral slopes close
to $-5/3$, while at large wavenumbers $kd_{i}>5$, the thermal pressure has a
plateau at approximately $10^{-7}$ while the magnetic pressure spectrum falls
down below $10^{-8}$.

Once again, the total pressure displays different behavior compared to its
constituent parts. The total pressure spectrum, in the inertial range displays
a spectral slope of $-7/3$ consistent with Eqn. 11. This is due to the
significant cross-spectrum correlation between the thermal and the magnetic
pressure, a direct consequence of the statistical pressure balance described
above.

At higher wavenumber, the total pressure spectrum steepens beyond the inertial
range $2\leq kd_{i}\leq 9$ with slope $-3$ before the hump on the thermal
pressure spectrum dominates the contribution.

The bottom panel of Fig. 4 shows the magnetic energy spectra along with its
components. In the inertial range ($kd_{i}\lesssim 1$) the magnetic spectra
approximate the familiar Kolmogorov $-5/3$ slope as expected in strong
turbulence. Near $kd_{i}=1$ the magnetic spectra steepen, as often reported
for plasma (Leamon et al., 1998), with a slope that varies in the range of
about $-7/3$ to $-8/3$ (Smith et al., 2006). In this range kinetic dissipation
and dispersion effects become important. Various steeper forms have been
proposed for the spectrum at still higher wavenumbers $kd_{e}>1$; an $11/3$
spectrum is shown here for reference. The spectrum of the out of plane
component $B_{z}$ is much flatter than the other two magnetic components, and
in the inertial range appears at a much lower amplitude, down by as much as
two orders of magnitude at the long wavelength end of the inertial range.
However at the intermediate kinetic range, the $B_{z}$ spectrum reaches
equipartition with the spectra of the other two components, and further along,
at sub electron scales, this component becomes dominant over the other two by
a modest factor.

We note in passing that the theoretical development in (Montgomery et al.,
1987) began with assumption that the magnetic spectrum exhibits a $-5/3$
spectrum and proceeded to compute the density and pressure spectra as a linear
response to the magnetic field. The conclusion was that in MHD the thermal
pressure (and density) spectra would take on a $-5/3$ power law as seen here.
The anticorrelation of thermal and magnetic pressure at inertial range scales
was also concluded. The statistical characterization of total pressure was not
included. The present findings are complementary to this important antecedent.

Figure 5: Probability density functions of pressure increments $\delta P$, for
magnetic pressure (blue circles), thermal pressure (red diamonds) and total
pressure (green boxes), each normalized by its standard deviation, at
$r\approx 1d_{e}=0.2d_{i}$ (left), and $r\approx 1d_{i}$ (right). The PDFs are
computed at $t=116.5\omega_{ci}^{-1}$. The black solid line is a normal
distribution for reference (see text for details).

Next, we show the probability density function (PDF) of the pressure
increments $(\delta P)$ defined in section 2.3 as

$\delta P\mathbf{(r})=P(\mathbf{x}+\mathbf{r})-P(\mathbf{x}),$ (12)

where $\mathbf{r}$ is the spatial lag. For simplicity, we compute the PDF
using lags along the $x$ and $y$ directions, and average the results. One can
also perform these calculations for lags along different directions and take
an angular average (Wang et al., 2022). The PDF’s are shown in Fig. 5 for two
different lags $r=1d_{e}\approx 0.2d_{i}$ (left panel) and $r\approx 1d_{i}$
(right panel).

At $r\approx d_{i}$, the PDF of both $\delta P_{mag}$ (blue circles) and
$\delta P_{th}$ (red diamonds) follow each other closely in contrast to the
$r\approx d_{e}$ case, however, displaying a departure from gaussianity in
both cases. The total pressure PDF (green box) for $r\approx d_{e}$ is very
close to a normal distribution (solid black line), drawn for reference, but is
distorted for $r\approx d_{i}$ implying better pressure balance at smaller
lags.

Figure 6: Structures (in black) contributing to regions I, II and III of the
probability density function (PDF) of the magnetic (left column), thermal
(middle column) and total pressure (right column) increments shown in Fig. 5
for a lag of $r\approx d_{e}$, averaged over $x$ and $y$ directions. Each row
corresponds to the region shown on the right side of the figure.

To examine the spatial distributions that originate from different sections of
these PDF’s, we divide the probability distribution function into regions I,
II and III (as in Greco et al. (2009)) and mask out all but one region in a
given plot. The results for small lag increments, at $r\approx d_{e}$, are
shown in Fig. 6. Here we have used a binary color representation where black
represents the structures contributing to a particular region, and white
represents the region with no contribution. It is immediately evident that the
total pressure increments at this scale are much more homogeneous than the
separate contributions from magnetic or thermal pressures. This is the case
for the region I core of the distributions, for the fine scale super-gaussian
tails representing region III, and for the intermediate subgaussian regions
II. This signature of pressure balance is particularly clear for the
supergaussian tail. However none of these regions is completely uniform. In
fact the spatial distribution of magnetic islands and their boundaries, i.e.,
what we sometimes called cellularization of the plasma (Servidio et al.,
2008), is evident in all three ranges of the PDF, as is evident in the panels
of Fig. 6. Quantitative diagnostics of this spatial structure leads us to a
discussion of intermittency.

Theories and classifications of intermittency provide a physical basis for
quantifying deviations from gaussianity. To understand the intermittency of
pressure fluctuations, we calculate the scale-dependent kurtosis and higher-
order structure function for pressure increments in the following sub-section.

###  3.2 Higher Order Statistics and Intermittency

####  3.2.1 Scale dependent kurtosis

The scale-dependent kurtosis $\kappa(\mathbf{r})$ for the pressure
fluctuations $\delta P(\mathbf{r})$ as a function of lag $\mathbf{r}$ is
defined as

$\kappa(\mathbf{r})=\frac{\langle\lvert
P(\mathbf{x}+\mathbf{r})-P(\mathbf{x})\rvert^{4}\rangle}{\langle\lvert
P(\mathbf{x}+\mathbf{r})-P(\mathbf{x})\rvert^{2}\rangle^{2}}.$ (13)

$\kappa(\mathbf{r})$ takes on a value of three for a scalar field having a
gaussian distribution. For a scalar quantity (or its increments) less randomly
distributed and more clustered in space, $\kappa(\mathbf{r})$ attains larger
values, and can be heuristically interpreted as the reciprocal of the spatial
filling factor. This property makes $\kappa(\mathbf{r})$ a good preliminary
indicator of the presence of spatial intermittency.

In Fig. 7 we plot $\kappa(r)$ for magnetic, thermal, and total pressure. It is
observed that $\kappa(r)$ is non-gaussian for all pressure forms at the
smaller lags $r\leq 1di$. This agrees with the PDF seen in Fig. 5. The scale-
dependent kurtosis for the magnetic field pressure peaks at $\approx
0.1d_{i}$, that for thermal pressure peaks at $\approx 0.3d_{i}$, while
$\kappa(r)$ of total pressure peaks at $\approx 0.5d_{i}$, approaching a value
closer to $3$ for larger lags. This implies that the fluctuations at smaller
lags are intermittent, and therefore the original Kolmogorov (1941) similarity
hypothesis is not valid and needs some discussion or clarification.

Figure 7: Scale dependent kurtosis $\kappa(\mathbf{r})$ for magnetic (dash-
dot), thermal (dashed) and total pressure (solid) as a function of lag. The
dotted line is a reference for the kurtosis of a gaussian distribution.

####  3.2.2 Structure functions

In this subsection, we discuss the higher order structure functions for
pressure fluctuations as measured by the increments $\delta P$, and compare
them with the higher order velocity and magnetic structure functions. The
structure function of order $n$ is defined as

$D_{P}^{n}(r)=\langle\lvert\delta P(r)\rvert^{n}\rangle.$ (14)

For Kolmogorov 1941 scaling ($K41$ from here on), which ignores intermittency,
(see section 2.3) the pressure fluctuation is expected to scale as $\delta
P\sim(\delta v)^{2}$, and therefore the $n^{th}$ order pressure structure
function $D_{P}^{n}(r)\sim(\delta v)^{2n}\sim r^{2n/3}$. However, this
idealization is not expected to be valid for higher-order structure functions
because of the non-gaussian behavior as seen in Fig. 5. This section explores
the lag dependence of the higher-order pressure structure functions in the
kinetic plasmas of interest.

The multi-fractal scaling of pressure fluctuations has not been explored as
much as the velocity and magnetic field. Some attempts have been made to
include intermittency corrections in the pressure scaling (Donzis et al.,
2012), but not for kinetic plasma turbulence as far as we are aware.

Figure 8:  Sequence of higher order structure functions for pressure. For
uniformity, we assume that the region between the dashed lines at
$d_{e}=0.2d_{i}$ and $d_{i}$ denotes the kinetic range (brown rectangle), and
the region between $4-15d_{i}$ represents the inertial range (olive rectangle)
(see Adhikari et al. (2021)).

For further analysis, first we assume the pressure structure function to take
the form

$D_{P}^{n}(r)\sim r^{\zeta_{n}},$ (15)

where $\zeta_{n}$ includes the intermittency correction to the $K41$ scaling
(Kolmogorov, 1962; Obukhov, 1962). We then compute the scaling exponents for
our system and compare that with a few intermittency models that are
formulated to describe the scaling of velocity structure functions.

One of the first descriptions of intermittency comes from Oboukhov’s
suggestion that motivated the Kolmogorov refined similarity hypothesis
(Oboukhov, 1962; Kolmogorov, 1962). Specifically Oboukhov reasoned that patchy
dissipation might be accounted for if the velocity fluctuation (i.e.,
increment) depends on the local value of dissipation and exhibits a log-normal
distribution. According to this model ($K62$ from here on), the scaling
exponent takes the form

$\zeta_{n}^{K62}=\frac{n}{3}+\frac{\mu}{18}(3n-n^{2}),$ (16)

where $\mu$ is the intermittency exponent, whose values typically ranges
between $0.2-0.3$ (Van Atta & Antonia, 1980; Sreenivasan & Kailasnath, 1993;
Davidson, 2015).

Another successful model used to describe intermittency is the She and Leveque
model (She & Leveque, 1994) which assumes a log-Poisson model distribution. In
this model, the exponent has the generalized form

$\zeta_{n}=\frac{n}{3}(1-\alpha)+C_{0}(1-\beta^{n/3})],$ (17)

where $C_{0}=\frac{\alpha}{1-\beta}$ is related to the dimension D of the
dissipative eddies through the relation $C_{0}=3-D$ (Biskamp, 2003) and
therefore called the co-dimension. The parameter $\alpha$ is the scaling
exponent, while $\beta$ is a characteristic of the intermittency. For
hydrodynamic turbulence, because of the quasi-$1$D structure of thin vortex
filaments, the co-dimension $C_{0}=2$. She and Leveque (She & Leveque, 1994)
further assumed that with Kolmogorov scaling, $\alpha=2/3$ and $\beta=2/3$,
which yields

$\zeta_{n}^{SL}=\frac{n}{9}+2\left[1-\left(\frac{2}{3}\right)^{n/3}\right].$
(18)

However, for magnetohydrodynamic turbulence the smallest dissipative
structures may be taken to be $2$D as we observe from small scale current
sheets. For that case, $C_{0}=1$. With the same $x=2/3$, we now have
$\beta=1/3$ that yields the scaling exponents (Politano & Pouquet, 1995;
Müller & Biskamp, 2000) of the form

$\zeta_{n}^{MHD}=\frac{n}{9}+1-\left(\frac{1}{3}\right)^{n/3}.$ (19)

Next, in Fig. 8 we show the higher-order structure functions $D_{P}^{n}$
versus lag $r$ up to the eighth order. Two vertical lines on the left denote
the kinetic range (between $d_{e}=0.2d_{i}$ and $1d_{i}$), while two vertical
lines on the right denote the inertial range (between $4d_{i}$ and $15d_{i}$)
(Adhikari et al., 2021). The magnitude of the higher-order structure functions
is smaller than the lower-order structure functions. The structure functions
are observed to be steeper in the kinetic range compared to the inertial
range. The steepness of the kinetic range increases sharply with higher-order
structure function compared to the inertial range because the fluctuations at
those lag scales are much smaller than those in the inertial range. Beyond the
inertial range, the structure functions display a similar slope. Next, we
calculate the power exponent of these structure functions. In Fig. 9 we plot
the scaling exponent of these structure functions as a function of the order
$n$, across the inertial (top) and kinetic (bottom) range and compare it with
those from the magnetic and velocity structure functions. The exponents are
obtained using power-law fits to the structure function curve at both the
inertial and kinetic ranges. It is instructive to compare the present results
with similar findings presented for kinetic plasma turbulence by (Wan et al.,
2016) and for the density in compressible MHD turbulence by (Yang et al.,
2017b).

The scaling exponents for the structure function follow the order velocity $<$
magnetic field $<$ pressure, with the ratio between pressure and velocity
structure function $\sim 2$. The value of $\zeta_{n}$ for the
velocity/magnetic and pressure structure functions decrease more and more from
the conventional (K$41$) value of $n/3$ and $2n/3$ respectively as $n$
increases. Interestingly, the exponent saturates to a value close to $1$ for
the velocity structure function.

Figure 9: Power index for different orders $n$ of the magnetic (blue squares),
velocity (green circles), and pressure (red diamond) structure functions in
the inertial (top), and kinetic (bottom) range as defined by the shaded
rectangles in Fig. 8. The scaling exponents predicted by different models in
the inertial range are drawn for reference (see text for details). In the
kinetic range, the scaling exponents follow a linear trend with a slope of
$0.927$ (velocity), $0.801$ (magnetic field), and $0.779$ (pressure).

We also plot the K$41$, SL, and K$62$ prediction for comparison. Clearly, none
of these models fully describe the exponents obtained from our simulations.
However, interestingly, the exponent for the magnetic structure functions is
in good agreement with the generalized She Leveque exponent (Eqn. 17) with
scaling exponent $x=3/4$, and $\beta=1/4$, such that the co-dimension
$C_{0}=1$. This curve is denoted by MSL (modified SL) in Fig. 9.

In the kinetic range (bottom panel Fig. 9), the exponents for all the
structure functions exhibits a linear relation with the order of the structure
function, contrary to the inertial range behavior. The values of the exponents
for lower-order $(n\lesssim 6)$ magnetic and velocity structure functions are
similar. However, for $n>6$ the difference in $\zeta_{n}$ starts to increase
and diverges with increasing $n$ (not shown). On the contrary, the power
exponent for the pressure structure functions is smaller in the kinetic range
suggesting that the previous description in the inertial range is not valid in
the kinetic range. More studies emphasizing the kinetic regime of the
structure functions will certainly shed some light on our understanding of the
structure functions.

Recent works on intermittency in the magnetic field provide guidance in this
issue, while it remains clear that some subtle distinctions may be needed. For
example, Leonardis et al. (2013) found multifractal scaling of the magnetic
field at subproton scales in simulations of kinetic reconnection, but also
found monofractal gaussian scaling at scales less than $d_{i}$ (Leonardis et
al., 2016). Other kinetic simulations of collisonless plasma turbulence
supported the finding that turbulence at proton scales and smaller is self-
similar, i.e., monofractal (Wan et al., 2016). Laboratory findings (e.g.,
(Schaffner & Brown, 2015)) and Parker Solar Probe data (Chhiber et al., 2021)
also find monofractal behavior at sub-proton scales. Here we have found very
similar results for intermittency in the total pressure. First we have found
that the pressure fluctuations are indeed intermittent. Second, the finding of
of multifractal scaling at scales larger than $d_{i}$, and monofractal scaling
at scales less than $d_{i}$ is in accordance with characterization of magnetic
field intermittency in a number of prior studies.

##  4 Conclusions

Pressure balance has long been an important element in space plasma physics
(Burlaga et al., 1990). The traditional viewpoint adopts explicitly an
assumption that fluctuations are small, leading to a menu of linear wave
solutions and small amplitude pressure fluctuations that emerge as special
static solutions of the MHD equations (Cowling, 1976). Here we have employed
particle-in-cell simulations of kinetic turbulence starting from conditions
that give rise to strong MHD scale fluctuations and subsequently the formation
of strong spatial intermittency, coherent structures, and nonuniform
dissipation.

In this paper, we emphasize the properties of pressure fluctuations that are
of inherently nonlinear and stochastic nature, and therefore the
characteristics we describe are attributable to turbulence and not a
collection of static or wavelike structures.

The key conclusions of this paper are listed below:

  1. 1.

Poisson’s equation for pressure (Eqn. 7) in terms of Elsässer’s variable was
presented. This equation explains the quasi-pressure balance obtained in
turbulent systems. Similarly, Poisson’s equation for the pressure fluctuation
is formulated in terms of Elsässer shear (Eqn. 9). The form of the equation is
the same as it would be for velocity replaced by Elsässer variables.

  2. 2.

In a PIC simulation of decaying turbulence, as turbulence is fully developed,
the volume averaged change in the thermal pressure behaves opposite to that of
the magnetic pressure keeping the total pressure nearly constant as seen in
Fig. 1. A 2D view of the thermal (ion plus electron) and magnetic pressure
indicates the suggested signatures of anticorrelation between the magnetic and
the thermal pressure locally, and results in the global anticorrelation (Fig.
2).

  3. 3.

The omnidirectional pressure spectrum has a characteristic slope of $-7/3$
(Fig. 4) in the inertial range while the magnetic and thermal pressure each
have a slope of $-5/3$, a less steep power spectrum than the total pressure
spectrum. The significant cross contribution between the two pressure terms
(magnetic and thermal) is a result of the anti-correlation between the two.

  4. 4.

The intermittent behavior of pressure fluctuations is discussed and contrasted
to the intermittency observed in the velocity and magnetic fields. The ratio
of the scaling exponents for the higher-order pressure structure function to
that of velocity structure function in the inertial range is $\sim 2$. None of
the previously described intermittency models fit the scaling exponent
obtained from our simulations; resolution of this issue requires detailed
study that is outside the scope of the present paper. The exponents for
velocity structure functions are much lower than the predicted value for
higher orders $n$ and saturate close to $1$. The scaling exponent for magnetic
structure functions does agree with the MHD description of the SL model but
with a different scaling exponent. Further, both monofractal and multifractal
scaling exponents are observed in the pressure, suggesting that traces of the
intermittent behavior are also observed in pressure fluctuations.

This study sheds light on the contribution of thermal and magnetic pressure to
the total pressure using pressure spectra and higher-order statistics. One can
use these results to characterize the fluctuations of mass density in many
astrophysical systems, exploring the contributions coming from the thermal and
magnetic pressure in the system (Lau et al., 2013). Therefore, studies of
pressure fluctuations in astrophysical systems can only lead to a better
understanding of the system.

Finally it is not yet clear if the findings presented here will extrapolate to
all the cases of two or three dimensions, including the expected anisotropy of
magnetic fields (Oughton et al., 1994). Given the great variation in physical
processes that occur in widely different parameter regimes of three
dimensional plasma turbulence, it is also not clear if the conclusions present
here will prove to be universal across all attainable parameters. Future
studies on these extensions are therefore warranted.

## Acknowledgements

S. A. acknowledges helpful discussions with Francesco Pecora. All the authors
would like to acknowledge the high-performance computing support from Cheyenne
(Computational and Information Systems Laboratory, 2017) provided by NCAR’s
Computational and Information Systems Laboratory, sponsored by the NSF. We
also thank NERSC resources, a U.S. DOE Office of Science User Facility
operated under Contract No. DE-AC02-05CH11231. S. A, and M. A. S. acknowledge
support from NASA LWS 80NSSC20K0198. P. A. C gratefully acknowledges the
hospitality from the University of Delaware and the Bartol Research Institute
during his sabbatical visit, and acknowledges support from NASA Grant No.
80NSSC19M0146, and the Bartol Research Institute. S. A. and P. A. C. are also
supported by DOE grant DE-SC0020294. W. H. M is supported by NSF DOE grant AGS
2108834 at the University of Delaware, the IMAP project (Princeton subcontract
SUB0000317)) and the NASA LWS program FST grant to New Mexico Consortium
(subcontract 655-001 to Delaware).

## Data Availability

The data that support the findings of this study are openly available in
Zenodo (Adhikari et al., 2023).

##  Appendix A System Overview

Fig. 10 shows the magnitude of the out-of-plane current ($j_{z}$) in the
system at the time of analysis.

Figure 10:  Magnitude of the out-of-plane current density $j_{z}$ at the time
of analysis $(t\omega_{ci}=116.5)$. The current sheets are the regions with
larger magnitudes of $j_{z}$.

## References

  * Adhikari et al. (2021) Adhikari S., Parashar T., Shay M., Matthaeus W., Pyakurel P., Fordin S., Stawarz J., Eastwood J., 2021, Physical Review E, 104, 065206
  * Adhikari et al. (2023) Adhikari S., Matthaeus W., Parashar T., Shay M., Cassak P., 2023, [Dataset], Zenodo
  * Angelinelli et al. (2020) Angelinelli M., Vazza F., Giocoli C., Ettori S., Jones T., Brunetti G., Brüggen M., Eckert D., 2020, Monthly Notices of the Royal Astronomical Society, 495, 864
  * Armstrong et al. (1981) Armstrong J., Cordes J., Rickett B., 1981, Nature, 291, 561
  * Balbus & Hawley (1998) Balbus S. A., Hawley J. F., 1998, Reviews of modern physics, 70, 1
  * Batchelor (1951) Batchelor G., 1951, in Mathematical Proceedings of the Cambridge Philosophical Society. pp 359–374
  * Bayly et al. (1992) Bayly B. J., Levermore C. D., Passot T., 1992, Physics of Fluids A: Fluid Dynamics, 4, 945
  * Biskamp (2003) Biskamp D., 2003, Magnetohydrodynamic turbulence. Cambridge University Press
  * Bruno & Carbone (2013) Bruno R., Carbone V., 2013, Living Reviews in Solar Physics, 10, 1
  * Burlaga et al. (1990) Burlaga L., Scudder J., Klein L., Isenberg P., 1990, Journal of Geophysical Research: Space Physics, 95, 2229
  * Carlberg et al. (1997) Carlberg R., et al., 1997, The Astrophysical Journal, 485, L13
  * Chhiber et al. (2021) Chhiber R., Matthaeus W. H., Bowen T. A., Bale S. D., 2021, The Astrophysical Journal Letters, 911, L7
  * Churazov et al. (2012) Churazov E., et al., 2012, Monthly Notices of the Royal Astronomical Society, 421, 1123
  * Cohen et al. (2009) Cohen I., Huang Y., Chen J., Benesty J., Benesty J., Chen J., Huang Y., Cohen I., 2009, Noise reduction in speech processing, pp 1–4
  * Coleman Jr (1968) Coleman Jr P. J., 1968, The Astrophysical Journal, 153, 371
  * Computational and Information Systems Laboratory (2017) Computational and Information Systems Laboratory 2017, Cheyenne: HPE/SGI ICE XA System (University Community Computing), doi:10.5065/D6RX99HX
  * Cowling (1976) Cowling T. G., 1976, Bristol, Adam Hilger, Ltd.
  * Davidson (2015) Davidson P. A., 2015, Turbulence: an introduction for scientists and engineers. Oxford university press
  * Donzis et al. (2012) Donzis D., Sreenivasan K., Yeung P., 2012, Physica D: Nonlinear Phenomena, 241, 164
  * George et al. (1984) George W. K., Beuther P. D., Arndt R. E., 1984, Journal of Fluid Mechanics, 148, 155
  * Goldstein et al. (1995) Goldstein M. L., Roberts D. A., Matthaeus W., 1995, Annual review of astronomy and astrophysics, 33, 283
  * Greco et al. (2009) Greco A., Matthaeus W. H., Servidio S., Chuychai P., Dmitruk P., 2009, apjl, 691, L111
  * Hellinger et al. (2013) Hellinger P., Trávníček P. M., Štverák Š., Matteini L., Velli M., 2013, Journal of Geophysical Research: Space Physics, 118, 1351
  * Higdon (1984) Higdon J. C., 1984, Astrophys. J., 285, 109
  * Hill & Boratav (1997) Hill R. J., Boratav O. N., 1997, Physical Review E, 56, R2363
  * Hill & Boratav (2001) Hill R. J., Boratav O. N., 2001, Physics of Fluids, 13, 276
  * Hill & Wilczak (1995) Hill R. J., Wilczak J. M., 1995, Journal of Fluid Mechanics, 296, 247
  * Khatri & Gaspari (2016) Khatri R., Gaspari M., 2016, Monthly Notices of the Royal Astronomical Society, 463, 655
  * Kim & Antonia (1993) Kim J., Antonia R., 1993, Journal of Fluid Mechanics, 251, 219
  * Klainerman & Majda (1981) Klainerman S., Majda A., 1981, Commun. Pure Appl. Math., 34, 481
  * Kolmogorov (1962) Kolmogorov A. N., 1962, Journal of Fluid Mechanics, 13, 82
  * Kolmogorov (1991) Kolmogorov A. N., 1991, Proceedings of the Royal Society of London. Series A: Mathematical and Physical Sciences, 434, 9
  * Kravtsov & Borgani (2012) Kravtsov A. V., Borgani S., 2012, Annual Review of Astronomy and Astrophysics, 50, 353
  * Lau et al. (2013) Lau E. T., Nagai D., Nelson K., 2013, The Astrophysical Journal, 777, 151
  * Leamon et al. (1998) Leamon R. J., Smith C. W., Ness N. F., Matthaeus W. H., Wong H. K., 1998, J. Geophys. Res., 103, 4775
  * Leonardis et al. (2013) Leonardis E., Chapman S. C., Daughton W., Roytershteyn V., Karimabadi H., 2013, Phys. Rev. Lett., 110, 205002
  * Leonardis et al. (2016) Leonardis E., Sorriso-Valvo L., Valentini F., Servidio S., Carbone F., Veltri P., 2016, Physics of Plasmas, 23, 022307
  * Markevitch et al. (2004) Markevitch M., Gonzalez A., Clowe D., Vikhlinin A., Forman W., Jones C., Murray S., Tucker W., 2004, The Astrophysical Journal, 606, 819
  * Matthaeus & Brown (1988a) Matthaeus W. H., Brown M. R., 1988a, Phys. Fluids, 31, 3634
  * Matthaeus & Brown (1988b) Matthaeus W. H., Brown M. R., 1988b, The Physics of fluids, 31, 3634
  * Matthaeus & Goldstein (1982) Matthaeus W. H., Goldstein M. L., 1982, Journal of Geophysical Research: Space Physics, 87, 6011
  * Matthaeus et al. (1991) Matthaeus W. H., Klein L. W., Ghosh S., Brown M. R., 1991, Journal of Geophysical Research: Space Physics, 96, 5421
  * Monin & Yaglom (1975) Monin A., Yaglom A., 1975, Statistical fluid mechanics: Mechanics of Turbulence. "" Vol. 2, MIT Press, Cambridge, MA
  * Montgomery et al. (1987) Montgomery D., Brown M. R., Matthaeus W., 1987, Journal of Geophysical Research: Space Physics, 92, 282
  * Müller & Biskamp (2000) Müller W.-C., Biskamp D., 2000, Physical Review Letters, 84, 475
  * Nelkin & Chen (1998) Nelkin M., Chen S., 1998, Physics of Fluids, 10, 2119
  * Oboukhov (1962) Oboukhov A. M., 1962, Journal of Fluid Mechanics, 13, 77
  * Obukhoff & Yaglom (1951) Obukhoff A., Yaglom A., 1951, NACA Tech. Mem, 1350
  * Obukhov (1962) Obukhov A. M., 1962, J. Geophys. Res., 67, 3011
  * Oughton et al. (1994) Oughton S., Priest E. R., Matthaeus W. H., 1994, J. Fluid Mech., 280, 95
  * Parashar et al. (2018) Parashar T. N., Matthaeus W. H., Shay M. A., 2018, The Astrophysical Journal Letters, 864, L21
  * Pessah (2010) Pessah M. E., 2010, The Astrophysical Journal, 716, 1012
  * Planelles et al. (2015) Planelles S., Schleicher D., Bykov A., 2015, Space Science Reviews, 188, 93
  * Politano & Pouquet (1995) Politano H., Pouquet A., 1995, Physical Review E, 52, 636
  * Pratt et al. (2019) Pratt G., Arnaud M., Biviano A., Eckert D., Ettori S., Nagai D., Okabe N., Reiprich T., 2019, Space Science Reviews, 215, 1
  * Pullin & Rogallo (1994) Pullin D., Rogallo R., 1994, Stanford Univ., Studying Turbulence Using Numerical Simulation Databases. 5: Proceedings of the 1994 Summer Program
  * Reisenfeld et al. (1999) Reisenfeld D. B., McComas D., Steinberg J., 1999, Geophysical research letters, 26, 1805
  * Ruffolo et al. (2021) Ruffolo D., Ngampoopun N., Bhora Y. R., Thepthong P., Pongkitiwanichakul P., Matthaeus W. H., Chhiber R., 2021, The Astrophysical Journal, 923, 158
  * Sahraoui et al. (2020) Sahraoui F., Hadid L., Huang S., 2020, Reviews of Modern Plasma Physics, 4, 1
  * Schaffner & Brown (2015) Schaffner D., Brown M. R., 2015, The Astrophysical Journal, 811, 61
  * Schuecker et al. (2004) Schuecker P., Finoguenov A., Miniati F., Böhringer H., Briel U., 2004, Astronomy & Astrophysics, 426, 387
  * Schumann & Patterson (1978) Schumann U., Patterson G., 1978, Journal of Fluid Mechanics, 88, 685
  * Servidio et al. (2008) Servidio S., Matthaeus W., Dmitruk P., 2008, Physical review letters, 100, 095005
  * She & Leveque (1994) She Z.-S., Leveque E., 1994, Physical review letters, 72, 336
  * Simionescu et al. (2019) Simionescu A., et al., 2019, Space Science Reviews, 215, 1
  * Smith et al. (2006) Smith C. W., Hamilton K., Vasquez B. J., Leamon R. J., 2006, Astrophys. J., 645, L85
  * Squire et al. (2023) Squire J., Kunz M. W., Arzamasskiy L., Johnston Z., Quataert E., Schekochihin A. A., 2023, arXiv preprint arXiv:2303.00468
  * Sreenivasan & Kailasnath (1993) Sreenivasan K., Kailasnath P., 1993, Physics of Fluids A: Fluid Dynamics, 5, 512
  * Tsuji & Ishihara (2003) Tsuji Y., Ishihara T., 2003, Physical Review E, 68, 026309
  * Uberoi (1953) Uberoi M. S., 1953, Journal of the Aeronautical Sciences, 20, 197
  * Van Atta & Antonia (1980) Van Atta C., Antonia R., 1980, The Physics of Fluids, 23, 252
  * Vellante & Lazarus (1987) Vellante M., Lazarus A., 1987, Journal of Geophysical Research: Space Physics, 92, 9893
  * Voigt & Fabian (2006) Voigt L., Fabian A., 2006, Monthly Notices of the Royal Astronomical Society, 368, 518
  * Wan et al. (2016) Wan M., Matthaeus W. H., Roytershteyn V., Parashar T. N., Wu P., Karimabadi H., 2016, Physics of Plasmas, 23, 042307
  * Wang & Richardson (2001) Wang C., Richardson J., 2001, Journal of Geophysical Research: Space Physics, 106, 29401
  * Wang et al. (2022) Wang Y., et al., 2022, The Astrophysical Journal, 937, 76
  * Yang et al. (2017a) Yang Y., et al., 2017a, Physics of Plasmas, 24, 072306
  * Yang et al. (2017b) Yang Y., Matthaeus W. H., Shi Y., Wan M., Chen S., 2017b, Physics of Fluids, 29, 035105
  * Zank & Matthaeus (1990) Zank G. P., Matthaeus W. H., 1990, Phys. Rev. Lett., 64, 1243
  * Zank & Matthaeus (1992) Zank G. P., Matthaeus W. H., 1992, J. Geophys. Res., 97, 17 189
  * Zank & Matthaeus (1993) Zank G. P., Matthaeus W. H., 1993, Phys. Fluids A, 5, 257
  * Zeiler et al. (2002) Zeiler A., Biskamp D., Drake J., Rogers B., Shay M., Scholer M., 2002, Journal of Geophysical Research: Space Physics, 107, SMP
  * Zhao et al. (2016) Zhao S., Cheng E., Qiu X., Burnett I., Liu J. C.-c., 2016, The Journal of the Acoustical Society of America, 140, 4178
  * Zhuravleva et al. (2014) Zhuravleva I., et al., 2014, Nature, 515, 85
  * Zhuravleva et al. (2019) Zhuravleva I., Churazov E., Schekochihin A., Allen S., Vikhlinin A., Werner N., 2019, Nature Astronomy, 3, 832