# Multihyperuniform Long-Range Order in Medium-Entropy Alloys

Duyu Chen<EMAIL_ADDRESS>Materials Research Laboratory, University
of California, Santa Barbara, California 93106, United States  Xinyu Jiang
Mechanical and Aerospace Engineering, Arizona State University, Tempe, AZ
85287  Duo Wang  Materials Science and Engineering, Arizona State University,
Tempe, AZ 85287  Houlong Zhuang<EMAIL_ADDRESS>Mechanical and Aerospace
Engineering, Arizona State University, Tempe, AZ 85287  Yang Jiao
<EMAIL_ADDRESS>Materials Science and Engineering, Arizona State
University, Tempe, AZ 85287  Department of Physics, Arizona State University,
Tempe, AZ 85287

###### Abstract

We provide strong numerical evidence for a hidden multihyperuniform long-range
order (MHLRO) in SiGeSn medium-entropy alloys (MEAs), in which the normalized
infinite-wavelength composition fluctuations for all three atomic species are
completely suppressed as in a perfect crystalline state. We show this MHLRO
naturally leads to the emergence of short-range order (SRO) recently
discovered in MEAs, which results in stable lower-energy states compared to
alloy models with random or special quasi-random structures (SQSs) possessing
no atomic SROs. The MHLRO MEAs approximately realize the Vegard’s law, which
offers a rule-of-mixture type predictions of the lattice constants and
electronic band gap, and thus can be considered as an ideal mixing state. The
MHLRO also directly gives rise to enhanced electronic band gaps and superior
thermal transport properties at low temperatures compared to random structures
and SQSs, which open up novel potential applications in optoelectronics and
thermoelectrics. Our analysis of the SiGeSn system leads to the formulation of
general organizing principles applicable in other medium- and high-entropy
alloys (HEAs), and a highly efficient computational model for rendering
realistic large-scale configurations of MEAs and HEAs.

Medium-entropy alloys (MEAs) are those composed of mixtures of equal or
relatively large proportions of three elements, and alloys containing five or
more elements are usually referred to as high-entropy alloys (HEAs) Zhang et
al. (2014). These MEAs and HEAs are distinctly different from traditional
metallic alloys which contain one or two major components with smaller amounts
of other elements. Studies have shown that MEAs and HEAs can possess superior
mechanical performance and resistance to corrosion compared to traditional
alloys Miracle and Senkov (2017), and can have unusual thermal and electronic
transport properties Mu et al. (2018), opening up their new device
applications.

Atomic structures of medium/high entropy alloys largely remain a mystery, due
to challenges in direct imaging of bulk structures in three dimensions (3D).
The preponderance of previous studies typically assumed MEAs and HEAs possess
a random or special quasi-random structure (SQS) Zunger et al. (1990)
consisting of atoms of different types distributed on the sites of a
crystalline lattice. However, a number of recent studies, in particular those
on the chromium-cobalt-nickel (CrCoNi) MEAs, have discovered significant
short-range order, often manifested as suppressed clustering of atoms of the
same type, which are shown to lead to energetically favorable states of the
alloys Ding et al. (2018); Zhang et al. (2020); Walsh et al. (2021). It is
interesting to note that a natural consequence of the existence of such short-
range order is the suppression of composition (or density) fluctuations. Since
the atoms prefer certain neighbor arrangements, it can be expected that as an
observation window moves randomly within the disordered system, the number of
atoms of a particular element that fall within the window would not fluctuate
significantly, especially for large window sizes.

A system in which the infinite-wavelength normalized density fluctuations are
completely suppressed is called hyperuniform Torquato and Stillinger (2003).
Among the hyperuniform systems, the disordered hyperuniform (DHU) varieties
are recently discovered exotic states of matter Torquato and Stillinger
(2003); Torquato (2018); Donev et al. (2005); Zachary et al. (2011); Jiao and
Torquato (2011); Jiao et al. (2014) that lie between a perfect crystal and
liquid. These systems are similar to liquids or glasses in that they are
statistically isotropic and possess no Bragg peaks, and yet they completely
suppress normalized large-scale density fluctuations like crystals and in this
sense possess a hidden long-range order Torquato and Stillinger (2003);
Zachary and Torquato (2009); Torquato (2018). Hyperuniformity is manifested as
vanishing static structure factor $S(k)$ in the zero-wavenumber limit, i.e.,
$\lim_{k\rightarrow 0}S(k)=0$ (where $k$ is the wavenumber), or equivalently
vanishing normalized number variance infinite-wavelength limit, i.e.,
$\lim_{R\rightarrow\infty}\sigma^{2}_{N}(R)/N(R)=0$, where $\sigma_{N}^{2}(R)$
and $N(R)$ are respectively the local number variance and average atom number
in a spherical observation window of radius $R$ Torquato and Stillinger
(2003); Torquato (2018). DHU states have been discovered in a variety of
equilibrium and non-equilibrium physical and biological systems, and appear to
endow such systems with desirable physical properties that cannot be achieved
in either the totally random or perfectly crystalline states Zachary et al.
(2011); Jiao and Torquato (2011); Chen et al. (2014); Zachary and Torquato
(2011); Torquato et al. (2015); Donev et al. (2005); Batten et al. (2009);
Kurita and Weeks (2011); Hunter and Weeks (2012); Dreyfus et al. (2015);
Hexner and Levine (2015); Jack et al. (2015); Weijs et al. (2015); Florescu et
al. (2009); Man et al. (2013); Jiao et al. (2014); Mayer et al. (2015); Hejna
et al. (2013); Klatt et al. (2019); Lei et al. (2019); Lei and Ni (2019); Chen
and Torquato (2018); Chremos and Douglas (2018); Rumi et al. (2019); Huang et
al. (2021); Torquato (2021); Jiao (2022).

We are thus motivated to ask the following fundamental question: Do MEAs (or
HEAs) with significant short-range order also possess the aforementioned
hidden hyperuniform long-range order? In this letter, using the ternary
silicon-germanium-tin (SiGeSn) MEA as an example, we provide strong numerical
evidence for a hidden multihyperuniform long-range order (MHLRO) in MEAs, in
which the normalized infinite-wavelength composition fluctuations for all
atomic species are completely suppressed as in a perfect crystalline state.
This is achieved by numerically rendering multihyperuniform realizations of
the alloy, which naturally leads to the emergence of short-range order (SRO)
quantitatively consistent with recently reported observations and results in
stable lower-energy states compared to alloys possessing random structure and
SQS without any atomic SRO Ding et al. (2018); Zhang et al. (2020); Walsh et
al. (2021). We note that multihyperuniform patterns have been discovered
previously in other physical and biological systems, e.g., chicken eyes Jiao
et al. (2014). The MHLRO MEAs approximately realize Vegard’s law Vegard
(1921), which assumes an ideal mixture, and predicts that for a multi-
component solid solution where the individual component in its pure form
shares the same crystal structure, the lattice parameter and electronic band
gap of the mixture should be the weighted averages of the constituents’
lattice parameters and band gaps. In particular, the MHLRO state minimizes the
lattice distortion compared to the prediction of Vegard’s law, and thus can be
considered as an ideal mixing state. We also demonstrate that the MHLRO
directly gives rise to enhanced electronic band gaps and superior thermal
transport properties at low temperatures compared to the random/SQS alloy
models, which open up many novel potential applications in optoelectronics and
thermoelectrics. Based on our numerical analysis of SiGeSn system, we
formulate general organizing principles predicting similar MHLRO in other MEAs
and HEAs.

Realizations and structural characteristics of multihyperuniform MEAs. We
devise a highly efficient general Fourier-space construction technique to
generate large samples of disordered multihyperuniform MEAs (or HEAs) as well
as those with SQSs. Specifically, we define a fictitious “energy” $E$ of the
system as the squared differences between the target and constructed structure
factors associated with different elements, i.e.,

$E=\sum_{i}\sum_{k}[S_{i}(k)-S_{i,0}(k)]^{2},$ (1)

where $S_{i,0}(k)$ and $S_{i}(k)$ are the angular-averaged structure factor
associated with type-$i$ atoms in the target and constructed structures. We
employ the simulated annealing procedure consisting of sequential Monte Carlo
(MC) moves that involve swapping particles of different types to obtain low-
energy states associated with this “energy” $E$, and drive the system under
construction towards one with targeted structure factors (see Supplemental
Materials for detail).

To obtain multihyperuniform MEAs or HEAs, we impose the strong condition of
stealthy hyperuniformity and require that the structure factor associated with
each atomic species is zero for a range of small $k$, i.e.,

$S_{i,0}(k)=0~{}~{}~{}~{}\textnormal{for}~{}~{}k<K_{i,0},$ (2)

where $K_{i,0}$ is the range of exclusion region in the Fourier space for the
target structure factor associated with element $i$. To obtain SQSs for MEAs
or HEAs that faithfully retain the salient ensemble statistics of random
mixtures, when considering type-$i$ atoms, we treat all the other types of
atoms as “voids”, and utilize a previously known formula for structure factor
of lattice structures with voids to obtain the target structure factor
associated with element $i$:

$S_{i,0}(k)=p_{i}+(1-p_{i})S_{l}(k)=1-c_{i}+c_{i}S_{l}(k),$ (3)

where the “void” concentration $p_{i}=1-c_{i}$, $c_{i}$ is the atomic fraction
of element $i$, and $S_{l}(k)$ is the structure factor associated with the
underlying lattice structure that the MEA or HEA is based on.

In this letter, as an example we consider equal-concentration ternary SiGeSn
alloys that were known to possess a diamond lattice structure, i.e.,
$c_{\mathrm{Si}}=c_{\mathrm{Ge}}=c_{\mathrm{Sn}}=1/3$. The SROs of such alloys
have not been systematically explored before, and for simplicity we assume
that $S_{i}(k)$ are the same for all three atom types, and set $K_{i,0}=K_{0}$
for all $i$ when generating multihyperuniform SiSeSn. This general setting
also makes our structural model immediately applicable to other MEAs. To
suppress the generation of periodic structures and at the same time ensure
multihyperuniformity, we further pick $K_{p}<K_{0}<K_{1}$, where $K_{p}$ is
the wavenumber associated with the smallest possible periodic unit compatible
with the prescribed atomic fractions, and $K_{1}$ is the location of the first
Bragg peak associated with the underlying lattice (see Supplemental Materials
for detailed explanation). We provide visualizations of representative SiGeSn
structures with MHLRO in Fig. 1, and structure factors associated with
different elements for $K_{0}a=3.6$ in Fig. 2, where $a$ is the lattice
constant of a diamond cubic unit. As expected, the structure factor $S_{i}(k)$
decreases as $k$ decreases at small $k$ and essentially approaches zero as $k$
goes to zero, indicating hyperuniformity of the atomic distribution for all
element types. However, the stealthiness is not realized, and $S_{i}(k)$
possesses small, but finite values at small $k<K_{0}$ in the exclusion region.
The lack of stealthiness is likely due to the frustration caused by the
presence of three element types as they simultaneously try to form a
hyperuniform pattern Jiao et al. (2014).

$\begin{array}[]{c}\\\ \includegraphics[width=216.81pt]{fig1.png}\end{array}$

Figure 1: (Color online) Two-dimensional slices along the (001) plane of
representative multihyperuniform (with an exclusion region of $K_{0}a=3.6$ for
targeted structure factor, where $a$ is the lattice constant of a diamond
cubic unit), random, and special quasirandom structures of SiGeSn alloys on a
9 $\times$ 9 $\times$ 9 diamond lattice. Figure 2: (Color online) Structure
factors (left column) and Warren-Cowley short-range order parameters
$\alpha_{ij}^{\nu}$ (right column) of multihyperuniform (with an exclusion
region of $K_{0}a=3.6$ for targeted structure factor), random, and special
quasirandom structures of SiGeSn alloys. The statistics are averaged over 5
configurations with $N=5832$ particles.

As a comparison, we also generate corresponding random mixtures and SQSs for
such alloys, as shown in Fig. 1. As expected, the structure factors $S_{i}(k)$
of random mixtures and SQSs are almost indistinguishable, and are essentially
flat and proportional to $1-c_{i}$ away from the Bragg peaks, although the
random mixtures possess larger statistical fluctuations, as shown in left
column of Fig. 2. To quantify the SRO of multihyperuniform, random, and SQS
alloys, we compute the Warren-Cowley SRO parameter $\alpha_{ij}^{\nu}$ Cowley
(1950) defined as

