Statistical model for the orientation of non-spherical particles settling in turbulence

The orientation of small anisotropic particles settling in a turbulent fluid determines some essential properties of the suspension. We show that the orientation distribution of small heavy spheroids settling through turbulence can be accurately predicted by a simple Gaussian statistical model that takes into account particle inertia and provides a quantitative understanding of the orientation distribution on the problem parameters when fluid inertia is negligible. Our results open the way to a parameterisation of the distribution of ice-crystals in clouds, and potentially leads to an improved understanding of radiation reflection, or particle aggregation through collisions in clouds.

The orientation of small anisotropic particles settling in a turbulent fluid determines some essential properties of the suspension. We show that the orientation distribution of small heavy spheroids settling through turbulence can be accurately predicted by a simple Gaussian statistical model that takes into account particle inertia and provides a quantitative understanding of the orientation distribution on the problem parameters when fluid inertia is negligible. Our results open the way to a parameterisation of the distribution of ice-crystals in clouds, and potentially leads to an improved understanding of radiation reflection, or particle aggregation through collisions in clouds. How non-spherical objects settle in a turbulent environment is a highly relevant question in several domains. An example is provided by very small ice crystals in clouds (size ∼ 100 µm), which grow through aggregation to form precipitation size particles (size ∼ 1 mm) [1][2][3][4]. The settling of plankton in the ocean [5][6][7] can induce patchiness of the population, therefore affecting mating, feeding and predation [8]. In these problems, the orientational degrees of freedom clearly affect not only settling and collision properties, but also light reflection [9]. As a prerequisite to a description of these effects, this Letter provides an understanding of the orientation statistics of small spheroids settling in a turbulent environment based on a statistical model, under the assumption that fluid inertia can be neglected.
The interaction between turbulence and settling leads to intriguing phenomena, even in the simpler case of spherical particles. Maxey found that turbulence increases the settling speed of a single small particle [10,11]. Substantial progress was recently achieved in understanding how two spherical particles settling together move relative to each other and collide [12][13][14][15][16].
In a fluid at rest the orientation dynamics of slowly settling non-spherical particles is determined by weak torques resulting from fluid inertia [17][18][19][20]. Turbulence affects the orientation of such particles through turbulent vorticity and strain. In the absence of settling this is well understood [7,[21][22][23][24][25][26][27][28][29]. Neglecting fluid inertia, the direct numerical simulations (DNS) of turbulence by Siewert et al. [30] demonstrated that settling induces a bias in the orientation distribution of the particles. The physical origin of this bias is not known, and it is not understood how the bias depends on the parameters of the problem: the turbulent Reynolds number, Re λ , the Stokes number (particle inertia), the gravitational acceleration, and the particle shape. Also, how significant are non-Gaussian, intermittent small-scale features of the turbulent flow [31], such as intense vortex tubes [2] in aligning the particles? To answer these questions we analyse a statistical model for the orientation of small heavy spheroids settling in homogeneous isotropic turbulence, for parameters relevant to cloud physics, and compare with results based on DNS of turbulence. Fig. 1 shows the predicted bias in the distribution of the vector n pointing along the particle symmetry axis. The statistical-model predictions agree very well with the DNS results. This shows that that non-Gaussian turbulent fluctuations are not important. The statistical model explains the sensitive parameter dependence of the DNS results. This is important because it allows us to parameterise the bias, to quantitatively understand the physical properties of the system. We analyse the model by an expansion in the 'Kubo number' Ku, a dimensionless correlation time of the flow [33]. Padé-Borel resummation yields excellent agreement with numerical simulations at Ku = 0.1, and qualitative agreement with DNS of turbulence. At larger Ku the theory fails to converge, but the model still explains qualitatively the underlying mechanisms. Last, we discuss possible effects of fluid inertia.
Formulation of the problem. The equations of motion for translation and rotation of a particle reads Here g is the gravitational acceleration (directionĝ), x is the position of the particle, n its symmetry vector, m its mass, ω its angular velocity, and J(n) is its inertia tensor in the lab frame. In the point-particle approximation, force f and torque T on a spheroid are [25,34,35]: In Eq. (2), v is the particle velocity, u(x, t) is the turbulent velocity field, Ω ≡ 1 2 ∇ ∧ u is half the turbulent vorticity, S is the strain-rate matrix, the symmetric part of the matrix A of fluid-velocity gradients (its antisymmetric part is called O), and M are translational and rotational resistance tensors: ⊥ )nn T , and M (r2) is a third-rank tensor. For a fore-aft symmetric particle, the equations of motion (1,2) are invariant under n → −n, so that only the magnitude n g ≡ |n ·ĝ| can play a role in the dynamics. The form of M (r2) and of the C and Kcoefficients are known for spheroidal particles, see Supplemental Material (SM) [32] and Ref. [36]. The parameter γ ≡ 9νρ f /(2a a ⊥ ρ p ) is Stokes constant, ν is the kinematic viscosity of the fluid, ρ f and ρ p are fluid and particle mass densities, 2a is the length of the particle symmetry axis, and 2a ⊥ is the particle diameter.
Our DNS of turbulence use the code described in [37] and in the SM [32]. The Kolmogorov scales u K , η K , and τ K are determined by the dissipation rate ε ≡ ν Tr AA T (the average is along steady-state Lagrangian trajectories), and by ν ≈1 × 10 −5 m 2 s −1 (air). The particle as-pect ratio is λ ≡ a /a ⊥ . The simulations were done for spheroids of varying λ and with max(a , a ⊥ ) =150 µm, much smaller than η K for values of ε pertaining to mixedphase clouds (DNS: ε ≈ 1, 16, and 256 cm 2 s −3 ). Particle inertia is measured by the Stokes number St K ≡ (γτ K ) −1 . The mass-density ratio is ρ p /ρ f ≈ 1000 (ice crystals in air), and the dimensionless gravity parameter is defined as Statistical model. The model is appropriate for particles smaller than η K . We approximate the universal [31] dissipative-range turbulent fluctuations by an incompressible, homogeneous, isotropic Gaussian random velocity field u(x, t) with zero mean, correlation length , correlation time τ , and rms speed u 0 [33] (details given in the SM [32]). In the persistent limit [33], for Ku ≡ u 0 τ / > 1, the model parameters St ≡ (γτ ) −1 and F ≡ gτ /u 0 map to St K = √ 5 Ku St and F K = [F/(5 Ku)] /η K . Here /η K is the ratio between the size of the dissipation range and the Kolmogorov length. In turbulence this ratio depends weakly on the Reynolds number Re λ [38], = cη K Re 1/2 λ . For the data shown in Fig. 1 we have Re λ = 95, and Fig. S1 in SM [32] shows results for other values of Re λ . We find good agreement between the statistical-model results at large Ku and the DNS for c ≈ 1.3. For Ku > 1, the model predictions depend on two parameter combinations only [33], Ku St and F/ Ku. In terms of the DNS parameters this means that the orientation bias depends only on St K and F K Re Perturbation theory. Eqs. (1,2) are solved by expansion in powers of Ku [33,39]. We outline the essential steps below, details are given in the SM [32]. We use dimensionless variables: t ≡ t/τ, r ≡ r/ , u ≡ u/u 0 , and drop the primes. To calculate the steady-state distribution of n g ≡ |n ·ĝ| we must evaluate the fluctuations of the fluid-velocity gradients along particle paths. This is achieved by an expansion in δx t ≡ x t − x ⊥ /(a a ⊥ (λ 2 + 1)). Eq. (4) is the lowest-order solution of Eqs. (1,2). The terms in Eq. (3) that do not involve δx t depend only on the history of the fluid-velocity gradients along the paths x (d) t ('history contribution'). The c-coefficients contain at most five powers of n 0 , and one must sum over all tensor products allowed by symmetry (Einstein convention). See SM [32].
The first integral shown in Eq. (3), by contrast, depends on δx t . It is therefore sensitive to how turbulence modifies the settling paths ('preferential sampling' [33]).
We determine the steady-state moments (n t ·ĝ) p ∞ by first calculating the moments conditional on the initial orientation n 0 , using Eq. (3) and the relation (n t ·ĝ) p n0 = (n 0 ·ĝ) p +p Ku(n 0 ·ĝ) p−1 n where n (i) t is the coefficient of Ku i in Eq. (3). Eq. (5) is valid to order Ku 2 . We average over the fluid-velocity fluctuations as described in Ref. [33]. The moments are independent of the initial position x 0 due to homogeneity of the flow. We expect that effects of the initial velocity v 0 and angular velocity ω 0 decay exponentially, so that they do not affect the steady state. We therefore set both to zero. Only the n 0 -dependence matters. In this way we obtain expressions for (n t ·ĝ) p n0 , which involve secular terms that increase linearly with time as t → ∞. But these terms must vanish since n t is a unit vector. This condition yields a recursion relation for the steady-state averages (n ·ĝ) p ∞ , independent of n 0 . This recursion is valid for arbitrary values of G ≡ Ku F St /C (t) ⊥ , and to order Ku 0 . Note that G can be large even if Ku is small. We solve the recursion by a series expansion in small G: . (6) The coefficients A (2i) j (St, λ) depend on the shape and inertia of the particle, but not on G or p. From Eq. (6) we obtain the Fourier transform of the probability distribution of n g = |n ·ĝ|. Inverse Fourier transformation yields the distribution. To order G 4 we find: The lowest-order term corresponds to a uniform distribution of n t . Let us examine the G 2 -term. It turns out that A 1 is negative for disks and positive for rods (see Fig. S2 in the SM [32]). This explains that the orientation of settling disks is biased: disks tend to fall edge on and rods settle tip first (as in Fig. 1).
Padé-Borel resummation. Now consider higher orders in the G-expansion. The series (6) is asymptotically divergent and must be resummed. Fig. 2 demonstrates that Padé-Borel resummation [33,40] of the series yields excellent results. Shown are results from a resummation of (6) to order G 34 (thick solid lines). These results agree very well with numerical simulations of the statistical model for Ku = 0.1 and St = 10 (symbols). The resummed theory works up to G = 10, and in this range the bias increases with increasing G. The resummed theory also predicts that the moments increase as St increases, for fixed G. A more detailed analysis of the recursion leading to Eq. (6) reveals, however, that the limit G → ∞ is delicate. Perfect alignment requires λ = ∞ [32].
In summary, perturbation theory in Ku shows that turbulence gives rise to an orientation bias ( Fig. 2), in excellent agreement with statistical-model simulations at Ku = 0.1 and in qualitative agreement with DNS ( Fig. 1).
The calculations leading to Eq. (6) reveal that each moment (n t ·ĝ) 2p ∞ is a sum of two contributions that stem from the 'preferential sampling' and 'history' terms in Eq. (3). For small Ku the history effect is dominant, the orientation bias is entirely determined by the history of fluid-velocity gradients along straight deterministic paths, Fig. 2a. Decomposing the leading-order contribution as A 1,hist. we find that |A 1,hist. | (Fig. S2b in the SM [32]). Fig. 3a leads to the same conclusion. It shows the distribution P (n g ) for Ku = 0.1. Also shown is P (n g ) computed for particles falling with constant velocity v = v s (n 0 ). We choose the squared initial orientation n 0 n T 0 in (4) as the steady state average n 0 n T 0 ∞ , evaluated using the small-Ku theory. This corresponds to keeping just the history contribution to n 2p g ∞ . We observe excellent agreement with the full statistical-model simulations. This shows that the history effect causes the orientation bias at small values of Ku.
Persistent limit. In the persistent limit we use numerical simulations with Ku = 10 to analyse the orientation bias in the same way as for small Ku. The result is shown in Fig. 3b (parameters correspond to two curves in Fig. 1(a). We plot the full statistical-model distribution and results for particles with a constant velocity (4) that neglects preferential sampling. For the data in Fig. 3b, the average n 0 n T 0 ∞ is computed using statistical-model simulations. We see that the history effect makes a substantial contribution to P (n g ). But since the distributions do not match, we infer that preferential sampling also contributes. This contribution is hatched in Fig. 3b. Limit of large settling speeds. Fig. 3c shows the moments n 2p g ∞ for p = 1, 2, 3 in the persistent limit as functions of the DNS Stokes number St K . Open symbols denote full statistical-model simulations, solid lines correspond to simulations based on straight deterministic paths. At intermediate Stokes numbers we see a clear difference between the two simulations, preferential sampling is important in this region.
As the Stokes number grows, however, the Figure demonstrates that preferential sampling ceases to play a role. In this limit the orientation bias is entirely caused by the history effect. The bias shown in Fig. 3c increases as St K increases. But as the perturbation theory indicates, the limit of large G is quite subtle. Statisticalmodel simulations for Ku = 1 show that the degree of alignment starts to decrease for very large G.
Conclusions. We analysed a statistical model for the orientational dynamics of small heavy spheroids settling in turbulence. The predictions of the model agree well with our own numerical results based on DNS of homogeneous isotropic turbulence (Fig. 1). Our statisticalmodel analysis shows that there are two distinct competing mechanisms causing the orientation bias: preferential sampling and the history effect. The latter dominates for large settling speeds, but it makes substantial contributions also in other parameter regimes. Preferential sampling dominates only when the bias is negligibly small. When the bias is significant, the history effect explains at least about 50% of the bias observed in Fig. 1.
We have shown that the orientation alignment depends on combinations of dimensionless numbers: St K and F K Re −1/2 λ . Our analysis shows that it is the smallscale properties of the flow that determine the orientation alignment. The Re λ -dependence arises only because it determines the ratio between the smooth scale to η K . We note that F K Re −1/2 λ equals the ratio of the settling velocity and the rms turbulent velocity fluctuations.
Our results pertain to small ice crystals settling in turbulent clouds, and allow us to model the sensitive depen-dence of the effect upon particle shape, size, and the turbulence intensity. This is important since turbulent dissipation rates vary widely in clouds. Our results predict strongly varying degrees of alignment. That the statistical model is in excellent agreement with the DNS opens a way to parameterise the orientation distribution of icecrystals in clouds. This potentially leads to an improved understanding of the radiative properties of clouds, and of particle aggregation through collisions in clouds.
The present work is based on the point-particle approximation of heavy particles, which neglects the effect of fluid inertia. This requires the particle Reynolds number Re p ≡ av c /ν to be small, where a = max(a , a ⊥ ). Estimating the slip velocity v c by the Stokes settling speed, we find that Re p is of order unity for the data shown in Fig. 1, so the condition is marginally satisfied. The shear Reynolds number, Re s , must also be small. Since Re s ≡ a 2 tr S 2 /ν ∼ (a/η K ) 2 [41], this condition is satisfied for small particles.
Lopez et al. [42] analysed the orientational dynamics of rods settling in a vortical flow. For small Re p they found a bi-modal distribution, with peaks at n g = 0 and 1. They explain the peak at n g = 0 by the effect of fluid inertia. Our results may explain the peak at n g = 1. These results, although not for a turbulent flow, indicate that turbulent and fluid-inertia torques compete in general. How to model this competition is an open question. For small Stokes numbers one may formulate an adhoc model by simply adding turbulent and fluid-inertia torques, along the lines suggested in Ref. [42]. But in general it remains a challenge to take into account effects due to fluid inertia from first principles, in a turbulent environment. Simulations resolving particle and fluid motion [43,44] and experiments [45][46][47][48] for micronsized particles in turbulence are needed to test the predictions, and to determine the orientational dynamics of larger particles where fluid inertia must matter [48]. Finally, how to extend the ideas developed here to particles lighter than the fluid remains a challenging task.

I. DETAILS ON EQUATIONS OF MOTION
We begin by recalling the form of the resistance tensors, defined by Eq. (2) in the main text. The tensor M (t) relates the force on the particle to the slip velocity, the difference between the fluid and the particle velocities. The superscript refers to the translational (t) degrees of freedom. The tensor is expressed as: For the rotational degrees of freedom, superscript (r), the resistance tensor is decomposed into two terms, M (r1) and M (r2) . The first term relates the torque to the angular slip velocity, the difference between the rotation rate of the fluid and that of the particle. The second term relates the torque to the local strain in the fluid. The expression for M (r1) is similar to Eq. (S1): The expression for M (r2) is of the form: Here ijl is the completely anti-symmetric Levi-Civita tensor, and the summation convention is used: repeated indices are summed from 1 to 3. For spheroidal particles, the coefficients in Eqs. (S1), (S2), and (S3) are known [1]: Here λ≡ a /as ⊥ is the aspect ratio of the spheroid, 2a is the length of the particle symmetry axis, and 2a ⊥ is the particle diameter. Prolate spheroids have λ > 1. Oblate spheroids correspond to λ < 1. In this case, Eq. (S4c) reduces to β = arccos(λ)/[λ √ 1 − λ 2 ]. The motion of the particle is fully characterised by the knowledge of v, the velocity of the particle, and the angular velocity of the particle in the laboratory frame, ω. The direction of the particle symmetry vector, n, is related to the angular velocity ω by the kinematic equation: The equation for the angular velocity reads in the laboratory frame: Here T is the torque on the particle, and J(n) is the inertia tensor of the particle. For a spheroid its elements are given in Ref. [1]. A difficulty is that the inertia tensor depends on the instantaneous orientation n(t). In the DNS we therefore express Eq. (S6) in the particle frame, at the cost of an extra term in the equations of motion, due to the rotation of the axes [2].

A. Method of resolution
The Navier-Stokes equations read: where u is the velocity field of the flow, which is assumed to be incompressible: p is the pressure field, Φ is a forcing term, ρ f is the density of the fluid, and ν is the kinematic viscosity. The equations are solved in a triply periodic domain of size L 3 , using a pseudo-spectral method. The code has been described elsewhere [3]. Briefly, the fluid is forced in the band of lowest wave-numbers by imposing an energy injection of ε. The value of ε chosen here correspond to cloud conditions, see subsection II B, in such a way that he Kolmogorov length, η K = (ν 3 /ε) 1/4 is of the order of the grid spacing. The code is fully dealiased, using the 2/3-rule. Specifically, if N is the number of grid points in each direction, the nonlinear term is computed by using only Fourier modes 0 ≤ n ≤ N/3. The equation of motion of the particles is implemented by interpolating the velocity and velocity gradients at the location of the particle, using tri-cubic schemes. The method has been carefully checked, by systematically comparing the results in the absence of motion u = 0, with the code, and with an elementary solver (Mathematica) of the equations of motion of the particle.

B. Choice of parameters
The values chosen in this study correspond to cloud conditions [4]. In physical units, the viscosity chosen here is ν = 0.113cm 2 /s. The number of points taken here was N = 384, with an energy injection rate of ε ≈ 1 cm 2 /s 3 , N = 784 (ε ≈ 16 cm 2 /s 3 ), and N = 1568 (ε ≈ 256 cm 2 /s 3 ). With the values of the parameters chosen here, the size L of the box is equal to 8π cm. ). The DNS data is compared to statistical-model simulations with large Ku. The statistical-model parameters St, and F are determined by the relations to the DNS parameters Re λ , St K , and F K described in the main text. The agreement is equally good as that observed in Fig. 1 in the main text.

