\documentclass{jfm}
\usepackage{graphicx}
\usepackage{epstopdf, epsfig, bm}
\usepackage{color}
\usepackage{amsmath}
\bibliographystyle{jfm}
%\usepackage{pgfplots}
%\graphicspath{{./fig/}{./png/}}
%\graphicspath{{./fig/}}
\graphicspath{{fig/}}
\def\red{\textcolor{red}}
\def\blue{\textcolor{blue}}

%\graphicspath{{./fig/}}
%\graphicspath{{/Users/nilshau/repositories/github/christer/fig/}}

\newcommand\nn{\nonumber}
\newcommand{\tr}{\rm Tr\,}
\newcommand\ve[1]{\boldsymbol{#1}}
\newcommand{\ma}[1]{\ensuremath{\amsbb{#1}}}
%\newcommand{\Reys}{\ensuremath{\text{Re}_{\textem{s}}}}
%\newcommand{\Reyp}{\ensuremath{\text{Re}_{\textem{p}}}}
%\newcommand{\St}{\ensuremath{\text{St}}}
%\newcommand{\Kn}{\ensuremath{\text{Kn}}}
%\newcommand{\Ma}{\ensuremath{\text{Ma}}}
%AB: changed; is this good enough?
%NILS: Looks good to me
\newcommand{\Reys}{Re_{s}}
\newcommand{\Reyp}{Re_{p}}
\newcommand{\Reyn}{Re}
\newcommand{\St}{St}
\newcommand{\Kn}{Kn}
\newcommand{\Ma}{Ma}
\newcommand{\Sl}{\ensuremath{\text{Sl}}}
\newcommand{\rD}{{\rm D}}
\def\epsK{\epsilon_{\rm K}}
\def\cs{\ensuremath{c_{\text{s}}}}
\def\kf{\ensuremath{k_{\text{f}}}}
\def\urms{\ensuremath{u_{\text{rms}}}}
\def\taus{\ensuremath{\tau_{\text{s}}}}
\def\taut{\ensuremath{\tau_{\text{t}}}}
\def\tauk{\ensuremath{\tau_{\text{k}}}}
\def\tauforce{\ensuremath{\tau_{\text{force}}}}
\def\tauK{\ensuremath{\tau_{\text{K}}}}
\def\nap{\ensuremath{n_{\text{a}}}}
\def\zon{\ensuremath{\phantom{0}}}
\def\ztw{\ensuremath{\phantom{00}}}
\newcommand{\tens}[1]{\boldsymbol{\mathsf{#1}}}
\renewcommand{\vec}[1]{\ensuremath{\mbox{\boldmath$#1$}}}
\newcommand{\nab}{{\bm{\nabla}}}

%\newcommand{\Tab}[1]{Table~\ref{#1}}
\newcommand{\Sec}[1]{Section~\ref{#1}}
\newcommand{\App}[1]{Appendix~\ref{#1}}
%\newcommand{\Fig}[1]{Figure~\ref{#1}}
%\newcommand{\Figs}[2]{Figures~\ref{#1} and \ref{#2}}

%swap these lines in/out only for purpose of editing
%\newcommand{\Fig}[1]{\textcolor{red}{figure~\ref{#1}}}
%\newcommand{\FFig}[1]{\textcolor{red}{Figure~\ref{#1}}}
%\newcommand{\Figs}[2]{\textcolor{red}{figures~\ref{#1} and \ref{#2}}}
%\newcommand{\Figss}[2]{\textcolor{red}{figures~\ref{#1}--\ref{#2}}}
%\newcommand{\Figp}[2]{\textcolor{red}{figure~\ref{#1}({\it #2})}}
%\newcommand{\Figsp}[3]{\textcolor{red}{figures~\ref{#1}({\it #2}) and ({\it #3})}}
%\newcommand{\Figssp}[3]{\textcolor{red}{figures~\ref{#1}({\it #2})--({\it #3})}}
%\newcommand{\Tab}[1]{\textcolor{magenta}{table~\ref{#1}}}
%\newcommand{\Tabs}[2]{\textcolor{red}{tables~\ref{#1} and \ref{#2}}}

\newcommand{\Fig}[1]{figure~\ref{#1}}
\newcommand{\FFig}[1]{Figure~\ref{#1}}
\newcommand{\Figs}[2]{figures~\ref{#1} and \ref{#2}}
\newcommand{\Figss}[2]{figures~\ref{#1}--\ref{#2}}
\newcommand{\Figp}[2]{figure~\ref{#1}({#2})}
\newcommand{\Figsp}[3]{figures~\ref{#1}({#2}) and ({#3})}
\newcommand{\Figssp}[3]{figures~\ref{#1}({#2})--({#3})}
\newcommand{\Tab}[1]{table~\ref{#1}}
\newcommand{\Tabs}[2]{tables~\ref{#1} and \ref{#2}}

\newcommand{\EEq}[1]{Equation~(\ref{#1})}
\newcommand{\Eq}[1]{equation~(\ref{#1})}
\newcommand{\Eqs}[2]{equations~(\ref{#1}) and~(\ref{#2})}
\newcommand{\Eqss}[2]{equations~(\ref{#1})--(\ref{#2})}
\newcommand{\bra}[1]{\langle #1\rangle}
\newcommand{\meanrho}{\overline{\rho}}
\newcommand{\dd}{{\rm d}}
\newcommand{\um}{\,\mu{\rm m}}
\newcommand{\cm}{\,\rm cm}
\newcommand{\m}{\,\rm m}
\newcommand{\s}{\,\rm s}
\def\FF{\bm{F}}
\def\ff{\bm{f}}
\def\vv{\bm{v}}
\def\uu{\bm{u}}
\def\xx{\bm{x}}
\definecolor{awesome}{rgb}{1.0, 0.13, 0.32}
\newcommand{\obs}[1]{{\color{awesome}#1}}
\definecolor{antiquefuchsia}{rgb}{0.57, 0.36, 0.51}
\newcommand{\obss}[1]{{\color{antiquefuchsia}#1}}
\newcommand{\new}[1]{{\color{blue}#1}}
\newcommand{\ComN}[1]{{\color{red}NILS: #1}}

\def\ie{\textit{i.e. }}
\newcommand*{\ten}[1]{\underline{\underline{#1}}}
\newcommand{\overbar}[1]{\mkern 1.5mu\overline{\mkern-1.5mu#1\mkern-1.5mu}\mkern 1.5mu}

%\input{symdef.inc}

\newtheorem{lemma}{Lemma}
\newtheorem{corollary}{Corollary}

\shorttitle{Spectral characterisation of inertial particle clustering}
\shortauthor{N.E.L. Haugen, et al.}

\title{Spectral characterisation of inertial particle clustering in turbulence}

\author{
N.\ E.\ L.\ Haugen\aff{1,2},
A.\ Brandenburg\aff{1,3},
C.\ Sandin\aff{1} \and
L.\ Mattsson\aff{1}
 }

\affiliation{\aff{1}Nordita,
    KTH Royal Institute of Technology and Stockholm University,
    Hannes Alfv\'ens v\"ag 12, SE-10691 Stockholm, Sweden
\aff{2}SINTEF Energi A.S., Sem Saelands vei 11, 7034 Trondheim, Norway
\aff{3}Department of Astronomy, Stockholm University, SE-10691 Stockholm, Sweden
}

\newcommand\aap{Astron. Astrophys.}  % Astronomy and Astrophysics
\let\astap=\aap                      % alternative shortcut
\newcommand\aapr{A\&ARv}             % Astronomy and Astrophysics Review (the)
\newcommand\aaps{A\&AS}              % Astronomy and Astrophysics Supplement Series
\newcommand\actaa{Acta Astron.}      % Acta Astronomica
\newcommand\afz{Afz}                 % Astrofizika
\newcommand\aj{AJ}                   % Astronomical Journal (the)
\newcommand\ao{Appl. Opt.}           % Applied Optics
\let\applopt=\ao                         % alternative shortcut
\newcommand\aplett{Astrophys.~Lett.} % Astrophysics Letters
\newcommand\apj{Astrophys. J.}       % Astrophysical Journal
\newcommand\apjl{ApJ}                % Astrophysical Journal, Letters
\let\apjlett=\apjl                       % alternative shortcut
\newcommand\apjs{ApJS}               % Astrophysical Journal, Supplement
\let\apjsupp=\apjs                       % alternative shortcut
\newcommand\apss{Ap\&SS}             % Astrophysics and Space Science
\newcommand\araa{ARA\&A}             % Annual Review of Astronomy and Astrophysics
\newcommand\arep{Astron. Rep.}       % Astronomy Reports
\newcommand\aspc{ASP Conf. Ser.}     % ASP Conference Series
\newcommand\azh{Azh}                 % Astronomicheskii Zhurnal
\newcommand\baas{BAAS}               % Bulletin of the American Astronomical Society
\newcommand\bac{Bull. Astron. Inst. Czechoslovakia} % Bulletin of the Astronomical Institutes of Czechoslovakia 
\newcommand\bain{Bull. Astron. Inst. Netherlands} % Bulletin Astronomical Institute of the Netherlands
\newcommand\caa{Chinese Astron. Astrophys.} % Chinese Astronomy and Astrophysics
\newcommand\cjaa{Chinese J.~Astron. Astrophys.} % Chinese Journal of Astronomy and Astrophysics
\newcommand\fcp{Fundamentals Cosmic Phys.}  % Fundamentals of Cosmic Physics
\newcommand\gca{Geochimica Cosmochimica Acta}   % Geochimica Cosmochimica Acta
\newcommand\grl{Geophys. Res. Lett.} % Geophysics Research Letters
\newcommand\iaucirc{IAU~Circ.}       % IAU Cirulars
\newcommand\icarus{Icarus}           % Icarus
\newcommand\japa{J.~Astrophys. Astron.} % Journal of Astrophysics and Astronomy
\newcommand\jcap{J.~Cosmology Astropart. Phys.} % Journal of Cosmology and Astroparticle Physics
\newcommand\jcp{J.~Chem.~Phys.}      % Journal of Chemical Physics
\newcommand\jgr{J.~Geophys.~Res.}    % Journal of Geophysics Research
\newcommand\jqsrt{J.~Quant. Spectrosc. Radiative Transfer} % Journal of Quantitiative Spectroscopy and Radiative Transfer
\newcommand\jrasc{J.~R.~Astron. Soc. Canada} % Journal of the RAS of Canada
\newcommand\memras{Mem.~RAS}         % Memoirs of the RAS
\newcommand\memsai{Mem. Soc. Astron. Italiana} % Memoire della Societa Astronomica Italiana
\newcommand\mnassa{MNASSA}           % Monthly Notes of the Astronomical Society of Southern Africa
\newcommand\mnras{Mon. Not. Roy. Astron. Soc.}             % Monthly Notices of the Royal Astronomical Society
\newcommand\na{New~Astron.}          % New Astronomy
\newcommand\nar{New~Astron.~Rev.}    % New Astronomy Review
\newcommand\nat{Nature}              % Nature
\newcommand\nphysa{Nuclear Phys.~A}  % Nuclear Physics A
\newcommand\pra{Phys. Rev.~A}        % Physical Review A: General Physics
\newcommand\prb{Phys. Rev.~B}        % Physical Review B: Solid State
\newcommand\prc{Phys. Rev.~C}        % Physical Review C
\newcommand\prd{Phys. Rev.~D}        % Physical Review D
\newcommand\pre{Phys. Rev.~E}        % Physical Review E
\newcommand\prl{Phys. Rev.~Lett.}    % Physical Review Letters
\newcommand\pasa{Publ. Astron. Soc. Australia}  % Publications of the Astronomical Society of Australia
\newcommand\pasp{PASP}               % Publications of the Astronomical Society of the Pacific
\newcommand\pasj{PASJ}               % Publications of the Astronomical Society of Japan
\newcommand\physrep{Phys.~Rep.}      % Physics Reports
\newcommand\physscr{Phys.~Scr.}      % Physica Scripta
\newcommand\planss{Planet. Space~Sci.} % Planetary Space Science
\newcommand\procspie{Proc.~SPIE}     % Proceedings of the Society of Photo-Optical Instrumentation Engineers
\newcommand\rmxaa{Rev. Mex. Astron. Astrofis.} % Revista Mexicana de Astronomia y Astrofisica
\newcommand\qjras{QJRAS}             % Quarterly Journal of the RAS
\newcommand\sci{Science}             % Science
\newcommand\skytel{Sky \& Telesc.}   % Sky and Telescope
\newcommand\solphys{Sol.~Phys.}      % Solar Physics
\newcommand\sovast{Soviet~Ast.}      % Soviet Astronomy (aka Astronomy Reports)
\newcommand\ssr{Space Sci. Rev.}     % Space Science Reviews
\newcommand\zap{Z.~Astrophys.}       % Zeitschrift fuer Astrophysik

%\def\pcode{\textsc{pencil} code}
\def\pcode{\textsc{Pencil code}}

\begin{document}
\date{\today}

\maketitle

\begin{abstract}
Clustering of inertial particles is important for many types of
astrophysical and geophysical turbulence, but it has been studied
predominately for incompressible flows.
Here we study compressible flows and
compare clustering in both compressively (irrotationally) and
vortically (solenoidally) forced turbulence.
Vortically and compressively forced flows are driven stochastically either by
solenoidal waves or by circular expansion waves, respectively.
For compressively forced flows, the power spectra of the density of inertial
particles are a particularly sensitive tool for displaying particle
clustering relative to the density enhancement.
We use both Lagrangian and Eulerian descriptions for the particles.
Particle clustering through shock interaction is found to
be particularly prominent in turbulence driven by spherical
expansion waves.
It manifests itself through a double-peaked distribution of spectral
power as a function of Stokes number.
The two peaks are associated with two distinct clustering mechanisms;
shock interaction for smaller Stokes numbers and the centrifugal
sling effect for larger values.
The clustering of inertial particles is associated with the formation of caustics.
Such caustics can only be captured in the Lagrangian description,
which allows us to assess the relative importance of caustics
in vortically and irrotationally forced turbulence.
We show that the statistical noise resulting from the limited number
of particles in the Lagrangian description
can be removed from the particle power spectra,
allowing us a more detailed comparison of the residual spectra.
We focus on the Epstein drag law relevant for rarefied gases,
but show that our findings apply also to the usual Stokes drag.
\end{abstract}

\begin{keywords}
particle/fluid flow, compressible flows, shock waves, compressible turbulence
\end{keywords}

\section{Introduction}

Turbulent gas flows often harbour particles
of different sizes, for example dust grains or water droplets.
Larger particles can have significant inertia and decouple from the flow.
This can lead to local enhancements of their density, which is
generally described as clustering.
This can be important for the coalescence of particles to larger ones.
This process eventually leads to the formation of planetesimals in
planetary accretion discs \citep{1980Icar...44..172W}
or to raindrops in clouds \citep{Shaw03,GW13}. 
Another situation where particle clustering is important is when 
reactive particles are ``competing'' for the same reactant gas.
The concentration of reactant gas may then be significantly lower
within a particle cluster than outside, which yields an overall
lower reaction rate \citep{kruger17a,2018JFM...836..932H,karchniwy19}.

In clouds, and in most industrial applications, the compressibility of air
is rather weak, because the turbulent flow speeds are strongly subsonic.
This tends to be quite different in astrophysical flows such as those
in accretion discs around young stars and the interstellar medium,
which is continuously being fed with dust from the outflows
of stars near the end of their lives.
The driving of turbulence is fundamentally different in the
meteorological and astrophysical contexts.
Cloud turbulence is driven by convection,
which is an intrinsically vortical flow.
There is a large variety of turbulent industrial flows, but for 
the vast majority, the
turbulence is typically driven by some sort of shear, which yields
vortical flows.
Turbulent flows with purely compressive driving are sometimes also
referred to as acoustic (irrotational) turbulence or wave turbulence.
Turbulence in the interstellar medium is driven predominately by
supernova explosions, which are intrinsically compressive flows
\citep[see, e.g.,][]{1999ApJ...514L..99K,2004RvMP...76..125M,
2005A&A...436..585D,Federrath+11,Gent13a,Gent13b,Gent20,Evirgen19}.
At larger flow speeds, however, vorticity can always
be produced by baroclinicity and shocks 
\citep{Federrath+11,DSB11,Porter+15}.

To isolate the distinctive properties of compressive and vortical
turbulence, we simulate isothermal turbulence by applying a stochastic
forcing that is either vortical or compressive.
In the latter case, the pressure enhancements, which are the result of
energy injection of the forcing, are completely irrotational.
It is only when the resulting spherical expansion waves interact with
each other that some vorticity can be produced -- especially for
large Mach number, which is the ratio of root-mean-square (rms) velocity
to the sound speed.
Likewise, the purely vortical driving can lead to significant
compression and finite flow divergences at larger Mach numbers
\citep[see, e.g.,][]{Federrath+11,Mattsson19a}.

For acoustic turbulence, the energy spectra drop slightly more rapidly
with wavenumber $k$ (like $k^{-2}$) compared to the
$k^{-5/3}$ Kolmogorov spectrum for vortical turbulence \citep{KP73}.
Acoustic turbulence does not necessarily imply large Mach numbers,
because the bulk speed may well be less than the wave speed.
Our primary interest lies in the clustering of particles in these two
types of flows for small and large Mach numbers.

There are two rather different approaches to simulating
non-interacting inertial particles.
One is to model them as a pressure-less fluid, and the
other is to model them as non-interacting point particles.
In both cases, the particles interact with the gas through friction.
We refer to these two approaches as Eulerian
and Lagrangian, respectively.
Each of them has advantages and disadvantages.
A Lagrangian description is ideal for dilute systems,
but it is susceptible to statistical noise,
especially at small length scales where there are fewer particles.
This is a disadvantage of the Lagrangian approach.
An important disadvantage of the fluid description is that it cannot
describe how particles of the same type can go past each other.
This is because in the fluid description one only considers the bulk flow,
which is the average velocity of all particles of a specific type in a
small local volume.
The bulk velocity is therefore a single-valued
function of position.
In the Lagrangian description, by contrast,
one does not average, and since there are usually
several particles in every small volume,
the flow velocity of the particle phase can be multi-valued.
In particular, when particles go past each other,
we have the formation of caustics.
This implies that particles of the same size can
have opposite velocities at the same location,
creating phase-space singularities.
In the Eulerian approach, this situation would
lead to the formation of shocks -- even for dilute
particle populations.
In the Lagrangian approach, by contrast, particle
populations can go through each other without any
interaction.

Caustic formation can be an important pathway to enhanced particle
interaction \citep{WM05,2010JFM...645..497B,Gustavsson+2012,2014JFM...749..841V}.
This applies mainly to particles large enough to decouple from the carrier
fluid and this phenomenon can be the main reason why such particles interact.
At high Mach numbers, it is mainly shock interaction and compression of
the carrier fluid when shocks meet that matters.
This is because the density increase due to compression is by far the
greatest effect.

A more commonly studied route to enhanced interaction rates is the
centrifugal effect of turbulent eddies,
which fling the particles to the edges of the eddies \citep{ML83,Maxey87,
Squires91,Eaton94,Bec03,Bec07b,Bhatnagar18,Yavuz18}.
This is because the particles do not experience the
confining pressure that keeps the gas on closed streamlines.
The relative importance of caustics and the centrifugal
effects is not well understood \citep{2014JFM...749..841V}.
One of the goals of the present study is therefore
to assess their roles separately for vortical and
acoustic turbulence.
This assessment will be based on the knowledge that caustics are
only captured by the Lagrangian approach while the centrifugal effect
is captured both by the Lagrangian and the Eulerian approaches.

The Lagrangian and Eulerian descriptions are complementary
and can also be used to gauge their respective regimes of validity.
This allows us to study, for example, when and at what
length scales statistical noise becomes important,
and when caustics formation becomes important.
We focus on three-dimensional simulations, but we
also use one-dimensional simulations, where caustic formation
can be studied in isolation.

In both types of approaches, we ignore the back-reaction
of particles on the flow.
This can become important at large mass loading parameters
and can lead to other interesting effects such as the
streaming instability \citep{JY07J} and the resonant drag instability 
\citep{Squire17}, which will not be addressed here.
We also neglect gravity and tidal forces.

\section{The model}
\label{sec:model}

We consider an isothermal gas where the pressure is proportional
to the density $\rho$ and is given by $\rho\cs^2$, with $\cs$ being
the isothermal sound speed.
The velocity of the gas $\bm{u}$ is governed by the Navier-Stokes
and continuity equations,
\begin{equation}
\label{momentum}
\frac{\partial\uu}{\partial t}+\uu\cdot\bm{\nabla}\uu=
-\cs^2\bm{\nabla}\ln\rho + \bm{f}
+\rho^{-1}\bm{\nabla}\cdot(2\rho\nu\bm{\mathsf{S}}+\rho\zeta_{\rm shock}\bm{\mathsf{I}}\nab{\bm\cdot}\uu),
\end{equation}
\begin{equation}
\frac{\partial\ln\rho}{\partial t}+\uu\cdot\bm{\nabla}\ln\rho
=-\bm{\nabla}\cdot\uu,
\end{equation}
where $\bm{f}$ is a stochastic forcing term,
$\nu$ is the kinematic viscosity,
$\bm{\mathsf{I}}$ is the unit matrix with indices ${\sf I}_{ij}=\delta_{ij}$, and
$\bm{\mathsf{S}}$ is the trace-less rate of strain tensor with the components
\begin{equation}
\mathsf{S}_{ij}={\textstyle\frac{1}{2}}(\partial u_i/\partial x_j+\partial u_j/\partial x_i)
-{\textstyle\frac{1}{3}}\delta_{ij}\bm{\nabla}\cdot\bm{u}.
\end{equation}
The forcing term consists either of random plane waves (vortical forcing)
that are $\delta$-correlated in time
\citep[see][for details]{2012PhFl...24g5106H},
or of localised pressure enhancements (compressive forcing)
with $\bm{f}=-\bm{\nabla}\phi$,
where $\phi$ is a Gaussian in space at new locations
in regular time intervals $\delta t_{\rm f}$ \citep{2006MNRAS.370..415M}.

In most of the simulations presented here, we perform direct numerical
simulations (DNS).
In those cases, $\zeta_{\rm shock}=0$.
However, to save resources, especially in astrophysics, one often uses
a shock-capturing viscosity \citep{vNR50}.
This broadens the shocks and allows one to resolve them on a coarser mesh.
To assess the effect of such an artificial treatment on the particle
clustering, in some cases we compare the results from the DNS
with runs where a coarser mesh is used together with a shock-capturing
viscosity.
We adopt the shock capturing viscosity of \cite{vNR50}, which corresponds to
a bulk viscosity with
\begin{equation}
\zeta_{\rm shock}=C_{\rm shock} \delta x^2 \langle-\nab{\bm \cdot}\uu\rangle_+.
\label{zeta_shock}
\end{equation}
Here, $\langle...\rangle_+$ denotes a running five point average over all positive
%arguments, corresponding to a compression.
%AB: we could refer either to our 2005 paper, or perhaps to an earlier preprint by Aake.
%AB: or even Caunt-Korpi?
%AB: Now referred to this paper
%NILS: Good, I think it makes sense to include even our paper with
%NILS: Tony. I have done that now. OK for you?
arguments, corresponding to a compression; see \cite{2001A&A...369..706C, 10.1111/j.1365-2966.2004.08127.x}
for a detailed description.
%AB.
In contrast to the DNS, we refer to those simulations as large eddy
simulations (LES).

In both the Lagrangian and Eulerian descriptions,
the velocity $\vv_p$ of the particle with
index $p$ couples to the gas through the friction force,
\begin{equation}
\label{drag}
\FF_p=-\frac{1}{\tau_p}\left(\vv_p-\uu\right).
\end{equation}
For the flows considered here,
the particle response time is given by a term
that is slightly different for dense and dilute gases.
When the mean-free path of the gas molecules
is short, the response time is given by the
Stokes time, modified by a Reynolds number-dependent factor of the form
\begin{equation}
\label{tau_St}
\tau_p^{\text{St}}=\frac{2}{9}\, \frac{\rho_p}{\rho} \, \frac{a_p^2}{\nu}
\left(1+0.15\,\Reyp^{0.687}\right)^{-1},
\end{equation}
where $a_p$ and $\rho_p$ are the radius and material density, respectively,
and Re$_{p} = a_{p}u_{\rm rms}/\nu$ is the particle Reynolds number.
For rarefied gases, the response time is based on the Epstein drag and is given by
\citep{Schaaf63,Kwok75,Draine79b,2019MNRAS.490.5788M}
\begin{equation}
\label{tau_Ep}
\tau_p^{\text{Ep}}=\sqrt{\frac{\pi}{8}}\frac{\rho_p}
{\rho}\,\frac{a_p}{\cs}\left(1+\frac{9\pi}{128}
\frac{\left|\vec{u}-\vec{v}_p\right|^2}{c_{\text{s}}^2}\right)^{-1/2}.
\end{equation}
The second term inside of the parenthesis becomes important
at large Mach numbers.
In the following, we assume that the physical size of particles is
always smaller than the Kolmogorov size of the turbulence.

In the Lagrangian description, the evolution of a particle $p$
with velocity $\vv_p$ at position $\xx_p$ is given by
\begin{equation}
\label{dvdt_dxdt}
\frac{d\vv_p}{dt}=\FF_p,\quad
\frac{d\xx_p}{dt}=\vv_p,
\end{equation}
where we have ignored the effects of Brownian motion,
which leads to the diffusion of particles.
In the associated Eulerian description, diffusion is included and,
instead of \Eq{dvdt_dxdt}, we solve instead
\begin{equation}
\label{momentump}
\frac{\partial\vv_p}{\partial t}+\vv_p\cdot\bm{\nabla}\vv_p=\FF_p
+\frac{2}{\rho_p}\bm{\nabla}\cdot(\rho_p\nu_p\bm{\mathsf{S}}_p),
\end{equation}
\begin{equation}
\label{densityp}
\frac{\partial n_p}{\partial t}+\vv_p\cdot\bm{\nabla} n_p
=-n_p\bm{\nabla}\cdot\vv_p+\kappa_p\nabla^2 n_p,
\end{equation}
where $n_p$ is the particle number density and
$\nu_p$ and $\kappa_p$ are artificial viscosity and diffusivity
for the particle fluid (denoted now by $p$ collectively for all particles).
Those terms are needed for reasons of numerical stability.

%We use a medium-resolution ($N_{\rm mesh}^3=256^3$ and $512^3$) in
%three-dimensional (3D) triply periodic cubic domains with side lengths $L$
%and the number of mesh points in each direction is $N_{\rm mesh}$.
%NILS: Some small changes
We use medium-resolution ($N_{\rm mesh}^3=256^3$ and $512^3$) in
three-dimensional (3D) triply periodic cubic domains with side lengths $L$
and the number of mesh points in each direction being $N_{\rm mesh}$.
%NILS.
The smallest wavenumber in the domain is then $k_1=2\pi/L$.
Unless otherwise specified, dust particles are
included as inertial particles in 5 size bins with $4\times10^6$ particles
in each (for the Lagrangian simulations).
We use the \pcode\ \citep{2020_JOSS}, which is a
high-order, finite-difference code (sixth order in space and third
order in time); see also \citet{2002CoPhC.147..471B} for details. 

\section{Diagnostic tools}
\subsection{Non-dimensional numbers}

The flow is characterised by the Reynolds and Mach numbers,
\begin{equation}
\label{eq:Reynolds}
\Reyn=\frac{\urms}{\nu\kf}
\quad\mbox{and}\quad
\Ma=\frac{u_{\rm rms}}{c_{\rm s}},
\end{equation}
respectively, where $\kf$ is the forcing wavenumber.
The behaviour of the particles is characterised by the Stokes numbers
\begin{equation}
\label{eq:stokes}
\St_{\rm int}=\tau_{p}/\tau_{\rm f}
\quad\mbox{and}\quad
\St_{\rm Kol}=\tau_{p}/\tau_{\rm Kol},
\end{equation}
where $\tau_p$ is the particle response time, 
$\tau_{\rm f}=\kf/u_{\rm rms}$ is the time scale related to the forcing scale,
and $\tau_{\rm Kol}=\sqrt{\nu/\epsK}$ is the Kolmogorov time,
where $\epsK=\bra{2\rho\nu\bm{\mathsf{S}}^2}$ is the energy dissipation rate.
Both variants of the Stokes number are related to particle clustering
due to particle inertia.
For small Mach numbers, $\rho$ is close to the mean density $\meanrho$,
allowing us to express $\epsK$ also in terms of the energy spectrum $E(k)$.
They are normalised such that $\int E(k)\,\dd k=\meanrho\bra{\uu^2}/2$.
We then have $\epsK=2\meanrho\nu\int k^2 E(k)\,\dd k$.
Also of interest is the associated Kolmogorov scale
$\ell_{\rm Kol}=(\nu^3/\epsK)^{1/4}$.

The Knudsen number of particles of a certain size is defined as the ratio
of the mean free path $\lambda$
of the gas molecules to the size of the particle, i.e.,
$\Kn=\lambda/d_{\text{p}}$.
The drag force on the particles is inversely proportional to the
particle response time; see Eq.~(\ref{drag}).
For a continuous fluid, where the Knudsen number is much smaller than unity,
the particle response
time is given by the Stokesian time; see Eq.~(\ref{tau_St}).
For rarefied gases, however, the mean free path is large compared to the
particle size and the response time is then given by the
Epstein time, as given in Eq.~(\ref{tau_Ep}).

The dimensionless particle size parameter
\citep[see][]{2016MNRAS.456.4174H} is defined as
\begin{equation}
    \alpha = {\rho_{p}\over \langle\rho\rangle} {a_{p}\over L_{\rm f}},
\end{equation}
where $L_{\rm f}$ is taken to be the physical forcing scale of the
turbulent flow.
For small values of $\Kn$, we find that the mean Stokes number is
\begin{equation}
    \langle\St_{\rm int}\rangle = {\langle \tau_p^{\rm St}\rangle\over \tau_{\rm f}}
    \approx {(2/9)\,\alpha\,Re_p\over 1+0.15\,Re_p^{0.687}} \sim \alpha\,Re_p,
\end{equation}
while in the Epstein limit we find
\begin{equation}
    \langle\St_{\rm int}\rangle = {\langle\tau_p^{\rm Ep}\rangle\over \tau_{\rm f}} \approx \sqrt{\pi\over 8}\,\alpha\,\mathcal{M}_{\rm rms}\, \left(1+\frac{9\pi}{128}
\frac{\left|\vec{u}-\vec{v}_p\right|^2}{c_{\text{s}}^2}\right)^{-1/2}\sim \alpha\,\mathcal{M}_{\rm rms}.
\end{equation}
From these two relations one can see that, while $\St_{\rm int}$
for a given particle size in the Epstein limit is mainly affected by
compression and essentially unaffected by viscosity of the carrier fluid,
it is inversely proportional to the viscosity in the Stokes limit.

\subsection{Power spectra of particle density}
\label{PowerSpectra}

To measure preferential clustering at all scales, from the smallest
scale resolved in the simulation to the size of the simulation box,
we compute power spectra of $n_p$ as \citep{2018JFM...836..932H}
\begin{equation}
   P_{n}(k)=\frac{1}{2}\!\!\!\!\!\!\!\!\!\!\!\!
\sum_{k-\delta k/2 \leq |\vec{k}|<k+\delta k/2} \!\!\!\!\!\!\!\!\!\!\!\!
\left|\hat{n}_p(\vec{k})\right|^2 \text{d}^3\vec{k},
\end{equation}
where $\hat{n}_p(\vec{k})={\cal{F}} \left( n_p(\vec{x}) \right)$ is the
Fourier transform of $n_p$, $\vec{k}=(k_x,k_y,k_z)$ is the wavevector,
and the integration is over concentric shells in wave number space.
From the above, we see that
\begin{equation}
\label{eq:int_P}
\int_{k_1}^{k_{\rm max}} P_{n}(k)\,\dd k=\frac{1}{2}\langle n_p^2\rangle,
\end{equation}
where $k_{\max}=k_1 N_{\rm mesh}/2$ is the Nyquist wavenumber, and angle
brackets, $\langle...\rangle$, represent spatial averaging, which means
that $\langle n_p^2\rangle^{1/2}$ is the root-mean-squared particle
number density.

If the particles were randomly distributed, $|\hat{n}_p(\vec{k})|$
would be on average independent of $k=|\vec{k}|$.
In three dimensions, however, the shell integration introduces an
additional $k^2$ factor, so we expect
\begin{equation}
\label{eq:P_n}
    P_{n}(k)=A k^2,
\end{equation}
which can be combined with \Eq{eq:int_P} to yield
\begin{equation}
\label{eq:np}
  \frac{1}{2} \langle n_p^2\rangle
  = \int_{k_1}^{k_{\rm max}} A k^2\, \dd k
  = \frac{A}{3}(k_{\rm max}^3-k_1^3)
  \approx \frac{A}{3} k_{\rm max}^3,
\end{equation}
where we have assumed $k_1\ll k_{\max}$, and $A$ is a constant
that we shall be concerned with later in \Sec{SpectralAnalysis}.
We improve on this description further below, when we analyse
concrete examples.

\section{Results}

In this paper, we are concerned with the differences in
particle clustering between vortical and compressive forcings.
To convey an impression of this phenomenon, we begin by showing in
\Fig{fig:projections} the projected particle number densities of snapshots
from DNSs with vortical and compressive forcings for Stokes numbers
around unity.
Evidently, the visual impressions for the two types of flow are rather
different.
Even though the typical length scales of the forcings are similar
in the two cases, the clustering phenomenon is markedly different.
For vortical forcing, the overall contrast between minimum and maximum
particle concentrations is much smaller than for compressive forcing.
In the vortical case, the particle concentrations take a more filamentary
and perhaps sheet-like structure, while in the compressive case, the
particle concentrations are more spherical in shape.

%isothermal_expwave_512_f4.0_dtf1.0_5adjagr_tsc_nu0.02_nushock05N_128proc
\begin{figure}
\centering\includegraphics[width=0.49\textwidth]{fig/W2_new.jpg}
\includegraphics[width=0.49\textwidth]{fig/E2_new.jpg}
\caption{Left: Projected particle number density in a snapshot from a DNS
%with vortical forcing. Right: the same from a DNS with compressive
%AB: Lars, is V2 correct?
with vortical forcing (later referred to as Run~V2).
%Right: the same from a DNS with compressive forcing.
%AB: Lars, is C2 correct?
Right: the same from a DNS with compressive forcing (later referred to as Run~C2).
Both cases correspond the particle size showing the most
%clustering ($\alpha = $, St~$\sim 1$).
%AB: Lars, could you say more about this run? What is Re and Ma, for example.
%AB: Is this a run that we discuss in this paper (V2 and C2)?
%AB: Skip alpha (don't know the value yet) and added values that are hopefully right.
clustering ($\St_{\rm int}=0.31$ and $0.36$, respectively).
\label{fig:projections}}
\end{figure}

\subsection{One-dimensional simulations}

To illustrate the effects specific to compressive clustering, let us
consider first a one-dimensional shock model as an illustrative example.
Here, we also compare the Lagrangian particle simulations with the ones
in the Eulerian description.

\subsubsection{Applicability of the Eulerian approach}
\label{Applicability}

Caustics formation in the particle distribution, 
which is evident from the presence of multi-valued particle velocities, 
is a phenomenon that cannot be described with the Eulerian approach.
For small enough particle inertia, i.e., for small Stokes
numbers, the Eulerian and Lagrangian approaches should agree with
each other.
However, there is also another source of discrepancy.
The Lagrangian approach is suitable for modelling dilute systems, but,
due to the finite number of particles,
it also has the disadvantage of suffering from statistical fluctuations when
the intention is to model non-dilute systems, where fluctuations should
be small.
Statistical noise does not occur in the Eulerian approach.
Thus, for dense systems and small Stokes numbers, the Eulerian approach
can be beneficial.
To determine the limits of applicability of the Eulerian approach
quantitatively, we consider a simple one-dimensional model.

We adopt a localised hump in the fluid density, which we model
by a Gaussian in $\ln\rho(x)$ at the position $x=0$ in a domain
of size $-\pi/2<x<3\pi/2$.
Thus, we take
\begin{equation}
    \ln\rho[(x)/\rho_0]={\cal A} \exp \left(-x^2/2\sigma_s^2 \right),
\end{equation}
where ${\cal A}$ is an amplitude factor,
$\rho_0$ is an overall normalisation coefficient,
the width is given by $\sigma=0.35$, and the ratio of 
the peak value over the background is 3.1.
For these simulations we use periodic boundary conditions.
The initial density profile launches an acoustic wave; see
\Fig{fig:pcomp_1d_initial} for plots of $u_x$ and $\rho$ at $t=0.1$
and $t=0.5$.
The front speed exceeds the sound speed when the gas speed approaches
a certain fraction of the sound speed; see Appendix~\label{pposition}
for an illustration.
We then use the gas velocity and gas density at $t=0.5$ as initial
condition for the particles by setting the particle velocity for all
particle sizes equal to the fluid velocity.
The particle number density of all particle sizes is set proportional
to the gas density.

\begin{figure}
\centering
\includegraphics[width=\textwidth]{fig/pcomp_1d_initial.eps}
\caption{
({\em a}) Gas velocity and ({\em b}) density at $t=0.1$ (upper row)
and $t=0.5$ (lower row) in panels ({\em c}) and ({\em d}).
The data for $t=0.5$ serve as initial condition for the particles.
\label{fig:pcomp_1d_initial}}
\end{figure}

\begin{figure}
\centering
\includegraphics[width=\textwidth]{fig/pcomp_1d_time.eps}
\caption{
Particle velocity in the Lagrangian simulation (black) and the Eulerian
one (red), together with the gas velocity (blue) at times $t=1$, 2,
%and 3 for $a_p=3\times10^{-3}$ (corresponding to $\St=1$, panels a, c, and e)
%AB: can/should we call it here St_int? We haven't defined anything else yet.
%NILS: I wonder if it is not better to just call it St (without
%NILS: any subscript), since St_int is defined based on the forcing
%NILS: scale, which is not defined in this 1D case.
%NILS: We could even look in the idl script to see exactly how
%NILS: we define the Stokes number here and then state the definition.
%NILS: What do you think?
and 3 for $a_p=3\times10^{-3}$ (corresponding to $\St_{\rm int}=1$, panels a, c, and e)
%and $10^{-1}$ (corresponding to $\St=30$, panels b, d, and f).
%AB: St_int
and $10^{-1}$ (corresponding to $\St_{\rm int}=30$, panels b, d, and f).
\label{fig:pcomp_1d_time}}
\end{figure}

In \Fig{fig:pcomp_1d_time} we show the particle velocities as a function
of position for Lagrangian and Eulerian models for particle radii
$a_p=3\times10^{-3}$ and $10^{-1}$ at different times.
We also show the gas velocity, which propagates at a speed slightly faster
than the sound speed, $\cs=1$; see \App{pposition} for a plot showing
the numerically obtained dependence of the front speed on the gas speed.
The lighter particles follow the fluid and are not shown, but the heavier
ones lag behind because they only inherit the speed of the gas at $t=0.5$,
and they are too heavy to get accelerated by the passing acoustic wave.
At early times ($t=1$), the velocities in the Lagrangian and Eulerian
models are close to each other, but at later times the particle velocities
in the Lagrangian description become multi-valued, which corresponds to
caustics formation.
This phenomenon becomes more prominent for the heavier particles
since they are not decelerated by the drag from the gas.
In the Eulerian description, we instead see the formation of a shock.
Away from the shock, the Lagrangian and Eulerian descriptions agree with
each other rather well, especially for the heavier particles.
To resolve the shock in the Eulerian simulations,
we must apply a certain amount
of artificial viscosity and diffusivity for the particle fluid.
If this artificial viscosity is too small, wiggles occur in the
downstream part of the shock, as can already be seen from the
profile of $v_x(x)$ in \Fig{fig:pcomp_1d_time}.
Including such artificial viscosity and diffusivity is a purely numerical
device to stabilise the solution, but it is likely to introduce errors
in the results.
By comparing with the Lagrangian approach, we will try to assess
the extent of such artifacts.

\begin{figure}
\centering
\includegraphics[width=\textwidth]{fig/pcomp_1d_dens.eps}
\caption{
Particle density in the Lagrangian simulation
(black) and the Eulerian one (red), together with
the gas density (blue) at times $t=1$, 2, and 3.
Two particle sizes are shown: $a_p=3\times10^{-3}$
(left column) and $10^{-1}$ (right column).
\label{fig:pcomp_1d_dens}}
\end{figure}

Snapshots of the particle number densities are shown in
\Fig{fig:pcomp_1d_dens} for the same times and the same two particle
sizes as in \Fig{fig:pcomp_1d_time}.
We see the development of extended structures
with two enhancements on their flanks, characteristic of caustics,
as is correctly reproduced with the Lagrangian approach.
The Eulerian approach, on the other hand, yields just a single albeit
very strong spike, 
which may cause an increase of the particle-interaction rate even without any caustics forming.

It is of interest to determine the Stokes number relevant for
caustics formation.
This is important for knowing the maximum Stokes number for which the
Eulerian approximation can still be used.
For smaller Stokes numbers,
no artificial particle viscosity and diffusivity are needed.
For larger Stokes numbers, however, the Eulerian approach can represent
the caustics only as shocks, which requires an increasing amount of
artificial viscosity and diffusivity to keep them numerically resolved.
The Stokes number is defined through \Eq{eq:stokes}.
For the solution shown in \Fig{fig:pcomp_1d_time}, we find that
the fluid travel time across the width of the front $\Delta x$ is
\begin{equation}
\tau_{\rm f}=\Delta x/\Delta u
\end{equation}
where $\Delta u$ 
is the fluid velocity at the front.
For the current experiments, $\Delta x$ is taken to be
the thickness of the front at the time 
when the particles were introduced in the simulation;
see Run~A in \Tab{tab:Stokes} for more details about the particles.

In the example discussed above, we have $\tau_{\rm f}=0.72/0.66\approx1.1$,
and we find caustics for $a_{p}\ge10^{-3}$,
which corresponds to $p=3$.
To compute the critical Stokes number $\St_{\rm crit}$
above which caustics formation occurs,
we need to know the $\rho_{p}/\rho$ ratio, where we take the fluid density
on the upstream side of the front, which is here $\rho\approx1.8$.
We used $\rho_{p}=10^3$, so, altogether, we have $\St_{\rm crit}=0.30$.
This means that, for simulations where the Stokes number is
larger than $\sim 0.3$, the Eulerian particle approach cannot be used.
It is important to realise that we are here talking about the Stokes
number based on the largest fluid scale, $\St_{\rm int}$, and not
the Kolmogorov scale.
As we shall see below, at the numerical resolutions accessible in our
three-dimensional simulations, the difference between $\St_{\rm int}$ and
$\St_{\rm Kol}$ is not very large.

\begin{table}%[]
    \centering
    \begin{tabular}{lcccccccc}
    & $p$       & 1        &       2          &      3    &      4          &     5      &       6          & 7 \\
Run & $k_2 a_p$ & $10^{-4}$ & $3\times10^{-4}$  &  $10^{-3}$& $3\times10^{-3}$ &  $10^{-2}$ &  $3\times10^{-2}$ &  $10^{-1}$ \\
    \hline            
A & $\St_p$     & $0.03 $   & $0.1$            &  $0.3$   & $1.0$             &  $3$       &  $10$          &  $30$ \\
B & $\St_p$     & $0.007$   & $0.02$           &  $0.07$  & $0.2$             &  $0.7$     &  $2.0$         &  $6.8$ \\
    \end{tabular}
    \caption{
    Stokes numbers for Runs~A (\Sec{Applicability}) and B (\Sec{sec:1D_B}).
}
    \label{tab:Stokes}
\end{table}

\subsubsection{A mechanism for irrotational supersonic clustering}
\label{sec:1D_B}

It is conceivable that irrotational supersonic flows can yield new ways
of clustering that would not occur in vortical subsonic flows.
One idea for such a mechanism is that particles of a suitable mass that
move toward each other on two colliding shocks, will be decelerated as the
shocks collide, but the particles are too heavy to become re-accelerated
as the shocks depart again immediately after the collision.
The particles will then be left behind after the shocks move away,
forming a cluster.
To test this idea, we use an experiment similar to that described in
\Sec{Applicability}, but with a stronger density enhancement.

In \Sec{Applicability}, we used the gas velocity and density at a certain
time to reinitialise the particles by setting the particle velocity for
all particle sizes to a value equal to the fluid velocity.
The initial width of the density distribution is again $0.35$.
However, the density enhancement of the gas is now so strong that its
distribution is so different from a Gaussian that it can no longer
be used for reinitialising the particles in a simple way.
We therefore reinitialise the fluid density equal to the density
of the lightest particles and then set the density and velocity of the
heavier particles to the same as the lightest particles.
The ratio of the peak value of the density over its background then
turns out to be 22.
In this experiment, we only use the Lagrangian approach.
See Run~B in \Tab{tab:Stokes} for more details about the particles.

\begin{figure}
\centering
\includegraphics[width=\columnwidth]{fig/pbut.eps}
\caption{$xt$ diagrams of ({\em a}) $\rho$, ({\em b}) $n_1$, ({\em c}) $n_4$,
({\em d}) $n_5$, ({\em e}) $n_6$, and ({\em f}) $n_7$.
}\label{fig:pbut}
\end{figure}

\begin{figure}
\centering
\includegraphics[width=\columnwidth]{fig/pprelvel}
\caption{
Velocity of particles with radius $a_6$ (black dots) and fluid velocity
(red lines), as well as the gas density (blue lines and axes on the right)
at $t=1.18$ (a), close to the time $t_*=1.16$ when the shocks
meet and the gas density develops a peak.
Panels (b) and (c) show the same at $t=t_*+\tau_6=1.22$ and
$t_*+2\tau_6=1.28$.
}\label{fig:pprelvel}
\end{figure}

The result of this experiment is shown in \Fig{fig:pbut}, where we plot
$xt$ diagrams of $\rho$ and $n_p$ for different particle sizes.
For small particles with small Stokes numbers, the particles follow
the gas.
This can be seen by the similarity between panels (a) and (b).
%To estimate $\St$ in this case, we used $\Delta u=2.6$ and
%AB: St_int
To estimate $\St_{\rm int}$ in this case, we used $\Delta u=2.6$ and
$\Delta x=0.08$, so $\tau_{\rm f}=0.08/2.6\approx0.03$.
For our largest particles with $a_{p}=0.1$, we find $\rho\approx300$;
see the blue lines in \Fig{fig:pprelvel}.
We again used $\rho_{p}=10^3$, which yields
$\rho_{p}/\rho=1000/300\approx3$.
Therefore, $\tau_{p}=\sqrt{\pi/8}\,(\rho_{p}/\rho)(a_{p}/\cs)\approx0.2$,
%so we have $\St\sim 7$.
%AB: St_int
so we have $\St_{\rm int}\sim 7$.

%For our largest particle size with $\St\approx7$,
%AB: St_int
For our largest particle size with $\St_{\rm int}\approx7$,
the two counter-streaming particle clouds associated with the two
opposing shocks tend to run through each other owing to their large inertia,
as we can see from \Figp{fig:pbut}{f}.
%For the intermediate size, where $\St\approx2$, however, a sizeable
%AB: St_int
For the intermediate size, where $\St_{\rm int}\approx2$, however, a sizeable
particle cloud is left behind at the original position of the collision
of the two shocks; see \Figp{fig:pbut}{e}.
This critical value is close to unity, as one might have expected.
We conclude from this that a particle cluster will
form after the collision of two shocks if the Stokes number is
around unity.
%AB: Lars to verify
Here, the Stokes number is based on the width and speed of the shock fronts.
This mechanism for particle clustering is fundamentally different from
the classical eddy mechanism of Maxey \& Riley and it operates only
for large enough Mach numbers.
For such flows, however, it may be the dominating mechanism.

%The biggest uncertainty in our estimate of $\St$ lies in the value of $\rho$.
%AB: St_int
The biggest uncertainty in our estimate of $\St_{\rm int}$ lies in the value of $\rho$.
In the following, we focus on the particles in bin~6,
and denote the corresponding velocity by $v_6$ and
the response time of those particles by $\tau_6$.
To determine the relevant value of $\rho$, we show in \Fig{fig:pprelvel} the
profiles of $\rho(x)$ together with those of $v_6(x)$ and $u(x)$ at times
$t_*+\tau_6/3$, $t_*+\tau_6$, and $t_*+2\tau_6$, where $t_*=1.16$ is the time
when the shocks meet and the gas density develops a peak.
We see that the peak in $\rho$ reaches values of around 400, but at the time
$t_*+\tau_6$, the particles will have slowed down considerably.
Therefore, the relevant density to be used is the temporally averaged
density at the peak until that time, which is below 300.

\subsection{Three-dimensional simulations}

%eulerian
%hel512f05a_nud1em3 crashed
%exp512f4_nud1em2var2 also crashed

%V1   isothermal_helical_256_f0.5_dtf1.0_5adjagr_tsc_nu1_noshocknu/"
%V2a  isothermal_helical_256_f4.0_dtf1.0_5adjagr_tsc_nu0.02_noshocknu'  %(t=262)
%V2b  isothermal_helical_512_f4.0_dtf1.0_5adjagr_tsc_nu0.005_noshocknu_512proc %(t=26)
%     isothermal_helical_512_f4.0_dtf1.0_5adjagr_tsc_nu0.005_noshocknu_512proc
%V2c  isothermal_helical_512_f4.0_dtf1.0_5adjagr_tsc_nu0.003_noshocknu_512proc %(crashed)
%V2_  isothermal_helical_512_f4.0_dtf1.0_5adjagr_tsc_nu0.01_noshocknu_512proc %(????)
%V2   isothermal_helical_256_f4.0_dtf1.0_5adjagr_tsc_nu0.01_noshocknu/"
%V3   isothermal_helical_512_f10_dtf1.0_5adjagr_tsc_nu0.05_noshocknu_512proc %(used to be called V3a)
%V3a  isothermal_helical_512_f10_dtf1.0_5adjagr_tsc_nu0.02_noshocknu_512proc %(used to be called V3)
%C1   isothermal_expwave_256_f0.5_dtf0.5_5adjagr_tsc_nu5e-3_noshocknu/"
%C1a  isothermal_expwave_512_f0.5_dtf0.5_5adjagr_tsc_nu2e-3_noshocknu
%     isothermal_expwave_256_f0.5_dtf0.5_5adjagr_tsc_nu2e-3_noshocknu
%C2   isothermal_expwave_256_f4.0_dtf1.0_5adjagr_tsc_nu0.05_noshocknu/" -> 0.02 at 512
%C2a  isothermal_expwave_512_f4.0_dtf1.0_5adjagr_tsc_nu0.02_noshocknu

%print,eta,dx: 0.0423458     0.0490874   (V1)
%print,eta,dx: 0.079597353   0.0245437 (V3)

% /cfs/klemming/nobackup/c/csandin/pencil-code/christer/isoth_expwave/isothermal_helical_256_f0.5_dtf1.0_5adjagr_tsc_nu1_noshocknu

\begin{table}%[]
    \centering
    \resizebox{\hsize}{!}{
    \begin{tabular}{ccllllcclccc}
    % 1      2        3       4        5                     6        7     8     9     10      11      12
     Case    &forcing& $f_0$&$\delta t_{\rm f}$& $\nu$ & $\Ma$ & $\Reyn$ &
     %$\tau_{\rm kin}$ & $\tau_{\rm Kol}$& $\St_{\rm int}$ & $\St_{\rm Kol}$ &
     %$\alpha$\\
%AB: don't have alpha yet
     $\tau_{\rm kin}$ & $\tau_{\rm Kol}$& $\St_{\rm int}$ & $\St_{\rm Kol}$ \\
    %$\alpha$\\
%AB: Nils, do we want to add the alpha values here?
%NILS: I dont believe we use alpha anywhere later in the paper, but
%NILS: since we mention alpha in section 3, we could do that.
%NILS: What do you prefer?
%NILS: Will it fit in the table?
\\
    \hline            
    V1  & vort.& 0.02 & --- & 0.001 & 0.15 &103 & 4.3 & 1.8  & $7.3\times10^{-4}$... 7.3 & $1.7\times10^{-3}$... 17 & \\
    V1S & vort.& 0.02 & --- & 0.001 & 0.15 &103 & 4.4 & 1.8  & $5.1\times10^{-4}$... 5.1 & $1.2\times10^{-3}$... 12 & \\    
    V2  & vort.& 0.2  & --- & 0.01~ & 0.67 & 45 & 1.0 & 0.56 & $3.1\times10^{-3}$... 31  & $5.5\times10^{-3}$... 55 & \\
    V2a & vort.& 0.2  & --- & 0.02~ & 0.59 & 20 & 1.13 & 0.81 & $2.8\times10^{-3}$... 28 & $3.9\times10^{-3}$... 39 & \\ % 
    V2b & vort.& 0.2  & --- & 0.005 & 0.71 & 95 & 0.94 & 0.39 & $3.3\times10^{-3}$... 33 & $8.0\times10^{-3}$... 80 & \\ %
    V3  & vort.& 0.5  & --- & 0.05  & 1.00 & 13 & 0.67 & 0.50 & $4.7\times10^{-3}$... 47 & $6.3\times10^{-3}$... 63 & \\ %used to be called V3a
    V3a & vort.& 0.5  & --- & 0.02  & 1.14 & 38 & 0.59 & 0.32 & $5.3\times10^{-3}$... 53 & $9.9\times10^{-3}$... 99 & \\ %used to be called V3
    C1  & comp. & 0.5  & 0.5 & 0.005 & 0.19 & 25 & 3.53& 1.89 & $8.9\times10^{-4}$... 8.9 & $1.7\times10^{-3}$... 17 & \\
    C1a & comp. & 0.5  & 0.5 & 0.002 & 0.19 & 64 & 3.47 & 1.22 & $9.0\times10^{-4}$... 9.0 & $2.6\times10^{-3}$... 26 & \\
    C1.5& comp. & 1.5  & 1.0 & 0.015 & 0.39 & 17 & 1.72 & 1.16 & $2.4\times10^{-2}$... 54  & $3.6\times10^{-2}$... 80 & \\
    C2  & comp. & 4.0  & 1.0 & 0.02~ & 0.76 & 25 & 0.88& 0.57 & $3.6\times10^{-3}$... 36  & $5.5\times10^{-3}$... 55 & \\
    \end{tabular}
    }
    \caption{Summary of our simulations; ``comp'' and ``vort'' refer to
    compressive and vortical forcings, respectively.
    \newline
    \label{tab:models}}
\end{table}
%dir="../isothermal_helical_256_f0.5_dtf1.0_5adjagr_tsc_nu1_noshocknu/"
%dir="../isothermal_helical_256_f4.0_dtf1.0_5adjagr_tsc_nu0.01_noshocknu/"
%dir="../isothermal_expwave_256_f0.5_dtf0.5_5adjagr_tsc_nu5e-3_noshocknu/"
%dir="../isothermal_expwave_256_f4.0_dtf1.0_5adjagr_tsc_nu0.05_noshocknu/"

In the previous sections, we studied several aspects of particle clustering in 
idealised one-dimensional simulations.
We will now proceed by turning our attention to fully three-dimensional turbulence
simulations.
As described in \Sec{sec:model}, two kinds of forcings are employed in this work.
In \Tab{tab:models}, run names starting with ``V'' use vortical forcings,
while those starting with ``C'' adopt spherical expansion wave forcing
(compressive forcing).
The numbers behind those letters indicate different forcing strengths,
which yield different Mach numbers.
Different Reynolds numbers are indicated by letters a and b.
The run with Stokes drag is denoted by the letter S at the end.
We also compare with corresponding Eulerian models
(\Tab{tab:models_eulerian}), where we have included
an artificial viscosity and diffusivity in order to
stabilise the simulations; see \Eqs{momentump}{densityp}.

% V1 5e-3 (hel256f05f variable)
% V2 1e-2 (hel512f05a_nud1em2)
% C1 1e-2 (exp512f05a_nud1em2 -> exp512f05a_nud1em2var [success])
% C2 5e-2 (exp256f4_a5em6 -> exp512f4_nud1em2var [crashed] -> exp512f4_nud1em2var2 [also crashed]
\begin{table}%[]
    \centering
    \begin{tabular}{ccllllccllll}
    % 1      2        3       4        5                     6        7     8     9     10      11      12
     Case    &forcing& $f_0$&$\delta t_{\rm f}$& $\nu$ & $\nu_{p}$ & $\Ma$ & $\Reyn$ & $N_{\rm mesh}$ \\
    \hline            
    V1  & vort.& 0.02 & --- & 0.001 & 0.005 & 0.15 &103 & 256 \\%hel256f05f
    V2  & vort.& 0.2  & --- & 0.01~ & 0.01~ & 0.67 & 45 & 512 \\%hel512f05a_nud1em2
    C1  & exp. & 0.5  & 0.05& 0.002 & 0.001 & 0.19 & 96 & 512 \\%exp512f05a_nud1em3b
    C2  & exp. & 4.0  & 1.0 & 0.05~ & 0.05~ & 0.72 & 14 & 512 \\%exp512f4_nud5em3
    \end{tabular}
    \caption{Summary of Eulerian runs.
    For all these runs, the artificial diffusivity ($\kappa_p$) equals the artificial 
    viscosity ($\nu_p$).
    \newline
    \label{tab:models_eulerian}}
\end{table}

\begin{figure}
\centering
\includegraphics[width=.8\textwidth]{fig/pnp_ap.eps}
\caption{
Contour plots of particle number density for (a) $\St_{\rm int}=$ 0.033,
(b) 0.33, (c) 3.3, and (d) 33 for case V2b. 
The particle number density has been integrated over the perpendicular
direction for four mesh zones.
\label{fig:pnp_ap}}
\end{figure}

\begin{figure}
\centering
\includegraphics[width=\textwidth]{fig/ppcoarse.eps}
\caption{Scatter plot of particle number density as a function of
fluid density for the three smallest particle sizes of case V2b. 
\label{fig:ppcoarse}}
\end{figure}

\begin{figure}
\centering
\includegraphics[width=.6\textwidth]{fig/fluid_spectra.eps}
\caption{$E(k)$ (solid lines) and $P_\rho(k)$ (dashed lines)
  for our main cases (V1, V2, C1, and C2).
\label{fig:fluid_spectra}}
\end{figure}

Contour plots of particle number density are shown in \Fig{fig:pnp_ap}
for Run~V2b.
It is clearly seen that the clustering is strongest for our intermediate
Stokes numbers, $\St_{\rm int}=0.33$ and $3.3$; see \Figsp{fig:pnp_ap}{b}{c}.
For the very smallest and largest Stokes numbers, we see almost no
clustering; see \Figsp{fig:pnp_ap}{a}{d}.

In \Fig{fig:ppcoarse}, we show scatter plots of the particle number density as 
a function of fluid density for all fluid grid cells in 
the domain. In the left panel, 
showing results for the smallest Stokes numbers, we see that there is
a strong correlation between the two. This is because these small 
particles follow the fluid almost perfectly, which means that,
when the fluid is compressed (high fluid density), the particle field is 
also compressed (high particle number density).
When the Stokes number is increased, but is still rather small, we see
in \Figp{fig:ppcoarse}{b} that the two fields are only weakly correlated.
Finally, for those Stokes numbers where we see the strongest clustering
in \Figp{fig:pnp_ap}{b}, we show in \Figp{fig:ppcoarse}{c} that there
is no correlation between particle and fluid densities.
We can also see the effect of the inertial clustering itself in that there is a large
number of grid cells without any particles ($n_p=0$), while there is
also a significant number of grid cells containing many particles ($n_p > 2$),
which is not the case for the smaller Stokes numbers.

\subsubsection{Kinetic energy and density power spectra}
\label{SpectralAnalysis}

In this paper, we make extensive use of particle power spectra.
In this context, it is useful to show first the relevant spectra for the gas.
In \Fig{fig:fluid_spectra}, we show kinetic energy spectra together with
power spectra of the gas density for cases V1, V2, C1, and C2.
The peak of the kinetic energy spectrum occurs at $k/k_1=3$ for the cases
with vortical forcings (V1 and V2).
This wavenumber indicates where the main power
of the external forcing is found.
For the cases with compressive forcing, however, the peak
in the kinetic energy spectrum is found for $k/k_1$ between 1 and 2.
The density power spectra follow the kinetic energy spectra fairly well,
although there is a vertical shift.
For the compressive forcing, we are driving strong flow
divergencies, which results in $P_\rho(k)$ being large compared to
$E_{\rm K}(k)$.
Since the vortical forcing is divergence-free, the corresponding density variation is small, as
seen through the smaller values of $P_\rho(k)$ compared to $E_{\rm K}(k)$.
We can also see that the extent of the vertical shift is larger for
smaller Mach numbers.
This is because the fluid is less compressed for smaller Mach numbers.

\subsubsection{Initial and tracer particle spectra}
\label{InitialAndTracer}

In the following, we study power spectra of particle number densities for particles
that are embedded in a fluid where turbulence is generated either with
vortical or compressive forcing (Runs~V1--V3 and Runs~C1--C2,
respectively).
Since particles are tracked in a Lagrangian fashion,
it is convenient to allocate each Lagrangian
particle to the nearest Eulerian neighbours in each direction in the fluid mesh.
In this way, for every particle size, we generate
a variable on the fluid mesh that contains the number of particles
that reside within or in the neighbourhood of a given grid cell.
These variables can now be used to calculate the particle power 
spectra for each particle size.
The size of the neighbourhood of grid cells that will be influenced by
a given particle depends on the interpolation scheme used, which
%in our case is a second order linear scheme.
%AB: we didn't refer to Johansen. Would this be the right paper?
%NILS: I can't find any reference to the interpolation scheme in 
%NILS: their Nature paper. I do not know which one that would be most 
%NILS: relevant. Since we are just talking about plain interpolation
%NILS: here I wonder if we actually need a reference?
in our case is a second order linear scheme \citep{2007Natur.448.1022J}.

Initially, the particles are randomly distributed
over the entire simulation box. This means that the initial power spectra
for the different particle sizes are just white noise, which corresponds to
a $k^2$ spectrum; see \Sec{PowerSpectra}.
If every particle is associated solely with the very nearest grid point
of the Eulerian mesh, the $k^2$ scaling would then be valid for the full
wavenumber range.
For the current work, however, the contribution from a particle is distributed
over several nearby grid points through a linear interpolation scheme.
This means that the $k^2$ scaling will not extend all
the way to the largest wavenumbers.
Instead, for $k>k_\ast$, the spectrum, hereafter $P_{{n},{\rm noise}}$,
becomes less steep and eventually reaches a maximum, before it goes
down towards the very end.

In analogy with \Eq{eq:np}, we now get
\begin{equation}
  \frac{1}{2} \langle n_p^2\rangle
  = \int_{k_1}^{{k}_{\rm max}} P_{{n},{\rm noise}}(k)\,\dd k
  = \int_{k_1}^{{k}_{\rm max}} A\tilde{k}^2\,\dd k
  \equiv \frac{A}{3}\tilde{k}_{\rm eff}^3,
\end{equation}
where $\tilde{k}=k/[1+(k/k_\ast)^3]$ has been substituted for $k$ in order 
to account for the departure from $k$
for $k\gtrsim k_\ast\equiv\kappa k_{\max}$, with $\kappa\approx0.789$ being
a parameter proportional to the position of the local maximum of $P_n(k)$.
This implies that $P_{{n},{\rm noise}}=A\tilde{k}^2$.
Defining therefore $\kappa=k_\ast/k_{\rm max}$, we find
$\tilde{k}_{\rm eff}^3=
(1+\kappa)(1-\kappa+\kappa^2)/\kappa^3k_{\rm max}^3\approx {(1.45\,k_{\max})^3}$.
By re-arranging the above equation, the constant $A$,
defined in \Sec{PowerSpectra}, is found to be
$A=3(\langle n_p^2\rangle/2)/\tilde{k}_{\rm eff}^3
\approx0.49\langle n_p^2\rangle/k_{\max}^3$.
Together with \Eq{eq:P_n} we then
obtain the initial power spectrum of the particles as
\begin{equation}
  \label{E_np}
   P_{{n},{\rm noise}}\approx
   0.49\,\frac{\langle n_p^2\rangle B^2}{k_{\rm max}^3}\tilde{k}^2,
\end{equation}
where the constant on the right-hand side is defined as
$B=\langle \rho\rangle/\langle n_p\rangle$ and is introduced
to compensate for the fact that the fluid and particle 
density fields do not have the same mean value.
This compensation is required in order to obtain \Eq{eq:P_model}.

Particles that are very small, having essentially vanishing Stokes
numbers, will follow the gas perfectly.
If the fluid is incompressible, particles will be re-shuffled owing
to turbulence, but their mean separation will be unchanged with time,
which means that there is no particle clustering.
Hence, the power spectrum of tracer particles in an incompressible
fluid would equal $P_{{n},{\rm noise}}$
for all times. If, however, the fluid is compressible,
the compression of the fluid may be so strong that the resulting fluctuations
in particle number density becomes larger than the white noise.
The particle power spectrum will then be the same as the one for the fluid density.
For tracer particles we therefore expect the power spectrum of the particles
to be given by
\begin{equation}
\label{eq:P_model}
  P_{n,{\rm model}}(k) = P_\rho(k) + P_{{n},{\rm noise}}(k).
\end{equation}
In \Fig{fig:power_tracer}({\it a}), we compare the calculated particle
spectra from the simulations $P_n(k)$ (solid lines) with the modelled
spectra $P_{n,{\rm model}}(k)$ (dashed lines) for the four runs V1, V2,
C1, and C2 for the lightest particles.
We see that the calculated and modelled spectra are remarkably similar
for all cases.
The $k^2$ part of the spectrum is a real physical effect and is a
result of having a finite number of particles in the simulation, such
a feature cannot be seen if the particles are tracked by the Eulerian
approach. Instead, one would then see that the particle spectra follow
the fluid density spectrum for all wavenumbers.
Unless the particle suspension consists of an infinite number of
particles, this behaviour is incorrect and is due to the fact that the
Eulerian approach treats the particle suspension as a continuous fluid
and not as a collection of a finite number of discrete particles.

\begin{figure}
\centering
\includegraphics[width=.49\textwidth]{fig/tracer.eps}
\includegraphics[width=.49\textwidth]{fig/tracer_E1.eps}
\caption{
Left:
Power spectra of particle number
density for the smallest Stokes numbers, which are essentially tracer
particles (solid lines) for cases V1, V2, C1, and C2.
The dashed lines denote the model spectrum as presented in 
\Eq{eq:P_model}.
Right:
Comparison of numerical power spectra (solid) and model
spectra using \Eq{eq:P_model} (dashed) for the smallest particles of
Run~C1 with $N_p=2.5\times10^6$ (blue),
$20\times10^6$ (red), and $160\times10^6$ (black).
\label{fig:power_tracer}}
\end{figure}

In \Fig{fig:power_tracer}({\it b}), we show the numerical particle spectra for 
the particles with the smallest Stokes number of
Run~C1, obtained with three different numbers of particles (solid lines).
The numerical results are then compared with the model data given by
\Eq{eq:P_model}. We see that the model results reproduce the simulation 
results for all particle numbers.
But, more importantly, we note that the effect of the finite number of 
particles becomes less dominant when the number of particles is increased, 
which is as expected. 
For the black line ($N_p=160\times 10^6)$,
10 times more particles than fluid mesh points were used. 
This highlights the difficulty in exploring weak clustering at large wavenumbers 
for large Reynolds numbers.

\subsubsection{Comparison of cases with Epstein and Stokesian drag}

For rarefied gases, the molecular mean free path is large compared to the
particle size, and the response time is then given by the
Epstein time, as presented in Eq.~(\ref{tau_Ep}).
For a dense (continuous) fluid, however, where $\Kn<1$,
the response time is given by the modified
Stokes time; see Eq.~(\ref{tau_St}).
To compare the effect of using these two response times,
we show in \Fig{fig:Ep_vs_Stokes}
power spectra of particle number densities for simulations using
Epstein and Stokes drags for approximate Stokes numbers between
$6\times 10^{-4}$ and 6.
We see that the two reflect rather similar trends.
At least part of the remaining discrepancies can probably be 
explained by
small differences in the actual Stokes numbers for the two drag laws.
Both for small and large Stokes numbers, we see a $k^2$ spectrum,
indicative of random particle distributions at high wavenumbers 
(small scales).
Initially, particles of all sizes were randomly distributed,
and hence they had a $k^2$ spectrum.
However, the reason that we still see a $k^2$ spectrum at later
times are different for the smallest and largest particles.
The largest particles are so heavy that their response times are far
too long for the particles to be able to react to the turbulent eddies
associated with the smallest spatial and temporal scales.
This is why the power spectra of the heavier particles
still show their initial $k^2$ spectrum at small scales.
The smallest particles, by contrast, are being re-shuffled by the turbulence
and therefore maintain a random particle distribution.
They do not have enough inertia to move from one fluid element to another,
which is a requirement for the particles to form inertia-based clusters.
They are fundamentally different from the short-lived
increase in the particle number density, which occurs always when the fluid
volume in which the tracer particles reside, is compressed -- for
example due to the passing of an acoustic wave or shock.

\begin{figure}
\centering
\includegraphics[width=.49\textwidth]{fig/Ep_vs_Stokes.eps}
\includegraphics[width=.49\textwidth]{fig/Ep_vs_Stokes_noise.eps}
\caption{
Comparison of power spectra of particle number densities 
for Run~V1 using Epstein (solid lines) and Stokes drag (dashed
lines) for $\St_{\rm int}=6\times10^{-4}$ (black), $6\times10^{-3}$ (red),
$6\times10^{-2}$ (orange), $0.6$ (green), and $6$
(blue) for (a) $P_n(k)$ and (b) $P_n(k)-P_{\rm noise}$.
\label{fig:Ep_vs_Stokes}}
\end{figure}

\subsubsection{Applicability of the Eulerian approach for particles}
\label{Applicability_Eulerian}

%pencil-code/christer/isoth_expwave/eulerian/idl
\begin{figure}\begin{center}
\includegraphics[width=.49\columnwidth]{fig/pcomp_spec_christer_V1}
\includegraphics[width=.49\columnwidth]{fig/pcomp_spec_christer_E1a_new}
\includegraphics[width=.49\columnwidth]{fig/pcomp_spec_christer_V2}
\includegraphics[width=.49\columnwidth]{fig/pcomp_spec_christer_E2}
\end{center}\caption[]{
Comparison of power spectra of particle number densities for
Runs~V1 ({\em a}), C1 ({\em b}), V2 ({\em c}), and C2 ({\em d}).
}\label{pcomp_spec_christer_helexp}\end{figure}

In \Fig{pcomp_spec_christer_helexp}, particle power spectra for
different Stokes numbers are shown for our four main runs, V1, V2, C1, and C2.
The solid lines correspond to results obtained with the Lagrangian approach
for the particles, while the Eulerian particle approach was used for
the simulations visualised by the dotted lines.
For V1 and V2, as shown in \Fig{pcomp_spec_christer_helexp}({\it a})
and ({\it c}), we see that the two approaches yield similar results for
the smaller Stokes number up to the wavenumbers where the $k^2$ scaling
commences in the Lagrangian simulations, which corresponds to $k \sim 4
k_1$ for the $5\times4$ million particles chosen for these simulations
with five different radii.
When the Stokes number is increased (different colours), the Lagrangian
and Eulerian spectra for $k/k_1 \leq 4$ are still comparable for V1,
but for V2 we see a clear difference for the largest Stokes numbers,
and not just for the largest wavenumbers.

We recall that in the 1-D simulations, the critical value of
$\St_{\rm int}$ for the occurrence of caustics was around 0.3;
see the discussion at the end of \Sec{Applicability_Eulerian}.
In the present case of 3-D turbulence, we also see a difference
between the Eulerian and Lagrangian simulations near
$\St_{\rm int}=0.3$, but only for the spherical expansion waves;
see \Figp{pcomp_spec_christer_helexp}{d}.
When the turbulence is driven compressively, as in Runs~C1 and C2, which are
shown in \Figsp{pcomp_spec_christer_helexp}{b}{d}, we see that 
the Eulerian and Lagrangian approaches yield comparable results only for the very 
smallest Stokes numbers.
Thus, the range of applicability of the Eulerian approach depends on
the type of forcing and is more restrictive in the compressive case.
For the other cases, the differences are rather small, except perhaps
for the heaviest particles.
This difference could also be caused by our usage of artificial
viscosity and diffusivity in the Eulerian simulations.
At large wavenumbers, on the other hand, there is always a large
difference, but this is mainly caused by the effect of noise.

For Run~C2 in \Figp{pcomp_spec_christer_helexp}{d}, we also show 
the power spectra from the Eulerian approach with the contribution
from $P_{{n},{\rm noise}}$ of \Eq{E_np} being added.
This models a power spectrum that is accounting for the noise from a
finite number of Lagrangian particles; see the dashed lines.
For the smallest Stokes number, we see that the solid and dashed
lines almost overlap, but for all other Stokes numbers the resemblance 
is poor. This means that, except for the very smallest Stokes
numbers, the differences between the particle power spectra obtained
with the Lagrangian and Eulerian particle approaches are \emph{not}
primarily due to the noise contribution of the Lagrangian approach.
Furthermore, from the results shown in \Fig{pcomp_spec_christer_helexp},
we can also conclude that the Eulerian approach should \emph{not} be
used to track particles unless the Stokes number is low.

As expected from the discussion above (\Sec{InitialAndTracer}),
we notice a $k^2$ behaviour for the smallest Stokes numbers; see
\Fig{pcomp_spec_christer_helexp}.
This applies, of course, only to the Lagrangian simulations
with a finite number of particles, and cannot be seen in the
corresponding Eulerian simulations.
However, there is a clear departure from the $k^2$ behaviour as the Stokes
number is increased.
The reason is that particle inertia now starts to have an effect
on the clustering. This clustering is not due to fluid compression,
but rather due to other inertia-based clustering 
mechanisms, such as the Maxey--Riley mechanism.

For intermediate Stokes numbers, the particle power spectra resemble
those presented in \cite{2018JFM...836..932H}.
In that paper, the authors speculate that the peak in the particle
power spectra is associated with the similarity in the characteristic
time scales of the turbulence and the response time of the particles.
However, such a connection could not be confirmed in their work owing
to their limited Reynolds number.
We see here the same trend as found by Haugen et al., namely that
the individual maxima of the spectra are insensitive to the Stokes
number.
%AB: could we say more here? It sounds like a boring end...

\subsubsection{Reynolds number dependence}

In order to investigate the nature of the inertia-based clustering
further, we would like to run simulations with much larger Reynolds numbers.
In the DNS, this becomes very costly when the Mach number is large and
shocks need to be resolved.
In \Fig{pcomp_spec_christer_Re_hel}, we show the power spectra for
different Stokes and Reynolds numbers.
For the smaller Stokes numbers, an increase in $\Reyn$ leads to an
increase in spectral power.
For the larger Stokes numbers, the trend is opposite.
Looking at the Kolmogorov-based Stokes number ($\St_{\rm Kol}$), however, we see that,
%AB: in many places (especially in the rest), we just write \St.
%AB: We could either add "int" in those cases (if that's what we mean)
%AB: or we could say clearly which one we mean. The former is better.
for a given value of $\Reyn$, we get more power and hence more clustering
when the Kolmogorov-based Stokes number is closer to unity.
For the Reynolds numbers obtained here, it therefore seems that it is the 
Kolmogorov based Stokes number that controls the strength of the
clustering, not the one based on the integral scale ($\St_{\rm int}$).
The same has also been found in other low Reynolds number studies \citep{Bec07b,Baker2017}.
As the Reynolds number is increased, however, one eventually
%reaches a point where particles with $\St$ around unity will be much slower than
%AB: St_int, right?
reaches a point where particles with $\St_{\rm int}$ around unity will be much slower than
the smallest turbulent eddies and they will therefore be totally decoupled from 
the Kolmogorov scale. Hence, the clustering cannot be determined by
$\St_{\rm Kol}$ in such cases.
The nature of this large-scale particle clustering still remains to
be understood, because much larger resolution would be needed.

%christer_dir1='isothermal_helical_256_f4.0_dtf1.0_5adjagr_tsc_nu0.01_noshocknu' ;(used in other Fig)
%christer_dir2='isothermal_helical_256_f4.0_dtf1.0_5adjagr_tsc_nu0.02_noshocknu'
\begin{figure}\begin{center}
\includegraphics[width=.49\columnwidth]{fig/pcomp_spec_christer_Re_hel}
\includegraphics[width=.49\columnwidth]{fig/pcomp_spec_christer_Re_noise_hel}
\end{center}\caption[]{
    Comparison of power spectra of particle number densities 
    for Run~V2 (dashed lines) with Run~V2a (dotted lines) and
    Run~V2b (solid lines) for (a) $P_n(k)$ and (b) $P_n(k)-P_{\rm noise}(k)$.
}\label{pcomp_spec_christer_Re_hel}\end{figure}

\begin{figure}\begin{center}
\includegraphics[width=.99\columnwidth]{fig/non-mono-spec.eps}
%AB: in this figure, we should use italics St and add "int"
%AB: Also the charactersin the legends are bigger than in Fig.15.
%AB: Also the subscript 1 on k should be roman.
\end{center}\caption[]{
$P_n(k)$ and $P_n(St_{\rm int})$ for
Run~C1 with more particle sizes.
}\label{non-mono-spec}\end{figure}

\subsubsection{Clustering mechanisms}

From Run~C1, shown in \Figp{pcomp_spec_christer_helexp}{b}, we see that
%the spectra for $\St=0.09$ and 0.9 are very similar.
%AB: St_int, right?
the spectra for $\St_{\rm int}=0.09$ and 0.9 are very similar.
This may be an indication of a non-monotonic behaviour of the
spectral evolution with Stokes number.
In order to investigate this further, we perform a new simulation
that is identical to Run~C1 except that we now include many more
closely spaced particle sizes (Stokes numbers).
The spectra for some of these Stokes numbers are shown in the left-hand
panel of \Fig{non-mono-spec}.
%From here we see that the spectral power increases from $\St=0.018$
%AB: here -> this, also added _{\rm int}
From this we see that the spectral power increases from $\St_{\rm int}=0.018$
%up until $\St=0.3$, before it decreases again as we move towards $\St=1.5$.
%AB: do we mean here St_int?
up until $\St_{\rm int}=0.3$, before it decreases again as we move towards $\St_{\rm int}=1.5$.
Finally, a clear increase is seen for larger Stokes numbers.
In the right-hand panel of \Fig{non-mono-spec} we plot the power for
five different wavenumbers as a function of the particle Stokes number
in order to see this non-monotonic behaviour more clearly.
Here we see that, for all wavenumbers shown, 
%the power reach two distinct
%AB: "reach" -> "spectra attain "
the power spectra attain two distinct
%maxima: the first is around Stokes number 0.3 while the other is 
%AB: abbreviation; should we add "int"?
maxima: the first is around $\St_{\rm int}=0.3$, while the other is 
%found between Stokes numbers 10 and 30.
%AB: abbreviation; should we add "int"?
found between $\St_{\rm int}=10$ and 30.
%It is not clear what is causing this non-monotonic behaviour,
%AB: added "immediately", because we now explain this, right?
It is not immediately clear what is causing this non-monotonic behaviour,
%but it could be due to a change in the relative importance
%AB: "it could be" -> "we argue here that it is"
but we argue here that it is due to a change in the relative importance
between two different particle clustering mechanisms.
Since compressive forcing was used for this run, the two competing
mechanisms are most likely (i) the shock-collision mechanism, as described
in Section~\ref{sec:1D_B}, and (ii) the classical Maxey--Riley
clustering mechanism.

If one of the peaks is due to the shock-collision mechanism, we would
expect it to be stronger as the Mach number is increased.
We therefore perform an intermediate simulation (Run~C1.5) where we
increase the Mach number to $\Ma=0.39$ to investigate this.
The results are shown in \Fig{non-mono-spec_E1p5}, where it is clearly
seen that the first peak has become substantially stronger, and also moved
somewhat to the right.
The second peak is, however, almost unchanged, except for a smaller
shift to the left.
This seems to indicate that it is the first peak that is due to the
shock-collision mechanism.
We now proceed by increasing the Mach number even further (Run~C2) and show the
results in \Fig{non-mono-spec_E2}. The first peak has now become so strong that
%the second peak is only visible as a weak shoulder for $\St\approx10$.
%AB: add "int"?
the second peak is only visible as a weak shoulder for $\St_{\rm int}\approx10$.
We will now continue by investigating the mechanism behind 
the second peak.

%FIG15
\begin{figure}\begin{center}
\includegraphics[width=.99\columnwidth]{fig/non-mono-spec_E1p5.eps}
\end{center}\caption[]{
Similar to \Fig{non-mono-spec}, but for
Run~C1.5 with more particle sizes.
}\label{non-mono-spec_E1p5}\end{figure}

%FIG16
\begin{figure}\begin{center}
\includegraphics[width=.99\columnwidth]{fig/non-mono-spec_E2.eps}
\end{center}\caption[]{
Similar to \Fig{non-mono-spec}, but for
Run~C2 with more particle sizes.
}\label{non-mono-spec_E2}\end{figure}

For the classical eddy mechanism of Maxey \& Riley, we expect the
clustering to depend primarily on the vortical part of the velocity field.
Hence, the Stokes number for the second peak should be
of the order of unity if the fluid time scale is calculated
based on the vortical part of the velocity field.
To check this, we have performed a Helmholtz decomposition of
the velocity field by computing the vector and scalar potential of
$\uu=\nab\times\psi+\nab\phi$.
The rms values of the corresponding velocity fields are given in
\Tab{tab:Helmholtz}.
We also list the estimated peak Stokes numbers from the simulations
at $k/k_1=32$ and those based on the vortical velocity field.
We see that the vortical peak Stokes numbers are 1.2 and 1.6 for Runs~C1
and C1.5, while for Run~C2, we only see an indication of a shoulder at
$\St_{\rm peak2}^{\rm vort}=1.6$, at least for smaller values of $k$.
This could be taken as an indication that the second peak is indeed
due to the classical eddy mechanism of Maxey \& Riley.

\begin{table}%[]
    \centering
    \begin{tabular}{llccccccc}
Run & $\urms$ & $\urms^{\rm vort}$ & $\urms^{\rm pot}$
& $\St_{\rm peak1}$ & $\St_{\rm peak2}$
& $\St_{\rm peak1}^{\rm vort}$ & $\St_{\rm peak2}^{\rm vort}$ \\
    \hline            
C1  &  0.170 & 0.017 & 0.169 & 0.15& 12 & 0.02 & 1.2 \\%E1_256_nrad15_more
C1.5&  0.352 & 0.11  & 0.334 & 0.6 &  5 & 0.19 & 1.6 \\%E1p5_256_nrad20
C2  &  0.865 & 0.200 & 0.842 &  2  & (7)& 0.46 &(1.6)\\%E2_256_nrad20
    \end{tabular}
    \caption{
Rms velocities for the Helmholtz decomposed velocities
together with the estimated peak Stokes numbers from the simulations
and those based on the vortical velocity field.
The numbers in parentheses are more uncertain. The reason for this is that,
for Run~C2, the second peak appears as a shoulder only, and no clear maximum can be identified.
}
    \label{tab:Helmholtz}
\end{table}


\begin{figure}
\centering
%/cfs/klemming/projects/snic/snic2020-4-12/brandenb/christer/isoth_expwave/E1p5_256_nrad20
\includegraphics[width=\textwidth]{fig/pnp_ap_E1p5.eps}
\caption{
Similar to \Fig{fig:pnp_ap}, but for Run~C1.5 showing
contour plots of particle number density for $\St_{\rm int}$
in the range from 0.3 to 57.
\label{fig:pnp_ap_E1p5}}
\end{figure}

%In \Fig{fig:pnp_ap_E1p5}, we have plotted the particle concentrations
%AB: have plotted -> show 
In \Fig{fig:pnp_ap_E1p5}, we show the particle concentrations
%for Run~C1.5 in order to investigate if we can see any trace of the two
%AB: investigate -> determine 
for Run~C1.5 in order to determine if we can see any trace of the two
different mechanisms behind the particle clustering.
We see that clustering is now apparent for a very broad range of Stokes
%numbers, ranging from $\St=0.3$ to $50$.
%AB: add "int"?
numbers, ranging from $\St_{\rm int}=0.3$ to $50$.
%Both for small and for large values of $\St$ do we see blob-like clusters,
%AB: skip second "for", add "int"?
Both for small and large values of $\St_{\rm int}$ do we see blob-like clusters,
while for intermediate values the structures are more sheet-like.
%Apart from that, there is no real difference in the morphology of 
%AB: Apart from -> Other than 
Other than that, there is no real difference in the morphology of 
structures between small and large Stokes numbers.

%In \Fig{non-mono-spec_V1} we show the same results as in
%AB: modified
In \Fig{non-mono-spec_V1} we present similar results as in
%\Figs{non-mono-spec}{non-mono-spec_E2} but for the low Mach-number
%AB: comma, and added "now"
\Figs{non-mono-spec}{non-mono-spec_E2}, but now for the low Mach-number
case with vortical forcing (Run~V1). Here there is no indication of
anything more than a single peak.
This peak, which is due to the classical eddy-clustering of Maxey--Riley, 
%is found  around $\St\approx1$, as expected.
%AB: added "to be" and "int"?
is found to be around $\St_{\rm int}\approx1$, as expected.

\begin{figure}\begin{center}
\includegraphics[width=.99\columnwidth]{fig/non-mono-spec_V1.eps}
%AB: Nils, you did this figure. Could you make the St italics and add int?
%AB: The characters should also be a bit smaller 0.65 is what I used.
%AB: My routines are also checked in.
\end{center}\caption[]{
%Run~V1 with more particle sizes.
%AB: inserted "Similar to \Fig{non-mono-spec}, but for"
Similar to \Fig{non-mono-spec}, but for Run~V1 with more particle sizes.
}\label{non-mono-spec_V1}\end{figure}

\begin{figure}\begin{center}
\includegraphics[width=.49\columnwidth]{fig/pcomp_spec_christer_Ma_hel}
\includegraphics[width=.49\columnwidth]{fig/pcomp_spec_christer_Ma_noise}
\end{center}\caption[]{
    Comparison of Run~V3a (solid lines),
    as well as Run~V3s01 (dashed-dotted lines) and Run~V3s02 (dotted lines)
    for (a) $P_n(k)$ and (b) $P_n(k)-P_{\rm noise}$.
}\label{pcomp_spec_christer_Ma_hel}\end{figure}

\subsubsection{Shock capturing viscosity}

Finally, let us discuss how well the clustering results can be modelled
using lower numerical resolution together with a shock viscosity to
stabilise the code.
To be able to do this, it is interesting to have an even larger value of $\Ma$.
%Therefore, we consider a simulation with $\Ma=1.14$ and compare with
%AB: simulation -> DNS, and added "(Run~V3a)"
Therefore, we consider a DNS with $\Ma=1.14$ (Run~V3a) and compare with
two LES with different values of $C_{\rm shock}$; see \Eq{zeta_shock}.
The result for the density spectra is shown in \Fig{pcomp_spec_christer_Ma_hel}.
%We see that for all values of $\St$, except for $\St=0.31$, 
%AB: add "int", changed 0.31 -> 0.53 (the table value was correct and legend wrong, but now corrected)
We see that for all values of $\St_{\rm int}$, except for $\St_{\rm int}=0.53$, 
the LES spectra (dashed and dashed-dotted lines) are close to the DNS
for $k/k_1<20$.
%For $\St=0.31$ (red lines), however, the agreement exists only up to
%AB: add "int"?
For $\St_{\rm int}=0.53$ (red lines), however, the agreement exists only up to
$k/k_1\approx6$.
Surprisingly, a similar departure is not seen in the kinetic energy and
fluid density spectra shown in \Fig{pspec_kine_comp}.
A major difference is, of course, that the LES do not resolve the
small length scales at all, which is also why their spectra are shorter.
The discrepancy in the particle density spectra between LES and DNS
therefore suggests that the clustering, which occurs mostly at those
%intermediate values of $\St$, depends on physical effects at
%AB: add "int"?
intermediate values of $\St_{\rm int}$, depends on physical effects at
the scale of the shocks, corresponding to high wavenumbers.
If this is indeed the case, this departure between DNS and LES may
become worse at larger values of $\Ma$.
%This is something for the future to verify.
%AB: modified
It will be interesting to revisit this question in future simulations
at higher resolution.

\begin{figure}\begin{center}
\includegraphics[width=\columnwidth]{fig/pspec_kine_comp}
\end{center}\caption[]{
$E_{\rm K}(k)$ (left) and $P_\rho(k)$ (right) for Runs~V3 and V3a
(DNS with $\nu=0.05$ and $0.02$, respectively) as well as Runs~V3s01,
V3s02, and V3s1 (LES with $C_\nu=0.1$, $0.2$, and, $1$ respectively).
}\label{pspec_kine_comp}\end{figure}

We should point out that the name LES is, in the present context,
somewhat of a misnomer, because here all the eddies are actually
resolved.
This is because the Kolmogorov scale in Run~V3 is about
three times as large as the mesh spacing; see \Tab{tab:eta}.
In the LES, the mesh spacing is four times larger than in the DNS,
which means that we are resolving down to almost the Kolmogorov scale.

%print,eta,dx: 0.0423458     0.0490874  (V1)
%print,eta,dx: 0.0443534     0.0245437 (V2b)
%print,eta,dx: 0.079597353   0.0245437 (V3a)

%V1   isothermal_helical_256_f0.5_dtf1.0_5adjagr_tsc_nu1_noshocknu/"
%V2a  isothermal_helical_256_f4.0_dtf1.0_5adjagr_tsc_nu0.02_noshocknu'  %(t=262)
%V2b  isothermal_helical_512_f4.0_dtf1.0_5adjagr_tsc_nu0.005_noshocknu_512proc %(t=26)
%     isothermal_helical_512_f4.0_dtf1.0_5adjagr_tsc_nu0.005_noshocknu_512proc
%     isothermal_helical_512_f4.0_dtf1.0_5adjagr_tsc_nu0.005_noshocknu
%V2c  isothermal_helical_512_f4.0_dtf1.0_5adjagr_tsc_nu0.003_noshocknu_512proc %(crashed)
%V2_  isothermal_helical_512_f4.0_dtf1.0_5adjagr_tsc_nu0.01_noshocknu_512proc %(????)
%V2   isothermal_helical_256_f4.0_dtf1.0_5adjagr_tsc_nu0.01_noshocknu/"
%V3   isothermal_helical_512_f10_dtf1.0_5adjagr_tsc_nu0.05_noshocknu_512proc %(used to be called V3a)
%V3a  isothermal_helical_512_f10_dtf1.0_5adjagr_tsc_nu0.02_noshocknu_512proc %(used to be called V3)
%C1   isothermal_expwave_256_f0.5_dtf0.5_5adjagr_tsc_nu5e-3_noshocknu/"
%C1a  isothermal_expwave_512_f0.5_dtf0.5_5adjagr_tsc_nu2e-3_noshocknu
%     isothermal_expwave_256_f0.5_dtf0.5_5adjagr_tsc_nu2e-3_noshocknu
%C2   isothermal_expwave_256_f4.0_dtf1.0_5adjagr_tsc_nu0.05_noshocknu/" -> 0.02 at 512
%C2a  isothermal_expwave_512_f4.0_dtf1.0_5adjagr_tsc_nu0.02_noshocknu

\begin{table}%[]
    \centering
    \begin{tabular}{lccccr}
     Run & $\ell_{\rm Kol}$ & $\delta x$ & $\epsK$ & $\Ma$ & $\Reyn$ \\
    \hline            
    V1  & 0.042 & 0.049 & 0.0003& 0.15 &101 \\%isothermal_helical_256_f0.5_dtf1.0_5adjagr_tsc_nu1_noshocknu
    V2b & 0.044 & 0.025 & 0.034 & 0.72 & 95 \\%isothermal_helical_512_f4.0_dtf1.0_5adjagr_tsc_nu0.005_noshocknu_512proc
    V3a & 0.080 & 0.025 & 0.199 & 1.14 & 38 \\%isothermal_helical_512_f10_dtf1.0_5adjagr_tsc_nu0.02_noshocknu_512proc
  V3s01 & 0.074 & 0.098 & 0.267 & 1.15 & 38 \\%isothermal_helical_128_f10_dtf1.0_5adjagr_tsc_nu0.02_nushock01N_128proc
  V3s02 & 0.076 & 0.098 & 0.234 & 1.14 & 38 \\%isothermal_helical_128_f10_dtf1.0_5adjagr_tsc_nu0.02_nushock02N_128proc
  C2    & 0.213 & 0.049 & 0.064 & 0.74 & 15 \\%isothermal_expwave_256_f4.0_dtf1.0_5adjagr_tsc_nu0.05_noshocknu
  C2a   & 0.107 & 0.025 & 0.061 & 0.76 & 38 \\%isothermal_expwave_512_f4.0_dtf1.0_5adjagr_tsc_nu0.02_noshocknu
  C2s05 & 0.106 & 0.098 & 0.065 & 0.75 & 37 \\%isothermal_expwave_512_f4.0_dtf1.0_5adjagr_tsc_nu0.02_nushock05N_128proc
  C2s1  & 0.106 & 0.098 & 0.063 & 0.74 & 37 \\%isothermal_expwave_512_f4.0_dtf1.0_5adjagr_tsc_nu0.02_nushock1N_128proc
    \end{tabular}
    \caption{
Summary of the Kolmogorov scale $\ell_{\rm Kol}=(\nu^3/\epsK)^{1/4}$,
mesh spacing $\delta x$, energy dissipation rate $\epsK$, as well as Mach
and Reynolds numbers.
    \newline
    \label{tab:eta}}
\end{table}

\begin{figure}\begin{center}
\includegraphics[width=.49\columnwidth]{fig/pcomp_spec_christer_Ma_C2}
\includegraphics[width=.49\columnwidth]{fig/pcomp_spec_christer_Ma_C2_noise}
\end{center}\caption[]{
    Comparison of Run~C2 (solid lines),
    as well as Run~C2s05 (dashed-dotted lines) and Run~C2s1 (dotted lines)
    for (a) $P_n(k)$ and (b) $P_n(k)-P_{\rm noise}$.
}\label{pcomp_spec_christer_Ma_hel_C2}\end{figure}

%AB: we forgot to describe our last figure. How about this?
The quality of LES is worse for runs with compressive forcing; see
\Fig{pcomp_spec_christer_Ma_hel_C2}, where we show the results for Run~C2
with $\Ma=0.76$.
We see that for $\St_{\rm int}=3.6$ (green lines), the agreement between
DNS and LES rather poor even for small values of $k/k_1$.
For $\St_{\rm int}=0.36$ (red lines), the agreement is slightly better,
but again only for small values of $k/k_1$.
%AB.

\section{Conclusions}

In this work we have used particle power spectra to investigate
particle clustering for compressible isotropic turbulence. 
We have studied the effect of using either the Epstein drag, which applies
%to $\Kn\gg1$, or the modified Stokesian drag, which applies to $\Kn\ll1$.
%AB: to -> for (in 2 cases)
for $\Kn\gg1$, or the modified Stokesian drag, which applies for $\Kn\ll1$.
%As long as the particle sizes are non-dimensionalised
%AB: sizes -> radii 
As long as the particle radii are non-dimensionalised
with the Stokes number, the power spectra resulting from the two drag
laws turned out to be similar.
This supports the general usefulness of the Stokes number -- even
for compressible flows and very diverse drag laws.

When using the concept of power spectra to analyse particle clustering
of Lagrangian particles, it is important to realise that the number of
particles used will have an effect on the power at large wavenumbers
(small scales).
%This is due to the fact that, if too few particles are used, it is
%impossible to identify clusters at small scales, because there are not
%enough particles to populate such small clusters.
%AB: reorganized. Please check.
This is due to the fact that, if too few particles are used, there are
not enough particles to populate such small clusters and it becomes
impossible to identify clusters at small scales.
%AB.
This effect is clearly seen through the presence of a $k^2$ contribution
%to the power spectra at small scales.
%AB: to -> in
in the power spectra at small scales.
The magnitude of the $k^2$ contribution is proportional to the inverse
of the total number of particles in the simulation.
If the Eulerian approach is used for the particles, it is implicitly
assumed that an infinite number of particles are involved.
%This means that the $k^2$ contribution to the power spectra is avoided.
%AB: means -> implies, avoided -> absent
This implies that the $k^2$ contribution to the power spectra is absent.
However, this only applies to cases with an infinite number of particles.

When using the Eulerian particle approach, multi-valued particle
velocities are not possible.
%This means that the Eulerian approach cannot represent caustics, which
%AB: means that -> is because  
This is because the Eulerian approach cannot represent caustics, which
implies that only very small Stokes numbers would be modelled correctly.
%For larger Stokes numbers, artificial diffusion and viscosity
%AB: start with "Furthermore"
Furthermore, for larger Stokes numbers, artificial diffusion and viscosity
must be used for the particle fluid in order to stabilise the simulations.
%This yields nonphysical results.
%AB: yields -> can yield 
This can yield nonphysical results.
%Hence, the Eulerian particle approach should only be used to simulate
%AB: should -> can
Hence, the Eulerian particle approach can only be used to simulate a
%very large number of particles that are so small that they behave
%AB: so small -> small enough so 
very large number of particles that are small enough so that they behave
almost as tracers.

\begin{figure}
\centering
%\includegraphics[width=.6\textwidth]{fig/clustering_drawing.pdf}
\includegraphics[width=.8\textwidth]{ppt/sketch.png}
\caption{
Sketch summarising the various clustering mechanisms discussed in this paper.
\label{fig:clustering}}
\end{figure}

We have seen significant differences in the clustering depending
on how the turbulence is generated.
For vortical forcing, the clustering peaks at
%a Stokes number of around unity.
%AB: would it be correct to write here "Kolmogorov-based"?
a Kolmogorov-based Stokes number of around unity.
However, when we use compressive forcing, we drive strong flow divergences.
In that case, we find that clustering peaks at
two different Stokes numbers, one somewhat below unity, and the other
at much larger Stokes numbers.
We argue that the first peak is explained by shock clustering,
while the second is the usual Maxey-Riley clustering
%(the centrifugal sling effect), 
%AB: inserted "based on" for better readability
(based on the centrifugal sling effect), 
and its Stokes number is found to be around unity if it is evaluated
based on the vortical part of the velocity field; see
\Fig{fig:clustering} for a summary.

In order to resolve shocks in high Mach number DNS, a very fine mesh
%is required. In many cases, the interesting physics is, however, not
is required.
In many cases, interesting physics is, however, not
%related to the internal structure of the shocks themselves. 
%AB: added "but to the flows outside the shocks"
related to the internal structure of the shocks themselves,
but to the flows outside the shocks.
%AB.
%It is therefore interesting to model the shocks instead
%AB: modified
It is therefore useful to model shocks through a shock viscosity instead
of resolving them, such that a coarser mesh can be used.
%We investigated the effect on the particle clustering of using a 
%AB: of -> by
We investigated the effect on the particle clustering by using a 
shock-capturing viscosity to broaden the shocks.
For the simulations performed here, the shock-broadening
%meant that the mesh could have four times less grid points in each
%AB: could -> was allowed to 
meant that the mesh was allowed to have four times less grid points in
each direction.
We found that the cases with shock-capturing viscosity reproduced the
%results of the fully resolved DNS simulations for the first decade of wavenumbers
%AB: the -> a, skip simulations 
results of a fully resolved DNS for the first decade of wavenumbers
rather accurately.
However, the relatively strong clustering at small scales that was found
%for DNS simulations of particles with Stokes numbers slightly less than
%AB: DNS simulations -> the DNS
for the DNS of particles with Stokes numbers slightly less than
unity were not reproduced.

%AB: Nils, do you know whether the appendix comes before or after the
%AB: acknowledgements? I've now moved the acknowledgements to here.
\section*{Acknowledgements}
This work was supported by the Knut and Alice Wallenberg Foundation
through the grant Dnr.\ KAW 2014.0048 on
``Bottlenecks for particle growth in turbulent aerosols.''
We thank the Swedish Research Council, grant numbers 2015-04505 (ML)
and 2019-04234 (AB), for partial financial support.
The simulations were performed using resources provided by
the Swedish National Infrastructure for Computing (SNIC)
at the Royal Institute of Technology in Stockholm.

\appendix

\section{Front speed versus gas speed}
\label{pposition}

In \Sec{Applicability}, we used counter-propagating acoustic waves to
drive inertial particles into each other.
The speed of these waves is equal to the sound speed when the speed is
small, but can become comparable to the gas speed for large velocities.
This is shown in \Fig{fig:pposition}.
Note that, even at subsonic gas speeds, the wave speed may exceed the
speed of sound.
Interestingly, the front speed is slightly faster than the associated
Doppler speed, $\cs+u_{\max}$, but slower than $(\cs^2+u_{\max}^2)^{1/2}$.

\begin{figure}
\centering
\includegraphics[width=.6\textwidth]{fig/pposition.eps}
\caption{
Front speed versus gas speed from the one-dimensional experiment described
in the text (plus signs and black line), compared with the associated
Doppler speed, $\cs+u_{\max}$ (blue) and $(\cs^2+u_{\max}^2)^{1/2}$ (red).
\label{fig:pposition}}
\end{figure}

\bibliography{Refs.bib}
\end{document}
