A modal analysis of the behavior of inertial particles in turbulence

TL;DR

This paper provides an analytical expression for the Lyapunov exponents predicting inertial particle clustering behavior in turbulence, explaining the non-monotonic trend and the transition from clustering to dispersion.

Contribution

The study introduces an analytical model for Lyapunov exponents that accurately predicts inertial particle clustering and crossing in turbulence, supported by simulations.

Findings

Lyapunov exponents predict clustering at St=O(1)

Trajectory crossing occurs in hyperbolic flows with Q<0

Clustering rate is bounded and highly nonlinear in hyperbolic flows

Abstract

The clustering of small heavy inertial particles subjected to Stokes drag in turbulence is known to be minimal at small and large Stokes number and substantial at $\rm St = \mathcal O(1)$. This non-monotonic trend, which has been shown computationally and experimentally, is yet to be explained analytically. In this study, we obtain an analytical expression for the Lyapunov exponents that quantitatively predicts this trend. The sum of the exponents, which is the normalized rate of change of the signed-volume of a small cloud of particles, is correctly predicted to be negative and positive at small and large Stokes numbers, respectively, asymptoting to $\tau Q$ as $\tau \to 0$ and $\tau^{-1/2} |Q|^{1/4}$ as $\tau \to \infty$, where $\tau$ is the particle relaxation time and $Q(\tau)$ is the difference between the norm of the rotation- and strain-rate tensors computed along the particle trajectory. Additionally, the trajectory crossing is predicted only in hyperbolic flows where $Q<0$ for sufficiently inertial particles with a $\tau$ that scales with $|Q|^{-1/2}$. Following the onset of crossovers, a transition from clustering to dispersion is predicted correctly. We show these behaviors are not unique to three-dimensional isotropic turbulence and can be reproduced closely by a one-dimensional mono-harmonic flow, which appears as a fundamental canonical problem in the study of particle clustering. Analysis of this one-dimensional canonical flow shows that the rate of clustering, quantified as the product of the Lyapunov exponent and particle relaxation time, is bounded by $-1/2$, behaving with extreme nonlinearity in the hyperbolic flows and always remaining positive in the elliptic flows. These findings, which are stemmed from our analysis, are corroborated by the direct numerical simulations.

Tables

Table 1. Table 1: The leading order terms in the asymptotic solution of Eq. ( 18 ). The solution of the first order expansion ( n = 1 𝑛 1 n=1 ) for λ 𝜆 \lambda is already provided in Eq. ( 30 ). For n > 1 𝑛 1 n>1 , ϵ i subscript bold-italic-ϵ 𝑖 \boldsymbol{\epsilon}_{i} and λ 𝜆 \lambda must be calculated iteratively. For arbitrary tensors 𝑨 𝑨 \boldsymbol{A} and 𝑩 𝑩 \boldsymbol{B} , ( 𝑨 ⋅ 𝑩 ) i ​ j subscript ⋅ 𝑨 𝑩 𝑖 𝑗 (\boldsymbol{A}\cdot\boldsymbol{B})_{ij} is defined as 1 2 ​ ( A i ​ k ​ B k ​ j + A j ​ k ​ B k ​ i ) 1 2 subscript 𝐴 𝑖 𝑘 subscript 𝐵 𝑘 𝑗 subscript 𝐴 𝑗 𝑘 subscript 𝐵 𝑘 𝑖 \frac{1}{2}\left(A_{ik}B_{kj}+A_{jk}B_{ki}\right) .

$n$	$i$	${ϵ}_i\left(\omega \right)$	$n\left({\lambda }^2+\lambda \right)$
1	1	$𝑮\left(\omega \right){\left(1+2\lambda +\hat{i}\omega \right)}^{-1}$	$-2{ϵ}_1\left(\omega \right):{ϵ}_1^T\left(-\omega \right)$
2	1	$\left(𝑮\left(\omega \right)-2{ϵ}_2\left(\omega \right)\cdot {ϵ}_1\left(-\omega \right)\right){\left(1+2\lambda +\hat{i}\omega \right)}^{-1}$
2	2	$-{ϵ}_1\left(\omega \right)\cdot {ϵ}_1\left(\omega \right){\left(1+2\lambda +2\hat{i}\omega \right)}^{-1}$	$-2{\sum }_{j=1}^n{ϵ}_j\left(\omega \right):{ϵ}_j^T\left(-\omega \right)$
3	1	$\left(𝑮\left(\omega \right)-2{ϵ}_2\left(\omega \right)\cdot {ϵ}_1\left(-\omega \right)-2{ϵ}_3\left(\omega \right)\cdot {ϵ}_2\left(-\omega \right)\right){\left(1+2\lambda +\hat{i}\omega \right)}^{-1}$
3	2	$\left(-2{ϵ}_3\left(\omega \right)\cdot {ϵ}_1\left(-\omega \right)-{ϵ}_1\left(\omega \right)\cdot {ϵ}_1\left(\omega \right)\right){\left(1+2\lambda +2\hat{i}\omega \right)}^{-1}$
3	3	$-2{ϵ}_2\left(\omega \right)\cdot {ϵ}_1\left(\omega \right){\left(1+2\lambda +3\hat{i}\omega \right)}^{-1}$	$-2{\sum }_{j=1}^n{ϵ}_j\left(\omega \right):{ϵ}_j^T\left(-\omega \right)$

Equations

$$\mathrm{St}:=\frac{\tau }{{\tau }_{\eta }}=\frac{1}{18}\left(\frac{{\rho }_{\mathrm{p}}}{{\rho }_{\mathrm{f}}}\right){\left(\frac{d_{\mathrm{p}}}{\eta }\right)}^2,$$

$$\ddot{\mathbf{x}}=\mathbf{u}(\mathbf{x},t)-\dot{\mathbf{x}},$$

$$\overline{\nabla \cdot \mathbf{u}}=0,$$

$$\mathcal{C}(\mathrm{St}):=\overline{\nabla \cdot \dot{\mathbf{x}}}\approx -\overline{\frac{\partial u_i}{\partial x_j}\frac{\partial u_j}{\partial x_i}}=\overline{\Vert \mathbf{\Omega }{\Vert }^2-\Vert \mathbf{S}{\Vert }^2}=Q(\mathrm{St}),$$

$$\mathcal{C}(\mathrm{St})={\int }_{-\infty }^{\infty }\frac{{\tilde{\rho }}^{\mathrm{Q}}(\omega ;\mathrm{St})}{1+{\omega }^2}\mathrm{d}\omega ,$$

$${\tilde{\rho }}^{\mathrm{Q}}(\omega ;\mathrm{St}):=\frac{1}{2\pi }{\int }_{-\infty }^{\infty }{\rho }^{\mathrm{Q}}(t;\mathrm{St})e^{-\hat{i}\omega t}\mathrm{d}t=\frac{1}{2\pi }{\int }_{-\infty }^{\infty }-\overline{\frac{\partial u_i(t^{\prime })}{\partial x_j}\frac{\partial u_j(t^{\prime }+t)}{\partial x_i}}{e}^{-\hat{i}\omega t}\mathrm{d}t.$$

$$x_i(\mathbf{X}+\delta \mathbf{X},t)=x_i(\mathbf{X},t)+\frac{\partial x_i}{\partial X_j}\delta X_j.$$

$$J_{ij}(t):=\frac{\partial x_i(t)}{\partial X_j},$$

$$V(t):=\det \left(x_i(\mathbf{X}+\delta {\mathbf{X}}^j,t)-x_i(\mathbf{X},t)\right),$$

$$x_i(\mathbf{X}+\delta {\mathbf{X}}^j,t)-x_i(\mathbf{X},t)=J_{ik}\delta X_k^j.$$

$$V(t)=V(0)\det \left(\mathbf{J}(t)\right),$$

$$V(t)=V(0)\exp ({\mathcal{C}}^{t}t),$$

$${\mathcal{C}}^t:=\frac{1}{t}\ln \left[\det \left(\mathbf{J}(t)\right)\right],$$

$$\mathcal{C}=\underset{t\to \infty }{\lim }{\mathcal{C}}^t.$$

$${\ddot{J}}_{ij}+{\dot{J}}_{ij}=\frac{\partial u_i}{\partial x_k}{J}_{kj},$$

$$F_{ij}:={\dot{J}}_{ik}{J}_{kj}^{-1},$$

$${\dot{F}}_{ij}={\ddot{J}}_{ik}{J}_{kj}^{-1}-F_{ik}{F}_{kj}$$

$${\dot{F}}_{ij}+F_{ik}{F}_{kj}+F_{ij}=\frac{\partial u_i}{\partial x_j},$$

$$\mathrm{tr}(\mathbf{F})=\det {\left(\mathbf{J}\right)}^{-1}\frac{\mathrm{d}\left[\det \left(\mathbf{J}\right)\right]}{\mathrm{d}t},$$

$$\mathcal{C}=\overline{\mathrm{tr}(\mathbf{F})},$$

$${\dot{F}}_{ij}+F_{ik}{F}_{kj}+F_{ij}=\sum _{\omega }{G}_{ij}(\omega )e^{\hat{i}\omega t},$$

$$\mathbf{G}(\omega ):=\frac{1}{2\pi }{\int }_{-\infty }^{\infty }\nabla \mathbf{u}(t)e^{-\hat{i}\omega t}\mathrm{d}t.$$

$$\mathbf{F}=\lambda \mathbf{I}+\sum _{\omega }\mathbf{\Psi }(\omega )e^{\hat{i}\omega t},$$

$$\mathcal{C}=n\lambda ,$$

$$(\lambda +{\lambda }^2)I_{ij}+\sum _{\omega }(1+2\lambda +\hat{i}\omega ){\Psi }_{ij}(\omega )e^{\hat{i}\omega t}+\sum _{{\omega }_n}\sum _{{\omega }_m}{\Psi }_{ik}({\omega }_n){\Psi }_{kj}({\omega }_m)e^{\hat{i}({\omega }_n+{\omega }_m)t}=\sum _{\omega }{G}_{ij}(\omega )e^{\hat{i}\omega t},$$

$$(\lambda +{\lambda }^2)I_{ij}+\sum _{\omega }(1+2\lambda +\hat{i}\omega ){\Psi }_{ij}(\omega )e^{\hat{i}\omega t}+\sum _{\omega }{\Psi }_{ik}(\omega ){\Psi }_{kj}(-\omega )=\sum _{\omega }{G}_{ij}(\omega )e^{\hat{i}\omega t}.$$

$${\Psi }_{ij}(\omega )=(1+2\lambda +\hat{i}\omega {)}^{-1}{G}_{ij}(\omega ).$$

$$\lambda +{\lambda }^2+\frac{1}{n}\sum _{\omega }{\Psi }_{ij}(\omega ){\Psi }_{ji}(-\omega )=0,$$

$$\lambda +{\lambda }^2+\frac{1}{n}\sum _{\omega }\frac{G_{ij}(\omega )G_{ji}(-\omega )}{(1+2\lambda {)}^2+{\omega }^2}=0.$$

$$\lambda +{\lambda }^2-\frac{1}{n}{\int }_{-\infty }^{\infty }\frac{{\tilde{\rho }}^{\mathrm{Q}}(\omega )}{(1+2\lambda {)}^2+{\omega }^2}\mathrm{d}\omega =0,$$

Full text

1 Introduction

The characterization, prediction, and design of a wide range of applications rely on understanding the dynamics of inertial particles in turbulent flows in general and particle clustering or dispersion in particular. The formation of planets and planetesimals in our early solar system is hypothesized to be a result of particle clustering [1, 2, 3]. The accurate prediction of weather relies on proper modeling of droplet coalescence in the clouds, a phenomenon with a direct connection to the problem of particle clustering [4, 5, 6]. The design of a particle-based solar receiver intended for the endothermic chemical reaction with a high-temperature requirement can be achieved through particles clustering [7, 8, 9, 10]. Apart from the few preceding examples, numerous areas in science and engineering, ranging from sediment and pollutant transport in environmental flows to pharmaceutical applications, also benefit from the better understanding of particle clustering and dispersion [11, 12].

