# Targeting high symmetry in structure predictions by biasing the potential
energy surface

Hannes Huber  Department of Physics, University of Basel, Klingelbergstrasse
82, CH-4056 Basel, Switzerland  Martin Sommer  Department of Physics,
University of Basel, Klingelbergstrasse 82, CH-4056 Basel, Switzerland  Moritz
Gubler  Department of Physics, University of Basel, Klingelbergstrasse 82,
CH-4056 Basel, Switzerland  Stefan Goedecker  Department of Physics,
University of Basel, Klingelbergstrasse 82, CH-4056 Basel, Switzerland

###### Abstract

Ground state structures found in nature are in many cases of high symmetry.
But structure prediction methods typically render only a small fraction of
high symmetry structures. Especially for large crystalline unit cells there
are many low energy defect structures. For this reason methods have been
developed where either preferentially high symmetry structures are used as
input or where the whole structural search is done within a certain symmetry
group. In both cases it is necessary to specify the correct symmetry group
beforehand. However it can in general not be predicted which symmetry group is
the correct one leading to the ground state. For this reason we introduce a
potential energy biasing scheme that favors symmetry and where it is not
necessary to specify any symmetry group beforehand. On this biased potential
energy surface, high symmetry structures will be found much faster than on an
unbiased surface and independently of the symmetry group to which they belong.
For our two test cases, a $C_{60}$ fullerene and bulk silicon carbide, we get
a speedups of 25 and 63. In our data we also find a clear correlation between
the similarity of the atomic environments and the energy. In low energy
structures all the atoms of a species tend to have similar environments.

Structure prediction methods are an important tool for the discovery of new
materials [1]. Such methods can not only be applied for materials at ambient
pressure but also under very high pressures that are relevant for geophysical
applications, but not accessible by experimental methods [2]. For this reason,
numerous methods such as simulated annealing [3], basin hopping [4], minima
hopping (MH) [5, 6, 7, 8], random structure searches [9], meta-dynamics [10]
and various variants of evolutionary genetic algorithms [11, 12, 13, 14, 15,
16] as implemented in the USPEX [17], CALYPSO [18] and XtalOpt [19] software
package have been developed. These advanced global geometry optimisation
methods have shown that they can efficiently explore[20, 21] the potential
energy surface (PES) of different clusters and bulk materials under a variety
of external conditions and predict new structures. However, such methods
require a high computational effort because the number of possible meta-stable
structures grows exponentially with respect to the number of atoms in the
system and the generation and relaxation of a single structure requires many
energy and force evaluations.

Unless defects of certain materials are studied explicitly, the ground state
and the lowest defect free meta-stable structures are of greatest interest
since they are the structures that can most likely be synthesized. However,
for large cells most structures that are found in a structure prediction
contain defects. Typically these defect structures represent local minima in a
funnel whose bottom corresponds to a defect-free, meta-stable or the global
minimum structure. To find this relatively small number of defect-free
structures, a structure search which visits a very large number of defect
structures, is inefficient.
To favor high symmetry, most crystal structure prediction methods use input
guess structures that are of high symmetry. If the correct symmetry is chosen,
the most similar low energy structure is found much more rapidly. For basin
hopping and genetic algorithms there exist also a versions where all the moves
of the atoms are constrained to conserve the desired symmetry [22, 4]. The
inconvenience in all these approaches is that there are more than 200 space
groups and it is a priori unknown which one will be adopted by the system.
Traditionally symmetry is defined by geometric operations such as rotations or
reflections that leave the structure invariant. We will use in this work an
alternative definition of symmetry. We will consider a system to be highly
symmetric if all atoms of the same element see only a small number of
different environments. Structures with a large number of environments are
actually unlikely to exist according Pauling’s rule of structural parsimony
originally established for ionic materials [23, 24]. In many cases we will
actually try to find systems where all atoms of the same element see the same
environment. Evidently this is true for many high symmetry structures such as
the $C_{60}$ fullerene or the diamond structure of silicon and carbon. A
structure is either invariant or not under certain symmetry operations. So
basing some penalty function on the number of possible symmetry operations
would give rise to a discontinuous function. Our definition of similarity is
however a continuous functions. It is zero if the environments are identical
and grows larger in a continuous way when the environments become more
different. Our definition is thus broader than the traditional one. We
classify a structure also as highly symmetric if there are a few distinct
environments which are however very similar.

The tendency that low energy structures have in general similar environments
has already been exploited to gain efficiency in the context of evolutionary
structure prediction algorithms. In this context mutation moves were
introduced that favor environments that have a similar radial distribution as
certain selected role model environments [25, 26]. Even in amorphous systems
it was observed that structures that had similar pair distribution functions
were also low in energy [27].

The basic idea of our approach is to perform a structure search on a biased
potential energy surface [28]

$E_{b}(R_{1},\ldots,R_{n})=E(R_{1},\ldots,R_{n})+\omega P(R_{1},\ldots,R_{n})$

where $E_{b}$ is the biased PES, $E$ is the physical PES, $P$ is the penalty
function and $\omega$ is the biasing weight. Since the number of environments
is larger for defective structures, the penalty part will push up these
defective structures on the biased PES. In this way the downhill barriers are
lowered compared to the uphill barriers and the PES becomes a stronger
structure seeker character which speeds up the search for the global minimum
and possibly other high symmetry structures at the bottom of other funnels.

We quantify the similarity of environments with the overlap matrix (OM)
fingerprint [29] based on s- and p-type orbitals, that was shown to be able to
detect in a highly reliable way different atomic environments [30]. In
particular this fingerprint has both radial and angular resolution. In the OM
method the eigenvalues of a localized overlap matrix, centered on the atom $i$
whose environment has to be characterized, are assembled into an atomic
environment fingerprint vector $\bf{f}_{i}$. This environment characterization
is done for all atoms in the system. If all atomic environments are identical,
the rank of the matrix $F$ formed by all these vectors $\bf{f}_{i}$ is one, if
there are two distinct elemental environments the rank is two, etc. The rank
can most easily be calculated from the eigenvalues $\lambda_{i}$ of the Gram
matrix $D=F^{T}F$, constructed from the these fingerprint vectors. The number
of the non-zero eigenvalues of this matrix gives the rank of the fingerprint
vectors. So the penalty function that favours one single environment for a
certain element is

$P_{1}(R_{1},\ldots,R_{N_{at}})=\sum_{i=2}^{N_{at}}\lambda_{i}=\operatorname{Tr}(D)-\lambda_{1},$
(1)

In case we want to allow for up to two environments, the penalty becomes

$P_{2}(R_{1},\ldots,R_{N_{at}})=\sum_{i=3}^{N_{at}}\lambda_{i}=\operatorname{Tr}(D)-\lambda_{1}-\lambda_{2},$
(2)

where $\operatorname{Tr}$ is the trace of the matrix, i.e. the sum over all
eigenvalues. As usual, we have assumed in all the above formulas that the
eigenvalues are sorted in decreasing order. For a multi-component system, each
element contributes its own penalty function and the total penalty function is
the sum of all the elemental contributions. For highly symmetric structures
where all local environments are equivalent, e.g. the ground state of C60, the
bias function will be exactly zero. If the environments get more distinct, the
bias function grows due to the positive semi-definiteness of the Gram matrix.
Since the Gram matrix gives essentially the effective dimension of the vector
space spanned by the local descriptor vectors, it is called dimensionality
matrix in this paper.

To test our method we selected two systems of quite different nature. The
first one, silicon carbide, is a crystalline system with two elements that
have to mix in the right way to find low energy structures and the second, the
C${60}$ fullerene, is a molecular clusters. Its global minimum is just one
structure out of a huge number of meta-stable structures with varying
structural motifs such as planar structures, chains and bowls.

For the exploration of the PES the minima hopping (MH) algorithm was used, but
our biasing scheme is in principle applicable to any structure prediction
method. The MH algorithm is not based on thermodynamic principles like
simulated annealing or basin hopping but uses a combination of molecular
dynamics, local geometry optimization and a history of previously found local
minima to escape quickly from already known regions and hence, efficiently
explore the entire PES. For the geometry optimization of C60 with free
boundary conditions the conjugate gradient method was used. In the case of
periodic boundary conditions (PBC) the highly efficient and stable vc-SQNM
method developed by Gubler et al. [31, 32] was used. The molecular dynamics
(MD) simulation was implemented using the velocity Verlet algorithm for the
non-periodic case and the variable cell shape MD [33] for PBC. This method
allows atoms as well as cell vectors to move dynamically during the MD
simulation for PBC.
For the calculation of the PES of C60 the transferable tight binding potential
for carbon from Xu et al. [34] was used. For the silicon carbide simulations
in PBC DFTB+ [35] was used with the Slater-Koster parameterisation set pbc-0-3
[36]. The MD and an initial local geometry optimisation are performed on the
biased PES followed by a local geometry optimisation on the unbiased PES. This
avoids falling in potentially existing spurious local minima on the biased
PES.

To obtain conservative forces of the biased PES the derivative of the symmetry
bias needs to be added to the physical forces. The same is true for the
derivative of the biased symmetry function with respect to the lattice vectors
which need to be added to the lattice derivatives in the case of PBC. The
derivations of these two quantities can be found in the supplementary
information.

Figure 1: Changes in the characteristics of the disconnectivity graphs of
silicon carbide (top row) and $C_{60}$ (bottom row) induced by a bias. The
left column shows the disconnectivity graphs of the PES without a bias and the
right column the disconnectivity graphs of the PES with a bias. The graph was
constructed with the disconnectionDPS software [37].

The character of a PES can best be deduced from the appearance of its
disconnectivity graph [38] Fig. 1. For a structure seeker [39], the downhill
barriers are much lower than the uphill barriers. As a consequence, any
algorithm that crosses preferentially lower barriers will experience some
driving force toward the minimum at the bottom of the funnel and find it
therefore faster. This driving force will of course depend on the strength
$\omega$ of the bias. While on the one hand it is desirable to choose a large
$\omega$, the penalty should on the other hand only induce some weak
perturbation that does not completely deform the physical PES. In particular
there should remain in most cases a one-to-one mapping between the local
minima on the physical and the biased PES. As already noted by Zwanzig in the
context of protein folding [40] a relatively small bias can have a large
effect on the dynamics of the system and reduce the folding time by several
orders of magnitude. We were indeed always able to find a range of values for
$\omega$ that did speed up the search for high symmetry structures
considerable without destroying the overall character of the PES. With our
weight the downhill barriers are typically twice as large as the uphill
barriers and the penalty difference between high and low symmetry structures
is a few times the difference of their physical energy. This later criterion
can be used to find suitable values of $\omega$.

Fig. 1 shows the differences of the disconnectivity graphs for the unbiased
and biased system. The changes in the appearance of the disconnectivity graphs
indicate that the biased PES has a much stronger structure seeker character
which should make the search for the lowest high symmetry structures
considerably faster.

Figure 2: Comparison of the distribution of the found minima on the unbiased
(a) and biased (b) PES. The symmetry, plotted along the x-axis, is measured by
the number symmetry operations that leave the structure invariant. The
coloring indicates the relative abundance of meta-stable structures for a
given symmetry and energy.

To investigate the effect of the symmetry bias on the speed of the global
geometry optimization in a systematic way, statistical tests were conducted
for $C_{60}$ and 16 atom silicon carbide cells. One hundred global geometry
optimizations were started from different initial configurations until the
ground state structure, or in the case of silicon carbide a polytype of the
ground state, was found on the unbiased and the biased PES. For $C_{60}$ the
carbon atoms were randomly placed on a plane and for silicon carbide the
carbon and silicon atoms were randomly placed in spatially separated sub-cells
that formed the crystalline cell. To avoid nonphysical structures a minimum
and maximum distance between the randomly placed atoms was enforced. Since
there are no phase separated low energy structures in a cell of this size, the
silicon and carbon atoms always had to mix to find the low energy structures.
As a measure for the computational cost of the runs we used the number of
required local geometry optimizations. As can be seen from Table 1 the biasing
reduces the average number of geometry optimizations by a factor of 25 for C60
and by a factor of 63 for silicon carbide. It can also be seen from Table 1
that other statistical markers like the standard deviation (std), quantiles
and the number of geometry optimizations for the fastest as well as the
slowest simulations decreased by about the same magnitude.

