Pattern formation by turbulent cascades – Nature

Direct numerical simulations of the Navier–Stokes equation with odd viscosity

Direct numerical simulations of the Navier–Stokes equation with odd viscosity (equation (1)) are performed in a cubic box of size L = 2π with periodic boundary conditions (such that the smallest wavenumber is 2π/L = 1). Our results can be reproduced with any Navier–Stokes solver by including a modified Coriolis term modulated by k² (or, equivalently, by a vector Laplacian for real-space-based methods) to account for odd viscosity. We use a pseudo-spectral method with Adams–Bashforth time-stepping and a 2/3-dealiasing rule⁵⁸. Both normal and odd viscosities are integrated exactly using integrating factors. The forcing f(t, k) acts on a band of wavenumbers k ∈ [k_in, k_in + 1] with random phases that are delta-correlated in space and time, ensuring a constant average energy injection rate ϵ = ⟨u ⋅ f⟩. It has a zero mean component ⟨f(t, k)⟩ = 0 and covariance \(\langle {\bf{f}}(t,{\bf{k}})\cdot {\bf{f}}({t}^{{\prime} },{{\bf{k}}}^{{\prime} })\rangle ={\epsilon }\delta (t-{t}^{{\prime} })\delta ({\bf{k}}-{{\bf{k}}}^{{\prime} })\). The time-step is chosen to resolve the fastest odd wave with frequency \({\tau }_{{\rm{odd}},\max }^{-1}={\nu }_{{\rm{odd}}}{k}_{\max }^{2}\), where \({k}_{\max }\) is the highest resolved wavenumber in the domain. We find that stable integration requires a time-step \(\Delta t\lesssim 0.1{\tau }_{{\rm{odd}},\max }\). A complete overview of the input parameters for the simulations in this work is provided in the Supplementary Information. Approximately 3 million CPU hours were required to perform the simulations underlying this work.

Effect of odd waves on the nonlinear energy transfer

In this section, we describe how the waves induced by odd viscosity (odd waves) affect the nonlinear energy transfer. Our analysis closely follows that of rotating turbulence^2,4,36,59.

Nonlinear energy transfer

Fourier-transforming the Navier–Stokes equation, multiplying with v^*(t, k) (where the asterisk denotes complex conjugation), and adding the complex conjugate, we find the energy balance equation^2,4

where ν = η/ρ is the kinematic viscosity, and in which

This term describes the nonlinear energy transfer between scales, whereas F = v^* ⋅ f corresponds to energy injection by the forcing term f. The term −2νk²E represents standard viscous dissipation. In equation (7), the sum runs on momenta p and q such that k + p + q = 0, and P_ij(k) = δ_ij − k_ik_j/k² is the projector on incompressible flows.

At first glance, equation (6) is left unchanged by odd viscosity, because of its non-dissipative nature. However, odd viscosity has indirect effects on the energy transfer (in the same way as the non-dissipative Coriolis force has an indirect effect on the energy transfer in rotating turbulence).

Odd waves

To see that, we first consider the linear and inviscid limit of the Navier–Stokes equation (1) (so we set ν = 0 and (u ⋅ ∇)u = 0). As detailed in the section ‘Linear stability of the fluid and odd waves’ of the Supplementary Information (in which we consider a more general odd viscosity tensor), this equation has wave solutions of the form

in which h^±(k) = e(k) × (k/k) ± ie(k) with \({\bf{e}}({\bf{k}})={\widehat{{\bf{e}}}}_{z}\times {\bf{k}}/\parallel {\widehat{{\bf{e}}}}_{z}\times {\bf{k}}\parallel \) (ref. ⁶⁰) with frequency

Taking into account normal viscosity leads to an additional exponential decay of the waves with the rate −νk² (Supplementary Information). In particular, we note that the linearized Navier–Stokes equation does not exhibit any linear instability. By construction, k ⋅ h^±(k) = 0, so these modes represent incompressible flows. Furthermore, \({({{\bf{h}}}^{+}({\bf{k}}))}^{* }\cdot {{\bf{h}}}^{-}({\bf{k}})=0\) and \({({{\bf{h}}}^{\pm }({\bf{k}}))}^{* }\cdot {{\bf{h}}}^{\pm }({\bf{k}})=2\). Hence, odd waves provide an orthonormal basis for incompressible flows. As k × h^± = −kh^±, the basis functions have a well-defined helicity ∓1.

