Coherent structures and extreme events in rotating multiphase turbulent flows

By using direct numerical simulations (DNS) at unprecedented resolution we study turbulence under rotation in the presence of simultaneous direct and inverse cascades. The accumulation of energy at large scale leads to the formation of vertical coherent regions with high vorticity oriented along the rotation axis. By seeding the flow with millions of inertial particles, we quantify -for the first time- the effects of those coherent vertical structures on the preferential concentration of light and heavy particles. Furthermore, we quantitatively show that extreme fluctuations, leading to deviations from a normal-distributed statistics, result from the entangled interaction of the vertical structures with the turbulent background. Finally, we present the first-ever measurement of the relative importance between Stokes drag, Coriolis force and centripetal forces along the trajectories of inertial particles. We discover that vortical coherent structures lead to unexpected diffusion properties for heavy and light particles in the directions parallel and perpendicular to the rotation axis.

The strength of rotation is measured by the Rossby number Ro = ( f k 2 f ) 1/3 /Ω, defined as the ratio of the rotation time, τ Ω = 1/Ω, and the flow time-scale, f k 2 f . Here f and k f are the input of energy and the wavenumber where the external forcing is applied (see table   1). The most striking phenomenon originated by the Coriolis force is the formation of intense and coherent columnar vortical structures (see Fig. 1), which has been observed in numerical simulations [15][16][17][18][19] and in experiments for rotating turbulence produced by an oscillating grid [9], for decaying turbulence [10][11][12], forced turbulence [14], and turbulent convection [24]. The appearance of these large-scale vortices is associated to a noticeable two-dimensionalization of the flow in the plane perpendicular to the rotation axis. Rotating turbulent dynamics with Rossby number O(1) is typical of many industrial and geophysical applications, but key fundamental questions are still open. These are mostly connected to the nature of the interaction between the two-dimensional vortical structures and the underlying fully three-dimensional anisotropic turbulent fluctuations, and to the way this impacts the Lagrangian dynamics of particles dispersed in the flow. In this paper, we empirically assess the Eulerian and Lagrangian statistical properties of rotating flow by using high resolution direct numerical simulations at unprecedented resolution. We present the first simulateneous study of Lagrangian and Eulerian properties, seeding the strongly rotating flow with billions of small-particles with and without inertia. In particular, we investigate statistical events much larger than the root mean squared fluctuations, measuring high order moments of velocity increments both along the rotation axis and in the perpendicular plane. To disentangle the statistical properties of the 2D structures from the underlying 3D turbulent background, we propose to decompose the velocity field on its instantaneous mean profile, obtained by averaging along the rotation axis, and on the fluctuations around it. We show that there exits a highly non-trivial entanglement among the vortical structures and the 3D background leading to a complex non-Gaussian distribution for both 2D and 3D components. Similarly, we quantify the singular role played by vortical structures for the preferential concentration of inertial particles' trajectories. We assess for the first time the properties of inertia in driving light and heavy particles advected by the rotating flow assessing the relative importance of the Centrifugal, Coriolis, added mass and Stokes forces and we show that rotation is extremely efficient in separating heavy from light particles, defeating the mixing properties of the underlying turbulent flow.
Eulerian fields. Rotation causes the generation of inertial waves in the flow [1].
Waves, and the associated instabilities, are of general interest given their fundamental character in atmospheric and oceanographic applications. The interplay between inertial waves and the two-dimensional three-components (2D3C) turbulent structures which develop in rotating turbulent flows is the subject of an active debate. Several authors [16,[25][26][27][28][29] have discussed the possibility to describe the dynamics of rapidly rotating 3D flows (limit of Rossby number much smaller than 1), in terms of wave turbulence triggered by triadic resonant interactions (for reviews on wave turbulence see, e.g., [30][31][32]). At the same time, experimental [13,33] and numerical studies [20,21] indicate that 2D turbulence provides an effective description of many aspects of rotating flows (for a recent review on 2D turbulence see [34]).
Theoretical studies [26,27], addressing inertial wave turbulence theory with a complete numerical solution in addition to the results of quasi-normal closures, and numerical simulations [16] have shown that the nonlinear wave interactions tend to concentrate energy in the wave-plane normal to the rotation axis, favoring the transfer of energy from the 3D fast modes toward the 2D slow manifold (see also [25] for a generalized quasi-normal approach, not restricted to the asymptotic limit and with quantitative comparisons to direct numerical simulation data). This has been proposed as a mechanism which creates the columnar vortices [23]. In particular, triadic wave interactions are able to capture the main part of the so called "spectral buffer layer", i.e., the spectral region close to the 2D slow manifold [27]. On the other hand, the leading resonant three-wave interactions cannot transfer energy directly to the 2D modes [16,20] and the wave approximation cannot be uniform as a function of the wavenumber. In other words, wave turbulence description ceases to be valid for very small wavenumbers in the direction of the rotation axis k || = k · Ω/Ω 0 and for very large wavenumbers in the perpendicular direction k ⊥ = k − k · Ω/Ω [29], see also [26] for a discussion about the decoupling of the 2D manifold. In such spectral regions the coupling of modes by near-resonant and non-resonant triads has been numerically investigated at moderate Rossby numbers by [36]. Previously, the decoupling of the 2D slow mode was questioned in [26], while a theoretical work [35] based on stability analysis of the 2D flow has shown the existence of a critical Rossby number, which depends on Reynolds, below which 3D rotating flow becomes exactly 2D in the long-time limit.
The scenario is complicated by the fact that the predictions obtained from the wave turbulence in infinite domains, with a continuous wavenumber space, could differ from the observations of numerical simulations and experiments, which necessarily deal with fluids confined in finite volumes. Note, in particular, that the exact decoupling of the 2D slow manifold from the inertial waves, due to resonant three-wave interactions, is not proven in the continuous case (see [26]). The discretization of the wavenumbers in finite volumes causes a gap between the 2D manifold and the 3D modes which could favor the decoupling of the 2D dynamics (see, e.g, [16] and references therein). The wave turbulence theory has been recently applied to the case of an infinite fluid layer confined between two solid boundaries [28]. In this case, the discretization of the k || allows to address the dynamics of the 2D manifold and its relationship with the wave-modes. In particular, it has been shown that the presence of a strong 2D mode might have a strong feedback on the waves dynamics, as inertial waves can be scattered by the vortices [28]. Along this line, recent experiments [37] and numerical simulations [38] have shown that a significant fraction of the kinetic energy is concentrated in the inertial waves whose period is shorter than the turnover time of the 2D structures, while waves with longer period are scrambled by the turbulent advection. Finally, recent numerical investigation of the rotating Taylor-Green flow [22]  Lagrangian particles. Lagrangian dynamics in rotating flows is at the core of many different physical and engineering problems, ranging from the dispersion and diffusion of pol-lutants, living species or mixing of chemical reagents, to cite just a few examples. However, the bulk of knowledge collected about the Eulerian properties of the flow has no counterpart in the Lagrangian framework. In the last decade a significant advance in the understanding of the dynamics of inertial particles suspended in turbulent flows has been achieved, notably for homogeneous and isotropic flows [39]. In the specific framework of particle dynamics in rotating flows, very few results are available. We mention a theoretical prediction for the spatial distribution of small, heavy particles in rotating turbulence [40], a prediction for the dispersion of fluid tracers in rotating turbulence [41], and the experimental study of the tracer-like particles acceleration statistics [42].
In this paper, we present the first attempt to assess the importance of Coriolis and centrifugal forces on the dynamical evolution and spatial dispersion of light and heavy particles, within the point-particles approximation. We show that the combined effect of inertia plus rotation leads to a singular behavior for the particles' statistics. In particular, the preferential sampling of high/low vorticity regions is strongly enhanced and characterized by anisotropic contributions on opposite directions: light particles tend to diffuse mainly vertically (i.e. along the rotation axis) while heavy particles are strongly confined in horizontal planes (see Fig. 2). As a result, the relative importance of Coriolis, centrifugal or addedmass forces might vary of order of magnitudes comparing light or heavy families. We suggest that, at any rotation rate of practical interest, both the 2D and the 3D turbulent structures are coupled together and that any attempt to separate them into a weak wave turbulence coupled with a quasi 2D slow-dynamics in the plane perpendicular to the rotation axis might fail to capture key properties for both Eulerian and Lagrangian statistics. This is an important remark for the phenomenology of Eulerian and Lagrangian rotating turbulence and to further improve its modelisation.
The paper is organized as follows. In section (II) we discuss the numerical set up concerning both Eulerian and Lagrangian properties. In section (III) we discuss the Eulerian statistical properties at changing both Rossby and Reynolds numbers, while in section (IV) we present the main results concerning the dispersion of light/heavy particles. Conclusions follow in section (V).

A. The equation of motion for the Eulerian flow and for the Lagrangian trajectories
The dynamics of an incompressible velocity field u in a rotating reference frame with angular frequency Ω is given by the three dimensional Navier-Stokes equations (NSE): Here ρ f and ν are the density and the kinematic viscosity of the fluid, respectively; 2Ω × u is the Coriolis force, and f is an external force. For an incompressible fluid, rotation breaks the statistical isotropy of the flow, but not its homogeneity. Note that the centrifugal force is the Zeman wavenumber [9,43,44] defined as the Fourier scale where the inertial turnover time, τ nl (k) = ε −1/3 k −2/3 becomes of the same order of τ Ω , i.e., k Ω ∼ (Ω 3 /ε) 1/2 , ε being the energy transfer rate. For Ro ≤ 1, the dynamics of the energy transfer will be largely influenced by rotation. Importantly enough, as soon as the Zeman wavenumber is larger than k f , an inverse energy transfer develops for k ≤ k f , characterized by a strong accumulation of the kinetic energy into 2D large-scale structures. As a result, for Ro ≤ 1, the system develops a forward cascade of energy, partially affected by the presence of rotation, and a simultaneous inverse energy cascade leading to a strong anisotropy. The need to resolve both interval of scales is the major bottleneck for direct numerical simulations.
In the reference frame rotating with angular frequency Ω, the equations for the trajectory r t and the velocity v(r t , t) of a small sphere of radius R and density ρ p suspended in the   fluid field u can be approximated as [45]: where r 0 is the position of the rotation axis. Within the point-particle model,the inertial dynamics is controlled by two non-dimensional parameters, the density ratio, β = 3ρ f /(ρ f + 2ρ p ), and the Stokes number, St = τ p /τ η , defined as the ratio between the particle relaxation time, τ p = R 2 /3βν, and the Kolmogorov time, τ η . The first term on the right hand side (rhs) of (3) is the fluid acceleration and results from an estimate of the added-mass and pressure gradients along the trajectory of the tracers. The second term is the Stokes drag. With respect to the case of homogeneous and isotropic flows two new forces appear in the rhs: the Coriolis and the centrifugal/centripetal, the third and fourth terms respectively.   To our knowledge, this is the first attempt to assess the effects of these two forces on the statistical and dynamical properties of inertial particles in turbulence. At variance with the NSE for the flow, the centrifugal force is present in the equation for the particle motion and it explicitly breaks homogeneity because of its dependency on the distance from the rotation axis. Its sign depends on the factor (β − 1): for heavy particles (0 ≤ β < 1) the force is centrifugal, while for light particles (1 < β ≤ 3) it is centripetal.
In equation (3), we have neglected the Basset history and gravity forces, and the Faxen corrections; moreover, we have approximated the material derivative along the inertial particle trajectories in terms of the material derivative along tracer paths. In the previous set up, tracer trajectories are evolved according to the equation:ṙ t = u(r t , t) .

B. Direct Numerical Simulations set-up
We performed a set of state-of-the-art high-resolution direct numerical simulations of the NSE in a periodic, cubic domain of size L = 2π with up to N 3 = 4096 3 collocation points. The rotation axis is in the x-direction, i.e, Ω = (Ω, 0, 0). The integration of Eqs.
(1) has been performed by means of a fully dealiased pseudo-spectral code, with second order Adams-Bashforth scheme with viscous term exactly integrated. The parameters of the Eulerian dynamics for the different runs are reported in Table I. The integration of eqs.
(2) is performed by interpolating the Eulerian velocity field and its derivatives with a 6 − th order B-spline algorithm on the particle position [46]. Parameters of Lagrangian dynamics can be found in Table II A. At high rotation, the presence of a simultaneous forward and inverse cascade asks for an optimized set-up, to minimize spurious finite-size effects. The critical Rossby number where energy starts to flow upscale is known/believed to depend on the way the system is forced [21,22,47] and on the aspect ratio of the volume where the flow is confined [16,28,48]. In presence of an inverse flux, it is crucial to force the system at intermediate wavenumbers to allow the large-scale flow to develop its own dynamics, without being directly influenced by the forcing. Moreover, an energy sink mechanism must be added to prevent the formation of a condensate at the lowest Fourier mode that could spoil the statistics at all scales.
To match the previous requirements, we adopted a stochastic isotropic Gaussian force, f , active on a narrow shell of wavenumbers at k f ∈ [4 : 6]. The amplitude vector of each forcing Fourier mode is obtained as f (k, t) = f 0 (i k × X(t)). The variables X i (t) are independent, identically distributed time-differentiable stochastic processes, solution of the following Ornstein-Uhlenbeck 2-nd order process: In the expression above, τ f is the correlation time of the process, and W i (t) is a Wiener process. It is important to stress that the above time-correlated process ensures the continuity of the Lagrangian acceleration of the tracers (see [49] for details). At low Rossby, to arrest the inverse cascade, we remove energy at large scales with a linear friction term, α∆ −1 u and acting on wavenumbers |k| ≤ 2 only. This term is added to evolve the Lagrangian particles on a stationary state without the need to over resolve the field in the infrared regime.
To understand the basic phenomenology of a rotating turbulent flow, it is useful to recall the different dynamical states that can be observed in the case of strong rotation, i.e. low Rossby number. In Figure (3), we show the temporal evolution of the total kinetic energy, E kin = dp |u p | 2 , starting from a fluid at rest and until a stationary regime is achieved.
In the early stage, the rotation rate Ω is zero, and the flow develops a 3D direct cascade.
Small-scale thermalization is indicated by the overshoot of the kinetic energy. This is the standard situation of stationary, non-rotating turbulent flows where the energy input is balanced by viscous dissipation. After this stage, we switch on the rotation and the inverse energy cascade starts to develop if Ω is large enough: this is indicated by the linear growth in time of the kinetic energy. Later, we switch on the damping term at large scale. Doing that, we end up with a statistically stationary regime for a strongly rotating turbulent flow.
In the inset of the same figure we show the presence of a simultaneous positive and negative spectral flux when the rotation rate is large enough, indicating the existence of a forward and inverse energy transfer for scales smaller and larger of the forcing scale, respectively. We remark that the spectral flux is here defined as the transfer of energy across a wavenumber k by the non-linear interactions N p of the Navier-Stokes equations [50]: Π(k) = |p|<k dp u * p · N p . The simultaneous presence of direct and inverse cascades is shown by the two plateaux in the spectral flux, in agreement with previous findings [16,20,48,[51][52][53][54].
Once the turbulent flow is stationary, we seed it with Lagrangian particles of different inertia, released with the same velocity of the underlying fluid. When Ro is small, the flow is characterized by the presence of few intense columnar cyclones, co-rotating with Ω. In the plane perpendicular to the rotation, the associated two-dimensional vortices are much slower than any other structure in the flow. Moreover, as shown in Figure (2), these cyclones strongly influence the distribution of the particles. Light particles are trapped inside, while heavy particles are ejected, leading to an extreme, singular preferential sampling of the underlying flow (see sec. IV). In all cases here investigated, the flow displays a few big cyclones well separated from each other. The breaking of the cyclone-anticyclone symmetry is a well known feature of rotating turbulence [22,[55][56][57][58][59]. It is also the indication that the formation of the vortical columnar structures cannot be entirely due to a 2D inverse cascade regime, because in the plane perpendicular to the rotation axis the symmetry is not broken.
Nevertheless, it is suggestive to interpret the presence of three long-living coherent columns in terms of the dynamics of point vortices, since a system of three equal-sign point-vortices is linearly stable in two dimensions [60]. We cannot exclude that the columnar vortices would eventually merge into a single cyclone, after long enough time. Considering that vortices with equal sign repel each other, and that their merging would cause the generation of an intense shear between them, it is arguable that the process of a collapse is unlikely to occur.

A. Fourier analysis
Rotation affects the spectral distribution of the kinetic energy on a wide range of scales.
For Fourier modes between the forcing and the Zeman wavenumbers, k f < k < k Ω , a standard phenomenological argument predicts for the energy spectrum: where the distance r is in the plane normal to Ω, and the versort is orthogonal to both Ω and r. Thus, we define the p-th order transverse structure function (TSF) as: where isotropy is assumed in the normal plane. In Figure (5 figure), rotation effects are always sub-leading. Here the statistics is in good agreement with the Kolmogorov K41 prediction [50]. At Ro << 1, the scaling laws are always spoiled by anisotropy and the only systematic way to disentangle scaling properties would be to resort to a decomposition in terms of eigenfunctions of the group of rotations [67,68]. Moreover, the flow is naturally bimodal, with a 2D3C dynamics superposed and entangled with the 3D turbulent fluctuations.
To better clarify the statistics scale-by-scale, we propose to decompose the velocity field into two components, one given by the 2D3C slow modes and the other associated to the 3D fast modes, u(x, y, z|t) = u 2D (y, z|t) + u (x, y, z|t) .
Here we have defined the two-dimensional field as the average of the velocity field in the direction of Ω: u 2D (y, z|t) = dx u(x, y, z|t). In Figure (6), we plot the second order transverse structure function, S ⊥ (r) measured for the undecomposed field, and for the two fields obtained by the above decomposition. The figure shows the existence of a scale, of the order of l Ω = 2π/k Ω , where the statistics changes from being 2D3C to 3D dominated. The background field follows quite closely the Kolmogorov scaling ∝ r 2/3 (not shown), while the 2D3C field has a scaling much smoother than the Zeman estimate. This does not necessarily contradicts the results of section (III A). Rather, it clearly shows that within the Eulerian statistics there are two different components that influence the physics at different scales.
It also suggests that any attempt to fit/predict scaling laws without a separation of the different contributions might lead to uncontrolled approximations.
In Figure (7), we plot the 4th-order (left panel) and 6th-order (right panel) flatness derived from the transverse structure functions, 3 for the undecomposed velocity field, the 2D projection, and the fluctuating part. Except for very large spatial increments, the curves are always far from the Gaussian limit. Consider the data for the whole field at large rotation rate, i.e., Ω = 10 (empty squares in both panels). The 4th-order flatness display a weak dependence on the analyzed scale in the  [33,65,69,70] is merely apparent and probably due to a non-trivial combination of effects induced by the coherent structures and contributions from the underlying 3D turbulent fluctuations. This is one of the main results of this paper.
Finally, let us notice that naively one would expect that the flatness of the fluctuating field for Ω = 10 should be equal to the flatness of the total field for Ω = 4 (filled squares). Our data show that this is not the case, meaning that rotation not only leads to the formation of the 2D columnar structures, but also modifies the 3D fluctuating turbulence, if the Rossby is small enough. We summarize in table III the results for the best fit to the flatness scaling exponents for the fluctuating components at high rotation and for the total component at small rotation rates.
A novel way to investigate the statistics of the columnar vortices is to exploit the peculiar features of the inertial particles in sampling the flow. It is known that light particles are attracted inside the vortices, while heavy particles are expelled out of them [39,71]. By studying the velocity statistics measured along the trajectories of particles with different inertia, it is possible to use their preferential concentration in specific flow regions, to enhance or deplete the contribution of the slow vortical modes with respect to the turbulent background.
Here, we start by analyzing the different contributions of the forces that influence inertial particles motion. In Figure ( When the Rossby number is small, i.e. for Ω = 10, the inertial particle dynamics does not always attain a statistically steady state. The relative importance of the forces is affected by two different reasons. The first one is purely kinematic, since both Coriolis and centripetal forces are proportional to the rotation rate. The second one is dynamical: the organization of the flow, with the formation of strong columnar vortices, competes with the kinematic effects. If particles are heavier than the fluid, the centrifugal force soon becomes dominant: particles not only tend to avoid coherent vortical structures, but also tend to spiral away from their rotation axis very efficiently (see also Fig. 2). This enhanced centrifugal action is balanced by the Stokes drag only. Comparing, a St rms and a Cp rms , it is interesting to note that this balance is very efficient, leading to a total acceleration, a tot rms , much smaller than the single contributions, and to an almost stationary statistics in the long time limit. In this regime, the dynamics of the heavy particles is uncorrelated with respect of the underlying fluid. Particles move away from their rotation axis with a spiral motion, whose radius grows exponentially in time, r(t) ∼ exp(Ω 2 τ p t). Since their velocity also increases exponentially over time, the particle Reynolds number might eventually become too large for the validity of the model equations (2) and (3), [45]. Hence, it would be crucial to perform a systematic comparison with experimental data, in order to understand the limitation of the point-like approach in the limit of very heavy particles.
For small Rossby number, we observe an opposite behavior in the case of light particles. The centripetal force attracts the light particles toward their original rotation axis, but its intensity vanishes as r t → r 0 . The overall effect is therefore to spatially confine the trajectories of light particles, depleting turbulent diffusion and preventing them from exploring regions far away from the rotation axis. Additionally, one needs to consider the dynamical attraction inside the coherent vortical structures. As visually shown in Fig. (2), preferential centripetal concentration is the leading effect and almost all light particles are trapped inside vortical structures. The leading term in light particles acceleration is the added mass, which is not balanced by other forces. We also notice that at long times the temporal behavior of the added mass term becomes noisy, in spite of the large number of particles used in computing the average. This is because eventually all light particles collapse into the cores of a few columnar vortices, thus reducing the effective statistics.
The singular role played by the presence of the coherent structures for Lagrangian statistics is better quantified in Fig. (9), where we plot the preferential sampling of specific flow regions, by measuring the average vertical vorticity at the particle positions normalized with the averaged vertical vorticity in the volume, At low rotation, Ω = 4, the preferential sampling by heavy or light particles is similar to what observed in homogeneous and isotropic turbulence, and quantitatively it is a O(1) effect with respect to the mean fluid vorticity. At large rotation rate, the situation is different: heavy particles, because of the sweeping due to the centrifugal forces, do not show preferential concentration, while light particles over-sample the intense vorticity regions with an effect which is a factor O(100) larger. In Fig. (9), we also show the Probability Distribution Function (PDF) of the vertical vorticity w x along particle trajectories for tracers and for one light and one heavy family. Notice the bimodal PDF for the light particles induced by the trapping in the vortex cores; for the heavy particles, the PDF is symmetric, because of the homogeneous sampling of the flow regions outside strong vortical structures.
Concerning absolute dispersion, the influence of the strong vortical structures will induce a systematic anisotropic effects for tracers [41]. On the other hand, since the vortical structures are fatal traps for the light particles and strong repellers for heavy particles, we expect to measure strong deviations in the single particle dispersion too. In Figure (10), we show the mean square absolute dispersion of the particles from their initial position as a function of along different directions i = (x, y, z), here normalized with the ones measured for tracers.
For the heavy particles, we find that the diffusion in the plane normal to the rotation axis Ω is enhanced, due to the centrifugal effect; while parallel diffusion is reduced. Moreover at fixed value of the density mismatch β, the effect is stronger for higher Stokes number.

V. CONCLUSIONS
Rotating turbulence is key for many industrial and geophysical applications. In many empirical set-up it is also key to control the dispersion and advection of particles.  [72] for a discussion of discreteness and resolution effects). Finite volume effects can be estimated by comparing the typical distance traveled by the waves during the duration of the simulation, estimated in terms of their group velocity. In our case, for the highest rotation case this distance is pretty small, of the order of 10% of the total volume. Using state-of-the-art highly resolved DNS, as done here, is crucial in order to reduce the spectral gap with the horizontal plane also. Here, for the highest resolved case we have an excellent resolution of the buffer layer near k || = 0, including wavenumbers with angle close to 0.04 degrees with the horizontal plane and thus reducing finite volume effects. However, improving angular resolution to better capture the spectral buffer layer also at small wavenumbers close the two-dimensional manifold is an issue [26]: further numerical investigations e.g., in slab geometries, permitting to obtain values of k || = 0 small enough are desirable to shed further light on the problem of 2D-3D modes coupling.
Other numerical approaches meant to understand the importance of different triadic interactions in Fourier-space, and to further clarify the nature of the inverse cascade in purely rotating turbulence, are possible. One important example is given by [73] where reduced Navier-Stokes equations including only near-resonant, non-resonant and near twodimensional triad interactions are considered. These numerical approaches are restricted to work on spectral space, with severe limitation in the number of modes that can be considered.  u(x, y, z) (filled squares). We also superpose the power law prediction ∝ r −0.15 obtained from independent measurements of isotropic turbulence without rotation and the best fit for the power law measured with intense rotation in the present data ∝ r −0.35 . Right: the same for the 6th-order flatness, K