# Pattern formation by turbulent cascades – Nature

Mar 20, 2024
Strange India

### 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 k2 (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 rule58. Both normal and odd viscosities are integrated exactly using integrating factors. The forcing f(t, k) acts on a band of wavenumbers k [kin, kin + 1] with random phases that are delta-correlated in space and time, ensuring a constant average energy injection rate ϵ = uf. 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 turbulence2,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 equation2,4

$${\partial }_{t}E=-2\nu {k}^{2}E-T+F$$

(6)

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

$$T({\bf{k}},t)={\rm{I}}{\rm{m}}\sum _{{\bf{k}}+{\bf{p}}+{\bf{q}}={\bf{0}}}{v}_{i}^{\ast }({\bf{k}},t){P}_{ij}({\bf{k}}){q}_{\ell }{v}_{\ell }^{\ast }({\bf{p}},t){v}_{j}^{\ast }({\bf{q}},t).$$

(7)

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νk2E represents standard viscous dissipation. In equation (7), the sum runs on momenta p and q such that k + p + q = 0, and Pij(k) = δij − kikj/k2 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

$${\bf{v}}(t,{\bf{x}})={{\bf{h}}}^{\pm }({\bf{k}})\,{{\rm{e}}}^{{\rm{i}}{\omega }_{\pm }({\bf{k}})t+{\rm{i}}{\bf{k}}\cdot {\bf{x}}}+\,{\rm{c.c.}}$$

(8)

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. 60) with frequency

$${\omega }_{\pm }({\bf{k}})=\pm {\nu }_{{\rm{odd}}}{k}_{z}| k| .$$

(9)

Taking into account normal viscosity leads to an additional exponential decay of the waves with the rate −νk2 (Supplementary Information). In particular, we note that the linearized Navier–Stokes equation does not exhibit any linear instability. By construction, kh±(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

$${\bf{v}}(t,{\bf{x}})=\sum _{{\bf{k}}}\sum _{s=\pm }{v}_{s}(t,{\bf{k}}){{\bf{h}}}^{s}({\bf{k}}){{\rm{e}}}^{{\rm{i}}\omega t+{\rm{i}}{\bf{k}}\cdot {\bf{x}}}$$

(10)

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

$${\partial }_{t}{v}_{{s}_{k}}=\sum _{\begin{array}{c}{\bf{k}}+{\bf{p}}+{\bf{q}}={\bf{0}}\\ {s}_{p},{s}_{q}=\pm \end{array}}{C}_{k| p,q}{{\rm{e}}}^{{\rm{i}}[\omega ({\bf{k}})+\omega ({\bf{p}})+\omega ({\bf{q}})]t}{v}_{{s}_{p}}^{* }{v}_{{s}_{q}}^{* }-\nu {k}^{2}{v}_{{s}_{k}}+{f}_{{s}_{k}}$$

(11)

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

$${C}_{k| p,q}=-\frac{1}{4}({s}_{p}p-{s}_{q}q){[({{\bf{h}}}^{{s}_{p}}({\bf{p}})\times {{\bf{h}}}^{{s}_{q}}({\bf{q}}))\cdot {{\bf{h}}}^{{s}_{k}}({\bf{k}})]}^{* }$$

(12)

satisfy Ckp,q = Ckq,p.

#### Helicity and energy conservation of inviscid invariants

In terms of the components v±(k), energy and helicity, respectively, read4,60

$$E=\sum _{{\bf{k}}}(| {v}_{+}({\bf{k}}){| }^{2}+| {v}_{-}({\bf{k}}){| }^{2})$$

(13)

$$H=\sum _{{\bf{k}}}k(| {v}_{+}({\bf{k}}){| }^{2}-| {v}_{-}({\bf{k}}){| }^{2}).$$

(14)

A direct calculation shows that36,60

$${C}_{k| p,q}+{C}_{p| q,k}+{C}_{q| k,p}=0$$

(15)

and

$${s}_{k}k{C}_{k| p,q}+{s}_{p}p{C}_{p| q,k}+{s}_{q}q{C}_{q| k,p}=0$$

(16)

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,

$${\partial }_{t}E(k)={v}_{{s}_{k}}^{* }{\partial }_{t}{v}_{{s}_{k}}+\,{\rm{c.c.}}$$

(17)

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

$${\partial }_{t}E(k)=\sum _{\begin{array}{c}{\bf{k}}+{\bf{p}}+{\bf{q}}={\bf{0}}\\ {s}_{p},{s}_{q}=\pm \end{array}}{C}_{k| p,q}{{\rm{e}}}^{{\rm{i}}[\omega ({\bf{k}})+\omega ({\bf{p}})+\omega ({\bf{q}})]t}{v}_{{s}_{k}}^{* }{v}_{{s}_{p}}^{* }{v}_{{s}_{q}}^{* }+\,{\rm{c.c.}}\,$$

(18)

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 kz = 0).

#### Resonant manifold