$\alpha_{ij}^{\nu}=1-\frac{p_{ij}^{\nu}}{c_{j}},$ (4)

where $p_{ij}^{\nu}$ is the probability of finding atomic species $j$ around
an atom of type $i$ in the $\nu$-th neighboring shell, and the results are
shown in right column of Fig. 2. Clearly, the multihyperuniform structures
possess substantial SRO, and atoms of the same type are unfavored in the first
three shells around any atom in the system, while random and SQS structures
possess essentially no SRO, which are expected. The natural emergence of SRO
in MEAs with hidden MHLRO suggest that these two types of order may be
coupled, and could be due to the fact that in order to achieve suppression of
normalized composition fluctuations on large length scales for each species
under lattice constraints, the system needs to suppress local clustering of
the same type of atoms. This rationale is not unique to the SiGeSn alloy, and
given the ubiquitous nature of SRO discovered previously in MEAs and HEAs,
strongly indicates the possibility of hidden MHLRO in other MEAs or HEAs.

MHLRO leads to lower-energy states and smaller lattice distortions in MEAs.
Conventional Si-Ge-Sn alloys were believed to possess a single solid solution
phase with amorphous atomic structures, consistent with the absence of any
binary or ternary intermetallic compounds in the Si-Ge-Sn phase diagrams
Olesinski and Abbaschian (1984a, b, c). However, it remains a long-standing
question whether any sort of atomic ordering exists in Si-Ge-Sn alloys
Mukherjee et al. (2017). For example, short-ranger order was recently
predicted in Ge-Sn alloys Cao et al. (2020), while previous Monte Carlo
simulations that only explore the vicinity of random mixtures/SQSs in the
configuration space fail to discover low-energy states with short-range order
Wang et al. (2019). Here we address this enduring issue by comparing the
energies of our DHU structures with those of random and SQS systems. In
particular, we perform molecular dynamics simulations with the Stillinger-
Weber potential parameterized for Si-Ge-Sn alloys Stillinger and Weber (1985);
Lee and Hwang (2012, 2017), which has been used to study thermal transport of
these alloys Wang et al. (2020), and compute the energy difference using the
DHU system’s energy as the reference. We perform calculations for different
supercell sizes ($N=5832$ and $N=46656$ atoms), and benchmark the energy
differences with the density functional theory (see the Supplemental Material
for calculation parameters) results of systems with $N=216$ atoms. Figure 3(a)
shows that although the magnitudes of the energy differences are both method
and size-dependent, all the results consistently show that (i) the energies of
the SQS and random systems are similar; more importantly, (ii) the
multihyperuniform system is the most energetically stable, which strongly
suggests the existence of long-range order in Si-Ge-Sn alloys.

We speculate that multihyperuniform structures with suppressed large-scale
composition fluctuations should be associated with a smaller degree of
deviation from Vegard’s law as well as smaller lattice distortion that is
intrinsic to any HEA or MEA system He and Yang (2018) compared to random/SQS
structures. We optimize the structures using DFT simulations and determine the
lattice constants of the optimized multihyperuniform, random, and SQS
structures with $3\times 3\times 3$ supercells to be 17.903, 17.883, and
17.883 Å, respectively. According to Vegard’s law, the lattice constant of an
equi-molar Si-Ge-Sn alloy is predicted to be 17.904 Å, the weighted average of
the lattice constants of pure Si (16.406 Å), Ge (17.346 Å), and Sn (19.960 Å)
crystals with $3\times 3\times 3$ supercells. These results confirm our
hypothesis that the multihyperuniform system exhibits an minimal amount of
deviation from Vegard’s law and reflects the nearly ideal mixing of the three
constituent elements in the multihyperuniform SiGeSn MEA. We also compute the
lattice distortion using the vector dissimilarity metric Zhuang (2021), i.e.,
the root mean squared displacement (RMSD) between the reference lattice (a
diamond lattice with a lattice constant as predicted by Vegard’s law) and the
optimized structure for the multihyperuniform, random, and SQS systems,
respectively. We find the RMSDs for the multihyperuniform, random, and SQS
systems are 2.14, 3.21, and 3.31 Å, respectively. Clearly, the
multihyperuniform system shows the least extent of lattice distortion, which
leads to the lowest ground-state energy of such system. The stability of the
multihyperuniform system can be attributed to the suppressed composition
fluctuations in such system. On the other hand, the RMSDs and energies for
random and SQS systems are nearly the same, as shown in Figure 3(a), which is
consistent of the fact that the individual SQSs by design should retain the
salient ensemble statistics of random mixtures.

Figure 3: (a) Histogram comparing the energies (on a log10 scale) of random,
SQS, and multihyperuniform structures using DFT and MD simulations with the
Stillinger-Weber potential Stillinger and Weber (1985) and cell sizes. Here we
use the multihyperuniform structure as the reference and set its energy to
zero; in the plot “DFT-216” refers to DFT results with a supercell of 216
atoms, and “SW-5832” and “SW-46656” refer to MD results with a supercell of
5832 and 46656 atoms, respectively. Band structures of the multihyperuniform
(b), SQS (c), and random (d) Wang et al. (2020) SiGeSn MEAs. (e) Phonon
density of states and partition ratio of SiGeSn MEA with the multihyperuniform
structure.

Larger electronic bandgaps. The effect of the presence of the long-range order
on the electronic structure and thermal transport properties of the SiGeSn MEA
remains elusive, which we address in our current work. We first calculate the
band structure of the SiGeSn MEA with the multihyperuniform and SQS structure
and use the previously computed band structure of the random structures for
comparison Wang et al. (2020). Figures 3(c), (d), and (e) reveals that all the
three systems have direct band gaps at the $\Gamma$ point and the
corresponding band gaps increase from 0.28, 0.36, to 0.50 eV for the random,
SQS, and multihyperuniform systems, respectively. For pure crystalline Si, Ge,
and Sn, the band gaps are 1.29, 0.33, and 0.15 eV, respectively. In the random
SiGeSn system, the clustering of Ge and Sn atoms significantly reduces the
band gap and the presence of MHLRO redistributes all atoms with suppressed
local clustering of atoms of the same type, leading to a band gap that is
nearly the average band gap (0.59 eV) of the three constituent elements, which
is another manifestation that the MHLRO SiGeSn alloy approximately realizes
Vegard’s law.

Enhanced low-temperature thermal conductivity. The thermal conductivity of the
multihyperuniform system at low temperatures are generally higher than that of
the random mixture. For example, the thermal conductivities are 1.53 and 0.87
W/m$\cdot$ K at 10 K for the multihyperuniform and random systems,
respectively. The improved thermal conductivity makes the multihyperuniform
SiGeSn MEA a better candidate for electronic device applications such as
computer chips that require a stricter management in heat dissipation. The
higher thermal conductivity of the multihyperuniform system once again can be
attributed to the smaller lattice distortions of such system compared to the
random mixture, since it is well known that in solid solutions such as alloys,
large lattice distortions would enhance the scattering of phonons and lead to
low thermal conductivity Fan et al. (2016). However, as the temperature
increases, the effect of phonon-phonon coupling becomes more prominent, which
leads to similar and decreased thermal conductivities for multihyperuniform
and random mixtures. For example, we find that the computed thermal
conductivities at 400 K for the multihyperuniform and random systems are 0.53
and 0.57 W/m$\cdot$ K, respectively. Figure 3(e) shows the phonon density of
states (PDOS) and partition ratio (PR) of the multihyperuniform mixture.
Although the overall shape of the PDOS is similar to that of the random system
Wang et al. (2020), we observe a gap almost opens around 280 cm-1
(corresponding to a temperature of about 400 K), leading to diminishing PR
modes in the vicinity of this frequency, and a reduction in the thermal
conductivity.

In summary, we present strong numerical evidence suggesting a new universal
MHLRO in MEA and HEA systems, exemplified by the SiGeSn MEA that is endowed
with a minimal lattice distortion and enhanced electronic band gap, both of
which matching the predictions from Vegard’s law that reflects the ideal
mixing behavior of atoms in the MEAs. We also find improved thermal
conductivity of the multihyperuniform alloy at low temperatures due to smaller
lattice distortions compared to random mixtures. Our analysis of the SiGeSn
system leads to the formulation of general organizing principles applicable in
other MEAs and HEAs, and a new type of highly efficient computational model
for rendering realistic large-scale configurations of MEAs and HEAs beyond the
conventional random mixture/SQS models. Our results will open up many novel
potential applications in optoelectronics, thermoelectrics, and quantum
devices.

###### Acknowledgements.

X. J., D. W., and H.Z. thank Arizona State University (ASU) for the start-up
funds. This research used computational resources of the Agave Research
Computer Cluster of ASU and the Texas Advanced Computing Center under Contract
TG-DMR170070.