IV. DETAILS ON THE STATISTICAL MODEL
In the simulations of the statistical model the smooth, incompressible, homogeneous, isotropic, threedimensional velocity field is given by u = ∇ ∧ Ψ/ √ 6, where Ψ is a vector potential [5]. Each component of Ψ consists of a spatial superposition of Fourier modes with random time-dependent prefactors. The statistics of the prefactors is chosen such that Ψ i (x, t) is Gaussian distributed with correlation function where . . . denotes an ensemble average. We choose the velocity scale as the rms speed of the flow, u 0 ≡ |u| 2 . The correlation length and correlation time τ are the scales at which Eq. (S9) decay in space and in time. These Eulerian scales constitute the Kubo number Ku = u 0 τ /η [5].

A. Implicit solution of equations of motion
We follow the method outlined in Ref. [5] and expand the solution of the equations of motion for the settling spheroid, Eqs. (1) and (2) in the main article, in powers of Ku: (S10) Here Λ ≡ (λ 2 − 1)/(λ 2 + 1), and we have decomposed the translational resistance tensor into a part that does not depend on the Kubo number, and a Ku-dependent part, Here ∆n is a rescaled orientation displacement vector The rotational resistance tensors are decomposed in a similar way. Eq. (S10) is an exact, but implicit, solution to Eqs. (1) and (2) in the main article, provided that the initial conditions of the particle velocity and angular velocity are set to zero. The solution (S10) depends implicitly on orientation n, the angular velocity ω, and on the centre-of-mass position x of the particle. If the Kubo number is small we can terminate the solution (S10) at some order in Ku and solve the resulting set of equations iteratively. Details are given in Ref. [5].

