Zhi Lin , , Sirui Zhu , Lingyun Ding

a School of Mathematical Sciences, Zhejiang University, Hangzhou, Zhejiang 310027, China

b Department of Mathematics, University of North Carolina, Chapel Hill, NC 27599, United States

Keywords: Mixing and transport Passive scalar Anisotropic diffusion

ABSTRACT We consider a fluid stirred by the locomotions of squirmers through it and generalize the stochastic hy- drodynamic model proposed by Thiffeault and Childress, Phys. Lett. A (2010) and Lin et al., J. Fluid Mech. (2011) to the case in which the swimmers move in anisotropically random directions. A non-diagonal effective diffusivity tensor is derived with which the diffusive preference of a passive particle along any given direction can be computed to provide more details of the phenomena beyond scalar statistics. We further identify a fraction from the orthogonal decomposition of the drift-induced particle displacement to distinguish the underlying nonlinear mixing mechanism for different types of swimmers. Numerical simulations verify the analytical results with explicit examples of prescribed, anisotropic stirring motions. We also connect our formulation to several measures used in clinical medical research such as diffusion tensor imaging where anisotropic diffusion has a significant consequence.

Scalar advection due to drift, namely, the permanent displace- ment of elements disturbed by the passage of a body through a bulk fluid, is a ubiquitous and important mechanism that enhances heat and mass transport in various natural and industrial pro- cesses. Recently, academic interests in drift-enhanced mixing has been sparked by their potential importance in geophysical and bi- ological contexts [1–4] .

A key feature in these fluidic environments is stratification that prohibits the mixing and transport of various quantities in certain directions [5,6] . Yet, field measurements suggest that this hindrance is somehow compensated by motions overcoming the cross-sectional potential energy gradient and stability between wa- ter layers [7,8] . In this paper we focus on the potential contribu- tion of the Darwin drift due to swimming bodies [9] as suggested to be a significant effect in the swimming-induced mixing in the ocean [10–12] .

To define and compute the effective diffusivity for the mass displacement due to a colony of submerged swimmers, Thiffeault et. al. followed the classical kinetic theory of diffusion and pro- posed an original stochastic hydrodynamic model that derives the effective diffusivity of a passive scalar drifted by random swim- mers [13,14] . In particular, the position of target particlexMafterMencounters withsquirmers, which are classical and accurate mod- els for potential or Stokesian swimmers in real life [15–19] , is com- puted as a discrete sum of sequential, independent drifts:

HereΔλ(ak,bk)is the displacement from each “kick”in the ran- dom direction ?rk,akandbkare decomposed impact parameters in a co-moving framework of reference. Under isotropic and non- interacting assumptions the flow-enhanced diffusivity is increased by about 100 times over molecular diffusion. The model has later been substantiated in observations from experiments [20–22] and from simulations [23–25] .

Among the extensive body of aforementioned research, a com- mon but oversimplifying assumption is isotropy. As mentioned ear- lier, a signature of the environments in which the swimmers oper- ate is the anisotropic distribution of important physical quantities and the resulting behavioral parameters [26–28] . Therefore in this paper we analyze and simulate the anisotropic squirming and the resulting mixing enhancements as a generalization to the isotropic model. Here, the full effective diffusivity tensor is computed from the mean squared particle displacement in two or three dimen- sions for an anisotropic probability distribution prescribed for the random swimming directions.

We study the mixing and transport of a passive scalar field stirred byNnon-interacting swimmers uniformly distributed in a two- or three-dimensional fluidic domainΩ. Each swimmer moves along a line segment of lengthλin a random direction at a con- stant speedU. The Darwin drift of a particle in an encounter with

one swimmer is therefore the vector

whereUis the random velocity vector of magnitudeUandx(s)is the particle trajectory starting at the pointη. It can be shown thatΔλ(η,U)is independent ofUfor largeλ[14,29] .

Consider the probability density function (PDF),pa(k), for the random swimming direction vectora=U/U. It was assumed to be isotropic in previous work, namely,pa(k)= 1/2 πin 2D orpa(k)= 1/4 πin 3D. To extend the applicability of the model by allow- ingpa(k)to be non-constant, we first note that the mean squared particle displacement due to one instance of Darwin drift over all swimming directions is

Furthermore, under the assumptions of diluteness and indepen- dent, identically distributed (i.i.d.) swimming directions, the total particle displacementdue toNswimmers/encounters is sim- ply the sumwith the displacementfrom theith encounter also being i.i.d. Therefore,