The 2D modes with kz = 0 form a so-called slow manifold, or resonant manifold, that contributes to most of the nonlinear energy transfer. Furthermore, isolated triads with kz ≠ 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 kz = 0, leading to an accumulation of energy in the slow manifold, enhancing the two-dimensionalization of the flow4,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 layers4,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. 39 (see ref. 42 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 τ3. 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 τ3(k). Therefore, using hypothesis 3, we write

$${\epsilon }=\varepsilon (k)=A{\tau }_{3}(k){k}^{\alpha }{E}^{\beta }$$

(19)

where A is a constant. The exponents are found using dimensional analysis (with [E] = L3T−2, [ϵ] = L2T−3, [k] = L−1, [τ3] = T), which yields α = 4 and β = 2.

In Fig. 2, we argue that the eddy turnover time τE is the relevant timescale when kkodd (so we set τ3 = τE in the expression above), whereas the frequency of odd waves ω is the relevant timescale when kkodd (so we set τ3 = ω−1). 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]/vk. As E(k) is the shell-average of vk, we have dimensionally $$E(k)\propto {v}_{k}^{2}/k$$, so vk = [kE(k)]1/2. Putting everything together, we end up with equations (3) and (4) of the main text.

#### Scaling relation for kc

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

$$\frac{E(k)}{{E}_{0}(k)}\propto \left\{\begin{array}{ll}1 & \,\mathrm{for}\,\,k\ll {k}_{{\rm{odd}}},\\ {\left(k/{k}_{{\rm{odd}}}\right)}^{s} & \,\mathrm{for}\,\,k\gg {k}_{{\rm{odd}}},\end{array}\right.$$

(20)

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

$$E(k)\propto \left\{\begin{array}{ll}{{\epsilon }}^{2/3}{k}^{-5/3} & \,\mathrm{for}\,\,k\ll {k}_{{\rm{odd}}},\\ {{\epsilon }}^{2/3}{k}^{-5/3}{\left(k/{k}_{{\rm{odd}}}\right)}^{s} & \,\mathrm{for}\,\,k\gg {k}_{{\rm{odd}}}.\end{array}\right.$$

(21)

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 kc from the balance between injection and dissipation. Neglecting contributions to the dissipation from wavenumbers k < kodd (where there is no meaningful change from Kolmogorov scaling), we obtain

$${\epsilon }\propto {\int }_{{k}_{{\rm{o}}{\rm{d}}{\rm{d}}}}^{{k}_{{\rm{c}}}}\nu {k}^{2}E(k){\rm{d}}k.$$

(22)

Assuming koddkc, this yields

$${\epsilon }\propto \nu \,{{\epsilon }}^{2/3}\,{k}_{{\rm{c}}}^{4/3}{({k}_{{\rm{c}}}/{k}_{{\rm{o}}{\rm{d}}{\rm{d}}})}^{s},$$

(23)

resulting in the scaling relation for the peak condensation

$${k}_{{\rm{c}}}\propto {\left({{\epsilon }}^{1/3}{\nu }^{-1}{k}_{{\rm{odd}}}^{s}\right)}^{\frac{1}{4/3+s}}\propto {\left({k}_{\nu }^{4/3}{k}_{{\rm{odd}}}^{s}\right)}^{\frac{1}{4/3+s}},$$

(24)

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

For s = 2/3, we find

$${k}_{{\rm{c}}}\propto {\left({k}_{\nu }^{4/3}{k}_{{\rm{odd}}}^{2/3}\right)}^{1/2}\propto {{\epsilon }}^{1/4}{\nu }^{-1/2}{\nu }_{{\rm{odd}}}^{-1/4}$$

(25)

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 effect62,63,64,65,66,67 generated by the usual viscosity.

A coarse estimate of the height of the peak in E(k)/E0(k) can be obtained by evaluating equation (4) (to get E(k)) and equation (3) (to get E0(k)) at k = kc 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 viscosity62,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π/kc 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 effects68 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 colloids69,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. 71 and references therein). We note, however, that finite energy fluxes have been reported in certain cases56,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 added71. 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 turbulence81,82,83,84,85,86 and laminar and turbulent patterns in wall-bounded shear flows87,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, kz), which distinguish the horizontal (perpendicular) directions from the vertical direction91,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 < kodd and then proceeds to condensate anisotropically into the low-kz 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 kc, which is thus primarily visible in the perpendicular directions because of the anisotropic condensation. The dominant vertical scale hence remains closer to kodd. This leads to a crude estimate for the aspect ratio γ of the features in the pattern produced by the odd viscosity as

$$\gamma =\frac{{k}_{c}}{{k}_{{\rm{odd}}}}\propto {\nu }^{-1/2}{\nu }_{{\rm{odd}}}^{1/2}.$$