B. Explicit solution of equations of motion for small values of Ku
To eliminate the implicit dependence on x through the flow velocity u(x(t), t) and through the flow velocity gradients A(x(t), t), we consider a decomposition of particle trajectories x(t) into two parts: (S13) The first part, x (d) (t), is a deterministic part that is obtained by setting the turbulent velocity u to zero in the particle equation of motion. The second part is the remainder, the flow-dependent fluctuating part, δx(t) ≡ x(t) − x (d) (t). The deterministic part of the trajectory is obtained by integrating the solution for the particle velocity in Eq. (S10) with u = 0 and A = 0. For large values of t, the corresponding particle velocity approaches the settling velocity v s (n 0 ) in a quiescent fluid The corresponding deterministic trajectory approaches x (d) (t) = x 0 + v s (n 0 )t. Here n 0 is the constant orientation at which the particle settles. In the steady state, the orientation n 0 must be chosen from a stationary distribution of particle orientations that must be determined self consistently.
For small values of Ku, the flow gives rise to small deviations δx from the deterministic settling trajectories x (d) . Following the procedure outlined in Ref. [5], we attempt a series expansion of the fluid velocity experienced by particles in terms of small δx and similarly for the fluid-velocity gradients. We insert these expansions into the implicit solutions (S10), and iterate the solution to find an expression for the orientation n valid for small values of Ku. To second order in Ku we find Eq. (3) in the main article. Writing out all terms explicitly this equation reads in vector notation: t , t). The coefficents are: where Λ = (λ 2 − 1)/(λ 2 + 1). To obtain the first term in the last row of Eq. (3) in the main text, for example, we must add the terms involving d 1 , d 2 , and d 10 .