In a physically realistic scenario, particle clustering is governed by multiple non-dimensional parameters such as the density ratio of particle to carrier flow, the particle volume fraction, the size of particles relative to the characteristic length scale of the flow, the Froude, Stokes, Reynolds, and Knudsen numbers, the effect of the particle shape as well as the anisotropy of the underlying flow. While a wide range of behaviors is observed in different regions of this multi-dimensional parameter space, particle clustering is primarily governed by a small subset of parameters. To make the problem analytically tractable, we only focus on that small subset of parameters and neglect the effect of other parameters that are of secondary importance. Namely, we neglect the effect of the particles on the flow as well as particle-particle interactions, both of which imply low volume and mass fractions. The effect of thermal fluctuations and body forces on the particle motion are neglected considering regimes with low Knudsen and high Froude numbers, respectively. In our derivation, we will be considering statistically isotropic flows. The finite-size effects are neglected by adopting a point-particle approach, which is to consider only particles that are much smaller than the smallest structure of the flow in the absence of particles. More specifically, we focus on regimes with $d_{\rm p}/\eta\ll 1$, where $d_{\rm p}$ is the particle diameter and $\eta$ is a length scale of the flow that characterizes curvature of the velocity field. Additionally, we assume that the particle Reynolds number is sufficiently small and that the drag is linearly proportional to the slip velocity. We also focus on the regimes at which the particles are much denser than the fluid, i.e. $\rho_{\rm p}/\rho_{\rm f}\gg 1$ with $\rho_{\rm p}$ and $\rho_{\rm f}$ denoting the particle density and the fluid density, respectively.

In such a simplified regime, the particle clustering becomes a function of the particle Stokes number primarily. The Stokes number represents the inertia of the particles and is defined as the ratio of the particle relaxation time $\tau:=\rho_{\rm p}d_{\rm p}^{2}/(18\rho_{\rm f}\nu)$, where $\nu$ is the fluid kinematic viscosity, to a flow time scale. For homogeneous isotropic turbulent flows, the flow time scale is often taken to be the Kolmogorov time scale $\tau_{\eta}:=\sqrt{(\nu/\epsilon)}$, in which $\epsilon$ is the mean volumetric dissipation rate, respectively. Since the Kolmogorov length scale $\eta:=(\nu^{3}/\epsilon)^{1/4}$ is also used for the flow length scale, the Stokes number can be defined as

[TABLE]

which depends on the density ratio $(\rho_{\rm p}/\rho_{\rm f})$ and the length scale ratio $(d_{\rm p}/\eta)$. For the regime under consideration, $(\rho_{\rm p}/\rho_{\rm f})$ is large and $(d_{\rm p}/\eta)$ is small, thus St is finite and can vary from 0 to $\infty$, a range that is investigated in this study. In practice, this regime translates to gas flows laden with small dense droplets or solid particles [7, 8].

In addition to the Stokes number, the underlying flow must be fully defined for the characterization of particle clustering. Among the wide range of possibilities, certain canonical flows are studied the most for their physical relevance or their fundamental importance. If we limit this investigation to stationary isotropic turbulence as a particular class of flows, the flow Reynolds number will appear as an additional non-dimensional parameter that affects particle clustering. The Reynolds number determines the range of time scales that a particle encounters along its trajectory.

Particle clustering, as a function of the particle Stokes number and the flow Reynolds number, has been studied in the past extensively. It has been observed experimentally [13, 14, 15, 16, 11, 17], simulated numerically [18, 19, 20, 21, 22, 23], and described analytically [24, 25, 26, 27, 28]. These studies have shown that the degree to which particles clustering mainly depends on the Stokes number. Depending on their Stokes number, particles may homogeneously disperse in space or preferentially concentrate in certain regions of the flow. For ${\rm St}\ll 1$, particles become neutral fluid tracers and experience minimum clustering. There is also minimal clustering at the limit of St $\gg 1$, where particles follow a ballistic trajectory uncorrelated with the underlying flow. The maximum clustering is achieved when ${\rm St}\approx 1$. This non-monotonic variation, which has been observed experimentally and numerically, is yet to be described analytically. In this study, we attempt to describe this trend using analytical tools while employing numerical simulations for verification purposes.

The particle motion subjected to the preceding simplifications is well characterized by the Stokes drag. Denoting the position of a particle by $\boldsymbol{x}(t)$ and the flow velocity at the particle location by $\boldsymbol{u}(\boldsymbol{x},t)$, the dimensionless equation of the motion of a particle subjected to the Stokes drag is

[TABLE]

where $\dot{(\bullet)}:={\rm d}(\bullet)/{\rm d}t$. We employed $\tau$, $L$, and $U=L/\tau$ as the time, length, and velocity scales to non-dimensionalize corresponding parameters in Eq. (2). All the following equations are also non-dimensionalized based on $\tau$ and $L$. $\tau$ is the particle relaxation time, defined above, and $L$ is a characteristic length scale. The following formulations are independent of the choice of $L$, hence its choice is arbitrary. Parameters normalized based on the flow time scale rather than the particle relaxation time are distinguished by the subscript $\eta$. The only exception is $\tau_{\eta}$ which denotes the dimensional Kolmogorov time scale rather than the particle relaxation time normalized by the flow time scale. Although Eq. (2) appears to be a simple ordinary differential equation, it behaves nonlinearly due to the dependence of $\boldsymbol{u}$ on $\boldsymbol{x}$. It is this nonlinear behavior that gives rise to particle clustering, which is a complex and nonlinear phenomenon.

One of the first analytical relations for quantifying particle clustering was derived in [25] and [24]. This relation, which hereafter is referred to as RM, is obtained by approximating the acceleration of particles with that of the fluid. One of the underlying assumptions of this analysis is

[TABLE]

where $\overline{(\bullet)}$ denotes Lagrangian time-averaging along the trajectory of a particle. This condition, which is less restrictive than $\nabla\cdot\boldsymbol{u}=0$, only requires the flow to be incompressible along particle trajectory over a long period. Thus, this condition can also be satisfied for compressible flows where the total mass of a given control volume remains the same over a long time interval. The central assumption of RM, however, is approximating the acceleration of particles with that of the flow via $\ddot{\boldsymbol{x}}\approx{\rm D}\boldsymbol{u}/{\rm D}t$. With these two assumptions, one can take the divergence of Eq. (2) to obtain

[TABLE]

where $\boldsymbol{\Omega}$ and $\boldsymbol{S}$ are the dimensionless rotation-rate and strain-rate and $Q$ is half of the Q-criterion, which is defined in the literature for identification of the vortical regions [29, 30]. The sign of $Q$ indicates the dominance of flow rotation-rate ($Q>0$) or strain-rate ($Q<0$), which are associated with elliptic and hyperbolic regions of the flow, respectively. Apart from the fact that $Q$ is normalized by $\tau$, its computation along the particle trajectory implies its dependence on St, i.e. $Q_{\eta}=Q_{\eta}({\rm St})$. The RM expression relates the degree of preferential concentration, here defined as $\mathcal{C}$, to the divergence of the particle velocity field. The sign of $\mathcal{C}$ determines the regime of particle clustering or dispersion. Dispersion occurs when $\mathcal{C}>0$, i.e. when particles move away from each other over time. Clustering occurs when $\mathcal{C}<0$, i.e. when particles get closer to each other over time. The Eulerian definition of $\mathcal{C}$ through divergence operation in Eq. (4) is harder to interpret when particle trajectories cross in which $\dot{\boldsymbol{x}}$ is not a well-defined function of the spatial co-ordinates [27]. Thus, in Section 2, we adopt a Lagrangian definition of $\mathcal{C}$, which mathematically corresponds to the above Eulerian definition, for quantification of the particle clustering.

Approximating the acceleration of particles with that of flow, i.e. $\ddot{\boldsymbol{x}}\approx{\rm D}\boldsymbol{u}/{\rm D}t$, limits the validity of RM to ${\rm St}\ll 1$. The accuracy of RM in predicting the first- and second-order statistics of $\mathcal{C}$ has been shown for ${\rm St}\ll 1$ in homogeneous turbulence [31, 32] as well as synthetic flows [33]. For St $\geq 1$, RM predicts an unbounded $\mathcal{C}$ proportional to St and fails to capture the non-monotonic behavior of $\mathcal{C}({\rm St})$. Additionally, Eq. (4) suggests that particles are repelled from the rotation-dominated regions, where $\|\boldsymbol{\Omega}\|>\|\boldsymbol{S}\|$, and preferentially concentrate in regions with higher strain-rate. These generally accepted qualitative assessments will be examined thoroughly in this study.

In an earlier attempt, we derived an alternative relationship for $\mathcal{C}$ via a first-order correction to RM [32]. In that study, we linearized Eq. (2) by expanding $\boldsymbol{u}(\boldsymbol{x},t)$ using the Taylor series in spatial co-ordinates and expressing its temporal variation in the Fourier space to obtain an eigenvalue problem for the Lyapunov exponents of pairs of inertial particles. We showed that the sum of these exponents (in three dimensions $\lambda_{1}+\lambda_{2}+\lambda_{3}$) is equal to the divergence of the particle velocity field and can be expressed as

[TABLE]

where

[TABLE]

In the last integral $t^{\prime}$ varies along the particle trajectory, thus $\overline{(\bullet)}$ operates on $t^{\prime}$. This definition implies $\rho^{\rm Q}(t)=\rho^{\rm\Omega}(t)-\rho^{\rm S}(t)$, where $\rho^{\rm\Omega}$ and $\rho^{\rm S}$ are the norm of the autocorrelation function of the rotation-rate and strain-rate tensors, respectively, and $\tilde{\rho}^{\rm Q}(\omega)$ is the Fourier transformation of $\rho^{\rm Q}(t)$.

The primary assumption associated with Eq. (5), which we refer to as SL (Small Lyapunov) hereafter, is that the Lyapunov exponents $\lambda_{i}$ are much smaller than 1. This assumption is valid at relatively small St in which $\lambda_{i}$ are also small. Comparing Eqs. (4) and (5), SL can be considered as the filtered RM which accounts for unresponsiveness of particles to high-frequency fluctuations. As St $\to 0$, SL exactly reproduces RM. At this limit, $\omega^{2}$ in the denominator of Eq. (5) can be neglected, thus reproducing Eq. (4) exactly. Both SL and RM are a linear function of the rotation-rate and strain-rate tensors, predicting a $\mathcal{C}$ that changes linearly with $Q$. Due to this linearity, both formulations predict dispersion ($\mathcal{C}>0$) only for rotating flows ($Q>0$) and clustering ($\mathcal{C}<0$) only for straining flows ($Q<0$). Neither of them predicts clustering for a rotating flow or dispersion for a straining flow. While these predictions agree with the numerical simulations at small St, they break down at higher Stokes numbers. Although SL remains bounded and provides a better prediction than RM, it still deviates from the reference for ${\rm St}>\mathcal{O}(1)$ due to the underlying assumption of $|\lambda_{i}|\ll 1$.

The objective of this study is to derive an analytical relationship for $\mathcal{C}({\rm St};\tilde{\rho}^{\rm Q})$ which is valid at both small and large St. In what follows, we present a step-by-step derivation of a solution that is extracted from Eq. (2) for particles experiencing statistically isotropic flow. Then, we compare the prediction of RM, SL, and our solution to the reference numerical results. For this purpose, we first consider a one-dimensional unimodal oscillatory flow. Through this case, we test the accuracy of the proposed solution and discuss its implications. Then, we extend our analysis to more complex multi-dimensional flows, oscillating at a continuous range of frequencies. We consider three-dimensional isotropic turbulence for this purpose, where we compare all the analyses at a wide range of Stokes numbers.

2 Analytical derivation

Consider a collection of four particles that are located at $\boldsymbol{X}$ and $\boldsymbol{X}+\delta\boldsymbol{X}^{i}$, $i=1,2,3$, at $t=0$ (Figure 1). We refer to this collection of inertial particles, which are within an infinitesimal distance from each other, a cloud. At time $t$, the particles in the cloud will move to new positions $\boldsymbol{x}(\boldsymbol{X},t)$ and $\boldsymbol{x}(\boldsymbol{X}+\delta\boldsymbol{X}^{i},t)$, respectively. Without loss of generality, we choose $\delta$ such that the particles are always at relatively small distances from each other. More specifically, $\|\boldsymbol{x}(\boldsymbol{X}+\delta\boldsymbol{X}^{i},t)-\boldsymbol{x}(\boldsymbol{X},t)\|\ll\eta$ is a sufficient condition to ensure particles in the cloud experience a linearly-varying velocity field. Given this condition, the trajectory of an arbitrary particle within the cloud that is initially located at $\boldsymbol{X}+\delta\boldsymbol{X}$ can be described as

[TABLE]

One may choose to consider a cloud as a collection of a larger number of particles than the one that is shown in Figure 1. Nevertheless, based on Eq. (7), the motion of all those particles can be represented as the linear combination of the motion of $n$ particles in an $n$-dimensional space. In other words, the motion of a cloud is fully characterized via Eq. (7) once $\boldsymbol{x}(\boldsymbol{X},t)$ and the linear deformation tensor

[TABLE]