Decomposition of the energy transfer on odd waves

Expanding the velocity field as a superposition of helical waves

in which \({v}_{s}^{* }(t,{\bf{k}})={v}_{s}(t,-{\bf{k}})\) to ensure the reality of v(t, x), the Navier–Stokes equation becomes

in which we have used the short \({v}_{{s}_{k}}\) for \({v}_{{s}_{k}}(t,{\bf{k}})\), the term \({f}_{{s}_{k}}({\bf{k}})\) corresponds to the forcing term, and

satisfy C_k∣p,q = C_k∣q,p.

Helicity and energy conservation of inviscid invariants

In terms of the components v_±(k), energy and helicity, respectively, read^4,60

A direct calculation shows that^36,60

and

from which we deduce that energy and helicity are conserved when normal viscosity and the forcing can be neglected (ν = 0 and f = 0), even if odd viscosities are present. In particular,

so using equation (11), we find (when ν = 0 and f = 0)

This equation shows that the nonlinear energy transfer T(k, t) in equation (7) is suppressed when averaged over long times compared to ω(k) + ω(p) + ω(q), unless this quantity vanishes exactly, as is the case for 2D modes (Fig. 2h, blue line, corresponding to modes with k_z = 0).

Resonant manifold

The 2D modes with k_z = 0 form a so-called slow manifold, or resonant manifold, that contributes to most of the nonlinear energy transfer. Furthermore, isolated triads with k_z ≠ 0 can also satisfy the resonance condition ω(p) + ω(q) + ω(k) = 0. In the case of rotating turbulence, resonant triads primarily transfer energy from the 3D modes to the quasi-2D slow manifold with k_z = 0, leading to an accumulation of energy in the slow manifold, enhancing the two-dimensionalization of the flow^{4,33,34,35,36}. We expect a similar phenomenon to occur in the case of fluids with odd viscosity owing to its similarity to rotating fluids, as is also suggested by the two-dimensionalization observed in our numerical simulations. As a consequence, the effective spatial dimension of the system depends on the scale at which it is observed (such as in rotating turbulence or thick layers^4,35,61). More insights may be obtained by developing a weak turbulence theory for odd waves, in the same spirit as for rotating flows (we refer to refs. ^26,27,28 for more details on wave turbulence).

Scaling relations and wavelength selection

Scaling relation for the energy spectrum

We first analyse the power spectrum, building on the phenomenological theory of ref. ³⁹ (see ref. ⁴² for a review). This theory relies on the following hypotheses: (1) energy is conserved away from injection and dissipative scales; (2) the cascade is local, which means that different length scales are coupled only locally (for example, very large scales are not directly coupled to very small scales); and (3) the rate of energy transfer ε(k) from scales higher than k to scales smaller than k is directly proportional to the triad correlation time τ₃. Because of hypotheses 1 and 2, the rate of energy transfer ε(k) is constant across the scales (that is, does not depend on k) and can be identified with the energy dissipation rate ϵ. Moreover, because of hypothesis 2, ε(k) should depend only on local quantities k and E(k), in addition to τ₃(k). Therefore, using hypothesis 3, we write

where A is a constant. The exponents are found using dimensional analysis (with [E] = L³T⁻², [ϵ] = L²T⁻³, [k] = L⁻¹, [τ₃] = T), which yields α = 4 and β = 2.

In Fig. 2, we argue that the eddy turnover time τ_E is the relevant timescale when k ≪ k_odd (so we set τ₃ = τ_E in the expression above), whereas the frequency of odd waves ω is the relevant timescale when k ≫ k_odd (so we set τ₃ = ω⁻¹). The dispersion relation of odd waves is computed in the Supplementary Information and given in equation (9). The eddy turnover time is τ_E(k) = [1/k]/v_k. As E(k) is the shell-average of v_k, we have dimensionally \(E(k)\propto {v}_{k}^{2}/k\), so v_k = [kE(k)]^1/2. Putting everything together, we end up with equations (3) and (4) of the main text.

Scaling relation for k
_c

For the condensation of the forward energy flux, the collapse of numerical results indicates that it can be described by a master scaling law