(26)

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 > kodd. In the region k < kodd, 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) k2 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 colloids32, magnetized graphene10 and magnetized polyatomic gases31,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 scale9,96,97. They have been experimentally measured in polyatomic gases under magnetic fields31,95, spinning colloids32 and magnetized electron fluids10. They have also been predicted in systems, including fluids under rotation98, magnetized plasma11,99,100, quantum fluids96,101,102,103, vortex matter104, sheared granular gases105, assemblies of spinning objects69,106,107,108,109,110,111,112,113,114,115,116,117 and circle swimming bacteria118,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 colloids32 and magnetized graphene)10. 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. 104), leading to an infinite value of νodd/ν. Kinetic theory calculations for magnetized plasma (ref. 99, section 19.44) predict ν = ν0/[1 + x2] and νodd = ν0x/[1 + x2] in which x = 2ωτ with τ is a collision time and ωB is a frequency proportional to the magnetic field B, whereas ν0 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 speed98. In electron gases in graphene, experiments have been performed that validate these theoretical calculations10 (with x = B/B0, where B0 is a reference magnetic field). This results in a ratio νodd/ν = xB. Likewise, in active fluids, theoretical works suggest that νodd is proportional to the rotation speed of the spinning particles106.

### 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) equation93,120,121,122,123,124,125,126

$${\partial }_{t}\omega +J(\psi ,\omega )+\beta {\partial }_{x}\psi =-\alpha \omega +\nu \Delta \omega +{f}_{\omega }$$

(27)

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 Dedalus128.

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 purposes13,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 droplets23,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 equation23,139,140,141 that can exhibit scale-invariant cascades, similar to that present in the Navier–Stokes equation142,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 sizes24,145,146, in which the distribution originates from complex mechanisms, including air turbulence and fluid fragmentation147,148,149,150,151,152.

In our toy model, we consider clusters Mn composed of 2n−1 monomers M1, with n = 1, …, N. This is reminiscent of what is done in shell models of turbulence52, 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 MN. 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

$$\varnothing \,\mathop{\longrightarrow }\limits^{{J}_{{\rm{i}}{\rm{n}}}}\,{M}_{1}$$

(28a)

$$2\,{M}_{n}\mathop{\mathop{\rightleftharpoons }\limits^{{k}_{n}^{+}}}\limits_{{k}_{n+1}^{-}}{M}_{n+1}$$

(28b)

$${M}_{N}\,\mathop{\longrightarrow }\limits^{{J}_{{\rm{o}}{\rm{u}}{\rm{t}}}}\,\varnothing$$

(28c)

in which Mn (n = 1, …, N) represents a cluster of size 2n−1 (M1 represents a monomer), and Jin, Jout and $${k}_{n}^{\pm }$$ are the rates of the corresponding reactions.

The number densities cn of clusters then follow the dynamical equation

$$\frac{{\rm{d}}{c}_{n}}{{\rm{d}}t}=2{k}_{n+1}^{-}{c}_{n+1}-{k}_{n}^{-}{c}_{n}+\frac{1}{2}{k}_{n-1}^{+}{c}_{n-1}^{2}-{k}_{n}^{+}{c}_{n}^{2}+{J}_{{\rm{ext}}}$$

(29)

where

$${J}_{{\rm{ext}}}={\delta }_{n,1}{J}_{{\rm{in}}}-{\delta }_{n,N}{J}_{{\rm{out}}}{c}_{N}$$

(30)

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

We can also consider the mass density ρn = 2n−1m0cn, in which m0 is the mass of a monomer. Multiplying equation (29) with 2n−1m0, we find that the terms with prefactors $${k}_{n}^{\pm }$$ cancel as in a telescoping series. This manifests that equation (29) with Jin = Jout = 0 conserves mass. It is therefore convenient to introduce the fluxes

$${J}^{+}(n)=-{\int }_{1}^{n}\,d{n}^{{\prime} }\left[\frac{1}{2}{k}_{{n}^{{\prime} }-1}^{+}{c}_{{n}^{{\prime} }-1}^{2}-{k}_{{n}^{{\prime} }}^{+}{c}_{{n}^{{\prime} }}^{2}\right]$$

(31)

and

$${J}^{-}(n)=-{\int }_{1}^{n}\,d{n}^{{\prime} }[2{k}_{{n}^{{\prime} }+1}^{-}\,{c}_{{n}^{{\prime} }+1}-{k}_{{n}^{{\prime} }}^{-}\,{c}_{{n}^{{\prime} }}]$$

(32)

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

$$\frac{{\rm{d}}{c}_{n}}{{\rm{d}}t}=-{\partial }_{n}[\,{J}^{+}(n)+{J}^{-}(n)]+{\delta }_{n,1}{J}_{{\rm{in}}}-{\delta }_{n,N}{J}_{{\rm{out}}}{c}_{N}.$$

(33)

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

$${k}_{n}^{+}={\kappa }_{0}^{+}+{\kappa }_{1}^{+}\frac{N-n}{N-1}\,{\rm{and}}\,{k}_{n}^{-}={\kappa }_{0}^{-}+{\kappa }_{1}^{-}\frac{n-1}{N-1}.$$

(34)

Equation (29) is then solved starting from the initial condition cn = 0 for all n using DifferentialEquations.jl (ref. 153) 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 cn 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 Jtot ≡ 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 points154,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 Jin = Jout = 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 system46,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 cascade159.

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 networks161,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−,ncn+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 cn obtained numerically (Extended Data Fig. 3f). The rate of entropy production vanishes when the system is isolated (Jin = Jout = 0), and increases as a function of the flux going through the system (which is equal to Jin as long as there is a stationary state).