are defined.

Considering the configuration of the cloud that is fully characterized based on $x_{i}(\boldsymbol{X}+\delta\boldsymbol{X}^{j},t)-x_{i}(\boldsymbol{X},t)$ (Figure 1), we define

[TABLE]

as the signed-volume of the cloud at any given time $t$, where ${\rm det}\left({\bullet}\right)$ is the determinant operator. Since $\boldsymbol{x}(\boldsymbol{X},0)=\boldsymbol{X}$, it directly follows from Eq. (9) that $V(0)={\rm det}\left({\delta X^{j}_{i}}\right)$. It is evident from the geometrical relationship between the volume of a tetrahedral and the determinant of its corresponding matrix, the actual volume of the cloud is $|V(t)|/6$, where factor of $1/6$ is specific to a tetrahedral cloud and eliminated from our definition for the sake of notation brevity. The sign of $V(t)$, on the other hand, denotes whether particles are oriented relative to each other following a right-hand rule. We choose the initial orientation of particles such that always $V(0)>0$. Thus, $V(t)<0$ only occurs when one of the particles passes through the plane constructed by the other three (or the line constructed by the other two in two dimensions or crossing over the other particle in one dimension), turning the cloud inside-out. Hereafter, we refer to such incidents as particle trajectory crossing.

To express $V(t)$ in terms of $\boldsymbol{J}$, we rearrange Eq. (7) for a collection of particles and use Eq. (8) to obtain

[TABLE]

Taking the determinant of Eq. (10) and substituting for $V$ based on Eq. (9) yields

[TABLE]

which shows that the determinant of $\boldsymbol{J}(t)$ is the signed-volume of the cloud, which has undergone a linear deformation characterized by $\boldsymbol{J}$, normalized by its initial volume. Although we derived Eq. (11) for a tetrahedral cloud composed of four particles (Figure 1), the result is generalizable to a cloud with an arbitrary number of particles by considering the entire cloud as a collection of tetrahedron elements and repeating the above procedure for each element.

Based one the definition of $V(t)$, we define the finite-time exponential rate of change of the signed-volume, $\mathcal{C}^{t}$, such that

[TABLE]

holds. Combining Eqs. (11) and (12) yields

[TABLE]

which is an alternative definition of $\mathcal{C}^{t}$. As shown in A, $\mathcal{C}^{t}$ is related to $\mathcal{C}$, which was defined earlier as the divergence of particle velocity field, by

[TABLE]

Before proceeding further, it is necessary to build a physical intuition around $\mathcal{C}$. From a Lagrangian perspective, $\mathcal{C}$ is the exponential rate at which $V(t)$ changes over a long period. Its real component, which only depends on $|V(t)|$, has an intuitive meaning. It is the exponential rate of change of the volume of the region in the flow that is contaminated by particles. Thus, ${\rm real}\left({\mathcal{C}}\right)>0$ signifies the regimes of particle dispersion, which is analogous to diffusion induced mixing in continuum flows. ${\rm real}\left({\mathcal{C}}\right)<0$, on the other hand, signifies the regimes of particle clustering in which the collection of particles come together, causing the cloud of particles to form fractals with a dimension smaller than the number of spatial dimensions. From Eq. (12), having a $\mathcal{C}$ with a nonzero imaginary component indicates the oscillatory behavior of $V(t)$. It is the particle trajectory crossing that changes the sign of $V(t)$, leading to ${\rm imag}\left({\mathcal{C}}\right)\neq 0$. According to Eq. (14), $\mathcal{C}^{t}$ is closely related to $\mathcal{C}$. $\mathcal{C}^{t}$ is to $\mathcal{C}$ what the finite-time Lyapunov is to the Lyapunov exponent. By tracking $\mathcal{C}^{t}$ in time, one may identify incidents of trajectory crossing at which $V(t)\to 0$ and ${\rm real}\left({\mathcal{C}^{t}}\right)\to-\infty$. We will use these concepts as we develop an expression for $\mathcal{C}$ for an arbitrary flow and $\rm St$.

Based on Eqs. (13) and (14), $\mathcal{C}$ can be computed once $\boldsymbol{J}$ is determined. To obtain a relationship for $\boldsymbol{J}$, we need the governing equation of motion of particles that is taken to be Eq. (2), when the assumptions enumerated in Section 1 are considered. Taking the derivative of Eq. (2) with respect to $\boldsymbol{X}$ and employing the chain rule, yields

[TABLE]

where $\dot{J}_{ij}=\partial\dot{x}_{i}/\partial X_{j}$ and $\ddot{J}_{ij}=\partial\ddot{x}_{i}/\partial X_{j}$, which hold true since $\boldsymbol{X}\neq\boldsymbol{X}(t)$. Also, $\nabla\boldsymbol{u}=\partial\boldsymbol{u}/\partial\boldsymbol{x}$ is expressed in terms of $\boldsymbol{x}$ and thus is tractable by computing the fluid velocity gradient along the trajectory of the cloud. The tensor $\nabla\boldsymbol{u}$ is a general function of time and as a result, Eq. (15) is not a constant coefficient ordinary differential equation (ODE) to be integrated directly. Nevertheless, analogous to an ODE with constant coefficients, its solution ${\rm det}\left({\boldsymbol{J}}\right)$ exponentially grows or decays indefinitely in time for $\mathcal{C}\neq 0$. To proceed further, we make use of

[TABLE]

a transformation that produces a more tractable constant coefficient ODE. Note $J_{kj}^{-1}=(\boldsymbol{J}^{-1})_{kj}$ is implied in Eq. (16) and what follows. Based on this transformation,

[TABLE]

and hence from Eq. (15)

[TABLE]

which has been derived by others as well [6, 34] and classifies as the Riccati equation in a tensorial form.

Equation (18), in contrast to Eq. (15), is nonlinear but has constant coefficients. It is expressed in terms of $\boldsymbol{F}$, which is an instantaneous rate of deformation of the cloud normalized by its size and is independent of the arbitrary choice of $\boldsymbol{X}$. As a result, its determinant ${\rm det}\left({\boldsymbol{F}}\right)$ is a statistically stationary variable for sufficiently long $t$ for flows reaching an equilibrium. To show that and establish a relationship between $\mathcal{C}$ and $\boldsymbol{F}$, we employ Jacobi’s formula,

[TABLE]

along with Eqs. (13) and (14) as outlined in B to obtain

[TABLE]

in which $\rm\operatorname{tr}(\bullet)$ is the trace operator. This simple relationship indicates that $\mathcal{C}$, which is the sum of the Lyapunov exponents, is the time average of the sum of the eigenvalues of $\boldsymbol{F}$. In other words, the eigenvalues of $\overline{\boldsymbol{F}}$ are the Lyapunov exponents associated with the pairs of inertial particles.

Next, we solve Eq. (18) for $\boldsymbol{F}$ to find an analytical estimate for $\mathcal{C}$ via Eq. (20). Given that the tensor $\nabla\boldsymbol{u}$ is a general function of time, we find a solution for Eq. (18) by expressing $\nabla\boldsymbol{u}$ as a set of harmonic functions using the Fourier transformation. Thus, Eq. (18) can be written as

[TABLE]

where

[TABLE]

In this equation and all that follows, $\omega$, which runs form $-\infty$ to $\infty$, can be considered as $2\pi k/T$ with $k\in\mathbb{Z}$ and $T\to\infty$ being the sampling period. We also express $\boldsymbol{F}$ by taking its Fourier transformation as

[TABLE]

where $\boldsymbol{I}$ is the identity tensor, $\lambda$ is a representative Lyapunov exponent, and tensor $\boldsymbol{\Psi}(\omega)$ is the oscillatory response of $\boldsymbol{F}$ to $\boldsymbol{G}$. In Eq. (23), the steady response, which is expressed as $\lambda\boldsymbol{I}$, is separated from $\boldsymbol{\Psi}$ to simplify the form of following expressions. The form of the steady term implies that the Lyapunov exponents are equal, which need not be true in general. The problems considered in this study, where the flow around the particles is assumed to be statistically isotropic, have no preferred directions. For such problems, all exponents are equal111Distinct $\lambda_{i}$ are reported in the literature even for isotropic flows, which are obtained by sorting each ensemble before averaging. If $\lambda_{i}$ correspondence with $i^{\rm th}$ co-ordinate is preserved – statistical uncertainty associated with the finite sampling aside – then the ensemble averages in all directions will be equal owing to the isotropy condition., and the form of the steady term in Eq. (23) is valid. Also, note that from Eqs. (20) and (23) $\mathcal{C}$ can be computed as

[TABLE]

where $n$ is the number of spatial dimensions.

Our goal is to determine $\lambda$ and $\boldsymbol{\Psi}$ using Eq. (21). Substituting for $\boldsymbol{F}$ in Eq. (21) yields

[TABLE]

To proceed further, we neglect the higher-order terms in Eq. (25) by assuming $\|\boldsymbol{\Psi}\|\ll 1$. We revisit this assumption in C, where the effect of higher order terms is analyzed for a unimodal excitation that is a model problem introduced in Section 3. Based on this assumption, all the time-dependent terms in the second summation can be neglected compared to those in the first summation. The remaining time-independent terms are retained as they might be comparable to the first term. Thus, keeping only the terms with $\omega_{n}=-\omega_{m}$ yields

[TABLE]

For Eq. (26) to hold,

[TABLE]

Additionally, $\sum_{\omega}\Psi_{ik}(\omega)\Psi_{kj}(-\omega)$ must be a diagonal matrix, which occurs in isotropic flows when the non-identical entries of $\boldsymbol{\Psi}$ are uncorrelated222While $\Psi_{ik}(\omega)\Psi_{kj}(-\omega)$ for $i\neq j$ is not necessarily smaller than that of $i=j$, $\sum_{\omega}\Psi_{ik}(\omega)\Psi_{kj}(-\omega)$ for $i\neq j$ will be negligible compared to that of $i=j$ when the sampling period $T$ is sufficiently long as is the case here.. Taking the trace of Eq. (26) yields

[TABLE]

where $n$, as defined above, is the number of spatial dimensions. Due to the long-term incompressibility condition in Eq. (3), $G_{ii}(0)$ is zero and does not appear in Eq. (28). From Eqs. (27) and (28)

[TABLE]

Using the convolution theorem, Eq. (29) is expressed in terms of a continuous integral as

[TABLE]

where $\tilde{\rho}^{\rm Q}$ is defined in Eq. (6).

Remarks on Eq. (30):

This equation, which will be further developed in Sections 3 and 4, is an integral equation. Evaluating the integral that appears in this equation requires a knowledge of $\lambda$, which itself is the solution. Hence, obtaining a closed-form explicit expression for $\lambda$ relies on further simplification and requires additional assumptions. 2. 2.

The expressions given by RM and SL are special forms of Eq. (30), which can be reproduced exactly by adopting further assumptions. Specifically, linearizing Eq. (30) at the limit of $|\lambda|\ll 1$ and using Eq. (24) to express it in terms of $\mathcal{C}$ exactly reproduces SL (Eq. (5)). Furthermore, neglecting $\omega^{2}$ in the denominator reduces the integral to $\rho^{\rm Q}(0)=\overline{\|\boldsymbol{\Omega}\|^{2}-\|\boldsymbol{S}\|^{2}}$, exactly reproducing RM (Eq. (4)). 3. 3.

No length scale appears in this equation, justifying our earlier arbitrary choice of $L$ in Section 1. The choice of time scale, on the other hand, affects terms with $\omega$, $\lambda$, and $\tilde{\rho}^{\rm Q}$. These parameters are all normalized by the particle relaxation time $\tau$ in Eq. (30). In most physical scenarios, $\tilde{\rho}^{\rm Q}$ is governed primarily by the underlying flow rather than $\tau$, and a natural choice for normalization of Eq. (30) is $\tau_{\eta}$. If normalized based on $\tau_{\eta}$, St will appear in the denominator of the integral, which is compatible with the notion that the large St particles barely respond to the flow fluctuations. From that re-normalization, one can show $\lambda_{\eta}\propto{\rm St}Q_{\eta}$ for $\rm St\ll 1$ and $\lambda_{\eta}\propto{\rm St}^{-1/2}|Q_{\eta}|^{1/4}$ for $\rm St\gg 1$, where subscript $\eta$ denotes normalization based on $\tau_{\eta}$ (refer to Section 3.1 for more detail). 4. 4.

Contrary to RM and SL, $\lambda$ or $\mathcal{C}$ are nonlinear functions of $\tilde{\rho}^{\rm Q}$. Scaling $\tilde{\rho}^{\rm Q}$ by a factor of $k$, clouds contract or expand $k$ times faster only if $k\to 0$. In the limit $k\to\infty$ and if $\lambda\neq-1/2$, $\lambda$ will scale as $k^{1/4}$. The numerical results of Section 3.2 confirm this asymptotic prediction. 5. 5.