## References

  * Zhang et al. (2014) Y. Zhang, T. T. Zuo, Z. Tang, M. C. Gao, K. A. Dahmen, P. K. Liaw, and Z. P. Lu, Prog. Mater. Sci. 61, 1 (2014).
  * Miracle and Senkov (2017) D. Miracle and O. Senkov, Acta Mater. 122, 448 (2017).
  * Mu et al. (2018) S. Mu, Z. Pei, X. Liu, and G. M. Stocks, J. Mater. Sci. 33, 2857â€“2880 (2018).
  * Zunger et al. (1990) A. Zunger, S.-H. Wei, L. G. Ferreira, and J. E. Bernard, Phys. Rev. Lett. 65, 353 (1990).
  * Ding et al. (2018) J. Ding, Q. Yu, M. Asta, and R. O. Ritchie, Proc. Natl. Acad. Sci. U.S.A. 115, 8919 (2018).
  * Zhang et al. (2020) R. Zhang, S. Zhao, J. Ding, Y. Chong, T. Jia, C. Ophus, M. Asta, R. O. Ritchie, and A. M. Minor, Nature 581, 283 (2020).
  * Walsh et al. (2021) F. Walsh, M. Asta, and R. O. Ritchie, Proc. Natl. Acad. Sci. U.S.A. 118 (2021).
  * Torquato and Stillinger (2003) S. Torquato and F. H. Stillinger, Phys. Rev. E 68, 041113 (2003).
  * Torquato (2018) S. Torquato, Phys. Rep. 745, 1 (2018).
  * Donev et al. (2005) A. Donev, F. H. Stillinger, and S. Torquato, Phys. Rev. Lett. 95, 090604 (2005).
  * Zachary et al. (2011) C. E. Zachary, Y. Jiao, and S. Torquato, Phys. Rev. Lett. 106, 178001 (2011).
  * Jiao and Torquato (2011) Y. Jiao and S. Torquato, Phys. Rev. E 84, 041309 (2011).
  * Jiao et al. (2014) Y. Jiao, T. Lau, H. Hatzikirou, M. Meyer-Hermann, J. C. Corbo, and S. Torquato, Phys. Rev. E 89, 022721 (2014).
  * Zachary and Torquato (2009) C. E. Zachary and S. Torquato, J. Stat. Mech. Theor. Exp. 2009, P12015 (2009).
  * Chen et al. (2014) D. Chen, Y. Jiao, and S. Torquato, J. Phys. Chem. B 118, 7981 (2014).
  * Zachary and Torquato (2011) C. E. Zachary and S. Torquato, Phys. Rev. E 83, 051133 (2011).
  * Torquato et al. (2015) S. Torquato, G. Zhang, and F. H. Stillinger, Phys. Rev. X 5, 021020 (2015).
  * Batten et al. (2009) R. D. Batten, F. H. Stillinger, and S. Torquato, Phys. Rev. Lett. 103, 050602 (2009).
  * Kurita and Weeks (2011) R. Kurita and E. R. Weeks, Phys. Rev. E 84, 030401(R) (2011).
  * Hunter and Weeks (2012) G. L. Hunter and E. R. Weeks, Rep. Prog. Phys. 75, 066501 (2012).
  * Dreyfus et al. (2015) R. Dreyfus, Y. Xu, T. Still, L. A. Hough, A. G. Yodh, and S. Torquato, Phys. Rev. E 91, 012302 (2015).
  * Hexner and Levine (2015) D. Hexner and D. Levine, Phys. Rev. Lett. 114, 110602 (2015).
  * Jack et al. (2015) R. L. Jack, I. R. Thompson, and P. Sollich, Phys. Rev. Lett. 114, 060601 (2015).
  * Weijs et al. (2015) J. H. Weijs, R. Jeanneret, R. Dreyfus, and D. Bartolo, Phys. Rev. Lett. 115, 108301 (2015).
  * Florescu et al. (2009) M. Florescu, S. Torquato, and P. J. Steinhardt, Proc. Natl. Acad. Sci. U.S.A. 106, 20658 (2009).
  * Man et al. (2013) W. Man, M. Florescu, E. P. Williamson, Y. He, S. R. Hashemizad, B. Y. C. Leung, D. R. Liner, S. Torquato, P. M. Chaikin, and P. J. Steinhardt, Proc. Natl. Acad. Sci. U.S.A. 110, 15886 (2013).
  * Mayer et al. (2015) A. Mayer, V. Balasubramanian, T. Mora, and A. M. Walczak, Proc. Natl. Acad. Sci. USA 112, 5950 (2015).
  * Hejna et al. (2013) M. Hejna, P. J. Steinhardt, and S. Torquato, Phys. Rev. B 87, 245204 (2013).
  * Klatt et al. (2019) M. A. Klatt, J. Lovrić, D. Chen, S. C. Kapfer, F. M. Schaller, P. W. A. Schönhöfer, B. S. Gardiner, A. Smith, G. E. Schröder-Turk, and S. Torquato, Nat. Commun. 10, 1 (2019).
  * Lei et al. (2019) Q.-L. Lei, M. P. Ciamarra, and R. Ni, Sci. Adv. 5, eaau7423 (2019).
  * Lei and Ni (2019) Q.-L. Lei and R. Ni, Proc. Natl. Acad. Sci. U.S.A. 116, 22983 (2019).
  * Chen and Torquato (2018) D. Chen and S. Torquato, Acta Mater. 142, 152 (2018).
  * Chremos and Douglas (2018) A. Chremos and J. F. Douglas, Phys. Rev. Lett. 121, 258002 (2018).
  * Rumi et al. (2019) G. Rumi, J. A. Sánchez, F. Elías, R. C. Maldonado, J. Puig, N. R. C. Bolecek, G. Nieva, M. Konczykowski, Y. Fasano, and A. B. Kolton, Phys. Rev. Res. 1, 033057 (2019).
  * Huang et al. (2021) M. Huang, W. Hu, S. Yang, Q.-X. Liu, and H. P. Zhang, Proc. Natl. Acad. Sci. U.S.A. 118, e2100493118 (2021).
  * Torquato (2021) S. Torquato, Proc. Natl. Acad. Sci. U.S.A. 118, e2107276118 (2021).
  * Jiao (2022) Y. Jiao, Phys. A: Stat. Mech. Appl. 585, 126435 (2022).
  * Vegard (1921) L. Vegard, Z. Phys. 5, 17 (1921).
  * Cowley (1950) J. M. Cowley, Phys. Rev. 77, 669 (1950).
  * Olesinski and Abbaschian (1984a) R. Olesinski and G. Abbaschian, Bull. Alloy Phase Diagr. 5, 273 (1984a).
  * Olesinski and Abbaschian (1984b) R. Olesinski and G. Abbaschian, Bull. Alloy Phase Diagr. 5, 265 (1984b).
  * Olesinski and Abbaschian (1984c) R. Olesinski and G. Abbaschian, Bull. Alloy Phase Diagr. 5, 180 (1984c).
  * Mukherjee et al. (2017) S. Mukherjee, N. Kodali, D. Isheim, S. Wirths, J. M. Hartmann, D. Buca, D. N. Seidman, and O. Moutanabbir, Phys. Rev. B 95, 161402 (2017).
  * Cao et al. (2020) B. Cao, S. Chen, X. Jin, J. Liu, and T. Li, ACS Appl. Mater. Interfaces 12, 57245 (2020).
  * Wang et al. (2019) D. Wang, L. Liu, W. Huang, and H. L. Zhuang, J. Appl. Phys. 126, 225703 (2019).
  * Stillinger and Weber (1985) F. H. Stillinger and T. A. Weber, Phys. Rev. B 31, 5262 (1985).
  * Lee and Hwang (2012) Y. Lee and G. S. Hwang, Phys. Rev. B 85, 125204 (2012).
  * Lee and Hwang (2017) Y. Lee and G. S. Hwang, J. Phys. D: Appl. Phys. 50, 494001 (2017).
  * Wang et al. (2020) D. Wang, L. Liu, M. Chen, and H. Zhuang, Acta Mater. 199, 443 (2020).
  * He and Yang (2018) Q. He and Y. Yang, Front. Mater. Sci. 5, 42 (2018).
  * Zhuang (2021) H. Zhuang, arXiv preprint arXiv:2110.03797 (2021).
  * Fan et al. (2016) Z. Fan, H. Wang, Y. Wu, X. J. Liu, and Z. P. Lu, Rsc Adv. 6, 52164 (2016).


# On the transport equation for probability density functions of turbulent
vorticity fields

Jiawei Li, Zhongmin Qian and Mingrui Zhou Department of Mathematical Sciences,
Carnegie Mellon University, Pittsburgh, PA 15213, United States. _Email_
<EMAIL_ADDRESS>Institute, University of Oxford, Oxford
OX2 6GG, United Kingdom, and Oxford Suzhou Center for Advanced Research.
_Email_<EMAIL_ADDRESS>Institute, University of
Oxford, Oxford OX2 6GG, United Kingdom. _Email_<EMAIL_ADDRESS>
###### Abstract

Vorticity random fields of turbulent flows (modeled over the vorticity
equation with random initial data for example) are singled out as the main
dynamic variables for the description of turbulence, and the evolution
equation of the probability density function (PDF) of the vorticity field has
been obtained. This PDF evolution equation is a mixed type partial
differential equation (PDE) of second order which depends only on the
conditional mean (which is a first order statistics) of the underlying
turbulent flow. This is in contrast with Reynolds’ mean flow equation which
relies on a quadratic statistics. The PDF PDE may provide new closure schemes
based on the first order conditional statistics, and some of them will be
described in the paper. We should mention that the PDF equation is interesting
by its own and is worthy of study as a PDE of second order.

_Keywords_ : Navier-Stokes equation, PDF method, turbulent flows, vorticity.

_MSC Classifications_ : 76F02, 76F55, 76D05, 76M35

##  1 Introduction

In statistical fluid mechanics (cf. [13]), the velocity $U(x,t)$ of a
turbulent flow is promoted to a random field, cf. [8], indexed by space
variable $x\in\mathbb{R}^{3}$ and time parameter $t$. From this point of view,
the turbulence problem, if there is one, seeks for a description of the
distribution of the velocity field. This distribution is rather complicated
and consists of all joint distributions of the velocity across over finitely
many locations and times, and therefore it is challenging to describe the
distributions of turbulent flows in general. As early as in 1950’s, E. Hopf
[6] (cf. [13] too) made an ambitious attempt and derived differential
equations for the distribution of a turbulent flow. Hopf’s differential
equations are however infinite dimensional and involve functional derivatives,
which are therefore too difficult to extract useful information about
turbulence. The dominant approach in the statistical theory of turbulence,
initiated in G. I. Taylor’s seminal work [18, 19], has been based on the
analysis of moment structure functions and relied on the spectral method (cf.
[1, 13] for example). In the past decades, attention has been paid to the one-
dimensional marginal distribution of the velocity field $U(x,t)$, whose
distribution is finite-dimensional. It is reasonable to assume that the
distribution of $U(x,t)$, with $(x,t)$ fixed, has a probability density
function (PDF) with respect to the Lebesgue measure on $\mathbb{R}^{3}$. The
PDF methods based on formal PDF transport equations have been developed
through a series of work by Pope and other researchers (cf. [14, 15]), which
become powerful tools for modeling turbulent flows. The PDF, by definition,
has to satisfy “the adjoint equation” of the fluid dynamic equations. If the
turbulent flow is an incompressible viscous fluid flow, then the PDF of the
velocity must satisfy the adjoint equation of the Navier-Stokes equation. The
transport equation for PDF known in the literature is still a formal adjoint
equation of the Navier-Stokes equation, and therefore only few features from
the transport equation may be used in modeling turbulent flows.

For incompressible fluid flows, the Navier-Stokes equation is equivalent to
the vorticity equation (see (2) below), and therefore it is natural to
consider the vorticity $W=\nabla\wedge U$ as the main fluid dynamic variable
for the study of turbulence [12]. There are good reasons why we should
concentrate on the vortex motion in turbulence. The velocity of a turbulent
flow is unlikely to be independently or possess conditional independence with
respect to the spatial variable, while vortex motions of many turbulent flows
observed in nature (such as vortex lines, vortex rings) acquire certain
conditional independence, in the sense that by focusing on the motion near a
fixed region, the future vortex motions evolve more or less independent of
what happens in other positions. In fact there are good evidences which
demonstrate that some sort of superposition property of vorticity may be
maintained, although not exactly due to highly non-linear and non-local nature
of turbulence. These observations are valuable in modeling turbulent flows via
the vorticity, which are already applied in vortex methods (cf. [3, 11]). In
this sense PDF methods based on the vorticity are valuable.

The main contribution of the present work is the partial differential equation
for the PDF of the vorticity (called PDF PDE, or PDF equation for short),
which will be derived in the main body of the paper. The PDF PDE is highly
non-linear, however the most striking aspect is that the PDF PDE for the
vorticity depends only on one single first-order statistical characteristic of
the turbulent flow. More precisely, we identify the main statistical
characteristic needed for the PDF PDE with the conditional mean function of
the increment of the vorticity given the current vorticity:

$\mu^{i}(x,y,w,t)=\mathbb{E}\left[\left.W^{i}(y,t)-W^{i}(x,t)\right|W(x,t)=w\right]$
(1)

for $i=1,2,3$, where $x,y,w\in\mathbb{R}^{3}$ and $t\geq 0$. The PDF of the
vorticity $W$ is a solution to the PDF PDE, which is a second-order partial
differential equation where coefficients appearing in the PDF PDE depend on
$\mu$ only. Besides its theoretical interest, this PDF PDE paves the way
towards practical modeling the statistics of turbulent flows based on the
vorticity PDF.

The paper is organized as the following. In Section 2, several notions and
notations together with several standard assumptions about fluid dynamic
random fields will be introduced, and the main result, namely the PDF PDE for
vorticity will be derived. In Sections 3 and 4, the theoretical foundation,
based on our PDF PDE for vorticity, will be laid for the purpose of modeling
various turbulent flows. In Section 5, we propose the most direct way of
modeling PDF of the vorticity by specifying the conditional mean function
$\mu^{i}$. The PDF PDE considered purely as a PDE theory is over-determined,
due to its non-linearlity in the sense that the coefficients appearing in the
PDF PDE are not independent of its solutions. Therefore care is needed to
ensure that the additional constraint is satisfied. In the last section, we
propose the heat flow method to model statistical quantities needed for
closing the PDF PDE and obtain concrete PDF examples.

_Conventions on notations_. The following conventions are used throughout the
paper. Firstly Einstein’s convention on summation on repeated indices through
their ranges is assumed, unless otherwise specified. If $A$ is a vector or a
vector field (in the space of dimension three) dependent on some parameters,
then its components are labeled with upper-script indices, i.e.
$A=(A^{i})=(A^{1},A^{2},A^{3})$. The same convention applies to coordinates
too. Partial derivatives of functions may be labeled with variables in sub-
scripts. For example, if $A(w;x,t)$ is a vector-valued function depending on
$w,x\in\mathbb{R}^{3}$ and $t$, then $\nabla_{w}A$ means the total derivatives
$(\frac{\partial}{\partial w^{i}}A^{j})$, $\Delta_{x}A$ means the vector
$(\frac{\partial^{2}}{\partial x^{k}\partial x^{k}}A^{i})$ of the Laplacians
of $A^{i}$. However, as a general rule, derivatives without sub-scripts mean
the derivatives with respect to the variable $x=(x^{i})$, unless otherwise
specified for avoiding possible confusion.

The velocity vector field will be denoted by $U$. $W=\nabla\wedge U$ is its
vorticity so that $W^{i}=\varepsilon^{ijk}\frac{\partial}{\partial
x^{j}}U^{k}$, where $\varepsilon^{ijk}$ are the Levi-Civita symbols.

##  2 PDF equation for the vorticity