Table 1: Table of statistical markers and biasing parameters of 100 global geometry optimisations started from randomly generated structures for C60 and 16 atom silicon carbide cells on the unbiased PES as well as on the biased PES. The statistical markers always relate to the required number of local geometry optimisations. The biasing parameters for C60 were $\omega=0.3$ with $\sigma_{c}=6.0$ [29] and for silicon carbide $\omega=3.5$ with $\sigma_{c}=4.5$. The numbers in parenthesis give the speedup with respect to the unbiased runs for corresponding quantities. All simulations were successfully carried out until the ground state structure or in the case of 16 atom silicon carbide a polytype of the ground state structure, was found.  | C60 unbiased | C60 biased | SiC unbiased | SiC biased
---|---|---|---|---
mean | 9254.71 | 370.87 (25) | 9715.36 | 153.84 (63)
std | 7716.87 | 246.74 (31) | 10181.81 | 97.92 (104)
min | 1184 | 89 (13) | 771 | 31 (25)
25% | 4048 | 180 (22) | 3506 | 87 (40)
50% | 6994 | 309 (23) | 5720 | 122 (47)
75% | 11765 | 520 (23) | 12170 | 202 (60)
max | 44136 | 1165 (38) | 56213 | 571 (98)

As expected and as shown in Fig. 2 the distribution of the found structures
with respect to their physical energy $E$ and their degree of symmetry is also
quite different. For the biased MH runs, the fraction of high symmetry
structures is considerably higher and the average physical energy of low
symmetry structures is higher since many low energy defects were not found.
Many of these high symmetry structures are quite interesting. Searching for
structures where all atoms of a certain species have the same environment, we
found for instance several SiC structures where all the carbon atoms were
3-fold coordinated, whereas all the silicon atoms are 4-fold coordinated. Such
a structure is shown in Fig. 3.

Figure 3: Two high symmetry structures found by a biased run. The structure
(a) is only 517 meV/atom higher in energy than the ground state structure,
even though the bonding character is completely different from the ground
state where all atoms are 4-fold coordinated. In this structure all carbon
atoms are 3-fold coordinated. Structure (b) is a 4-fold coordinated silicon
carbide structure with an energy of 6 meV/atom above the ground state. Four
different environments exist for each carbon and silicon atom, but the
environments are so similar that the differences can not be detected by eye.
Carbon atoms are displayed by smaller spheres than silicon atoms for which
each environment has its own colour.

Since our penalty function goes smoothly to zero when the environments get
more similar, we actually also found most low energy structures of silicon
carbide with up to 4 different environments per atom with a penalty function
that favours a single environment. It turned out that in these cases, the
environments tend to be quite similar and result thus in a small instead of a
strictly zero penalty function (see Fig. 3). This finding is related to a
strong correlation between the structural environment diversity as measured by
our penalty function and the total energy as shown in Fig. 4.

Figure 4: Correlation between our penalty function $P_{1}$, that is a measure
of the structural diversity and the total energy $E$. Only if all atoms in the
C60 (a) or silicon carbide (b) structure have similar environments, the energy
of the structure will be low.