Similar to RM and SL, the only term that appears in this equation which relates $\lambda$ to the underlying flow is $\tilde{\rho}^{\rm Q}$, which is closely related to the second invariant of the velocity gradient tensor. As expected, this equation is invariant under Galilean transformation (translation and rotation). The absence of other flow-related parameters suggests that particle clustering, although complex, solely depends on the difference between the norms of the autocorrelation function of the rotation- and strain-rate tensors.

To better understand the implication of Eq. (30), it is simplified in the next section to find an explicit expression for $\lambda$ for a case in which the underlying flow oscillates at a single frequency.

3 One-dimensional unimodal excitation

Earlier, we argued that the clustering phenomenon is complex. It exhibits a multiscale behavior (an aspect which is not discussed here) and follows a non-monotonic trend versus Stokes number. One may associate this complexity in its entirety to the complexity of the background flow. Such an assessment might stem from the fact that turbulent flows, which are often employed for studying clustering, are also very complex. On the other hand, one may hypothesize that this complexity is partially a manifestation of the fundamental response of Eq. (2) to what could be considered a far simpler flow than turbulence. To distinguish the role of the two and discern the mechanism by which some complexities in the clustering phenomenon arises, we will consider a simple model problem in this section. Later on, we will show that this simple model problem is fundamentally important in the study of particle clustering as it explains much of our observation when we consider a three-dimensional turbulent flow.

Consider a two-dimensional pure straining flow that oscillates in time with a given frequency and amplitude. That is

[TABLE]

where $\Phi$ and $\omega$ are the root square of amplitude and frequency of the oscillations, respectively, and are independent parameters that are defined based on the underlying flow. Suppose $\Phi<0$ for now, as $\Phi>0$ is a more complex case to be considered later. The factor of $-2$ and square are incorporated into this definition to simplify the notation of the following expressions and allow for a one-to-one correspondence of $\Phi$ and $\tilde{\rho}^{\rm Q}$. For this flow, the velocity gradient tensor is diagonal with $|\partial u_{1}/\partial x_{1}|=|\partial u_{2}/\partial x_{2}|$. The particle motion in the $x$-direction is independent of $u_{2}$ as is the particle motion in the $y$-direction is independent of $u_{1}$. Owing to the diagonal structure of $\nabla\boldsymbol{u}$ and the equivalence of the rate of strain in $x$- and $y$-directions, one can analyze the motion of a particle in one direction and apply the results to both directions. Thus, instead of analyzing two identical problems, we analyze a one-dimensional flow with

[TABLE]

which, in other words, corresponds to the flow around particles located on the $x$-axis in the two-dimensional setting described above.

Both straining and rotating flow can be represented in this one-dimensional setting using $\Phi<0$ and $\Phi>0$, respectively. Interpreting Eq. (32) is simple in the straining regime, since the underlying $u(x,t)$ is real and varies linearly along the line with all fluid particles moving away or toward a single point. Its interpretation in the rotation regime, on the other hand, requires imagining the line to have an oscillatory rotational motion with all fluid particles oscillating along concentric arches. From Eq. (18), one can prove that the rotating regime described in Eq. (32) with $\Phi>0$ produces a Lyapunov exponent that is identical to that of a forced vortex defined by $\boldsymbol{u}=[x_{2},-x_{1}]\sqrt{2\Phi}\cos(\omega t)$ in two dimensions. For this two-dimensional rotating flow, the eigenvalues of the velocity gradient tensor are imaginary and correspond to the single eigenvalue of the one-dimensional model problem. In contrast to the one-dimensional case, however, $\nabla\boldsymbol{u}$ is real in the two-dimensional case, allowing one to carry out the computations in the real plane.

Pertaining to what was discussed in Section 2, the long-term rate at which the distance between two particles changes in a cloud is exponential. To show this in a one-dimensional setting, consider two particles separated by $\delta x$. The difference between the fluid velocity at the location of these particles is $(\partial u/\partial x)\delta x$. Since $\partial u/\partial x$ does not scale with $\delta x$, the rate at which particle move relative to each other scales with their distance, leading to an exponential change of $\delta x$. The exponential rate at which these particles converge or diverge to each other as $t\to\infty$ is by definition the Lyapunov exponent $\lambda$. Depending on the value of $\lambda$, four scenarios may occur. As shown in Figure 2, particles move relative to each other and may diverge (${\rm real}\left({\lambda}\right)>0$) or converge (${\rm real}\left({\lambda}\right)<0$), representing regimes at which particles contaminate a larger or smaller space over time, respectively. The convergence or divergence of particles may occur while their trajectories cross (${\rm imag}\left({\lambda}\right)\neq 0$) or do not cross (${\rm imag}\left({\lambda}\right)=0$). Similar scenarios can be hypothesized for a rotating flow. Some of these scenarios, however, may never occur in reality (e.g., ${\rm real}\left({\lambda}\right)<0$ for heavy particles in a rotating flow). Our goal in this section is to find an analytical relationship for $\lambda$ that firstly delineates between the above scenarios and secondly, provides a quantitative estimate of $\lambda$ when Eq. (32) holds.

The one-dimensional flow described above is characterized using Eq. (6) and Eq. (32) as

[TABLE]

where $\delta$ is the Dirac delta function. The sign of $\tilde{\rho}^{\rm Q}$ reaffirms our earlier distinction between the rotating ($\Phi>0$) and straining ($\Phi<0$) regimes. Substituting $\tilde{\rho}^{\rm Q}$ into Eq. (30) and taking $n=1$ simplifies it to

[TABLE]

In total, there are two non-dimensional parameters that appear in Eq. (34), which are $\Phi$ and $\omega$. In Section 1, we employed particle relaxation time $\tau$ to normalize all parameters. Since $\Phi$ and $\omega$ are normalized by $\tau^{2}$ and $\tau$, respectively, the effect of Stokes number, viz. $\tau_{\eta}/\tau$, is embedded in both parameters. Equation (33) and consequently Eq. (34) also represent multi-dimensional isotropic flows as long as $\tilde{\rho}^{\rm Q}$ contains only a single frequency.

Next, we find a relationship for $\lambda=\lambda(\Phi,\omega)$ using Eq. (34). According to Eq. (34), $\lambda$ is one of the roots of a fourth-order polynomial. Closer examination shows that this polynomial has two real roots for $\Phi>-\omega^{2}/4$ or $\Phi>0$, no real roots for $\Phi<-(\omega^{2}+1)^{2}/16$ or $-(\omega^{2}+1)^{2}/16<\Phi<-\omega^{2}/4$, and four real roots otherwise (Figure 3). We show below that out of four roots of this polynomial, only one is physically relevant.

The roots of Eq. (34) can be analytically computed by converting it to a depressed quadratic form and taking

[TABLE]

With this change of variable, Eq. (34) becomes

[TABLE]

Therefore,

[TABLE]

No contraction or expansion is expected in the absence of flow. Thus, the condition $\lambda=0$, which corresponds to $\gamma=1$, must be satisfied when $\Phi=0$. Therefore, from the two possible solutions in Eq. (37), only the root corresponding to the plus sign is admissible. Computing $\lambda$ from Eqs. (35) and (37) and imposing the condition $\lambda(\Phi=0,\omega)=0$ for one more time gives

[TABLE]

which is valid for $\omega$ and $\Phi\in\mathbb{R}$. Based on this equation $\lambda$ can be complex for a specific combination of $\Phi$ and $\omega$, pertaining to the occurrence of particle crossovers. The magnitude of ${\rm imag}\left({\lambda}\right)$ is crossover frequency and proportional to the number of times particles exchange side in a unit time. The real and imaginary part of $\lambda(\Phi,\omega)$ are plotted in Figure 4.

Remarks on Eq. (38):

${\rm real}\left({\lambda}\right)$ represents the rate of expansion (${\rm real}\left({\lambda}\right)>0$) or contraction (${\rm real}\left({\lambda}\right)<0$). Since the second term in Eq. (38) is always positive, the strongest possible rate of contraction is $-1/2$. For $\omega<1$, the minimum value of ${\rm real}\left({\lambda}\right)$ occurs at the discriminant curve $\Phi=-(\omega^{2}+1)^{2}/16$, at which ${\rm real}\left({\lambda}\right)$ is not differentiable. For $\omega>1$, the minimum occurs in a region enclosed between $-(\omega^{2}+1)^{2}/16<\Phi<-\omega^{2}/4$, where ${\rm real}\left({\lambda(\Phi,\omega)}\right)$ is constant and equal to the global minimum $-1/2$ (Figure 4). 2. 2.

No expansion or contraction is predicted for two cases. The first is the trivial case in which $\Phi=0$. The second case occurs at $\Phi=-(\omega^{4}+10\omega^{2}+9)/16$. In the latter case only ${\rm real}\left({\lambda}\right)=0$ and ${\rm imag}\left({\lambda}\right)\neq 0$, indicating pure oscillation of the volume of the clouds between positive and negative values (Figure 2). 3. 3.

In the rotation dominated regimes where the flow is elliptical and $\Phi>0$, the clouds of particles always experience pure expansion with ${\rm real}\left({\lambda}\right)>0$. No crossovers occurs in this regime (Figure 4). 4. 4.

Contrary to the contraction-rate, the expansion-rate is unbounded. In general ${\rm real}\left({\lambda}\right)$ is proportional to $\Phi$ for $|\Phi|\ll 1$ and to $|\Phi|^{1/4}$ for $|\Phi|\gg 1$. 5. 5.

As discussed above, ${\rm imag}\left({\lambda}\right)\neq 0$ is the byproduct of ${\rm det}\left({\boldsymbol{J}}\right)<0$, occurring in strong-straining flows as the cloud fully collapses and turns inside-out. In the one-dimensional setting, ${\rm imag}\left({\lambda}\right)\neq 0$ occurs when two particles exchange sides on the line. 6. 6.

Particle crossover occurs when $\Phi<\max\left[-(\omega^{2}+1)^{2}/16,\min\left(-\omega^{2}/4,-1/4\right)\right]$ (blue region in Figure 3). The maximum of ${\rm imag}\left({\lambda}\right)$ occurs at $\Phi=-(\omega^{2}+1)^{2}/16$ for $\omega>1$, at which ${\rm imag}\left({\lambda}\right)$ is not differentiable. 7. 7.

Out of four scenarios in Figure 2, only three are predicted to occur in the one-dimensional straining flow and one in the one-dimensional rotating flow. In a straining regime $\{{\rm real}\left({\lambda}\right)>0,\;{\rm imag}\left({\lambda}\right)=0\}$ never occurs while in a rotating regime only $\{{\rm real}\left({\lambda}\right)>0,\;{\rm imag}\left({\lambda}\right)=0\}$ occurs. In other words, particles have to cross each other to diverge in a straining flow, whereas they may or may not cross when they converge. Thus, in a straining flow, particle dispersion only occurs if particle trajectories cross, whereas they can cluster regardless of the occurrence of crossovers. In a rotating flow, they always diverge without crossing over each other. 8. 8.

For a straining flow oscillating at sufficiently small $\omega$, the onset of crossover occurs at $\Phi=-1/16$ whereas the onset of expansion occurs at $\Phi=-9/16$. Therefore, particles that disperse in a straining flow are at least several times more inertial than those which do not cross each other.

Equation (38) is the exact solution of Eq. (30) when the flow is one-dimensional and excited at a single frequency. Thus, the underlying assumption of Eq. (38) is that of Eq. (30). Namely, the high-order oscillatory terms, which originally appeared in Eq. (25), are neglected by assuming $\|\boldsymbol{\Psi}\|\ll 1$. To evaluate the significance of this assumption, we compare the prediction of Eq. (38) to the numerical results in Section 3.2. Before that, however, we discuss the interpretation of these findings when they are expressed as a function of the Stokes number.

3.1 St-dependent re-normalization of Eqs. (30) and (38)

All the derived equations so far have been normalized based on the particle relaxation time $\tau$. The choice of $\tau$ for the time scale led to compact equations in which only an oscillation amplitude and frequency appear. This normalization translates to experiments in which a single class of particles is reused in a variety of flows. Often in practice, however, we encounter multiple classes of particles in one particular flow. To predict the trends observed in the latter case, we normalize the previous results by a flow-dependent time scale rather than $\tau$. In this section, by changing the normalization parameter, we show the effect of St on $\mathcal{C}$ in a one-dimensional pure straining and rotating flows as well as a general arbitrary multi-dimensional flow.

