Nonequilibrium steady solutions of the Boltzmann equation
Davide Proment, Miguel Onorato, Pietro Asinari, Sergey Nazarenko
NNonequilibrium steady solutions of the Boltzmann equation
Davide Proment,
1, 2, ∗ Miguel Onorato, Pietro Asinari, and Sergey Nazarenko Dipartimento di Fisica Generale, Universit`a di Torino,Via Pietro Giuria 1, 10125 Torino, Italy INFN, Sezione di Torino, Via Pietro Giuria 1, 10125 Torino, Italy Dipartimento di Energetica, Politecnico di Torino,Corso Duca degli Abruzzi 24, 10129 Torino, Italy Mathematics Institute, The University of Warwick, Coventry, CV4-7AL, UK (Dated: September 20, 2018)
Abstract
We report a study of the homogeneous isotropic Boltzmann equation for an open system. Weseek for nonequilibrium steady solutions in presence of forcing and dissipation. Using the languageof weak turbulence theory, we analyze the possibility to observe Kolmogorov-Zakharov steady dis-tributions. We derive a differential approximation model and we find that the expected nonequilib-rium steady solutions have always the form of warm cascades. We propose an analytical predictionfor relation between the forcing and dissipation and the thermodynamic quantities of the system.Specifically, we find that the temperature of the system is independent of the forcing amplitudeand determined only by the forcing and dissipation scales. Finally, we perform direct numericalsimulations of the Boltzmann equation finding consistent results with our theoretical predictions.
PACS numbers: 47.27.Gs, 05.70.Ln, 47.70.Nd ∗ Electronic address: [email protected]; URL: a r X i v : . [ n li n . C D ] J a n . INTRODUCTION Systems in a steady state are characterized by observables that do not change in time;they can be either in equilibrium or out of equilibrium. Systems in nonequilibrium steadystates have net currents (fluxes): examples of nonequilibrium steady-state systems includean object in contact with two thermal sources at different temperatures, for which the cur-rent is a heat flux; a resistor with electric current flowing across it; the kinesin-microtubulesystem, for which kinesin motion is the current. Most biological systems, including molecu-lar machines and even whole cells, are in nonequilibrium states [1]. In particular, biologicalsystems rely on a continuous flux of energy and/or particles supplied by some proper envi-ronmental reservoirs.In statistical mechanics, investigating the general properties of a system in contact withreservoirs, namely an open system, is a long lasting problem (e.g. see the second problemdiscussed by E.H. Lieb on the occasion of the award of the Boltzmann medal [2]), eventhough these theoretical challenges are sometimes neglected in applied engineering at large.The difficulties arise from the fact that finding the large deviation functional for a stationarystate with fluxes is still an open problem (see [3] and references therein). In the present work,for focusing our attention and considering an affordable goal, we consider the kinetic theoryof gases. In particular, we consider a system composed of a large number of interactingparticles, comparable to the Avogadro number. The Boltzmann kinetic equation (BKE)describes the time evolution of the single-particle distribution function, which provides astatistical description of the positions and velocities (momenta) of the gas molecules. Thisintegro-differential kinetic equation, proposed by Boltzmann at the end of the XIX century,has been derived starting from the phase-space Liouville equation, assuming the stosszahlansatz [4]. Its equilibrium state, which maximizes the entropy measure, is the Maxwell-Boltzmann distribution. In case of small deviations from the local equilibrium, it is possibleto systematically derive hydrodynamic equations for macroscopic quantities of the system;e.g., in the lowest order approximation for small departures from equilibrium, the Navier-Stokes equations [4].Kinetic equations have also been studied in the framework of wave turbulence theory25] where it has been shown that other solutions with respect to thermodynamic solutionscan be stationary states of the system, in case of external forcing and dissipation. Thesedistributions, which have usually the form of power-laws in momentum space, are calledKolmogorov-Zakharov (KZ) and they represent constant flux of conserved quantities similarto the Kolmogorov energy cascade in strong Navier-Stokes turbulence [6, 7]. These solu-tions, named cascade solutions, become important when considering an open system, i.e.with forcing and dissipation terms. They have been studied for a great variety of weaklynonlinear dispersive models: examples can be found in water waves [8–10], internal waves[11], nonlinear optics [12], Bose-Einstein condensation [13–15], magnetohydrodynamics [16].An out of equilibrium description of the Boltzmann equation using the KZ solutions wasfirst devised in [17] considering different types of interaction potential between particles.Problems of interaction locality scale-by-scale and wrong flux direction were pointed out. Inparticular in [18] Kats showed that for all realistic physical situations the direction of thecascades in the system is always in the wrong orientation with respect to the one predictedby the Fjørtoft theorem [33]. When a formal KZ solution has a flux direction contradictingwith the Fjørtoft theorem, this spectrum (even if local) cannot be established because itcannot be matched to any physical forcing and dissipation at the ends of the inertial range.For example in [12], the particle cascade KZ solution was found to be of this type in the two-dimensional nonlinear Schr¨odinger equation model the authors argued that in this case theKZ solution is not achievable and a mixed state, with both a cascade and a thermodynamiccomponents were proposed. Another example of mixed cascade-thermodynamic states canbe found in the context of three-dimensional Navier-Stokes turbulence [19], where such mixedstates were called warm cascades [34].The present manuscript will focus on warm cascades found in the homogenous isotropicBoltzmann equation (HIBE) and in particular it will answer to the following importantquestions. • What is precisely the relation between the conserved quantity fluxes and the thermo-dynamics quantities of the system? • How does this relation depends on the forcing and dissipation rates and acting scales?To answer the above questions we will perform numerical simulations of the homogeneousisotropic Boltzmann equation with forcing and dissipation. We will then use a diffusion3pproximation model (DAM) to derive analytical predictions on how the thermodynamicquantities, temperature and chemical potential, are related to fluxes, forcing and dissipativescales. We will then test these predictions by numerically simulating both DAM and thecomplete homogenous isotropic Boltzmann equation.The work is organized as follows: in Section II we review the properties of the Boltzmannequation for the homogeneous isotropic case; in Section III we introduce DAM and we derivethe analytical predictions; Section IV is dedicated to numerical results of DAM and HIBE;in Section V we draw the conclusions. A set of Appendixes also provide detailed calculationsof those results which are briefly reported in the main text.
II. THE BOLTZMANN KINETIC EQUATION
The Boltzmann kinetic equation describes the time evolution of the single-particle dis-tribution function, which provides a statistical description for the positions and momentaof the gas molecules: the function n ( x , k , t ) express a probability density function in theone-particle phase space R d x × R d k with respect to time, where d is the dimension. Note thatwe denote the momentum variable with the letter k instead of the conventional p to followthe common notation of wave turbulence [5]. The Boltzmann equation takes the followingform: ∂n∂t ( x , k , t ) + k m · ∂n∂ x ( x , k , t ) = I coll ( x , k , t ) , (1)where I coll = (cid:90) + ∞−∞ W [ n ( x , k , t ) n ( x , k , t ) − n ( x , k , t ) n ( x , k , t )] d k d k d k (2)sums the effect of the two-body collisions of particles with all possible values of momenta.The form of the collision integral we are reporting is equivalent to the standard one andcorresponds to Eq. (4.18), page 64 in Cercignani’s book [4]. Here W describes syntheticallythe scattering amplitude transition 2 → W is W = Γ δ ( k + k − k − k ) δ ( | k | + | k | − | k | − | k | ) , (3)where δ -functions assure conservation of the total momentum and the total kinetic energy(which is proportional to | k | ) of incoming and outgoing particles. The collision probability,4xpressed by Γ ≡ Γ( k , k | k , k ) ≥
0, is invariant under permutations { , } → { , } , { , } → { , } , and { , } → { , } . In the present paper we will consider the case ofthree-dimensional rigid spheres with diameters σ and mass m , for which Γ simply results inΓ = 2 σ /m [4]. For other interaction potentials, as Coulomb or Born approximation, referto [18, 20].For the purposes of our work, we consider a homogeneous and isotropic (in physical space R d x ) system with the one-particle probability density function independent of x and its mo-mentum dependency coming only via the modulus k = | k | , so n ( x , k , t ) → n ( k, t ). It is usefulto express the distributions in the energy space ω i = | k i | where we use again the notation ω for the energy in analogy with wave turbulence. Then, the particle density in ω -space satisfiesthe relation (cid:82) N ( ω, t ) dω = (cid:82) n ( k, t ) d k or, in the other words, N ( ω, t ) = n ( ω, t ) Ω ω d − (cid:12)(cid:12) dkdω (cid:12)(cid:12) ,where Ω is the solid angle. After these considerations Boltzmann equation (1) simplifies tothe homogeneous isotropic Boltzmann equation (HIBE): ∂N ∂t = (cid:90) ∞ S ( n n − n n ) δ ( ω + ω − ω − ω ) dω dω dω , (4)where we denote for brevity N i = N ( ω i , t ) and n i = n ( ω i , t ), and the functional S = ( ω ω ω ω ) d − (cid:12)(cid:12)(cid:12)(cid:12) dk dω (cid:12)(cid:12)(cid:12)(cid:12) (cid:12)(cid:12)(cid:12)(cid:12) dk dω (cid:12)(cid:12)(cid:12)(cid:12) (cid:12)(cid:12)(cid:12)(cid:12) dk dω (cid:12)(cid:12)(cid:12)(cid:12) (cid:12)(cid:12)(cid:12)(cid:12) dk dω (cid:12)(cid:12)(cid:12)(cid:12) (cid:90) Ω Γ δ ( k + k − k − k ) d Ω d Ω d Ω d Ω (5)takes into account the change of coordinates and the average over solid angles. Here-after, we always consider a three-dimensional gas of hard-sphere particles in a non-dimensional form with m = 1 and σ = 8 . Then the functional simply results in S = 2 π min (cid:2) √ ω , √ ω , √ ω , √ ω (cid:3) (see Appendix A for details of the angular integration).The HIBE has two conserved quantities, the mass and energy densities, ρ M = (cid:90) + ∞−∞ n ( ω, t ) d k = 2 π (cid:90) + ∞ n ( ω, t ) √ ωdω,ρ E = (cid:90) + ∞−∞ n ( ω, t ) k d k = 2 π (cid:90) + ∞ n ( ω, t ) ω dω. (6)Note that ρ M and ρ E are always constant in time for any distribution n and interactionpotential, due to the fact that collisions are 2 → . Steady solutions
1. Equilibrium in a closed system
The HIBE (4) is an integro-differential equation with no general analytic solution. Itis easy, however, to look for steady (time independent) solutions. In closed system, i.e.without forcing and/or dissipation mechanisms, the only steady solution corresponds to thethermodynamic equilibrium described by the Maxwell-Boltzmann (MB) distribution, n MB ( ω ) = e − ω + µT = Ae − ωT , (7)where A = e − µT and constants µ and T have the meaning of the chemical potential andthe temperature respectively (we consider the natural unit system, where the Boltzmannconstant is one). Validation is trivial by plugging (7) into (4): for any value of T , µ and theinteraction potential S , the δ -function assures that the integrand is zero. Moreover, thetotal mass density of the system is ρ M = A ( πT ) , the total energy density is ρ E = Aπ T ,and any other moment of ω , due to the bi-parametric nature of the MB distribution, is afunction of ρ M and ρ E . The H theorem states that in a closed system any out of equilibriumdistribution with defined mass and energy densities will always relax to the MB distributionhaving same ρ M and ρ E .In Fig. 1 we show a numerical simulation of the HIBE with initial condition given bya Gaussian function centered around a particular value of energy; as it is clear from thefigure, the initial condition relaxes to the MB distribution. The numerical algorithm usedto perform this simple example will be discussed in Section IV B. We can observe that theinitial condition evolves reaching an equilibrium MB distribution: the exponential behaviorbecome evident by observing the inset where we plot it lin-log plot scale. Moreover by fittingthe results with the MB function we can find the thermodynamic quantities A and T : thosecorrespond exactly to ones expected knowing initial mass and energy densities (note thatnow integrals (6) are evaluated from 0 to a finite value of ω due to numerical finiteness of ω -space). 6 .5 ⋅ -4 ⋅ -4
0 20 40 60 80 100 n ( ω ) ω initialdistribution A=1.8e-04,T=2.6e+01 FIG. 1: Numerical computation of HIBE with an initial Gaussian shaped distribution (continuosblack line): intermediate states are shown with gray lines and final steady distribution with dashedblack line. The latter has a MB behavior (7) with fitted parameters A and T printed in figure.The inset shows the same plot in lin-log scale.
2. Nonequilibrium steady states
Now, what can we expect in an open system driven by external forcing and dissipationmechanisms? We will answer this question keeping in mind the main results of the waveturbulence theory. Part of this theory is dedicated to study steady solutions to kineticequations in the power-law form, n ( ω ) ∼ ω − x , where the constant x assumes different valuesdepending on the considered wave system. It is sometimes possible to find the so-calledKolmogorov-Zakharov (KZ) solutions n KZ ( ω ) ∼ ω − x which correspond to constant fluxesof conserved quantities through scales. The KZ distribution always appears in a range ofscales, known as inertial range, between the forcing and dissipation were the source and sinkare located.As already mentioned, the HIBE conserves the number of particles and the energy, andso one could expect to observe two turbulent KZ cascades. The KZ exponent x can beevaluated by applying the standard Zakharov transformations [5], by dimensional analysis[21], or by using the method (equivalent to Zakharov transformation) proposed by Balk[22]. We have chosen the last one and the complete analytical calculations are presented inAppendix B. The KZ exponents depend on the scaling behavior of the scattering term Γ d of the particle system. For the particular case of three-dimensionalhard spheres we haveconstant particle flux η = ⇒ n KZ ( ω ) ∼ ω − , constant energy flux (cid:15) = ⇒ n KZ ( ω ) ∼ ω − . (8)The simplest way to mimic an open system where steady nonequilibrium distributions ofthe form of turbulent KZ solutions can be establish is to consider a forced-dissipated HIBE ∂N∂t ( ω , t ) = I coll ( ω , t ) + F ( ω ) − D ( ω ) N ( ω ) . (9)The forcing F is constant in time and very narrow near a particular energy value ω f : with thischoice the incoming fluxes of particles η and energy (cid:15) roughly satisfy relation (cid:15) = ω f η . Thedissipation term D is implemented as a filter which removes, at each iteration time, energyand particles outside of the domain ω ∈ ( ω min , ω max ). Further details on the numericalscheme are explained in Section IV B. What happens if we try to solve numerically suchforced/damped integro-differential equation?In Fig. 2 and Fig. 3 we plot the nonequilibrium steady states obtained with numericalsimulations of the HIBE with forcing and dissipation; the initial conditions are characterizedby n ( ω, t = 0) = 0 . The parameters in the simulations are ω min = 5, ω max = 195, the forcingrate F = 10 − . In Fig. 2 forcing is located at ω f = 22 and in Fig. 3 at ω f = 182. Nopower-law distributions, and so no KZ solutions (8), are observed (note that both plots arein lin-log scales), but instead one can see weakly perturbed exponential curves. We canattempt to measure the quantities T and A in (7) by fitting our numerical curves; however,those are not perfect straight lines (in the lin-log plot) and left and right branches withrespect to forcing scale may give different results. For such reason we will denote by ( · ) L the quantities evaluate on the left brach and with ( · ) R the right ones.Another example we analyze is the case where we fix the forcing and dissipative scalesand change the forcing rate. Numerical results for final steady states evaluated for threedifferent forcing amplitudes, F = 10 − , − , − , are presented in Fig. 4. The effect ofincreasing the amplitude F results in an upward shift of the curves. Therefore, qualitatively,the temperature appears to be the same for each value of the flux. The only difference is8 -7 -6 -5 -4 -3
0 50 100 150 200 n ( ! ) ! forcing scaleA L =1.6e-04, T L =4.3e+01A R =2.0e-04, T R =2.8e+01 FIG. 2: An example of HIBE (9) steady state shown with black line in lin-log scale. Simulationparameters are: ω min = 5, ω f = 22, ω max = 195 and F = 10 − , and ω cutoff = 200. The dashedand point/dashed lines are left and right branch best fits obtained with the MB distribution (7):fitting parameters are reported in label. -7 -6 -5 -4 -3
0 50 100 150 200 n ( ! ) ! forcing scaleA L =7.7e-05, T L =5.6e+01A R =3.7e-04, T R =3.7e+01 FIG. 3: An example of HIBE (9) steady state shown with black line in lin-log scale. Simulationparameters are: ω min = 5, ω f = 182, ω max = 195 and F = 10 − , and ω cutoff = 200. The dashedand point/dashed lines are left and right branch best fits obtained with the MB distribution (7):fitting parameters are reported in label. -7 -5 -3
0 20 40 60 80 100 n ( ω ) ω F=10 -4 F=10 -5 F=10 -6 FIG. 4: Steady states of HIBE (9) in lin-log scale obtained for different values of the forcing rate F . Parameters are ω cutoff = 200, ω min = 5, ω f = 21 and ω max = 95. ρ E ( t ) t F=10 -4 F=10 -5 F=10 -6 FIG. 5: Total energy densities ρ E ( t ) in function of time for different forcing rates F . For the systemparameters refer to ones in Fig. 4. the speed at which the system, initially empty, reaches its steady state. In fig. 5 we showthe energy density evolution (same line styles corresponds to same systems).After these preliminary numerical results, a lot of questions can be posed. Why noKZ constant flux solutions are observed but just small deviations from MB distributions?What happens when forcing or dissipation scales are changed? What is in general the10elation between physical quantities such as fluxes, forcing and dissipation scales and theMB parameters? The aim of this manuscript is to provide explanations to such phenomenaand answer these questions. B. Locality of interactions
For the KZ spectra to be valid mathematical (and therefore physically relevant) solutions,it is necessary that they satisfy the locality condition. A spectrum is local when the collisionintegral converges. In other words, non-locality means that the collision integral is notweighted scale by scale but most of the contributions come from the limits of integrationcorresponding to the ends of the inertial range. Physically, the non-locality is in contradictionwith the assumption that the flux of the relevant conserved quantity in the inertial range iscarried only by the nearest scales. Mathematically, locality guaranties that the KZ spectrumis a valid solution in an infinite inertial range, which is not guarantied a priori becauseZakharov transformation is not an identity transformation and could, therefore, lead tospurious solutions.For the HIBE case, locality depends on the particular interaction potential, which affectsthe scaling of Γ , and on the dimensionality of the system - for detailed calculations seeAppendix B. Locality is not always found for both KZ solutions: for example for the Coulombpotential only the energy cascade is local, as shown in [17]. In the case of three-dimensionalhard spheres considered in the present work, the criterion of locality is never satisfied forany of the two KZ solutions, which means that these solutions are un-physical and irrelevantin this model.
C. The flux directions
Besides locality, another important requirement for establishment of the KZ spectra isthe correctness of the flux directions for the respective conserved quantities. In a systemwhere two quantities are conserved, the following Fjørtoft-type argument is used to establishwhich quantity must have a direct or an inverse cascade.11 . The Fjørtoft argument
Consider an open system where forcing scale ω f is widely separated from a low- ω dissipa-tion scale ω min and a high- ω dissipation frequency ω max , thus ω min (cid:28) ω f (cid:28) ω max . Becausethe energy density in the ω -space is different from the particle density by factor ω , theforcing rate of the energy (cid:15) is related with the forcing rate of the particles η as (cid:15) ∼ ω f η .Suppose that some energy is dissipated at the low scale ω min at a rate comparable withthe forcing rate (cid:15) . But then the particles would have to be dissipated at this scale at therate proportional to (cid:15)/ω min ∼ η ω f /ω min (cid:29) η , which is impossible in steady state becausethe dissipation cannot exceed the forcing. Thus we conclude that in the steady state theenergy must dissipate only at ω max . By a symmetric contradiction argument one can eas-ily show that the only place where the particles can be dissipated in such systems is ω min .This means that energy must have a direct cascade (positive flux direction) and particles aninverse cascade (negative flux direction).
2. Flux directions in the HIBE
It has been proved in [18], see also Appendix B for details, that fluxes of the KZ solutionsfor all types of the interaction coefficient Γ have always the wrong directions with respect tothe Fjørtoft argument requirements (in the case x > n ( ω, t ) ∼ ω − x and plotting them as functions of x for a fixed ( ω, t ), see Fig. 6. Threeexponents x correspond to steady solutions of HIBE: the particle equipartition x eq = 0, theKZ particle cascade x η and the KZ energy cascade x (cid:15) . As shown in Appendix B, we knowthat (cid:15) = 0 on the particle cascade, η = 0 on the energy cascade, whereas in the equipartitionboth fluxes are zero, i.e. (cid:15) = η = 0. We also know that for large negative x (large positiveslope) both fluxes must be negative, as such a steep unsteady spectrum would evolve tobecome less steep, toward equipartition. Note that always x η < x (cid:15) , when x >
0. Now wecan sketch the particle and energy fluxes as function of the exponent x as it is done in Fig.6. From this sketch, it can be easily understood that whenever the condition x η < x (cid:15) isvalid, the particle flux will be positive and the energy flux will be negative, contradicting theFjørtoft argument. This means that one cannot match these formal KZ solutions, obtained12
0 0 η ( x ) , ε ( x ) x η (x) ε (x) FIG. 6: Energy and particle fluxes on power-law solutions n ( ω, t ) ∼ ω − x as functions of x for thethree-dimensional hard sphere model. for an infinite inertial range, to any physical forcing or dissipation at the ends of a large(but finite) inertial range.What is then happening when fluxes have wrong direction? It has been observed inoptical wave turbulence [12] that the pure KZ spectra are not established in these casesand one has to expect a mixed solution where both a flux and a thermal components arepresent. Such mixed states are quite common for turbulent systems of different kinds,including strong Navier-Stokes turbulence and have been named warm cascades [19]. Suchcascades were obtained within the Leith model (which belongs to the class of the differentialapproximation models) as exact analytical solutions. III. DIFFERENTIAL APPROXIMATION MODEL
Numerical integration of the Boltzmann collision integral is very challenging because thenumber of degrees of freedom grows as a polynomial. A great simplification comes fromthe isotopic assumption, which reduces the degrees of freedom from N to N ( N is thenumber of points needed to describe the distribution). However spanning a large numberof momentum scales is still difficult. For those reasons, some approximations to the kineticequations were proposed in order to increase the range of modeled scales, see for example[23]. 13 great simplification is to replace the collision integral operator of the kinetic equation bya nonlinear differential operator which mimics the basic scalings of the original one and yieldsthe same steady solutions. The HIBE then results in a nonlinear partial differential equationcalled the differential approximation model (DAM). Such models have been proposed tosimulate turbulence in different research fields: for example in water waves [24], in nonlinearoptics [12], in strong Navier-Stokes turbulence [25, 26], in Kelvin quantum turbulence [27],in astrophysics (Kompaneets equation) [28], in semiconductors [29]. Replacing the integraloperator by a differential one amounts to assuming locality of the scale interactions, whichmeans the relevant distributions must be local for DAM to have a good predictive power.We mentioned in Section II that for hard sphere Boltzmann equation the pure KZ spectraare non-local and so no DAM would be advisable. However, we observed in some examples(Fig. 2 and Fig. 3) that the relevant solutions in this case are not pure KZ spectra butdistributions which are close to MB, warm cascades, which appear to be local. Thus, we usethe DAM for describing this system, after which we will validate our results by computingthe full HIBE.For the dual cascade systems, such as gravity water waves [30], nonlinear Schr¨odingerequation [12], two-dimensional hydrodynamic turbulence [26], Kelvin waves [27] or HIBEconsidered here, DAM has always the form of a dual conservation law, ∂ t N ( ω, t ) = ∂ ωω R [ n ( ω, t )] , (10)where R is a nonlinear second-order differential term whose details depend on the particularmodel. This equation can be written as a continuity equation for the particle invariant, ∂ t N ( ω, t ) + ∂ ω η ( ω, t ) = 0 , with the particle flux η ( ω, t ) = − ∂ ω R [ n ( ω, t )] . (11)Moreover, equation (10) can be written as a continuity equation for the energy [12], ∂ t [ N ( ω, t ) ω ] + ∂ ω (cid:15) ( ω, t ) = 0 , with the energy flux (cid:15) ( ω, t ) = R [ n ( ω, t )] − ω∂ ω R [ n ( ω, t )] . (12)14e are now able to find the functional R by requiring it to yield the MB distribution (7)and the KZ spectra (8) as steady state solutions of DAM (10). These constraints lead to R [ n ( ω, t )] = − S ω n ( ω, t ) ∂ ωω log n ( ω, t ) , (13)where S is a constant. A formal derivation starting from the kinetic equation can be obtainedfollowing [12, 29]. It is trivial to verify by substitution that KZ solutions (8) correspondto constant fluxes through scales. Namely, the KZ particle cascade has a constant particleflux and zero energy flux while the KZ energy cascade viceversa. Let us again consider thethe flux directions on the KZ distributions, but now using DAM. Substituting power-lawspectra n = c ω − x into (13), equations (11) and (12) yield η = c S x (9 / − x ) ω / − x (cid:15) = c S x (11 / − x ) ω / − x . (14)By plotting η and (cid:15) as functions of the exponent x at fixed ω , we arrive again at Fig. 6.Note that it is by using DAM such plot was obtained. Once again we note that the particleand the energy fluxes on the respective KZ solutions ( x = 7 / x = 9 /
4) have wrongdirections with respect to the Fjørtoft argument.The beauty of the DAMs is the possibility to solve numerically the system for widefrequency ranges and, therefore, to find clear scalings. In particular, such models are veryefficient for finding constant steady flux solutions because they become simple ordinarydifferential equations (ODEs). In the following we will present some analytical results forsuch steady states.
A. Constant energy flux: direct cascade
We will now find an ODE that describes a constant direct energy cascade (cid:15) with no fluxof particles, which we call ODE- (cid:15) . According to Fjørtoft argument, this implies a largedirect-cascade inertial range. Putting η = 0 in (11) and (12), we haveconstant energy flux = ⇒ (cid:15) = R ( ω, t ) = const. (15)15sing (13), we arrive at the following Cauchy problem (cid:15) = − S ω n ( ω ) ∂ ωω log n ( ω ) ,n ( ω ) = n ,∂ ω n ( ω ) = n (cid:48) , (16)where we have chosen the boundary conditions fixing the values of the distribution and itsderivative at the same point ω (e.g. at the forcing scale) for ease of numerical solution.If we solve numerically in ω -forward the ODE- (cid:15) for different values of the energy flux wefind curves presented in Fig. 7. Here we do not want to discuss the details (it will be donewidely in Section IV), but just remark that the solutions follow the MB distribution andsuddenly change behavior going very fast to a zero value of the distribution. We will callthis rapid change a front solution .
1. Compact front behavior
It is possible to find a front solution for the equation (15) describing the behavior nearthe dissipation scale. Let us seek for a front solution which in the vicinity of a certain point ω max behaves like n ( ω ) = B ( ω max − ω ) σ . If we plug this expression into (15) and take thelimit ω → ω max we find that to satisfy this equation in the leading order in ( ω max − ω ) wemust have σ = 1 B = (cid:113) (cid:15)S ω / = ⇒ n ( ω ) = (cid:114) (cid:15)S ω / ( ω max − ω ) . (17)Thus, the front solution is linear in the vicinity of ω max with a slope depending on thedissipation scale ω max and the value of the energy flux (cid:15) . Note that the compact frontbehavior at the dissipation scale is typical for DAM. We will soon discover that ω max is avery useful physical parameter which allows us to find a link between the temperature, thechemical potential and the energy flux in the forced-dissipated system.
2. Kats-Kontorovich correction
Lets summarize our preliminary observations. We expect a warm cascade, that is adistribution which contains both the flux and the thermal components. We have also foundthat the solution has a compact front which arrests the cascade at the dissipation scale16 max . We will now assume (verifying it later) that in the most of the inertial range thewarm cascade solution is close to the thermodynamic MB distribution and the correctiondue to finite flux is small. We then perform a qualitative matching of the flux-correctedMB distribution to the compact front, and thereby obtain a relation between ω max , T and A in (7). To find the warm cascade solution in the inertial range, we consider the Kats-Kontorovich (KK) correction to the Maxwell-Boltzmann distribution: n ( ω ) = n MB (1 + ˜ n ) = (1 + ˜ n ) A e − ωT , (18)where ˜ n is small, ˜ n (cid:28)
1. By plugging this solution into (15) and linearizing in ˜ n we end upwith the following ODE- (cid:15) for the correction (cid:15) ω − A − e ωT = − S ∂ ωω ˜ n. (19)
3. Matching
We will now match the KK correction to the front solution. The basic idea is to forcethe KK solution to satisfy the n ( ω max ) = 0 and to have at ω max the same slope as the frontsolution. Detailed calculation is presented in Appendix C. The prediction results in: (cid:15) = S ω max A e − ω max T . (20)This relation is very important because it gives an analytical relation between the thermo-dynamic quantities T and A in terms of the energy flux (cid:15) and the dissipation scale ω max .However we note that our matching is only qualitative, because the KK correction is sup-posed to be small which is not the case near the front. Thus, the relation (20) is approximateand we do not expect it to hold precisely.
4. Alternative approach to find ω max Another simple way to find a prediction for the value of ω max is the following. As weexpect to observe a warm cascade, we can ask what will be the range where the thermalcomponent will dominate the dynamics. We can simply assume that in most of the inertialrange we will have a distribution n (cid:39) n MB . Note that the MB distribution always has a17ositive concavity, ∂ ωω n ≥
0. On the other hand, we note that our ODE- (cid:15) can be re-writtenas ∂ ωω n = 1 n (cid:20) ( ∂ ω n ) − (cid:15)S ω (cid:21) , (21)from which it is clear that ∂ ωω n may change sign. The point at which ∂ ωω n = 0 can beconsidered as s boundary separating the MB range (with negligible flux correction) and thefront solution (with large flux correction). This boundary can be estimated by a simplesubstitution of the MB distribution to the r.h.s. of (21), which gives (cid:15) = A S ω e − ωT T = g (cid:15) ( ω, A, T ) . (22)As this relation contains the exponential factor which decays very fast (for ω max (cid:29) T , seeAppendix C), it is natural to think that the range at which (cid:15) becomes important appearsvery sharply and is very near to the point ω max . Thus we arrive at the following estimate, (cid:15) = A S ω max e − ω max T T . (23) B. Constant particle flux: inverse cascade
In analogy of what has been done for the direct cascade, we now look for predictionsin the inverse particle cascade η with no flux of energy. The ODE- η that describes such acascade is simple to obtain: by integrating equation (11) once and putting (cid:15) = 0 in (12), wehave: constant particles flux = ⇒ η = − R ( ω, t ) ω = const. (24)This yields the following Cauchy problem, η = S ω n ( ω ) ∂ ωω log n ( ω ) ,n ( ω ) = n ,∂ ω n ( ω ) = n (cid:48) . (25)This problem is most naturally solved backwards in the ω -space, as we are interested inthe inverse cascade. We seek for a solution having a particle flux going from high to lowfrequencies, i.e. η < η → −| η | inequation (25). The Cauchy problem (25) is very similar to (16) with the only difference inthe ω -scaling. Thus we will use the same approach for studying it.18 . Compact front behavior Let us find a front solution for the equation (24). We now expect the front to be onthe left edge of the (inverse cascade) inertial range, i.e. in the vicinity of a certain point ω min < ω f . By plugging n ( ω ) = B ( ω − ω min ) σ expression into (24) and taking the limit ω → ω min , in the leading order in ( ω − ω min ) we have σ = 1 B = (cid:114) | η | S ω / = ⇒ n ( ω ) = (cid:115) | η | S ω / ( ω − ω min ) . (26)Thus, the front solution for the inverse particle cascade is also linear in the vicinity of ω min ,with a slope depending on ω min and the value of the particle flux η .
2. Kats-Kontorovich correction
As previously supposed for the direct energy cascade, we expect in the most of the inverse-cascade range a corrected thermodynamic spectrum and a front solution behavior at the leftend of this range. Let us evaluate the Kats-Kontorovich correction (18), and after thatmatch it to the front solution. By plugging the expression (18) into (24) and linearizing in˜ n we obtain the following ODE- η for the correction, | η | ω − A − e ωT = − S ∂ ωω ˜ n. (27)
3. Matching
Again, we want to match the KK correction to the front solution. The idea is very similarto the previously used for the direct cascade, except for the fact that now the limit taken is ω min (cid:28) T ; for details refer to Appendix D. This results with the following condition on theflux, | η | = S (cid:18) (cid:19) A ω min . (28)19 . Alternative estimate of ω min Again, we can obtain an alternative estimate for predicting the range of the warm cascade.Let us rewrite the ODE- η as ∂ ωω n = 1 n (cid:20) ( ∂ ω n ) − | η | S ω (cid:21) . (29)Keeping in mind that the MB distribution is always characterized by a positive concavity,i.e. ∂ ωω n ≥
0, and considering the hypothesis ∂ ω n (cid:39) ∂ ω n MB we find | η | = A S ω e − ωT T = g η ( ω, A, T ) . (30)Similarly to what we have done for the inverse cascade, we now can suggest that the changeof concavity occurs near ω min . This results in | η | = A S ω min e − ω min T T . (31)However, we do not expect a good prediction as before because in this case the exponentialterm is not a rapidly varying function near ω min . C. Double cascade
We have now all tools to study the double cascade process. Let us force at ω f , dissipateat ω max and ω min , and consider the case ω min (cid:28) ω f (cid:28) ω max . If the forcing range is narrow,the simple relation (cid:15) = η ω f holds for the fluxes. Using this relation, and combining (20)and (28), we can estimate T and A in the system: T = 2 ω max72 ln ω max ω min + ln ω max ω f − ,A = 29 (cid:115) | η | S ω / , (32)and, therefore, the chemical potential µ = T (cid:32)
12 ln
S ω / | η | + ln 92 (cid:33) . (33)Note that the temperature appears to be independent of the fluxes and is completely con-trolled by the forcing and the dissipation scales. This means that increasing the forcingstrength without moving ω f simply adds more particles into the system with the energy perparticle remaining the same. 20 V. NUMERICAL RESULTS
In this Section we present the numerical results obtained by using the DAM and byintegrating, at lower resolution, the HIBE. Our aim is to compare results for the warmcascade solutions of DAM, which has been devised as a local approximation of the integralcollision operator, with direct numerical simulation of the full integro-differential equation(9).
A. DAM resutls
We will first present some numerical experiments on integration of the Cauchy problems(16) and (25) in which we take for simplicity S = 1. Note that all numerical simulationscan be performed without any loss of generality starting with a particular value ω becauseof re-scaling properties described in Appendix E.
1. Constant direct energy cascade
In Fig. 7 we show the results obtained by integrating equation (16) with ω = 3 . (cid:15) . As initial conditions, we choose the values of the spectrum n and its slope n (cid:48) from the MB distribution having A = 1 and T = 1. The solutions followthe thermodynamic solution (shown as a continuous line) until they rapidly deviate and reachthe front in the vicinity of particular values of ω max . This numerical experiment exhibits twoimportant facts always observed in simulations performed with different initial conditions:the presence of a long transient in which the flux correction is negligible with respect to thethermodynamic MB distribution and the presence of a particular value ω max at which n ( ω )goes to zero. A lin-log plot of the function g (cid:15) ( ω, , (cid:15) = 1, (cid:15) = 10 − and (cid:15) = 10 − marksthe predicted cut-off frequencies for the respective flux values. Agreement with the behaviorin Fig. 7 is evident: the values of ω max obtained with equation (22) and Fig. 8 coincidewith the observed values in Fig. 7 within 5%. Note that the peak of g (cid:15) ( ω, , , ) is around ω = 3 .
5: this is why we set this value as initial condition ω .21 -8 -4
4 6 8 10 12 14 16 18 20 n ( ω ) ω ε =10 ε =10 -2 ε =10 -4 A e - ω /T FIG. 7: DAM simulations of (16) for different constant energy flux (cid:15) starting with the same initialcondition at ω = 3 . n MB ( ω ) = Ae − ωT where A = 1 and T = 1(continuos line). -4 -2
4 6 8 10 12 14 16 ε ω g ε ( ω ) FIG. 8: Plot in lin-log scale of the function g (cid:15) ( ω, , In Fig. 9 we present the results for a particular case with flux (cid:15) = 1. We can appreciate thepresence of warm cascade and the front solution near ω max . The linear behavior of the front isevident in the zoom near ω max showed in the inset. Numerically we are able to measure ω max and so evaluate B from equation (17). The theoretical prediction agrees with the measured22 -10 -8 -6 -4 -2 n ( ω ) ω A e - ω /T -8 ω max =5.751B err =0.997% FIG. 9: DAM simulation of ODE- (cid:15) (16) with constant energy flux (cid:15) = 1 (dashed line). The initialconditions in ω = 3 . T = 1 and A = 1 (continuos line). Theinset shows a zoom of numerical n ( ω ) (dots) in the vicinity of the point ω max where a linear fit isshown by continuos line. slope with the error B err = 0 . B err = | B meas − B est | /B meas where B meas is the measured linear coefficient and B est is the one taken form relation (17).In all other simulations performed with different values of (cid:15) or different initial conditions, B err is always within 5%.We now check numerically the validity of the matching prediction (20) by taking differentinitial condition n MB ( ω = 3 .
5) varying T and keeping A = 1 and (cid:15) = 1: results are plottedin Fig. 10. It is evident from the figure that the predicted temperature (continuous blackline) evaluated from relation (20) is an overestimation of the numerical results (dots) andthe error is around 10%. Finally prediction for the alternative temperature relation (23) isplotted with gray dashed line: it appears to give a better estimation than relation (20).
2. Constant particle cascade
We now investigate the inverse particle cascade by solving Cauchy problem (25) going ω -backward. In Fig. 11 we show numerical results obtained by taking initial conditions at ω from MB distribution n MB ( ω ) = Ae − ωT with T = 1, A = 1. As in the case of constant energy23 T ( ! m a x ) ! max T KK T g FIG. 10: Checking of predictions in DAM constant energy flux cascade (16): the points representthe temperature of the initial condition T with respect to the measured ω max . Solid line is thematching relation (20) while dashed one is obtained from (23). flux, here the warm cascade range is wider for smaller flux values. We also observe fronts invicinities of cutoff points ω min . In Fig. 12 we show the function g η ( ω, ,
1) which representsthe prediction of the thermodynamic range (30). Qualitative front values of results in Fig.11 show poor agreement with this na¨ıve estimation.The front solution is analysed in detail in Fig. 13 where we choose the particular casewith η = −
1. The linear behavior is demonstrated in the inset. Moreover a numericalestimation of ω min lets us evaluate B , see equation (26). The error B err is presented in thefigure; for all other simulations we have performed B err remained within 4%.Finally we check KK matching prediction for the thermodynamic quantity A with re-spect to ω min presented in equation (28): results are showed in Fig. 14. In this case the KK analytical prediction (continuous line) underestimates the numerical data while theestimation (31) is completely out of range (dashed line). However the scaling A ∼ ω − / ofKK prediction tends to be reached for small values of ω min , where ω min (cid:28) T .24 -6 -3 -2 -1 n ( ω ) ω η =-10 η =-10 -2 η =-10 -4 A e - ω /T FIG. 11: DAM simulations of (25) for different constant inverse particle flux η starting with thesame initial condition at ω = 3 . n MB ( ω ) = Ae − ωT where A = 1 and T = 1 (plotted with continuos line). -4 -2 -2 -1 | η | ω g η ( ω ) FIG. 12: Plot of the function g η ( ω, ,
3. Double cascade
An example of double cascade is presented in Fig. 15 where we set the forcing at ω f = ω = 3 .
5. We show here three cases where the particle fluxes are respectively η = − η = − − and η = − − . Measuring ω min and ω max for each case we are able to estimate the25 -6 -3 n ( ω ) ω A e - ω /T -5 ω min =1.942B err =0.839% FIG. 13: DAM simulation of (25) with η = − ω = 3 . n MB ( ω ) = Ae − ωT where A = 1 and T = 1 (plotted withcontinuos line). Inset: lin-lin scale zoom in the vicinity ω min (dots) where the best linear fit ispresented with a continuous line. temperature T est from prediction (32). Results do not agree with the expected temperature(the initial conditions set it at T = 1) but they approach this value for bigger ranges,i.e. when the condition ω min (cid:28) ω f (cid:28) ω max is better satisfied (see for example the case η = − − ). B. HIBE results
We now to present results of the direct simulation of HIBE with the full Boltzmann colli-sion integral and compare them with predictions obtained by DAM. As we have mentionedabove, the evaluation of (4) is numerically challenging and it is nowadays practically impos-sible to simulate such wide ω -space ranges as we have done using the DAM. In the presentwork, we will always use a low resolution of 101 points by considering ω ∈ [0 , ω cutoff ] andtaking a uniform distribution with ∆ ω = ω cutoff . We have checked that the numerical solu-tions are mesh independent by taking a finer mesh, 201 points, and comparing the solutionof one critical case.The δ -function in (4) defines a resonant manifold over which the integrand need to beevaluated; numerically it is a set of discrete resonant conditions M = { ω , ω , ω , ω } / ω +26 A ( ω m i n ) ω min A KK A g -1 -2 -1 FIG. 14: Checking predictions in DAM constant inverse flux cascade (25): the points represent thethermodynamic quantity A with respect to the measured ω min . Continuos line is the KK matchingprediction given in equation (28) while dashed one is obtained from (31). ω = ω + ω which can be pre-computed. Note that the dissipation at high wave numbers ischosen to satisfy ω max ≤ ω cutoff /
1. Direct cascade study
We first analyze the direct energy cascade by putting the forcing scale near the low- ω dissipation scale in order to have a wider direct inertial range. Numerical results forthese final steady states were previously presented as examples in Fig. 2 and Fig. 4.We concentrate now only on the last one: here we kept fixed ω min = 5, ω f = 21 and ω max = 95 and varied the forcing coefficient, i.e. the fluxes η and (cid:15) . We were claimingthat the temperature of the systems is the same because qualitatively the distributionshave identical slopes. Moreover we observed in all the examples that left and right branchchemical potentials and temperatures can be defined by the forcing scale.With these previous DAM results in mind we have measured A and T in three examplespresented in Fig. 4: the results are shown in Fig. 16 and are compared to analytical27 -10 -6 -2 -2 -1 n ( ! ) ! " =10 T est =18.5 " =10 -2 T est =1.5 " =10 -4 T est =1.3 energyparticles FIG. 15: DAM double cascade simulations of equations (16) and (25) for three different values ofthe particle flux η (and consequently of the energy flux (cid:15) = η ω f ). The initial condition are takenat ω f = 3 . T = 1 and A = 1. Measuring ω min and ω max in eachcase we estimate of the temperature T est from prediction (32). predictions (32). As expected the quantity A ∼ √ η but the line (in log-log plot) is shiftedwith respect to the interval between A L and A R , represented respectively with filled andempty circles. However, the theoretical prediction is much closer to A R , which is naturalbecause the right inertial interval is wider than the left one. In fact, the agreement of A R with the theory is quite good considering the presence of the undefined constant S inthe theoretical prediction. The temperature is shown in Fig. 17: even though T L and T R are different they both appear to be forcing independent, as predicted. The temperatureevaluated from relation (32): temperature (dashed line) stands in between of these values,and closer to T R , which, again, is natural because the right inertial interval is wider.We have also analyzed sensitivity of the temperature to varying the high- ω dissipationrange and results are presented in Fig. 18. Keeping the forcing constant and changing thevalue of ω max the system reaches steady states characterized by different temperatures T L (filled circles) and T R (empty circles). The prediction (32), shown by the continuous line, isin between of the two temperatures and is closer to T R - again due to the wider right range.28 -4 -3 -2 -4 -3 -2 A ( | η | ) | η | FIG. 16: Results for the fitted values of the thermodynamic quantity A = e − µT plotted against theforcing levels, as obtained in the simulations shown in Fig. 4: the values of A L are shown by filledcircles while A R - by empty ones, the blue line refers to prediction (32) with S = 1. -3 -2 -1 T ( ε ) ε FIG. 17: Measured temperatures T L (filled triangles) and T R (empty ones) obtained in the simu-lations shown in Fig. 4. The dashed line is the analytical prediction (32).
2. Inverse cascade study
Finally, we have performed some simulations putting the forcing scale near the dissipationat high ω ’s in order to study the inverse cascade process. In this case too, as reported in Fig.3, we observe two different values of thermodynamic quantities on the left and on the right29
15 30 45 60 95 130 165 200 T ( ! m a x ) ! max FIG. 18: Temperature in different steady state keeping the forcing constant and varying the dissi-pation scale ω max : big empty circles correspond to the temperature T R on the right of the forcingscale, whereas small filled circles to the left side, T L . The continuous line is the prediction (32). T ( ω m i n ) ω min FIG. 19: Temperature in different steady state keeping the forcing constant and varying the dissi-pation scale ω min : big empty circles correspond to the temperature T R on the right of the forcingscale, whereas small filled circles to the left side, T L . The continuous line is the prediction (32). from the forcing. Here we are able to study the scaling of the thermodynamic quantities T and A with respect to changes of the small- ω dissipation scale ω min . Results for T areshown in Fig. 19 and for A in Fig. 20, with the “left” quantities shown by filled circles30 -6 -4 -2
1 10 100 A ( ω m i n ) ω min FIG. 20: Thermodynamic amplitude A for different ω min and fixed forcing in log-log scales innumerical simulations of HIBE. Filled circles corresponds to A L while empty circles - A R . Thecontinuos line is formula (32) with S = 1. and the “right” ones by empty circles. There is a reasonably good agreement of T with theprediction (32) for small ω min . This is natural because smaller ω min corresponds to largerinverse cascade inertial range and also because the prediction is valid when ω min (cid:28) T .On the other hand, for A the prediction (32) is in better agreement with the data at large ω min with ω min ∼ ω f . This is due to two possible reasons. First, we underline that theagreement can be made more suitable since the analytical prediction contains the undefinedorder-one parameter S which could be adjusted to better fit the numerical results. Second,the particle flux which defines A in relation (32) can be smaller due to finite range effects.Indeed, following [32], the ratio of the leftward particle flux to the total particle productionrate is estimated in η L = η ω max − ω f ω max − ω min . (34)This equation, in addition to other corresponding to rightward fluxes in the cited paper,states that for the particle flux to be mostly to the left inertial ranges in both directionsmust be large (note that this is also the condition of validity of the Fjørtoft argument).Fig. 21 shows the behavior of the normalised left measured flux η L /η (empty triangles) withrespect to ω min . We can clearly see that the measured particle flux is indeed much smaller31 .00.51.0 1 10 100 η L ( ω m i n ) / η ω min FIG. 21: The plot shows in log-linear scale the ratio of the measured inverse particle flux over thetotal one η L /η (empty triangles) with respect to ω min ; the continuous line correspond to the totalflux while the dashed one follows the finite range correction (34). than the one imposed by the forcing term (continuos line), around one third of it. This is inquite good agreement with the finite range prediction (34) plotted with dashed line. Similarreasoning can be made for corrections on A ( | η | ) in the case of direct cascade example in Fig.16. V. CONCLUSIONS
In the present paper we investigated stationary turbulent states in the isotropic Boltz-mann kinetic equation for hard spheres. This was done by looking for steady nonequilibriumstates in open systems, that is when forcing and dissipation mechanisms are present. Analo-gies with similar results of wave turbulence theory suggest the manifestation of a warmcascade , i.e. a constant direct flux of energy and inverse flux of particles on backgroundof thermodynamic Maxwell-Boltzmann distribution. This is a consequence of wrong fluxdirections in KZ solutions with respect to the Fjørtoft argument.We have built an ad-hoc differential approximation model to easily simulate the cascadeprocesses. Indeed, this simplification allowed us to reach a wide range of scales inaccessible bysolving the isotropic Boltzmann kinetic equation directly. Simulations show the presence ofa warm cascade with approximately the MB shape followed by a sharp front for both energy32nd particle cascades. We have physically interpreted ω min and ω max as intrinsic dissipationscales at low and high ω ’s which are necessary to establish the steady state. Moreover,we have found analytical predictions relating the particle and energy fluxes, forcing anddissipations scales to the thermodynamic quantities of the system. In particular we haveshown that the temperature is independent of the amplitude of the fluxes but only dependson the forcing and dissipation scales.We have then compared the theoretical predictions and the numerical results obtainedwith the differential approximation model with simulations of the complete isotropic Boltz-mann kinetic equation. Even though the resolution for the latter was limited by the availablecomputational power, the results are comparable and in good agreement with the analyticalpredictions. In particular we have verified that the steady state is characterized by a warmcascade where a fitted thermodynamic Maxwell-Boltzmann distribution has been used tomeasure temperature and chemical potential of the system. We observe, in agreement withour analytical predictions, that the temperature is completely defined by the forcing anddissipation scales and does not depend on the fluxes.We hope that this work may open some perspectives towards understanding nonequi-librium steady states and their net currents (fluxes) by cross-fertilization with the weakturbulence theory. Acknowledgments
We would like to thank Guido Boffetta, Colm Connaughton, Filippo De Lillo, StefanoMusacchio, Al Osborne, and Arturo Viero for fruitful discussions. Simulations were per-formed on computational resources founded by the Office of Naval Research (ONR). Finally,we are grateful to the Gnu Scientific Library (GSL) developers for providing free softwarewhich has been used for simulations.
Appendix A: Three-dimensional δ -function angular average The angular average of the four-wave linear momentum conservation δ ( k ) = δ ( k + k − k − k ) is evaluated by splitting it into two δ -functions of three particle collision. This33esults in (cid:90) Ω δ ( k ) d Ω = (cid:90) Ω (cid:90) k max k min δ ( k + k − k ) δ ( k + k − k ) d k d Ω = (cid:90) k max k min (cid:20)(cid:90) Ω δ ( k + k − k ) d Ω (cid:21) (cid:20)(cid:90) Ω δ ( k + k − k ) d Ω (cid:21) k d − dkd Ω= 4 π (cid:90) k max k min kk k kk k k dk = 2 πk k k k min( k , k , k , k ) , (A1)where geometrically k min = | k − k | = | k − k | and k max = k + k = k + k . For detailsabout the integration of three particle δ -function see Appendices in [5]. Appendix B: Kolmogorov-Zakharov solutions for general HIBE
The Boltzmann collision integral I coll is defined as I coll ( x , k , t ) = (cid:90) ∞−∞ Γ [ n ( x , k , t ) n ( x , k , t ) − n ( x , k , t ) n ( x , k , t )] × δ ( k + k − k − k ) δ ( | k | + | k | − | k | − | k | ) d k , (B1)where the two δ -functions assure the conservation of the linear momentum and kinetic energy.In the isotropic case it is convenient to move in the energy domain ω i = | k i | ∈ [0 , + ∞ ) andso the HIBE results in I ( ω ) = (cid:90) ∞ S ( n n − n n ) δ ( ω ) dω , (B2)where I ( ω ) = Ω I coll ( x , ω , t ) ω d − (cid:12)(cid:12)(cid:12) dk dω (cid:12)(cid:12)(cid:12) and we use for brevity n i = n ( ω i ) = n ( x , | k i | , t ),and δ ( ω ) = δ ( ω + ω − ω − ω ). The functional S is S = 116 ( ω ω ω ω ) d − (cid:104) Γ δ ( k + k − k − k ) (cid:105) Ω (B3)and the operator (cid:104)·(cid:105) Ω states for the integration over solid angles. It is important for thefollowing to estimate the homogeneity degree of S . Supposing that the collisional kernelscales as Γ λ (34) λ (12) = λ β Γ , we have S λ (34) λ (12) = λ ( d − ) +2 β − d S = λ d +2 β − S . (B4)Moreover its behavior at the boundaries of integration islim ω i → + ∞ S ∼ ω d − τ i lim ω i → + S ∼ ω d − τ i (B5)34f we assume that lim ω i → + ∞ (cid:104) Γ δ ( k + k − k − k ) (cid:105) Ω ∼ ω τ i (B6)and lim ω i → + (cid:104) Γ δ ( k + k − k − k ) (cid:105) Ω ∼ ω τ i (B7)(note that for ω i → ∞ also another ω j must go to infinity due to the δ -function).In the following we will suppose that the particle distribution function follows the power-law distribution n ( ω ) = A ω − ν and so I ( ω ) = A (cid:90) ∞ S − (cid:2) ω − ν ω − ν − ω − ν ( ω + ω − ω ) − ν (cid:3) Θ( ω + ω − ω ) dω , (B8)where Θ is the Heaviside step function.
1. Kolmogorov-Zakharov solutions
We will present the Kolmogorov-Zakharov solutions of the collision integral using themethod presented by Balk in [22]. The collision integral, without any loss of generality, canbe rewritten as I ( ω ) = A ω − − µ (cid:90) ∞ S ( ω − ν ω − ν − ω − ν ω − ν ) ( ω ω ω ω ) ω µ δ ( ω ) dω ω dω ω dω ω (B9)where the exponent µ = 2 ν + 1 − β − d dω i ω i ). If the integral converges, Balk proved that is possible to interchange thethree integration index in the integrand with the fourth one, ω . Thanks to the symmetricproperties of the collision kernel we can write I ( ω ) = A ω − − µ (cid:90) ∞ S ( ω − ν ω − ν − ω − ν ω − ν ) ( ω ω ω ω ) × ( ω µ + ω µ − ω µ − ω µ ) δ ( ω ) dω ω dω ω dω ω , (B11)which clearly vanishes for µ = 0 or µ = 1. This corresponds to the condition on the exponent ν = ν | µ =0 = 3 d −
24 + βν = ν | µ =1 = 3 d β. (B12)Note that first KZ solution for HIBE were presented in [17].35 . Convergence of the integral (locality condition) The locality of interactions is guaranteed by the convergence of the collision integral. Wethen investigate the possible values of ν which assure the convergence around the integrandsingularities. a. Limit ω → ∞ In the limit of ω → ∞ we can approximate ( ω + ω − ω ) − ν = ω − ν − νω − ν − ( ω − ω ) + O ( ω − ν − ) at the second order. The argument in the square brackets of (B8) results in[ ... ] (cid:39) ω − ν (cid:2) ω − ν − ω − ν + νω − ν ω − ( ω − ω ) (cid:3) . (B13)As a consequence, when ν >
0, the integrand for large ω goes like ω − ν − ω − ν ω ν − d +2 − τ and so theconvergence condition is ν > d − τ . (B14) b. Limit ω → + In the limit of ω → + we can approximate ( ω + ω − ω ) − ν = ( ω − ω ) − ν − νω ( ω − ω ) − ν − + O ( ω ) at the second order. The argument in the square brackets of (B8) resultsin [ ... ] = ω − ν (cid:2) ω − ν − ω − ν ω ν ( ω − ω ) − ν + νω − ν ω ν +13 ( ω − ω ) − ν − (cid:3) . (B15)So, when ν >
0, the integrand for small ω goes like ω − ν ω ν − d − τ and so the convergencecondition is ν < d τ . (B16)Analogue condition holds for the singularity ( ω + ω − ω ) − ν → + .
3. Constant fluxes
The solutions n ( ω ) = A ω − ν and n ( ω ) = A ω − ν correspond, respectively, to constantflux of particle and energy. To demonstrate this fact we perform the substitution ω i = ω ξ i i (cid:54) = 1 in the equation (B11) which results, recalling the homogeneity of the integrandfunction, in I ( ω ) = A ω − − µ (cid:90) ∆ S ξ ξ ξ ( ξ − ν ξ − ν − ξ − ν ) ( ξ ξ ξ ) × (1 + ξ µ − ξ µ − ξ µ ) δ (1 + ξ − ξ − ξ ) dξ = A ω − − µ U ( µ ) (B17)The integral U ( µ ) is now performed over the triangle ∆ in the ξ × ξ space satisfying theconditions 0 ≤ ξ i ≤ ξ ≥ − ξ , without any dependence on ω . a. Flux of particles The flux of particles is defined as Q ( ω ) = − (cid:90) ω I ( ω ) dω = − A U ( µ )4 (cid:90) ω ω − − µ dω = A U ( µ ) ω − µ µ . (B18)If µ = 1 the flux is zero while in the case µ = 0 it is indeterminate. By applying the Del’Hˆopital rule in the latter case we find Q ( ω ) | µ =0 = A (cid:90) ∆ S ξ ξ ξ ( ξ − ν ξ − ν − ξ − ν ) ( ξ ξ ξ ) ln (cid:18) ξ ξ ξ (cid:19) δ (1 + ξ − ξ − ξ ) dξ (B19)The integrand, and so the sign of the particle flux, is always negative for ν >
0. This isclear by looking at the sign of every factors in the integral: all are trivially positive except( ξ − ν ξ − ν − ξ − ν ) and ln (cid:16) ξ ξ ξ (cid:17) . Recalling that (1 − ξ )(1 − ξ ) ≥ ξ = ξ + ξ − ≤ (1 − ξ )(1 − ξ ) = ξ ξ − ξ − ξ + 1 = ξ ξ − ξ = ⇒ ξ ξ ≥ ξ , (B20)which leads to ln (cid:16) ξ ξ ξ (cid:17) ≤ ξ − ν ξ − ν − ξ − ν ) ≤ ν ). As a consequence Q ( ω ) ≥
0, that is the particle flux goes from low to high frequencies. b. Flux of energy
The flux of energy is P ( ω ) = − (cid:90) ω I ( ω ) ω dω = − A U ( µ )4 (cid:90) ω ω − µ dω = − A U ( µ ) ω − µ − µ ) (B21)37nd is null when µ = 0 while indeterminate in the case µ = 1. Again applying the Del’Hˆopital rule we have P ( ω ) | µ =1 = A (cid:90) ∆ S ξ ξ ξ ( ξ − ν ξ − ν − ξ − ν ) ( ξ ξ ξ ) (B22) × [ ξ ln( ξ ) − ξ ln( ξ ) − ξ ln( ξ )] δ (1 + ξ − ξ − ξ ) dξ As previously discussed, the term ( ξ − ν ξ − ν − ξ − ν ) ≤ ν >
0. Differently, thefactor [ ξ ln( ξ ) − ξ ln( ξ ) − ξ ln( ξ )] is always positive but here the demonstration is notso trivial as in the previous case and for a complete discussion see [18]. So P ( ω ) ≤
0, whichmeans that the energy flux goes from high to low frequencies.
Appendix C: Matching Kats-Kontorovich to ω max front solution We will here find the match between the KK correction and the front solution for theODE- (cid:15) . We make the hypothesis that the front occurs for ω max (cid:29) T and so it is reasonableto think that the term ω − in equation (19) it is slowly varying with respect to e ωT . So byintegrating twice in ω (19) and match to the front we get the Cauchy problem (cid:15)ω − T A − e ωT = − S ˜ n + c ( ω − ω max ) + c ˜ n ( ω max ) = − ∂ ω ˜ n ( ω max ) = BAe ω max T (C1)where the first condition assures that n ( ω max ) = 0 and the second that the front behavior islinear with slope B = − (cid:15) ω − max found in equation (17). The integration constants are then c = − S + (cid:15) ω − max T A − e ω max T c = S − (cid:15) ω − max A − e ω max T − (cid:15) ω − max T A − e ω max T . (C2)We will now match this solution in the regime where T (cid:28) ω max and the flux is negligiblewith respect to the thermodynamic solution. In this regime, by assuming the scaling relation (cid:15) ∼ ω max A e − ω max T , the coefficients results in c (cid:39) − Sc (cid:39) S − (cid:15) ω − max A − e ω max T (C3)and so the smallness of the correction reads as˜ n ( ω ) = 0 = − (cid:15)ω − T A − e ωT + ( ω − ω max ) S − (cid:15) ω − max A − e ω max T − S. (C4)38inally, considering that ω (cid:28) ω max , we recover and validate the relation (cid:15) = S ω max A e − ω max T . (C5) Appendix D: Matching Kats-Kontorovich to ω min front solution The KK correction for that case is given by equation (27). We will consider the limit ω (cid:28) T and so e ωT (cid:39)
1. By integrating twice in ω we get the Cauchy problem | η | ω − A − = − S ˜ n + c ( ω − ω min ) + c ˜ n ( ω min ) = − ∂ ω ˜ n ( ω min ) = | η | S ω − min A − (D1)with the initial conditions chosen in order to match with the front solution. The integrationconstants result in c = − S + | η | ω − min 463 A − c = − S | η | ω − min A − + | η | ω − min 29 A − . (D2)We assume and guess that the particles flux scales as | η | ∼ ω min A . Now, in the regime ω (cid:29) ω min were the correction is negligible we have˜ n ( ω ) = 0 = −| η | ω − A − + c ( ω − ω min ) + c (cid:39) ω c . (D3)So, finally, we impose that c = 0 to get the condition on the flux | η | = S (cid:18) (cid:19) A ω min . (D4) Appendix E: Scaling properties of DAM
The constant energy flux DAM (15) and the constant particles flux DAM (24) can begenerally written as c = − Sω p n ( ω ) ∂ ωω log n ( ω ) (E1)where the exponent is respectively p = 13 / p = 11 / c is a constant that representthe considered flux. Lets now analyze the rescaling properties of that equation by thefollowing change of variables c = λ α ¯ cω = λ β ¯ ωn = λ γ ¯ n. (E2)39fter some easy algebra we find that the system is invariant if α = ( p − β + 2 γ. (E3)As a consequence we can establish how the thermodynamic quantities defined by theMaxwell-Boltzmann distribution n MB ( ω ) = Ae − ωT = e − ω + µT vary: the temperature T scalesas ω and so T = λ β ¯ T , while the chemical potential µ scales as µ = λ β γ ¯ µ . [1] C. Bustamante, J. Liphardt, and F. Ritort, Physics Today , 43 (2005), URL http://link.aip.org/link/?PTO/58/43/1 .[2] E. H. Lieb, Physica A: Statistical Mechanics and its Applications , 491 (1999),ISSN 0378-4371, proceedings of the 20th IUPAP International Conference on StatisticalPhysics, URL .[3] T. S. Komatsu, N. Nakagawa, S.-i. Sasa, and H. Tasaki, Phys. Rev. Lett. , 230602 (2008).[4] C. Cercignani, The Boltzmann equation and its applications (Springer, 1988).[5] V. Zakharov, V. L’vov, and Falkovich,
Kolmogorov Spectra of Turbulence 1: Wave Turbulence (Springer-Verlag, 1992).[6] A. Kolmogorov, in
Proceedings (Doklady) Academy of Sciences, USSR (1941), vol. 30, pp.301–305.[7] U. Frisch,
Turbulence: the legacy of AN Kolmogorov (Cambridge University Press, 1995).[8] P. Janssen,
The Interaction of Ocean Waves and Wind (Cambridge University Press, 2004).[9] M. Onorato, A. Osborne, M. Serio, D. Resio, A. Pushkarev, V. Zakharov, and C. Brandini,Physical Review Letters , 144501 (2002).[10] A. Dyachenko, A. Korotkevich, and V. Zakharov, Journal of Experimental and TheoreticalPhysics Letters , 546 (2003).[11] Y. Lvov, K. Polzin, and E. Tabak, Physical Review Letters , 128501 (2004).[12] S. Dyachenko, A. C. Newell, A. Pushkarev, and V. E. Zakharov, Physica D: Nonlin-ear Phenomena , 96 (1992), URL .
13] N. Berloff and B. Svistunov, Physical Review A , 13603 (2002).[14] S. Nazarenko and M. Onorato, Physica D: Nonlinear Phenomena , 1 (2006).[15] D. Proment, S. Nazarenko, and M. Onorato, Physical Review A (Atomic, Molecular, andOptical Physics) , 051603 (pages 4) (2009), URL http://link.aps.org/abstract/PRA/v80/e051603 .[16] S. Galtier, S. Nazarenko, A. Newell, and A. Pouquet, Journal of Plasma Physics , 447(2000).[17] A. Kats, V. Kontorovich, S. Moiseev, and V. Novikov, ZhETF Pis ma Redaktsiiu , 13(1975).[18] A. Kats, Soviet Journal of Experimental and Theoretical Physics , 1106 (1976).[19] C. Connaughton and S. Nazarenko, Phys. Rev. Lett. , 044501 (2004).[20] V. Karas, S. Moiseev, and V. Novikov, Zhurnal Eksperimental’noi i Teoreticheskoi Fiziki ,1421 (1976).[21] C. Connaughton, S. Nazarenko, and A. C. Newell, Physica D: Nonlinear Phenomena ,86 (2003), URL .[22] A. Balk, Physica D: Nonlinear Phenomena , 137 (2000).[23] C. Connaughton, Physica D: Nonlinear Phenomena , 2282 (2009), ISSN0167-2789, URL .[24] S. Hasselmann, K. Hasselmann, J. Allender, and T. Barnett, J. Phys. Oceanogr , 1378(1985).[25] C. E. Leith, Physics of Fluids , 1409 (1967), URL http://link.aip.org/link/?PFL/10/1409/1 .[26] V. L’vov and S. Nazarenko, JETP letters , 541 (2006).[27] G. Boffetta, A. Celani, D. Dezzani, J. Laurie, and S. Nazarenko, Journal of Low TemperaturePhysics , 193 (2009).[28] J. Peacock, Cosmological physics (Cambridge Univ Pr, 1999).[29] Y. Lvov, R. Binder, and A. Newell, Physica D: Nonlinear Phenomena , 317 (1998), ISSN0167-2789.[30] V. Zakharov and A. Pushkarev, Nonlinear Processes in Geophysics , 1 (1999).
31] P. Asinari, Computer Physics Communications , 1776 (2010), ISSN 0010-4655, URL .[32] A. C. Newell, S. Nazarenko, and L. Biven, Physica D: Nonlinear Phenomena ,520 (2001), ISSN 0167-2789, URL .[33] This theorem, originally put forward by Fjørtoft in 1953 for the 2D turbulence, says that thatthe integral whose density grows fastest with the wavenumber/momentum must cascade fromlow to high wavenumbers/momenta. The other integral must cascade inversely, from high tolow wavenumbers/momenta. For the classical particles, this means that the energy flux mustbe from low to high momenta, and the flux of particles must be toward low momenta; seeSection II C.[34] In Navier-Stokes the warm cascades correspond to so called bottleneck phenomenon whicharises in numerics due to an energy flux stagnation near the maximum wave-number..[33] This theorem, originally put forward by Fjørtoft in 1953 for the 2D turbulence, says that thatthe integral whose density grows fastest with the wavenumber/momentum must cascade fromlow to high wavenumbers/momenta. The other integral must cascade inversely, from high tolow wavenumbers/momenta. For the classical particles, this means that the energy flux mustbe from low to high momenta, and the flux of particles must be toward low momenta; seeSection II C.[34] In Navier-Stokes the warm cascades correspond to so called bottleneck phenomenon whicharises in numerics due to an energy flux stagnation near the maximum wave-number.