wheren=N/Vis the number density of the swimmers. In the largeλlimit as we will show in later sections, this average is fur- ther reduced tonregardless of whether the distribu- tionpa(k)is constant and the effects of anisotropy only emerge in the full diffusivity tensor while its trace, the total diffusivity, is in- variant. Therefore the total diffusivity is only a first order descrip- tion of the effective diffusion in question and more details can only be recovered with further analysis.

In this section, a standard procedure in statistical mechanics will be adopted to compute the full diffusivity tensor for analysis and for simulation.

Following the classical It?approach the effectively particle dif- fusion in question can be characterized by a second-order tensor which we denote ashenceforth. It is derived from statistical av- erages of asymptotic particle displacements as the following,

To facilitate future derivation, we now orthogonally decompose the displacement vectorΔλ(η,U)in theco-moving frameworkof reference aligned with the swimming path. For example, in two dimensions we have

whereais the swimming direction as we described before and the unit vectorbis perpendicular toa. Combined with Eq. (3) for the particle displacements, the full effective diffusivity tensor in Eq. (4) is then reduced to

An encounter between the particle and the swimmer is illustrated in Fig. 1 .

Fig. 1. Schematic diagram of the 2D decomposition of particle displacement in a single encounter.

Fig. 2. Spatially averaged cross-stream fraction of particle displacements, ξ, as a function of swimming length, λ, for three types of swimmers.

Fig. 3. Polar plots of two-dimensional probability density functions f i (θ) , i = 1 , 2 , 3 for different swimming preferences and their induced directional effective diffusivities κf (θ) .

And we further replace the probability distribution functionpa(k)for the random swimming direction with an equivalent but simpler form, namely,

Then the diffusivity tensor Eq. (4) in 2D can be readily computed as

where the parameters

are Fourier components of the PDF and

which is essentially the spatially averaged cross-stream fraction of particle displacements. This fraction is specified by the fluid dy- namical properties of the flow field generated by one individual swimmer and is also dependent on the swimming lengthλ.

Similarly, in a 3D spherical coordinate system a function of two angular variables,h(θ,?),θ∈ [0,2 π),?∈ [0,π] , should be used to describe the probability distribution of random swimming direc- tions. In particular, whenhis separable, namely,

the three-dimensional version of the diffusivity tensor Eq. (4) is

where the parameters are computed by

For a non-separableh(θ,?), a closed-form expression of the tensor Eq. (4) is not available to our knowledge due to the compli- cated structure in the flow field. But as we will show in the next section, we can decompose the distribution as a linear combination of separable components and identify the anisotropic signatures in each component. Furthermore, with the above formulation the nu- merical computation for the tensor given any general distribution is then straightforward by superposition.

To understand how directionally preferred stirring translates into anisotropic diffusion, we now establish a quantitative connec- tion between the diffusivity tensor Eq. (6) and the swimming mo- tion by introducing an analogous swimmer mobility tensor as

and again,ais the unit vector in the swimming direction. It is then obvious that the streamwise component of the diffusivity Eq. (6) is a constant multiple of. To evaluate the significance of the trans- verse contribution we rewrite Eq. (6) as

with the average cross-stream fraction of particle displacements,ξ, defined in Eq. (10) .

Figure 2 demonstrates the behavior ofξfor three types of typ- ical swimmers: a cylinder in a potential flow (2D potential), a sphere in a potential flow (3D potential) and a Stokesian squirmer with the stresslet intensity set to 0.5 [14] . It can be seen that for potential swimmers,ξdecreases rapidly asλ/lincreases which in- dicates that the overall transverse contribution to the effective dif- fusion is negligible. In this case the diffusivity tensoris trivially proportional to the swimmer mobility tensor. However, for the Stokesian case,ξhas an asymptote of a significantly nonzero value and now we recover a nonlinear relation between the two as a result of indispensable cross-stream effects in the Darwinian drift. This distinction in terms of mixing mechanism can be viewed as a direct consequence of the difference between the mechanical and dynamical features of potential and Stokes flows.

In this section, we will see the explicit results for several pre- scribed distributions,pa(k)in both 2D and 3D cases. Moreover, we will introduce a scalar function to obtain a more straightforward, geometrical interpretation of the anisotropic effective diffusion.

Notice that in two dimensions we can always apply a coordi- nate rotation so that the diffusivity tensor Eq. (8) is equivalent to the diagonal form

in which the diagonal elements are the diffusivities along two coordinate axes, respectively. This linear transformation has been widely used by clinical and animal studies and the ratioκx/κywas used to characterize anisotropic diffusion [30] . In this spirit, we now reduce the rank-2 tensor to a single scalar function to have a geometrical and physical understanding to our problem. That is,