In this paper we regard the vorticity of an incompressible turbulent flow as
the main dynamic variable. The goal of this section is to derive the evolution
equation for the PDF of the vorticity of an incompressible turbulent flow.

###  2.1 Prelims and assumptions

Let $U(x,t)$ be the velocity of an incompressible turbulent flow in
$\mathbb{R}^{3}$ with viscosity $\nu$. Suppose there is no external force
supplied to the turbulence. Being a random field though, $U(x,t)$ satisfies
the Navier-Stokes equations

$\frac{\partial U^{i}}{\partial t}+U^{j}\frac{\partial U^{i}}{\partial
x^{j}}=\nu\Delta U^{i}-\frac{\partial P}{\partial x^{i}}$

and

$\frac{\partial U^{j}}{\partial x^{j}}=0,$

where $i=1,2,3$, and $P$ is the pressure, subject to initial conditions which
are random. The motion equations of the vorticity $W=\nabla\wedge U$ are the
vorticity equations

$\frac{\partial W^{i}}{\partial t}+U^{j}\frac{\partial W^{i}}{\partial
x^{j}}=\nu\Delta W^{i}+W^{j}\frac{\partial U^{i}}{\partial x^{j}}$ (2)

for $i=1,2,3$.

We make the following technical assumptions.

First, we assume that both $U(x,t)$ and $W(x,t)$ have derivatives in $x$ and
$t$ of any order, and theses derivatives decay to zero at infinity
sufficiently fast, so that possible boundary terms arising in applications of
the Stokes’ formula have no contributions in computations below. Therefore,
since $U$ is divergence-free, $\Delta U=-\nabla\wedge W$, according to Green’s
formula

$U^{i}(x,t)=\int_{\mathbb{R}^{3}}\frac{1}{4\pi|x-y|}\varepsilon^{ijk}\frac{\partial}{\partial
y^{j}}W^{k}(y,t)\mathrm{d}y,$

which yields the Biot-Savart law:

$U^{i}(x,t)=-\int_{\mathbb{R}^{3}}\varepsilon^{ijk}\frac{x^{j}-y^{j}}{4\pi|x-y|^{3}}W^{k}(y,t)\mathrm{d}y.$
(3)

for $i=1,2,3$.

Our second technical assumption is to impose certain regularity on the
distribution of the vorticity. At each point $x$ and instance $t\geq 0$,
$W(x,t)$ is a random variable defined on a probability space
$(\Omega,\mathcal{F},\mathbb{P})$, taking values in $\mathbb{R}^{3}$, and
therefore its law (or called distribution) is a probability measure on
$(\mathbb{R}^{3},\mathcal{B}(\mathbb{R}^{3}))$. We assume that $W(x,t)$ has a
positive and smooth probability density function (PDF), denoted by $f(w;x,t)$,
in the sense that

$\mathbb{P}\left(W(x,t)\in A\right)=\int_{A}f(w;x,t)\mathrm{d}w$

for every Borel measurable subset $A$. The joint distribution of $W(x,t)$ and
$W(y,t)$ for any pair $x\neq y$ has a positive and smooth PDF, denoted by
$f_{2}(w_{1},w_{2};x,y,t)$.

To state our third assumption, we first notice that

$f_{2|1}(w;w_{1},x,y,t)=\frac{f_{2}(w_{1},w;x,y,t)}{f(w_{1};x,t)}$

is the PDF of the conditional law of $W(y,t)$ given $W(x,t)=w_{1}$. The
conditional mean function $\mu^{i}(x,y,w,t)$ is defined by

$\displaystyle\mu^{i}(x,y,w,t)$
$\displaystyle=\mathbb{E}\left[\left.W^{i}(y,t)-W^{i}(x,t)\right|W(x,t)=w\right]$
$\displaystyle=\int_{\mathbb{R}^{3}}(w_{1}^{i}-w^{i})f_{2|1}(w_{1};w,x,y,t)\mathrm{d}w_{1}$
(4)

for $i=1,2,3$, which will play a dominant role in the sequel.

The third technical assumption is about the regularity of the conditional
average function $\mu=(\mu^{1},\mu^{2},\mu^{3})$. It is assumed that the
derivatives of $\mu$ of any order exist and decay to zero sufficiently fast at
infinity. Moreover it is assumed that $\mu$ has an asymptotic expansion

$\mu^{i}(x,y,w,t)=a_{k}^{i}(x,w,t)(y^{k}-x^{k})+b_{jk}^{i}(x,w,t)(y^{k}-x^{k})(y^{j}-x^{j})+o\left(|y-x|^{2}\right)$
(5)

as $|y-x|\rightarrow 0$, where $a_{k}^{i}$’s and $b_{jk}^{i}$’s are assumed to
be continuous with respect to all of their arguments. We denote
$b^{i}=b_{kk}^{i}$ for $i=1,2,3$. It is clear that

$a_{k}^{i}(x,w,t)=\mathbb{E}\left[\left.\frac{\partial W^{i}}{\partial
x^{k}}(x,t)\right|W(x,t)=w\right]=\left.\frac{\partial}{\partial
y^{k}}\mu^{i}(x,y,w,t)\right|_{y=x}$ (6)

which represents the local rate of change in the vortex motion over the
turbulent region, and

$b^{i}(x,w,t)=\frac{1}{2}\mathbb{E}\left[\Delta_{x}W^{i}(x,t)|W(x,t)=w\right]=\frac{1}{2}\left.\Delta_{y}\mu^{i}(x,y,w,t)\right|_{y=x}.$
(7)

In the remainder of the paper, we will work with a turbulent flow for which
the three assumptions listed above are satisfied.

###  2.2 PDF equation and its derivation

In this section, we derive the main result of the paper, that is, a partial
differential equation which the PDF of the vorticity must satisfy.

######  Theorem 1.

Under the assumptions and notations established in Sec. 2.1, suppose the PDF
$f(w;x,t)$ of the vorticity is smooth in $(w,x,t)$ and has finite moments,
that is $\int_{\mathbb{R}^{3}}|w|^{n}f(w;x,t)dw<\infty$ for every
$n=1,2,\cdots$. Then $f$ satisfies the following partial differential
equation:

$\left(\frac{\partial}{\partial t}+\frac{\partial B^{i}}{\partial
x^{i}}+B^{i}\frac{\partial}{\partial
x^{i}}-\nu\Delta_{x}\right)f=\nu\frac{\partial}{\partial
w^{i}}\left(\frac{\partial}{\partial
x^{k}}\left(fa_{k}^{i}\right)-2b^{i}f\right)+\frac{\partial}{\partial
w^{i}}\left(fD^{i}\right),$ (8)

where $\Delta_{x}$ is the Laplacian with respect to the space variable $x$,
$a_{k}$ and $b$ are given as in (6) and (7),

$B^{i}(x,w,t)=\int_{\mathbb{R}^{3}}\frac{1}{4\pi|y-x|}\varepsilon^{ijk}\frac{\partial}{\partial
y^{j}}\mu^{k}(x,y,w,t)\mathrm{d}y,$ (9)

and

$D^{i}(x,w,t)=w^{l}\int_{\mathbb{R}^{3}}\frac{y^{l}-x^{l}}{4\pi|y-x|^{3}}\varepsilon^{ijk}\frac{\partial}{\partial
y^{j}}\mu^{k}(x,y,w,t)\mathrm{d}y.$ (10)

###### Proof.

Let $F$ be a smooth function on $\mathbb{R}^{3}$ with a compact support. By
the definition of PDF one has

$\int_{\mathbb{R}^{3}}F(w)\frac{\partial}{\partial
t}f(w;x,t)\mathrm{d}w=\mathbb{E}\left[\frac{\partial}{\partial
t}F(W(x,t))\right].$

We are going to calculate the right-hand side expectation in terms of the PDF
$f(w;x,t)$ and other statistical characteristics of the turbulent flow. This
will be done by using the following equation

$\frac{\partial}{\partial
t}F(W(x,t))=\nu\Delta_{x}F(W)-\frac{\partial\left(U^{i}F(W)\right)}{\partial
x^{i}}+F_{j}(W)W^{i}\frac{\partial U^{j}}{\partial x^{i}}-\nu\frac{\partial
F_{j}(W)}{\partial x^{k}}\frac{\partial W^{j}}{\partial x^{k}},$

where, for simplicity, $F_{j}$ denote the partial derivatives $\frac{\partial
F}{\partial x^{j}}$ of $F$, $j=1,2,3$. This equation follows directly from the
vorticity equations. From the previous equation we obtain that

$\displaystyle\mathbb{E}\left[\frac{\partial}{\partial t}F(W(x,t))\right]$
$\displaystyle=\nu\mathbb{E}\left[\Delta_{x}F(W)\right]-\mathbb{E}\left[\frac{\partial\left(U^{i}F(W)\right)}{\partial
x^{i}}\right]$ $\displaystyle+\mathbb{E}\left[F_{j}(W)W^{i}\frac{\partial
U^{j}}{\partial x^{i}}\right]-\mathbb{E}\left[\nu\frac{\partial
F_{j}(W)}{\partial x^{k}}\frac{\partial W^{j}}{\partial x^{k}}\right]$
$\displaystyle:=I_{1}+I_{2}+I_{3}+I_{4}.$ (11)

Let us calculate $I_{i}$’s on the right-hand side of Eq. (11). First, it is
easy to see that by definition

$\displaystyle I_{1}$
$\displaystyle=\nu\Delta_{x}\mathbb{E}\left[F(W)\right]=\int_{\mathbb{R}^{3}}F(w)\nu\Delta_{x}f(w;x,t)\mathrm{d}w.$
(12)

For computing $I_{2}$, we shall use Eq. (3) and obtain that

$\displaystyle I_{2}$
$\displaystyle=-\mathbb{E}\left[\frac{\partial(U^{i}F(W))}{\partial
x^{i}}\right]=-\frac{\partial}{\partial
x^{i}}\mathbb{E}\left[F(W)U^{i}\right]$
$\displaystyle=-\frac{\partial}{\partial
x^{i}}\mathbb{E}\left[F(W)\varepsilon^{ijk}\int_{\mathbb{R}^{3}}\frac{1}{4\pi|y-x|}\frac{\partial}{\partial
y^{j}}W^{k}(y,t)\mathrm{d}y\right].$

To work out the expectation on the right-hand side, we may rewrite the partial
derivative as a limit:

$\frac{\partial}{\partial y^{j}}W^{k}(y,t)=\lim_{h\rightarrow
0}\frac{1}{h}\left(W^{k}(y+he^{(j)},t)-W^{k}(y,t)\right),$

where $e^{(j)}$ represents the unit vector with $j$-th component equal to $1$,
and the rest two components $0$, so that we can rewrite $I_{2}$ as the
following limit

$I_{2}=-\frac{\partial}{\partial x^{i}}\lim_{h\rightarrow
0}\frac{1}{h}\mathbb{E}\left[F(W)\varepsilon^{ijk}\int_{\mathbb{R}^{3}}\frac{1}{4\pi|y-x|}\left(W^{k}(y+he^{(j)},t)-W^{k}(y,t)\right)\mathrm{d}y\right].$
(13)

The expectation in this expression may be written in terms of two-point joint
distributions of $W$ as the following

$\int_{\mathbb{R}^{3}}F(w)\left[\int\frac{\varepsilon^{ijk}}{4\pi|y-x|}\left(\int_{\mathbb{R}^{3}}w_{1}^{k}\left(f_{2}(w,w_{1};x,y+he^{(j)},t)-f_{2}(w,w_{1};x,y,t)\right)\mathrm{d}w_{1}\right)\mathrm{d}y\right]\mathrm{d}w.$

The inner integral against the variable $w_{1}$ equals

$f(w;x,t)\left(\mu(x,y+he^{(j)},w,t)-\mu(x,y,w,t)\right).$

After substituting this in Eq. (13) and sending $h\rightarrow 0$, we obtain
that

$I_{2}=\int_{\mathbb{R}^{3}}F(w)\left[-\frac{\partial}{\partial
x^{i}}\left(f(w;x,t)\int_{\mathbb{R}^{3}}\frac{\varepsilon^{ijk}}{4\pi|y-x|}\frac{\partial}{\partial
y^{j}}\mu^{k}(x,y,w,t)\mathrm{d}y\right)\right]\mathrm{d}w.$ (14)

Next we deal with $I_{3}$. Again, using (3), we may write