Using the Kolmogorov spectrum for the case without odd viscosity, E₀(k) ∝ ϵ^2/3k^−5/3, we find for the energy spectrum

Using the scaling argument of the previous section (see equations (3) and (4)), we find s = 2/3, which is compatible with the numerical results. This scaling continues until dissipation saturates the condensation. We can thus estimate the location of the condensation peak k_c from the balance between injection and dissipation. Neglecting contributions to the dissipation from wavenumbers k < k_odd (where there is no meaningful change from Kolmogorov scaling), we obtain

Assuming k_odd ≪ k_c, this yields

resulting in the scaling relation for the peak condensation

where in the last relation, we substituted the normal Kolmogorov wavenumber k_ν ∝ ϵ^1/4ν^−3/4.

For s = 2/3, we find

as quoted in the main text.

Estimation of the height for the peak

The mechanism of non-dissipative arrest analysed in this work is reminiscent of but distinct from the bottleneck effect^{62,63,64,65,66,67} generated by the usual viscosity.

A coarse estimate of the height of the peak in E(k)/E₀(k) can be obtained by evaluating equation (4) (to get E(k)) and equation (3) (to get E₀(k)) at k = k_c given by equation (5), yielding \(h\equiv E({k}_{{\rm{c}}})/{E}_{0}({k}_{{\rm{c}}})\,\propto \)\({({\nu }_{{\rm{odd}}}/\nu )}^{1/3}\) (see Extended Data Fig. 1d for a comparison with numerical data). Notably, this suggests that h depends on only the ratio of odd to normal viscosity. We also note that h increases as normal viscosity ν decreases (that is, when the Reynolds number increases), in contrast with the bottleneck effect due to dissipative viscosity^{62,63,64,65,66,67} in which the magnitude of the effect decreases as viscosity decreases.

Wavelength selection

In Extended Data Fig. 1c, we plot an estimate of the power spectrum of the vorticity, evidencing wavelength selection in the vorticity. This suggests that the characteristic wavelength 2π/k_c should be directly visible in snapshots of the vorticity field. This can be seen in Fig. 3b. The width of the peak leads to a wide distribution of structure sizes in the image.

We expect the wavelength selection mechanism due to the arrested cascade to persist at arbitrarily long times and to resist small perturbations, in contrast with metastable patterns arising from kinetic effects⁶⁸ in which the system resides in metastable states for long but finite periods (see Supplementary Information for convergence plots).

The wavelength selection mechanism we have described can be compared with that in active turbulence, for instance in bacterial suspensions and self-propelled colloids^{69,70,71,72,73,74}. In active turbulence, however, it has been reported that there is no energy transfer across scales (and hence no cascade): energy is typically dissipated at the same scale as it is injected, and it is believed that the wavelength selection is the result of a scale-by-scale balance (see, for instance, Figs. 3d and 4g and sections 3.2.2 and 4.2.3 in ref. ⁷¹ and references therein). We note, however, that finite energy fluxes have been reported in certain cases^{56,75,76,77,78,79,80}.

In these systems, wavelength selection has been described as the result of a Swift–Hohenberg-type term included in the stress tensor (leading to a finite-wavelength linear instability), to which noise is added⁷¹. By contrast, cascade-induced pattern formation cannot be directly traced to a linear instability of Navier–Stokes equation (1) (see section ‘Effect of odd waves on the nonlinear energy transfer’ as well as Supplementary Information section ‘Linear stability of the fluid and odd waves’ for a linear stability analysis). The linear stability analysis does not predict any instability, neither to a stable branch with a particular wavelength nor to an unstable branch that could itself bifurcate to the state of interest as part of a subcritical bifurcation.

An analogy with similar situations such as Rossby and drift wave turbulence^{81,82,83,84,85,86} and laminar and turbulent patterns in wall-bounded shear flows^87,88,89,90 suggests that the wavelength selection may be described by considering the linear stability of the statistically averaged Navier–Stokes equation, for instance, using an appropriate turbulence closure model.

Anisotropic energy spectra