and it can be viewed as the diffusivity along any given direction specified by the angleθ.

Here the constants

and the dynamical parameterξ, the probabilistic parametersαfandβfhave been defined in Eqs. (9) and (10) .

Now we look at the results for the following angular functions describing the swimming preference in 2D as indicated in Eq. (7) :

Figure 3 is a comparison between the distributions Eq. (19) and the corresponding directional diffusivitiesκfi(θ),i= 1,2,3 . Here we have a uniform distribution of circular swimmers with num- ber densityn= 10?5each of which generates a 2D potential flow with parametersU=l= 1 andλ= 20 .

It can be observed that although the swimming preferences are vastly different and anisotropic, their diffusive consequences are similar indicated by the mild angular variations inκfi. In particular,κf2is essentially a rotated version ofκf1andf3 even implies an isotropic diffusivity. This should have been expected because diffu- sion is a smoothing process and the sharp variations infiwould be attenuated. On the other hand, we have seen in Eq. (9) that the controlling components infin terms of diffusive impact are sin(2θ)and cos(2θ)throughthe parametersαfandβf. Thereforef1andf2produce similar results. In contrast, the modes sin(3θ)and cos(3θ)inf2andf3have no significance in the effective diffu- sivity, a second-order statistics and can only contribute to higher- order statistics of the scalar field.

Alternatively a physical perspective of the above results can be provided by visualizing the evolution of an initial scalar concentra- tion being expanded under the effective diffusivity tensor Eq. (8) . With the same dynamical settings used in Fig. 3 , we numerically solve the effective heat equation in free 2D space

forfi,i= 1,2,3 and summarize the results in Fig. 4 . The initial, isotropic distribution of the passive scalar fieldθis

In column (a) of Fig. 4 ,θ0is plotted in grayscale superimposed with concentric contours. In columns (b) and (c), we can clearly see that the scalar field develops anisotropic structures forf1andf2astimeprogresses, highlighted by the elliptical contours. How- ever, forf3the plots are indistinguishable from isotropic diffusion since the effective diffusivity tensor in this case is merely a con- stant multiple of the identity matrix. These are all consistent with our previous analyses, especially with the plots forκfishown in Fig. 3 .

Fig. 4. Numerical simulations of the effective heat Eq. (20) . In the i th row, i = 1 , 2 , 3 , three snapshots of the evolution of the scalar concentration are shown for t = 0 , 50 and 100 respectively, with the diffusivity tensor Eq. (8) computed for the prescribed distribution f i (θ) given in Eq. (19) .

The diffusivity tensor in three dimensions can be computed with Eq. (12) if the swimming preference distribution is separable as

For a general, non-separable distribution, we may decompose it as

with the coefficientsak= 1,2,···chosen to satisfy renormalization conditions for a proper probability distribution. Consequently, the diffusion tensor would also be a linear combination of terms in the form of Eq. (12) . That is,

and

with probabilistic parameters

We now focus on one such component without the loss of gen- erality and drop the subscriptskhenceforth. Following a similar derivation, we can define a three-dimensional, directional diffusion coefficient as

for given anglesθ∈ [0,2 π),?∈ [0,π] where parametersci,i= 1,2,3,4 are constants connected to those in Eq. (25) by

For example, for distribution densities the corresponding directional diffusivities are

Figure 5 illustrates the density functions Eq. (28) and the di- rectional diffusivities Eq. (29) with similar features as in two- dimensional cases shown in Fig. 3 , especially when 2D projections onto thex?yplane are considered. Here the surfaces are rendered withx=rcosθcos?,y=rsinθcos?andz=rsin?where the ra- dial distanceris specified by the function value ofh(θ,?)or byκh(θ,?).

To validate the analytical results in previous sections we im- plement direct particle simulations. Similar numerical experiments as documented in Ref. [14] are performed with the same param- eter values used, with the only difference being that the swim- mers here are swimming in anisotropically random directions. For each scenario, the particle starts at the origin with 105realiza- tions in the domains(x,y)∈ [ ?50 0,50 0]2in two dimensions and(x,y,z)∈ [ ?50 0,50 0]3in three dimensions, respectively.

Figure 6 compares the mean squared displacements of the par- ticle stirred by anisotropic swimmers each of which generates a potential flow around a dipole for the probability distributions Eq. (19) . In each case, the displacement and its projections onto thex?andy?axis all grow linearly in time which indicate the dif- fusive nature of the process. In each panel, the black dashed line segment slightly below the curve of total displacement 〈R2〉 has a slope equal to the trace of the theoretical formula for diffusivity tensor Eq. (8) with numerical parameter values. The agreement be- tween simulation and analysis verifies that the modified stochas- tic hydrodynamic model we propose is an accurate model for the anisotropic stirring we started to investigate.