The Lyapunov exponent $\lambda$ is related to the underlying flow in Eq. (30) through $\tilde{\rho}^{\rm Q}$. Its dimensional counterpart, $\tilde{\rho}^{\rm Q_{\rm d}}$ is computed based on $\overline{\|\boldsymbol{S}_{\rm d}\|}$ and $\overline{\|\boldsymbol{\Omega}_{\rm d}\|}$, where subscript d denotes dimensional variables. The Kolmogorov time scale $\tau_{\eta}=\overline{\|\boldsymbol{S}_{\rm d}\|}^{-1}$ thus emerges as the most natural choice for normalization of $\tilde{\rho}^{\rm Q_{\rm d}}$. Employing subscript $\eta$ to distinguish parameters that are normalized based on $\tau_{\eta}$, we have $\tilde{\rho}^{\rm Q}={\rm St}\tilde{\rho}^{\rm Q_{\eta}}$, $\lambda={\rm St}\lambda_{\eta}$, and $\mathcal{C}={\rm St}\mathcal{C}_{\eta}$ and $\Phi={\rm St}^{2}\Phi_{\eta}$. Thus, re-nondimensionalization of Eq. (38) based on $\tau_{\eta}$ yields

[TABLE]

$\lambda_{\eta}$, as oppose to $\lambda$, is often reported in the literature when analyzing the clustering of inertial particles (e.g.,  [35]). Similar to the remarks at the end of the last section, one can analyze the behavior of $\lambda_{\eta}=\lambda_{\eta}({\rm St},\Phi_{\eta},\omega_{\eta})$ under various circumstances. For the sake of brevity, we have condensed this information in a schematic (Figure 5) for cases in which $\omega_{\eta}\approx 0$ and postpone its discussion to Section 4.3, where we compare the results of our one-dimensional model problem to that of a turbulent flow.

The asymptotic predictions in Figure 5, which are obtained from Eq. (39) for a one-dimensional unimodal excitation case, can be generalized to an arbitrary multi-dimensional flow by taking $\Phi_{\eta}=Q_{\eta}/n$. To show this, consider Eq. (30) that holds for the general case. After re-normalization based the Kolmogorov time scale $\tau_{\eta}$

[TABLE]

To obtain the asymptotic behavior of $\lambda_{\eta}$, suppose $\lambda_{\eta}\propto{\rm St}^{p}$ where $p$ is an exponent to be determined. Assuming $p>-1$ for ${\rm St}\ll 1$, then the second term in Eq. (40) can be neglected compared to the first term, and 1 will be the leading order term in the denominator of the integrant. Thus, the integral reduces to $Q_{\eta}$, which implies $\lambda_{\eta}={\rm St}Q_{\eta}/n$ and $p=1$. The analysis for ${\rm St}\gg 1$ is similar. Again assuming $p>-1$, this time the first term can be neglected compared to the second term. Assuming the energetic frequencies of $\tilde{\rho}^{\rm Q_{\eta}}$ are at low frequencies333This assumption, which is in agreement with the results obtained from a turbulent flow in Figure 11-b, implies very heavy inertial particles only respond to the largest structures of the flow since $\omega\sim\dot{x}/L$ with $\dot{x}$ and $L$ being the typical velocity of particle and the size of a flow structure that significantly influences the particle motion, respectively., then the second term in the denominator of the integrant can be neglected. Thus, the integral scales as $({\rm St}\lambda_{\eta})^{-2}Q_{\eta}/n$, implying $\lambda_{\eta}\propto{\rm St}^{-1/2}|Q_{\eta}/n|^{1/4}$ and $p=-1/2$. These results, which are consistent with those of the one-dimensional analysis in Figure 5, are also confirmed by our numerical results in Section 4.

3.2 Numerical validation

To validate the present analysis and compute $\lambda$, we ideally need to compute $\boldsymbol{J}(t)$ as $t\to\infty$ provided that $\mathcal{C}$ (and thus $\lambda$) is defined by Eq. (14). In practice, however, ${\rm det}\left({\boldsymbol{J}}\right)$ grows exponentially in time, producing an ill-conditioned system for a long integration period. To overcome this shortcoming, we infer the long-term response of $\boldsymbol{J}$ to excitation $\partial\boldsymbol{u}/\partial\boldsymbol{x}$ from a single cycle of excitation, as outlined below.

For either the one-dimensional unimodal case discussed in Section 3 or the three-dimensional turbulence that will be discussed in Section 4, $\partial\boldsymbol{u}/\partial\boldsymbol{x}$ will be a harmonic function with period $T$. Evidently, $T=2\pi/\omega$ for the one-dimensional flow. For the turbulent flow, which is chaotic, $T\to\infty$. However, our numerical result shows that if we take $T\gg 1$ for the turbulent flow, i.e., when the dimensional sampling period is much larger than particle relaxation time, the results become independent of $T$. Now, if we denote $\boldsymbol{J}(t=NT)$ by $\boldsymbol{J}^{(N)}$, our goal is to compute ${\rm det}\left({\boldsymbol{J}^{(N)}}\right)$ as the number of cycles $N\to\infty$. A brute force approach of integrating Eq. (15) for many cycles will produce a large numerical error. A better alternative is to compute a transformation matrix $\boldsymbol{A}\in\mathbb{R}^{2n\times 2n}$ that satisfies

[TABLE]

Then, since ${\rm det}\left({\boldsymbol{J}^{(0)}}\right)=1$, ${\rm det}\left({\boldsymbol{J}^{(N)}}\right)$ can be computed as $(\lambda_{1}^{A}\lambda_{2}^{A}\cdots\lambda_{n}^{A})^{N}$, where $\lambda_{1}^{A}$, $\lambda_{2}^{A}$, $\cdots$, $\lambda_{n}^{A}$ are the largest $n$ eigenvalues of $\boldsymbol{A}$. The remaining eigenvalues of $\boldsymbol{A}$ are associated with the decay of the initial conditions and are thus excluded from our calculations. Thus, $\mathcal{C}$, which is $\ln\left({\rm det}\left({\boldsymbol{J}^{(N)}}\right)\right)/(NT)$ as $N\to\infty$, can be computed as $\ln(\lambda_{1}^{A}\lambda_{2}^{A}\cdots\lambda_{n}^{A})/T$.

In practice, we compute $\boldsymbol{A}$ column-by-column by selecting $2n$ linearly independent initial conditions for $[\begin{matrix}\boldsymbol{J}&\boldsymbol{\dot{J}}\end{matrix}]^{\rm T}$ and calculating its time evolution for $T$ using Eq. (15). $\partial\boldsymbol{u}/\partial\boldsymbol{x}$ in Eq. (15) is prescribed based on Eq. (32) for the unimodal excitation case. As discussed in detail in Section 4, it is numerically extracted from the direct numerical simulation for the three-dimensional turbulence by recording its temporal variation along particle trajectories. All calculations are performed in the complex plane using the fourth-order Runge-Kutta time integration scheme. The time step size, which has been verified to be sufficiently small, is $2\pi\times 10^{-6}/\omega$ for the one-dimensional and a fraction of $\tau_{\eta}$ for the three-dimensional calculations. We refer to $\lambda$ (or $\mathcal{C}$) obtained from this procedure as the reference numerical result and use it to validate our analysis.

For the one-dimensional case, the behavior of $\ln({\rm det}\left({\boldsymbol{J}}\right))$ as a function of time varies depending on $\omega$ and $\Phi$ significantly (Figure 6). Based on Eq. (20), the long-term slope of these curves is equivalent to $\lambda$ (or $\mathcal{C}$ since $n=1$). These one-dimensional calculations are also repeated for $n=2$ with non-diagonal $\nabla\boldsymbol{u}$ representing a forced vortex. The results are identical to that of the one-dimensional rotating flow, showing the applicability of Eq. (33) to the higher-dimensional unimodal isotropic flows. Spikes in Figure 6 corresponds to the occurrence of particle crossover when the distance between two particles and ${\rm det}\left({\boldsymbol{J}}\right)$ becomes zero and changes sign. The occurrence of these incidences can be readily deduced from the homogeneous form of Eq. (18). When $F<-1$, then ${\rm det}\left({F}\right)=-F^{2}-F<0$, leading to unstable growth of $F$ to $-\infty$ that indicates ${\rm det}\left({J}\right)=0$. These crossover incidents, which has been called the sling effect in the literature, have been predicted theoretically and shown experimentally in the past [6, 36].

Following the above procedure, $\lambda$ is computed for $\omega\in(0,10]$ and $\Phi\in[-10,10]$ on a $1024\times 1024$ discrete parameter space. These numerical calculations show that $\lambda$ is a highly nonlinear function of $\omega$ and $\Phi$ (Figure 7). For a certain combination of $\omega$ and $\Phi$, $\lambda$ is not differentiable or has a very sharp gradient. These results confirm the possibility of expansion (dispersion) in straining flows, where $\Phi<0$ and ${\rm real}\left({\lambda}\right)>0$. They also confirm the existence of a bound on the contraction-rate that never drops below $-1/2$. The asymptotic behaviors $\lambda\propto|\Phi|^{1/4}$ as $|\Phi|\to\infty$ and $\lambda\propto\Phi$ when $\Phi\ll 1$ are also confirmed by the numerical results (Figure 5). These numerical observations are in full agreement with the prediction of our analysis (Figures 4 and 7).

The numerically computed ${\rm imag}\left({\lambda}\right)$ is also nonzero for a certain combination of $\omega$ and $\Phi$. In those cases, ${\rm imag}\left({\lambda}\right)$ can be either positive or negative depending on the initial condition. Its magnitude, however, is unique and depends on the crossover frequency. According to Figures 4 and 7, there is an agreement between the numerical and analytical results in terms of $\omega$ and $\Phi$ at which ${\rm imag}\left({\lambda}\right)\neq 0$, viz. the blue region in Figure 3. This agreement shows the ability of our analysis in predicting the regimes of particle crossover correctly.

Discrepancies between the numerical and analytical results are also observed, specifically for $\Phi\lesssim-1$ and $\omega\lesssim 1$ where ${\rm imag}\left({\lambda}\right)\neq 0$ (blue region in Figure 3). A higher degree of non-linearity is observed in the numerical result for $\Phi<-(\omega^{4}+10\omega^{2}+9)/16$. Additionally for $\omega>1$, while our analysis predicts ${\rm real}\left({\lambda}\right)=-1/2$ for $-(\omega^{2}+1)^{2}/16<\Phi<-\omega^{2}/4$, the numerical result shows a broader range of $\Phi$ producing ${\rm real}\left({\lambda}\right)=-1/2$. Fitting a curve to the numerical result shows ${\rm real}\left({\lambda}\right)=-1/2$ for $-(\omega^{2}+1)(|\omega|+1/\sqrt{2})^{2}/16<\Phi<-\omega^{2}/4$, which includes $\Phi$ slightly lower than the analytical prediction (Figures 8). In overall, however, the discriminant curves obtained from our analysis provide a good approximation for the values at which the numerical solution changes sign or reaches a plateau.

$\Phi$ and $\omega$ are normalized by the particle relaxation time $\tau$ and thus, are proportional to St. This correspondence implies that $\Phi\gg 1$ and $\omega\gg 1$ represent particles with a high St. As a result, much of the depicted parameter space in Figures 4 and 7 is relevant to high St regime. Our analysis captures the general trend of $\lambda(\Phi,\omega)$ in this part of the parameter space. However, extremely narrow valleys of parameter space with ${\rm real}\left({\lambda}\right)<0$ (e.g. $\Phi=-10$ and $\omega=1$) are missing in the analytical result (Figure 8). The pattern produced by these valleys resembles a fractal structure. The distance between valleys reduces, and they become shallower as $\omega\to 0$. For a given $\omega<1$, these valleys appear only at $\Phi<-(\omega^{2}+1)(|\omega|+1/\sqrt{2})^{2}/16$ (Figure 8). These valleys may override the extremum of ${\rm real}\left({\lambda}\right)$ predicted at $\Phi=-(\omega^{2}+1)^{2}/16$ by our analysis. As discussed in detail in Section 4.3, the predicted extremum agrees with the numerical result when $\omega\to 0$. Neglecting the local extremum associated with the valleys at finite $\omega<1$, the remainder of ${\rm real}\left({\lambda(\omega,\Phi)}\right)$ is a smooth envelope that has an extremum at the predicted $\Phi=-(\omega^{2}+1)^{2}/16$.