In line with the inherent symmetry of the system, we now consider cylindrically averaged energy spectra E(k_⊥, k_z), which distinguish the horizontal (perpendicular) directions from the vertical direction^91,92,93,94. To reveal in which part of the k-space the energetic condensation occurs, we compute the cylindrically averaged spectrum of the cases with odd viscosity normalized by the spectrum of the reference case without odd viscosity (Extended Data Fig. 1b). Starting with the direct cascading case in Extended Data Fig. 1b (top panel), we see that indeed the flow remains mostly 3D isotropic for k < k_odd and then proceeds to condensate anisotropically into the low-k_z manifold because of the quasi-2-dimensionalization effect of the odd viscosity. As detailed in the main text, the condensation is saturated by dissipation, leading to a peak condensation wavelength k_c, which is thus primarily visible in the perpendicular directions because of the anisotropic condensation. The dominant vertical scale hence remains closer to k_odd. This leads to a crude estimate for the aspect ratio γ of the features in the pattern produced by the odd viscosity as

For the case presented in Fig. 3b, this leads to an aspect ratio γ ≃ 3.

For the inverse cascading case (Extended Data Fig. 1b, bottom panel), we again observe anisotropic condensation in the region k > k_odd. In the region k < k_odd, however, the kinetic energy for the case with odd viscosity is larger than the case without odd viscosity, as indicated in dark orange. This is because in this range, we expect the same diffusive equipartitioned scaling E(k) ∝ k² for both cases with and without odd viscosity, and there is no active dissipative mechanism to deplete the excess energy that has accumulated at higher wavenumbers in the case with odd viscosity.

Experimental considerations

In this section, we discuss the conditions required to observe the wavelength selection described in the main text in a fluid with odd viscosity. In short, we expect this effect to occur, for instance, in a fluid of self-spinning particles large enough to be inertial (not overdamped).

First, the Reynolds number Re = UL/ν has to be large enough. This puts constraints on the viscosity ν of the fluid, the details of which depend on the experimental setup considered. The current experimental systems we are aware of in which explicit measurements of odd viscosities were reported (active spinning colloids³², magnetized graphene¹⁰ and magnetized polyatomic gases^31,95) are all in a regime in which the nonlinear advective term in the Navier–Stokes equation can be neglected, either because ν is large enough or for geometric reasons; effectively, Re ≪ 1. Note also that experimental instances of (especially 2D) odd fluids may include a substrate, on top of which the active particles move. This can lead to the addition of an effective drag force −γv in the Navier–Stokes equation describing the odd fluid made of these particles. If such a term is large, it would prevent the existence of an inertial regime, and probably spoil the phenomenology discussed here.

Second, the ratio ν^odd/ν has to be large enough for the effect to be visible. When ν^odd ≲ ν, energy is dissipated as soon as, or before any effect of odd waves can arise. Henceforth, observing the effects of odd waves on turbulence would require ν^odd > ν. Odd viscosities (ν^odd ≠ 0) typically arise in systems breaking time-reversal and inversion symmetry at the microscopic scale^9,96,97. They have been experimentally measured in polyatomic gases under magnetic fields^31,95, spinning colloids³² and magnetized electron fluids¹⁰. They have also been predicted in systems, including fluids under rotation⁹⁸, magnetized plasma^11,99,100, quantum fluids^{96,101,102,103}, vortex matter¹⁰⁴, sheared granular gases¹⁰⁵, assemblies of spinning objects^{69,106,107,108,109,110,111,112,113,114,115,116,117} and circle swimming bacteria^118,119. In the systems mentioned above, in which experimental measurements of odd viscosity have been reported, ν_odd/ν reaches at most 1/3 (in active spinning colloids³² and magnetized graphene)¹⁰. From a theoretical point of view, the ratio ν_odd/ν is expected to increase linearly with the time-reversal breaking field. For instance, ideal vortex fluids are predicted to have a finite ν_odd but a vanishing ν (ref. ¹⁰⁴), leading to an infinite value of ν_odd/ν. Kinetic theory calculations for magnetized plasma (ref. ⁹⁹, section 19.44) predict ν = ν₀/[1 + x²] and ν_odd = ν₀x/[1 + x²] in which x = 2ωτ with τ is a collision time and ω ∝ B is a frequency proportional to the magnetic field B, whereas ν₀ is the value of normal shear viscosity when B = 0. Similarly, kinetic theory in rotating gases leads to an identical result in which x ∝ Ω is proportional to the rotation speed⁹⁸. In electron gases in graphene, experiments have been performed that validate these theoretical calculations¹⁰ (with x = B/B₀, where B₀ is a reference magnetic field). This results in a ratio ν_odd/ν = x ∝ B. Likewise, in active fluids, theoretical works suggest that ν_odd is proportional to the rotation speed of the spinning particles¹⁰⁶.

