aa r X i v : . [ phy s i c s . b i o - ph ] M a y Model flocks in a steady vortical flow
A. W. Baggaley
1, 2 School of Mathematics and Statistics, Newcastle University, Newcastle upon Tyne, NE1 7RU, UK Joint Quantum Centre Durham-Newcastle
We modify the standard Vicsek model to clearly distinguish between intrinsic noise due to imperfectalignment between organisms, and extrinsic noise due to fluid motion. We then consider the effectof a steady vortical flow, the Taylor Green vortex, on the dynamics of the flock, for various flowspeeds, with a fixed intrinsic particle speed. We pay particular attention to the morphology of theflow, and quantify its filamentarity. Strikingly, above a critical flow speed there is a pronouncedincrease in the filamentarity of the flock, when compared to the zero flow case. This is due to thefact that particles appear confined to areas of low vorticity; a familiar phenomena, commonly seenin the clustering of inertial particles in vortical flows. Hence the cooperative motion of the particlesgives them an effective inertia, which is seen to have a profound effect on the morphology of theflock, in the presence of external fluid motion. Finally we investigate the angle between the flowand the particles direction of movement and find it follows a power law distribution.
PACS numbers: 47.32.Ef,47.63.-b,87.18.-h
I. INTODUCTION
Flocks of birds, schools of fish or swarms of insects are allexamples of self-organised aggregations, offering some ofthe most spectacular displays of motion in nature. Thecomplexity of the animals moving together, seemingly asone, has intrigued observers for centuries. Indeed, theability of a flock of birds to change direction apparentlysimultaneously and spontaneously led Edmund Selous,who published numerous books on natural history, to be-lieve birds to be telepathic [1].To mathematically model animal aggregations it is typ-ical to use models based on self-propelled particles, suchas the Vicsek model [2]. Here the direction of motion istypically based on the idea that flocking animals tend toalign themselves to others within a local “radius of in-teraction” or based on a number of nearest-neighbours,as observed in starling flocks [3]. More complex mod-els with additional short range repulsion have also beeninvestigated [4], and a recent study extended the Vicsekmodel to include decision making, needed for example bybirds to make a synchronised landing [5].We note that this is not the only approach with a num-ber of interesting studies based on nonlinear partial dif-ferential equations [6]. In addition a recent study byBallerini et al. [7] questioned the validity of approachesbased on the standard Vicsek model, arguing that the ra-dius of interaction should be based on a topological dis-tance, rather than the metric distance used in [2]. How-ever we believe that the standard Vicsek model [2] isthe ideal starting point to investigate the effect of ex-ternal fluid motions on flock dynamics and morphology,and note that the behaviour of marching locusts couldbe modelled using a similar approach [8]. There are alsointeresting extensions to model human crowd stampedes[9].One of the drawbacks of the standard Vicsek model isthat all sources of noise in the system arise in a singleterm, which is both temporally and spatially uncorre- lated. It is perhaps more natural to distinguish betweentwo sources of noise, which we would expect to arise insuch systems. The first is intrinsic noise, due to the factthat animals will never perfectly align. Here we see noissue with approach used in [2] where uncorrelated noiseseems an appropriate approximation. However animalsmove in fluid environments, and so will be subject to ex-ternal noise, due to the motion of the fluid. In almostall natural flows a combination of large lengthscales andsmall viscosities mean that the flow is turbulent [10].Hence the external forcing will be subject to complexspatio-temporal correlations which are not present in thestandard Vicsek model.A recent study by Khurana & Ouellette [11] soughtto address this issue by considering a slight modificationto the Vicsek model, in an unsteady flow made up of asummation of random Fourier modes, with a prescribedEnergy spectrum consistent with the Kolmogorov spec-trum of classical turbulence. The focus of [11] was thestability of flocks, particularly in the limit of a small ra-dius of interaction, where they showed an external forcingwith spatio-temporal correlations was far more efficientin destabilising the flock. Here our interest is in the effectof fluid motion on coherent flocks, with larger values forthe radius of interaction.One may be tempted, therefore, to surmise that it isthen necessary to consider model flocks in fully turbu-lent flows. However, implementing a full direct numeri-cal simulation of the Navier-Stokes equation, necessary toproduce a turbulent velocity field is a formidable compu-tational challenge [12]. Moreover interpretation of resultsis not necessarily straightforward or unambiguous, giventhe strongly nonlinear nature of fully developed turbu-lence. Even the ‘reduced’ KS model considered in [11]is beyond what we wish to consider here, and it fails toreproduce many of the important features of true turbu-lence, such as intense vortical structures [13]. There isa long tradition in fluid dynamics of using simple ‘toy’models, to probe the underlying physics of a problem, x -3 -2 -1 0 1 2 3 y -3-2-10123 x -3 -2 -1 0 1 2 3 y -3-2-10123 FIG. 1: Snapshots of the Vicsek model in the statistically steady-state (left η = 0 . , R = 0 .
2, right η = 0 . , R = 0 . with a certain level of control. This is the approach wetake here, turning to the Taylor-Green vortex [14] as anideal toy flow. Indeed this flow has been applied to fieldsas diverse as swimming microorganisms [15], the settlingof inertial particles [16, 17] and Magnetohydrodynamics[18].The aim of this paper is to investigate the effect ofa simple vortical flow on a Vicsek flock, paying partic-ular attention to the flock morphology. We note thatthe shapes of schools of fish and flocks of birds have re-ceived both theoretical [19, 20] and experimental atten-tion [19, 21]. However less attention has focused on therole that external fluid motion could play in affectingflock morphology.The plan of the paper is as follows: In Section II wedescribe the standard Vicsek model. We also introducestatistical measures of the shape and structure of theflock, namely Ripley’s K function and Minkowski func-tionals. In Section III we extend the model to includea smooth external forcing due to a steady vortical flow.Detailed comparisons of the statistical measures of shapeand structure are made and we also study the angle be-tween the particles and the flow. We close with a discus-sion of the results in Section IV and highlight possibledirections for future studies. II. THE VICSEK MODEL
We begin by introducing the Vicsek model [2], a well-known analogue of the Ising model of ferromagnetism.We consider N particles in a two-dimensional square box(with sides of size L ), with periodic boundary condi-tions. Whilst the systems we are interested in are three-dimensional, by working in a low dimensional space we r h ˆ L ( r ) i -0.4-0.200.20.40.60.8 R = 0 . R = 0 . R = 0 . R = 0 . R = 1 . FIG. 2: (Color online) The time averaged Besag’s function h ˆ L ( r ) i , Eq. (5) plotted as a function of spatial scale r forsimulations with η = 0 .
2. ˆ L ( r ) > L ( r ) < gain a huge computational speedup which allows for athorough examination of parameter space. A follow upstudy to investigate flocking in three-dimensional timedependent flows is planned in the near future. Each par-ticle has a position x i ( t ) and an intrinsic, self-driven,velocity v i ( t ). All particles are assumed to move withthe same speed, V , and a particles intrinsic velocity isdetermined by v i = ( V cos( θ i ) , V sin( θ i )) , (1)where θ i determines the direction the particle moves in.Key to the model is that θ i is periodically (at each timeincrement) determined from the average of the particle’sown direction, plus the directions of its neighbours withina critical radius, R , such that [22] θ i = atan2 N X | x i − x j | 1] and η is the intensity of the noise. From these sim-ple rules Vicsek et al. [2] observed the emergence of col-lective flocking or swarming of the particles. The globalorder of the system can be characterised by computing, ψ η,R ( t ) = 1 N V (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) N X i =1 v i (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) , (3)which is a convenient measure to establish if the systemhas reached a statistically steady state. We follow [2] andmodel a system of N = 500 self-propelled particles in adomain of size L = 2 π with a speed V = 1. Particles areevolved according to an Euler scheme such that x i ( t + ∆ t ) = x i ( t ) + ∆ t v i ( t )where at each tilmestep v i ( t ) is updated according toEqns. (1) & (2). Note ζ i is drawn randomly at eachtimestep and we take ∆ t = 0 . 1. For large value of η co-operative motion is not observed, as the system is domi-nated by noise, hence we restrict our parameter values to η ∈ [0 , . R ∈ [0 , k -d tree [23].Simulations are evolved for 200 time units. Within ashort period of time ( ≈ 20 time units) ψ η,R saturatesto a statistically steady value, which for the low noiseintensity considered here is close to or equal to unity,indicating substantial global alignment of the particles.Snapshots of the particles for simulations with η = 0 . R = 0 . , . R plays in determining the structure of the flock is clear.In order to study this effect more quantitatively we useRipley’s K function, a statistical pattern analysis methodused as a measure of spatial clustering. For N particlesthis is defined asˆ K ( r ) = L N X i = j I ( | x i − x j | < r ) N − , (4)where I is the indicator function (1 if the Euclidean sep-aration is less than r , 0 otherwise), and r ∈ [0 , L ]. Notethat this can easily be adapted to account for the periodicboundary conditions used in this study. Particle cluster-ing results in ˆ K ( r ) increasing faster than if particles aredistributed in a spatially random manner, that is, if theyfollow a Poisson distribution. Ripley’s K function for aPoisson-distributed data set takes the form ˆ K ( r ) = πr . x -1 0 1 2 y -2.5-2-1.5-1-0.500.511.5 FIG. 3: (Color online) An example of a flock formed from thestandard Vicsek model. Red (filled) points indicate wherethe separation between that point and its nearest neighbouris less than ℓ/ 2, where ℓ = L/ √ N is the expected separationbetween random points. Blue (open) points are those whichdo not satisfy this criteria. The black structure shows the α -shape [24] constructed from the red points, which is usedto quantify the morphology of the flock. For a linear scaling of Poisson-distributed data, it is com-mon to normalise ˆ K ( r ) and define Besag’s function asˆ L ( r ) = p K ( r ) /π − r. (5)The advantage of Besag’s function is that it has a simpleinterpretation, it is zero for randomly distributed points,takes positive values for particles clustered over the spa-cial scale r and is negative if the particles are dispersedover a given scale. To increase our statistics we applytemporal averaging, and use angled brackets hi to de-note temporal averaging in the statistically steady state(where ψ η,R has saturated); this is done every 10 time-steps. Figure 2 shows the averaged Besag’s functions forfive different values of interaction radius R ( η = 0 . h ˆ L i to the rightis observed, as well as a broadening of the values of r where h ˆ L i > 0. Both are consistent with clustering ofparticles over larger scales with increasing R . We wantto emphasise that this is not a particularly novel result,however it is important to establish ‘baseline’ statisticsof the model, before studying the effects of flocking inthe presence of a structured external flow.It is also instructive to introduce a statistic which mea-sures the characteristic size of the flock. Whilst there isno unique way to do this, we find it convenient to mea-sure D = arg max r h ˆ L ( r ) i , (6)the location of the maximum value of the time averagedBesag function. The variation of h D i will be investigatedlater when we introduce a modified Vicsek model to ac-count for an external flow. x -3 -2 -1 0 1 2 3 y -3-2-10123 ω -0.8-0.6-0.4-0.200.20.40.60.8 FIG. 4: (Color online) The Taylor-Green vortex flow. Thepseudocolour plot displays the flow’s vorticity, Eq. (9), blackarrows indicate the velocity field, Eq. (8). We also study the morphology of the flock, and we arenot aware of any study which has previously done so.To do this one must first compute the shape of the flock.Again there is no unique way to do this and it is clear thatresults will be somewhat dependent on the method thatis used. Here our goal is not to perform a detailed analy-sis of flock morphology, for example to compare with thatobserved in bird flocks, insect swarms etc. But insteadto quantify the change in morphology when an externalflow acts on the particles. Hence we choose a compu-tationally efficient method, and apply this consistentlyin all simulations. We use the particle positions to de-fine an α -shape [24] (a generalisation of the convex hull).Our procedure is the following: Given N particles in adomain of area L , the expected separation between ran-dom points is ℓ = L/ √ N . Hence we define a subset ofpoints y i = { x i : | x i − x i ′ | < γℓ } , where i ′ denotes theclosest particle to x i . We take γ = 0 . 5, which leaves theset of points whose distance to the nearest neighbour is nogreater than half the expected separation value. This isto eliminate the contribution of the few points which arenot part of the coherent flock (which does occur for thelarger values of η used in this study). Note our results arenot sensitive to the choice of γ (within a sensible range),indeed results are virtually identical if γ ∈ [0 . , . α -shape is subsequently constructed from the points y i , using a standard package in MATLAB . The resultsof this process can be seen in Fig. 3. This is then used tocompute the morphological properties of the flock usingMinkowski functionals, which we discuss next.Minkowski functionals measure the topological and ge-ometrical features of a given shape or structure. Had-wiger’s theorem [25] states that in a d -dimensional space,there exist d + 1 quantities which completely describethe morphology of a given structure. In 2D, the threeMinkowski functionals of a closed contour are its enclosedarea S , perimeter P , and the Euler characteristic χ [26], but for the purposes of this study it is most convenientto analyse the dimensionless filamentarity F = ( P − πS ) / ( P + 4 πS ) (7)By definition, F ∈ [0 , F = 0 for a circle and F = 1 for a (not necessarily straight) line segment. Both S and P can be efficiently computed from the α -shape,and so we are able to compute the evolution of F intime, and report an average value for the filamentarity, h F i , as we did for Besag’s functions. Again we refrainfrom reporting results here, and present the bulk of re-sults once we have discussed a modification of the model.Note that the boundary ∂ S α of the α -shape is a subsetof the Delaunay triangulation of S , hence it is plausi-ble that we over estimate P and hence F due to smalltriangular artefacts as can be seen in Fig. 3. Again weemphasise that we are not attempting a detailed analysisof flock morphology, to compare with natural flocks. Butto quantify the change in morphology when an externalflow acts on the particles using a consistent methodologyacross all simulations. III. SELF-PROPELLED PARTICLES IN THETAYLOR GREEN FLOW We now turn our attention to Vicsek flocks in the pres-ence of external noise, due to fluid motion. As discussedin section I, we forgo the computational expense, andcomplexity of a fully turbulent flow. Instead we turn toa steady vortical flow, the Taylor Green (TG) vortex [14],defined as v f ( x ) = ( u f , v f ) = V f (sin( x ) cos( y ) , − cos( x ) sin( y )) , (8)where x = ( x, y ). The vorticity field is given by ω = ∇ × v f = 2 V f sin( x ) sin( y ) , (9)and the flow is incompressible ( ∇· v f = 0). V f is a scalingparameter which can be adjusted to modify the relativeintrinsic particle speed to that of the background flow.Figure 4, displays the vorticity field of the TG flow withthe corresponding velocity field.A perhaps overly simplistic approach to modelling thesystem, would be to simply allow the particles to be ad-vected by the flow, such that the equation of motion forthe particles is modified to becomed x i d t = ( V cos( θ i ) + u f ( x i ) , V sin( θ i ) + v f ( x i )) . (10)However, it is not clear to us how Eq. (2) can still re-main valid. We would assume that ‘particles’ rather thanorienting themselves to match the direction the otherparticles are oriented towards, they would instead orientthemselves to the direction of motion of nearby particles. x -3 -2 -1 0 1 2 3 y -3-2-10123 -0.8-0.6-0.4-0.200.20.40.60.8 x -3 -2 -1 0 1 2 3 y -3-2-10123 -0.8-0.6-0.4-0.200.20.40.60.8 FIG. 5: (Color online) Snapshots of the Vicsek flock with external fluid forcing, in the statistically steady-state (left η =0 . , R = 0 . 2, right η = 0 . , R = 0 . 8) with V f = 0 . 5. A pseudocolour plot (made slightly translucent to improve visualisation ofthe particles) displays the flow’s vorticity ω , Eq. (9). Arrows indicate the particles direction of motion. Note the concentrationof particles in regions where ω ≃ r h ˆ L ( r ) i -0.2-0.100.10.20.30.40.5 R = 0 . R = 0 . R = 0 . R = 0 . R = 1 . FIG. 6: (Color online) The time averaged Besag’s function h ˆ L ( r ) i , Eq. (5) plotted as a function of spatial scale r forsimulations with V f = 0 . 5. ˆ L ( r ) > L ( r ) < Such an approach was also taken in [11]. Hence we arriveat a modified form for the equation used to update θ , θ i = atan2 N X | x i − x j | 1. Unlikefor the standard model, we use a 3 rd -order Runge-Kuttatime-stepping scheme, to allow for the additional com-plexity of the flow, however θ is only updated at thebeginning of a timestep. Again simulations ( η ∈ [0 , . R ∈ [0 , ψ η,R to ensure statistics are only computed in a statisti-cally steady value. We consider three values of V f = 0 . . . 5, with a fixed value V = 1 for the particlespeed. Snapshots of the particles for simulations with η = 0 . R = 0 . , . V f = 0 . R , that the morphology of the flock is dramatically af-fected by the influence of the flow. We note also thatparticles tend to be found in areas of low vorticity, thisis not seen (and would not be expected) if we considertracer particles V = 0. We shall discuss this feature laterin the text.In all simulations we find the statistical properties ofthe flocks are nearly independent of magnitude of theintrinsic noise, and so for all quantities reported belowwe perform additional averaging over η . Figure 6 showsthe averaged Besag’s functions (Eq. (5)) for five differ-ent values of interaction radius R , with V f = 0 . 5. Againa noticeable (although less pronounced when comparedwith Fig. 2) migration of the peak of h ˆ L i to the right is R < D > V f = 0 V f = 0 . V f = 0 . V f = 0 . R < F > V f = 0 V f = 0 . V f = 0 . V f = 0 . FIG. 7: (Color online) (left) The time averaged characteristic size of the flock h D i (Eq. (6)), plotted as a function of the radiusof interaction R , for varying flow velocities V f . (right) The corresponding plot of the variation of the flock filamentarity h F i (Eq. (7)). Note a profound change is observed in the flock morphology ( R ≥ . 6) when V f ≥ . φ P D F ( φ ) -2 -1 φ π P D F ( φ ) φ − . FIG. 8: (Color online) The probability density function(PDF) of φ , the angle between the flow direction and theparticles intrinsic velocity, for the simulation with V f = 0 . η = 0 . R = 0 . φ ) ∼ φ − . displayed as a red (dashed)line. observed. However due to the periodic boundary condi-tions a secondary peak where h ˆ L ( r ) i > r > π canbe seen for all values of R . This indicates that the flowis forcing the flocks to form an elongated shape, tightlyclustered in one direction, but coherent over the scale ofthe box, as we see in Fig. 5.This is also visible if we inspect plots of h D i (Eq. (6))against R for the values of V f considered, Fig. 7 (left).A clear reduction in the typical distance of clustering h D i is seen (at least for larger values of R ) as the flowvelocity increases. The effect of the flow on the flocksmorphology, namely h F i (Eq. (7)), can be seen in Fig. 7(right). At low flow speeds ( V f = 0 , . 1) increasing the radius of interaction leads to a dramatic decrease in thefilamentarity, as the flock tends to become close to cir-cular. However above some critical value of V f , there isa dramatic transition, and even for large values of R thestructures exhibit large filamentarity. It is important tonote that this is visible even when the flow speed is onlyone-fifth of the particles intrinsic swimming speed.Finally we turn our attention to the relationship be-tween the particles intrinsic direction of motion, and thedirection of the local flow. We define φ i = cos − (cid:18) v i · v f ( x i ) V | v f ( x i ) | (cid:19) , (12)from this we use kernel density estimation [27] to esti-mate the probability density function (PDF) of φ . Thisis displayed in Fig. 8. We observe empirically thatthere is a power law relationship between the angle be-tween the flow and the particles intrinsic velocity, i.e.PDF( φ ) ∼ φ α . Of course it is then natural to ask if thispower law is universal or depends on η , R and/or V f .We address this question directly in Fig. 9, for the twolarger flow velocities V f = 0 . , . 5. It is clear that φ isindependent of the intrinsic noise η , but depends on boththe flow speed V f and the radius of interaction R . Thebehaviour is as one may naively expect with a much shal-lower power law relationship for large values of R , andsmaller values of V f , where one may expect more scat-ter between the flow direction and the particles intrinsicdirection. IV. DISCUSSION The main goal of this work was to understand the ef-fect of an external forcing flow, made up of steady coun-terrotating vortices, on the standard Vicsek model. Todo this we made a slight modification to the model and R α -3.2-3-2.8-2.6-2.4-2.2-2-1.8-1.6-1.4 η = 0 . η = 0 . η = 0 . η = 0 . η = 0 . R α -4-3.5-3-2.5-2-1.5 η = 0 . η = 0 . η = 0 . η = 0 . η = 0 . FIG. 9: (Color online) The power law scaling, α , of the probability density function (PDF) of φ (PDF( φ ) = φ α ) where φ is theangle between the flow direction and the particles intrinsic velocity. (left) V f = 0 . 2, (right) V f = 0 . 5, plotted for varying radiiof interaction R . The slope is clearly dependent on both R and V f , but appears independent of the intrinsic noise η . x -3 -2 -1 0 1 2 3 y -3-2-10123 -2-1.5-1-0.500.511.52 FIG. 10: (Color online) Snapshots of the Vicsek flock withexternal fluid forcing given by Eq. (13), with η = 0 . R = 0 . V f = 0 . 5. A pseudocolour plot (made slightly translucentto improve visualisation of the particles) displays the flow’svorticity ω . Arrows indicate the particles direction of motion. then compared statistically properties of the flocks, pay-ing particular attention to morphological properties. Oneof the striking features of the flocks with an external flow V f ≥ . v f ( x ) = V f ν X m =1 ( A m × k m cos φ m + B m × k m sin φ m ) , (13)where φ m = k m · x , k m are wave vectors, A m and B m are random unit vectors (orthogonal to k m ) and ν is ascaling parameter chosen such that the root mean square(RMS) value of v f is V f . Taking k m = m creates therandom large scale (solenoidal) flow whose vorticity fieldis displayed in Fig. 10. A simulation of the Vicsek flockin the presence of this random fluid flow is performedtaking η = 0 . R = 0 . V f = 0 . 5; a snapshot of thesystem in the steady-state is also displayed in Fig. 10.The elongated structure of the flock visible in the sim-ulations using the TG flow is also found in this randomflow, which gives some indication of the robustness of theresults presented here.Of course it is an open question as to the extentto which these results apply to real animal flocks andswarms. However it is both interesting and importantto understand how structures in the flow can affect boththe shape and fluctuations (which we have not consid-ered here) of the shape of animal flocks. As mentionedearlier in the text we have neglected how rotation of theparticles affects their dynamics, which could provide aninteresting extension to the model. Natural next stepsare to understand the role that temporal fluctuations of the flow play, before turning attention to more realisticflows, such as those generated by direct numerical simu-lation of the Navier Stokes equation. [1] E. Selous. Thought-transference (or What!) in Birds.Constable Limited, 1931.[2] T. Vicsek, A. Czir´ok, E. Ben-Jacob, I. Cohen, andO. Shochet. Novel type of phase transition in a systemof self-driven particles. Phys. Rev. Lett., 75:1226–1229,Aug 1995.[3] G. F. Young, L. Scardovi, A. Cavagna, I. Giardina, andN. E. Leonard. Starling flock networks manage uncer-tainty in consensus at low cost. PLoS Comput. Biol.,9:e1002894, 01 2013.[4] I. D. Couzin, J. Krause, R. James, G. D. Ruxton, andN. R. Franks. Collective memory and spatial sorting inanimal groups. J. Theor. Biol., 218(1):1 – 11, 2002.[5] Kunal Bhattacharya and Tams Vicsek. Collective deci-sion making in cohesive flocks. New J. Phys., 12, 2010.[6] J. K. Parrish and L. Edelstein-Keshet. Complexity, pat-tern, and evolutionary trade-offs in animal aggregation.Science, 284(5411):99–101, 1999.[7] M. Ballerini, N. Cabibbo, R. Candelier, A. Cavagna,E. Cisbani, I. Giardina, V. Lecomte, A. Orlandi,G. Parisi, A. Procaccini, M. Viale, and V. Zdravkovic.Interaction ruling animal collective behavior depends ontopological rather than metric distance: Evidence from afield study. Proc. Natl. Acad. Sci. U.S.A., 105(4):1232–1237, 2008.[8] J. Buhl, D. J. T. Sumpter, I. D. Couzin, J. J. Hale, E. De-spland, E. R. Miller, and S. J. Simpson. From disorder toorder in marching locusts. Science, 312(5778):1402–1406,2006.[9] D. Helbing, I. Farkas, and T. Vicsek. Simulating dynam-ical features of escape panic. Nature, 407:487–490, 2000.[10] F. Uriel. Turbulence: The Legacy of A. N. Kolmogorov.Cambridge University Press, 1995.[11] N. Khurana and N. T. Ouellette. Stability of model flocksin turbulent-like flow. New J. Phys., 15(9):095015, 2013.[12] T. Ishihara, T. Gotoh, and Y. Kaneda. Study of high–reynolds number isotropic turbulence by direct numeri-cal simulation. Annu. Rev. Fluid Mech., 41(1):165–180,2009.[13] Zhen-Su She, Eric Jackson, and Steven A. Orszag. Inter-mittent vortex structures in homogeneous isotropic tur-bulence. Nature, 344(6263):226–228, 03 1990.[14] G. I. Taylor. Lxxv. on the decay of vortices in a viscousfluid. Philos. Mag., 46(274):671–674, 1923.[15] W. M. Durham, E. Climent, and R. Stocker. Gyrotaxisin a steady vortical flow. Phys. Rev. Lett., 106:238102,Jun 2011. [16] M. R. Maxey and S. Corrsin. Gravitational settling ofaerosol particles in randomly oriented cellular flow fields.J. Atmos. Sci., 43(11):1112–1134, 2015/03/01 1986.[17] L. Bergougnoux, G. Bouchet, D. Lopez, and E. Guazzelli.The motion of solid spherical particles falling in a cellularflow field at low stokes number. Phys. Fluids, 26(9), 2014.[18] C. Nore, M. E. Brachet, H. Politano, and A. Pouquet.Dynamo action in the taylor–green vortex near threshold.Phys. Plasmas, 4(1):1–3, 1997.[19] C. K. Hemelrijk, H. Hildenbrandt, J. Reinders, and E. J.Stamhuis. Emergence of oblong school shape: Modelsand empirical data of fish. Ethology, 116(11):1099–1112,2010.[20] C. K. Hemelrijk and H. Hildenbrandt. Schools of fishand flocks of birds: their shape and internal structure byself-organization. Interface Focus, 2012.[21] B. L. Partridge, T. Pitcher, J. M. Cullen, and J. Wilson.The three-dimensional structure of fish schools. Behav.Ecol. Sociobiol., 6(4):277–288, 1980.[22] Note the use of the atan2 function, the four-quadrantinverse tangent tan − ( y, xy, x