The discovery of these narrow valleys has an important implication in designing new hydrodynamic particle separators with high precision. One can generate a pure oscillatory straining flow laden with particles of various density or size and adjust the oscillation frequency such that some particles cluster and others disperse in time. Note that $\omega$ and $\Phi$ in Figures 7 and 8 are normalized by $\tau$ and $\tau^{2}$, respectively. Hence, the dispersion or clustering of particles with different $\tau$ in a given oscillatory pure straining flow is determined based on the value of ${\rm real}\left({\lambda}\right)$ extracted from a parabolic curve in Figures 7 and 8. This parabolic curve can intersect with the valleys, which occur in a part of parameter space where $\Phi$ is a quartic function of $\omega$. By adjusting the oscillation frequency of the flow, the intersection point can be changed such that $(\omega,\Phi)$ of a pre-specified particle class coincides with a valley. As a result, while that pre-specified particle class clusters over time, particles with slightly higher or lower $\tau$ will disperse as their $(\omega,\Phi)$ is outside of the valley. To demonstrate this idea, we simulated a 2D pure oscillatory straining flow laden with two class of particles. Particles in these two classes are different in size merely by 1%. The flow parameters are adjusted according to the above procedure to have $(\omega,\Phi)$ of particles inside and outside of the valley that passes through $(\omega,\Phi)=(1,-10)$. An animation in the supplementary material shows the response of these two classes of particles (red and black) to the flow. Only after a few oscillations, one class of particles clusters and accumulates toward the flow stagnation point whereas the other class disperses and gets removed from the oscillatory region of the flow. In practice, neither the drag law or the flow field is perfectly linear. The particles have a finite size and density and may interact with each other or the solid boundaries of the apparatus. Although our preliminary numerical results show that the existence of these narrow valleys is robust against these imperfections, an experimental study is still needed to confirm the feasibility of this novel concept for particle separation in the future.

3.3 Comparison against other models

The present analysis is a more general form of our former [32] and Maxey’s [24] analysis. For the case of one-dimensional unimodal excitation, SL is simplified by combining Eqs. (5), (24), and (33) as

[TABLE]

and similarly, RM is simplified by combining Eqs. (4), (24), and (33) as

[TABLE]

These two relationships can be derived from Eq. (38) or Eq. (36) under the assumptions of $|\lambda|\ll 1$ and $|\omega|\ll 1$ as well, therefore they can be considered a special case of the present analysis.

RM and SL are valid only within the regimes compatible with their underlying assumptions (Figure 9). RM is derived for small St, which translates to $\omega\ll 1$ and $|\Phi|\ll 1$ (unless $|\Phi|\approx 0$, where there is no flow and $\lambda\approx 0$). SL is derived for small ${\rm real}\left({\lambda}\right)$, which translates to $\omega\gg 1$ or $|\Phi|\ll 1$. The present analysis, on the other hand, provides a reasonable approximation at all $\omega$ and $\Phi$. To highlight these differences, we have compared the predictions of all models against the reference numerical results in Figures 9 and 10. Figure 10(c) confirms RM and SL predictions are invalid at small $\omega$, except for $\Phi\ll 1$. For example, these two relations predict $\lambda=-5$ at $\omega\approx 0$ and $\Phi=-5$, whereas $\lambda=0.3$ and 0.29 from the numerical calculations and our analysis, respectively. At higher $\omega$, SL prediction approaches the numerical solution for a wide range of $|\Phi|$, while RM prediction is wrong everywhere except $|\Phi|\approx 0$ (Figure 10(e) and Figure 9).

The nonlinear behavior of ${\rm real}\left({\lambda}\right)$ at high St regime is captured only by our analysis. According to Eqs. (42) and (43), RM and SL are both linear functions of $\Phi$. Hence, they both fail to predict the nonlinear behavior of ${\rm real}\left({\lambda}\right)$ at higher $|\Phi|/\omega$ (Figures 7 and 10). Neither RM or SL predicts the lower bound on the contraction rate ($\min(\lambda)\geq-0.5$), possibility of expansion in straining flows ($\lambda>0$ for $\Phi\ll-\omega^{4}$), and the asymptotic variation of expansion or contraction rate at high amplitude oscillations ($\lambda\propto|\Phi|^{1/4}$ for $|\Phi|\gg 1$). These two models also fail to predict particle crossovers that occur in strong-straining flows where ${\rm imag}\left({\lambda}\right)\neq 0$.

4 Three-dimensional isotropic turbulence: extension to multimodal excitation

In a physically realistic turbulent flow, excitation is not at a single frequency but involves a continuous range of frequencies. To analyze a multimodal excitation, all modes that appear in Eq. (30) must be retained. With the transformation introduced in Eq. (35), Eq. (30) can be expressed as

[TABLE]

which provides a generic relationship for $\gamma^{2}$ in the form of an eigenvalue problem. An iterative approach can be adopted for computing $\gamma^{2}$ from this equation, as the evaluation of the integral requires prior knowledge of $\gamma^{2}$.

Under certain conditions, it is possible to find an explicit closed-form relationship for $\gamma$ based on Eq. (44). One such scenario is when the explicit form of $\tilde{\rho}^{\rm Q}(\omega;{\rm St})$ is known, as was the case in Section 3. If the analytical form of $\tilde{\rho}^{\rm Q}(\omega)$ were to be available and $\tilde{\rho}^{\rm Q}\to 0$ as $|\omega|\to\infty$, Cauchy’s integral formula could be employed to express

[TABLE]

Depending on the form of $\tilde{\rho}^{\rm Q}(\hat{i}\gamma;{\rm St})$, an explicit relationship for $\gamma$ can be obtained from this expression.

Another possible scenario that may arise is a design problem formulated as finding a specific $\tilde{\rho}^{\rm Q}(\omega;{\rm St})$ when a desirable $\gamma({\rm St})$ is given. Designing a hydrodynamic particle separator by enhancing the clustering of a particular class of particles in a polydisperse distribution is an instance of such a scenario. The present formulation can be instrumental in solving this inverse problem by expressing Eq. (44) as a Fredholm integral equation of the first kind with a kernel function $(\gamma^{2}+\omega^{2})^{-1}$ [37].

4.1 $\tilde{\rho}^{\rm Q}$ in a turbulent flow

In a nutshell, $\tilde{\rho}^{\rm Q}$ determines at what frequencies different gradients in the fluid velocity field oscillate as seen by the particle. Analytical exploration of $\tilde{\rho}^{\rm Q}$ is a non-trivial task for an arbitrary flow. In general, fluid flows are too complex to be solved analytically and expressed in a closed-form solution. This statement is particularly true about turbulence that is typically studied experimentally or numerically when a detailed solution is needed. What makes the analytical calculation of $\tilde{\rho}^{\rm Q}$ even more challenging in such flows is its dependence on the particle trajectory. Besides having an explicit relationship for the velocity gradient at all points in space and time, one needs to know how particles preferentially sample the flow to compute $\tilde{\rho}^{\rm Q}$. Despite all these complexities, analytical modeling of $\tilde{\rho}^{\rm Q}$ in an approximate form could be a feasible task to be accomplished by future investigations. An essential component in such an investigation would be a model for the velocity gradients along particle trajectories. Such an effort has been undertaken in the context of large-eddy simulations [38] for tracers and needs to be extended to the inertial particles. In this study, we do not attempt to model $\tilde{\rho}^{\rm Q}$, but rather compute it directly from a numerical simulation. This way, the accuracy of the present analysis can be evaluated independently of the accuracy of the model used for $\tilde{\rho}^{\rm Q}$.

Although the results presented in this section are based on the accurate quantification of $\tilde{\rho}^{\rm Q}$, in Section 4.3 we will show that such an accurate knowledge is not necessary for a qualitative prediction of $\mathcal{C}$. Even with a minimal knowledge of the flow, i.e., the fact that particles preferentially sample the straining region of the flow or $\langle Q_{\eta}\rangle<0$, the present analysis can predict the existance of a dip in $\langle\mathcal{C}_{\eta}\rangle$ at $\rm{St}\approx 1$ and $\langle\mathcal{C}_{\eta}\rangle\to 0$ as $\rm{St}\to 0$ or $\infty$. This lack of sensitivity to the form of $\tilde{\rho}^{\rm Q}$ shows that the clustering phenomenon is rather universal, reaffirming our intuition that the clustering ought to occur primarily at $\rm{St}\approx 1$ regardless of the details of the flow field.

To generate the background flow, we perform direct numerical simulation of a triply periodic homogeneous isotropic turbulence using an in-house solver with a specialized linear solver [39]. A second-order spatial discretization on a $256^{3}$ grid and 4th order Runge-Kutta time integration scheme are employed. Stationary turbulence is maintained by adding a forcing term to the momentum equation that is proportional to the velocity [40]. The forcing term is dynamically computed at each time step to prevent fluctuation of $\tau_{\eta}$ and thus St [32, 41]. The maximum deviation of $\tau_{\eta}$ from the target Kolmogorov time scale is 0.3%. The Reynolds number based on the Taylor micro-scale is Re${}_{\lambda}=100$. Special care has been taken in interpolating quantities at the location of particles from the Eulerian grid. In particular, the interpolation scheme is designed to correctly translate the incompressibility condition to the Lagrangian velocity gradient tensor. Additionally, the Lagrangian gradients are kept $C^{0}$ continuous by interpolating from a pre-constructed continuous Eulerian field [32].

The particles trajectory is computed using Eq. (2). A dilute mixture is considered so that the effect of particles on the fluid can be neglected, thereby eliminating potential errors that can be introduced by two-way coupling forces [42, 43, 44]. 113 classes of particles are considered in total with ${\rm St}=2^{p/8}$, $p\in\{-32,\ldots,80\}$. At each Stokes number, approximately $10^{4}$ randomly seeded particles were simulated for several large eddy turnover time to allow development of clusters. Starting with this time-evolved distribution, we record the velocity gradient tensor at the position of each particle for 700$\tau_{\eta}$ with $0.1\tau_{\eta}$ intervals. Based on $\nabla\boldsymbol{u}_{\eta}(t)$, $\tilde{\rho}^{\rm Q_{\eta}}$ is computed using Eq. (6) at each St. The number of particles and the integration period are verified to be sufficient for achieving statistical convergence. The results of these calculations at few Stokes numbers are shown in Figure 11.

In a turbulent flow, in contrast to the previous case study, $\nabla\boldsymbol{u}_{\eta}$ has full rank, and $\|\boldsymbol{S}_{\eta}\|$ and $\|\boldsymbol{\Omega}_{\eta}\|$ are nonzero simultaneously. These parameters, computed along the trajectory of particles, exhibit complex behaviors due to their dependence on St, caused by the preferential sampling of the flow field by the particles. These dependencies are briefly mentioned here for the sake of completeness and discussed in detail in [32]. In the homogeneous turbulence under consideration, $\rho^{\rm S_{\eta}}(t;{\rm St})$ and $\rho^{\rm\Omega_{\eta}}(t;{\rm St})$ (the norm of the autocorrelation function of the strain- and rotation-rate tensors, respective) are both strictly positive. $\langle\rho^{\rm S_{\eta}}\rangle$ is relatively independent of St and exponentially decays with time, which is analogous to what can be observed with particles in a random straining flow. $\langle\rho^{\rm\Omega_{\eta}}\rangle$, on the other hand, varies significantly versus St. As a result, $\langle\rho^{\rm Q_{\eta}}\rangle=\langle\rho^{\rm\Omega_{\eta}}\rangle-\langle\rho^{\rm S_{\eta}}\rangle$ strongly depends on St.

For ${\rm St}\lesssim 1$, $\langle\rho^{\rm Q_{\eta}}\rangle$ undergoes an increasing-decreasing trend in time. $\langle\rho^{\rm Q_{\eta}}\rangle$ being negative at $t_{\eta}\ll 1$ and ${\rm St}\ll 1$ is due to the smaller value of $\langle\rho^{\rm\Omega_{\eta}}\rangle$. Hence, particles with small St tend to centrifuge out of rotational regions with the short time constant and follow slow vortical features since $\langle\tilde{\rho}^{\rm Q_{\eta}}\rangle>0$ at $\omega_{\eta}\ll 1$ and ${\rm St}<1$ (Figure 11).

For ${\rm St}\gg 1$, particles are not responsive to the velocity fluctuations and follow a trajectory that is uncorrelated with the flow. As a result, particles distribute uniformly in space and the Lagrangian and Eulerian statistics become almost identical. Additionally, one can show that the Eulerian strain-rate and rotation-rate autocorrelation functions are equal in a periodic domain. Therefore, $\langle\rho^{\rm S_{\eta}}\rangle$ and $\langle\rho^{\rm\Omega_{\eta}}\rangle$ converge to the same value, leading to $\langle\rho^{\rm Q_{\eta}}\rangle\to 0$ and $\langle\tilde{\rho}^{\rm Q_{\eta}}\rangle\to 0$ as ${\rm St}\to\infty$.