Rossby and drift wave turbulence

Extended Data Fig. 2 shows examples of simulations of the Rossby and drift wave turbulence mentioned in Fig. 4b. A brief review is contained in the Supplementary Information, and we refer the reader to refs. ^{13,14,15,16,17,18,19,93,120,121,122,123,124,125,126,127} for more details. In the figure, we simulate the Charney–Hasegawa–Mima (CHM) equation^{93,120,121,122,123,124,125,126}

in which J(a, b) = (∂_xa)(∂_yb) − (∂_ya)(∂_xb), ω = Δψ and ψ is the stream function, defined such that the velocity field is v = −ϵ ⋅ ∇ψ (ϵ is the 2D Levi-Civita symbol). The parameter β represents the gradient of the Coriolis force in a β-plane approximation; α represents large-scale friction and ν is viscosity, whereas f_ω is a vorticity forcing. Simulations are performed using the open-source pseudo-spectral solver Dedalus¹²⁸.

Note that in Rossby wave turbulence, the only exact inviscid invariants are energy and helicity. However, it has been established that a quantity dubbed zonostrophy evolves slowly enough to be considered as an invariant for practical purposes^{13,93,129,130,131}. This raises the question of whether such an adiabatic invariant may exist for odd turbulence, and whether it can predict the direction of the cascades (see refs. ^4,132 for discussions of the relation between inviscid invariants and the direction of turbulent cascades).

Minimal model of mass cascade with scale selection

In this section, we consider a simple model of the mass cascade that exhibits wavelength selection.

Mass cascades can, for instance, occur in the pulverization of objects into debris or in the coalescence and breakup of droplets^{23,133,134,135,136,137,138}. These processes can be modelled by the aggregation and fragmentation of clusters composed of monomers linked together: two clusters that collide may merge into a larger cluster; and a given cluster may split into smaller ones, spontaneously or on collision. The mean-field kinetics of these processes is described by a population balance equation generalizing the so-called Smoluchowski equation^{23,139,140,141} that can exhibit scale-invariant cascades, similar to that present in the Navier–Stokes equation^142,143,144. This kinetic equation may describe two classes of situations: (1) closed systems in which mass is conserved and (2) open systems in which particles are injected and removed from the system. Case 1 may somehow be compared with freely decaying turbulence, whereas case 2 may be compared with driven turbulence in which energy is injected and dissipated.

We expect that the balance between aggregation and fragmentation will lead to a preferred size if large clusters tend to break up, whereas small clusters tend to coalesce. Such a preferred size should manifest as a peak in the distribution of aggregate sizes. Such a peak has been reported, for instance, in the case of raindrop sizes^24,145,146, in which the distribution originates from complex mechanisms, including air turbulence and fluid fragmentation^{147,148,149,150,151,152}.

In our toy model, we consider clusters M_n composed of 2ⁿ⁻¹ monomers M₁, with n = 1, …, N. This is reminiscent of what is done in shell models of turbulence⁵², in which the wavenumbers are chosen in geometric progression. We assume that (1) there are interactions only between clusters of the same size and (2) there is a maximum cluster size N. The first assumption ensures that the mass fluxes are local, and the second enables us to consider a finite number of equations. We include a constant source of monomers, and a sink that removes the largest clusters M_N. In the case of raindrops in a cloud, for instance, the source may describe the condensation of droplets from vapour, and the sink may describe the precipitation of large droplets out of the cloud. The model is summarized by the set of reactions

in which M_n (n = 1, …, N) represents a cluster of size 2ⁿ⁻¹ (M₁ represents a monomer), and J_in, J_out and \({k}_{n}^{\pm }\) are the rates of the corresponding reactions.

The number densities c_n of clusters then follow the dynamical equation

where

in which it is implied that c_n ≡ 0 for n < 1 and n > N.

We can also consider the mass density ρ_n = 2ⁿ⁻¹m₀c_n, in which m₀ is the mass of a monomer. Multiplying equation (29) with 2ⁿ⁻¹m₀, we find that the terms with prefactors \({k}_{n}^{\pm }\) cancel as in a telescoping series. This manifests that equation (29) with J_in = J_out = 0 conserves mass. It is therefore convenient to introduce the fluxes