In summary, based on a non-conventional measure of symmetry, that is motivated
by Paulings rule of parsimony, we construct a penalty function that measures
the dissimilarity between different atomic environments in a structure. To
construct the penalty function no guesses of which symmetry will be adopted by
the system are required. Adding this penalty function to the physical PES
gives a biased PES where disordered structures are pushed up in energy
relative to high symmetry structures. This leads to a lowering of the downhill
barriers compared to the uphill barriers. This stronger structure seeker
property of the biased potential energy surface allows for much faster
searches for high symmetry ground states.
The penalty function also allows us to find high symmetry structures of higher
energy rapidly. This feature opens the way to perform structure prediction on
a PES that was constructed with a cheap but not very accurate method to find
high symmetry structures in low as well as moderately higher energy regions.
These high symmetry structures can then be reranked by calculating their
energies with a more accurate but also more expensive electronic structure
method. In this way high energy structures that were higher in energy with the
cheap method can become low in energy with the accurate method. This procedure
would not be possible without the bias because in this case the overwhelming
majority of higher energy structures are typically all defect structures,
which are unlikely to become low energy structures when reranked.
Our results also clearly show the general validity of Paulings rule that in
low energy structures the variability of the atomic environments is quite
limited.
Financial support from SNF and computing time from CSCS ( project s963) and
sciCORE (http://scicore.unibas.ch/) are acknowledged. We thank Prof. Alireza
Ghasemi and Prof. Andris Gulans for useful comments on the manuscript.

## References

  * Oganov  _et al._ [2019] A. R. Oganov, C. J. Pickard, Q. Zhu, and R. J. Needs, Nature Reviews Materials 4, 331 (2019).
  * Hilleke  _et al._ [2022] K. P. Hilleke, T. Bi, and E. Zurek, APPLIED PHYSICS A-MATERIALS SCIENCE & PROCESSING 128 (2022), 10.1007/s00339-022-05576-z.
  * Kirkpatrick  _et al._ [1983] S. Kirkpatrick, C. D. Gelatt, and M. P. Vecchi, Science (New York, N.Y.) 220, 671 (1983).
  * Wales and Doye [1997] D. J. Wales and J. P. K. Doye, The Journal of Physical Chemistry A 101, 5111 (1997).
  * Goedecker [2004] S. Goedecker, The Journal of Chemical Physics 120, 9911 (2004), https://doi.org/10.1063/1.1724816 .
  * Amsler and Goedecker [2010] M. Amsler and S. Goedecker, The Journal of chemical physics 133, 224104 (2010).
  * Sicher  _et al._ [2011] M. Sicher, S. Mohr, and S. Goedecker, The Journal of chemical physics 134, 044106 (2011).
  * Roy  _et al._ [2008] S. Roy, S. Goedecker, and V. Hellmann, Physical Review E 77, 056707 (2008).
  * Pickard and Needs [2006] C. J. Pickard and R. J. Needs, Phys. Rev. Lett. 97, 045504 (2006).
  * Martoňák  _et al._ [2007] R. Martoňák, A. R. Oganov, and C. W. Glass, Phase Transitions 80, 277 (2007), https://doi.org/10.1080/01411590701228398 .
  * Johnston [2003] R. L. Johnston, Dalton Transactions , 4193 (2003).
  * Oganov and Glass [2006] A. R. Oganov and C. W. Glass, The Journal of chemical physics 124, 244704 (2006).
  * Bhattacharya  _et al._ [2013] S. Bhattacharya, S. V. Levchenko, L. M. Ghiringhelli, and M. Scheffler, Phys. Rev. Lett. 111, 135501 (2013).
  * Vilhelmsen and Hammer [2014] L. B. Vilhelmsen and B. Hammer, The Journal of chemical physics 141, 044711 (2014).
  * Zhu  _et al._ [2015] Q. Zhu, A. R. Oganov, A. O. Lyakhov, and X. Yu, Physical Review B 92 (2015), https://doi.org/10.1103/PhysRevB.92.024106.
  * Curtis  _et al._ [2018] F. Curtis, X. Li, T. Rose, Á. Vázquez-Mayagoitia, S. Bhattacharya, L. M. Ghiringhelli, and N. Marom, Journal of chemical theory and computation 14, 2246 (2018).
  * Glass  _et al._ [2006] C. W. Glass, A. R. Oganov, and N. Hansen, Computer Physics Communications 175, 713 (2006).
  * Wang  _et al._ [2012] Y. Wang, J. Lv, L. Zhu, and Y. Ma, Comput. Phys. Commun. 183, 2063 (2012).
  * Falls  _et al._ [2021] Z. Falls, P. Avery, X. Wang, K. P. Hilleke, and E. Zurek, The Journal of Physical Chemistry C 125, 1601 (2021).
  * Jørgensen  _et al._ [2018] M. S. Jørgensen, U. F. Larsen, K. W. Jacobsen, and B. Hammer, The journal of physical chemistry. A 122, 1504 (2018).
  * Yamashita  _et al._ [2018] T. Yamashita, N. Sato, H. Kino, T. Miyake, K. Tsuda, and T. Oguchi, Physical Review Materials 2 (2018), https://doi.org/10.1103/PhysRevMaterials.2.013803.
  * Wheeler  _et al._ [2007] S. E. Wheeler, P. v. R. Schleyer, and H. F. Schaefer, The Journal of Chemical Physics 126, 104104 (2007), https://doi.org/10.1063/1.2646940 .
  * Pauling [1929] L. Pauling, Journal of the American Chemical Society 51, 1010 (1929).
  * Pauling [1989] L. Pauling,  _The nature of the chemical bond and the structure of molecules and crystals: An introduction to modern structural chemistry_ , 3rd ed. (Cornell Univ. Pr, Ithaca, NY, 1989).
  * Sørensen  _et al._ [2018] K. H. Sørensen, M. S. Jørgensen, A. Bruix, and B. Hammer, The Journal of Chemical Physics 148, 241734 (2018), https://doi.org/10.1063/1.5023671 .
  * Chiriki  _et al._ [2019] S. Chiriki, M.-P. V. Christiansen, and B. Hammer, Phys. Rev. B 100, 235436 (2019).
  * Cliffe  _et al._ [2017] M. J. Cliffe, A. P. Bartók, R. N. Kerber, C. P. Grey, G. Csányi, and A. L. Goodwin, Phys. Rev. B 95, 224108 (2017).
  * De  _et al._ [2019] D. S. De, M. Krummenacher, B. Schaefer, and S. Goedecker, Physical review letters 123, 206102 (2019).
  * Zhu  _et al._ [2016] L. Zhu, M. Amsler, T. Fuhrer, B. Schaefer, S. Faraji, S. Rostami, S. A. Ghasemi, A. Sadeghi, M. Grauzinyte, C. Wolverton, _et al._ , The Journal of chemical physics 144, 034203 (2016).
  * Parsaeifard  _et al._ [2020] B. Parsaeifard, D. Tomerini, D. S. De, and S. Goedecker, The Journal of chemical physics 153, 214104 (2020).
  * Gubler  _et al._ [2022] M. Gubler, M. Krummenacher, H. Huber, and S. Goedecker, “Efficient variable cell shape geometry optimization,”  (2022).
  * Schaefer  _et al._ [2015] B. Schaefer, S. Alireza Ghasemi, S. Roy, and S. Goedecker, The Journal of chemical physics 142, 034112 (2015).
  * Parrinello and Rahman [1981] M. Parrinello and A. Rahman, Journal of Applied Physics 52, 7182 (1981), https://doi.org/10.1063/1.328693 .
  * Xu and Scuseria [1994] C. Xu and G. E. Scuseria, Physical review letters 72, 669 (1994).
  * Hourahine  _et al._ [2020] B. Hourahine, B. Aradi, V. Blum, F. Bonafé, A. Buccheri, C. Camacho, C. Cevallos, M. Y. Deshaye, T. Dumitrică, A. Dominguez, S. Ehlert, M. Elstner, T. van der Heide, J. Hermann, S. Irle, J. J. Kranz, C. Köhler, T. Kowalczyk, T. Kubař, I. S. Lee, V. Lutsker, R. J. Maurer, S. K. Min, I. Mitchell, C. Negre, T. A. Niehaus, A. M. N. Niklasson, A. J. Page, A. Pecchia, G. Penazzi, M. P. Persson, J. Řezáč, C. G. Sánchez, M. Sternberg, M. Stöhr, F. Stuckenberg, A. Tkatchenko, V. W.-Z. Yu, and T. Frauenheim, The Journal of chemical physics 152, 124101 (2020).
  * Sieck [2000] A. Sieck, “Structure and physical properties of silicon clusters and of vacancy clusters in bulk silicon,” (2000).
  * Miller  _et al._ [2022] M. Miller, D. Wales, and V. de Souza, “disconnectionDPS: Fortran program to generate disconnectivity graphs from stationary point databases.” http://www-wales.ch.cam.ac.uk/software.html (2022).
  * Becker and Karplus [1997] O. M. Becker and M. Karplus, J. Chem. Phys. 106, 1495 (1997).
  * Wales [2004] D. Wales,  _Energy Landscapes: Applications to Clusters, Biomolecules and Glasses_ , 1st ed. (Cambridge University Press, 2004).
  * Zwanzig  _et al._ [1992] R. Zwanzig, A. Szabo, and B. Bagchi, Proceedings of the National Academy of Sciences 89, 20 (1992), https://www.pnas.org/content/89/1/20.full.pdf .

##  Appendix A Symmetry Bias

The goal of the symmetry bias is to find a measure for the structural symmetry
of the system and to use it as a bias on the PES to drive the system faster to
the global minimum during a minima hopping simulation. As a measure for the
structural symmetry of a system we quantify the differences between the
individual atomic environments.
First a matrix $F$ is formed containing the overlap matrix (OM) fingerprints
[29] of every atom $k$ of the cluster or cell as vector $V_{k}$ in columns. In
the OM fingerprint method the eigenvalues of a localized overlap matrix are
assembled into a vector. All entries of each fingerprint vectors $V_{k}$ need
to be sorted before forming the matrix $F$.

$F=\left(\begin{array}[]{ccccc}V_{1}(1)&V_{2}(1)&V_{3}(1)&\cdots&V_{nat}(1)\\\
V_{1}(2)&V_{2}(2)&V_{3}(2)&\cdots&V_{nat}(2)\\\
\vdots&\vdots&\vdots&\ddots&\vdots\\\
V_{1}(l_{fp})&V_{2}(l_{fp})&V_{3}(l_{fp})&\cdots&V_{nat}(l_{fp})\\\
\end{array}\right)$

with $V_{k}(j)$ being the j-th entry of the OM fingerprint and $l_{fp}$ is the
length of the fingerprint vectors.

The Gram matrix

$D=F^{T}F$ (3)

can now be formed.

If all atomic environments are identical, the rank of the matrix $D$ formed by
the OM fingerprint vectors is one, when there are only two different
environments the rank is two, etc. The rank can most easily be calculated from
the eigenvalues $\lambda_{i}$ of the matrix $D$, constructed from the OM
fingerprint vectors. The eigenvalues $\lambda_{i}$ of matrix $D$ are sorted in
descending order, i.e. $\lambda_{1}$ is the largest eigenvalue. The matrix
elements $D_{i,j}$ equals $<V_{i}|V_{j}>$ for the atom pair $(i,j)$.

The number of the non-zero eigenvalues of this matrix gives the rank of the
fingerprint vectors. So the penalty function that favours one single
environment for a certain element is

$P_{1}(R_{1},\ldots,R_{N_{at}})=\sum_{i=2}^{N_{at}}\lambda_{i}=\operatorname{Tr}(D)-\lambda_{1},$
(4)

In case we want to allow for up to $l$ environments, the penalty becomes

$P_{l}(R_{1},\ldots,R_{N_{at}})=\sum_{i=l+1}^{N_{at}}\lambda_{i}=\operatorname{Tr}(D)-\sum_{i=1}^{l}\lambda_{i}$
(5)

where $R_{i}$ is the position of the atom $i$ in the system in Cartesian
coordinates, $N_{at}$ equals the number of atoms in the system and
$\operatorname{Tr}(D)$ is the trace of matrix D.
For a multi-component system, each element contributes its own penalty
function and the total penalty function is the sum of all the elemental
contributions.

##  Appendix B Symmetry bias derivatives

To obtain conservative forces of the biased PES the derivative of the symmetry
bias needs to be added to the physical forces. The same is true for the
derivative of the biased symmetry function with respect to the lattice vectors
which need to be added to the lattice derivatives in the case of PBC.

###  B.1 Symmetry bias derivative with respect to atomic coordinates

The derivative of the symmetry bias function with respect to the atomic
coordinates is

$\frac{\partial P}{\partial
R}=\frac{\partial(\operatorname{Tr}(D)-\lambda_{1})}{\partial
R}=\sum_{i}\left(\frac{\partial D_{i,i}}{\partial
R}\right)-\frac{\partial\lambda_{1}}{\partial R}.$ (6)

For the term $\frac{\partial\lambda_{1}}{\partial R}$ the Hellman Feynman
Theorem is used.

$\frac{\partial\lambda_{1}}{\partial R}=\braket{X_{1}}{\frac{\partial
D}{\partial R}}{X_{1}}$

with $X_{1}$ being the eigenvector belonging to the largest eigenvalue
$\lambda_{1}$ of matrix $D$. The derivative of the dimensionality matrix $D$
$\frac{\partial D}{\partial R}$ depends on the derivatives of the OM
fingerprints.

$\frac{\partial D_{i,j}}{\partial R}=\sum_{l}\frac{\partial V_{i}(l)}{\partial
R}V_{j}(l)+\frac{\partial V_{j}(l)}{\partial R}V_{i}(l)$ (7)

with $l$ counting over all entries in the atomic fingerprint eigenvectors
$V_{i}$. The derivative $\frac{\partial V_{i}}{\partial R}$ is formed with the
help of publication [29].
Now the negative gradient of the derivative can be added to the physical
forces to obtain the biased forces belonging to the biased PES.

###  B.2 Symmetry bias derivative with respect to lattice vectors

Analog to the derivative with respect to the atom positions we can find the
derivative with respect to the lattice vectors

$\frac{\partial{P}}{\partial{{\mathbf{h}}}}=\frac{\partial{(Tr(D)-\lambda_{1})}}{\partial{{\mathbf{h}}}}=\sum_{i}\left(\frac{\partial{D_{i,i}}}{\partial{\mathbf{h}}}\right)-\frac{\partial{\lambda_{1}}}{\partial{\mathbf{h}}}$

with h being the lattice vector matrix

$\mathbf{h}=\begin{bmatrix}h_{1}(1)&h_{2}(1)&h_{3}(1)\\\
h_{1}(2)&h_{2}(2)&h_{3}(2)\\\ h_{1}(3)&h_{2}(3)&h_{3}(3)\\\ \end{bmatrix}$

with $h_{i}$ being the lattice vectors. Like before we can use the the Hellman
Feynman Theorem for the term $\frac{\partial\lambda_{1}}{\partial\mathbf{h}}$.
This results in

$\frac{\partial\lambda_{1}}{\partial\mathbf{h}}=\braket{X_{1}}{\frac{\partial
D}{\partial\mathbf{h}}}{X_{1}}$

The derivative of the matrix entries $D_{i,j}$ with respect to the lattice
vectors is

$\frac{\partial{D_{i,j}}}{\partial{\mathbf{h}}}=\sum^{l_{fp}}_{l}\frac{\partial{V_{i}(l)}}{\partial{\mathbf{h}}}V_{j}(l)+\frac{\partial{V_{j}(l)}}{\partial{\mathbf{h}}}V_{i}(l)$

where we now need the derivative of the OM fingerprints $V_{i}$ with respect
to the lattice vectors.

To calculate the derivative of the OM fingerprints $V_{k}$ with respect to the
lattice vectors we can apply the chain rule so that we can use the already
known derivation of the OM fingerprint with respect to the atomic positions
$\frac{\partial V_{k}}{\partial R}$. It is important to note that the OM
fingerprints $V_{k}$ for atom $k$ in the system is formed by putting Gaussian
type orbitals only on all atoms within a given cutoff radius around the
central atom $k$ and then forming an overlap matrix from them. Therefore, we
only need to consider the atomic positions $\tilde{R}_{j}^{k}$ of all atoms
$j$ in the sphere around the central atom $k$. This leads to the fact that we
now have two counting schemes, one for the atoms in the sphere and one for the
atoms in the main cell. To deal with this we introduce a function index$(i,k)$
that maps atom number $i$ from the sphere counting scheme of the central atom
$k$ to the main cell counting scheme that gives back the index of the
corresponding atom in the main cell counting scheme. This results in

$\frac{\partial{V_{k}(l)}}{\partial{{h}}}=\sum^{N_{sp}}_{j}\frac{\partial{V_{k}(l)}}{\partial{{\tilde{R}_{j}^{k}}}}\frac{\partial{\tilde{R}_{j}^{k}}}{\partial{{h}}}$

with $N_{sp}$ being the number of atoms in a sphere that is formed by the
cutoff radius around the central atom $k$. Since the derivative of the overlap
matrix fingerprint is invariant under the change of the counting scheme we get

$\frac{\partial V_{k}}{\partial\tilde{R}_{j}^{k}}=\frac{\partial
V_{k}}{\partial R_{\text{index}(j,k)}}.$

One needs to be aware of the fact that in periodic boundary conditions it is
possible that multiple images of the same atom from the main cell can be
inside the cutoff radius. The position of the atoms in the sphere around atom
$k$ can then be described as

$\tilde{R}^{k}_{i}=R_{\text{index}(i,k)}+\mathbf{h}\cdot\mathbf{n}_{i,k}=\mathbf{h}\cdot
R_{\text{index}(i,k)}^{\text{frac}}+\mathbf{h}\cdot\mathbf{n}_{i,k}$

with

$\mathbf{n}_{i,k}=\begin{bmatrix}a_{i,k}\\\ b_{i,k}\\\ c_{i,k}\\\
\end{bmatrix}$

being the multiplier for atoms outside the periodic cell and

$R_{\text{index}(i,k)}^{\text{frac}}=\begin{bmatrix}\alpha_{\text{index}(i,k)}\\\
\beta_{\text{index}(i,k)}\\\ \gamma_{\text{index}(i,k)}\\\ \end{bmatrix}$

the atomic coordinates in fractional form of the atom belonging to
$\text{index}(i,k)$ in the main cell.

This results in

$\frac{\partial{\tilde{R}_{i}^{k}}}{\partial{{h}}}=\begin{bmatrix}\alpha_{\text{index}(i,k)}+a_{i,k}&\beta_{\text{index}(i,k)}+b_{i,k}&\gamma_{\text{index}(i,k)}+c_{i,k}\\\
\alpha_{\text{index}(i,k)}+a_{i,k}&\beta_{\text{index}(i,k)}+b_{i,k}&\gamma_{\text{index}(i,k)}+c_{i,k}\\\
\alpha_{\text{index}(i,k)}+a_{i,k}&\beta_{\text{index}(i,k)}+b_{i,k}&\gamma_{\text{index}(i,k)}+c_{i,k}\\\
\end{bmatrix}$

for atom $i$ in the sphere counting scheme around the central atom $k$ in the
main cell counting scheme.


# Single inertial particle statistics in turbulent flows from Lagrangian
velocity models

Jan Friedrich  Institute of Physics and For Wind, University of Oldenburg,
26129 Oldenburg, Germany  Univ. Lyon, ENS de Lyon, Univ. Claude Bernard, CNRS,
Laboratoire de Physique, F-69342, Lyon, France  Bianca Viggiano  Univ. Lyon,
ENS de Lyon, Univ. Claude Bernard, CNRS, Laboratoire de Physique, F-69342,
Lyon, France  Department of Mechanical and Materials Engineering, Portland
State University, Portland, Oregon, USA  Mickael Bourgoin  Univ. Lyon, ENS de
Lyon, Univ. Claude Bernard, CNRS, Laboratoire de Physique, F-69342, Lyon,
France  Raúl Bayoán Cal  Department of Mechanical and Materials Engineering,
Portland State University, Portland, Oregon, USA  Univ. Lyon, ENS de Lyon,
Univ. Claude Bernard, CNRS, Laboratoire de Physique, F-69342, Lyon, France
Laurent Chevillard  Univ. Lyon, ENS de Lyon, Univ. Claude Bernard, CNRS,
Laboratoire de Physique, F-69342, Lyon, France

###### Abstract

We present the extension of a modeling technique for Lagrangian tracer
particles [B. Viggiano et al., J. Fluid Mech.(2020), vol. 900, A27] which
accounts for the effects of particle inertia. Thereby, the particle velocity
for several Stokes numbers is modeled directly by a multi-layered Ornstein-
Uhlenbeck process and a comparison of key statistical quantities (second-order
velocity structure function, acceleration correlation function, and root mean
square acceleration) to expressions derived from Batchelor’s model as well as
to direct numerical simulations (DNS) is performed. In both approaches,
Stokes’ drag is treated by an approximate “linear filter” which replaces the
particle position entering the fluid velocity field by the corresponding ideal
tracer position. Effects of preferential concentration of inertial particles
are taken into account in terms of an effective Stokes number that is
determined from the zero crossing of the acceleration correlation function
from DNS. This approximation thus allows the modeling of inertial particle
statistics through stochastic methods and models for the Lagrangian velocity;
the particle velocity is effectively decoupled from the particle position. In
contrast to the ordinary filtering technique [Cencini et al., J. Turbul.
(2006), 7, N36], our method captures the effects of preferential concentration
of particles at low Stokes numbers which manifest themselves for instance by a
sharp decrease of the acceleration variance for increasing Stokes numbers.

##  I Introduction

Transport processes of particles in turbulent flows play an important role in
turbulence theory [1, 2, 3] and are intimately related to the problems of
turbulent mixing and turbulent diffusion. Recent advances in particle tracking
[4, 5] allow for the accurate quantification of single-particle statistics and
yield important insights into the behavior of marked fluid particles (so-
called Lagrangian tracers) as well as heavier particles, i.e., particles with
finite inertia [6, 7]. For such inertial particles, for instance, numerical
and experimental evidence suggests the segregation of particles into clusters
[8, 9, 10, 11]. This phenomenon, generally known as _preferential
concentration_ , is sometimes explained by particles evading dominant vortical
structures and accumulating in flow regions of high strain [12] although other
mechanisms have been proposed [13, 14]. A comprehensive understanding of these
mechanisms, however, is complicated by the random and multiscale structure of
the turbulent velocity field. On the other hand, such particle-laden flows are
commonly encountered in many industrial and environmental processes, which
creates an abundance of research devoted to characterizing their dynamical
properties [15, 16]. Despite numerous experimental, numerical, and theoretical
works [4, 11, 14], simple generic models which are capable to quantitatively
reproduce basic dynamical Lagrangian properties (e.g., velocity and
acceleration statistics) as well as the combined role of particle inertia and
preferential concentration are still lacking.

Turbulent flow properties can be described either in a fixed frame of
reference, the so-called Eulerian description, or alternatively by the
Lagrangian description, where evolution of the flow is observed temporally
through the trajectories of point particles. The Lagrangian description in
turbulence provides a complete view of particle transport and dispersion which
can be traced back to the seminal works by G.I. Taylor who set the diffusion
problem in the context of fluid element trajectories [17]. Thereby, the
trajectory of an ideal tracer, the so-called Lagrangian path, can be
determined from the first-order ordinary differential equation (ODE)

$\mathbf{\dot{X}}(\mathbf{y},t)\equiv\mathbf{v}(\mathbf{y},t)=\mathbf{u}(\mathbf{X}(\mathbf{y},t),t)\;.$
(1)

where $\mathbf{v}(\mathbf{y},t)$ is the Lagrangian velocity, with the initial
condition $\mathbf{X}(\mathbf{y},0)=\mathbf{y}$. The statistical description
of ideal tracer particles in turbulent flows, which is governed by Eq. (1), is
one of the main concerns of turbulence theory. On the basis of the self-
similar theory of Kolmogorov, moments of the Lagrangian velocity increments
$\delta_{\tau}\mathbf{v}(\mathbf{y},t)=\mathbf{v}(\mathbf{y},t+\tau)-\mathbf{v}(\mathbf{y},t)$
are supposed to scale as $\langle(\delta_{\tau}v)^{n}\rangle\sim\tau^{n/2}$.
Hence, Kolmogorov’s theory suggests that a tracer particle exhibits a
diffusive process (ordinary Brownian motion) in velocity space [1, 18] at
inertial scales. The corresponding Ornstein-Uhlenbeck (OU) process for the
Lagrangian velocity has been extended by Sawford in order to account for
finite Reynolds number effects [19] (see also [20] for further references).
Recently, an infinitely differentiable causal random walk has been proposed
[21]. It can be considered as a layered OU processes, where each layer of
order $n$ ensures the existence of the derivative of order $n$ and causality.
Nonetheless, experimental and numerical evidence suggests that moments exhibit
multifractal scaling $\langle(\delta_{\tau}v)^{n}\rangle\sim\tau^{\zeta_{n}}$
which manifests itself in form of a nonlinear dependence of the scaling
exponents $\zeta_{n}$ on order $n$ [22, 23, 24, 25]. The latter feature is a
direct consequence of the phenomenon of intermittency, reminiscent of the non-
self-similarity of turbulence, which entails strongly non-Gaussian behavior of
the Lagrangian velocity increment probability density function at small time
separations $\tau$. Precisely those multiscaling (non-Gaussian) features
impose a major challenge for the modeling of particle transport in turbulent
flows [26, 21].

Once a particle gains inertia, either from its size or its density in
comparison to the surrounding fluid with kinematic viscosity $\nu$, the
dynamics of the particle is governed by additional forces [27, 12]. For
example, there are now effects due to the Stokes’ drag force based on the
relative velocities of the fluid and the particles and an added mass force
which is caused by the displacement of the fluid from the motion of the
particle, to name a few. If we only retain Stokes’ drag and disregard these
other terms (which is typically assumed to be the case when particles are
rather heavy and small), the particle’s velocity is determined by the
following set of equations [27, 12, 28]

$\displaystyle\mathbf{\dot{X}}_{p}(\mathbf{y},t)=$
$\displaystyle\mathbf{v}_{p}(\mathbf{y},t)\;,$ (2)
$\displaystyle\mathbf{\dot{v}}_{p}(\mathbf{y},t)=$
$\displaystyle\frac{1}{\tau_{p}}\left[\mathbf{u}(\mathbf{X}_{p}(\mathbf{y},t),t)-\mathbf{v}_{p}(\mathbf{y},t)\right]\;,$
(3)

where $\mathbf{X}_{p}(\mathbf{y},t)$ is the particle position,
$\mathbf{v}_{p}(\mathbf{y},t)$ is the particle velocity, and
$\mathbf{u}(\mathbf{x},t)$ is the fluid velocity field. Initial conditions
fully determine the future evolution of the particle position and we assume
them to be given by $\mathbf{X}_{p}(\mathbf{y},0)=\mathbf{y}$ and
$\mathbf{v}_{p}(\mathbf{y},0)=\mathbf{u}(\mathbf{y},0)$, where the matching of
tracer and inertial particle position at $t=0$ has been imposed arbitrarily.
The particle response time, $\tau_{p}$, is defined as
$\tau_{p}=m_{p}/3\pi\mu_{f}d_{p}$ where $m_{p}$ and $d_{p}$ denote the mass
and the diameter of the particle, respectively, and $\mu_{f}$ is the viscosity
of the fluid.

In this simplified model, the Stokes number, defined as St$=\tau_{p}/\tau_{K}$
(where $\tau_{K}=\left(\nu/\langle\varepsilon\rangle\right)^{1/2}$ is the
Kolmogorov time scale and $\langle\varepsilon\rangle$ denotes the averaged
local energy dissipation rate) is the relevant parameter to characterize
particles inertia. The Stokes number thus describes the ratio of the particle
response time $\tau_{p}$ to the Kolmogorov time scale $\tau_{K}$ at which
viscous forces dominate the flow. In particular, for St$=0$ (i.e., for the
case of an overdamped inertial particle) we recover Eq. (1) and the particle
behaves as an ideal Lagrangian tracer.

The main purpose of this paper is to generalize recent stochastic and
multifractal modeling approaches [21] \- which were initially devised for the
Lagrangian velocity (1) - in order to include finite inertia effects. The
simplest approach in this direction is to investigate the particle dynamics
under the assumption that each inertial particle samples homogeneously the
carrier flow, and hence the explored Lagrangian flow velocity
$\mathbf{u}(\mathbf{X}_{p}(\mathbf{y},t),t)$ in Eq. (3) is representative of
the complete background turbulence regardless of any eventual preferential
concentration mechanism. This approximation, determines the particle velocity
$\mathbf{v}_{p}(\mathbf{y},t)$ by a _linear filter_ of the Lagrangian velocity
$\mathbf{v}(\mathbf{y},t)$ (see Eq. (3)) and has been invoked for the first
time by Tchen [29] and Hinze [30]. In particular, under the assumption of an
exponential decay of the Lagrangian velocity correlation function for the
tracers of the carrier flow (what would correspond for instance to a simple OU
dynamics for the tracers), the Tchen-Hinze theory relates the variance of the
particle velocity $\sigma_{\textbf{v}_{p}}^{2}$ directly to that of the ideal
tracer $\sigma_{\textbf{v}}^{2}$ according to

$\sigma_{\textbf{v}_{p}}^{2}=\frac{T_{L}^{2}}{T_{L}^{2}-\tau_{p}^{2}}\sigma_{\textbf{v}}^{2}\;,$
(4)

where $T_{L}$ denotes the Lagrangian integral time scale. An extension of this
linear filter approximation has been suggested by Deutsch and Simonin [31, 32]
who argued that fluid quantities in Eq. (4) do not take into account the fact
that the Lagrangian dynamics of the tracers probed at the position of the
inertial particles may deviate from the global Lagrangian dynamics of the
carrier flow, due to preferential concentration. Hence, they propose to
replace the Lagrangian integral time scale $T_{L}$ and variance
$\sigma_{\textbf{v}}^{2}$ by quantities determined _at the position of
inertial particles_ ($T_{L@p}$ and $\sigma_{\textbf{v}@p}^{2}$). This
refinement improves the agreement of Eq. (4) with direct numerical simulations
[33], in particular for particles with small Stokes numbers, which seem to be
more sensitive to preferential sampling effects. From a practical point of
view, such a refinement requires however to determine $T_{L@p}$ and
$\sigma_{\textbf{v}@p}^{2}$ on the basis of statistics of trajectories of
individual inertial particle and its coincident Lagrangian particle, which is
a quite challenging task, hardly feasible in experiments. Furthermore, the
framework by Tchen-Hinze as well as the refined theory [31, 32] does not
reproduce the observed behavior of statistical quantities at small time lags.
This applies in particular to the particle’s acceleration which remains non-
differentiable under the assumption of a simple exponential correlation
function for the Lagrangian velocity, in contrast to empirical evidence which
shows smooth behavior. These limitations have also been discussed recently in
the broader context of particle settling in turbulent flows [34].

In this context, the purpose of the present work thus is twofold: _i.)_ we
generalize the Tchen-Hinze theory to a stochastic process which is infinitely
differentiable (similarly to the framework which has been devised for the
Lagrangian velocity [21]) and, hence, contains important information on small-
scale fluctuations, and, _ii.)_ we take into account the effects of
preferential concentration (e.g., the correct amplitudes of particle
accelerations [35]) by the introduction of an effective Stokes number, which
is determined on the basis of the zero crossing of the empirically determined
acceleration correlation function.