C. Moments of alignment
Eq. (S16) relates the orientation dynamics to the fluid velocity and its gradients. This expression makes it possible to calculate the moments (n ·ĝ) 2p n0 conditional on n 0 (odd-order moments vanish because spheroids are fore-aft symmetric). By raising Eq. (S16) to the power 2p, terminating the resulting expressions at order Ku 2 , and taking the average using the known stationary statistics of the fluid velocity and its gradients, we obtain the moments of n ·ĝ conditional on initial alignment n 0 ·ĝ.
However, the resulting expressions contain secular terms. These terms are proportional to Ku 2 and grow linearly with time. We know that such secular terms must vanish, because the vector n must remain normalised to unity. Its components cannot continue to grow. We therefore demand that the initial alignment n 0 ·ĝ is distributed according to a stationary distribution of alignments that causes all secular terms to vanish. This procedure gives rise to one condition per value of p, and this is sufficient to determine the desired stationary distribution of alignments.
The dependence of the secular terms upon (n 0 ·ĝ) 2p is not only through algebraic powers, but also functional through exponentials and error functions. As a consequence, the conditions that make sure that the secular terms vanish are hard to solve in general. We therefore seek a perturbative solution and expand the secular conditions and the moments m p ≡ (n 0 ·ĝ) 2p in terms of the parameter G ≡ Ku F/C (t) ⊥ , for instance: To lowest order in G we find the recursion Using the boundary condition m (0) 0 = 1 we find These moments correspond to a uniform distribution of n 0 . Higher-order contributions in G to the moments can be recursively determined by series expansion of the secular terms. We find that moments m . (S20) The coefficients A (2i) j (St, λ) are quite lengthy in general. We therefore only give the lowest-order expression here. To order G 2 we have .
Decomposing A 1 (St, λ) in one contribution due to preferential sampling and one due to the history effect, we have The preferential-sampling contribution reads: The contribution due to the history effect is 30Λ 2 − 49Λ + 49 + 4C and C (r) = 10 3 to simplify. Eq. (S4) relates C (t) ⊥ and C (t) to λ, so that A 1 is a function of St and λ only. Fig. S2 shows how the expressions for A (2) 1 and for the relative preferential-sampling contribution |A (2) 1,pref. /A (2) 1 | depend on the aspect ratio λ, for different values of St. We see that the preferential contribution is small compared to the contribution due to the history effect for all parameter values.

D. Large-G limit
In the previous Section we summarised how the moments of the orientation distribution can be calculated by using perturbation expansions in the parameter G. Now consider the opposite limit of large values of G. In this limit we obtain a recursion equation for the moments (n ·ĝ) 2p that is independent of preferential effects. This is expected because the particles fall rapidly through the turbulence in this limit, experiencing the turbulent gradients approximately as a white-noise signal. The history contribution is obtaind from a fourth-order recursion that is hard to solve in general.
One might expect that the orientation of rod-like particles approaches perfect alignment with gravity, that is (n ·ĝ) 2p = 1 for all values of p as G → ∞. However, when we insert this limiting expression into the fourth-order recursion equation, we obtain a condition that is only satisfied when λ = ∞. This shows that perfect alignment is only possible when λ = ∞. For other values of λ the moments must approach limiting values as G → ∞. We have not yet managed to calculate these limits from our statistical-model theory.