and

corresponding to the reactions with rates \({k}_{n}^{\pm }\), and such that

To induce wavelength selection, we choose particular forms for \({k}_{n}^{\pm }\). The basic idea is the forward flux \({k}_{n}^{+}\) should decrease with n, whereas the backward flux \({k}_{n}^{-}\) should increase with n. Experimentation suggests that various strictly increasing functions of (N − n)/(N − 1) and (n − 1)/(N − 1), respectively, lead to similar results. We choose

Equation (29) is then solved starting from the initial condition c_n = 0 for all n using DifferentialEquations.jl (ref. ¹⁵³) with a fourth-order A-stable stiffly stable Rosenbrock method (Rodas4P) until a steady state is reached. The resulting steady state is shown in Extended Data Fig. 3. In Extended Data Fig. 3c, we observe that the density c_n is peaked at an intermediate value \({n}_{c}^{* }\) (pink dashed line), which is neither the maximum cluster size N, nor the monomer size 1, demonstrating wavelength selection. Similarly, Extended Data Fig. 3d shows that the mass density ρ_n is peaked around a (different) scale \({n}_{\rho }^{* }\) (red dashed line). As we have considered a mean-field description that does not take space into account, there is no proper pattern-only wavelength selection.

We observe in Extended Data Fig. 3e that the flux J_tot ≡ J⁺ + J⁻ (black curve in inset) is constant and nonzero for 1 < n < N. In 1D, the existence of a steady state is equivalent to a constant flux. (Note that certain models of aggregation–fragmentation may exhibit oscillations, that is, limit cycles instead of fixed points^154,155). The total flux can be decomposed into the forward flux J⁺ associated with reactions with rates \({k}_{n}^{+}\) and the backward flux J⁻ associated with reactions with rates \({k}_{n}^{-}\), respectively, defined in equations (31) and (32), and plotted in Extended Data Fig. 3e (red and blue curves, respectively).

In Extended Data Fig. 3h, we analyse the initial value problem obtained by setting J_in = J_out = 0 in equation (29). An exact solution of this model is given in the Supplementary Information. Wavelength selection may occur, although there is no net flux. This can be compared with the arrest of coarsening that can arise in mixtures and similar mass-conserving systems, even if the mass is not injected and removed from the system^{46,47,48,49,50,51,156,157,158}. We also observe that wavelength selection occurs only when the total number of monomers is large enough, which is reminiscent of what happens in so-called beam self-cleaning in optics, in which light in an optical waveguide at sufficiently high power may undergo a nonlinear redistribution of the mode powers that favours the fundamental, similar to an inverse cascade¹⁵⁹.

Equation (29) describes the mean-field dynamics of the reactions (28). To check whether the effect is still present beyond mean field, we solve the corresponding Doob–Gillespie kinetic Monte Carlo problem using the package Catalyst.jl (refs. ^153,160). The result of the simulation is shown in Extended Data Fig. 3g, and compared with mean-field simulations, with excellent agreement.

Finally, we discuss the rate of entropy production in the system. To do so, it is convenient to introduce the rates \({k}^{+,n}={k}_{n}^{+}/2\) and \({k}^{-,n}={k}_{n+1}^{-}\) to match the notations used in the literature on chemical reaction networks^{161,162,163,164}. We identify the forward and backward fluxes corresponding to the reaction with rates k^±,n as \({J}^{+,n}={k}^{+,n}{c}_{n}^{2}\) and J^−,n = k^−,nc_n+1. The rate of entropy production corresponding to the reaction is then \({\mathop{\sigma }\limits^{.}}_{n}=({J}^{+,n}-{J}^{-,n})\log ({J}^{+,n}/{J}^{-,n})\) (refs. ^{161,162,163,164}). We can then evaluate this quantity and the total rate of entropy production \(\mathop{\sigma }\limits^{.}={\sum }_{n}{\mathop{\sigma }\limits^{.}}_{n}\) from the steady-state distributions c_n obtained numerically (Extended Data Fig. 3f). The rate of entropy production vanishes when the system is isolated (J_in = J_out = 0), and increases as a function of the flux going through the system (which is equal to J_in as long as there is a stationary state).