The paper is organized as follows: Section II discusses the implications of
the linear filter approximation for the particle velocity correlation
function. Moreover, we provide a comparison between trajectories obtained
under this approximation, and the ones obtained from direct numerical
simulations (DNS) of turbulence. We present then models for particles
correlation functions based on filtered versions of _i.)_ an infinitely
differentiable causal random walk in Section III and _ii.)_ by the Batchelor
model in Section IV. These predictions are directly compared to DNS in Section
V and a summary of our results is included in Section VI.

##  II Particle response based on the linear filter approximation

Inertial particle motion at low Stokes numbers, as mentioned in the
introduction, is determined by the system of first order ODEs (2-3) and
requires the knowledge of the full spatiotemporal (Eulerian) fluid velocity
field in Eq. (3). In the following, an approximation will be discussed, which
replaces the Eulerian fluid velocity at the position of the inertial particle
$\mathbf{X}_{p}(\mathbf{y},t)$ by its value at the position of the tracer,
i.e. by the Lagrangian velocity. In this approximation, particle velocity is
obtained by the linear filtering of tracer velocity.

###  II.1 Linear filtering of Lagrangian velocity

In order to model inertial particle statistics on the basis of the Lagrangian
velocity (1), an approximation of the coupled system of first-order ODEs is
invoked which can be termed “linear filtering of the particle velocity” due to
its analogy to methods from signal processing [35]. In this approximation, the
particle position $\mathbf{X}_{p}(\mathbf{y},t)$ that enters the fluid
velocity field in Eq. (3) is approximated as the position of the ideal tracer
$\mathbf{X}(\mathbf{y},t)$ whose temporal evolution is governed by equation
(1). Therefore, Eqs. (2-3) are approximated by

$\displaystyle\mathbf{\dot{X}}_{p}(\mathbf{y},t)=$
$\displaystyle\mathbf{v}_{p}(\mathbf{y},t)\;,$ (5)
$\displaystyle\mathbf{\dot{v}}_{p}(\mathbf{y},t)=$
$\displaystyle\frac{1}{\tau_{p}}\left[\mathbf{u}(\mathbf{X}(\mathbf{y},t),t)-\mathbf{v}_{p}(\mathbf{y},t)\right]\;.$
(6)

Hence, in this linear filter approximation, the particle velocity is
effectively decoupled from the particle position. Therefore, the temporal
evolution of the particle velocity is solely determined by the Lagrangian
velocity $\mathbf{v}(\mathbf{y},t)$ along the tracer trajectory starting from
the initial position of the inertial particle
$\mathbf{X}_{p}(\mathbf{y},0)=\mathbf{y}$. To some extent, the linear filter
neglects the spatiotemporal organization of the fluid velocity, and thus the
segregation of inertial particles in regions of low vorticity, as suggested by
the phenomenon of preferential concentration. Nevertheless, under this
approximation, the evolution equation for the particle velocity (6) can be
solved according to

$\mathbf{v}_{p}(\mathbf{y},t)=\mathbf{v}(\mathbf{y},0)e^{-t/\tau_{p}}+\frac{1}{\tau_{p}}\int_{0}^{t}\textrm{d}t^{\prime}e^{-(t-t^{\prime})/\tau_{p}}\mathbf{v}(\mathbf{y},t^{\prime})\;.$
(7)