$\frac{\partial}{\partial
x^{j}}U^{i}(x,t)=\int_{\mathbb{R}^{3}}\frac{y^{j}-x^{j}}{4\pi|y-x|^{3}}\varepsilon^{ilk}\frac{\partial}{\partial
y^{l}}W^{k}(y,t)\mathrm{d}y,$

so that

$\displaystyle I_{3}$
$\displaystyle=\mathbb{E}\left[F_{i}(W)W^{j}\frac{\partial U^{i}}{\partial
x^{j}}\right]$
$\displaystyle=\mathbb{E}\left[F_{i}(W)W^{j}\int_{\mathbb{R}^{3}}\frac{y^{j}-x^{j}}{4\pi|y-x|^{3}}\varepsilon^{ilk}\frac{\partial}{\partial
y^{l}}W^{k}(y,t)\mathrm{d}y\right],$

which can be evaluated by using the two-point joint PDF. Indeed, we may repeat
the same idea as in the computation of $I_{2}$, to obtain that

$I_{3}=\int_{\mathbb{R}^{3}}F(w)\frac{\partial}{\partial
w^{i}}\left[-f(w;x,t)\varepsilon^{ilk}w^{j}\int_{\mathbb{R}^{3}}\frac{y^{j}-x^{j}}{4\pi|y-x|^{3}}\frac{\partial}{\partial
y^{l}}\mu^{k}(x,y,w,t)\mathrm{d}y\right]\mathrm{d}w.$ (15)

Similarly, for the last term $I_{4}$, we write

$\displaystyle I_{4}$
$\displaystyle=-\nu\sum_{k=1}^{3}\mathbb{E}\left[\frac{\partial
F_{i}(W)}{\partial x^{k}}\frac{\partial W^{i}}{\partial x^{k}}\right]$
$\displaystyle=-\nu\sum_{k=1}^{3}\lim_{h\rightarrow
0}\frac{1}{h^{2}}\mathbb{E}\left[\left(F_{i}(W(x+he^{(k)},t))-F_{i}(W(x,t))\right)\left(W^{i}(x+he^{(k)},t)-W^{i}(x,t)\right)\right]$
$\displaystyle=-\nu\sum_{k=1}^{3}\lim_{h\rightarrow
0}\frac{1}{h^{2}}\iint_{\mathbb{R}^{3}\times\mathbb{R}^{3}}\left(F_{i}(w_{2})-F_{i}(w_{1})\right)\left(w_{2}^{i}-w_{1}^{i}\right)f_{2}(w_{1},w_{2},x,x+he^{(k)},t)\mathrm{d}w_{1}\mathrm{d}w_{2}$
$\displaystyle=\nu\int_{\mathbb{R}^{3}}F_{i}(w)\lim_{h\rightarrow
0}\frac{1}{h^{2}}\sum_{k=1}^{3}\left\\{f(w;x+he^{(k)},t)\int_{\mathbb{R}^{3}}\left(w_{1}^{i}-w^{i}\right)f_{2|1}(w_{1};w,x+he^{(k)},x,t)\mathrm{d}w_{1}\right.$
$\displaystyle+\left.f(w;x,t)\int_{\mathbb{R}^{3}}\left(w_{1}^{i}-w^{i}\right)f_{2|1}(w_{1};w,x,x+he^{(k)},t)\mathrm{d}w_{1}\right\\}\mathrm{d}w$
$\displaystyle=\nu\int_{\mathbb{R}^{3}}F_{i}(w)\lim_{h\rightarrow
0}\frac{1}{h^{2}}I_{4}^{h}dw,$

where for simplicity we have introduced the notation

$\displaystyle I_{4}^{h}:=$
$\displaystyle\sum_{k=1}^{3}\left\\{\left(f(w;x+he^{(k)},t)-f(w;x,t)\right)\mu^{i}(x+he^{(k)},x,w,t)\right.$
$\displaystyle+\left.f(w;x,t)\left(\mu^{i}(x,x+he^{(k)},w,t)+\mu^{i}(x+he^{(k)},x,w,t)\right)\right\\}.$

$I_{4}^{h}$ can be evaluated by using (5) when $h$ is sufficiently small, by
which we mean that

$\mu^{i}(x,x+he^{(k)},w,t)=a_{k}^{i}(x,w,t)h+b_{kk}^{i}(x,w,t)h^{2}+o\left(h^{2}\right)$

and

$\mu^{i}(x+he^{(k)},x,w,t)=-a_{k}^{i}(x+he^{(k)},w,t)h+b_{kk}^{i}(x+he^{(k)},w,t)h^{2}+o\left(h^{2}\right).$

Consequently, we have

$\displaystyle\mu^{i}(x,x+he^{(k)},w,t)+\mu^{i}(x+he^{(k)},x,w,t)$
$\displaystyle=-\left(a_{k}^{i}(x+he^{(k)},w,t)-a_{k}^{i}(x,w,t)\right)h$
$\displaystyle+\left(b_{kk}^{i}(x,w,t)+b_{kk}^{i}(x+he^{(k)},w,t)\right)h^{2}+o\left(h^{2}\right),$

which yields that

$\lim_{h\rightarrow 0}\frac{1}{h^{2}}I_{4}^{h}=-\frac{\partial}{\partial
x^{k}}\left(fa_{k}^{i}\right)+2f\sum_{k}b_{kk}^{i}.$

and we may conclude that

$\displaystyle I_{4}$
$\displaystyle=-\nu\int_{\mathbb{R}^{3}}F_{i}(w)\left[\frac{\partial}{\partial
x^{k}}\left(fa_{k}^{i}\right)-2f\sum_{k}b_{kk}^{i}\right]\mathrm{d}w$
$\displaystyle=\int_{\mathbb{R}^{3}}F(w)\frac{\partial}{\partial
w^{i}}\left[\nu\frac{\partial}{\partial x^{k}}\left(fa_{k}^{i}\right)-2\nu
f\sum_{k}b_{kk}^{i}\right]\mathrm{d}w.$ (16)

The PDF equation (8) now follows by substituting (12, 14, 15, 16) into (11). ∎

######  Remark 2.

The partial differential equation (8), although it appears linear, is in fact
highly non-linear, and more importantly, the coefficients $a,b,B$ and $D$ that
define the partial differential equation cannot be in general determined by
the PDF $f$ alone. We have seen that these functions are functionals of the
conditional average function $\mu$, hence the PDF $f$ cannot be obtained
solely through the PDF equation. However, the PDF PDE provides an appealing
and new method for modeling the PDF $f$ based on the modeling of $\mu$ alone,
we will provide in this work several results on modeling the PDF $f$ based on
the PDF equation we just derived.

######  Remark 3.

If all coefficients $a,b,B$ and $D$ are considered as given, then we can pose
the initial value problem for solving the PDF PDE (8). The good news is that
then the PDF PDE is a scalar linear partial differential equation of second
order, although in general it is a mixed type of partial differential
equations of second order on 6 dimensional space. To the best knowledge of the
authors of the present paper, this kind of partial differential equations has
not been studied systemically in the existing literature.

###  2.3 The distribution of turbulence

As an application of the PDF PDE obtained in the previous sub-section, we
address a long-standing question in turbulence about the distribution of
turbulent flows. It has been conjectured and verified by measurements over
many years (cf. [1] and [4]) that the distribution of a genuine turbulent flow
cannot be Gaussian, while a mathematical proof for this statement, to the best
knowledge of the present authors, is not available yet. By using the PDF PDE
we are able to prove the following

######  Theorem 4.

Consider an incompressible viscous (with viscosity $\nu>0$) turbulent flow
with vorticity $W(x,t)$. Suppose

1) the mean vorticity is constant at any instance,

2) $\\{W(x,t)\\}$ is weakly isotropic in the sense that

$\mathbb{E}\left[\nabla W(x,t)|W(x,t)=w\right]=0$

for all $x,w$ and $t\geq 0$, and

3) the distribution of $\\{W(x,t)\\}$ is Gaussian.

Then the distribution of $W(x,t)$ is independent of $x$.

###### Proof.

Since $\left\\{W(x,t):x\in\mathbb{R}^{d}\textrm{ and }t\geq 0\right\\}$ is a
Gaussian random field, for $x\neq y$, the joint law of $W(x,t)$ and $W(y,t)$
is a normal distribution on $\mathbb{R}^{6}$ whose co-variance matrix
$\Sigma(x,y)$ may be decomposed into blocks

$\Theta(x,y)=\left(\begin{array}[]{cc}\Sigma_{x}&\Sigma_{x,y}\\\
\Sigma_{y,x}&\Sigma_{y}\end{array}\right),$

where the dependence on $t$ is suppressed for simplicity. The PDF of $W(x,t)$
is given by

$f(w;x,t)=\frac{1}{(2\pi)^{3/2}\sqrt{\det\Sigma_{x}}}\exp\left(-\frac{1}{2}(w-m)^{T}\Sigma_{x}^{-1}(w-m)\right),$

where $m$ denotes the mean vorticity, so that

$\ln
f(w;x,t)=-\frac{3}{2}\ln(2\pi)-\frac{1}{2}\ln\det\Sigma_{x}-\frac{1}{2}(w-m)^{T}\Sigma_{x}^{-1}(w-m).$
(17)

Since under the assumption $a_{k}^{i}=0$, $\psi=\ln f$ according to PDF PDE
(8) must satisfy the following equation:

$\displaystyle\frac{\partial\psi}{\partial t}$
$\displaystyle=\nu\Delta_{x}\psi-B^{i}\frac{\partial\psi}{\partial
x^{i}}+\nu|\nabla_{x}\psi|^{2}$ $\displaystyle+\left(D^{i}-2\nu
b^{i}\right)\frac{\partial\psi}{\partial w^{i}}+\frac{\partial\left(D^{i}-2\nu
b^{i}\right)}{\partial w^{i}}-\nabla_{x}\cdot B$

(which is the PDF PDE for $\psi=\ln f$ in the case where $a_{k}^{i}$ vanish).
Using (17) and by a lengthy but completely elementary computation, we deduce
that

$\displaystyle-2\frac{\partial\psi}{\partial t}$
$\displaystyle=\nu\Delta_{x}\ln\det\Sigma_{x}-\frac{1}{2}\nu\left|\nabla_{x}\ln\det\Sigma_{x}\right|^{2}+4\nu\frac{\partial
b^{i}}{\partial w^{i}}$
$\displaystyle-B^{i}\frac{\partial\ln\det\Sigma_{x}}{\partial
x^{i}}-2\frac{\partial D^{i}}{\partial w^{i}}+2\nabla_{x}\cdot B+\nu
u^{T}\Delta_{x}\Sigma_{x}^{-1}u$ $\displaystyle-\nu
u^{T}\left(\sum_{i}\frac{\partial\ln\det\Sigma_{x}}{\partial
x^{i}}\frac{\partial\Sigma_{x}^{-1}}{\partial x^{i}}\right)u-4\nu
b^{i}(\Sigma_{x}^{-1})_{i\beta}u^{\beta}$
$\displaystyle+2D^{i}(\Sigma_{x}^{-1})_{i\beta}u^{\beta}-B^{i}u^{T}\frac{\partial\Sigma_{x}^{-1}}{\partial
x^{i}}u-\frac{1}{2}\nu\left|u^{T}\nabla_{x}\Sigma_{x}^{-1}u\right|^{2},$

where $u=w-m$. The last equation looks complicated, but the important point in
our argument is the observation that, the left-hand side
$-2\frac{\partial\psi}{\partial t}$ is a quadratic function of $u$, while the
right-hand side is a polynomial in $u$ of degree $4$, and only the last term
on the right-hand side has order 4 in $u$. Therefore

$\frac{1}{2}\nu\left|u^{T}\nabla_{x}\Sigma_{x}^{-1}u\right|^{2}=0$

for every $u$, which yields that $\nabla_{x}\Sigma_{x}^{-1}=0$, so that
$\Sigma_{x}$ is independent of $x$. This completes the proof. ∎

##  3 Inviscid fluid flows

There is a significant simplification in PDF PDE for an inviscid “turbulent”
flow, although such turbulence may not exist in nature. For inviscid fluid
flows, the velocity $U(x,t)$ satisfies the Euler equations

$\frac{\partial}{\partial t}U^{i}+U^{j}\frac{\partial U^{i}}{\partial
x^{j}}=-\nabla P,\quad\frac{\partial U^{i}}{\partial x^{i}}=0.$