We notice that although the total mean squared displacement 〈R2〉 has the same growth rate for all three distributions given in Eq. (19) , the evolution of its axial projections varies. Forf2(θ)andf3(θ)the projected displacementsandare essen- tially equal. This similarity suggests that the second-order statis- tics alone does not sufficiently identify the signatures of the anisotropic dynamics and one needs a more detailed perspective such as the full diffusivity tensor, or the directional diffusivity function Eq. (17) . In contrast, a clear distinction betweenandcan be observed forf1(θ)attributed to the anisotropic dif- fusivity tensor Eq. (8) in this case. Three dimensional simulations provide similar insights for the distributions Eq. (28) and confirm the analytical calculations shown in Fig. 5 .

Fig. 5. 3D surfaces specified by the direction distributions: (a) h 1 (θ, ?) = sin ?, (b) h 2 (θ, ?) = and by their corresponding direc- tional diffusivities κh 1 (θ, ?) and κh 2 (θ, ?) in panels (c) and (d), respectively.

Fig. 6. Mean squared displacements of a passive particle under 2D anisotropic stirring as a function of time for three distributions: (a) , (b) f 2 (θ) =

In previous sections, we have derived the full diffusivity ten- sor and introduced an angular function to characterize the overall diffusive impact by a colony of swimmers with directional prefer- ences. In biological and clinical research, anisotropic diffusion has also been documented in a great body of literature [31,32] and this phenomenon allows technologies like diffusion tensor imag- ing (DTI) to focus on the diffusivity tensor and to uncover even more microstructural details in traditional MRI data for neu- ropathological diagnosis and treatment [33] . For practical purposes in these areas, the tensoris often further reduced to various scalar indices for measuring the extent of anisotropic structures with technical treatments aimed to avoid estimation biases [34].

Most of these measures are computed from the eigenvalues of the full diffusivity tensor, denoted asλi,i= 1,2,···,dinddi- mensions. Eigenvalues and their associated eigenvectors represent the effective diffusivities in the principle directions, respectively and are invariant under coordinate transforms. Common examples include fractional anisotropy (FA), relative anisotropy (RA) and vol- ume ratio (VR) which are defined as follows:

Table 1 (2D cases) and Table 2 (3D cases) list the values of these rotation-invariant measures for anisotropy for the distribu- tion functions Eqs. (19) and (28) . Similar to what was demon- strated in Figs. 3 and 4 , in two dimensions we recover the identical quantifications of anisotropy forf1 andf2 , whereasf3 is indistin- guishable from the isotropic case. Furthermore, the indices FA and RA seem to suggest a larger deviation from isotropy than VR does. This reinforces the applicability of the framework and analytics introduced in previous sections and suggests their potential abil- ity to reveal more details of dynamical features and mechanisms beyond such scalar measures as well as general second-order statistics.

Table 1 Comparison of scalar measures Eq. (30) for diffusive anisotropy for 2D distributions Eq. (19) .

Table 2 Comparison of scalar measures Eq. (30) for diffusive anisotropy for 3D distributions Eq. (28) .

In this paper, we propose an anisotropic extension to the stochastic hydrodynamic model for biogenic mixing motivated by the non-uniform spatial preference in the swimming motion of marine species, or any dilute colony of active swimmers in a flu- idic environment in general. In particular, the diffusivity tensor of a submerged passive scalar field is calculated and discussed in de- tail based on classical results in statistical fluid mechanics connect- ing mean squared particle displacements and Darwin drifts due to swimming bodies. These analyses and results provide a simpli- fied, theoretical framework and reference to explain the underlying mechanisms in various mixing phenomena observed in nature and in engineering.

Future work will be dedicated to taking into account more realistic factors to apply the model to broader practical settings studied in the literature. For example, we can further random- ize the parameter for the swimming pathλ[35] . In addition, preliminary results show that the anisotropic initial distribution of the swimmers, which is ubiquitous in nature [36] , also has a great impact on the mean-squared displacement of the particle. Moreover, we will strive to find a good quantification to char- acterize the anisotropy in three dimensions analogous to its 2D counterpart discussed, due to the challenges posed by the com- plicated interactions between the entries in the diffusion tensor Eq. (24) .

Declaration of Competing Interest

The authors declare that they have no known competing finan- cial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgement

This work was supported by the National Natural Science Foun- dation of China (Grant No. 12071429).