The particle position is thus determined as

$\mathbf{X}_{p}(\mathbf{y},t)=\mathbf{y}+\tau_{p}\mathbf{v}(\mathbf{y},0)(1-e^{-t/\tau_{p}})+\frac{1}{\tau_{p}}\int_{0}^{t}\textrm{d}t^{\prime}\int_{0}^{t^{\prime}}\textrm{d}t^{\prime\prime}e^{-(t^{\prime}-t^{\prime\prime})/\tau_{p}}\mathbf{v}(\mathbf{y},t^{\prime\prime})\;.$
(8)

The implications of this linear filter approximation for particle velocity and
acceleration correlation functions will be discussed in the following section.

###  II.2 Linear filtering of correlation functions

In this section, the linear filter approximation (6) is applied in order to
establish relations between inertial particle velocity correlation functions
and the Lagrangian velocity correlation function. To this end, without loss of
generality, we denote $v_{p}$ as of the component of the particle velocity
vector, and correspondingly $a_{p}$ any component of the particle acceleration
vector. Furthermore, we assume that the particle velocity has reached a
statistically stationary state, such that we can neglect the dependence on its
initial position. Eq. (6) can then be integrated according to

$\displaystyle
v_{p}(t)=\frac{1}{\tau_{p}}\int_{-\infty}^{t}\text{d}t^{\prime}e^{-(t-t^{\prime})/\tau_{p}}u(\bm{X}(t^{\prime}),t^{\prime})=\frac{1}{\tau_{p}}\int_{-\infty}^{t}\text{d}t^{\prime}e^{-(t-t^{\prime})/\tau_{p}}v(t^{\prime})=\frac{1}{\tau_{p}}\int_{-\infty}^{+\infty}\text{d}t^{\prime}g_{\tau_{p}}(t-t^{\prime})v(t^{\prime})\;,$
(9)

where $g_{\tau_{p}}(t)=e^{-t/\tau_{p}}1_{t\geq 0}$. The linear filter
approximation for inertial particle velocities based on the trajectory of
individual tracer particles can also be applied to the correlation functions
of velocity and acceleration. To do so, $\mathcal{C}_{v}(\tau)=\langle
v(t)v(t+\tau)\rangle$ defines the correlation function of the velocity of
tracers where $v(t)$ is any component of the tracer velocity. Similarly, the
correlation function of the inertial particle velocity $v_{p}(t)$ is defined
as $\mathcal{C}_{v_{p}}(\tau)=\langle v_{p}(t)v_{p}(t+\tau)\rangle$ and we
obtain, in the statistically stationary range,

$\displaystyle\mathcal{C}_{v_{p}}(\tau)$
$\displaystyle=\frac{1}{\tau_{p}^{2}}(G_{\tau_{p}}\star\mathcal{C}_{v})(\tau)=\frac{1}{\tau_{p}^{2}}\int_{-\infty}^{+\infty}\text{d}tG_{\tau_{p}}(\tau+t)\mathcal{C}_{v}(t)$
$\displaystyle=\frac{1}{\tau_{p}^{2}}\int_{0}^{+\infty}\text{d}t\left[G_{\tau_{p}}(\tau+t)+G_{\tau_{p}}(\tau-t)\right]\mathcal{C}_{v}(t)=\frac{1}{2\tau_{p}}\int_{0}^{+\infty}\textrm{d}t\left[e^{-|\tau+t|/\tau_{p}}+e^{-|\tau-t|/\tau_{p}}\right]\mathcal{C}_{v}(t)\;,$
(10)

where the kernel $G_{\tau_{p}}(t)=(g_{\tau_{p}}\star
g_{\tau_{p}})(t)=\frac{\tau_{p}}{2}e^{-|t|/\tau_{p}}$ was introduced using the
parity of correlation functions. Furthermore, the convolution product is
defined as

$(g_{1}\star
g_{2})(\tau)=\int_{-\infty}^{\infty}\textrm{d}tg_{1}(t)g_{2}(t-\tau)\;.$ (11)

The Fourier representation of Eq. (10) has been proposed for the first time by
Tchen [29]. The determination for the inertial particle acceleration
correlation function $\mathcal{C}_{a_{p}}(\tau)$ based on the filtering of the
tracer acceleration $\mathcal{C}_{a}(\tau)$ is done in a similar fashion,
resulting in

$\mathcal{C}_{a_{p}}(\tau)=\frac{1}{2\tau_{p}}\int_{0}^{+\infty}\text{d}t\left[e^{-|\tau+t|/\tau_{p}}+e^{-|\tau-t|/\tau_{p}}\right]\mathcal{C}_{a}(t).$
(12)

The linear filter approximation thus allows for direct assessment of velocity
and acceleration correlation functions via the corresponding tracer
correlation functions $\mathcal{C}_{v}(t)$ and $\mathcal{C}_{a}(t)$. The
implications of this approximation will be addressed in the following section
where ideal tracer, filtered inertial particles, and true inertial particles
will be compared in DNS of turbulence.

###  II.3 Direct numerical simulations of Lagrangian trajectories, inertial
particles, and comparison to linear filter approach