The PDF $f(w;x,t)$ of $W=\nabla\wedge U$ satisfies the (non-linear) transport
differential equation

$\frac{\partial f}{\partial t}+\frac{\partial}{\partial
x^{i}}\left(fB^{i}\right)=\frac{\partial}{\partial w^{i}}\left(fD^{i}\right),$
(18)

where $B$ and $D$ are given as in (9) and (10) respectively.

The following theorem provides a mathematical tool for modeling PDF of an
inviscid turbulent flow.

######  Theorem 5.

Consider PDE (18) where the data $B$ and $D$ are assumed as given. Assume that
$B$ and $D$ are Lipschitz continuous in $(x,w)$ uniformly in $t>0$. Suppose
$f(w;x,t)$ is a solution to (18) with continuous initial data $f(w;x,0)$, and
suppose $fD$ decays to zero sufficiently fast as $|w|\to\infty$. Then we have
the followings:

1) $f(w;x,0)\geq 0$ for all $x\in\mathbb{R}^{3}$, then $f(w;x;t)\geq 0$ for
all $x\in\mathbb{R}^{3}$ and $t>0$.

2) Suppose $\int_{\mathbb{R}^{3}}f(w;x,0)\mathrm{d}w=1$ for all
$x\in\mathbb{R}^{3}$, then $\int_{\mathbb{R}^{3}}f(w;x,t)\mathrm{d}w=1$ for
all $x\in\mathbb{R}^{3}$ and $t>0$, if and only if

$\frac{\partial}{\partial
x^{i}}\int_{\mathbb{R}^{3}}B^{i}(x,w,t)f(w;x,t)\mathrm{d}w=0$ (19)

for all $x\in\mathbb{R}^{3}$ and $t>0$.

###### Proof.

Let us define the integral curves $X$ and $Y$ by the following ordinary
differential equation system:

$\begin{cases}\frac{\mathrm{d}}{\mathrm{d}s}X^{i}(t,s)=-B^{i}(X(t,s),Y(t,s),t-s)),&X(t,0)=x,\\\
\frac{\mathrm{d}}{\mathrm{d}s}Y^{i}(t,s)=D^{i}(X(t,s),Y(t,s),t-s)),&Y(t,0)=w.\end{cases}$

Define

$h(s)=f(Y(t,s);X(t,s),t-s)e^{N(w;x,s)}$

for all $s\in[0,t]$, where

$N(w;x,s):=\int_{0}^{s}\left(\frac{\partial D^{i}}{\partial
w^{i}}-\frac{\partial B^{i}}{\partial
x^{i}}\right)(X(t,r);Y(t,r),t-r)\mathrm{d}r.$

Then clearly $h(0)=f(w;x,t)$ and

$h(t)=f(Y(t,t);X(t,t),0)e^{N(w;x,t)}.$

It is clear that

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}s}h(s)=$ $\displaystyle
e^{N}\frac{\partial f}{\partial
w^{i}}\frac{\mathrm{d}}{\mathrm{d}s}Y^{i}(t,s)+e^{N}\frac{\partial f}{\partial
x^{i}}\frac{\mathrm{d}}{\mathrm{d}s}X^{i}(t,s)$
$\displaystyle-e^{N}\frac{\partial}{\partial t}f+e^{N}f\frac{\partial
D^{i}}{\partial w^{i}}-e^{N}f\frac{\partial B^{i}}{\partial x^{i}}$
$\displaystyle=$ $\displaystyle e^{N}\left\\{\frac{\partial(fD^{i})}{\partial
w^{i}}-\frac{\partial(B^{i}f)}{\partial x^{i}}-\frac{\partial}{\partial
t}f\right\\}$ $\displaystyle=$ $\displaystyle 0.$

Integrating with respect to the variable $s$ over $[0,t]$, we may conclude
that

$f(w;x,t)=f(Y(t,t);X(t,t),0)\exp\left[\int_{0}^{t}\left(\frac{\partial
D^{i}}{\partial w^{i}}-\frac{\partial B^{i}}{\partial
x^{i}}\right)(X(t,s);Y(t,s),t-s)\mathrm{d}s\right].$

The conservation of positivity property 1) follows immediately from this
representation. To show 2) we consider the function
$u(x,t)=\int_{\mathbb{R}^{3}}f(w;x,t)\mathrm{d}w$ for $x\in\mathbb{R}^{3}$ and
$t\geq 0$. If $u(x,t)=1$ for all $x$ and $t>0$, then by integrating (18) over
$\mathbb{R}^{3}$ we obtain (19). Conversely, if (19) holds, by integrating
(18) over $\mathbb{R}^{3}$ then $\frac{\partial}{\partial t}u(x,t)=0,$which
yields that $u(x,t)=1$. ∎

##  4 Weakly homogeneous and weakly isotropic flows

It has been pointed out in Secs. 2 and 3, the coefficients appearing in the
PDF PDE are determined by the conditional mean function $\mu(x,y,w,t)$ defined
in (1). The significance of this statistical quantity of a turbulent flow has
been demonstrated in the derivation of the PDF PDE. By definition the
conditional mean describes the deviation from the local isotropicity of the
turbulent flow, and therefore this statistical characteristic can be used in
the classification of turbulent flows. In this section, we define weak
homogeneous and weak isotropic turbulent flows, then study the PDF equation
for these turbulent flows.

The homogeneity and the isotropicity can be defined in general for random
fields indexed by a space variable $x\in\mathbb{R}^{d}$, which has been
introduced into the study of turbulence by G. I. Taylor [19]. The local
homogeneous and local isotropic flows were introduced by Kolmogorov for
formulating K41 theory (and its improved version K62 theory). According to
Kolmogorov [9, 8], a random field $\\{Z(x,t)\\}$ is locally homogeneous if the
conditional distribution of $Z(y,t)-Z(x,s)$ given $Z(x,s)=z$ is independent of
$(z,x,s)$ for any $t\geq s$, and further it is locally isotropic if the
conditional distribution of $Z(y,t)-Z(x,s)$ are invariant under reflections
and rotations. The usage of conditional mean function makes it possible to
generalize these concepts to their weak versions.

We may simplify the PDF PDE for turbulent flows when the vorticity random
field $W(x,t)$ is _weakly homogeneous_ in the sense that for every pair
$(y,x)$ and $t>0$, the conditional mean of $W^{i}(y,t)-W^{i}(x,t)$ given
$W(x,t)=w$ depends only on $y-x$, $w$ and $t>0$ but independent of $x$, so
that we may write

$\mu^{i}(x,y,w,t)=\beta^{i}(y-x,w,t),$

where $\beta^{i}(0,w,t)=0$. Then

$a_{k}^{i}(x,w,t)=\frac{\partial\beta^{i}}{\partial x^{k}}(0,w,t),$

and

$b^{i}(x,w,t)=\Delta\beta^{i}(0,w,t),$

which are still denoted by $a_{k}^{i}(w,t)$ and $b^{i}(w,t)$, though they are
independent of $x$. Similarly

$B^{i}(w,t)=\int_{\mathbb{R}^{3}}\frac{1}{4\pi|y|}\varepsilon^{ijk}\frac{\partial\beta^{k}}{\partial
y^{j}}(y,w,t)\mathrm{d}y,$

and

$D^{i}(w,t)=w^{l}\int_{\mathbb{R}^{3}}\frac{y^{l}}{4\pi|y|^{3}}\varepsilon^{ijk}\frac{\partial\beta^{k}}{\partial
y^{j}}(y,w,t)\mathrm{d}y,$

which again are functions of $(w,t)$ only. Therefore the PDF PDE can be
simplified to be the following mixed type PDE

$\left(\frac{\partial}{\partial t}+B^{i}\frac{\partial}{\partial
x^{i}}-\nu\Delta\right)f=\nu\frac{\partial}{\partial
w^{i}}\left(a_{k}^{i}\frac{\partial f}{\partial
x^{k}}\right)+\frac{\partial}{\partial w^{i}}\left(fA^{i}\right),$

where $A^{i}=D^{i}-2\nu b^{i}$ and $i=1,2,3$.

We say that the vorticity $W$ is _weakly isotropic_ , if

$\mathbb{E}\left[\nabla W(x,t)|W(x,t)=w\right]=0$

for all $x,w$ and $t$. By definition, $W$ is locally isotropic in Kolmogorov’s
sense, then it is weakly isotropic.

If $W$ is weakly homogeneous and weakly isotropic, then $a=0$, $B$ and $A$
depend on $(w,t)$ only, so for this case, the PDF PDE becomes a parabolic-
transport equation

$\left(\frac{\partial}{\partial t}+B^{i}\frac{\partial}{\partial
x^{i}}-\nu\Delta\right)f=\frac{\partial}{\partial w^{i}}\left(fA^{i}\right).$
(20)

After we have deduced the PDE (20), there is no need to assume that $B$ and
$A$ are independent of $x$. Therefore we may consider the following PDE

$\left(\frac{\partial}{\partial t}+B^{i}\frac{\partial}{\partial
x^{i}}+\frac{\partial B^{i}}{\partial
x^{i}}-\nu\Delta\right)f=\frac{\partial}{\partial w^{i}}\left(fA^{i}\right).$
(21)

where both $B$ and $A$ are functions of three variables $x$, $w$ and $t$.

The following theorem provides the foundation for modeling weakly isotropic
turbulent flows based on the vorticity PDF.

######  Theorem 6.

Consider the parabolic-transport differential equation (21), where $A$ and $B$
are considered as given. Assume that $A(x,w,t)$ and $B(x,w,t)$ are Lipschitz
continuous in space variable $(x,w)$ uniformly in $t$. Suppose $f(w;x,t)$ is a
solution to (21) such that $fA$ is continuous and decays to zero sufficiently
fast as $|w|\to\infty$.

1) If $f(w;x,0)\geq 0$ for all $x$ and $w$, then $f(w;x,t)\geq 0$ for all
$t\geq 0$, $x$ and $w$.

2) If $f(w;x,0)$ is a PDF for every $x,$ i.e. it is non-negative and
$\int_{\mathbb{R}^{3}}f(w;x,0)\mathrm{d}w=1$ for every $x$, then so are
$f(w;x,t)$ for all $x$ and $t>0$ if and only if

$\frac{\partial}{\partial
x^{i}}\int_{\mathbb{R}^{3}}B^{i}(x,w,t)f(w;x,t)\mathrm{d}w=0$ (22)

for all $x$ and $t>0$.

3) Suppose in addition $B$ is independent of $w$ and $\frac{\partial
B^{i}}{\partial x^{i}}=0$, then if $f(w;x,0)$ is a PDF for every $x$,
$f(w;x,t)$ is also a PDF for every $x$ and for every $t>0$.

###### Proof.

The proof is rather similar to that of Theorem 5, so we only provide an
outline. Let $(M_{t}^{i})_{t\geq 0}$ ($i=1,2,3$) be a three dimensional
standard Brownian motion on a probability space
$(\varOmega,\mathcal{F},\mathbb{P})$. For every fixed $t>0$ we run the
following SDE system:

$\mathrm{d}X^{i}(t,s)=\sqrt{2\nu}\mathrm{d}M_{s}^{i}-B^{i}(X(t,s),Y(t,s),t-s)\mathrm{d}s,\quad
X(t,0)=x,$ $\mathrm{d}Y^{i}(t,s)=A^{i}(X(t,s),Y(t,s),t-s)\mathrm{d}s,\quad
Y(t,0)=w,$

which has a unique solution running up to $t$. Apply Itô’s formula to

$H_{s}=f(Y(t,s);X(t,s),t-s)e^{P(s)},$

where

$P(s)=\int_{0}^{s}\frac{\partial A^{i}}{\partial
w^{i}}(X(t,r),Y(t,r),t-r)\mathrm{d}r.$

Then

$\displaystyle\mathrm{d}H_{s}$ $\displaystyle=e^{P}\left(\frac{\partial
f}{\partial w^{i}}\mathrm{d}Y_{s}^{i}+\frac{\partial f}{\partial
x^{i}}\mathrm{d}X_{s}^{i}+f\frac{\partial A^{i}}{\partial
w^{i}}-\frac{\partial f}{\partial s}\mathrm{d}s+\nu\Delta f\mathrm{d}s\right)$
$\displaystyle=\sqrt{2\nu}e^{P}\frac{\partial f}{\partial
x^{i}}\mathrm{d}M_{s}^{i}$

so that

