Passive directors in turbulence

L. Zhao,1,2,* K. Gustavsson,3,* R. Ni,4 S. Kramel,5 G. A. Voth,5 H. I. Andersson,2 and B. Mehlig3 1Department of Engineering Mechanics, Tsinghua University, 100084 Beijing, China 2Department of Energy and Process Engineering, NTNU, 7491 Trondheim, Norway 3Department of Physics, Gothenburg University, 41296 Gothenburg, Sweden 4Department of Mechanical Engineering, Johns Hopkins University, Baltimore, Maryland 21218, USA 5Department of Physics, Wesleyan University, Middletown, Connecticut 06459, USA


I. INTRODUCTION
Suspensions of small particles in turbulence determine the physics and chemistry of many natural processes. The analysis of the underlying highly nonlinear and multiscale dynamics poses formidable challenges, because any description of the problem must refer to the turbulence that the particles experience as they move through the fluid. Experiments resolving the particle dynamics have only recently become possible, and direct numerical simulations (DNSs) of such systems are still immensely difficult. Recently there has been substantial progress in understanding the dynamics of spherical particles in turbulence by means of statistical models [1,2].
Yet most solid particles we encounter in nature and engineering are not spherical, such as ice crystals in turbulent clouds [3], plankton in the turbulent ocean [4][5][6][7], and turbulent fiber flows in industrial processing [8]. Therefore, it is necessary to understand how nonspherical particles translate and rotate in turbulence [9]. For very small particles, inertial effects are negligible [10], and the disturbance caused by the particles can be treated in the Stokes approximation [9,11,12]. Understanding the turbulent angular dynamics of nonspherical particles in this limit is a question of great current interest [9,[13][14][15][16][17][18][19]. However, even the angular dynamics of a single small rod in turbulence is quite intricate: Rods tend to align with the local vorticity of the flow [14,18,20]. Vorticity in turn aligns with the second eigenvector of the turbulent strain-rate matrix [21] and picks up that turbulence breaks time-reversal invariance [22,23]. Polymers tend to align with the main stretching direction of the flow [24].
Very little is known about how nonspherical particles orient relative to each other in a turbulent flow, even in the limit of very small particles whose center-of-mass position x simply follows the flow (passive particles). How do nearby nonspherical particles align with each other? The simplest case is that of axisymmetric particles with fore-aft symmetry. The symmetry axes n of an ensemble of such particles in turbulence form a spatial field n(x, t ). At every point in space and time n is normalized to unity. Our goal is to determine the geometrical properties of this field. For foreaft symmetric particles the problem is invariant under n → −n, so n(x, t ) is in effect a field of directors. It is plausible that turbulent strains align the particles as they approach, and in this case one expects the passive director field n(x, t ) to be a smooth function of particle position. However, in this paper we show that the spatial passive-director patterns of nonspherical particles in turbulence are not smooth in general, not even at the smallest scales where the turbulent fluid velocities are smooth functions of position. We show that angles between the symmetry axes of nearby particles are anomalously large. Our results of DNSs and statistical-model simulations allow us to conclude, first, that the attractor determining the director patterns is fractal, in general, and second, the steadystate distribution of angles between nearby particles has power-law tails. We derive a theory based on diffusion approximations that can at least qualitatively explain these observations and we relate the power-law tails to steps of different widths that occur in the director patterns.

II. FORMULATION OF THE PROBLEM
The center-of-mass positions x of small particles simply follow the flow u(x, t ), if their spatial diffusion is neglected. In this case one says that the center of mass is advected by the flow [1]. The directors of small axisymmetric particles with fore-aft symmetry obey Jeffery's equation [11] d dt n = On + Sn − (n · Sn)n.
Here O is the antisymmetric part of the matrix A of fluid-velocity gradients at the particle position and S is its symmetric part. Inertial effects [10] and angular diffusion [25] are neglected in Eq. (2). The parameter parametrizes particle shape [11,26]: = 0 for spheres, =−1 for thin disks, and = 1 for slender rods. For spheroids where κ is the aspect ratio of the spheroidal particle [11]. In the following we consider prolate particles with 0. We use two different kinds of simulations to analyze the director patterns n(x, t ). First, we employ DNS of a turbulent channel flow to obtain a turbulent velocity field and integrate Eqs. (1) and (2) numerically to determine the director patterns. Second, we perform simulations of Eqs. (1) and (2) using a statistical model, representing the small-scale turbulent velocities as a random Gaussian incompressible, homogeneous, and isotropic field with correlation length η, correlation time τ , and rms speed u 0 . The theory employs diffusion approximations that are valid in the limit of small Kubo number Ku ≡ u 0 τ /η and that have yielded important insights into the dynamics of small spherical particles in turbulence [2], just like Kraichnan's diffusion model for passive-scalar advection [1]. Our experiments measure the angles between symmetry axes of fibers advected in a turbulent flow between oscillating grids.
To characterize the director patterns n(x, t ), we measure the statistics of δn(x, t ) ≡ n(x + R, t ) ± n(x, t ) at small distances R ≡ |R|. The configurations ±n(x, t ) correspond to identical physical situations, and we choose the sign so that |δn| is minimal. We define the angular structure functions 054602-2 Here · · · R is a steady-state average over particle pairs conditional on their center-of-mass distance R. The order p of these moments need not be an integer, and it can also assume negative values. We consider small separations R between the particles, in the dissipation range of turbulence. In this range the second-order longitudinal velocity structure function δu 2 L R of the turbulent flow scales as [27][28][29][30][31][32] Here ·R is the velocity increment of the turbulent velocity in the directionR of separation R. This means that the turbulent velocity field is smooth in the dissipation range. It can be Taylor expanded to give δu L ∝ R. If the director field n(x, t ) were smooth too, then S p (R) ∝ R p for small R, as in Eq. (5). However, our results show that this is usually not the case. In our DNSs and statistical-model simulations we find instead anomalous scaling of the director field with |ξ p | |p|. Since |ξ p | |p| this implies that angles between the symmetry axes of nearby particles are anomalously enhanced at small separations. This means that the spatial director field of nonspherical particles is not smooth, in general, so angles between the symmetry axes of nearby nonspherical particles in turbulence are larger than expected. This is our main result. The remainder of this paper is organized as follows. In Sec. III we give details about our experiments, DNS, and the statistical model that we use to analyze the director patterns. Section IV summarizes the results of our experiments, DNS, and statistical-model analysis. Section V contains a discussion, conclusions, and an outlook.
In the Re λ = 140 flow the rod length is 2.3η K , which is small enough that the particles are in the tracer limit. In the Re λ = 277 flow, the rod length is 5.2η K . For particles of this size, tumbling rates are still roughly in the tracer limit [19], but finite-size effects start to become important. The Stokes numbers of the rods, defined as the ratio of the rod response time to the Kolmogorov time, are 0.002 and 0.01, so the particles behave like neutrally buoyant tracers even though the fluid density ρ f = 1.00 g cm −3 is slightly lower than the particle density ρ p = 1.12 g cm −3 .
The two experiments used different imaging setups. The lower-Reynolds-number data were taken with a laser scanning system and three cameras recording 5000 frames per second (fps) [32]. This data set has 3.8 × 10 4 frames. Each frame imaged one of eight slabs that were scanned sequentially. A typical image contains 60 rods within a slab with dimensions 1.5 × 1.5 × 0.3 cm 3 . The higher-Reynolds-number data were taken with the imaging system using volume illumination and four cameras at 450 fps [19]. This data set has 1.5 × 10 6 frames, each with typically eight rods in view in an imaging volume with dimensions 2 × 2 × 3 cm 3 . The seeding densities chosen are a compromise between obtaining sufficient numbers of rod-angle differences at small distances and minimizing the overlap of rods in the two-dimensional (2D) images. When rods overlap, it is difficult to separate 054602-3 them and measure their 3D positions. We discard such samples. Because rods with large-angle differences are more likely to overlap in the images, this introduces a sampling bias at small R. We therefore only report data for R greater than the rod length, where the bias is not large.
The experimental raw data are analyzed as follows. The camera images [ Fig. 1(a)] are first segmented, identifying clusters of bright pixels. The two-dimensional centers of mass of clusters are then stereo matched and tracked over time using the camera calibration data and a predictive tracking algorithm [33]. Rod angles are extracted from multiple images using methods described previously [15,34]. Figure 1(b) shows the three-dimensional reconstruction of the locations and angles between a pair of nearby rods [same as the pair shown in Fig. 1(a)]. The rods used in these experiments have an aspect ratio of more than four times that of previous experiments on tracer rods [15], resulting in a smaller uncertainty of the relative angle, compared to previous measurements.

B. Direct numerical simulations
We perform DNS of a turbulent channel flow using one-way coupling for the particle dynamics: Given the fluid-velocity field and its gradient, spheroids with shape parameter move according to Eqs. (1) and (2). The turbulent channel flow is characterized by the Reynolds number Re * = u * h/ν based on the wall-friction velocity u * and the half-channel height h. The wall-friction velocity is determined by the wall stress and the fluid density. Since the channel flow is inhomogeneous, Re λ varies throughout the channel cross section. We take our statistics near the channel center, in a region of linear size 2η K , where the turbulent vorticity is approximately homogeneous and isotropic [35]. In this region we estimate Re λ = u rms λ/ν using the local rms turbulent velocity u rms = |u | 2 /3 1/2 and the Taylor scale λ = u rms √ 15ν/ . The prime denotes the fluctuating part of the fluid velocity obtained by Reynolds decomposition. The dissipation rate is calculated from the local turbulent velocity gradients = ν TrA T A . We choose Re * = 180. Near the channel center this gives u + rms = 0.686, + = 5.4 × 10 −3 , λ + = 36.2, Re λ = 24.8, η + K = 3.68, τ + K = 13.6, and u + K = 0.27. All nondimensional quantities are quoted in wall units, expressed in terms of u * and ν. The simulation domain is 12h × 6h × 2h in the streamwise, spanwise, and wall-normal directions. We apply periodic boundary conditions in the spanwise and streamwise directions and no-slip boundary conditions at the two walls. We use a pseudospectral method in the periodic directions and a second-order central finite-difference scheme [36] in the wall-normal direction. For time integration we use an explicit second-order Adams-Bashforth scheme.

C. Statistical model
In two spatial dimensions we use a stream function (x, t ) to represent a smooth, incompressible, homogeneous, isotropic random Gaussian velocity field u: We take [2] The stream function is constructed as a superposition of Fourier modes with Gaussian random timedependent coefficients. The coefficients are chosen such that (x, t ) has zero mean and correlations This correlation function defines the Eulerian scales of the flow, namely, the correlation length η and the correlation time τ . The typical speed is u 0 ≡ u 2 .
In three spatial dimensions, the velocity field u(x, t ) is constructed as the rotation of a vector field A(x, t ) [2]: The three components of A(x, t ) are independent Gaussian random functions with the same statistics as (x, t ). Further details of the statistical model are described in Ref. [2].

IV. RESULTS
This section is organized into two parts: a summary of the experimental and DNS data (Sec. IV A) and a summary of the statistical-model results (Sec. IV B).
In Sec. IV A we show our results for the angular structure functions S p (R) that illustrate our key conclusion, that the probability of observing large angles between nearby particles is anomalously large. Then we describe our analysis of the anomalous exponents ξ p in Eq. (6): for DNSs and statistical-model simulations in two and three spatial dimensions.
The statistical-model analysis (Sec. IV B) yields three main results. First, the exponents ξ p saturate for large p because the distribution of relative angles has power-law tails (Sec. IV B 1). Second, these tails result from a steady-state distribution of steps of different widths in the director patterns that occur when turbulent strains act on particles aligned orthogonal to the main stretching direction (Sec. IV B 2). We discuss how these steps are related to scar-line singularities observed in 2D random [37], deterministic chaotic [38], and quasiperiodic [39] maps. Third, we show that the attractor determining the steady-state director patterns is fractal for small and thus not smooth (Sec. IV B 3). There is a phase transition: For < c the director patterns are fractal, but for > c they become locally smooth. In Sec. IV B 4 we demonstrate that patterns in three spatial dimensions are similar to the patterns analyzed for the two-dimensional statistical-model system.

Angular structure functions
Figures 2(a)-2(c) shows our experimental results for the angular structure functions S p (R) (symbols). In our experimental apparatus the dissipation range extends up to roughly 10η K , where η K is the Kolmogorov length. More specifically, Fig. 2(d) shows that δu 2 L R ∝ R 2 in this range. The experimental data in Fig. 2(d) are taken from Fig. 4 in Ref. [31], obtained using the same measurement apparatus as in this paper, but at a slightly larger Reynolds number Re λ = 285. At larger separations an inertial-range power law [27][28][29][30][31][32] emerges (the data are roughly consistent with δu 2 L R ∝ R 2/3 ). The DNS data in Fig. 2(d) exhibit dissipation-range scaling δu 2 L R ∝ R 2 up to approximately 10η K . There is no inertial-range scaling for the DNS data because the Reynolds number is quite small, Re λ = 24.8.  Figures 2(a)-2(c) demonstrate that the angular structure functions decay much more slowly than R p in the dissipation range. Also shown are the results of our DNS (solid lines). We see that experimental and DNS results agree in the range where we have both experimental and DNS data. In the experiment the range of spatial scales in the dissipative range is too small to extract reliable values for the scaling exponents, but the DNS results exhibit clear power-law scaling with anomalous exponents ξ p p for p = 1 2 , 1, and 2.

Figures 3(a)-3(c)
show DNS results for the exponents ξ p , as well as results of 3D statisticalmodel simulations. The DNS results shown in Fig. 3(a) are obtained by fitting the DNS data for log S p (R) to A + ξ p log R. We consider two fitting ranges 0.05 R/η K 0.5 and 0.02 R/η K 0.2. The resulting estimates for ξ p are almost the same; the largest discrepancy is 5%. The data displayed in Fig. 3(a) correspond to 0.05 R/η K 0.5. The statistical-model results are obtained in a similar way. The values of ξ p are obtained by a linear least-squares fit to log S p (R) = A + ξ p log R. The data are fitted in two ranges 0.005η < R < 0.01η and 0.01η < R < 0.05η. In both ranges the values of C( ) and ξ ∞ converge to the same values, except for ∼ 1, where differences of order 0.1 are observed in ξ p . The data displayed are for the range 0.01η < R < 0.05η. Figure 3(a) exhibits good qualitative agreement between these simulations. For small values of |p|, the exponent ξ p is proportional to p, To estimate the constant of proportionality, we fit ξ p = C( )p to the DNS data in the range −0.1 p 0.1. The corresponding statistical-model results are obtained from the average value of ξ p /p in the same range, −0.1 < p < 0.1. Figure 3(b) shows the results, how C( ) depends on . The values observed in the DNSs and the statistical-model simulations are slightly different, but in both cases we find C( ) 1 for small values of . This observation is explained by the fact that the director patterns are fractal, as we show below, and this is one source of the large differences between angles of nearby rods. As increases, C( ) grows. At large values of we see that C( ) approaches unity; however, as long as C( ) < 1 the director field is fractal. For large values of p, by contrast, ξ p tends to a constant. Both the DNS and the statistical-model simulation show this behavior, characteristic of highly intermittent fields [30] ξ p → ξ ∞ ( ) as p → ∞.
We demonstrate below that this saturation is caused by narrow steps in the director field across which rods rotate by π . This is a second source of large-angle differences at small R. Figure 3(c) shows how ξ ∞ ( ) depends on . The plateau values ξ ∞ ( ) are obtained as averages of ξ p in the interval 4 < p < 5. For large , ξ ∞ ( ) is very close to unity. In summary, the scaling exponents are well approximated by ξ p ≈ min{p, 1} in the limit of large .
Note that the range of in Figs. 3(b) and 3(c) exceeds the physical limit for slender rods, = 1. Equation (3) shows that values of > 1 correspond to spheroids with imaginary aspect ratios κ (however, one can construct particles that have > 1 [26]). It is nevertheless instructive to consider the limit of large : We show in Sec. IV B 1 that the problem admits an exact solution for ξ p in the limit 1 and Ku 1. There we also discuss which insights this solution gives and how it fails at small values of .
B. Statistical-model analysis The results are very similar to the 3D case. We therefore begin by analyzing the 2D model. The problem is simplest to analyze in the white-noise limit Ku → 0, although turbulence corresponds to large Kubo numbers [2]. We will argue in this section that the white-noise limit nevertheless yields important insights, just like white-noise approximations for heavy-particle dynamics in turbulence [2] or Kraichnan's model for passive-scalar advection [1].

Large-limit
The reason why the white-noise limit is simpler to analyze is that diffusion approximations can be applied. To make progress in our problem we must currently assume that 1 and Ku 1. This means essentially that strain dominates over vorticity in aligning the particles. In this limit 054602-7 we can compute the steady-state form of the distribution P(R, δψ ) of center-of-mass distances R and angle differences δψ between particle pairs. Details are given in Appendix A. The result is of power-law form The factor N is a normalization constant. Figure 4(a) shows the results of statistical-model simulations at Ku = 0.01 and = 21 for the joint distribution P(R, δψ ) of angle differences and separations. We see that Eq. (12) is an excellent approximation at small Kubo numbers and large values of . Evaluating Eq. (4) with the distribution (12) gives for small R and p = 1, with coefficients a p = 2 −p 3 p/2 cos(pπ/2) and b p = 2 −p+1 √ 3π p−2 /(p − 1) (Appendix A). Comparison with Eq. (6) shows that for 1 and Ku 1. The saturation of ξ p for large p is a consequence of the power-law tail of P(R, δψ ). Equation (14) is consistent with the large-numerical results: Figs. 3(b), 3(c), 3(e) and 3(f) show that C ≈ 1 and ξ ∞ ≈ 1 down to = O(1). This indicates that the director patterns are smooth for 1 and Ku 1. This is no longer true for small . In Sec. IV B 3 we show that the director patterns become fractal at small values of .
The large-approximation discussed above fails to account for the fractality observed at small . However, numerical simulations demonstrate that the qualitative conclusions remain unchanged. Figure 4(b) shows that the distribution P(R, δψ ) still has power-law tails, albeit now with exponents different from −2 (the exponent is approximately equal to −1.5 for = 2 3 ). For small values of R, these power-law tails in δψ are cut off at δψ ∼ π/2, independent of R. As a consequence, the exponents ξ p saturate for p 1, but now at a constant smaller than unity.

Scar lines
What causes the power-law tails in P(R, δψ )? Consider a simple model in which strain is constant in space and time and vorticity is zero. We take Thus all initial angles converge toward the extensional strain direction ψ = ±π/2, except for ψ (x 0 , 0) = 0, marking the location of steps of height π in the director pattern. These steps are related to singularities observed in computer simulations of slender rods rotated by 2D random [37], deterministic chaotic [38], and quasiperiodic [39] maps. These singularities occur where the extensional eigenvector of S is orthogonal to the initial director pattern n(x 0 , 0) [37,40], just as the steps in the example above. In two dimensions, this constraint is satisfied on lines [37][38][39], termed scar lines in Ref. [37]. In turbulence A(x, t ) changes as a function of space and time, so new steps are continuously created, old steps sharpen, and their height approaches π . Since the problem is invariant under ψ → ψ + π , the steps effectively disappear as they sharpen. This is illustrated in Fig. 5(a) for = 1 where the attractor is smooth. Older steps leave thinner traces because they are less likely to be sampled by particles. We conclude that a steady-state distribution of steps of different widths develops, independently of the initial condition.
How are these steps related to the power-law tails in P(R, δψ ) and the saturation of the exponents ξ p at large p? For large , where the director patterns are smooth, we can estimate the width w s of a step in the x 1 direction as w s ∼ π/|∂ 1 ψ|, where ∂ 1 ψ is the derivative of ψ with respect to x 1 . In the limit of large and small Ku we find, using diffusion approximations (Appendix B), that P(∂ 1 ψ ) ∼ |∂ 1 ψ| −2 for large values of |∂ 1 ψ|. To obtain the step contributions to the angular structure function for large p, we note that a step of width w s contributes with weight w s /R to |δψ| p R for w s < 2R. Wide steps with w s > 2R give a smooth contribution of order (R/w s ) p−1 . Upon integrating over the distribution of w s up to w s 2R we find that |δψ| p R ∼ R for large p, establishing the connection to Eq. (13).

Fractal director patterns
The discussion in Secs. IV B 1 and IV B 2 applies only in the limit of large when the director patterns are smooth. Now we show that the director patterns are fractal at smaller . Consider first two dimensions. In this case the phase space of Eqs. (1) and (2) is three dimensional, spanned by the two components x 1 and x 2 of the center-of-mass position x and by the angle ψ of n with the x 1 axis, say. We find that the scatter of points in this phase space is not smooth, but fractal. To characterize the fractal we compute the Lyapunov dimension D L [2,41] to order Ku 2 (Appendix C). The result is . Equation (16) shows that the phase-space attractor is fractal for 0 < < c because D L is not an integer in this range. There is a phase transition at c and the fractal dimension D L equals 2 for > c , indicating that the attractor is smooth in this range. Figure 6 shows numerical results from statistical-model simulations for D L , in two spatial dimensions. We observe good agreement with Eq. (16) for small Ku. We also see that D L depends only very weakly on Ku [note that Eq. (16) does not apply for Ku = 1]. This indicates that preferential-sampling effects [2] are weak.
For a Gaussian random function f (x) with power-law spatial correlations, Orey [42] derived a relation between the increments δ f ≡ f (x + R) − f (x) and the fractal Hausdorff dimension D of the set of points embedded in the (d + 1)-dimensional space with coordinates x and f : To lowest order in Ku, Eq. (16) gives D L ∼ 3 − 2 2 . Setting D = D L in Eq. (17) yields C( ) = 2 2 , roughly consistent with the numerical results in Fig. 3(e). A more quantitative comparison would have to take into account that the Hausdorff and Lyapunov dimensions assume, in general, different values. Figure 5(b) illustrates steps in the director patterns for = 2/3 where the attractor is fractal. Steplike structures are still present, but less distinct because the attractor is fractal. In inviscid one-dimensional Burgers turbulence [43,44] the turbulent velocity structure functions |δu| p R obey |δu| p R ∼ R ζ p at small R and ζ p → 1 at large p. Bifractal theory [30] relates the saturation of the exponent ζ p at large p to steps (shocks) in the velocity field. This provides an insightful analogy: The saturation of the scaling exponents is caused by steps in the spatial field, as in our problem. However, there are fundamental differences: The velocity field in inviscid Burgers turbulence exhibits sharp jumps on a fractal set of positions, but the director field is itself fractal, in addition to exhibiting jumps. How to calculate the exponents in our system for 0 < < 1 is an open problem, even in the diffusive limit Ku → 0.

Three spatial dimensions
The argument leading to Eq. (15) generalizes to three dimensions. We can therefore conclude that steps form also in three dimensions. In this case phase space is five dimensional: three centerof-mass dimensions plus two Euler angles for the azimuthal (ϕ) and polar (θ ) degrees of freedom. Figure 5(c) shows a director pattern from 3D statistical-model simulations. The value of ϕ is color coded and plotted as a function of x 1 , x 2 , and x 3 . The spatial pattern is consistent with steps: In a given x 2 -x 3 plane we see sharp transition lines where ϕ jumps by π , and these lines appear also in neighboring x 2 -x 3 planes, at different values of x 1 .
For smaller we expect that the director patterns are fractal, as in two dimensions. The numerical results from 3D statistical-model simulations shown in Fig. 7(a) demonstrate that this is the case. The results indicate that there is a phase transition, as in two dimensions. Since phase space is five dimensional, D L changes from 5 at = 0 to 3 at large values of . For Ku = 1 the critical shape parameter is approximately c ≈ 1. Also shown are data for Ku = 0.1. The results for Ku = 1 and 0.1 are slightly different, unlike in two spatial dimensions. This could be due to numerical errors: The Ku = 0.1 data for D L in three spatial dimensions are the most difficult to obtain among the displayed statistical-model data. However, we cannot exclude that there is a Ku dependence in three spatial dimensions. This would indicate that preferential effects [2] or large time correlations matter.
For the turbulent channel flow our conclusions are qualitatively similar, although we could not determine the Lyapunov dimension reliably because the flow is inhomogeneous and the long trajectories needed to estimate D L do not remain in the center of the channel where we must take the statistics. Therefore, we have numerically computed the correlation dimension D 2 . It is defined as [45] P(|δw w w|) ∼ |δw w w| D 2 −1 as |δw w w| → 0, where δw w w = [δx 1 , δx 2 , δx 3 , δθ, δϕ] T and δθ and δϕ are differences of particle azimuthal angles θ and polar angles ϕ. The value of the power-law exponent in (18) is obtained by fitting Eq. (18) to the DNS data, in the range 0.02 < |δw w w| < 0.2. We also tested a slightly lower range, from 0.01 to 0.1. This makes only a small difference in the results, but to quantify the error it would be necessary to measure at substantially smaller values of |δw w w|. Our present DNS data do not permit this. 054602-11 Figure 7(b) shows that the fractal correlation dimension D 2 ranges from 5 at = 0 to approximately 3 at large values of , like the Lyapunov dimension for the 3D statistical model. However, for the DNS it is difficult to estimate c precisely. Figure 7(b) indicates that the critical shape parameter is of order unity, so the steady-state attractor is fractal for generic axisymmetric particles with fore-aft symmetry, except possibly for very slender rods. Figure 2 shows experimental and DNS data that illustrate our key result: Angles between the symmetry axes of nearby particles are anomalously large. We quantified this phenomenon by analyzing angular structure functions S p (R) at small spatial scales R. These functions measure the moments of relative angles between nearby particles. Direct numerical simulation and statisticalmodel simulation data demonstrate that the angular structure functions S p (R) exhibit power laws at small R, with anomalous scaling exponents ξ p . Using diffusion approximations in the limit 1 and Ku 1, we derived a theory that qualitatively explains the dependence of ξ p upon p. The exponents saturate at large p. This is caused by turbulent strains that form a steady-state distribution of steps in the director patterns. In the limit of 1 the director patterns are spatially smooth. However, for small , the director patterns are fractal, and this explains the small-p limit of the scaling exponents ξ p . At large p the exponents saturate, but no longer to unity. The saturation is still a consequence of steps in the director patterns, but now the steps inherit their properties from the fractal nature of the attractor.

V. DISCUSSION AND CONCLUSIONS
The experimental angular structure functions agree very well with the DNS results (Fig. 2), down to separations of the order of the rod length, although the DNS is performed at much smaller Reynolds number. There is excellent agreement between the DNS results and our statistical-model predictions. These facts indicate that the mechanisms causing anomalously large angles are robust, so our theory provides a foundational framework for understanding the relative alignment of particles in turbulence.
Our problem is related to, yet fundamentally different from, passive-scalar and vector problems [1,30,[46][47][48][49][50]. Instead of the magnitude of a scalar or a vector, we analyze the spatial director field n, normalized to unity and invariant under n → −n. Moreover, we consider finite particle distances in the dissipative range of turbulence. This is a two-particle problem, more general than the question of how spatial gradients of angles evolve, advected by turbulence. The latter dynamics is a singleparticle problem; it refers to an initially smooth manifold in phase space, as does the question of how the curvature of a material surface evolves in turbulence [51].
An open question is to find the form of the joint distribution P(R, δψ ) of separations R and angles δψ for finite . Our numerical results indicate that the distribution has algebraic tails too and that the power-law exponents vary as a function of particle shape. A related open problem is to derive the distribution in three spatial dimensions.
In the experiment it is difficult to separate rods that overlap in the camera images. In Fig. 2 we therefore show the experimental data only for distances larger than a rod length. It would be of great interest to obtain precise data at smaller distances, to systematically study the effect of hydrodynamic interactions [52,53] in the presence of nontrivial fluid flow.
A more far-reaching and more difficult problem is to understand how inertial effects change the patterns formed by larger particles. Particle inertia is relatively straightforward to take into account using the techniques reviewed in Ref. [2]. We anticipate that caustics [54][55][56][57] in the angular dynamics may increase the probability of large angles between nearby particles. For larger particles it also matters how the particles accelerate the surrounding fluid. There are corrections due to turbulent shears [10], and convective fluid inertia must matter for rapidly settling particles.
The relative alignment of particles in turbulence is a critical question in many scientific and engineering problems (including scattering of electromagnetic radiation from icy clouds [3], the dynamics of fiber suspensions [8], and plankton ecology [4]). Specifically, the relative alignment of approaching particles affects the rate at which they collide and possibly also collision outcomes. An 054602-12 ambitious long-term goal is to derive a theory for the collision rate between nonspherical particles in turbulence, with industrial fiber flows in mind [8], but also to model encounter rates of motile microorganisms [4,6] and of organic matter [58] in the turbulent ocean. In this context we must also ask how breaking of fore-aft symmetry changes the results summarized here. A first step in this direction is to use the present theory to determine the statistics of differences between the angular velocities ω = + n ∧ Sn of nearby particles. We expect that fractal steps in the director patterns are important, because the angular velocity depends explicitly on n and upon the shape parameter . However, this is just a starting point. The general problem is a very difficult one, yet important because of its wide range of applications. Our results show that the analysis of statistical models is a promising way of approaching this impactful question.  1) and (2)] for two nearby directors of separation R ≡ x 2 − x 1 and relative angle δψ ≡ ψ 2 − ψ 1 , we find, in two spatial dimensions, Here we have represented n as [cos ψ, sin ψ] T and we have expanded the fluid velocity u and the angular velocity ω in terms of small separations. Equation (A1) is expressed in dimensionless variables, x = x/η, t = t/τ , u = u/u 0 , and we have dropped the primes. Further, B = O + S, repeated indices are summed (Einstein summation convention), and jl is the two-dimensional Levi-Cività symbol. For small Kubo numbers, the fluid-velocity gradients and second derivatives in Eq. (A1) fluctuate rapidly and can be approximated by white noise. In this limit, we approximate the dynamics (A1) by a four-dimensional diffusion process. For the diffusion approximation to hold we must require not only that Ku → 0 (white-noise limit), but also that the change in n during one correlation time τ is much smaller than the magnitude |n| = 1. We find that this change is of order Ku n · Sn for large . Since the magnitude of n · Sn is of order unity in the dimensionless variables adopted in Eq. (A1), we must require that | | Ku 1 for the diffusion approximation to hold.
We calculate the drift and diffusion coefficients using an expansion in the Kubo number [2]. We use the dimensionless variables in Eq. (A1). We find that the drift coefficients vanish and that the diffusion coefficients are given by In the resulting diffusion equation, we change to polar coordinates R ≡ R[cos β, sin β] T . Since the flow is isotropic, we can remove one angular degree of freedom, leaving a three-dimensional diffusion equation in terms of R, δψ, and δβ ≡ β − ψ. This three-dimensional equation is hard to solve in general, but in the limit of large we have found the solution. Keeping only terms to highest order in , the diffusion equation simplifies considerably. In this limit, the steady-state diffusion equation reads Here f is a function of R, δψ, and δβ. It is related to the probability distribution P(R, δψ, δβ) by f = P/R. To solve Eq. (A3) we attempt a Fourier expansion This ansatz results in The solution of this equation has power-law tails for large δψ/R. The large-δψ/R asymptote of the solution can be written as a combination of two independent power laws in δψ/R: We must require that the tails are integrable to |δψ/R| = ∞. This implies that a m (R) = 0 for all values of m and that b m (R) = 0 for m = 0. Since the center of mass of the particles is advected in an incompressible flow, the marginal (spatial) distribution must be uniform. Therefore, we we must require that for small values of R. We now match the general m = 0 solution of (A5), to the asymptote (A7). As concluded above, we have a 0 (R) = 0, so we must match b 0 (R) = 3N /4R for small values of R, where N is a normalization factor. For the distribution P = f 0 R we thus obtain Eq. (12). Evaluating the moments |δψ| p R with this distribution yields Eq. (13) for small R and p = 1, and the coefficients are given by a p = 2 −p 3 p/2 cos(pπ/2) and b p = 2 −p+1 √ 3π p−2 /(p − 1) for p = 1. For p = 1 we find logarithmic corrections to power-law scaling, S 1 (R) ∼ R log R.
Using c 0 = 0 and d 0 = 1 the normalized probability distribution of angle gradients takes the form for large . The distribution of Y i is obtained by noting that the joint distribution P(Y, δα, δψ ) is uniform in δα and δψ. Therefore, the distribution P(Y, α) is uniform in α. It follows that P(Y 1 , This distribution has the P(Y 1 ) ∼ Y −2 1 tails mentioned in the main text.

APPENDIX C: CALCULATION OF THE LYAPUNOV DIMENSION D L
The Lyapunov dimension is a measure of the fractal dimension of the attractor in phase space. In two spatial dimensions, phase space is three dimensional (x 1 , x 2 , and ψ), so there are three Lyapunov exponents σ 1 σ 2 σ 3 . The sign of the maximal Lyapunov exponent σ 1 determines whether small separations grow (positive sign) or shrink (negative sign) exponentially. The signs of partial sums of the n upper Lyapunov exponents σ 1 + σ 2 + · · · + σ n determine whether n-dimensional subvolumes of phase space grow or shrink exponentially. Since the underlying flow is incompressible, the fractal dimension cannot be smaller than 2 (where σ 1 + σ 2 = 0 and σ 3 < 0) and it cannot exceed 3, the dimensionality of phase space (where σ 1 + σ 2 + σ 3 = 0). The Lyapunov dimension is defined as the linear interpolation between these limits [2,41]: To evaluate the Lyapunov exponents in Eq. (C1), we first note that two phase-space Lyapunov dimensions are given by the spatial Lyapunov exponents σ spatial 1 because the spatial dynamics is not influenced by the angular dynamics. Also σ spatial 2 = −σ spatial 1 . This equality follows from incompressibility of the flow.