Twenty-one datasets containing particle trajectories for varying Stokes
numbers have been generated through JHTDB (http://turbulence.pha.jhu.edu)
[36]. The spatiotemporal evolution of the fluid velocity has been obtained by
solving the Navier-Stokes equations in a periodic box with a resolution of
10243. The equation of motion of inertial particles (2-3) are solved by a
second order Runge-Kutta scheme and each DNS subset of data contains $32^{3}$
trajectories. The relevant turbulence quantities are summarized in Tab. 1.

Reλ | $\nu$ | $\langle\varepsilon\rangle$ | d$t$ | $\tau_{K}$ | $\tau_{p}$ | $T_{L}$ | $T_{tot}$ | $N_{p}$ | $N^{3}$
---|---|---|---|---|---|---|---|---|---
418 | $1.85\times 10^{-4}$ | 0.103 | 0.002 | 0.0424 | St$\times\tau_{K}$ | 1.3003 | 7.692$T_{L}$ | $32^{3}$ | $1024^{3}$
Table 1: Characteristic parameters of the direct numerical simulations of of
inertial particles : Taylor-Reynolds number $\textrm{Re}_{\lambda}$, kinematic
viscosity $\nu$, averaged kinetic energy dissipation rate
$\langle\varepsilon\rangle$, time step of simulation , dissipation time
$\tau_{K}=\left(\frac{\nu}{\langle\varepsilon\rangle}\right)^{1/2}$, particle
response time $\tau_{p}$ with St $=[0,0.1,\ldots,1.9,2]$, Lagrangian integral
time $T_{L}$, total time of simulation $T_{tot}$, number of particles in each
simulation $N_{p}$, and resolution $N$ of the periodic simulation domain.

The temporal evolution of particle position $\mathbf{X}_{p}(t)$, a velocity
component $v_{p}(t)$, and the corresponding acceleration $a_{p}(t)$ for
St$=0.2$ are presented in Figure 1 in red for a time span of $\approx
110\tau_{K}$. For comparison, Figure 1(b), shows the trajectory of an ideal
tracer $\mathbf{X}(\mathbf{y},t)$ (blue) starting from the same initial
condition $\mathbf{X}_{p}(0)=\mathbf{y}$. Whereas the tracer particle’s
velocity (blue) exhibits several strong oscillations, the inertial particle’s
velocity (red) seems not so much affected. Indeed, from Figure 1(b) one can
deduce that - already at such low Stokes numbers - the trajectory of the
inertial particle follows a substantially different path. In the context of
preferential concentration [35, 37], one could interpret this in terms of the
inertial particle evading strong vortical flow structures. This might also be
supported by the evolution of the particle’s acceleration whose amplitude (and
statistically speaking its variance) is significantly decreased in comparison
to the tracer particle.

Figure 1: (a) Velocity $v_{x}(t)$ (upper panel) and acceleration $a_{x}(t)$
(lower panel) of tracers (blue), filtered tracers (orange and green) and true
inertial particles for $\textrm{St}=0.2$. (b) Trajectories of tracer, filtered
tracer, and inertial particles.

The linear filter approximation from Eqs. (5-6) with $\tau_{p}=0.2\tau_{K}$
(St$=0.2$, orange curve) remains very close to its determining tracer
trajectory (blue). Especially the temporal evolution of the acceleration
suggests that the linear filter (orange) overestimates the acceleration
variance in comparison to the ordinary inertial particle. Cencini $et\,al.$
[35] applied the same technique to compare root-mean-squared accelerations and
observed large discrepancies between filtered and true $a_{rms}$ at low St.
Nonetheless, for larger Stokes numbers, true and filtered $a_{rms}$-values
seemed to approach one another. Hence, this study suggests that non-trivial
effects of the preferential concentration, when inertia is introduced to the
particles, creates the gap between true and filtered results for
$\textrm{St}<$1\. On that account, Figure 2(a) depicts the acceleration auto-
correlation function which has been obtained by averaging over all
$N_{p}=32^{3}$ particle trajectories. Here, the dash-dotted curves correspond
to the linear filter approximation (12) for $\tau_{p}=\text{St}\tau_{K}$ with
St$=[0.1,0.2,0.5,1,2]$ from left to right. Thereby, the fluid velocity
correlation function $\mathcal{C}_{a}(\tau)$ was integrated over the entire
time range of the simulation $T$ (see also Tab. 1). In agreement with the
findings by Cencini $et\,al.$ [35], the modeling breaks down initially,
drastically underestimating the correlation of the particle for
$\textrm{St}\leq$1, but recovers to generate meaningful statistics for St=2.
Especially, the regime of low Stokes numbers, e.g., for St=0.1, suggests that
the filter approximation (first dash-dotted curve from left) is nearly
identical to the ideal tracer St=0 (not shown in plot) and decorrelates much
faster than the true DNS for St=0.1 (orange curve).

These profound changes between inertial particle and tracers (or filtered
tracers) acceleration properties at low Stokes also manifest themselves in the
root mean square of acceleration itself. Figure 3 depicts the root mean square
values of acceleration $a_{rms}$ as a function of St for inertial particles
(blue) and filtered tracer particles (orange) according to Eq. (9). The
$a_{rms}$-values of inertial particles decreases much faster than their
filtered counterparts. These strong discrepancies were interpreted by Cencini
$et\,al.$ [35] in terms of inertial particles which preferentially sample
regions of low turbulence intensity (or depleted vorticity regions) whereas
the filtered tracer particles are still impacted by strong acceleration events
of tracer trajectories trapped in vortical structures. In other words, by
restricting itself to individual tracer trajectories
$\mathbf{X}(\mathbf{y},t)$, the filtering approach bears no information on the
spatial organization of the surrounding fluid velocity field, which apparently
is crucial for a better understanding of the dynamics of inertial particles.

Figure 2: (a) Acceleration auto-correlation from DNS for ideal tracer (blue)
and inertial particles with Stokes numbers $\textrm{St}=0.1,0.2,0.5,1.,2.$ The
dash-dotted lines correspond to the linear filter approximation (5-6) with
$\tau_{p}=\textrm{St}\tau_{K}$ and underestimate the correlations at small St
in comparison to their DNS counterparts. The linear filter approximation
approaches the DNS acceleration correlation function for larger
$\textrm{St}=2$. (b) Same as in (a) but for a filter with an effective Stokes
number $\textrm{St}^{*}=\tau_{p}^{*}/\tau_{K}$ which has been determined from
the zero crossing of the acceleration correlation from DNS (dash-dotted curves
and DNS curves now possess the same zero crossing). The agreement between DNS
and filtered curves is more accurate than in (a) for small St.

In the following, we suggest a modification of the linear filter which is in
better agreement with DNS but at the same time still is solely based on the
Lagrangian velocity $\mathbf{v}(\mathbf{y},t)$ of an individual tracer. To
this end, we introduce an effective particle response time $\tau_{p}^{*}$ in
the linear filter approximation (6) which does not necessarily obey the usual
relation $\tau_{p}=\text{St}\tau_{K}$ but henceforward is considered as a free
model parameter. In the spirit of recent stochastic models for Lagrangian
tracers [21], which identified the zero crossing of the acceleration
correlation as a crucial ingredient for model calibration, we proceed in a
similar way: First, we determine the zero crossing $\tau_{0}(\textrm{St})$ of
the acceleration auto-correlation functions from DNS in Figure 2.
Subsequently, with the help of a suitable roots-finding algorithm, we
determine the effective particle response time $\tau_{p}^{*}(\textrm{St})$ in
Eq. (12) in such a way that $\mathcal{C}_{a_{p}}(\tau_{0}(\textrm{St}))=0$.

Figure 2(b) depicts the acceleration auto-correlation functions (dash-dotted
curves) after this calibration. Hence, by virtue of an effective particle
response time $\tau_{p}^{*}$ based on matching zero crossing of the
acceleration correlation functions, a better agreement with DNS has been
achieved, particularly for low Stokes numbers. The corresponding filtered
trajectories with St${}^{*}=\tau_{p}^{*}/\tau_{K}$ are also included in Figure
1 and correspond to the green curves. Due to the increased damping
($\tau_{p}^{*}>\tau_{p}$) in the filter, $\tau_{p}^{*}$-filtered accelerations
(green) are closer to their DNS counterparts (red). In conclusion, the
introduction of an effective particle response time $\tau_{p}^{*}$ based on
the zero crossing of the acceleration correlation function reinstates the
effects of preferential sampling of inertial particles to some extend. This
notion will be further assessed in the following section by example of
stochastic models for inertial particle statistics.

Figure 3: Root mean square values of acceleration $a_{rms}$ as a function of
the Stokes number for inertial particles in DNS (blue) and filtered tracer
particles (orange) according to Eq. (6). The $a_{rms}$-values from DNS exhibit
a rapid decrease for small but finite St whereas the linear filter
approximation decays much slower.

##  III Process for inertial particles subject to Gaussian infinitely
differentiable fluid velocity

This section is devoted to a generalization of a recently developed model for
the Lagrangian velocity [21] which will be extended to take into account
finite particle inertia. The simplest model for the Lagrangian velocity $v(t)$
consists of an Ornstein-Uhlenbeck process which obeys the following Langevin
equation

$\textrm{d}v(t)=-\frac{1}{T}v(t)\textrm{d}t+\sqrt{\frac{2\sigma^{2}}{T}}\textrm{d}W(t)\;,$
(13)

where $\textrm{d}W(t)$ denote the increments of a Wiener process $W(t)$ with
zero mean and variance $\textrm{d}t$, and $\sigma^{2}$ denotes the variance of
the process $\langle v^{2}\rangle$. Furthermore, $T$ is a large turbulence
time scale. In this framework, velocity is non-differentiable and thus
acceleration only has a meaning in a distributional sense, which is at odds
with experimental and numerical observations. This model has been extended by
Sawford [19] in order to allow for the differentiability of the velocity,
however, the acceleration of the process remains non-differentiable itself.
Therefore, a multilayered Ornstein-Uhlenbeck process has been proposed where
the velocity itself is infinitely differentiable when the number of layers
tends to infinity [21]. In making use of the linear filter approximation
discussed in the previous section, the latter model can readily be extended in
order to allow for the effects of finite particle inertia. Following the
approach of [21], we consider the system of coupled stochastic differential
equations

$\displaystyle\dot{v}_{p}=$
$\displaystyle\frac{1}{\tau_{p}^{*}}\left(v-v_{p}\right)\;,$ (14)
$\displaystyle\dot{v}=$ $\displaystyle-\frac{1}{T}v+f_{\tau_{\eta}}\;,$ (15)

where $\tau_{p}^{*}$ denotes the particle response time of the model and
$\tau_{\eta}$ is a small-scale turbulence time scale of the order of the
Kolmogorov dissipative time scale $\tau_{K}$ defined in the caption of table
1. The forcing term $f_{\tau_{\eta}}$ is a zero-average Gaussian random force,
which is fully determined by its covariance in the statistically stationary
regime. Following the procedure outlined in [21], the covariance in both,
physical and spectral space, reads

$\displaystyle\mathcal{C}_{f_{\tau_{\eta}}}(\tau)\equiv\langle
f_{\tau_{\eta}}(t)f_{\tau_{\eta}}(t+\tau)\rangle=\frac{\sigma^{2}e^{-\tau_{\eta}^{2}/T^{2}}}{\sqrt{\pi}T\tau_{\eta}\operatorname{erfc}\left(\tau_{\eta}/T\right)}e^{-\frac{\tau^{2}}{4\tau_{\eta}^{2}}}=\frac{2\sigma^{2}e^{-\tau_{\eta}^{2}/T^{2}}}{T\operatorname{erfc}\left(\tau_{\eta}/T\right)}\int\textrm{d}\omega\;e^{2i\pi\omega\tau}e^{-4\pi^{2}\tau_{\eta}^{2}\omega^{2}}\;,$
(16)

where we have introduced $\operatorname{erfc}(x)=1-\operatorname{erf}(x)$ as
the complementary of the error function $\operatorname{erf}(x)$. As we can see
from the structure of its covariance (16), the Gaussian forcing term
$f_{\tau_{\eta}}$ is correlated over the dissipative time scale $\tau_{\eta}$,
the correlation function being itself a Gaussian function. Moreover, this
covariance structure shows that the random process $f_{\tau_{\eta}}(t)$ is
infinitely differentiable, which is a direct consequence of the smoothness of
its Gaussian shape. As developed in Ref. [21], the remaining parameter
$\sigma^{2}$ entering Eq. 16 corresponds to the variance of the Lagrangian
velocity $v$, that is the variance of the solution of Eq. (15) in the
statistically steady range

$\displaystyle\sigma^{2}=\langle v^{2}\rangle.$ (17)

At this stage, causality of the forcing term $f_{\tau_{\eta}}$ is not obvious,
and cannot be deduced from the covariance
$\mathcal{C}_{f_{\tau_{\eta}}}(\tau)$ (Eq. 16). It is indeed shown in Ref.
[21] that this forcing term can be obtained as an infinite sequence of causal
equations, i.e., an infinite number of embedded Ornstein-Uhlenbeck processes,
making the asymptotic process $f_{\tau_{\eta}}$ itself causal. Finally, the
covariance function of the particle velocity $v_{p}$, defined as successive
linear operations (Eqs. 14 and 15) on the forcing term $f_{\tau_{\eta}}$, can
be obtained as

$\mathcal{C}_{v_{p}}(\tau)\equiv\langle
v_{p}(t)v_{p}(t+\tau)\rangle=\frac{2\sigma^{2}e^{-\tau_{\eta}^{2}/T^{2}}}{T\operatorname{erfc}\left(\tau_{\eta}/T\right)}\int\textrm{d}\omega\;e^{2i\pi\omega\tau}\frac{1}{1+4\pi^{2}(\tau_{p}^{*})^{2}\omega^{2}}\frac{T^{2}}{1+4\pi^{2}T^{2}\omega^{2}}e^{-4\pi^{2}\tau_{\eta}^{2}\omega^{2}}\;,$
(18)

which can be calculated explicitly in physical space according to

$\displaystyle\mathcal{C}_{v_{p}}(\tau)=$
$\displaystyle\frac{\sigma^{2}T}{2(T^{2}-(\tau_{p}^{*})^{2})\operatorname{erfc}(\tau_{\eta}/T)}\Bigg{\\{}Te^{-|\tau|/T}\left[\operatorname{erfc}\left(\frac{\tau_{\eta}}{T}-\frac{|\tau|}{2\tau_{\eta}}\right)+e^{2|\tau|/T}\operatorname{erfc}\left(\frac{|\tau|}{2\tau_{\eta}}+\frac{\tau_{\eta}}{T}\right)\right]$
$\displaystyle-\tau_{p}^{*}e^{-|\tau|/\tau_{p}^{*}}e^{\frac{\tau_{\eta}^{2}}{(\tau_{p}^{*})^{2}}-\frac{\tau_{\eta}^{2}}{T^{2}}}\left[\operatorname{erfc}\left(\frac{\tau_{\eta}}{\tau_{p}^{*}}-\frac{|\tau|}{2\tau_{\eta}}\right)+e^{2|\tau|/\tau_{p}^{*}}\operatorname{erfc}\left(\frac{|\tau|}{2\tau_{\eta}}+\frac{\tau_{\eta}}{\tau_{p}^{*}}\right)\right]\Bigg{\\}}.$
(19)

Herein, the dependence of the variance of the particle velocity on
$\tau_{p}^{*}$ can be determined. Setting $\tau=0$ in equation (19) yields

$\sigma_{v_{p}}^{2}=\frac{\sigma^{2}T^{2}}{T^{2}-(\tau_{p}^{*})^{2}}\left[1-\frac{\tau_{p}^{*}}{T}e^{\frac{\tau_{\eta}^{2}}{(\tau_{p}^{*})^{2}}-\frac{\tau_{\eta}^{2}}{T^{2}}}\frac{\operatorname{erfc}\left(\frac{\tau_{\eta}}{\tau_{p}^{*}}\right)}{\operatorname{erfc}\left(\frac{\tau_{\eta}}{T}\right)}\right]\;.$
(20)

Hence, in the limit $\tau_{p}^{*}\rightarrow 0$, the variance of the
Lagrangian model, $\lim_{\tau_{p}^{*}\rightarrow
0}\sigma_{v_{p}}^{2}=\sigma^{2}$, is recovered [21]. Furthermore, the
prediction of the “single-layered” Tchen-Hinze theory (4) can be recovered by
discarding the angular brackets in Eq. (20).
The acceleration correlation function can be calculated from Eq. (19)
according to

$\mathcal{C}_{a_{p}}(\tau)=-\frac{\textrm{d}^{2}\mathcal{C}_{v_{p}}(\tau)}{\textrm{d}\tau^{2}},$
(21)

which yields

$\displaystyle\mathcal{C}_{a_{p}}(\tau)=$
$\displaystyle\frac{\sigma^{2}T}{2(T^{2}-(\tau_{p}^{*})^{2})\operatorname{erfc}(\tau_{\eta}/T)}\Bigg{\\{}\frac{1}{T}e^{-|\tau|/T}\left[\frac{2T}{\sqrt{\pi}\tau_{\eta}}e^{-\left(\frac{\tau_{\eta}}{T}-\frac{|\tau|}{2\tau_{\eta}}\right)^{2}}-\operatorname{erfc}\left(\frac{\tau_{\eta}}{T}-\frac{|\tau|}{2\tau_{\eta}}\right)-e^{2|\tau|/T}\operatorname{erfc}\left(\frac{|\tau|}{2\tau_{\eta}}+\frac{\tau_{\eta}}{T}\right)\right]$
$\displaystyle-\frac{1}{\tau_{p}^{*}}e^{-|\tau|/\tau_{p}^{*}}e^{\frac{\tau_{\eta}^{2}}{(\tau_{p}^{*})^{2}}-\frac{\tau_{\eta}^{2}}{T^{2}}}\left[\frac{2\tau_{p}^{*}}{\sqrt{\pi}\tau_{\eta}}e^{-\left(\frac{\tau_{\eta}}{\tau_{p}^{*}}-\frac{|\tau|}{2\tau_{\eta}}\right)^{2}}-\operatorname{erfc}\left(\frac{\tau_{\eta}}{\tau_{p}^{*}}-\frac{|\tau|}{2\tau_{\eta}}\right)-e^{2|\tau|/\tau_{p}^{*}}\operatorname{erfc}\left(\frac{|\tau|}{2\tau_{\eta}}+\frac{\tau_{\eta}}{\tau_{p}^{*}}\right)\right]\Bigg{\\}}.$
(22)

The variance of the acceleration can be determined according to

$\sigma_{a_{p}}^{2}=\frac{\sigma^{2}}{T^{2}-(\tau_{p}^{*})^{2}}\left[\frac{T}{\tau_{p}^{*}}\frac{e^{\frac{\tau_{\eta}^{2}}{(\tau_{p}^{*})^{2}}}}{e^{\frac{\tau_{\eta}^{2}}{T^{2}}}}\frac{\operatorname{erfc}\left(\frac{\tau_{\eta}}{\tau_{p}^{*}}\right)}{\operatorname{erfc}\left(\frac{\tau_{\eta}}{T}\right)}-1\right]\;.$
(23)

In the limit of $\tau_{p}^{*}\rightarrow 0$, the reduced tracer model
discussed in [21] is recovered and we get

$\lim_{\tau_{p}^{*}\rightarrow
0}\sigma_{a_{p}}^{2}=\frac{\sigma^{2}}{T^{2}}\left[\frac{T}{\sqrt{\pi}\tau_{\eta}}\frac{e^{-\frac{\tau_{\eta}^{2}}{T^{2}}}}{\operatorname{erfc}\left(\frac{\tau_{\eta}}{T}\right)}-1\right]\;.$
(24)

So far, the proposed model possesses Gaussian properties, i.e., the particle
velocity statistics is fully determined by the correlation function (19). The
inclusion of non-Gaussian properties into the stochastic process can be
achieved in the same manner discussed in [21] and will be the subject of
future work. In the present model, the fluid properties enter through the
turnover time $T$ and the small time scale $\tau_{\eta}$, which can be
parameterized in order to match certain characteristics of Lagrangian tracers.
In [21], for instance, $T$ and $\tau_{\eta}$ were determined in order to match
the integral Lagrangian time scale

$T_{L}=\int_{0}^{\infty}\textrm{d}\tau\frac{\mathcal{C}_{v}(\tau)}{\mathcal{C}_{v}(0)}\;,$
(25)

and the zero crossing of the empirically determined acceleration correlation
function $\mathcal{C}_{a}(\tau)$. In addition to this parameterization, here,
we also have to determine the particle response time $\tau_{p}^{*}$. In the
spirit of the effective particle response time of the linear filter
approximation discussed in Section II, $\tau_{p}^{*}$ will be determined in
order to match the zero crossing of the acceleration correlation function for
different Stokes numbers. However, before we explicitly carry out the model
calibration, we want to discuss the linear filter approximation in the context
of the Batchelor model for Lagrangian tracers in turbulence. The latter
modeling technique operates directly at the level of the second order
structure function $\mathcal{S}_{2}(\tau)=2\sigma^{2}-2C_{v}(\tau)$ and
differentiability of the Lagrangian velocity is ensured by introducing the
correct dissipation range behavior at small $\tau$. By contrast to the present
stochastic model, which was based on a Langevin equation (14) at its basic
level, the Batchelor model will not be discussed in form of a stochastic
process, but rather in the context of the so-called multifractal formalism
which operates directly on the structure function level [24].

##  IV Inertial particle statistics from Batchelor’s model for Lagrangian
tracers

In the seminal work of Batchelor [38] an interpolation between dissipation and
inertial range was proposed in the Eulerian frame of reference. The
interpolation was capable to reproduce the inertial range behavior of the
second order structure function $\langle(\delta_{r}u)\rangle\sim r^{2/3}$ and,
at small $r$, the dissipation range prediction
$\langle(\delta_{r}u)\rangle\sim r^{2}$ by a simple polynomial interpolation
[39]. Similarly, an interpolation for the Lagrangian second order structure
function $\left\langle(\delta_{\tau}v)^{2}\right\rangle$, which covers both
inertial $\left\langle(\delta_{\tau}v)^{2}\right\rangle\sim\tau$ and
dissipation range $\tau^{2}$, can be deduced [22, 25, 24, 40]. In its simplest
form, the second order Lagrangian structure function reads

$S_{2}(\tau)=\left\langle(\delta_{\tau}v)^{2}\right\rangle=2\sigma^{2}\frac{\frac{\tau}{T}}{\left[1+\left(\frac{\tau}{\tau_{\eta}}\right)^{-\delta}\right]^{\frac{1}{\delta}}}\;,$
(26)

where $\delta$ determines the transition between dissipation and inertial
range. Within the multifractal formalism, which proposes a more complete
modeling of the differential action of viscosity not only up to second order
statistics as it is done here, but also for higher order structure functions,
Lagrangian statistics seem to be well reproduced by choosing $\delta=4$.
Nevertheless, in what follows, where we completely neglect the implications of
intermittent corrections on the width and shape of the so-called intermediate
dissipative range, we will restrict ourselves to the case $\delta=2$, as
initially chosen by Batchelor to reproduce the behavior of the second order
structure function of Eulerian velocity, which leads to an acceptable
agreement with DNS data, as it will be shown in the following Section.
Furthermore, because the proposed parametrization for $S_{2}$ (equation 26 is
understood as a model for scales $\tau\ll T$, we need to come up with a proper
cutoff at large scales if we intend to use such a form for $\tau\gg T$).
Following the proposition of Ref. [21], we include in this simple picture an
additional exponential decrease over the characteristic large time scale $T$,
and consider the alternative form

$S_{2}(\tau)=2\sigma^{2}\frac{1-e^{-\frac{\tau}{T}}}{\left[1+\left(\frac{\tau}{\tau_{\eta}}\right)^{-\delta}\right]^{\frac{1}{\delta}}}\;,$
(27)

which entails a finite Lagrangian integral time scale. For our case, i.e.,
$\delta=2$, this time scale can be calculated explicitly

$T_{L}=\int_{0}^{\infty}\textrm{d}\tau\frac{\mathcal{C}_{v}(\tau)}{\mathcal{C}_{v}(0)}=\tau_{\eta}\left[1-\frac{\pi}{2}Y_{1}\left(\frac{\tau_{\eta}}{T}\right)-\frac{\pi}{2}H_{-1}\left(\frac{\tau_{\eta}}{T}\right)\right]\;,$
(28)

where $Y_{n}(z)$ denotes the Bessel function of the second kind and $H_{n}(z)$
the Struve function, as it is provided by a symbolic calculation software. In
this framework, the acceleration correlation function for the tracer particle
can be calculated from the second order structure function (26) by
$\mathcal{C}_{a}(\tau)=\frac{1}{2}\frac{\textrm{d}^{2}\mathcal{S}_{2}(\tau)}{\textrm{d}\tau^{2}}$.

Finite particle inertia can again be included in using the linear filter
approximation discussed in Section II. To this end, Eq. (10) has to be
evaluated from the Lagrangian velocity correlation
$\mathcal{C}_{v}(\tau)=\sigma^{2}-\frac{S_{2}(\tau)}{2}$. We could not obtain
a closed expression for the integral and thus have restricted ourselves to a
numerical evaluation, which will be further elaborated upon in the next
section.

##  V Comparison to DNS

In this section, we will apply the presented modeling techniques to the DNS
discussed in Section II.3. In order to connect these approaches to the
simulated data, parameters of the DNS must be defined to properly calibrate
the models, namely the integral length scale, $T$, the dissipative scale,
$\tau_{\eta}$, and the effective particle response time $\tau_{p}^{*}$.
Detailed explanations of the calibration technique can be found in [21]. In
summary, two conditions have to be imposed for the Lagrangian stochastic model
($\tau_{p}^{*}=0$), namely the matching of the Lagrangian time scale of the
stochastic model (25) to $T_{L}$ of the DNS as well as the matching of the
zero crossing of the acceleration correlation function (22) to the zero
crossing from the DNS data. This calibration suggests the values
$(\tau_{\eta}/\tau_{K},T/T_{L})=(0.5759,0.9791)$. For the Batchelor model, the
same calibration has to be carried out, where now the Lagrangian integral time
scale is given by Eq. (28). Furthermore, by matching the zero crossing, we
obtain the values $(\tau_{\eta}/\tau_{K},T/T_{L})=(1.7956,0.9941)$. It can be
noted that $\tau_{\eta}$ from the Batchelor model is larger than $\tau_{K}$
whereas the stochastic model exhibits a $\tau_{\eta}$ smaller than $\tau_{K}$
after calibration.

With the inclusion of St, an additional free parameter of the models is
available for calibration, the effective particle response time
$\tau_{p}^{*}$. As discussed, the acceleration correlation function for the
stochastic Gaussian process is given by Eq. (22), from which $\tau_{p}^{*}$
can be extracted based on the zero crossing of the DNS data for each Stokes
number. In a similar fashion, the new model parameter $\tau_{p}^{*}$ is
obtained from the linear filtering of the acceleration correlation function of
the Batchelor model derived from Eq. (27).

Figures 4(a) and 4(b) depict the comparison of second order structure function
$S_{2}(\tau)=\left\langle(\delta_{\tau}v_{p})^{2}\right\rangle$ stochastic and
Batchelor model, respectively (dash-dotted curves), to DNS for St=0-1 (solid
curves). For the tracers, St=0, the original stochastic model [21] and
Batchelor model are implemented.

For $\textrm{St}>$0, linear filtering of the velocity correlation function for
tracer particles agrees well with $S_{2}(\tau)$ from DNS for both models. The
stochastic approach, Figure 4(a) shows agreement between the model and DNS at
small scales and deviates slightly as the time lag $\tau$ increases. Notably,
slight deviations appear in the inertial range and might be attributed to the
Gaussianity of the stochastic model, which thus neglects intermittency
corrections. The application of the our filtering technique to the statistics
of the Batchelor model, Figure 4(b), shows similar tendencies. At small scales
the model coincides with the DNS profiles for the given St presented. At
$\tau/\tau_{K}>$1, again a deviation occurs where the model begins to
overestimate the structure function of its corresponding DNS curve, the near-
dissipative range seems to extend further than the one present in the
simulated data. These deviations slightly increase with St.

Figure 4: Comparison of the second-order structure function to DNS from (a)
the stochastic model (dash-dotted lines) and (b) the Batchelor model for 5
different St parameters in between 0 and 1. The structure functions have been
shifted vertically by multiplication of $10^{-n}$ for $n=0-4$ from top to
bottom.

The application of the models to the acceleration correlation function is
presented in Figure 5. Here, the discrepancy between modeled correlations and
those obtained from the simulated data can be observed at all scales. The
stochastic modeling of the inertial particle correlations is presented in
Figure 5(a) for the considered St range. The model over-predicts correlations
in the dissipative range. For increasing Stokes number, linear filtering of
the model improves the ability to accurately describe small scale correlations
and still the large scale variations between the two curves are minimal, for
example for St=0.2. As St increases further, quickly the filtered model and
the DNS show increased discrepancies between the profiles, as the
decorrelation of the acceleration occurs more rapidly than the model predicts.
Filtering of the Batchelor model, Figure 5(b), shows similar tendencies of the
results of the stochastic approach, but with even greater variation. The
correlation of tracer velocity at St=0 is slightly over-predicted at small
scales and slightly under-predicted at large scales. Comparable dissimilarity
is observed between the DNS and model at St=0.1. As the Stokes number
increases, the linear filter of the model breaks down and the predictions
decorrelate slower when compared to the simulated data responses.

Figure 5: Comparison of acceleration auto-correlation function from (a) the
stochastic model (dash-dotted lines) and (b) the Batchelor model (dash-dotted
lines) to DNS for 5 different St parameters in between 0 and 1.

Figure 6: (a) Zero crossing $\tau_{0}$ of the acceleration correlation
functions from DNS (as partially shown in Figure 2) as a function of St. The
inset shows a double-logarithmic representation with the black line indicating
$\tau_{0}\sim\textrm{St}^{0.35}$. (b) Recalibrated particle response time
$\tau_{p}^{*}$ expressed as $\textrm{St}^{*}=\tau_{p}^{*}/\tau_{K}$ based on
the zero crossing $\tau_{0}$ (see main text for more details on the
calibration process) as a function of St.

Due to the fact that the model calibration for finite St is based on the zero
crossing of the acceleration correlation, it is worth studying the St-
dependence of this quantity as well. Therefore, Figure 6(a) depicts the zero
crossing $\tau_{0}$ from DNS as a function of the Stokes number St. For the
case of Lagrangian tracers $\textrm{St}=0$, the zero crossing is at around
$2.2\tau_{K}$ and increases sharply for $\textrm{St}>0$. For
$\textrm{St}\approx 1$, a new quantitative behavior sets in, and the zero
crossing exhibits a slower increase. Furthermore, the zero crossing becomes
rather noisy, therefore, it is not entirely clear whether the zero crossing
would saturate at even higher St. The inset of Figure 6(a) shows a double-
logarithmic representation of $\tau_{0}$. For $\textrm{St}<1$, the zero
crossing appears to be a power law, whereas deviations from this power law
appear at $\textrm{St}\approx 1$. For comparison, the black line shows a power
law $\sim\textrm{St}^{0.35}$. However, at this point, no clear
phenomenological description could be provided that would allow for the
explanation of such a power law of the zeros of acceleration for inertial
particle motion.

Figure 6(b) depicts the calibrated effective Stokes number
$\textrm{St}^{*}=\frac{\tau_{p}^{*}}{\tau_{K}}$ as a function of the DNS
Stokes number St for both the stochastic (red) and the Batchelor model
(violet). Interestingly, the curves strongly resemble Figure 6(a) which
suggests a nearly linear relation between the zero crossing $\tau_{0}$ and the
effective (calibrated) particle response time $\tau_{p}^{*}$. The $a_{rms}$
values as a function of the Stokes number for all modeling techniques are
included in Figure 7 for direct comparison of the statistic. The filtered
acceleration, found directly from Eq. (12) (orange), and the stochastic model,
based on $\tau_{p}$ (green), quickly deviate from the DNS $a_{rms}$ curve
while the two models with the updated $\tau_{p}^{*}$ for the stochastic
process (red) and Batchelor representation (purple) show improved agreement at
all St and good agreement between $0.2\leq\textrm{St}\leq 1.1$. Cencini
$et\,al.$ [35] present a similar comparison and suggest that the deviation
between the profiles at small Stokes numbers is due to preferential
concentration captured which is not captured by the linear filter
approximation. The inclusion of the effective particle response time
$\tau_{p}^{*}$ in our model counteracts this discrepancy, providing accurate
representations of acceleration statistics for the presented Stokes numbers.

Figure 7: (a) Comparison of root mean square values of acceleration for DNS
and the different models. (b) Semi-logarithmic representation of (a).

##  VI Conclusion

We have presented a modeling technique for single inertial particle statistics
based on a filtering approach for the Lagrangian fluid velocity. The
introduction of an effective particle response time $\tau_{p}^{*}$ in the
linear filter is motivated by the strong discrepancies in the dynamics of
tracer-like and true inertial particles at low Stokes numbers. In particular,
the effective particle response time accounts for the effects of preferential
sampling of depleted vorticity regions by inertial particles, which manifests
itself, for instance, in a strong decrease of the variance of the particle
acceleration as suggested by Figure 7. Both the stochastic and the Batchelor
model are capable of reproducing this characteristic feature of inertial
particles. In contrast to the refinement of the linear filtering approximation
proposed by Deutsch and Simonin [31, 32], here, we do not operate on the level
of fluid quantities as “seen” by inertial particles, but rather modify the
“response of particles to fluid structures”. This has the advantage that the
calibration of our model equations is rather simple (i.e., it suffices to
determine the empirically observed zero crossing of the acceleration
correlation) in comparison to evaluating fluid quantities on the basis of
individual inertial particle trajectories. It will be a task for the future,
to further investigate the peculiar notion of preferential
sampling/preferential concentration by extending the modeling techniques and
possibly combine them with generalized vortex models similar to those in Refs.
[41, 42, 43, 44] . Furthermore - as we currently operate at the level of
second order statistics - much more could be understood by focusing for
instance on the behavior of the flatness of velocity increment which
highlights the implication of the existence of intermittent corrections. To
include this behavior, we would first need to derive the action of the linear
filtering at this higher level of statistics, and then to generalize both the
stochastic process as well as the parametrization of Batchelor to include
intermittent corrections, as it was done for Lagrangian velocity in Ref. [21].
Furthermore, the proposed model might be generalized to particle pairs which
opens up the way for investigating particle collisions [28] on the basis of
simple stochastic models or the multifractal formalism.

###### Acknowledgements.

J.F. acknowledges funding from the Humboldt Foundation within a Feodor-Lynen
fellowship and also benefited from financial support of the Project IDEXLYON
of the University of Lyon in the framework of the French program “Programme
Investissements d’Avenir” (ANR-16-IDEX-0005). B.V. and R.B.C. are supported by
U.S. National Science Foundation grant (NSF-GEO-1756259). R.B.C. is also
thankful for the support provided through the Fulbright Scholar Program.

## References

  * Monin and Yaglom [2007] A. S. Monin and A. M. Yaglom,  _Statistical Fluid Mechanics: Mechanics of Turbulence_ (Courier Dover Publications, 2007).
  * Frisch [1995] U. Frisch,  _Turbulence_ (Cambridge University Press, 1995).
  * Pope [2000] S. B. Pope,  _Turbulent Flows_ (Cambridge University Press, 2000).
  * Toschi and Bodenschatz [2009] F. Toschi and E. Bodenschatz, Lagrangian properties of particles in turbulence, Annu. Rev. Fluid Mech. 41, 375 (2009).
  * Machicoane  _et al._ [2019] N. Machicoane, P. D. Huck, A. Clark, A. Aliseda, R. Volk, and M. Bourgoin, Recent developments in particle tracking diagnostics for turbulence research, in  _Flowing Matter_ (Springer, 2019) pp. 177–209.
  * Qureshi  _et al._ [2008] N. M. Qureshi, U. Arrieta, C. Baudet, A. Cartellier, Y. Gagne, and M. Bourgoin, Acceleration statistics of inertial particles in turbulent flow, Eur. Phys. J. B 66, 531 (2008).
  * Ayyalasomayajula  _et al._ [2006] S. Ayyalasomayajula, A. Gylfason, L. R. Collins, E. Bodenschatz, and Z. Warhaft, Lagrangian measurements of inertial particle accelerations in grid generated wind tunnel turbulence, Phys. Rev. Lett. 97, 144507 (2006).
  * Squires and Eaton [1991] K. D. Squires and J. K. Eaton, Preferential concentration of particles by turbulence, Phys. Fluids 3, 1169 (1991).
  * Eaton and Fessler [1994] J. K. Eaton and J. R. Fessler, Preferential concentration of particles by turbulence, Int. J. Multiph. Flow 20, 169 (1994).
  * Monchaux  _et al._ [2012] R. Monchaux, M. Bourgoin, and A. Cartellier, Analyzing preferential concentration and clustering of inertial particles in turbulence, Int. J. Multiph. Flow 40, 1 (2012).
  * Bourgoin and Xu [2014] M. Bourgoin and H. Xu, Focus on dynamics of particles in turbulence, New J. Phys. 16, 085010 (2014).
  * Maxey and Riley [1983] M. R. Maxey and J. J. Riley, Equation of motion for a small rigid sphere in a nonuniform flow, Phys. Fluids 26, 883 (1983).
  * Chen  _et al._ [2006] L. Chen, S. Goto, and J. C. Vassilicos, Turbulent clustering of stagnation points and inertial particles, J. Fluid Mech. 553, 143 (2006).
  * Gustavsson and Mehlig [2016] K. Gustavsson and B. Mehlig, Statistical models for spatial patterns of heavy particles in turbulence, Advances in Physics 65, 1 (2016).
  * Bec  _et al._ [2006] J. Bec, L. Biferale, G. Boffetta, A. Celani, M. Cencini, A. Lanotte, S. Musacchio, and F. Toschi, Acceleration statistics of heavy particles in turbulence, J. Fluid Mech. 550, 349 (2006).
  * Bourgoin  _et al._ [2006] M. Bourgoin, N. T. Ouellette, H. Xu, J. Berg, and E. Bodenschatz, The role of pair dispersion in turbulent flow, Science 311, 835 (2006).
  * Taylor [1922] G. I. Taylor, Diffusion by Continuous Movements, Proc. London Math. Soc. s2-20, 196 (1922).
  * Obukhov [1959] A. Obukhov, Description of turbulence in terms of lagrangian variables (Elsevier, 1959) pp. 113–116.
  * Sawford [1991] B. L. Sawford, Reynolds number effects in Lagrangian stochastic models of turbulent dispersion, Phys. Fluids 3, 1577 (1991).
  * [20] J.-F. Pinton and B. L. Sawford, A lagrangian view of turbulent dispersion and mixing, in  _Ten Chapters in Turbulence_, edited by P. A. Davidson, Y. Kaneda, and K. R. Sreenivasan (Cambridge University Press) pp. 132–175.
  * Viggiano  _et al._ [2020] B. Viggiano, J. Friedrich, R. Volk, M. Bourgoin, R. B. Cal, and L. Chevillard, Modelling Lagrangian velocity and acceleration in turbulent flows as infinitely differentiable stochastic processes, J. Fluid Mech. 900, A27 (2020).
  * Chevillard  _et al._ [2003] L. Chevillard, S. Roux, E. Lévêque, N. Mordant, J.-F. Pinton, and A. Arnéodo, Lagrangian velocity statistics in turbulent flows: Effects of dissipation, Phys. Rev. Lett. 91, 214502 (2003).
  * Biferale  _et al._ [2004] L. Biferale, G. Boffetta, A. Celani, B. Devenish, A. Lanotte, and F. Toschi, Multifractal statistics of lagrangian velocity and acceleration in turbulence, Phys. Rev. Lett. 93, 064502 (2004).
  * Arnéodo et al. [2008] A. Arnéodo et al., Universal intermittent properties of particle trajectories in highly turbulent flows, Phys. Rev. Lett. 100, 254504 (2008).
  * Chevillard  _et al._ [2012] L. Chevillard, B. Castaing, A. Arneodo, E. Lévêque, J.-F. Pinton, and S. G. Roux, A phenomenological theory of Eulerian and Lagrangian velocity fluctuations in turbulent flows, C R Phys 13, 899 (2012).
  * Friedrich [2003] R. Friedrich, Statistics of lagrangian velocities in turbulent flows, Phys. Rev. Lett. 90, 84501 (2003).
  * Gatignol [1983] R. Gatignol, The faxen formulas for a rigid particle in an unsteady non-uniform Stokes flow, J. Mec. Theor. Appl, 1, 143 (1983).
  * Pumir and Wilkinson [2016a] A. Pumir and M. Wilkinson, Collisional Aggregation Due to Turbulence, Annu. Rev. Condens. Matter Phys. 7, 141 (2016a).
  * Tchen [2013] C.-M. Tchen,  _Mean value and correlation problems connected with the motion of small particles suspended in a turbulent fluid_ (Springer, 2013).
  * Hinze [1959] J. Hinze, Turbulence mcgraw-hill, New York 421, 3 (1959).
  * Deutsch [1992] E. Deutsch,  _Dispersion de particules dans une turbulence homogène isotrope stationnaire calculée par simulation numérique directe des grandes échelles_ , Ph.D. thesis, Ecully, Ecole centrale de Lyon (1992).
  * Deutsch and Simonin [1991] E. Deutsch and O. Simonin, Large eddy simulation applied to the modelling of particulate transport coefficients in turbulent two-phase flows, in  _8th Symposium on Turbulent Shear Flows, Volume 1_ , Vol. 1 (1991) pp. 10_1_1–10_1_6.
  * Février [2000] P. Février,  _Etude numérique des effets de concentration préférentielle et de corrélation spatiale entre vitesses de particules solides en turbulence homogene isotrope stationnaire_ , Ph.D. thesis, Toulouse, INPT (2000).
  * Berk and Coletti [2021] T. Berk and F. Coletti, Dynamics of small heavy particles in homogeneous turbulence: a lagrangian experimental study, J. Fluid Mech. 917 (2021).
  * Cencini  _et al._ [2006] M. Cencini, J. Bec, L. Biferale, G. Boffetta, A. Celani, A. S. Lanotte, S. Musacchio, and F. Toschi, Dynamics and statistics of heavy particles in turbulent flows, J. Turbul. , N36 (2006).
  * Yu  _et al._ [2012] H. Yu, K. Kanov, E. Perlman, J. Graham, E. Frederix, R. Burns, A. Szalay, G. Eyink, and C. Meneveau, Studying Lagrangian dynamics of turbulence using on-demand fluid particle tracking in a public turbulence database, J. Turbul. , N12 (2012).
  * Pumir and Wilkinson [2016b] A. Pumir and M. Wilkinson, Collisional Aggregation Due to Turbulence, Annu. Rev. Condens. Matter Phys. 7, 141 (2016b).
  * Batchelor [1951] G. Batchelor, Pressure fluctuations in isotropic turbulence, in  _Math. Proc. Cambridge Philos. Soc._ , Vol. 47 (Cambridge University Press, 1951) pp. 359–374.
  * Meneveau [1996] C. Meneveau, Transition between viscous and inertial-range scaling of turbulence structure functions, Phys. Rev. E 54, 3657 (1996).
  * Benzi  _et al._ [2010] R. Benzi, L. Biferale, R. Fisher, D. Q. Lamb, and F. Toschi, Inertial range Eulerian and Lagrangian statistics from numerical simulations of isotropic turbulence, J. Fluid Mech. 653, 221 (2010).
  * Friedrich and Friedrich [2013] J. Friedrich and R. Friedrich, Generalized vortex model for the inverse cascade of two-dimensional turbulence, Phys. Rev. E 88, 1 (2013).
  * Marcu  _et al._ [1995] B. Marcu, E. Meiburg, and P. K. Newton, Dynamics of heavy particles in a Burgers vortex, Phys. Fluids 7, 400 (1995).
  * Köhler  _et al._ [2016] J. Köhler, J. Friedrich, A. Ostendorf, and E. L. Gurevich, Characterization of azimuthal and radial velocity fields induced by rotors in flows with a low Reynolds number, Phys. Rev. E 93, 23108 (2016).
  * Ayyalasomayajula  _et al._ [2008] S. Ayyalasomayajula, Z. Warhaft, and L. R. Collins, Modeling inertial particle acceleration statistics in isotropic turbulence, Phys. Fluids 20, 95104 (2008).