$H_{0}=H_{t}-\sqrt{2\nu}e^{P}\frac{\partial f}{\partial
x^{i}}\mathrm{d}M_{s}^{i}.$

By taking expectation, one may deduce that the solution can be expressed as

$f(w;x,t)=\mathbb{E}\left[f(Y(t,t);X(t,t),0)\exp\left(\int_{0}^{t}\frac{\partial
A^{i}}{\partial w^{i}}(X(t,s),Y(t,s),t-s)\mathrm{d}s\right)\right].$

Statements 1) - 3) can be easily shown using this expression and the argument
in the proof of Theorem 5. If $B$ is independent of $w$, then, by integrating
(20),

$\left(\frac{\partial}{\partial t}+B^{i}\frac{\partial}{\partial
x^{i}}-\nu\Delta\right)u(x,t)=0,$

where $u(x,t)=\int_{\mathbb{R}^{3}}f(w;x,t)\mathrm{d}w$. Thus 3) follows from
the uniqueness of the solution of the previous parabolic equation. ∎

##  5 Modeling PDF of weakly isotropic flows

In this and the next sections, several simple models based on our PDF equation
for modeling vorticity distributions in turbulence are discussed. Clearly, the
most straightforward way for modeling the distribution of vorticity is to
assign the conditional mean function $\mu(x,y,w,t)$. The other parameters in
the PDF PDE may be determined accordingly.

For practical reasons which will be clarified in our computations below, the
simplest yet not trivial model for $\mu(x,y,w,t)$ should be a function of
$|y-x|^{2}$ only and $\mu$ decays sufficiently fast at infinity. For such a
model, the parameters $B^{i}$ and $D^{i}$ appearing in the PDF PDE (8) vanish
identically. The other parameters $a_{k}^{i}$ and $b_{jk}^{i}$ are determined
by the asymptotic condition as $|y-x|\rightarrow 0$. For this reason, we
assign

$\mu^{i}(x,y,w,t)=\begin{cases}-\frac{C}{6\nu}|y-x|^{2},&\textrm{ if
}|x-y|<\delta,\\\ h(|y-x|^{2}),&\textrm{ if }|x-y|\geq\delta,\end{cases}$ (23)

for $i=1,2,3$, where $\delta>0$ and $h$ is a $\mathcal{C}^{1}$-function such
that $h$ and its derivative decay sufficiently fast at infinity. $C\geq 0$ is
a model parameter which should be reflecting certain physical properties of
the underlying flows. The other parameters can be read out from this model
easily: $a_{k}^{i}=0$, $b_{jk}^{i}=-\delta_{jk}\frac{C}{6\nu}$ and therefore
$b^{i}=-\frac{C}{2\nu}$. With this model of $\mu$, the PDF PDE is simplified
to be the following parabolic-transport equation:

$\left(\frac{\partial}{\partial t}-\nu\Delta\right)f=C\sum_{i}\frac{\partial
f}{\partial w^{i}},$ (24)

The solution to (24) has a nice probabilistic representation which may be read
out from the proof of Theorem 6. In fact

$f(w;x,t)=\mathbb{E}\left[f(Y(t,t);X(t,t),0)\right]$ (25)

where $f(w;x,0)$ is the PDF of the initial vorticity, and $X$ and $Y$ are
solutions to the SDE:

$\begin{cases}\mathrm{d}X^{i}(t,s)=\sqrt{2\nu}\mathrm{d}M_{s}^{i},&X(t,0)=x,\\\
\mathrm{d}Y^{i}(t,s)=C\mathrm{d}s,&Y(t,0)=w,\end{cases}$

where $M$ is a standard 3D Brownian motion. The solutions when $s=t$ are given
by

$Y^{i}(t,t)=w^{i}+Ct$

and

$X(t,t)=x+\sqrt{2\nu}M_{t}.$

Substituting these into (25) we obtain

$\displaystyle f(w;x,t)$
$\displaystyle=\mathbb{E}\left[f(w+Ct;x+\sqrt{2\nu}M_{t},0)\right]$
$\displaystyle=\mathbb{E}\left[f(w+Ct;x+\sqrt{2\nu t}\xi,0)\right],$ (26)

where $\xi$ is a random vector with the standard 3D normal distribution
$N(0,I_{3})$. Clearly the representation (26) may also be written in terms of
Gaussian density

$\displaystyle f(w;x,t)$ $\displaystyle=\frac{1}{(4\pi\nu
t)^{3/2}}\int_{\mathbb{R}^{3}}f\left(w+Ct,z,0\right)e^{-\frac{|z-x|^{2}}{4\nu
t}}\mathrm{d}z$
$\displaystyle=\frac{1}{(2\pi)^{3/2}}\int_{\mathbb{R}^{3}}f\left(w+Ct,x+\sqrt{2\nu
t}z,0\right)e^{-\frac{|z|^{2}}{2}}\mathrm{d}z.$ (27)

Both representations (26) and (27) may be used to evaluate $f(w;x,t)$ by using
for example Monte-Carlo scheme by sampling Gaussian random variables.

The initial PDF $f(w;x,0)$ is determined by the initial vorticity
distribution. For a turbulent flow, the initial vorticity may be written as a
sum:

$W(x,0)=\omega_{0}(x)+\varepsilon(x),$

where $\omega_{0}(x)$ is the initial mean vorticity at the location $x$, and
$\varepsilon(x)$ represents a small random perturbation. It is reasonable,
therefore, to assume that $\varepsilon(x)$ has a normal distribution
$N(0,\sigma(x)^{2}I_{3})$, where the variance $\sigma(x)$ may or may not
depend on the location $x$. However, if $\varepsilon(x)$ is allowed to depend
on the location $x$, then the noise $\varepsilon(x)$ has to satisfy the
divergence-free condition as well, a technical issue we will not address here
in detail. It is reasonable to assume that the initial vorticity mean
$\omega_{0}$ is distributed in a small region of the space. In particular
$\omega_{0}$ decays to zero sufficiently fast near infinity. Under this
assumption,

$f(w;x,0)=\frac{1}{\left(2\pi\sigma(x)^{2}\right)^{3/2}}\exp\left[-\frac{|w-\omega_{0}(x)|^{2}}{2\sigma(x)^{2}}\right].$
(28)

By substituting this into the representation (26), we obtain

$f(w;x,t)=\mathbb{E}\left[\frac{1}{\left(2\pi\sigma\left(x+\sqrt{2\nu
t}\xi\right)^{2}\right)^{3/2}}e^{-\frac{\left|w+Ct-\omega_{0}(x+\sqrt{2\nu
t}\xi)\right|^{2}}{2\sigma\left(x+\sqrt{2\nu t}\xi\right)^{2}}}\right].$ (29)

It is easy to see that $f(w;x,t)$ is no longer Gaussian except for special
cases which we would like to discuss below.

If $\omega_{0}$ vanishes and if $\sigma$ is a positive small constant, the
initial distribution is homogeneous and equation (29) leads to a simple
expression

$f(w;x,t)=\frac{1}{\left(2\pi\sigma^{2}\right)^{3/2}}\exp\left(-\frac{\left|w+Ct\right|^{2}}{2\sigma^{2}}\right),$

which is independent of the location $x$ and remains a Gaussian density. The
interesting feature about this model is that the variance stays as the
constant $\sigma^{2}$ but new mean vorticity $-Ct$ is created evenly after
duration $t$. This is the case of a turbulence with a small constant random
perturbation. For turbulent flows observed in nature, the initial vorticity
mean $\omega_{0}$ does exist and does not vanish, and the random noise
$\varepsilon(x)$, for simplicity, may be modeled by a Gaussian random variable
independent of $x$. Therefore when $\sigma>0$ is a small constant, we have

$f(w;x,t)=\frac{1}{\left(2\pi\sigma^{2}\right)^{3/2}}\mathbb{E}\left[e^{-\frac{|w+Ct-\omega_{0}\left(x+\sqrt{2\nu
t}\xi\right)|^{2}}{2\sigma^{2}}}\right],$ (30)

where $\xi\sim N(0,I_{3})$.

For this case the mean vorticity $\omega(x,t)$ at $(x,t)$ can be evaluated.
Indeed

$\displaystyle\omega(x,t)$
$\displaystyle=\int_{\mathbb{R}^{3}}\frac{w}{\left(2\pi\sigma^{2}\right)^{3/2}}\mathbb{E}\left[e^{-\frac{|w+Ct-\omega_{0}\left(x+\sqrt{2\nu
t}\xi\right)|^{2}}{2\sigma^{2}}}\right]\mathrm{d}w$
$\displaystyle=\mathbb{E}\left[\int_{\mathbb{R}^{3}}\frac{w}{\left(2\pi\sigma^{2}\right)^{3/2}}e^{-\frac{|w+Ct-\omega_{0}\left(x+\sqrt{2\nu
t}\xi\right)|^{2}}{2\sigma^{2}}}\mathrm{d}w\right]$
$\displaystyle=\mathbb{E}\left[\omega_{0}\left(x+\sqrt{2\nu
t}\xi\right)\right]-Ct,$

and therefore

$\omega^{i}(x,t)=\int_{\mathbb{R}^{3}}\frac{\omega_{0}^{i}(y)}{\left(4\pi\nu
t\right)^{3/2}}e^{-\frac{\left|y-x\right|^{2}}{4\nu t}}dy-Ct.$ (31)

This equality shows that the vorticity mean $\omega(x,t)$ under this simple
model is independent of the noise parameter $\sigma^{2}$, and $\omega(x,t)$
evolves according to the heat type equation

$\frac{\partial}{\partial t}\omega(x,t)-C=\nu\Delta\omega(x,t),\textrm{
}\omega(x,0)=\omega_{0}(x),$

which is a rather crude approximation to the mean vorticity equation.

##  6 Heat flow method

In this section we propose another model for the PDF of the vorticity, based
on the heat flow method, in which the conditional mean function $\mu$ is
generated by a random field $R(x,\tau)$. The random field $R(x,\tau)$ evolves
according to the heat flow

$\left(\frac{\partial}{\partial\tau}-\nu\Delta\right)R^{i}(x,\tau)=0,$ (32)

where $R(x,\tau)=(R^{1}(x,\tau),R^{2}(x,\tau),R^{3}(x,\tau))$. The initial
value $R(x,0)=\xi(x)=(\xi^{i}(x))$ is a centered Gaussian noise white in space
in the sense that the co-variance of $\xi^{i}$ at two locations $u$ and $v$ in
$\mathbb{R}^{3}$ is given by

$\mathbb{E}\left[\xi^{i}(u)\xi^{i}(v)\right]=\delta(u-v).$ (33)

Also, we assume that $\xi^{i}$’s are i.i.d. random variables. For simplicity,
we assume that $R^{i}$ is centered here, but our argument can definitely be
generalized to the case when it is non-centered.

The solution to equation (32) is given by

$R^{i}(x,\tau)=\int_{\mathbb{R}^{3}}\frac{1}{(4\pi\nu\tau)^{\frac{3}{2}}}e^{-\frac{|x-y|^{2}}{4\nu\tau}}\xi^{i}(y)\mathrm{d}y,$
(34)

so that $R(x,\tau)$ is a centered Gaussian random field with its co-variance

$\displaystyle\sigma_{\tau}(x,y)=$
$\displaystyle\int_{\mathbb{R}^{3}}\int_{\mathbb{R}^{3}}\frac{1}{(4\pi\nu\tau)^{3}}e^{-\frac{|x-u|^{2}+|y-v|^{2}}{4\nu\tau}}\mathbb{E}\left[\xi^{i}(u)\xi^{i}(v)\right]\mathrm{d}u\mathrm{d}v$
$\displaystyle=$
$\displaystyle\int_{\mathbb{R}^{3}}\int_{\mathbb{R}^{3}}\frac{1}{(4\pi\nu\tau)^{3}}e^{-\frac{|x-u|^{2}+|y-v|^{2}}{4\nu\tau}}\delta(u-v)\mathrm{d}u\mathrm{d}v,$
$\displaystyle=$ $\displaystyle
e^{-\frac{|x-y|^{2}}{8\nu\tau}}\int_{\mathbb{R}^{3}}\frac{1}{(4\pi\nu\tau)^{3}}e^{-\frac{|u-\frac{(x+y)}{2}|^{2}}{2\nu\tau}}\mathrm{d}u$
$\displaystyle=$
$\displaystyle\frac{1}{8(2\pi\nu\tau)^{\frac{3}{2}}}e^{-\frac{|x-y|^{2}}{8\nu\tau}}.$