In Section 3.2, we demonstrated that the particle clouds only contract in a straining regime as a rotating regime only leads to cloud expansion. Thus, in turbulence, where regions of higher rotation-rate and strain-rate coexist in space, particles tend to accumulate in regions of higher strain-rate. The preferential concentration of particles in the straining regions, which occurs at all St, is supported by the dominance of $\langle\rho^{\rm S_{\eta}}\rangle(t=0)$ over $\langle\rho^{\rm\Omega_{\eta}}\rangle(t=0)$. This dominance is most noticeable at ${\rm St}\approx 1$, where $\langle\rho^{\rm Q_{\eta}}\rangle(0)$ has an extremum. At large and small Stokes numbers, on the other hand, $\langle\rho^{\rm Q_{\eta}}\rangle(0)$ asymptotes to zero as $\langle\rho^{\rm S_{\eta}}\rangle(t=0)$ and $\langle\rho^{\rm\Omega_{\eta}}\rangle(t=0)$ reach a balance.

4.2 Clustering statistics

From $\tilde{\rho}^{\rm Q_{\eta}}$, $\mathcal{C}^{t}_{\eta}$ is computed for RM, SL, and the present analysis. The method described in Section 3.2 is employed to compute the reference quantities. To improve accuracy, computations are performed with 12 time steps between consecutive records of $\nabla\boldsymbol{u}_{\eta}$, i.e. $\Delta t=\tau_{\eta}/120$. To prevent ill-conditioning of $\boldsymbol{A}$ in Eq. (41), the time integration period is limited to 12.5$\tau$ when $\tau<56$, whereas the entire interval of 700$\tau_{\eta}$ is considered for $\tau\geq 56$. This division to sub-intervals is uniformly applied to the computation of the input to the models and also the calculation of the reference results to ensure a one-to-one comparison between the two. To construct the PDF of $\mathcal{C}^{t}_{\eta}$, each sub-interval associated with each $10^{4}$ particles is treated as an independent ensemble (Figure 12). The overall trend is similar to what was observed in Section 3.3 for the one-dimensional case with the present analysis being the closest to the reference followed by SL and RM. The only exception is ${\rm St}=1$ where the PDF from the present analysis is skewed and shows a second unphysical peak at $\mathcal{C}^{t}_{\eta}\approx-0.5$.

The ensemble-averaged of $\mathcal{C}^{t}_{\eta}$ is computed from the PDFs (Figure 13). Since RM and SL are a linear function of $\langle\tilde{\rho}^{\rm Q_{\eta}}\rangle$, these ensemble-averaged quantities can be computed directly from the results shown in Figure 11. For the present analysis, however, $\langle\mathcal{C}_{\eta}(\tilde{\rho}^{\rm Q_{\eta}})\rangle$ is slightly different from $\mathcal{C}_{\eta}(\langle\tilde{\rho}^{\rm Q_{\eta}}\rangle)$. Their difference depends on the integration period $t_{\eta}$ and asymptotes to zero as $t_{\eta}\to\infty$.

The accuracy of all three models in predicting the reference results is similar to what was observed for the one-dimensional cases. All models collapse with the reference for ${\rm St}\ll 1$, whereas for ${\rm St}\gg 1$, their prediction widely varies. Among the three models, only the present analysis captures the expansion of clouds at ${\rm St}\gg 1$. ${\rm real}\left({\langle\mathcal{C}_{\eta}\rangle}\right)>0$ is predicted at ${\rm St}\gg 1$ despite the fact that $\langle\tilde{\rho}^{\rm Q_{\eta}}\rangle<0$ at all frequencies (Figure 11). Prediction of expansion in a straining regime stems from the nonlinear behavior of $\lambda(\Phi)$ that appeared as ${\rm real}\left({\lambda}\right)>0$ for $\Phi<-(\omega^{4}+10\omega^{2}+0)/16$ in Figure 4.

Despite capturing the overall trend, the present analysis is not in full quantitative agreement with the reference at high St. The disagreement can be attributed to the assumption of single $\lambda$ in Eq. (23), where a full rank matrix was replaced with a diagonal matrix444Our earlier argument that all $\lambda_{i}$ are the same in an isotropic flow relies on a sufficiently long sampling period to diminish statistical differences between $\lambda_{i}$. Here, however, the sampling period is limited to 12.5$\tau$ to prevent ill-conditioning of $\boldsymbol{A}$, resulting in 3 distinct $\lambda_{i}$ for each ensemble, which often are significantly different., and excluding sub-harmonics from our analysis when deriving Eq. (26).

There are two predictions that are unique to the present analysis and in agreement with the reference results. The first is the asymptotic behavior of $\langle\mathcal{C}_{\eta}\rangle$ as ${\rm St}\to\infty$ and the second is the sequence of the onset of crossover and dispersion. At very large Stokes number, $\langle\mathcal{C}_{\eta}\rangle\propto\rm{St}^{-1/2}$ (the right inset in Figure 13), which is in agreement with the prediction of the present analysis (Figure 5). Additionally, the present model predicted crossovers to occur at a Stokes number that precedes the Stokes number at which clustering is transitioned to dispersion regime. Based on Figure 13, the onset of crossover and dispersion occur at $\rm St\approx 0.18$ and $\rm St\approx 1.6$, respectively, following the predicted sequence. The ratio of two Stokes number, however, does not quantitatively agree with our model ($\approx 9$ versus 3).

To show the effect of turbulence intermittency on the particle clustering, we computed the second moment of $\mathcal{C}^{t}_{\eta}$ (Figure (14)). The plotted moments are normalized by $\sqrt{t_{\eta}}$ to ensure their independence from sampling period $t$. Among the available models, the present analysis provides the best estimation for $(\mathcal{C}_{\eta}^{t})^{\prime}$. All models collapse with the direct computations at the limit of small St and deviate from it as St increases. For ${\rm St}>1$, RM and SL predictions linearly increase with St while the present analysis remains bounded. A slope of $-1/3$ is observed in this regime of Stokes number, indicating a reduction in sensitivity of particle to turbulence fluctuations, which can be explained by their higher inertia and smoother trajectories [45].

4.3 Relevance of the one-dimensional model problem

In section 3, we considered a very simple flow that was one-dimensional and oscillated at a single frequency. The results of that analysis were primarily presented as a function $\Phi$ and $\omega$ that characterize the amplitude and frequency of oscillations, respectively, when normalized based on the particle relation time $\tau$. In Section 3.1, we showed that those results become $\rm St$-dependent if they are re-normalized based on the flow time-scale $\tau_{\eta}$. Later in Section 4.2, we applied a more general version of the same analysis to a three-dimensional isotropic turbulent flow and presented the results that were also normalized based on the flow time-scale $\tau_{\eta}$. Our goal in this section is to find how the result of these two cases compare in terms of St and whether the simple one-dimensional flow can adequately lead to an understanding of the clustering phenomenon in the three-dimensional turbulent flow.

Let us revisit the one-dimensional flow of Section 3 in the regime of $\omega\ll 1$. The condition of $\omega\ll 1$ corresponds to a flow oscillating at a frequency much lower than the inverse of particle relaxation time $\tau$. Based on the discussion of Section 3.1, $\lambda_{\eta}$ for this flow can be computed from Eq. (39) with $\omega_{\eta}\approx 0$. Since $\lambda_{\eta}$ can be replaced by $\mathcal{C}_{\eta}$ for a one-dimensional flow, Eq. (39) can be written as

[TABLE]

where $\Phi_{\eta}=+1$ and $-1$ represent a rotating and straining flow, respectively555Note that $\Phi_{\eta}=\overline{\|\boldsymbol{\Omega}_{\eta}\|^{2}-\|\boldsymbol{S}_{\eta}\|^{2}}$ (see Eq. (32)). Thus, $\Phi_{\eta}=\pm 1$ since $\tau_{\eta}$ is taken as $\|\boldsymbol{S}_{\rm d}\|^{-1}$ or $\|\boldsymbol{\Omega}_{\rm d}\|^{-1}$, where $\boldsymbol{S}_{\rm d}$ and $\boldsymbol{\Omega}_{\rm d}$ are the dimensional strain- and rotation-rate tensors, respectively.. For this one-dimensional flow in which $\omega_{\eta}\approx 0$, the prediction of SL collapses with that of RM to $\mathcal{C}_{\eta}=\pm{\rm St}$. The reference results can also be obtained by following the procedure of Section 3.2 using $\omega_{\eta}=10^{-4}$. The result of all these calculation is plotted in Figure 15, where ${\rm imag}\left({\mathcal{C}_{\eta}}\right)$ is also shown as the measure of particle crossover frequency.

As mentioned earlier, the accumulation of particles in the straining region of a turbulent flow leads to $\langle\|\boldsymbol{S}_{\eta}\|\rangle>\langle\|\boldsymbol{\Omega}_{\eta}\|\rangle$. As a result, the behavior of $\mathcal{C}_{\eta}({\rm St})$ in the isotropic turbulence (Figure 13) should be compared against the one-dimensional straining flow (black curves in Figure 15) rather than the rotating flow. One can verify that there is a one-to-one correspondence between the two cases. Increasing St from zero, ${\rm real}\left({\mathcal{C}_{\eta}}\right)$ decreases till the onset of crossover. The trend is reversed once ${\rm imag}(\mathcal{C}_{\eta})\neq 0$ up to a Stokes number at which $\mathcal{C}_{\eta}>0$. For larger St, ${\rm real}\left({\mathcal{C}_{\eta}}\right)$ changes non-monotonically with $\mathcal{C}_{\eta}\propto{\rm St}^{-1/2}$ as ${\rm St}\to\infty$ (not shown in Figure 15 explicitly). Such a close similarity between one- and three-dimensional cases provides a window from which one can view how clustering phenomenon unfolds as the Stokes number changes in a flow. We provide this description for the simpler one-dimensional case that is easier to understand, however, one can extend the following arguments to the turbulent flow that in essence is three oscillatory one-dimensional straining or rotating flows acting on a particle cloud (along the principal directions of the velocity gradient tensor) as it traverses the flow.

To better understand the behavior of particles in the straining flow (black curves in Figure 15), it is necessary to distinguish between two regimes in which ${\rm real}\left({\mathcal{C}_{\eta}}\right)$ decreases and increases with St. In the first regime occurring at ${\rm St}<\sqrt{1/32}$666This transition is predicted to occur at ${\rm St}\approx\sqrt{1/16}$ rather than $\sqrt{1/32}$ by the present analysis despite the fact that the corresponding extremum of ${\rm real}\left({\lambda}\right)$ was correctly predicted in Section 4. A closer examination shows that the extremum of ${\rm real}\left({\mathcal{C}_{\eta}}\right)$ is shifted to higher St due to the re-normalization of Eq. (46) and under-prediction of ${\rm real}\left({\lambda}\right)$ by our analysis (note the extremum of ${\rm real}\left({\mathcal{C}_{\eta}({\rm St})}\right)$ occurs at ${\rm St}=\sqrt{-\Phi}$ with $\Phi$ satisfying ${\rm real}\left({\lambda(\Phi)}\right)=2\Phi[{\rm d}({\rm real}\left({\lambda}\right))/{\rm d}\Phi]$)., ${\rm imag}\left({\mathcal{C}_{\eta}}\right)=0$ and ${\rm real}\left({\mathcal{C}_{\eta}}\right)$ decrease almost linearly versus St. No particle crossover occurs in this regime and the increase in the particle inertia is met with proportionally stronger slippage, leading to a faster rate at which particles get close to each other ($\{{\rm real}\left({\lambda}\right)<0,\;{\rm imag}\left({\lambda}\right)=0\}$ in Figure 2). Further increase in St leads to the second regime, in which particles have enough inertia to cross over each other (${\rm imag}\left({\lambda}\right)\neq 0$ in Figure 2). The maximum rate of convergence is obtained at the onset of crossovers before particles begin to overshoot each other. For ${\rm St}>\sqrt{1/32}$, the relative velocity of particles at the moment of crossover increases with St, such that at ${\rm St}\geq 3/4$ their mean distance, rather than decreasing, begins to increase over time. As St is increased beyond ${\rm St}\approx 3.6$, the rate at which particles diverge decreases, i.e., ${\rm d}({\rm real}\left({\mathcal{C}_{\eta}}\right))/{\rm dSt}<0$ for ${\rm St}>3.6$. Due to the high inertia of particles at this limit, particles hardly respond to the oscillations of the underlying flow and as a result tend to maintain their initial position, lowering their divergence rate. Another consequence of having very inertial particles is fewer incidents of crossovers that leads to a decrease in ${\rm imag}\left({\mathcal{C}_{\eta}}\right)$.

