{'fraction_non_alphanumeric': 0.03185136031851361, 'fraction_numerical': 0.009584900095849001, 'mean_word_length': 4.442750133761369, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Nonequilibrium steady solutions of the Boltzmann equation

We report a study of the homogeneous isotropic Boltzmann equation for an open system. We seek for nonequilibrium steady solutions in presence of forcing and dissipation. Using the language of weak turbulence theory, we analyze the possibility to observe Kolmogorov-Zakharov steady distributions. We derive a differential approximation model and we find that the expected nonequilibrium steady solutions have always the form of warm cascades. We propose an analytical prediction for 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 amplitude and determined only by the forcing and dissipation scales. Finally, we perform direct numerical simulations of the Boltzmann equation finding consistent results with our theoretical predictions.


I. 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 steady states have net currents (fluxes): examples of nonequilibrium steady-state systems include an object in contact with two thermal sources at different temperatures, for which the current is a heat flux; a resistor with electric current flowing across it; the kinesin-microtubule system, for which kinesin motion is the current. Most biological systems, including molecular machines and even whole cells, are in nonequilibrium states [1]. In particular, biological systems rely on a continuous flux of energy and/or particles supplied by some proper environmental reservoirs.
In statistical mechanics, investigating the general properties of a system in contact with reservoirs, namely an open system, is a long lasting problem (e.g. see the second problem discussed by E.H. Lieb on the occasion of the award of the Boltzmann medal [2]), even though 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 stationary state 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 theory of gases. In particular, we consider a system composed of a large number of interacting particles, comparable to the Avogadro number. The Boltzmann kinetic equation (BKE) describes the time evolution of the single-particle distribution function, which provides a statistical description of the positions and velocities (momenta) of the gas molecules. This integro-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 stosszahl ansatz [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 possible to 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 theory [5] where it has been shown that other solutions with respect to thermodynamic solutions can be stationary states of the system, in case of external forcing and dissipation. These distributions, which have usually the form of power-laws in momentum space, are called Kolmogorov-Zakharov (KZ) and they represent constant flux of conserved quantities similar to the Kolmogorov energy cascade in strong Navier-Stokes turbulence [6,7]. These solutions, 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 weakly nonlinear dispersive models: examples can be found in water waves [8][9][10], internal waves [11], nonlinear optics [12], Bose-Einstein condensation [13][14][15], magnetohydrodynamics [16].
An out of equilibrium description of the Boltzmann equation using the KZ solutions was first 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. In particular in [18] Kats showed that for all realistic physical situations the direction of the cascades in the system is always in the wrong orientation with respect to the one predicted by the Fjørtoft theorem [33]. When a formal KZ solution has a flux direction contradicting with the Fjørtoft theorem, this spectrum (even if local) cannot be established because it cannot 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 twodimensional nonlinear Schrödinger equation model the authors argued that in this case the KZ solution is not achievable and a mixed state, with both a cascade and a thermodynamic components were proposed. Another example of mixed cascade-thermodynamic states can be found in the context of three-dimensional Navier-Stokes turbulence [19], where such mixed states were called warm cascades [34].
The present manuscript will focus on warm cascades found in the homogenous isotropic Boltzmann equation (HIBE) and in particular it will answer to the following important questions.
• What is precisely the relation between the conserved quantity fluxes and the thermodynamics 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 homogeneous isotropic Boltzmann equation with forcing and dissipation. We will then use a diffusion approximation model (DAM) to derive analytical predictions on how the thermodynamic quantities, temperature and chemical potential, are related to fluxes, forcing and dissipative scales. We will then test these predictions by numerically simulating both DAM and the complete homogenous isotropic Boltzmann equation.
The work is organized as follows: in Section II we review the properties of the Boltzmann equation for the homogeneous isotropic case; in Section III we introduce DAM and we derive the 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 calculations of 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 distribution function, which provides a statistical description for the positions and momenta of the gas molecules: the function n(x, k, t) express a probability density function in the one-particle phase space R d x × R d k with respect to time, where d is the dimension. Note that we denote the momentum variable with the letter k instead of the conventional p to follow the common notation of wave turbulence [5]. The Boltzmann equation takes the following form: where 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 and corresponds to Eq. (4.18), page 64 in Cercignani's book [4]. Here W describes synthetically the scattering amplitude transition 2 → 2 as a function of the momenta of the interacting particles. As we consider elastic collisions, the general way to express W is where δ-functions assure conservation of the total momentum and the total kinetic energy (which is proportional to |k| 2 ) of incoming and outgoing particles. The collision probability, expressed by Γ 34 12 In the present paper we will consider the case of three-dimensional rigid spheres with diameters σ and mass m, for which Γ simply results in Γ 34 12 = 2σ 2 /m [4]. For other interaction potentials, as Coulomb or Born approximation, refer to [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 momentum dependency coming only via the modulus k = |k|, so n(x, k, t) → n(k, t). It is useful to express the distributions in the energy space ω i = |k i | 2 where we use again the notation ω for the energy in analogy with wave turbulence. Then, the particle density in ω-space satisfies the relation N (ω, t)dω = n(k, t)dk or, in the other words, N (ω, t) = n(ω, t) Ω ω where Ω is the solid angle. After these considerations Boltzmann equation (1) where we denote for brevity N i = N (ω i , t) and n i = n(ω i , t), and the functional takes into account the change of coordinates and the average over solid angles. Hereafter, we always consider a three-dimensional gas of hard-sphere particles in a nondimensional form with m = 1 and σ 2 = 8 . Then the functional simply results in √ ω 4 (see Appendix A for details of the angular integration).
The HIBE has two conserved quantities, the mass and energy densities, Note that ρ M and ρ E are always constant in time for any distribution n and interaction potential, due to the fact that collisions are 2 → 2 and elastic. This is evident by evaluating their time derivatives using equation (4): the symmetries with respect to the integration indices immediately show that these quantities are zero.
A. Steady solutions

Equilibrium in a closed system
The HIBE (4) is an integro-differential equation with no general analytic solution. It is 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 the thermodynamic equilibrium described by the Maxwell-Boltzmann (MB) distribution, where A = e − µ T and constants µ and T have the meaning of the chemical potential and the temperature respectively (we consider the natural unit system, where the Boltzmann constant is one). Validation is trivial by plugging (7) into (4): for any value of T , µ and the interaction potential S 34 12 , the δ-function assures that the integrand is zero. Moreover, the total mass density of the system is ρ M = A (πT ) 3 2 , the total energy density is ρ E = 3 2 Aπ and any other moment of ω, due to the bi-parametric nature of the MB distribution, is a function of ρ M and ρ E . The H theorem states that in a closed system any out of equilibrium distribution with defined mass and energy densities will always relax to the MB distribution having same ρ M and ρ E .
In Fig. 1 we show a numerical simulation of the HIBE with initial condition given by a Gaussian function centered around a particular value of energy; as it is clear from the figure, the initial condition relaxes to the MB distribution. The numerical algorithm used to perform this simple example will be discussed in Section IV B. We can observe that the initial condition evolves reaching an equilibrium MB distribution: the exponential behavior become evident by observing the inset where we plot it lin-log plot scale. Moreover by fitting the results with the MB function we can find the thermodynamic quantities A and T : those correspond exactly to ones expected knowing initial mass and energy densities (note that now integrals (6) are evaluated from 0 to a finite value of ω due to numerical finiteness of ω-space). The inset shows the same plot in lin-log scale.

Nonequilibrium steady states
Now, what can we expect in an open system driven by external forcing and dissipation mechanisms? We will answer this question keeping in mind the main results of the wave turbulence theory. Part of this theory is dedicated to study steady solutions to kinetic equations in the power-law form, n(ω) ∼ ω −x , where the constant x assumes different values depending on the considered wave system. It is sometimes possible to find the so-called Kolmogorov-Zakharov (KZ) solutions n KZ (ω) ∼ ω −x which correspond to constant fluxes of conserved quantities through scales. The KZ distribution always appears in a range of scales, known as inertial range, between the forcing and dissipation were the source and sink are located.
As already mentioned, the HIBE conserves the number of particles and the energy, and so one could expect to observe two turbulent KZ cascades. The KZ exponent x can be evaluated 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 in Appendix B. The KZ exponents depend on the scaling behavior of the scattering term Γ 34 12 and on the dimension d of the particle system. For the particular case of three-dimensional hard spheres we have The simplest way to mimic an open system where steady nonequilibrium distributions of the form of turbulent KZ solutions can be establish is to consider a forced-dissipated HIBE The forcing F is constant in time and very narrow near a particular energy value ω f : with this choice the incoming fluxes of particles η and energy roughly satisfy relation = ω f η. The dissipation term D is implemented as a filter which removes, at each iteration time, energy and particles outside of the domain ω ∈ (ω min , ω max ). Further details on the numerical scheme are explained in Section IV B. What happens if we try to solve numerically such forced/damped integro-differential equation?
In Fig. 2 and Fig. 3 we plot the nonequilibrium steady states obtained with numerical simulations of the HIBE with forcing and dissipation; the initial conditions are characterized by n(ω, t = 0) = 0 . The parameters in the simulations are ω min = 5, ω max = 195, the forcing rate F = 10 −5 . In Fig. 2 forcing is located at ω f = 22 and in Fig. 3 at ω f = 182. No power-law distributions, and so no KZ solutions (8), are observed (note that both plots are in lin-log scales), but instead one can see weakly perturbed exponential curves. We can attempt 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 with respect 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 scales and change the forcing rate. Numerical results for final steady states evaluated for three different forcing amplitudes, F = 10 −4 , 10 −5 , 10 −6 , are presented in Fig. 4. The effect of increasing 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 is   the speed at which the system, initially empty, reaches its steady state. In fig. 5 we show the energy density evolution (same line styles corresponds to same systems).
After these preliminary numerical results, a lot of questions can be posed. Why no KZ 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 the relation between physical quantities such as fluxes, forcing and dissipation scales and the MB parameters? The aim of this manuscript is to provide explanations to such phenomena and 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 collision integral converges. In other words, non-locality means that the collision integral is not weighted scale by scale but most of the contributions come from the limits of integration corresponding to the ends of the inertial range. Physically, the non-locality is in contradiction with the assumption that the flux of the relevant conserved quantity in the inertial range is carried only by the nearest scales. Mathematically, locality guaranties that the KZ spectrum is a valid solution in an infinite inertial range, which is not guarantied a priori because Zakharov transformation is not an identity transformation and could, therefore, lead to spurious solutions.
For the HIBE case, locality depends on the particular interaction potential, which affects the scaling of Γ 34

The Fjørtoft argument
Consider an open system where forcing scale ω f is widely separated from a low-ω dissipation scale ω min and a high-ω dissipation frequency ω max , thus ω min ω f ω max . Because the energy density in the ω-space is different from the particle density by factor ω, the forcing rate of the energy is related with the forcing rate of the particles η as ∼ ω f η.
Suppose that some energy is dissipated at the low scale ω min at a rate comparable with the forcing rate . But then the particles would have to be dissipated at this scale at the rate proportional to /ω min ∼ η ω f /ω min η, which is impossible in steady state because the dissipation cannot exceed the forcing. Thus we conclude that in the steady state the energy must dissipate only at ω max . By a symmetric contradiction argument one can easily 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 an inverse cascade (negative flux direction).

Flux directions in the HIBE
It has been proved in [18], see also Appendix B for details, that fluxes of the KZ solutions for all types of the interaction coefficient Γ have always the wrong directions with respect to the Fjørtoft argument requirements (in the case x > 0). An alternative way for finding the sign of the fluxes is considering them for general (not necessarily steady) power-law spectra n(ω, t) ∼ ω −x and plotting them as functions of x for a fixed (ω, t), see Fig. 6. Three exponents x correspond to steady solutions of HIBE: the particle equipartition x eq = 0, the KZ particle cascade x η and the KZ energy cascade x . As shown in Appendix B, we know that = 0 on the particle cascade, η = 0 on the energy cascade, whereas in the equipartition both fluxes are zero, i.e. = η = 0. We also know that for large negative x (large positive slope) both fluxes must be negative, as such a steep unsteady spectrum would evolve to become less steep, toward equipartition. Note that always x η < x , when x > 0. Now we can 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 is valid, the particle flux will be positive and the energy flux will be negative, contradicting the Fjørtoft argument. This means that one cannot match these formal KZ solutions, obtained 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 in optical wave turbulence [12] that the pure KZ spectra are not established in these cases and one has to expect a mixed solution where both a flux and a thermal components are present. Such mixed states are quite common for turbulent systems of different kinds, including strong Navier-Stokes turbulence and have been named warm cascades [19]. Such cascades were obtained within the Leith model (which belongs to the class of the differential approximation models) as exact analytical solutions.

III. DIFFERENTIAL APPROXIMATION MODEL
Numerical integration of the Boltzmann collision integral is very challenging because the number of degrees of freedom grows as a polynomial. A great simplification comes from the isotopic assumption, which reduces the degrees of freedom from N 8 to N 2 (N is the number of points needed to describe the distribution). However spanning a large number of momentum scales is still difficult. For those reasons, some approximations to the kinetic equations were proposed in order to increase the range of modeled scales, see for example [23].
A great simplification is to replace the collision integral operator of the kinetic equation by a nonlinear differential operator which mimics the basic scalings of the original one and yields the same steady solutions. The HIBE then results in a nonlinear partial differential equation called the differential approximation model (DAM). Such models have been proposed to simulate turbulence in different research fields: for example in water waves [24], in nonlinear optics [12], in strong Navier-Stokes turbulence [25,26], in Kelvin quantum turbulence [27], in astrophysics (Kompaneets equation) [28], in semiconductors [29]. Replacing the integral operator by a differential one amounts to assuming locality of the scale interactions, which means 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 spectra are 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 but distributions which are close to MB, warm cascades, which appear to be local. Thus, we use the DAM for describing this system, after which we will validate our results by computing the full HIBE.
For the dual cascade systems, such as gravity water waves [30], nonlinear Schrödinger equation [12], two-dimensional hydrodynamic turbulence [26], Kelvin waves [27] or HIBE considered here, DAM has always the form of a dual conservation law, where R is a nonlinear second-order differential term whose details depend on the particular model. This equation can be written as a continuity equation for the particle invariant, Moreover, equation (10) can be written as a continuity equation for the energy [12], We 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 where S is a constant. A formal derivation starting from the kinetic equation can be obtained following [12,29]. It is trivial to verify by substitution that KZ solutions (8) correspond to constant fluxes through scales. Namely, the KZ particle cascade has a constant particle flux and zero energy flux while the KZ energy cascade viceversa. Let us again consider the the flux directions on the KZ distributions, but now using DAM. Substituting power-law spectra n = c ω −x into (13), equations (11) and (12) yield By plotting η and 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 particle and the energy fluxes on the respective KZ solutions (x = 7/4 and x = 9/4) have wrong directions with respect to the Fjørtoft argument.
The beauty of the DAMs is the possibility to solve numerically the system for wide frequency ranges and, therefore, to find clear scalings. In particular, such models are very efficient for finding constant steady flux solutions because they become simple ordinary differential equations (ODEs). In the following we will present some analytical results for such steady states.

A. Constant energy flux: direct cascade
We will now find an ODE that describes a constant direct energy cascade with no flux of particles, which we call ODE-. According to Fjørtoft argument, this implies a large direct-cascade inertial range. Putting η = 0 in (11) and (12), we have constant energy flux =⇒ = R(ω, t) = const.
Using (13), we arrive at the following Cauchy problem where we have chosen the boundary conditions fixing the values of the distribution and its derivative at the same point ω 0 (e.g. at the forcing scale) for ease of numerical solution.
If we solve numerically in ω-forward the ODE-for different values of the energy flux we find curves presented in Fig. 7. Here we do not want to discuss the details (it will be done widely in Section IV), but just remark that the solutions follow the MB distribution and suddenly change behavior going very fast to a zero value of the distribution. We will call this rapid change a front solution.

Compact front behavior
It is possible to find a front solution for the equation (15) Thus, the front solution is linear in the vicinity of ω max with a slope depending on the dissipation scale ω max and the value of the energy flux . Note that the compact front behavior at the dissipation scale is typical for DAM. We will soon discover that ω max is a very useful physical parameter which allows us to find a link between the temperature, the chemical potential and the energy flux in the forced-dissipated system.

Kats-Kontorovich correction
Lets summarize our preliminary observations. We expect a warm cascade, that is a distribution which contains both the flux and the thermal components. We have also found that the solution has a compact front which arrests the cascade at the dissipation scale ω max . We will now assume (verifying it later) that in the most of the inertial range the warm cascade solution is close to the thermodynamic MB distribution and the correction due to finite flux is small. We then perform a qualitative matching of the flux-corrected MB 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: whereñ is small,ñ 1. By plugging this solution into (15) and linearizing inñ we end up with the following ODE-for the correction

Matching
We will now match the KK correction to the front solution. The basic idea is to force the KK solution to satisfy the n(ω max ) = 0 and to have at ω max the same slope as the front solution. Detailed calculation is presented in Appendix C. The prediction results in: This relation is very important because it gives an analytical relation between the thermodynamic quantities T and A in terms of the energy flux and the dissipation scale ω max .
However we note that our matching is only qualitative, because the KK correction is supposed to be small which is not the case near the front. Thus, the relation (20) is approximate and we do not expect it to hold precisely.

Alternative approach to find ω max
Another simple way to find a prediction for the value of ω max is the following. As we expect to observe a warm cascade, we can ask what will be the range where the thermal component will dominate the dynamics. We can simply assume that in most of the inertial range we will have a distribution n n M B . Note that the MB distribution always has a positive concavity, ∂ ωω n ≥ 0. On the other hand, we note that our ODE-can be re-written as ∂ ωω n = 1 n (∂ ω n) 2 − S ω 13 2 , from which it is clear that ∂ ωω n may change sign. The point at which ∂ ωω n = 0 can be considered as s boundary separating the MB range (with negligible flux correction) and the front solution (with large flux correction). This boundary can be estimated by a simple substitution of the MB distribution to the r.h.s. of (21), which gives As this relation contains the exponential factor which decays very fast (for ω max T , see Appendix C), it is natural to think that the range at which becomes important appears very sharply and is very near to the point ω max . Thus we arrive at the following estimate,

B. Constant particle flux: inverse cascade
In analogy of what has been done for the direct cascade, we now look for predictions in the inverse particle cascade η with no flux of energy. The ODE-η that describes such a cascade is simple to obtain: by integrating equation (11) once and putting = 0 in (12), we have: This yields the following Cauchy problem, This problem is most naturally solved backwards in the ω-space, as we are interested in the inverse cascade. We seek for a solution having a particle flux going from high to low frequencies, i.e. η < 0 and for convenience we will make the substitution η → −|η| in equation (25). The Cauchy problem (25) is very similar to (16) with the only difference in the ω-scaling. Thus we will use the same approach for studying it.

Compact front behavior
Let us find a front solution for the equation (24). We now expect the front to be on the 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 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 η.

Kats-Kontorovich correction
As previously supposed for the direct energy cascade, we expect in the most of the inversecascade range a corrected thermodynamic spectrum and a front solution behavior at the left end of this range. Let us evaluate the Kats-Kontorovich correction (18), and after that match it to the front solution. By plugging the expression (18) into (24) and linearizing iñ n we obtain the following ODE-η for the correction,

Matching
Again, we want to match the KK correction to the front solution. The idea is very similar to the previously used for the direct cascade, except for the fact that now the limit taken is

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 Keeping in mind that the MB distribution is always characterized by a positive concavity, i.e. ∂ ωω n ≥ 0, and considering the hypothesis ∂ ω n ∂ ω n M B we find Similarly to what we have done for the inverse cascade, we now can suggest that the change of concavity occurs near ω min . This results in However, we do not expect a good prediction as before because in this case the exponential term 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 , dissipate at ω max and ω min , and consider the case ω min ω f ω max . If the forcing range is narrow, the simple relation = η ω f holds for the fluxes. Using this relation, and combining (20) and (28), we can estimate T and A in the system: and, therefore, the chemical potential Note that the temperature appears to be independent of the fluxes and is completely controlled by the forcing and the dissipation scales. This means that increasing the forcing strength without moving ω f simply adds more particles into the system with the energy per particle remaining the same.

IV. NUMERICAL RESULTS
In this Section we present the numerical results obtained by using the DAM and by integrating, at lower resolution, the HIBE. Our aim is to compare results for the warm cascade solutions of DAM, which has been devised as a local approximation of the integral collision 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 simulations can be performed without any loss of generality starting with a particular value ω 0 because of re-scaling properties described in Appendix E.

Constant direct energy cascade
In Fig. 7 we show the results obtained by integrating equation (16) Fig. 7 within 5%. Note that the peak of g (ω, 1, 1, ) is around ω = 3.5: this is why we set this value as initial condition ω 0 .  In Fig. 9 we present the results for a particular case with flux = 1. We can appreciate the presence of warm cascade and the front solution near ω max . The linear behavior of the front is evident 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 measured  (20) is an overestimation of the numerical results (dots) and the error is around 10%. Finally prediction for the alternative temperature relation (23) is plotted with gray dashed line: it appears to give a better estimation than relation (20).

Constant particle cascade
We now investigate the inverse particle cascade by solving Cauchy problem (25) going ω-backward. In Fig. 11  flux, here the warm cascade range is wider for smaller flux values. We also observe fronts in vicinities of cutoff points ω min . In Fig. 12 we show the function g η (ω, 1, 1) which represents the 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

Double cascade
An example of double cascade is presented in Fig. 15 where we set the forcing at ω f = ω 0 = 3.5. We show here three cases where the particle fluxes are respectively η = −1, η = −10 −2 and η = −10 −4 . Measuring ω min and ω max for each case we are able to estimate the 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 ω f ω max is better satisfied (see for example the case η = −10 −4 ).

B. HIBE results
We now to present results of the direct simulation of HIBE with the full Boltzmann collision integral and compare them with predictions obtained by DAM. As we have mentioned above, the evaluation of (4) is numerically challenging and it is nowadays practically impossible to simulate such wide ω-space ranges as we have done using the DAM. In the present work, we will always use a low resolution of 101 points by considering ω ∈ [0, ω cutof f ] and taking a uniform distribution with ∆ω = ω 2 = ω 3 + ω 4 which can be pre-computed. Note that the dissipation at high wave numbers is chosen to satisfy ω max ≤ ω cutof f /2 in order to prevent ultraviolet bottleneck effects. The time evolution is performed by using the Euler scheme. Further details on numerical methods for solving the HIBE and a simple code can be found in [31].

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 for these 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 predictions (32). As expected the quantity A ∼ √ η but the line (in log-log plot) is shifted with respect to the interval between A L and A R , represented respectively with filled and empty circles. However, the theoretical prediction is much closer to A R , which is natural because 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 in the 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 temperature evaluated 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-ω dissipation range and results are presented in Fig. 18. Keeping the forcing constant and changing the value 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, is in between of the two temperatures and is closer to T R -again due to the wider right range.  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 are shown in Fig. 19 and for A in Fig. 20, with the "left" quantities shown by filled circles and the "right" ones by empty circles. There is a reasonably good agreement of T with the prediction (32) for small ω min . This is natural because smaller ω min corresponds to larger inverse cascade inertial range and also because the prediction is valid when ω min 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 the agreement can be made more suitable since the analytical prediction contains the undefined order-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 production rate is estimated in 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 directions must be large (note that this is also the condition of validity of the Fjørtoft argument).  We have built an ad-hoc differential approximation model to easily simulate the cascade processes. Indeed, this simplification allowed us to reach a wide range of scales inaccessible by solving the isotropic Boltzmann kinetic equation directly. Simulations show the presence of a warm cascade with approximately the MB shape followed by a sharp front for both energy and particle cascades. We have physically interpreted ω min and ω max as intrinsic dissipation scales 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 and dissipations scales to the thermodynamic quantities of the system. In particular we have shown that the temperature is independent of the amplitude of the fluxes but only depends on the forcing and dissipation scales.
We have then compared the theoretical predictions and the numerical results obtained with the differential approximation model with simulations of the complete isotropic Boltz- where geometrically k min = |k 1 − k 2 | = |k 3 − k 4 | and k max = k 1 + k 2 = k 3 + k 4 . For details about 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 where the two δ-functions assure the conservation of the linear momentum and kinetic energy.
In the following we will suppose that the particle distribution function follows the powerlaw distribution n(ω) = A ω −ν and so where Θ is the Heaviside step function.

Kolmogorov-Zakharov solutions
We will present the Kolmogorov-Zakharov solutions of the collision integral using the method presented by Balk in [22]. The collision integral, without any loss of generality, can be rewritten as where the exponent is chosen in order to have zero as homogeneity coefficient of the integrand (excluding the differentials dω i ω i ). If the integral converges, Balk proved that is possible to interchange the three integration index in the integrand with the fourth one, ω 1 . Thanks to the symmetric properties of the collision kernel we can write which clearly vanishes for µ = 0 or µ = 1. This corresponds to the condition on the exponent ν 0 = ν| µ=0 = 3d − 2 4 + β ν 1 = ν| µ=1 = 3d 4 + β. (B12) Note that first KZ solution for HIBE were presented in [17].
∀ i = 1 in the equation (B11) which results, recalling the homogeneity of the integrand function, in The integral U (µ) is now performed over the triangle ∆ in the ξ 3 × ξ 4 space satisfying the conditions 0 ≤ ξ i ≤ 1 and ξ 4 ≥ 1 − ξ 3 , without any dependence on ω 1 .
Finally, considering that ω ω max , we recover and validate the relation = S ω 9 2 max A 2 e − 2ωmax T . (C5) After some easy algebra we find that the system is invariant if α = (p − 2)β + 2γ.
As a consequence we can establish how the thermodynamic quantities defined by the Maxwell-Boltzmann distribution n M B (ω) = Ae − ω T = e − ω+µ T vary: the temperature T scales as ω and so T = λ βT , while the chemical potential µ scales as µ = λ β γμ.