It follows that the conditional distribution of $R(y,\tau)$ given
$R(x,\tau)=w$ has a normal distribution with mean
$e^{-\frac{|x-y|^{2}}{8\nu\tau}}w$ and co-variance matrix

$\frac{1}{8(2\pi\nu\tau)^{\frac{3}{2}}}\left(1-e^{-\frac{|x-y|^{2}}{4\nu\tau}}\right)I_{3}.$

Therefore the conditional mean of $R(y,\tau)-R(x,\tau)$ given $R(x,\tau)=w$
can be easily found to be

$-\left(1-e^{-\frac{|x-y|^{2}}{8\nu\tau}}\right)w,$

which will be our $\mu(x,y,w,t)$ with $\tau=\varphi(t)$ a reparametrization as
part of the model. For simplicity let us consider the power law model, that
is,

$\tau=\lambda_{1}(t+\lambda_{2})^{\alpha},$ (35)

where $\lambda_{1}$ and $\lambda_{2}$ are two positive numbers which then
become our model parameters. Since

$\mu(x,y,w,t)=-\frac{|x-y|^{2}}{8\nu\tau}w+o(|x-y|^{2})w,$

therefore $a_{k}^{i}=0$ and $b_{jk}^{i}=-\delta_{jk}\frac{1}{8\nu\tau}w^{i}$.
In particular,

$b^{i}(w,t)=-\frac{3w^{i}}{8\nu\lambda_{1}(t+\lambda_{2})^{\alpha}}$ (36)

for $i=1,2,3$. Since the conditional mean function $\mu$ depends only on
$|x-y|^{2}$ and decays exponentially fast at infinity, $B^{i}=0$ and
$D^{i}=0$, which implies that

$A^{i}(w,t)=\frac{3w^{i}}{4\lambda_{1}(t+\lambda_{2})^{\alpha}}.$ (37)

Thus the PDF PDE with these parameters is reduced to the simple parabolic-
transport equation

$\left(\frac{\partial}{\partial
t}-\nu\Delta_{x}\right)f=\frac{\partial}{\partial w^{i}}\left(A^{i}f\right).$

The divergence of $A$ is given by

$\frac{\partial A^{i}}{\partial
w^{i}}(x,w,t)=\frac{9}{4\lambda_{1}(t+\lambda_{2})^{\alpha}}.$

If $\alpha\neq 1$ but $\alpha>0$, by using the Feynman-Kac formula we have

$f(w;x,t)=\theta(\alpha,t)^{3}\mathbb{E}\left[f(Y(t,t),X(t,t),0)\right],$ (38)

where we have introduced

$\theta(\alpha,t)=\exp\left[\frac{3\left((t+\lambda_{2})^{1-\alpha}-\lambda_{2}^{1-\alpha}\right)}{4\lambda_{1}(1-\alpha)}\right]$
(39)

for simplicity, and $X$ and $Y$ are solutions to SDE

$\begin{cases}\mathrm{d}X^{i}(t,s)=\sqrt{2\nu}\mathrm{d}M_{s}^{i},&X(t,0)=x,\\\
\mathrm{d}Y^{i}(t,s)=\frac{3Y^{i}(t,s)}{4\lambda_{1}(t-s+\lambda_{2})^{\alpha}}\mathrm{d}s,&Y(t,0)=w.\end{cases}$

If $\alpha=1$ then

$f(w;x,t)=\lambda(t)^{3}\mathbb{E}\left[f(Y(t,t),X(t,t),0)\right],$ (40)

where

$\lambda(t)=\left(\frac{t+\lambda_{2}}{\lambda_{2}}\right)^{\frac{3}{4\lambda_{1}}}.$
(41)

The previous SDEs have explicit solutions:

$X(t,t)=x+\sqrt{2\nu}M_{t}$

and

$Y(t,t)=w\theta(\alpha,t)$

if $\alpha\neq 1$, and

$Y(t,t)=\lambda(t)w$

if $\alpha=1$.

Therefore, if $\alpha>0$ and $\alpha\neq 1$, by plugging these explicit
solutions $X$ and $Y$ into (38) we have

$\displaystyle f(w;x,t)$
$\displaystyle=\theta(\alpha,t)^{3}\mathbb{E}\left[f\left(w\theta(\alpha,t),x+\sqrt{2\nu}M_{t},0\right)\right]$
$\displaystyle=\theta(\alpha,t)^{3}\mathbb{E}\left[f\left(w\theta(\alpha,t),x+\sqrt{2\nu
t}\xi,0\right)\right]$ (42)

where $\xi\sim N(0,I_{3})$.

As in the previous section, the initial vorticity is assumed to be of the form

$W(x,0)=\omega(x)+\varepsilon(x),$

where $\varepsilon(x)$ has a normal distribution $N(0,\sigma^{2})$ and
$\omega(x)$ is the mean vorticity at $x$, so that

$f(w;x,0)=\frac{1}{\left(2\pi\sigma^{2}\right)^{3/2}}\exp\left[-\frac{\left|w-\omega(x)\right|^{2}}{2\sigma^{2}}\right]$

and therefore, according to (42),

$f(w;x,t)=\theta(\alpha,t)^{3}\mathbb{E}\left[\frac{1}{\left(2\pi\sigma(x+\sqrt{2\nu
t}\xi)^{2}\right)^{3/2}}e^{-\frac{\left|w\theta(\alpha,t)-\omega(x+\sqrt{2\nu
t}\xi)\right|^{2}}{2\sigma(x+\sqrt{2\nu t}\xi)^{2}}}\right],$

where $\xi\sim N(0,I_{3})$. Hence $f(w;x,t)$ is no longer Gaussian density
unless $\sigma(x)$ is independent of $x$ and $\omega=0$ identically.

Let us discuss a special case where $\omega=0$ identically and $\sigma>0$ is
constant. Then

$f(w;x,t)=\frac{\theta(\alpha,t)^{3}}{\left(2\pi\sigma^{2}\right)^{3/2}}\exp\left[-\frac{\left|w\theta(\alpha,t)\right|^{2}}{2\sigma^{2}}\right]$

is a Gaussian density with mean $0$ and variance $\rho^{2}I_{3}$, where

$\rho^{2}(x,t)=\frac{1}{\sqrt{\theta(\alpha,t)}}\sigma^{2}.$ (43)

If $\sigma>0$ is a constant and $\omega$ does not vanish identically, then
$f(w;x,t)$ is not Gaussian, but has a nice representation:

$f(w;x,t)=\frac{\theta(\alpha,t)^{3}}{\left(2\pi\sigma^{2}\right)^{3/2}}\mathbb{E}\left[e^{-\frac{\left|w\theta(\alpha,t)-\omega(x+\sqrt{2\nu
t}\xi)\right|^{2}}{2\sigma^{2}}}\right],$ (44)

where $\xi\sim N(0,I_{3})$, and $\theta(\alpha,t)$ is defined by (39).

A similar discussion applies to the model where $\alpha=1$, which is certainly
interesting too. For this model the PDF $f(w;x,t)$ is given by

$f(w;x,t)=\lambda(t)^{3}\mathbb{E}\left[f\left(w\lambda(t),x+\sqrt{2\nu
t}\xi,0\right)\right]$ (45)

where $\xi\sim N(0,I_{3})$ and $\lambda(t)$ is defined by (41).

Suppose again the initial vorticity

$W(x,0)=\omega(x)+\varepsilon(x),$

where $\varepsilon(x)$ has a normal distribution $N(0,\sigma^{2}(x))$ and
$\omega(x)$ is the mean. Then

$f(w;x,t)=\mathbb{E}\left[\frac{\lambda(t)^{3}}{\left(2\pi\sigma(x+\sqrt{2\nu
t}\xi)^{2}\right)^{3/2}}e^{-\frac{\left|w\lambda(t)-\omega\left(x+\sqrt{2\nu
t}\xi\right)\right|^{2}}{2\sigma(x+\sqrt{2\nu t}\xi)^{2}}}\right],$ (46)

which is not Gaussian except for the following special case.

If $\sigma(x)=\sigma>0$ is a constant and $\omega(x)=0$, then

$f(w;x,t)=\frac{\lambda(t)^{3}}{\left(2\pi\sigma^{2}\right)^{3/2}}e^{-\frac{\left|w\lambda(t)\right|^{2}}{2\sigma^{2}}}$

is Gaussian with mean zero and variance $\rho^{2}I_{3}$, where

$\rho^{2}(x,t)=\sigma^{2}\left(\frac{\lambda_{2}}{t+\lambda_{2}}\right)^{\frac{3}{2\lambda_{1}}}.$

The most interesting case for the purpose of modeling turbulent flows is the
model where $\sigma>0$ is a small parameter and $\omega$ is not zero, so that

$f(w;x,t)=\frac{\lambda(t)^{3}}{(2\pi\sigma^{2})^{3/2}}\mathbb{E}\left[e^{-\frac{\left|w\lambda(t)-\omega\left(x+\sqrt{2\nu
t}\xi\right)\right|^{2}}{2\sigma^{2}}}\right]$

where $\xi\sim N(0,I_{3})$. Although it is no longer Gaussian, some of its
features can be extracted by doing Monte-Carlo simulations.

We may conclude that this model provides some nice features of the propagation
of the vorticity which may be helpful for the understanding of the energy
dissipation in turbulence, and we will explore this in a separate work.

## References

  * [1] G. K. Batchelor.  The Theory of Homogeneous Turbulence.  Cambridge University Press, 1953.
  * [2] A. J. Chorin.  Vorticity and Turbulence.  Springer Science & Business Media, 2013\.
  * [3] G.-H. Cottet and P. D. Koumoutsakos.  Vortex Methods: Theory and Practice.  Cambridge University Press, 2000.
  * [4] P. A. Davidson.  Turbulence: An Introduction for Scientists and Engineers.  Oxford University Press, 1995.
  * [5] U. Frisch.  Turbulence: The Legacy of A. N. Kolmogorov.  Cambridge University Press, 1995.
  * [6] E. Hopf.  Statistical hydromechanics and functional calculus.  Journal of Rational Mechanics and Analysis, $\mathbf{1}$:87–123, 1952.
  * [7] A. N. Kolmogorov.  A refinement of previous hypotheses concerning the local structure of turbulence in a viscous incompressible fluid at high Reynolds number.  Journal of Fluid Mechanics, $\mathbf{13}$(1):82–85, 1962.
  * [8] A. N. Kolmogorov.  The local structure of turbulence in incompressible viscous fluid for very large Reynolds numbers.  Proceedings of the Royal Society of London. Series A: Mathematical and Physical Sciences, $\mathbf{434}$(1890):9–13, 1991\.
  * [9] A. N. Kolmogorov.  Dissipation of energy in the locally isotropic turbulence.  Proceedings of the Royal Society of London. Series A: Mathematical and Physical Sciences, $\mathbf{434}$(1890):15–17, 1991.
  * [10] M. Lesieur.  Turbulence in Fluids.  Springer Science & Business Media, 1997.
  * [11] A. J. Majda and A. L. Bertozzi.  Vorticity and Incompressible Flow.  Cambridge University Press, 2002.
  * [12] H. K. Moffatt.  A brief introduction to vortex dynamics and turbulence.  In Environmental Hazards: The Fluid Dynamics and Geophysics of Extreme Events, pages 1–27. World Scientific, 2011.
  * [13] A. Monin and A. Yaglom.  Statistical Fluid Mechanics: Mechanics of Turbulence, volume 1 and 2.  MIT Press, 1971 and 1975.
  * [14] S. B. Pope.  PDF methods for turbulent reactive flows.  Progress in Energy and Combustion Science, $\mathbf{11}$(2):119–192, 1985.
  * [15] S. B. Pope.  Turbulent Flows.  Cambridge University Press, 2000.
  * [16] Y. A. Rozanov.  Markov Random Fields. Springer, 1982.
  * [17] P. G. Saffman.  Vortex Dynamics.  Cambridge University Press, 1992.
  * [18] G. I. Taylor.  Diffusion by continuous movements.  Proceedings of the London Mathematical Society, $\mathbf{s2}$-$\mathbf{20}$(1):196–212, 1921.
  * [19] G. I. Taylor.  Statistical theory of turbulence, Parts 1-5.  Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, $\mathbf{151}$(873):421–478, 1935\.