In contrast to the one-dimensional straining flow, which produced a clustering phenomenon analogous to the isotropic turbulent flow, it is hard to find a real-world example that corresponds to the one-dimensional rotating flow (red curves in Figure 15). The reason is that the most commonly studied spatially and temporally oscillating particle-laden flows are turbulent. These flows contain regions in space where $\|\boldsymbol{S}\|<\|\boldsymbol{\Omega}\|$ and $\|\boldsymbol{S}\|>\|\boldsymbol{\Omega}\|$. Inertial particles in these flow tend to accumulate in regions where $\|\boldsymbol{S}\|>\|\boldsymbol{\Omega}\|$, thus leading to a behavior that is similar to that of the straining flow. If one were to identify a spatially and temporally oscillating flow that on average exposes particles to a higher rotation- than strain-rate (i.e., $Q>0$), then we may observe a behavior similar to that of the one-dimensional rotating flow in Figure 15. The behavior of particles in such rotating flow will be less complicated as there is no particle crossover. The distance between particles always increases over time at $\mathcal{O}({\rm St})$ rate at ${\rm St}\ll 1$ and $\mathcal{O}({\rm St}^{-1/2})$ at ${\rm St}\gg 1$ (if dependence of $Q_{\eta}$ on St is neglected). The asymptotic behavior of particles at these two limits will remain similar to the straining flow. At small St, particles follow fluid tracers and their distance barely changes over time. At high St, particles barely respond to the underlying flow oscillations and maintain their initial position777Since ${\rm real}\left({\mathcal{C}_{\eta}({\rm St})}\right)$ is linear at low St, the magnitude of ${\rm real}\left({\mathcal{C}_{\eta}}\right)$ for straining and rotating flows is equal up to the leading order term. At large St, however, the leading order terms have a similar exponent, i.e. ${\rm real}\left({\mathcal{C}_{\eta}({\rm St})}\right)\propto{\rm St^{-1/2}}$, but different magnitude with the rotating flow leading to a larger magnitude (Figure 15)..

There are also some differences between the results of the one-dimensional straining flow and the three-dimensional turbulence. In the turbulent flow, the maximum contraction rate occurs at a higher St and is less significant (Figures 13 and 15). This weaker clustering can be attributed to the presence of rotating regions that are absent in the one-dimensional case. Additionally, the asymptotic behavior of ${\rm real}\left({\mathcal{C}_{\eta}}\right)$ at ${\rm St}\ll 1$ is linear for the one-dimensional case, whereas it is superlinear for the turbulent flow. Since ${\rm real}\left({\langle\mathcal{C}_{\eta}\rangle}\right)\approx{\rm St}\langle Q_{\eta}\rangle$ at small St, this difference is a result of the behavior of $\langle Q_{\eta}\rangle$ versus St. $\langle Q_{\eta}\rangle$ was independent of St in the one-dimensional case, hence the linear rate, whereas it is proportional to St in the turbulent case888The rate at which $\langle Q_{\eta}\rangle$ grows versus St in an isotropic turbulent flow can be an artifact of the periodic boundary condition imposed for computational consideration. This potential artifact caused by the spatial confinement must be removed in the future for a more realistic assessment of the asymptotic behavior of ${\rm real}\left({\langle\mathcal{C}_{\eta}\rangle}\right)$ at small St in a turbulent flow., hence the superlinear rate. Note the differences in the magnitude of ${\rm imag}\left({\mathcal{C}_{\eta}}\right)$ is immaterial since the three-dimensional results are normalized by their maxima. Finally, the Stokes at which clustering is maximized coincides with the onset of crossovers for the one-dimensional case, whereas it occurs after the onset of crossovers for the turbulence case. This difference could be a result of the turbulence intermittency as ${\rm imag}\left({\mathcal{C}_{\eta}}\right)\neq 0$ for a fraction of ensembles leads to ${\rm imag}\left({\langle\mathcal{C}_{\eta}\rangle}\right)\neq 0$, thus shifting the onset to smaller St where only the trajectories of a few particles cross.

In overall, the one-to-one correspondence between the results of the one- and three-dimensional cases is remarkable (Figure 13 versus black curves in Figure 15). The former is a simple flow expressed by Eq. (32) oscillating at a single frequency, while the latter involves a multiscale three-dimensional chaotic flow with gradients oscillating at a continuous spectrum of frequencies. Such a close comparison between the two is a testimony to the fundamental significance and relevance of the one-dimensional model problem. It confirms our earlier hypothesis that the clustering phenomenon and how it varies versus St is primarily determined based on the governing equations of the motion of particles rather than the detailed structure of the underlying flow.

5 Conclusions

We derived a solution (Eq. (30)) for the Lyapunov exponents of inertial particles subjected to oscillatory fluid motion. Our analysis is aimed at predicting the rate of expansion or contraction of clouds of inertial particles, and also their crossovers frequency. We employed the sum of the Lyapunov exponents, i.e., the rate of change of volume of a cloud of particles in three dimensions, to characterize regimes of preferential concentration. We showed that our solution is more general and reproduces the pre-established models in the literature ([24] and [32]). Consistent with the previous models, the only flow-related parameter that appears in our model is the difference between the spectrum of rotation and strain rate tensors norm, viz. a closely related parameter to the Q-criterion, underscoring its fundamental role in clustering of inertial particles. We employed a canonical setup with unimodal excitation to investigate the behavior of the Lyapunov exponent under a wide range of flow conditions. Only the expansion with no crossovers was observed in a rotating regime, whereas both the contraction and expansion with the possibility of crossovers was observed in a straining regime. In a straining regime, the expansion and crossover occur for a sufficiently large oscillation amplitude. Additionally, a $-1/2$ bound on the rate of contraction (normalized by the particle relaxation time) was found. Our analysis also showed the Lyapunov exponent is linearly proportional to the Q-criterion at low oscillation amplitude and its power of $1/4$ at high oscillation amplitude. These observations, which are confirmed by the numerical simulations, are not captured by the other models. Other available models capture only the linear regime, where the oscillation amplitude is small. Discrepancies were also observed between our analysis and the reference results. Neglecting the contribution of the sub- and super-harmonics in our formulation are deemed to be the primary sources of discrepancies. Accounting for higher-order terms and extending the present analysis to statistically anisotropic flows remains as topics for future studies.

Following this canonical setting, we extended our analysis to complex multi-dimensional flows, in which a continuous range of frequencies is present. We considered a three-dimensional isotropic forced turbulence for validation of our analysis. Despite the added complexity, this case produced results analogous to that of the one-dimensional straining regime with unimodal excitation. In all cases, the contraction rate was proportional to St at small St with an extremum around the onset of crossovers. For larger St, the rate of contraction decreases till net expansion is observed at ${\rm St}\approx 1$. Beyond that St, the rate of expansion reaches a maximum and then asymptotes to zero proportional to $\rm St^{-1/2}$ as ${\rm St}\to\infty$. While all models correctly capture the linear trend at low St, only the present analysis provides a good prediction of the subsequent nonlinear behaviors at higher Stokes numbers. Additionally, only our analysis captures the occurrence of particle crossovers at high St and also the non-monotonic variation of the standard deviation of the rate of expansion or contraction versus St.

acknowledgments

We gratefully thank Prof. Andreas Acrivos for his extensive effort in reading this manuscript and providing detailed comments, which have been quite instrumental in clarification of the introduced concepts. This work was supported by the United States Department of Energy under the Predictive Science Academic Alliance Program 2 (PSAAP2) at Stanford University.

Appendix A Derivation of Eq. (14)

To relate $\mathcal{C}$ to $\mathcal{C}^{t}$, the Eulerian form of Eq. (4) must be expressed in terms of Lagrangian quantities. Since $\nabla$ operator can be expressed as $\partial/\partial x_{i}=(\partial/\partial X_{j})(\partial X_{j}/\partial x_{i})$, we have

[TABLE]

From Jacobi’s formula

[TABLE]

and as a result

[TABLE]

From Eqs. (49) and (13)

[TABLE]

which in combination with Eq. (4) gives

[TABLE]

completing the derivation.

Appendix B Derivation of Eq. (20)

To relate $\mathcal{C}$ to $\boldsymbol{F}$, we start from the Jacobi formula in Eq. (19) that can expressed as

[TABLE]

Averaging Eq. (52) over time produces

[TABLE]

Thus

[TABLE]

From Eq. (13) and ${\rm det}\left({\boldsymbol{J}(0)}\right)=0$, the above equation can be simplified to

[TABLE]

Based on Eq. (14), the term on the right-hand side of Eq. (55) is $\mathcal{C}$, thus completing the proof.

Appendix C Higher order expansions

To explain the difference between the analytical prediction and the numerical result observed in Section 3.2, we need to revisit the underlying assumptions of Eq. (38). Equation (38) is an exact solution of Eq. (30) that was obtained from Eq. (18) when the higher order terms in Eq. (25) are neglected. To account for those higher order effects, we consider

[TABLE]

as an asymptotic form of $\boldsymbol{F}$, where $\boldsymbol{\epsilon}_{1}$ corresponds to $\boldsymbol{\Psi}$, which is the first order solution obtained in Section 2. In this expression, the higher order terms, i.e. $\boldsymbol{\epsilon}_{i}$ for $i>1$, can be computed such that the second summation in Eq. (25) is represented more accurately. The addition of higher order terms modifies $\lambda$ in two ways. The first is to interact with the lower order terms and alter their amplitude. The second is to directly contribute to $\lambda$ through contraction of $\boldsymbol{\epsilon}_{i}(\omega)$ and $\boldsymbol{\epsilon}_{i}(-\omega)$. In either case, $\|\boldsymbol{\epsilon}_{i}\|$ provides a measure of the significance of those higher order terms. Therefore, we evaluate the importance of these higher order terms by investigating the decay rate of $\|\boldsymbol{\epsilon}_{i}\|$ versus $i$.

The first three terms of the asymptotic expansion are derived and provided in Table 1. As the order of expansion $n$ increases, additional terms appear in the lower order terms as a result of the interaction between higher order terms. All these additional terms are produced by the nonlinear term in Eq. (18) that turns into a convolution in the Fourier space. Further analysis of these terms shows that the decay rate of $\|\boldsymbol{\epsilon}_{i}\|$ depends primarily on $|\Phi|/\omega$ in the one-dimensional unimodal excitation setting. In general $\|\boldsymbol{\epsilon}_{i}\|\propto(|\Phi|/\omega)^{i}$ if $|\Phi|/\omega<1$ and in the worst-case scenario $\|\boldsymbol{\epsilon}_{1}\|\approx\|\boldsymbol{\epsilon}_{2}\|\approx\cdots\approx\|\boldsymbol{\epsilon}_{n}\|$. In the latter case, the asymptotic form in Eq. (23) will not converge. To show this behavior in practice, $\boldsymbol{\epsilon}_{i}$ is derived by continuing Table 1 beyond $n=3$. $\boldsymbol{\epsilon}_{i}$ and $\lambda$ are then calculated iteratively for several values of $\omega$ and $\Phi$. The decay rate of $\|\boldsymbol{\epsilon}_{i}\|$ versus $i$ is then computed and shown in Figure 16. This figure confirms that $\|\boldsymbol{\epsilon}_{i}\|$ may not decay monotonically if $|\Phi|/\omega\gg 1$, as is the case of $\Phi=-10$ and $\omega=0.1$. Non-converging $\|\boldsymbol{\epsilon}_{i}(\omega,\Phi)\|$ is particularly observed where the prediction of Eq. (30) disagrees with the reference numerical result.

Neglecting the higher order terms accounts for only some of the observed discrepancies. For a certain combination of $\omega$ and $\Phi$, including higher order terms does not produce a better estimate for $\lambda$. In these cases, $\|\boldsymbol{\epsilon}_{i}\|$ may not even decay for $i>1$, suggesting that the assumed form in Eq. (23) is incomplete. This form assumes that excitation at $\omega$ would generate a solution oscillating at $\omega,2\omega,3\omega,\cdots$, producing only super-harmonics. The numerical simulation shows, however, that the solution may oscillate at lower frequencies, producing sub-harmonics. A closer examination shows that the amplitude of sub-harmonics can far exceed that of super-harmonics. It is these sub-harmonics that cause highly nonlinear behaviors that are missing in the prediction of the present analysis. Identifying the nontrivial patterns of these sub-harmonics and thereby accounting for their contribution to $\lambda$ remains as a topic for future studies.